How is a counterpoise corrected geometry optimization done?

Question

I understand the problem of basis set superposition error (BSSE) and I know how the counterpoise correction for single point energies is calculated.

Today I found out that many software packages allow for counterpoise correction during optimization calculations, but how does this actually work, especially for methods where analyical gradients are used for optimization?

During geometry optimization we calculate the first derivative of the energy to get energy gradients which we follow down to our minima. I understand that I could use counterpoise correction if I calculate those gradients numerically, which is quite easy to understand but very expensive to do, but it seems that counterpoise can also be used in combination with analytical gradients. How is the counterpoise implemented to get counterpoise corrected gradients?

Answer 1

Background

For a system consisting of two molecules (monomers or fragments are also used) X and Y, the binding energy is

$$ \Delta E_{\text{bind}} = E^{\ce{XY}}(\ce{XY}) - [E^{\ce{X}}(\ce{X}) + E^{\ce{Y}}(\ce{Y})] \label{eq:sherrill-1} \tag{Sherrill 1} $$

where the letters in the parentheses refer to the atoms present in the calculation and the letters in the superscript refer to the (atomic orbital, AO) basis present in the calculation. The first term is the energy calculated for the combined X + Y complex (the dimer) with basis functions, and the next two terms are energy calculations for each isolated monomer with only their respective basis functions. The remainder of this discussion will make more sense if the complex geometry is used for each monomer, rather than the isolated fragment geometry.

The counterpoise-corrected (CP-corrected) binding energy [1] to correct for basis set superposition error (BSSE) [2] is defined as

$$ \Delta E_{\text{bind}}^{\text{CP}} = E^{\ce{XY}}(\ce{XY}) - [E^{\ce{XY}}(\ce{X}) + E^{\ce{XY}}(\ce{Y})] \label{eq:sherrill-3} \tag{Sherrill 3} $$

where the monomer calculations are now performed in the dimer/complex basis. Let's explicitly state how this works for the $E_{\ce{XY}}(\ce{X})$ term. The first molecule X contributes nuclei with charges, basis functions (AOs) centered on those nuclei, and electrons that will count to the final occupied molecular orbital (MO) index into the MO coefficient array. There is no reason why additional AOs that are not centered on atoms can't be added to a calculation. Depending on their spatial location, if they're close enough to have significant overlap, they may combine with atom-centered MOs, increasing the variational flexibility of the calculation and lowering the overall energy. Put another way, place the AOs that would correspond to molecule Y at their correct positions, but don't put the nuclei there, and don't consider the number of electrons they would contribute to the total number of occupied orbitals. This means that for the full electronic Hamiltonian

$$ \hat{H}_{\text{elec}} = \hat{T}_{e} + \hat{V}_{eN} + \hat{V}_{ee} $$

calculating the electron-nuclear attraction $\hat{V}_{eN}$ term is now different. Considered explicitly in matrix form in the AO basis,

$$ \begin{align*} V_{\mu\nu} &= \int \mathop{d\mathbf{r}_{i}} \chi_{\mu}(\mathbf{r}_{i}) \left( \sum_{A}^{N_{\text{atoms}}} \frac{Z_{A}}{|\mathbf{r}_{i} - \mathbf{R}_{A}|} \right) \chi_{\nu}(\mathbf{r}_{i}) \\ &=\sum_{A}^{N_{\text{atoms}}} Z_{A} \left< \chi_{\mu} \middle| \frac{1}{r_{A}} \middle| \chi_{\nu} \right> \end{align*} $$

there are now fewer terms in the summation, since the nuclear charges from molecule Y are zero (the atoms just aren't there), but the number of $\mu\nu$ are the same as for the XY complex. This and the $\hat{T}_{e}, \hat{V}_{ee}$ terms aren't really mathematically or functionally different then; this is more to show where the additional basis functions enter, or to show where nuclei appear in the equations [3].

These atoms that don't have nuclei or electrons, only basis functions, are called ghost atoms. Sometimes you also see the term ghost functions, ghost basis, or ghost {something} calculation. Adding the basis of monomer Y to make the full "dimer basis" means taking the monomer X and including basis functions at the nuclear positions for Y.

Geometry optimization

Now to calculate the molecular gradient, that is, the derivative of the energy with respect to the $3N$ nuclear coordinates. This is the central quantity in any geometry optimization. For the sake of simplicity, consider a steepest descent-type update of the nuclear coordinates $$ R_{A,x}^{(n+1)} = R_{A,x}^{(n)} - \alpha \frac{\partial E_{\text{total}}^{(n)}}{\partial R_{A,x}} \label{eq:steepest-descent} \tag{Steepest Descent} $$ where $n$ is the optimization iteration number, $\alpha$ is some small step size with units [length²][energy], and the last term is the derivative of the total (not just electronic) energy with respect to a change in atom $A$'s $x$-coordinate. Even Newton-Raphson-type updates with approximate Hessians (2nd derivative of the energy with respect to nuclear coordinates, rather than the 1st) need the gradient, so we must formulate it.

Formulation of the energy

We're in a bit of trouble, because we want to replace $E_{\text{total}}$ in the gradient with $E_{\text{total}}^{\text{CP}}$, but all we have is $\Delta E_{\text{bind}}^{\text{CP}}$. The concept of CP correction can still be applied to a total energy, but the BSSE must be removed from each monomer. The BSSE correction itself for each monomer is $$ \begin{split} E_{\text{BSSE}}(\ce{X}) &= E^{\ce{XY}}(\ce{X}) - E^{\ce{X}}(\ce{X}), \\ E_{\text{BSSE}}(\ce{Y}) &= E^{\ce{XY}}(\ce{Y}) - E^{\ce{Y}}(\ce{Y}), \end{split} \label{eq:2} $$ which, when subtracted from $\eqref{eq:sherrill-1}$, gives $\eqref{eq:sherrill-3}$. More correctly, considering that the geometry for each step is at the final cluster geometry and not the isolated geometry, the above is [4] $$ \begin{split} E_{\text{BSSE}}(\ce{X}) &= E_{\ce{XY}}^{\ce{XY}}(\ce{X}) - E_{\ce{XY}}^{\ce{X}}(\ce{X}), \\ E_{\text{BSSE}}(\ce{Y}) &= E_{\ce{XY}}^{\ce{XY}}(\ce{Y}) - E_{\ce{XY}}^{\ce{Y}}(\ce{Y}). \end{split} \label{eq:sherrill-10} \tag{Sherrill 10} $$

The CP-corrected total energy is the full dimer energy with BSSE removed from each monomer is then $$ \begin{split} E_{\text{tot}, \ce{\widetilde{XY}}}^{\text{CP}} &= E_{\ce{\widetilde{XY}}}^{\ce{XY}}(\ce{XY}) - E_{\text{BSSE}}(\ce{X}) - E_{\text{BSSE}}(\ce{Y}), \\ &= E_{\ce{\widetilde{XY}}}^{\ce{XY}}(\ce{XY}) - \left[ E_{\ce{\widetilde{XY}}}^{\ce{XY}}(\ce{X}) - E_{\ce{\widetilde{XY}}}^{\ce{X}}(\ce{X}) \right] - \left[ E_{\ce{\widetilde{XY}}}^{\ce{XY}}(\ce{Y}) - E_{\ce{\widetilde{XY}}}^{\ce{Y}}(\ce{Y}) \right]. \end{split} \label{eq:sherrill-15} \tag{Sherrill 15} $$ Note that I have modified which geometry is used for each monomer in $\eqref{eq:sherrill-15}$. All monomers are calculated at the supermolecule geometry. This is convenient for two reasons: 1. We are only interested in removing the BSSE, not the effect of monomer deformation, and 2. a isolated monomer geometry without deformation doesn't make sense in the context of a geometry optimization. I also added the tilde to signify that the supermolecular/dimer geometry used may not be the final or minimum-energy geometry, as would be the case during a geometry optimization. We simply extract all structures consistently from a given geometry iteration. Perhaps $\ce{XY}(n)$ would be better notation.

Formulation of the gradient

As Pedro correctly states, the differentiation operator is a linear operator. Because there are no products in $\eqref{eq:sherrill-15}$, the total gradient needed for $\eqref{eq:steepest-descent}$ will be a sum of gradients [5]: $$ \frac{\partial E_{\text{tot}, \ce{\widetilde{XY}}}^{\text{CP}}}{\partial R_{A,x}} = \frac{\partial E_{\ce{\widetilde{XY}}}^{\ce{XY}}(\ce{XY})}{\partial R_{A,x}} - \left[ \frac{\partial E_{\ce{\widetilde{XY}}}^{\ce{XY}}(\ce{X})}{\partial R_{A,x}} - \frac{\partial E_{\ce{\widetilde{XY}}}^{\ce{X}}(\ce{X})}{\partial R_{A,x}} \right] - \left[ \frac{\partial E_{\ce{\widetilde{XY}}}^{\ce{XY}}(\ce{Y})}{\partial R_{A,x}} - \frac{\partial E_{\ce{\widetilde{XY}}}^{\ce{Y}}(\ce{Y})}{\partial R_{A,x}} \right], $$ so each step of a CP-corrected geometry optimization will require 5 gradient calculations rather than 1. Note that the nuclear gradient should be included for each term as well, which is a trivial calculation.

Extension to other molecular properties

Although not commonly done, counterpoise correction can be applied to any molecular property, not just energies or gradients. Simply replace $E$ or $\partial E/\partial R$ with the property of interest. For example, the CP-corrected polarizability $\alpha$ of two fragments is $$ \alpha_{\text{tot}, \ce{\widetilde{XY}}}^{\text{CP}} = \alpha_{\ce{\widetilde{XY}}}^{\ce{XY}}(\ce{XY}) - \left[ \alpha_{\ce{\widetilde{XY}}}^{\ce{XY}}(\ce{X}) - \alpha_{\ce{\widetilde{XY}}}^{\ce{X}}(\ce{X}) \right] - \left[ \alpha_{\ce{\widetilde{XY}}}^{\ce{XY}}(\ce{Y}) - \alpha_{\ce{\widetilde{XY}}}^{\ce{Y}}(\ce{Y}) \right] $$ where I believe it now makes even less sense to have each individual fragment calculation not be at the cluster geometry. In papers that calculate CP-corrected properties, no mention is usually made of which geometry the individual calculations are performed at for this reason.

---

References

1. Boys, S. Francis; Bernardi, F. The calculation of small molecular interactions by the differences of separate total energies. Some procedures with reduced errors. Mol. Phys. 1970, 19, 553-566.
2. Sherrill, C. David. Counterpoise Correction and Basis Set Superposition Error. 2010, 1-6.
3. One implementation note: Most common quantum chemistry packages should allow for the usage of ghost atoms in energy and gradient calculations. However, as Sherrill states, they do not properly allow for composing the full gradient expression to perform CP-corrected geometry optimizations. Gaussian can, and Psi4 may. For programs that can calculate gradients with ghost atoms, Cuby can be used to drive the calculation of CP-corrected geometries and frequencies.
4. There is a typo in the Sherrill paper; the subscripts for all 4 energy terms should be $AB$, which here are $\ce{XY}$.
5. Simon, S.; Bertran, J.; Sodupe, M. Effect of Counterpoise Correction on the Geometries and Vibrational Frequencies of Hydrogen Bonded Systems. J. Chem. Phys. A 2001, 105,, 4359-4364.

Answer 2

The CP-correction can be seen as a sum of a correction term to the total energy of the complex. That correction term is essentialy the diference between the energies of the fragments with their basis set and with the whole basis set at the complex's geometry. The total CP-corrected energy is thus a sum of different complex and monomer's energies. Then, as the derivative is a linear operator, any derivative of the total CP-corrected energy is a sum of derivatives of the individual energies.. as simple as that... You can obtain CP-corrected 1st derivatives wrt nuclear coordinates and get optimized structures, but also CP-corrected armonic (and anarmonic) frequencies from 2nd (and 3rd) derivatives, dipole moments, etc...