Møller-Plesset perturbation theory to second order
Contents
Møller-Plesset perturbation theory to second order¶
Theory refresher¶
Formal Rayleigh-Schrödinger perturbation theory¶
We want to solve the time-independent Schrödinger equation:
In perturbation theory, the Hamiltonian is partitioned into zeroth-order and perturbation terms:
with \(\lambda\) the coupling strength of the perturbation. We assume that the complete spectrum of the zeroth-order Hamiltonian is known:
The energy and wavefunction are expanded in a formal power series in \(\lambda\) and terms to the same order are collected on the left- and right-hand sides:
Thus:
Energy corrections and perturbative expansion coefficients, \(a_{kn}^{(m)} \equiv \langle \Psi_{k}^{(0)} | \Psi_{n}^{(m)} \rangle\) , for the wavefunction are obtained by projection onto the set of known zeroth-order eigenfunctions:
and:
Møller-Plesset partitioning and second order energy¶
In general, the partitioning of the Hamiltonian can be achieved in a number of ways. In molecular electronic structure theory, \(H\) is the Born-Oppenheimer, many-body, molecular electronic Hamiltonian. The Møller-Plesset partitioning is probably the most common, since it follows quite naturally when first approximating the desired ground state in terms of a single reference determinant. The Hamiltonian is in fact rewritten as:
where \(F\) is the Fock operator and \(\Phi\) is the so-called fluctuation potential. \(F\) is a one-body operator whose spectrum consists of all determinants that can be built from excitation of the reference \(| 0 \rangle\):
where the zeroth-order energy and the orbital energy denominators are:
The fluctuation potential is a two-body operator. With this partinioning, the zeroth-order energy is the sum of orbital energies. The first-order correction is:
that is, the energy of the reference single determinant is correct throught first order in the perturbative series: \(E_{\mathrm{ref}} = E_{0}^{(0)} + E_0^{(1)}\).
The first-order wavefunction is obtained from the general RSPT expression. In a basis of molecular spin-orbitals:
where the orbital energy denominator is: \(\varepsilon_{ij}^{ab} = \varepsilon_{i} + \varepsilon_{j} -\varepsilon_{a} - \varepsilon_{b} \). The second order energy correction follows:
For a closed-shell, restricted reference using real MOs:
which we can further rearrange into to two terms, opposite-spin and same-spin:
Implementation¶
To compute \(E_{\mathrm{MP2}}\) we need to:
Obtain the reference closed-shell determinant from a Hartree-Fock calculation.
Transform the AO basis ERI tensor to MO basis.
Assemble the energy denominators.
Combine the results of steps 2 and 3 to form the perturbative correction.
We start with the declaration of the usual water molecule and its basis set. We also perform the SCF calculation with the ScfRestrictedDriver.
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
"""
mol = vlx.Molecule.from_xyz_string(h2o_xyz)
basis = vlx.MolecularBasis.read(mol, "cc-pvdz")
scfdrv = vlx.ScfRestrictedDriver()
scfdrv.compute(mol, basis)
* Warning * Environment variable OMP_NUM_THREADS not set.
* Warning * Setting OMP_NUM_THREADS to 8.
Self Consistent Field Driver Setup
====================================
Wave Function Model : Spin-Restricted Hartree-Fock
Initial Guess Model : Superposition of Atomic Densities
Convergence Accelerator : Two Level Direct Inversion of Iterative Subspace
Max. Number of Iterations : 50
Max. Number of Error Vectors : 10
Convergence Threshold : 1.0e-06
ERI Screening Scheme : Cauchy Schwarz + Density
ERI Screening Mode : Dynamic
ERI Screening Threshold : 1.0e-12
Linear Dependence Threshold : 1.0e-06
* Info * Nuclear repulsion energy: 9.3436381580 au
* Info * Overlap matrix computed in 0.05 sec.
* Info * Kinetic energy matrix computed in 0.03 sec.
* Info * Nuclear potential matrix computed in 0.01 sec.
* Info * Orthogonalization matrix computed in 0.10 sec.
* Info * SAD initial guess computed in 0.06 sec.
* Info * Starting Reduced Basis SCF calculation...
* Info * ...done. SCF energy in reduced basis set: -75.979046359568 au. Time: 1.66 sec.
* Info * Overlap matrix computed in 0.08 sec.
* Info * Kinetic energy matrix computed in 0.00 sec.
* Info * Nuclear potential matrix computed in 0.01 sec.
* Info * Orthogonalization matrix computed in 0.00 sec.
Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change
--------------------------------------------------------------------------------------------
1 -76.025869744699 0.0000000000 0.09892066 0.01121867 0.00000000
2 -76.026932827150 -0.0010630825 0.01920269 0.00294570 0.02594960
3 -76.026981056596 -0.0000482294 0.00412658 0.00076652 0.00507499
4 -76.026983601528 -0.0000025449 0.00247905 0.00043111 0.00166621
5 -76.026984175646 -0.0000005741 0.00021766 0.00003639 0.00052311
6 -76.026984187023 -0.0000000114 0.00003206 0.00000540 0.00012084
7 -76.026984187254 -0.0000000002 0.00000208 0.00000029 0.00001663
8 -76.026984187255 -0.0000000000 0.00000053 0.00000008 0.00000087
*** SCF converged in 8 iterations. Time: 4.48 sec.
Spin-Restricted Hartree-Fock:
-----------------------------
Total Energy : -76.0269841873 au
Electronic Energy : -85.3706223452 au
Nuclear Repulsion Energy : 9.3436381580 au
------------------------------------
Gradient Norm : 0.0000005316 au
Ground State Information
------------------------
Charge of Molecule : 0.0
Multiplicity (2S+1) : 1.0
Magnetic Quantum Number (M_S) : 0.0
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 1:
--------------------------
Occupation: 2.0 Energy: -20.54819 au
( 1 O 1s : 1.00)
Molecular Orbital No. 2:
--------------------------
Occupation: 2.0 Energy: -1.34520 au
( 1 O 2s : -0.44) ( 1 O 3s : -0.36) ( 2 H 1s : -0.20)
( 3 H 1s : -0.20)
Molecular Orbital No. 3:
--------------------------
Occupation: 2.0 Energy: -0.70585 au
( 1 O 1p-1: -0.49) ( 1 O 2p-1: -0.22) ( 2 H 1s : -0.33)
( 3 H 1s : 0.33)
Molecular Orbital No. 4:
--------------------------
Occupation: 2.0 Energy: -0.57109 au
( 1 O 2s : 0.15) ( 1 O 3s : 0.36) ( 1 O 1p0 : -0.55)
( 1 O 2p0 : -0.36) ( 2 H 1s : -0.21) ( 3 H 1s : -0.21)
Molecular Orbital No. 5:
--------------------------
Occupation: 2.0 Energy: -0.49457 au
( 1 O 1p+1: -0.63) ( 1 O 2p+1: -0.49)
Molecular Orbital No. 6:
--------------------------
Occupation: 0.0 Energy: 0.18787 au
( 1 O 3s : 1.02) ( 1 O 1p0 : 0.19) ( 1 O 2p0 : 0.33)
( 2 H 2s : -0.83) ( 3 H 2s : -0.83)
Molecular Orbital No. 7:
--------------------------
Occupation: 0.0 Energy: 0.25852 au
( 1 O 1p-1: -0.28) ( 1 O 2p-1: -0.67) ( 2 H 2s : 1.48)
( 3 H 2s : -1.48)
Molecular Orbital No. 8:
--------------------------
Occupation: 0.0 Energy: 0.79749 au
( 1 O 1p-1: 0.26) ( 1 O 2p-1: 0.49) ( 2 H 1s : -0.95)
( 2 H 2s : 0.66) ( 2 H 1p0 : -0.16) ( 3 H 1s : 0.95)
( 3 H 2s : -0.66) ( 3 H 1p0 : 0.16)
Molecular Orbital No. 9:
--------------------------
Occupation: 0.0 Energy: 0.87271 au
( 1 O 2s : 0.25) ( 1 O 3s : -0.34) ( 1 O 1p0 : 0.34)
( 2 H 1s : -0.77) ( 2 H 2s : 0.54) ( 2 H 1p-1: -0.31)
( 3 H 1s : -0.77) ( 3 H 2s : 0.54) ( 3 H 1p-1: 0.31)
Molecular Orbital No. 10:
--------------------------
Occupation: 0.0 Energy: 1.16315 au
( 1 O 3s : -0.78) ( 1 O 1p0 : 0.74) ( 1 O 2p0 : -1.31)
( 2 H 1s : 0.59) ( 2 H 1p0 : -0.24) ( 3 H 1s : 0.59)
( 3 H 1p0 : -0.24)
We can now access orbital energies and MO coefficients from the driver:
epsilon = scfdrv.scf_tensors["E"]
C = scfdrv.scf_tensors["C"]
The Integral transformation¶
We compute the MP2 energy correction with the ERI expressed in MO basis: we need to transform the ERI tensor from AO basis. The transformation reads:
with the MO integrals in physicists’ notation. The transformation requires \(O(N^{8})\) operation count. However, we can perform it more efficiently as a stepwise contraction:
We should also note that we do not need the full ERI tensor in MO basis, but rather the OOVV class of integrals, which involve two occupied and two virtual MO indices:
eridrv = vlx.ElectronRepulsionIntegralsDriver()
mknl = eridrv.compute_in_mem(mol, basis)
import numpy as np
N_O = mol.number_of_electrons() // 2
N_V = scfdrv.mol_orbs.number_mos() - N_O
# lambda -> b transformation
mknb = np.einsum("mknl,lB->mknB", mknl, C[:, N_O:])
print(f"{mknb.shape=}")
# TODO kappa -> a transformation
mnab =
print(f"{mnab.shape=}")
# TODO nu -> j transformation
mjab =
print(f"{mjab.shape=}")
# TODO mu -> i transformation
ijab =
print(f"{ijab.shape=}")
Let’s compare our OOVV ERI tensor with the one computed by using VeloxChem’s own MOIntegralsDriver:
moeridrv = vlx.MOIntegralsDriver()
moeri = moeridrv.compute_in_mem(mol, basis, mol_orbs=scfdrv.mol_orbs, mints_type="OOVV")
np.testing.assert_allclose(ijab, moeri, atol=1.e-10)
The MP2 energy correction¶
We now have all the ingrediente to compute the opposite-spin and same-spin components of the MP2 energy correction:
e_mp2_ss = 0.0
e_mp2_os = 0.0
# TODO extract the occupied subset of the orbital energies using NumPy slicing
e_ij =
# TODO extract the virtual subset of the orbital energies using NumPy slicing
e_ab =
# loop over occupied orbitals
for i in range(N_O):
# loop over occupied orbitals
for j in range(N_O):
# loop over virtual orbitals
for a in range(N_V):
# loop over virtual orbitals
for b in range(N_V):
# TODO form enegy denominator for given i, j, a, b quadruplet of indices
e_ijab =
# TODO update opposite-spin component of the energy
e_mp2_os -=
# TODO update same-spin component of the energy
e_mp2_ss -=
print(f"Opposite-spin MP2 energy: {e_mp2_os:20.12f}")
print(f"Same-spin MP2 energy: {e_mp2_ss:20.12f}")
print(f"MP2 energy: {e_mp2_os + e_mp2_ss:20.12f}")
VeloxChem has its own implementation of the MP2 energy correction. We can check our result against it.
mp2drv = vlx.Mp2Driver()
mp2drv.compute_conventional(mol, basis, scfdrv.mol_orbs)
np.testing.assert_allclose(e_mp2_os + e_mp2_ss, mp2drv.e_mp2, atol=1e-9)
An alternative implmentation using np.einsum¶
We can compute the MP2 energy terms using np.einsum to handle the tensor contractions. First, we need to define the energy denominators as a 4-index tensor with np.reshape.
e_ijab = 1 / (e_ab.reshape(1, 1, N_V, 1) + e_ab.reshape(1, 1, 1, N_V) - e_ij.reshape(N_O, 1, 1, 1) - e_ij.reshape(1, N_O, 1, 1))
The same-spin and opposite-spin are the contraction of the ERI and denominators tensors:
mp2_os = - np.einsum("ijab,ijab,ijab->", ijab, ijab, e_ijab)
mp2_ss = - np.einsum("ijab,ijab,ijab->", ijab, ijab - ijab.swapaxes(2, 3), e_ijab)
The results are unchanged with respect to the quadruple-loop approach:
print(f"Opposite-spin MP2 energy: {mp2_os:20.12f}")
print(f"Same-spin MP2 energy: {mp2_ss:20.12f}")
print(f"MP2 energy: {mp2_os + mp2_ss:20.12f}")
# compare with MP2 energy components computed with a quadruple loop
np.testing.assert_allclose(mp2_os, e_mp2_os, atol=1e-9)
np.testing.assert_allclose(mp2_ss, e_mp2_ss, atol=1e-9)
# compare with MP2 energy computed by VeloxChem
np.testing.assert_allclose(mp2_os + mp2_ss, mp2drv.e_mp2, atol=1e-9)
Bonus: the MP2 one-particle density matrix¶
The computation of expectation values of one-electron observables, such as the electric dipole moment, requires the one-particle density matrix (OPDM). In restricted Hartree-Fock theory, the MO basis OPDM is quite simple:
Let us compute the electronic component of the electric dipole moment. First, we obtain the necessary integrals in AO basis and transform them to MO basis. Note that the dipole moment has a 3 Cartesian components: the ElectricDipoleIntegralsDriver will compute all of them and the to_numpy method thus returns a list of AO basis integral matrices.
mudrv = vlx.ElectricDipoleIntegralsDriver()
mu_ao = mudrv.compute(mol, basis).to_numpy()
# TODO obtain the MO basis electric dipole moment integrals
mu_mo =
Next, we construct the RHF density matrix and compute the dipole moment.
# TODO form the OPDM of the RHF reference determinant
Dpq =
# TODO compute the RHF electronic component of the electric dipole moment
mu_rhf =
print(mu_rhf)
Lagrangian, amplitudes, and multipliers¶
Computing the MP2 correction to the density matrix requires not only the MP1 amplitudes (wavefunction parameters), but also the MP2 multipliers. The MP2 Lagrangian is:
At its stationary points, the Lagrangian is the MP2 energy. The amplitudes and multipliers that make the Lagrangian stationary are the MP1 amplitudes, \(t_{ij}^{ab\,(1)}\), and multipliers, \(\lambda^{ij\,(1)}_{ab}\):
For a closed-shell system, the amplitudes and multipliers are explicitly given as:
The unrelaxed MP2 density matrix¶
The expectation value of a one-electron operator \(X\) is then:
where we identifyied, after some algebra, the MP2 density matrix as:
Note that:
These expressions are not manifestly symmetric. In practice, we will symmetrize the MP2 density matrix before computing expectation values.
This is the so-called unrelaxed density matrix. We are using fixed molecular orbitals, unaffected by electron correlation.
Let us use np.einsum to form all these quantities and compute the unrelaxed MP2 electronic component of the electric dipole moment.
# TODO form the doubles amplitudes, i.e. the MP1 wavefunction coefficients, from the ERI and denominator tensors
t2 =
# TODO form the doubles multiplitiers, from the ERI and denominator tensors
l2 =
We now form the OO and VV blocks of the MP2 OPDM, following the equations above, and update the reference MO density matrix with these contributions:
# TODO form the OO block of the MP2 OPDM, don't forget to symmetrize it!
Dij =
# TODO form VV block of the MP2 OPDM, don't forget to symmetrize it!
Dab =
# add the MP2 contributions to the OO and VV blocks
Dpq[:N_O, :N_O] += Dij
Dpq[N_O:, N_O:] += Dab
# test that the total density matrix is symmetric
np.testing.assert_allclose(Dpq, np.transpose(Dpq))
Finally, the dipole moment:
# TODO compute the MP2 unrelaxed electronic component of the electric dipole moment
mu_mp2 =
print(mu_mp2)
Natural occupations and natural orbitals¶
The eigenvectors and eigenvalues of the one-particle density matrix provide the natural orbitals and the natural-orbital occupation numbers. These are useful, for example, in setting up the active space in multireference calculations.
# TODO diagonalize the density matrix
omega, NOs =
# eigenvalues and eigenvectors are in *ascending* order, so we resort them in *descending* order
idx = omega.argsort()[::-1]
omega = omega[idx]
NOs = NOs[:,idx]
One way to check that our MP2 unrelaxed OPDM is correct is to sum the natural occupation numbers: this should give us the number of electrons in the system.
print(sum(omega))
We can also double-check against reference values for the occupations from another code.
# reference natural orbital occupation numbers from DALTON
ref_omega = np.array([
1.99990540, 1.98720752, 1.97426785, 1.97108868, 1.96924405,
0.02241866, 0.02020351, 0.01713431, 0.01024357, 0.00551830,
0.00517755, 0.00472951, 0.00414944, 0.00404548, 0.00090056,
0.00086293, 0.00060545, 0.00051955, 0.00046726, 0.00045319,
0.00039740, 0.00037924, 0.00004262, 0.00003795
])
np.testing.assert_allclose(omega, ref_omega, atol=1e-3)
References¶
Shavitt, I.; Bartlett, R. J. Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory Cambridge Molecular Science; Cambridge University Press, 2009.
Helgaker, T.; Jørgensen, P.; Olsen, J. Molecular Electronic-Structure Theory Wiley, 2000.