Clarification / Advice flag `with_complexes`

[pultar]

First of all, thank you for providing autode, it really is an amazingly useful tool!

I was wondering if you could comment on a couple of questions I had regarding its usage:

Why does the option with_complexes for the calculate_reaction_profile function produce a second maximum? I could not find any files associated with a second transition state in the folder structure.

WARNINGS: Reaction coordinate ∆E‡ not calculated for reactant_complex, barrierless reaction assumed.

also suggest there is no second transition state for associating reactant A with reactant B.

So my question is, if a barrierless reaction is assumed, how come I get very specific number in my plots suggesting the presence of a second transition state?

In general, do you recommend using this flag or omitting it and just have autode start from two separate xyz files, also having in mind the requirement to set the discouraged allow_association_complex_G flag?

I was also wondering about the structures found in output: Do I understand correctly that v1yjY7_reactant.xyz would be the reactant reactive complex and v1yjY7_product.xyz the product reactive complex (before relaxation / dissociating)? It is especially unclear to me what would be the meaning of a reactive product complex if there is only one product.

Maybe you could comment on the identity of each of the five transitionary points in my plot? Even after studying your paper and the documentation, I don't understand how autode identifies these points.

To Reproduce I attached some files of the simple reaction shown here:

test-complexes.tar.gz

import autode as ade
ade.Config.n_cores = 20
ade.Config.max_core = 4000
ade.Config.lcode = "xtb"
ade.Config.hcode = "orca"

if not ade.methods.ORCA().is_available and ade.methods.XTB().is_available:
    exit("This example requires an ORCA and XTB install")

ade.Config.allow_association_complex_G = True

ethylene = ade.Reactant(name='ethylene', smiles='C=C')
butadiene = ade.Reactant(name='butadiene', smiles='C=CC=C')
cyclohexene = ade.Product(name='cyclohexene', smiles='C1CCC=CC1')
rxn = ade.Reaction(ethylene, butadiene, cyclohexene, name='test', solvent_name='toluene')

rxn.calculate_reaction_profile(with_complexes=True, free_energy=True)
print("∆E_r (kcal mol-1) = ", rxn.delta("E").to("kcal mol-1"))
print("∆E‡ (kcal mol-1)  = ", rxn.delta("E‡").to("kcal mol-1"))
print("∆H_r (kcal mol-1) = ", rxn.delta("H").to("kcal mol-1"))
print("∆H‡ (kcal mol-1)  = ", rxn.delta("H‡").to("kcal mol-1"))
print("∆G_r (kcal mol-1) = ", rxn.delta("G").to("kcal mol-1"))
print("∆G‡ (kcal mol-1)  = ", rxn.delta("G‡").to("kcal mol-1"))
print("TS imaginary freq = ", rxn.ts.imaginary_frequencies[0])

Environment

• Operating System: Linux CentOS 7
• Python version: 3.10.6
• autodE version: 1.3.2

[pultar]

I think I can answer parts of my question:

https://github.com/duartegroup/autodE/blob/6813afd7523091352634f37339e75c9cf95f232f/autode/reactions/reaction.py#L371

For reaction complexes, I assume this is where energy estimated are made?

[t-young31]

Hi – sorry for taking so long to reply – too many side projects and not enough time!

Why does the option with_complexes for the calculate_reaction_profile function produce a second maximum? I could not find any files associated with a second transition state in the folder structure.

As you found it's a effective barrier rather than a real value. It probably should be temperature and solvent dependent but probably not a huge effect. Agree that there being a number on a plot make it rather unclear – perhaps it should be ommited if estimated?

In general, do you recommend using this flag or omitting it and just have autode start from two separate xyz files, also having in mind the requirement to set the discouraged allow_association_complex_G flag?

Unless there's a possibility the potential energy difference to the association complex is strongly negative (i.e more than the ~10 kcal mol-1 entropy barrier) then you'll get an equally accurate answer with with_complexes=False.

I was also wondering about the structures found in output: Do I understand correctly that v1yjY7_reactant.xyz would be the reactant reactive complex and v1yjY7_product.xyz the product reactive complex (before relaxation / dissociating)? It is especially unclear to me what would be the meaning of a reactive product complex if there is only one product.

You're totally right on this one. The product complex is the same as the product in this case and sadly the code isn't smart enough to omit printing it.

Maybe you could comment on the identity of each of the five transitionary points in my plot?

r0 = ethylene r1 = butadiene p = cyclohexene

r0 + r1 -> [no structure]‡ -> r0.r1 -> [TS]‡ -> p

[pultar]

Thank you so much! That explanation was very helpful!