# 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?