Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add torsiondrive as a potential torsion scan tool. #43

Merged
merged 52 commits into from
Apr 6, 2022

Conversation

xiki-tempula
Copy link
Contributor

@xiki-tempula xiki-tempula commented Mar 25, 2022

@selimsami Do you mind help me on this?
I have allowed the torsiondrive to be used as torsion scan tool but it seems that the charge is not read in correctly.

https://github.com/Exscientia/qforce/blob/0fded8348d59a9967048f525c2fb766580162669/qforce/qm/qm.py#L122
I'm loading the charge from the torsion scan to be used as charge for the fitting.

But I found that the charge that I feed in to the Qforce (propane_qforce/fragments/CC_H8C3_d91b46644317dee9c2b868166c66a18c~1.charges is not the same as the topology that Qforce generates (Qforce (propane_qforce/fragments/CC_H8C3_d91b46644317dee9c2b868166c66a18c~1/step00/propane_qforce.itp) during the fitting.

Please find the example.
Archive.zip

@xiki-tempula
Copy link
Contributor Author

@selimsami Sorry, my fault. Forget to set

[ff]
charge_scaling = 1.0

@xiki-tempula xiki-tempula marked this pull request as draft March 28, 2022 10:54
@xiki-tempula
Copy link
Contributor Author

@selimsami The code is done now and I'm testing the robustness. Do you mind having a look and tell me your thought, please? Thanks.
Sorry that I need to change the sys.exit to raise SystemExit as otherwise my pycharm won't work https://youtrack.jetbrains.com/issue/PY-53625.
Sorry, I have to use some bash script again as torsiondrive command is a bit complicated and the output is an xyz file.

Copy link
Owner

@selimsami selimsami left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks a lot for this contribution, looks very nice. Sorry for the delay, had a really busy teaching week.

I made a few minor comments. And I have one more general one here: I think it would be nice to move xtb-torsiondrive reader/writer to a different module named torsiondrive or torsiondrive_xtb

qforce/qm/qm.py Outdated
Comment on lines 146 to 151
try:
os.makedirs(f'{dir}/{scan_id}_torsiondrive')
except FileExistsError:
warn(f'Folder {dir}/{scan_id}_torsiondrive exists, remove the folder.')
shutil.rmtree(f'{dir}/{scan_id}_torsiondrive')
os.makedirs(f'{dir}/{scan_id}_torsiondrive')
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why not just remove the directory if it exists?

if os.path.exists(name):
    shutil.rmtree(name)
os.makedirs(name)

qforce/qm/qm.py Outdated
os.makedirs(f'{dir}/{scan_id}_torsiondrive')

# Create the coordinate input file
with open(f'{dir}/{scan_id}_torsiondrive/input.xyz', 'w') as f:
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

perhaps you can use ASE here again to write the XYZ file?

See an example here.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@selimsami Thanks. I was thinking of the same thing but I couldn't get the mol.info to write what I want. It is nice that I could just put it in the comment.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@selimsami Sorry about more ase question. I wonder how do I create an ase mol object given ele_id and pos? see the input here https://github.com/selimsami/qforce/pull/43/files#diff-f579a382b4ca6ded95ab6837635635e6329d9fe38329927c5f4e2d8a4ff2b243R119

Copy link
Owner

@selimsami selimsami Apr 5, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

from ase import Atoms
mol = Atoms(positions=pos, numbers=ele_id)

@xiki-tempula
Copy link
Contributor Author

@selimsami Thanks for the advice on ase. Once I have tested that it is robust, I will ping you and we can have this merged.

@xiki-tempula xiki-tempula marked this pull request as ready for review April 6, 2022 09:32
@xiki-tempula
Copy link
Contributor Author

@selimsami I have tested this routine on 30 molecules and it seems fine. If you are happy, we can have this merged.

@selimsami selimsami merged commit a9b23a7 into selimsami:master Apr 6, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
2 participants