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-

  • 3
    $\begingroup$ 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. $\endgroup$ Feb 20, 2021 at 23:47
  • $\begingroup$ 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! $\endgroup$ Dec 25, 2021 at 21:18
  • $\begingroup$ @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). $\endgroup$
    – David CY
    Dec 27, 2021 at 3:13
  • $\begingroup$ @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. $\endgroup$ Dec 27, 2021 at 3:15
  • $\begingroup$ I also think a brief summary of the paper would make for an interesting answer. $\endgroup$
    – Tyberius
    Jan 26, 2022 at 4:52

1 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.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .