Problem with decarboxylation reactions

[Geert38]

Hi,

I am new to autodE. Many reactions work really smoothly sur as dehydration, hydration, deamination and ring closures. Others are more tricky. I am trying to get decarboxylation reactions to work with little success. I have tried many different. Here is my code for the decarboxylation of acetic acid. Now this reaction is supposed to be difficult with a high activation energy. I suspect that there is something more going on than just a difficult reaction. You say you are interested in a test case, so here is one. Is is possible that there is a thing going on with CO2?

import autode as ade
ade.Config.n_cores = 1
AceticAcid = ade.Reactant(smiles='CC(=O)O')
Methane = ade.Product(smiles='C')
CO2 = ade.Product(smiles='O=C=O')
rxn = ade.Reaction(AceticAcid, Methane, CO2, solvent_name='water', name='AA-Decarboxilation')
rxn.calculate_reaction_profile(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("∆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])

The energy E is nicely calculated but it does not calculate the free energy G. Actually it gives a value of 4 for dGts. The graph shows a warning that it assumes a barrierless reaction. When trying to read the value of the imaginary frequency it crashes, as it probably does not exist. I suspect that it has not located the transition state.
I have tested a variety of techniques, Orca, Gaussian, B3LYP, MP2, HF, solvent or not, but it never works.
The same is true for decarboxilation of formic acid (half of the WGS reaction), glycine, malonic acid and that of glutamic acid. Malonic acid is supposed to be easy to decarboxylate.

I work with
Python 3.9 on Windows10
Orca 5.04, Gaussian09
autodE version 1.4.1
rdkit version 223.9.5

Any thoughts ? It looks like autodE does not like CO2 as most other reactions work fine. CO2 probably does not react well to the search of confomers. Would this be a problem in the search for a TS?

Thanks in advance.

[t-young31]

Hi @Geert38

I've just tried it out and indeed autodE doesn't do a great job. Part of the problem is this early-stopping criteria on the initial path

autode.log.log: WARNING  Have a peak, products not made on isomorphism, but at least one of the distances is final. Assuming the peak is suitable    

which picks up a proton transfer across the carboxyl. We can force matters with a little manual intervention by creating reactant and products in the correct conformation

import autode as ade
ade.Config.lcode = ade.Config.hcode = "orca"

r = ade.Reactant("r.xyz")
p = ade.Product("p.xyz")
rxn = ade.Reaction(r, p, name='neb')
rxn.optimise_reacs_prods()
rxn.locate_transition_state()
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"))

8
r.xyz
C         -0.94699       -0.07829       -0.14189
C          0.49717        0.30811       -0.03089
O          0.96924        1.36670       -0.34881
O          1.23865       -0.68632        0.48331
H         -1.51414        0.74203       -0.59422
H         -1.34656       -0.30848        0.85676
H         -1.04552       -0.99068       -0.74768
H          0.53828       -1.44002        0.72219

8
p.xyz
C         -1.44659       -0.74815        0.05192
C          1.87467        0.88434       -0.02008
O          1.21899        1.79465       -0.43515
O          2.53012       -0.02628        0.39507
H         -1.21104        0.26014       -0.34547
H         -2.18369       -0.65644        0.87580
H         -1.87258       -1.37341       -0.75929
H         -0.51975       -1.22181        0.43597

A relevant aside: I'd be very surprised if this sort of decarboxylation is elementary with innocent solvent (consistent with the >60 kcal mol-1 barrier at DFT). I'd imagine some other acid/base or at least water mediated pathway is much more feasible.

[Geert38]

Dear Tom, Thanks for that. I will test this and post a response, it may help others. Yes of course this reaction is not so easy. The thing is that I can readily do this in the lab at 300-350 °C and get confirmation of the results. The gas is easily analysed and the water phase should only have acetic acid and so a simple COT measurement will give me the conversion without laborious HPLC analysis. I like simple things. I was suspecting an intermediate reaction but it is not so easy to find one. So far autodE did a great job in finding pathways that work, so I was surprised that it did not work on this one and on any of the same family. The library suggests that it is possible to make quantum chemistry simple, but this is deceptive, it will never be simple. I did actually test the following, but I suspect this is also too complex: Acetate = ade.Reactant(smiles='CC(=O)[O-]') Oxonium = ade.Reactant(smiles='[OH3+]') Methane = ade.Product(smiles='C') CO2 = ade.Product(smiles='O=C=O') Water = ade.Product(smiles='O') rxn = ade.Reaction(Acetate, Oxonium, Methane, CO2, Water, solvent_name= 'water', name= 'Acetate' ) rxn.calculate_reaction_profile(with_complexes=True) I joined the resulting pdf file. It shows that this may not be a simple reaction. I will go to work with your solution and I will let you know. The thing is that the documentation is a bit on the light side for a complex tool. I asked for access to the slack some time ago without response. By digging through old Q&A discussions we sometimes find answers to our questions. Have a good weekend, Geert

…

On Fri, 29 Mar 2024 at 13:40, Tom Young ***@***.***> wrote: Hi @Geert38 <https://github.com/Geert38> I've just tried it out and indeed autodE doesn't do a great job. Part of the problem is this early-stopping criteria on the initial path autode.log.log: WARNING Have a peak, products not made on isomorphism, but at least one of the distances is final. Assuming the peak is suitable which picks up a proton transfer across the carboxyl. We can force matters with a little manual intervention by creating reactant and products in the correct conformation import autode as adeade.Config.lcode = ade.Config.hcode = "orca" r = ade.Reactant("r.xyz")p = ade.Product("p.xyz")rxn = ade.Reaction(r, p, name='neb')rxn.optimise_reacs_prods()rxn.locate_transition_state()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")) 8r.xyz C -0.94699 -0.07829 -0.14189 C 0.49717 0.30811 -0.03089 O 0.96924 1.36670 -0.34881 O 1.23865 -0.68632 0.48331 H -1.51414 0.74203 -0.59422 H -1.34656 -0.30848 0.85676 H -1.04552 -0.99068 -0.74768 H 0.53828 -1.44002 0.72219 8p.xyz C -1.44659 -0.74815 0.05192 C 1.87467 0.88434 -0.02008 O 1.21899 1.79465 -0.43515 O 2.53012 -0.02628 0.39507 H -1.21104 0.26014 -0.34547 H -2.18369 -0.65644 0.87580 H -1.87258 -1.37341 -0.75929 H -0.51975 -1.22181 0.43597 A relevant aside: I'd be *very* surprised if this sort of decarboxylation is elementary with innocent solvent (consistent with the >60 kcal mol-1 barrier at DFT). I'd imagine some other acid/base or at least water mediated pathway is much more feasible. — Reply to this email directly, view it on GitHub <#331 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BG22KLWJ5HOAD3QY4GVL2LDY2VHKTAVCNFSM6AAAAABFMXTQP2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDAMRXGE4DQOJUGU> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[t-young31]

The library suggests that it is possible to make quantum chemistry simple,
but this is deceptive, it will never be simple.

Great quote. Should get it on a T-shirt!

I did actually test the following, but I suspect this is also too complex:

I'd suspect that too. A nano-reactor approach (e.g. 10.1021/acs.jctc.9b00143) is probably the best

The thing is that the documentation is a bit on the light side for a
complex tool

Fair comment. Unfortunately writing documentation never reaches the top of the priority list. Contributions are very welcome!

I asked for access to the slack some time ago without response

Sorry for missing this. I'll dig it out

[Geert38]

Hi, I have been running quite a few tests. Combining the products in one file does help convergence. None of these are the magic bullet, but at least now I am convinced that CO2 is not the orgin of the problem. I was somewhat surprised to find that the code you proposed: rxn.optimise_reacs_prods() rxn.locate_transition_state() Produces different results from my original code: rxn.calculate_reaction_profile() The results are of course very close, but not exactly the same. Sometimes one fares better, sometimes the other, there is no clear winner. Also it appears that there are quite a few algorithms that it can use. The files and directory it creates during the search for the transition state vary considerably between cases. Thanks anyway for the help.

…

On Sun, 31 Mar 2024 at 18:21, Tom Young ***@***.***> wrote: The library suggests that it is possible to make quantum chemistry simple, but this is deceptive, it will never be simple. Great quote. Should get it on a T-shirt! I did actually test the following, but I suspect this is also too complex: I'd suspect that too. A nano-reactor approach (e.g. 10.1021/acs.jctc.9b00143) is probably the best The thing is that the documentation is a bit on the light side for a complex tool Fair comment. Unfortunately writing documentation never reaches the top of the priority list. Contributions are very welcome! I asked for access to the slack some time ago without response Sorry for missing this. I'll dig it out — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you were mentioned.Message ID: ***@***.***>