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

Inconsistent molecular structure before and after CREST #8

Closed
jungsdao opened this issue Aug 1, 2020 · 3 comments
Closed

Inconsistent molecular structure before and after CREST #8

jungsdao opened this issue Aug 1, 2020 · 3 comments

Comments

@jungsdao
Copy link

jungsdao commented Aug 1, 2020

Hello, thank you for developing such a useful tool like gfn2-xtb and crest.
I'm using crest version 2.9 along with gfn2-xtb version 6.2.3.

While using crest for conformer search of a hydrocarbon chain, I've discovered that original molecular structure has not been preserved after crest job.
Following is an example.
"cal.xyz" is initial molecule before conformer search and "crest_best.xyz" is the output of crest job.
As can be found in the attached xyz file, weird CO2 molecule has been produced after crest.
Also following is my input command to run the job.

"crest cal.xyz -uhf 1 "
I want molecular structure (molecular connectivity) to be preserved but it seems to give false structure.
Probably the output structure might be thermodynamically more stable compared to input structure but this result defeats the purpose of conformer search.
Could anyone explain the reason or how I can circumvent this?
Many thanks in advance!

Crest_job.zip

@jungsdao jungsdao changed the title Inconsistenty molecular structure before and after CREST Inconsistent molecular structure before and after CREST Aug 1, 2020
@awvwgk awvwgk transferred this issue from grimme-lab/xtb Aug 1, 2020
@pprcht
Copy link
Contributor

pprcht commented Aug 3, 2020

Hi,
if you use the latest version from here, there is a command line flag called -cbonds that puts a constraint on each bond. This can help to maintain the connectivity and prevent dissociation. The -cbonds command accepts a force constant (in Eh) as an optional second argument, i.e., -cbonds <force constant>, you want to keep that as small as possible. Typical values are between 0.1 and 0.01 Eh.

In your calculation the reason for the dissociation could be that it is an open shell system (-uhf 1) in combination with a too large force due to the metadynamics potential. There might have been unwanted rearrangements induced by the RMSD potential at the respective open shell position, which lead to the dissociation.

I hope this will help you.

@jungsdao
Copy link
Author

jungsdao commented Aug 12, 2020

Thank you for your comment!
I'm glad to notice that recent features have included bond constraints.
I have updated version of crest and xtb to 2.10.2 and 6.3.2 respectively, and tried -cbonds option with 0.01 as criteria.
(with this command : crest cal.xyz -uhf 1 -cbonds 0.01)
However, this produces another error like forming a new "unwanted" bond.
In the attached *.xyz structure file, weird cyclic ring has been produced after crest job.

As you mentioned it seems -uhf 1 option induces such problem during meta-dynamics simulation.
It seems original bonds are constrained but newly formed bonds are not prohibited.
Would be there any possible way to prevent this from happening?

Best regards,
Hyunwook

structures.zip

@pprcht
Copy link
Contributor

pprcht commented Aug 13, 2020

In principle, all constraint that can be defined for a xtb calculation can also be applied in the conformational search. The details are described in the documentation here and here. The -cbonds option also just writes bond constraints in this format and applies them to the calculation. However, preventing something from froming new bonds while retaining an overall flexibility is one of the more difficult things to achieve with constraints. It is maybe more save to first try out constraints in xtb calculations and only afterwards -if they worked- include them in the conformational sampling.

If you are only interested in the structures (and not so much the energies) you could also try to conduct the search at the GFN-FF level with the -gfnff option and omit the uhf option. Because GFN-FF is a classical force field, bonds are defined by definition and no new bond formation should happen. The calculation will also be a lot faster.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants