Cusp correction

When Gaussian basis sets are used, molecular orbitals generally fail to satisfy the electron–nucleus cusp condition at nuclear positions. The casino.Cusp class implements the correction scheme of Ma, Towler, Drummond and Needs [12].

Physical motivation

The exact wavefunction must satisfy the Kato cusp condition [8] at every nucleus:

\[\left.\frac{\partial \langle\Psi\rangle}{\partial r_{iI}}\right|_{r_{iI}=0} = -Z_I \left.\langle\Psi\rangle\right|_{r_{iI}=0}\]

where \(Z_I\) is the nuclear charge, \(r_{iI}\) is the electron–nucleus distance and \(\langle\Psi\rangle\) is the spherical average of the wavefunction about the nuclear position (the spherical average was later shown to be redundant [15]). For a determinant of orbitals it is sufficient that every orbital which is nonzero at a given nucleus obeys the cusp condition there [12]. Gaussian functions have zero derivative at the nucleus (they are smooth at the origin), so a linear combination of Gaussians cannot satisfy this condition directly: the local energy then diverges as \(-Z_I/r_{iI}\) when an electron approaches a nucleus and shows wild oscillations close to it, which increase the VMC variance and can lead to severe bias or even numerical instabilities in DMC [12].

The Kato condition relates only the value and the first radial derivative of the wavefunction at the coalescence point; higher-order coalescence conditions, constraining the second and higher derivatives, are also known [13] [14] [16].

Correction scheme

Each molecular orbital \(\psi\) expanded in a Gaussian basis can be written as:

\[\psi = \phi + \eta\]

where \(\phi\) is the contribution from s-type Gaussians centered on the nucleus in question and \(\eta\) contains all other contributions. The corrected orbital is:

\[\tilde{\psi} = \tilde{\phi} + \eta\]

Inside a cusp correction radius \(r_c\), the s-part is replaced by:

\[\tilde{\phi} = C + \text{sgn}[\tilde{\phi}(0)] \exp(p(r)) \equiv C + R(r)\]

where \(C\) is a shift chosen so that \(\tilde{\phi} - C\) is of one sign within \(r_c\), and \(p(r)\) is a quartic polynomial:

\[p(r) = \alpha_0 + \alpha_1 r + \alpha_2 r^2 + \alpha_3 r^3 + \alpha_4 r^4\]

Since \(\eta\) is cuspless (it arises from Gaussians of nonzero angular momentum centered on the nucleus in question and from tails of Gaussians centered on other sites), the cusp condition for the corrected orbital reads [12]:

\[\left.\frac{d\tilde{\phi}}{dr}\right|_{r=0} = -Z\left(\tilde{\phi}(0) + \eta(0)\right)\]

The five coefficients \(\alpha_i\) are fixed by five constraints [12]: continuity of \(\tilde{\phi}\) and of its first and second derivatives at \(r = r_c\), the cusp condition at \(r = 0\), and a chosen value of \(\tilde{\phi}(0)\). With these constraints the equations for \(\alpha_i\) are solved analytically, and \(\tilde{\phi}(0)\) is then varied to minimise the maximum square deviation of the effective one-electron local energy

\[E_L^s(r) = \tilde{\phi}^{-1}\left[-\frac{1}{2}\nabla^2 - \frac{Z_{\text{eff}}}{r}\right] \tilde{\phi}, \qquad Z_{\text{eff}} = Z\left(1 + \frac{\eta(0)}{C + R(0)}\right)\]

from an “ideal” curve \(E_L^{\text{ideal}}(r) = Z^2(\beta_0 + \beta_1 r^2 + \beta_2 r^3 + \ldots + \beta_7 r^7)\), whose coefficients \(\beta_1 \ldots \beta_7\) were obtained by fitting to the 1s orbital of the carbon atom (\(\beta_0\) depends on the particular atom and its environment) [12].

Choice of the correction radius

The maximum possible cusp correction radius is \(r_{c,\max} = 1/Z\). The actual \(r_c\) for each orbital and nucleus is set to the largest radius less than \(r_{c,\max}\) at which the deviation of the effective one-electron local energy of the uncorrected orbital from the ideal curve exceeds \(Z^2/c_c\), where \(c_c\) is a universal parameter with a default value of 50 [12].

Enabling cusp correction

Cusp correction is applied automatically when atom_basis_type : gaussian is set in the input file (it is disabled for Slater-type orbitals, which can satisfy the cusp condition exactly).

To obtain the cusp polynomial coefficients from CASINO for comparison, add to your input:

cusp_info : T

and set POLYPRINT=.true. in gaussians.f90 before compiling CASINO.

Initialisation

casino.Cusp is created by casino.cusp.CuspFactory when cusp correction is enabled and passed to casino.Slater. It is not normally instantiated directly by the user. The correction is applied transparently inside slater.value_matrix, slater.gradient_matrix, and higher derivatives.

Properties

orbital_sign

Array of shape (natom, norbitals) storing the sign of each corrected orbital at each nuclear position.

shift

Array of shape (natom, norbitals) storing the constant shift \(C\) for each orbital–nucleus pair.

rc

Array of shape (natom, norbitals) storing the cusp correction radii.

alpha

Array of shape (5, natom, norbitals) storing the quartic polynomial coefficients \(\alpha_0 \ldots \alpha_4\).

References

[8]

T. Kato, On the eigenfunctions of many-particle systems in quantum mechanics, Commun. Pure Appl. Math. 10, 151 (1957).

[12] (1,2,3,4,5,6,7)

A. Ma, M. D. Towler, N. D. Drummond, and R. J. Needs, Scheme for adding electron–nucleus cusps to Gaussian orbitals, J. Chem. Phys. 122, 224322 (2005).

[13]

V. A. Rassolov and D. M. Chipman, Behavior of electronic wave functions near cusps, J. Chem. Phys. 104, 9908 (1996).

[14]

Y. I. Kurokawa, H. Nakashima, and H. Nakatsuji, General coalescence conditions for the exact wave functions: higher-order relations for Coulombic and non-Coulombic systems, Adv. Quantum Chem. 73, 59 (2016).

[15]

P. V. Tóth, Boundary conditions for many-electron systems, arXiv:1010.2700 [quant-ph] (2011).

[16]

J. Karwowski and A. Savin, Two-particle coalescence conditions revisited, arXiv:2204.09897 [physics.chem-ph] (2022).