# Two-stage unbinding & binding energies with DFT

I'm trying to calculate the stability of the product of a coordinate bonding process of molecules A and B in the polarizable continuum solvent model (PCM, solvent set to water).
I'm getting energies that (at first sight) don't agree with experiment, so I'd like to check things with the readers here to make sure I didn't make any obvious mistakes.

In the molecular process, a part of A is removed (say A - A1 → A0) before having B bind to what remains (A0) where the fragment A1 was removed. Let's call the final supermolecule C.
The fragments A1 (once removed) and B belong to an ambient "reservoir" -- we assume no interaction between them.

I understood this to mean that I only need to calculate the changes of energy through two reactions:

• A → A0 + A1
• A0 + B → C

Here are the specifics of the calculations I did. I used Gaussian 16 for DFT with the settings PBE-D3(BJ)/Def2SVP and SCRF(PCM) to optimize the geometries of A, A0, A1, B, and C, all separately.

Next, I defined the fragments A={A0, A1} and C={A0, B} to additionally calculate the counterpoise corrections for BSSE. Since counterpoise and SCRF can't go together (at least in g16), I calculated the corrections (single-point, no further optimization) using only the geometries of A and C that were optimized from the last step while not using the SCRF(PCM) command in this step.

After all this, the relative energies came out as

• A + (192.5 kJ/mol) → A0 + A1
• A0 + B → C + (162.1 kJ/mol)

which seem to mean that A is 30.4 kJ/mol more stable than C.

But my experiment somehow suggests that it favors the formation of C after having started from A.

So before I go into theorizing about any intermediate processes that I may have overlooked, can anyone point out any obvious mistake (either with concepts or software) that I may have made?

• As a brief note, PBE-D3(BJ)/def2-SVP is not a particularly accurate level of theory. The basis set is small, and there are likely better functionals for the task you're interested in modeling. The level of theory could very well explain the difference you're seeing. Feb 20, 2021 at 23:47
• David, have you had any luck or update with this question? It's been almost a year in the unanswered queue so I'd like to see some update if possible! Dec 25, 2021 at 21:18
• @NikeDattani it looks like there were no obvious mistakes and I concluded that there indeed was an overlooked intermediate process. The trend didn't change even across a few other higher level functionals that I tested, but ended up publishing the PBE data for better reproducibility (or ease of verification). Dec 27, 2021 at 3:13
• @DavidCY Would you be able to write a self-answer with a link to that paper and a summary of your findings? The question can still remain open, it would just be nice to have that information in an answer since comments are temporary and can be deleted without us being notified. Dec 27, 2021 at 3:15
• I also think a brief summary of the paper would make for an interesting answer.
– Tyberius
Jan 26, 2022 at 4:52

PS. Though not stated in the paper, changing PBE-D3(BJ)/Def2SVP to M06L-D3/Def2TZVPP didn't change our result.