duartegroup / autodE

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

How to get coordinates of reactant and product before conducting NEB? #284

Closed KS-spec closed 1 year ago

KS-spec commented 1 year ago

Dear all,

I'd like to get coordinates of reactant and product before conducting NEB. For example, if I type the below command and search for conformer and then perform NEB calculation, how can I obtain the xyz files of the reactant, which contains two molecules, and the product, which contains one molecule, before conducting the NEB calculation?

(ex.) rxn = ade.Reaction('CC=O.[H]N(C)C>>OC(C)(N(C)C)[H]', name='enamine_1') rxn.calculate_reaction_profile()

t-young31 commented 1 year ago

Hi @KS-spec – you should be able to do that with the following:

import autode as ade 

rxn = ade.Reaction("CC=O.[H]N(C)C>>OC(C)(N(C)C)[H]", name="enamine_1") 
rxn.find_lowest_energy_conformers()
rnx.optimise_reacs_prods()

for molecule in rxn.reacs + rnx.prods: 
    molecule.print_xyz_file()
KS-spec commented 1 year ago

Thank you very much for your kind reply.

The method you provided gave me the coordinates of each of the three molecules separately. I'd like to get the first and last coordinate of a trajectory such as the below coordinate (1 and 2) before conducting NEB. Could you please tell me if there is any way to do so? image

shoubhikraj commented 1 year ago

@KS-spec

NEB requires exact atom-mapping, and obtaining atom-mapped reactant and product geometries (especially when there are multiple molecules in the reactant i.e. a complex) is non-trivial (requires matching stereochemistry, RMSD etc.). AutodE actually side-steps the problem by only using the bond distances of the forming and breaking bond, and then performing "adaptive" constrained optimisations starting from the reactant complex.

The adaptive path is then refined with a NEB if required. The best you can do, at this stage, is to use the first and last point of the adaptive path trajectory. You can find this in the "transition_states" sub-folder which is generated when you call reaction.locate_transition_state() or reaction.calculate_reaction_profile(). The file is usually named <hash>_ll_ad_<bond_rearr>_path.xyz or <hash>_hl_ad_<bond_rearr>_path.xyz. @t-young31 pls correct me in case I missed something.

KS-spec commented 1 year ago

Thank you very much for your reply.