The McMurchie-Davidson integral scheme

Since a two-center Gaussian overlap distribution, can be rewritten as a one-center Gaussian, we might wonder if we can use the Hermite Gaussian recurrence relations to simplify the calculation of molecular integrals. Since the overlap distribution \(\Omega_{ij}(x)\) is in essence a Gaussian times a polynomial in \(x\) with degree \(i+j\), we may express it as a linear combination of Hermite Gaussians: \[\begin{equation} \require{physics} \Omega_{ij}(x) = \sum_t^{i+j} E_t^{ij} \Lambda_t(x_P) \thinspace , \end{equation}\] where \(E_t^{ij}\) are are the expansion coefficients and we must have \[\begin{equation} t >= 0 \qquad \text{and} \qquad t <= i + j \thinspace . \end{equation}\] We may then establish the following recurrence relations for the expansion coefficients of the overlap distribution: \[\begin{align} &E^{i+1,j}_t = \frac{1}{2p} E^{ij}_{t-1} + (t+1) E^{ij}_{t+1} - \frac{\beta}{p} \Delta_{x, KL} E^{ij}_t \\ % &E^{i,j+1}_t = \frac{1}{2p} E^{ij}_{t-1} + (t+1) E^{ij}_{t+1} + \frac{\alpha}{p} \Delta_{x, KL} E^{ij}_t \thinspace , \end{align}\] where the first coefficient can be calculated as \[\begin{equation} E_0^{00} = C_x^{KL} \thinspace . \end{equation}\]

Overlap integrals

The overlap integrals over the Cartesian Gaussian primitives \(G_{ikm}\) and \(G_{jln}\) can be factored: \[\begin{equation} \braket{G_{ikm}}{G_{jln}} = \braket{G_i}{G_j} \braket{G_k}{G_l} \braket{G_m}{G_n} \thinspace , \end{equation}\] where each of the factors is an integral over Gaussian overlap distribution. For example, we can work out the following one-dimensional overlap: \[\begin{align} \braket{G_i}{G_j} &= \int_{-\infty}^{+\infty} \dd{x} \Omega_{ij}(x) \\ &= \sum_t^{i+j} E^{ij}_t \int_{-\infty}^{+\infty} \dd{x} \Lambda_t(x_P) \\ &= E^{ij}_0 \sqrt{\frac{\pi}{p}} \thinspace , \end{align}\] where \(p\) is the total exponent of the overlap distribution and \(\vb{P}\) is its center of mass. The total overlap integral thus becomes: \[\begin{equation} \braket{G_{ikm}}{G_{jln}} = E^{ij}_0 E^{kl}_0 E^{mn}_0 \qty( \frac{\pi}{p} )^{3/2} \thinspace . \end{equation}\]

Linear momentum integrals

For the \(x\)-component of the linear momentum operator \(\vb{p}^c\), we find: \[\begin{equation} -i \matrixel{G_{ikm}}{\qty(\pdv{x})}{G_{jln}} = -i \matrixel{G_i}{\qty(\pdv{x})}{G_j} S_{kl} S_{mn} \thinspace , \end{equation}\] where the one-dimensional linear momentum integral is given by: \[\begin{equation} -i \matrixel{G_i}{\qty(\pdv{x})}{G_j} = 2 i \beta S_{i, j+1} - i j S_{i,j-1} \thinspace , \end{equation}\] where the prefactor \(i\) is the complex unit and should not be confused with the Cartesian exponent in \(G_i\).

Kinetic energy integrals

The kinetic energy integrals may be factored. We find: \[\begin{equation} \matrixel{G_{ikm}}{T^c}{G_{jln}} = T_{ij} S_{kl} S_{mn} + S_{ij} T_{kl} S_{mn} + S_{ij} S_{kl} T_{mn} \thinspace , \end{equation}\] where \(T_{ij}\) is a one-dimensional kinetic energy integral: \[\begin{equation} T_{ij} = -\frac{1}{2} \matrixel{G_i}{ \qty(\pdv[2]{x}) }{G_j} = -2 \beta^2 S_{i,j+2} + \beta (2j + 1) S_{ij} - \frac{1}{2} j(j-1) S_{i,j-2} \thinspace . \end{equation}\]

Dipole/position integrals

The dipole integrals may also be factored. If we want to calculate the \(x\)-component of the dipole integral, with respect to the dipole origin \(\vb{O}\), we calculate: \[\begin{equation} - \matrixel{G_{ikm}}{x_O}{G_{jln}} = - \matrixel{G_i}{x_O}{G_j} S_{kl} S_{mn} \thinspace , \end{equation}\] with the position integral: \[\begin{equation} \matrixel{G_i}{x_O}{G_j} = \sqrt{\frac{\pi}{p}} \qty( E^{ij}_1 + \Delta_{x, PO} E^{ij}_0 ) \thinspace . \end{equation}\]

Angular momentum integrals

For the integrals over the \(x\)-component of the angular momentum operator \(\vb{L}^c(\vb{O})\) about the reference point \(\vb{O}\), we find: \[\begin{equation} -i \matrixel{G_{ikm}}{ \qty[\vb{r}_O \cross \grad]_x }{G_{jln}} = -i S_{ij} \qty[ \matrixel{G_k}{y_O}{G_l} \matrixel{G_m}{ \qty(\pdv{z}) }{G_n} - \matrixel{G_k}{ \qty(\pdv{y}) }{G_l} \matrixel{G_m}{z_O}{G_n} ] \thinspace , \end{equation}\] so we can see that these angular momentum integrals are combinations of dipole/position integrals and linear momentum integrals. The other integrals (over the other components) can be found analogously by cyclic permutation.