duartegroup / autodE

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

Defining solvents for transition states inconsistent with molecules #164

Closed ffmulks closed 2 years ago

ffmulks commented 2 years ago

Describe the bug Adding solvent variables seems to function differently in transition states than in other parts of the code. I understand that this usually will take its variables from the provided reactant and product but it is possible to use the script to run a standard TS optimisation straight from an xyz-file without reactant and product complexes. This is helpful for enabling workarounds for reactions that autodE cannot do at the moment (i.e. I am doing the rotation barrier calculations that I mentioned previously via manually guessed structures or via PES-scan guesses). Neither TSguess() nor TransitionState() eats "solvent_name". TSguess() does take, however, "charge" and "mult" so it appears inconsistent with how Molecule() works that "solvent_name" cannot be used. TransitionState() does take TransitionState().solvent, though (which is then the only part that was defined on this function and not on TSguess).

To Reproduce Assume we have a molecule defined.

tsguess = TSguess(atoms=molecule.atoms,molecule.mult
                                   charge=molecule.charge,
                                   mult=molecule.mult,
                                   solvent_name=molecule.solvent
                                   )
molecule = TransitionState(tsguess)
molecule.optimise(method=hmethod)

This throws an unexpected keyword.

tsguess = TSguess(atoms=molecule.atoms,molecule.mult
                                   charge=molecule.charge,
                                   mult=molecule.mult
                                   )
molecule = TransitionState(tsguess, solvent_name = molecule.solvent)
molecule.optimise(method=hmethod)

This throws an unexpected keyword.

tsguess = TSguess(atoms=molecule.atoms,molecule.mult
                                   charge=molecule.charge,
                                   mult=molecule.mult
                                   )
molecule = TransitionState(tsguess)
molecule.solvent = get_solvent(molecule.solvent, kind='implicit')
molecule.optimise(method=hmethod)

This works.

tsguess = TSguess(atoms=molecule.atoms,molecule.mult
                                   charge=molecule.charge,
                                   mult=molecule.mult
                                   )
tsguess.solvent = get_solvent(molecule.solvent, kind='implicit')
molecule = TransitionState(tsguess)
molecule.optimise(method=hmethod)

This one is dangerous, it did not throw an exception but the CPCM(acetonitrile) is nowhere to be found in the generated ORCA inputs. When my reaction energy was roughly 700 kJ/mol wrong, I got a little suspicious here.

Expected behavior Solvents should be possible to define where charge and mult is allowed as well.

t-young31 commented 2 years ago

You are completely correct – I've been meaning to fix this, so now maybe the time! Expect to see this in v1.3.0 😄


tsguess.solvent = get_solvent(molecule.solvent, kind='implicit')
molecule = TransitionState(tsguess)

This doesn't work because the solvent gets inherited from the optional 'reactant' argument to TransitionState rather than the ts_guess, which I very much agree is confusing and bad!

t-young31 commented 2 years ago

Hi @ffmulks – this has been fixed in #167. If you have time I'd very much appreciate if you could try it out and see if the behaviour matches what you'd expect – thanks 😄


working example:

from autode import Atom
from autode.transition_states import TransitionState, TSguess

ts_guess = TSguess(mult=2, solvent_name='water',
                   atoms=[Atom("H", 4.52368, -0.00030,   0.00065),
                          Atom("H", 5.47129,  0.00009,  -0.00023),
                          Atom("H", 6.41149,  0.00045,  -0.00109)])

ts = TransitionState(ts_guess)
assert ts.solvent is not None

ts.optimise()
t-young31 commented 2 years ago

see if the behaviour matches what you'd expect

I'm going to hope it has because i'm going to merge #167 now to try and get a release out shortly

ffmulks commented 2 years ago

Weird, that works for me whereas my script in SLURM failed to do the very same thing. I suppose I may just have an environment mixup somewhere.

t-young31 commented 2 years ago

👻