Møller–Plesset

Møller–Plesset#

Møller–Plesset partitioning#

Using Rayleigh–Schrödinger perturbation theory (RSPT), the Hamiltonian can be partitioned according to the Møller–Plesset (MP) scheme [HJO14, SB09, SO12]

\hat{H} = \hat{F} + \hat{Φ}

where the zeroth-order Hamiltonian, $\hat{F}$ , is the Fock operator and the perturbation, $\hat{Φ}$ , is known as the fluctuation potential.

The exact solutions to the zeroth-order Hamiltonian are found from the time-independent Schrödinger equation

\begin{aligned} \hat{F} | Φ_{n} ⟩ & = E_{n} | Φ_{n} ⟩ \end{aligned}

The exact state solutions to this equation are given by the Hartree–Fock state and its excited determinants

\begin{aligned} | Φ_{0} ⟩ & = | Ψ_{HF} ⟩ \\ | Φ_{n} ⟩ & = | Ψ_{i j \dots}^{a b \dots} ⟩; n > 0 \end{aligned}

and, in terms of orbital energies, the associated exact energies are equal to

\begin{aligned} E_{0} & = \sum_{i}^{occ} ε_{i} \\ E_{i j \dots}^{a b \dots} & = ε_{a} + ε_{b} + \dots - ε_{i} - ε_{j} - \dots \end{aligned}

Perturbation expansion in the fluctuation potential#

The exact ground state is a solution to

\hat{H} | Ψ ⟩ = E | Ψ ⟩

where $| Ψ ⟩$ and $E$ are expanded in orders of the perturbation

| Ψ ⟩ = \sum_{k = 0}^{\infty} | Ψ^{(k)} ⟩; E = \sum_{k = 0}^{\infty} E^{(k)}

From RSPT, we get

\begin{aligned} E^{(0)} & = E_{0} \\ E^{(1)} & = ⟨ Ψ_{HF} | \hat{Φ} | Ψ_{HF} ⟩ = - \frac{1}{2} \sum_{i, j}^{occ} ⟨ i j ‖ i j ⟩ \end{aligned}

so that

E_{HF} = E^{(0)} + E^{(1)}

The first-order correction to the wave function is obtained from the RSPT master equation

| Ψ^{(1)} ⟩ = - \frac{1}{4} \sum_{i, j}^{occ} \sum_{a, b}^{virt} t_{i j}^{a b} | Ψ_{i j}^{a b} ⟩; t_{i j}^{a b} = \frac{⟨ a b ‖ i j ⟩}{E_{i j}^{a b}}

and, in terms of the introduced $t$ -amplitudes, we thereby get a second-order energy correction equal to

E^{(2)} = - \frac{1}{4} \sum_{i, j}^{occ} \sum_{a, b}^{virt} t_{i j}^{a b} ⟨ i j | | a b ⟩

such that

E_{MP 2} = E_{HF} + E^{(2)}

Closed-shell reference state#

For a closed-shell Hartree–Fock state, the MP2 energy correction can be written in terms of sums referring to molecular orbitals (MOs)

E^{(2)} = E_{os}^{(2)} + E_{ss}^{(2)}

where

\begin{aligned} E_{os}^{(2)} & = - \sum_{i, j}^{occ} \sum_{a, b}^{virt} \frac{(i a | j b) (i a | j b)}{E_{i j}^{a b}} \\ E_{ss}^{(2)} & = - \sum_{i, j}^{occ} \sum_{a, b}^{virt} \frac{(i a | j b)}{E_{i j}^{a b}} [(i a | j b) - (i b | j a)] \end{aligned}

The chemist’s notation has here been used for the integrals.

Illustration#

Let us illustrate the theory with a calculation of the MP2 energy of water.

import veloxchem as vlx

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, "cc-pvdz", ostream=None)

nocc = molecule.number_of_alpha_electrons()
nvirt = basis.get_dimension_of_basis() - nocc

print("Number of occupied orbitals:", nocc)
print("Number of unoccupied orbitals:", nvirt)

Number of occupied orbitals: 5
Number of unoccupied orbitals: 19

molecule.show()

Hartree–Fock reference state#

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

We obtain the orbital energies.

E = scf_results["E_alpha"]

Integrals in the MO basis#

We create an instance of the MOIntegralsDriver class to get the electron repulsion integrals in the MO basis and in the chemist’s notation.

erimo_drv = vlx.MOIntegralsDriver()

ovov = erimo_drv.compute_in_memory(molecule, basis, scf_drv.mol_orbs, "chem_ovov")
print("(ov|ov):", ovov.shape)

(ov|ov): (5, 19, 5, 19)

The MP2 energy correction#

We now have all the ingredients to compute the opposite-spin and same-spin components of the MP2 energy correction.

e_mp2_ss = 0.0
e_mp2_os = 0.0

# extract the occupied subset of the orbital energies
e_o = E[:nocc]
# extract the virtual subset of the orbital energies
e_v = E[nocc:]

for i in range(nocc):
    for j in range(nocc):
        for a in range(nvirt):
            for b in range(nvirt):
                # enegy denominators
                e_ijab = e_v[a] + e_v[b] - e_o[i] - e_o[j]

                # update opposite-spin component of the energy
                e_mp2_os -= (ovov[i, a, j, b] * ovov[i, a, j, b]) / e_ijab

                # update same-spin component of the energy
                e_mp2_ss -= (
                    ovov[i, a, j, b] * (ovov[i, a, j, b] - ovov[i, b, j, a]) / e_ijab
                )

print(f"Opposite-spin MP2 energy: {e_mp2_os:12.8f}")
print(f"Same-spin MP2 energy    : {e_mp2_ss:12.8f}")
print(f"MP2 energy              : {e_mp2_os + e_mp2_ss:12.8f}")

Opposite-spin MP2 energy:  -0.15163083
Same-spin MP2 energy    :  -0.05138187
MP2 energy              :  -0.20301270

Let us compare this result with the MP2 energy obtained with VeloxChem.

mp2_drv = vlx.Mp2Driver()
mp2_drv.ostream.mute()

mp2_results = mp2_drv.compute_conventional(molecule, basis, scf_drv.mol_orbs)

print(f"Energy difference: {mp2_results['mp2_energy']:12.8f}")

Energy difference:  -0.20301270

We note that the two results are in perfect agreement.

Size consistency#

We see that to second order, Møller-Plesset perturbation theory only involves up to double excitations from the HF reference. It would thus be natural to consider it an approximation to CISD, and expect it to suffer from the same issue, namely a lack of size consistency. Yet this is not the case as discussed in the RSPT section on this topic.

dimer_xyz = """6
2 water 100 Å apart                                                                                                            
O    0.000000000000        0.000000000000        0.000000000000                         
H    0.000000000000        0.740848095288        0.582094932012                         
H    0.000000000000       -0.740848095288        0.582094932012
O  100.000000000000        0.000000000000        0.000000000000                         
H  100.000000000000        0.740848095288        0.582094932012                         
H  100.000000000000       -0.740848095288        0.582094932012
"""

molecule = vlx.Molecule.read_xyz_string(dimer_xyz)
basis = vlx.MolecularBasis.read(molecule, "cc-pvdz", ostream=None)

molecule.show()

scf_results = scf_drv.compute(molecule, basis)

mp2_dimer_results = mp2_drv.compute_conventional(molecule, basis, scf_drv.mol_orbs)

print("=" * 30)
print("MP2 energies")
print("-" * 30)
print(f"Dimer       : {mp2_dimer_results['mp2_energy']:16.10f}")
print(f"Two monomers: {2 * mp2_results['mp2_energy']:16.10f}")
print("=" * 30)

==============================
MP2 energies
------------------------------
Dimer       :    -0.4060254140
Two monomers:    -0.4060254051
==============================

We can see that the two energies match. MP2 is size consistent! So why is it behaving better than CISD in this aspect?

The key is that in MP2, the coefficients of the excited determinants are independent of the system size. Thus a molecule would have the same MP2 energy correction regardless of the presence or not of another, non-interacting, molecule. By contrast, in CISD, the coefficients depend on the system size through normalization, lowering the weight of these determinants as the size of the system increases.

Taylor expansion of the energy#

Using RSPT, it is a common practice to use an order parameter $λ$ for the perturbation

\hat{H} (λ) = \hat{H_{0}} + λ \hat{V}

Each order of the perturbation then corresponds to a power of $λ$

E (λ) = E^{(0)} + λ E^{(1)} + λ^{2} E^{(2)} + \dots

The RSPT equations are found by collecting terms in powers of $λ$ and, at the end, set $λ$ 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_{FCI} (λ) = \sum_{i = 0}^{\infty} λ^{i} E_{MP}^{(i)}

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

E_{MP}^{(n)} = \frac{1}{n!} {\frac{d^{n} E_{FCI} (λ)}{d λ^{n}} |}_{λ = 0}

The total MP energy to order $n$ is a sum of the corrections

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

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

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

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})_{λ = 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()

<frozen importlib._bootstrap>:488: RuntimeWarning: mpi4py.MPI.Session size changed, may indicate binary incompatibility. Expected 32 from C header, got 40 from PyObject

Determine $E_{FCI} (λ)$ 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_MP2 = 0.5 * (E_FCI[0] + E_FCI[2] - 2 * E_FCI[1]) / 0.01**2

print(f"MP2 energy correction:{E_MP2:12.8f}")

MP2 energy correction: -0.12747116

Perform a reference calculation using the MP2 module in VeloxChem.

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.