search

duartegroup / autodE

automated reaction profile generation
https://duartegroup.github.io/autodE/
MIT License
165 stars 51 forks source link

Implementation IRC Scan #259

Open pultar opened 1 year ago

pultar commented 1 year ago

Describe the feature

autode is very good in automating lots of the workflow for transition state finding. It does require some manual checking afterwards though to verify found transition states are indeed meaningful. I would find it very cool if autode could provide more help also at this stage by following the intrinsic reaction coordinate (IRC, https://onlinelibrary.wiley.com/doi/full/10.1002/qua.24757) and check that the transition state found indeed connects starting material and product. I have had some unfortunate cases in which a manual scan by me revealed that autode had found a saddle point instead.

What do you think?

t-young31 commented 1 year ago

Thanks for the suggestion!

There is a check that the TS connects reactants and products using a quick reaction coordinate, which was the preference over an IRC because (a) it's easter to implement and (b) computationally cheaper. That does however only get executed if there isn't a conclusive set of displacements in the imaginary mode (i.e. a breaking bond lengthens by x Å), which sadly does rely on some magic numbers.. Would be interested to know what breaks in the cases you've seen. Would perhaps be nice to have an option to always do a QRC and be a little more conservative – always a balance between speed and accuracy!

pultar commented 1 year ago

Thanks for the answer and sorry for the very late reply. The calculation that gave me trouble comes from this paper:

https://onlinelibrary.wiley.com/doi/epdf/10.1002/anie.201711372

In brief, the authors tried to access the skeleton of a natural product using a Diels-Alder reaction and four different substrates. The experimental observation was that only 2/4 reactions gave product.

In all four cases, when I was using autode (xtb, orca, RDKit backend), only 1 bond is formed during the transition state (see files attached).

nanini.zip

It is also worrying that during one reaction (1d) an isomerization of the olefin occurred (reactant vs transition state). I believe this might be due to a missing check in the RDKit conformer generator, which should verify that all olefin geometries are preserved during conformer search. I talked to the authors of the RDKit and it is a known bug already fixed in master.

I have compiled the RDKit from source and will try to verify if this bug fix improves the outcome of the calculation if the next RDKit release does not happen earlier. I just thought you guys wanted to know about this problem either way.