Two-stage unbinding & binding energies with DFT

Question

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?

Thanks in advance-

Answer 1

As per request from some interested readers, I'm leaving a summary of what we ended up concluding and publishing.

In our study, "A0" is the cluster of Zr and O atoms in the metal-organic framework named PCN-224, and "B" is one of the two drug molecules (nalidixic acid "NDX" or ofloxacin "OFL") that we're trying to bind to the Zr-cluster. The experiment is done with water as the solvent, though the solvent was acetic acid in earlier stages of the synthesis of PCN-224.

The Zr atoms in the cluster has some "exposed" sites that allow coordinate bonding with other molecules. In water, these sites are usually occupied by H2O...OH ("A1"), together making "A". When pitted against B, A1 is supposed to detach from A0 before allowing B to bind to A0, forming "C". As stated in my original question, removing A1 requires energy, which is supplied thermally from the experimental setup and is observed to further proceed to the formation of C. So we expected the DFT calculation to show C having lower energy than A, but got the opposite result.

We then asked ourselves whether the acetate molecules (conjugate base of acetic acid) from the earlier stages of the experiment still remain attached to the Zr-cluster. When the calculation was redone with A1 set to acetate, the energy to remove A1 from A turned out to be smaller and gave C an energy lower than A. So the whole process was energetically explained.

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