Other partitionings

Taylor expansion of the energy

Using RSPT, it is a common practice to use an order parameter \(\lambda\) for the perturbation

\[ \hat{H}(\lambda) = \hat{H_0} + \lambda \hat{V} \]

Each order of the perturbation then corresponds to a power of \(\lambda\)

\[ E(\lambda) = E^{(0)} + \lambda E^{(1)} + \lambda^2 E^{(2)} + \ldots \]

The RSPT equations are found by collecting terms in powers of \(\lambda\) and, at the end, set \(\lambda\) equal to 1 to recover the correct Hamiltonian. Taking this approach, it is seen that RSPT relates to Taylor expansions. We will here demonstrate that

\[ E_\mathrm{FCI}(\lambda) = \sum_{i=0}^\infty \lambda^i E_{\mathrm{MP}}^{(i)} \]

or, in other words, the \(n\)th-order MP energy correction is given by

\[ E_{\mathrm{MP}}^{(n)} = \frac{1}{n!} \left. \frac{d^n E_\mathrm{FCI}(\lambda)}{d \lambda^n} \right|_{\lambda=0} \]

The total MP energy to order \(n\) is a sum of the corrections

\[ E_{\mathrm{MPn}} = \sum_{i=0}^n E_{\mathrm{MP}}^{(i)} \]

We will explicitly create \(\hat{H}(\lambda)\) and determine the FCI energy. This allows us to determine the MP2 energy by means of numerical differentiation to second order with respect to \(\lambda\).

First, we determine the Hartree–Fock state of water in the 6-31G basis set.

import veloxchem as vlx
import numpy as np

h2o_xyz = """3
water                                                                                                                          
O    0.000000000000        0.000000000000        0.000000000000                         
H    0.000000000000        0.740848095288        0.582094932012                         
H    0.000000000000       -0.740848095288        0.582094932012
"""

molecule = vlx.Molecule.read_xyz_string(h2o_xyz)
basis = vlx.MolecularBasis.read(molecule, "6-31G", ostream=None)

scf_drv = vlx.ScfRestrictedDriver()
scf_drv.ostream.mute()
scf_results = scf_drv.compute(molecule, basis)

# MO coefficients
C = scf_results["C_alpha"]

# orbital energies
E = scf_results["E_alpha"]

Second, we create \(\hat{H_0} = \hat{F}\) and \(\hat{V} = (\hat{H} - \hat{F})_{\lambda=1}\) in the MO basis. The perturbation has a one-electron as well as a two-electron part that we will keep separate

\[ \hat{V} = \hat{V}_1 + \hat{V}_2 \]

import numpy as np

F = np.diag(E)

# Get the core Hamiltonian in MO basis
T_ao = vlx.compute_kinetic_energy_integrals(molecule, basis)

V_nucpot_ao = vlx.compute_nuclear_potential_integrals(molecule, basis)

h = np.einsum("ai, ab, bj -> ij", C, T_ao + V_nucpot_ao, C)

# Compute the 2-electron integrals
fock_drv = vlx.FockDriver()
g_ao = fock_drv.compute_eri(molecule, basis)

g_1 = np.einsum("ds, abcd -> abcs", C, g_ao)
g_2 = np.einsum("cr, abcs -> abrs", C, g_1)
g_3 = np.einsum("bq, abrs -> aqrs", C, g_2)

g = np.einsum("ap, aqrs -> pqrs", C, g_3)

H0 = F

V_1 = h - F
V_2 = g

Third, we use MultiPsi to determine the FCI energies.

import multipsi as mtp

ci_drv = mtp.CIDriver()
ci_drv.ostream.mute()

space = mtp.OrbSpace(molecule, scf_drv.mol_orbs)
space.fci()

Determine \(E_\mathrm{FCI}(\lambda)\) and calculate the second derivative with use of a 3-point stencil.

E_FCI = []

for lambda_val in [-0.01, 0, 0.01]:

    ci_drv._update_integrals(
        molecule.nuclear_repulsion_energy(), H0 + lambda_val * V_1, lambda_val * V_2
    )

    ci_results = ci_drv.compute(molecule, basis, space, 1)

    E_FCI.append(ci_results["energies"][0])

E_Taylor2 = 0.5 * (E_FCI[0] + E_FCI[2] - 2 * E_FCI[1]) / 0.01**2

print(f"Numerical 2nd order Taylor energy correction:{E_Taylor2:12.8f}")

Numerical 2nd order Taylor energy correction: -0.12747116

Perform a reference calculation using the MP2 module in VeloxChem.

mp2_drv = vlx.Mp2Driver()
mp2_results = mp2_drv.compute(molecule, basis, scf_drv.mol_orbs)

print(f"MP2 energy correction:{mp2_results['mp2_energy']:12.8f}")

MP2 energy correction: -0.12747066

We note the the two results are in perfect agreement. The small discrepancy is due to the use of numerical differentiation in the former case.

Convergence of the perturbation series

While Møller–Plesset is often used to second order, it forms a series that one would expect to converge to the full CI energy. However, like any Taylor expansion, there is a convergence radius. This means if the perturbation is too large, the series may not converge to the desired outcome, even diverging. This can easily be illustrated on a two electrons in two orbital basis like H\(_2\) in a minimal basis, for which we can compute the exact result easily.

In such a system, there are in principle four determinants, but symmetries reduce this number to two relevant ones for the ground state.

\[ | 0 \rangle = | \sigma_g \bar{\sigma_g} \rangle \]

\[ | 1 \rangle = | \sigma_u \bar{\sigma_u} \rangle \]

with energies \(\langle 0 | \hat{H} | 0 \rangle = E_0 \) and \(\langle 1 | \hat{H} | 1 \rangle = E_1 \) and interaction term \(\langle 0 | \hat{H} | 1 \rangle = V_{01}\).

The Hamiltonian is thus

\[\begin{split} \begin{pmatrix} E_0 & V_{01}\\ V_{01} & E_1 \end{pmatrix} \end{split}\]

or if we shift everything by the HF ground state energy \(E_0\)

\[\begin{split} \begin{pmatrix} 0 & V_{01}\\ V_{01} & E_1 - E_0 \end{pmatrix} \end{split}\]

The eigenvalues of such a 2x2 matrix are known, as it corresponds to the roots of a second order polynomial:

\[ E_\mathrm{corr} = \frac{\Delta E \pm \sqrt{\Delta E^2 + 4V_{01}^2}}{2} \]

with \(\Delta E = E_1 - E_0 \geq 0\). Assuming \(\Delta E \neq 0\) we can write this as

\[ E_\mathrm{corr} = \frac{\Delta E \pm \Delta E \sqrt{1 + 4V_{01}^2/\Delta E^2}}{2} \]

Since we want the ground state, we keep only the solution with a minus sign. Now we can use the known Taylor expansion \(\sqrt{1+x} = 1 - \frac{x}{2} + o(x) \)

\[ E_\mathrm{corr} \approx - \frac{V_{01}}{\Delta E}\]

We recognize here the second order perturbation theory result. Specifically, here, this is not the Epstein-Nesbet perturbation theory result, but similar reasoning could be done using Møller–Plesset.

The important thing is that the Taylor series \(\sqrt{1+x} = 1 - \frac{x}{2} + o(x) \) is only convergent for \(|x| \leq 1\), which in our case means for \(2|V_{01}| \leq \Delta E \). In other words, if the energy difference between the ground and excited determinants become too small, the Taylor series -and thus the perturbation expansion- would not converge.

This regime corresponds to the strong correlation regime and requires more sophisticated electronic structure methods.