Møller–Plesset#

Møller–Plesset partitioning#

In molecular electronic structure theory, \(\hat{H}\) is the Born-Oppenheimer, many-body, molecular electronic Hamiltonian. However, in general, the partitioning of the Hamiltonian, i.e., the definition of \(\hat{H}_0\), can be achieved in a number of ways.

The most common partitioning in quantum chemistry is probably the Møller–Plesset (MP) partitioning [HJO14, SO12], 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:

\[\begin{equation*} \hat{H} = \hat{F} + \hat{\Phi} \end{equation*}\]

where \(\hat{F}\) is the Fock operator and \(\hat{\Phi}\) is the so-called fluctuation potential. \(\hat{F}\) is a one-body operator whose spectrum consists of all determinants that can be built from excitation of the reference \(| 0 \rangle\):

\[\begin{equation*} \hat{F} | 0 \rangle = E_{0}^{(0)}| 0 \rangle, \quad \hat{F} \left|_{ij\ldots} ^{ab\ldots} \right\rangle = \left(E_{0}^{(0)} + \varepsilon^{ab\cdots}_{ij\cdots} \right) \left|_{ij\ldots} ^{ab\ldots} \right\rangle, \end{equation*}\]

where the zeroth-order energy and the orbital-energy denominators are:

\[\begin{equation*} E_{0}^{(0)} = \sum_{i} \varepsilon_{i},\quad \varepsilon^{ab\cdots}_{ij\cdots} = \varepsilon_{a} + \varepsilon_{b} + \ldots - \varepsilon_{i} - \varepsilon_{j} - \ldots \end{equation*}\]

Here and in the following, the indices \(i,j,\ldots\) refer to occupied orbitals, whereas \(a,b,\ldots\) refer to unoccupied (virtual) ones. 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:

\[\begin{equation*} E_0^{(1)} = \left\langle 0 \left| \hat{\Phi} \right| 0 \right\rangle = -\frac{1}{2} \sum_{ij} \langle ij \| ij \rangle \end{equation*}\]

that is, the energy of the reference single determinant is correct through first order in the perturbative series: \(E_{\mathrm{ref}} = E_{\mathrm{HF}} = E_{0}^{(0)} + E_0^{(1)}\).

The first-order correction to the wave function is obtained from the general RSPT expression. In a basis of molecular spin-orbitals, it is given by

\[\begin{equation*} | \Psi_0^{(1)} \rangle = -\frac{1}{4}\sum_{ijab} \frac{\langle ab \| ij \rangle}{\varepsilon_{ij}^{ab}} |_{ij}^{ab}\rangle \, , \end{equation*}\]

where the orbital-energy denominator is \(\varepsilon_{ij}^{ab} = \varepsilon_{a} + \varepsilon_{b} - \varepsilon_{i} - \varepsilon_{j}\), and the expansion coefficients are identified as the so-called \(t\) amplitudes with

\[\begin{equation*} t_{ijab} = \frac{\langle ab \| ij \rangle}{\varepsilon_{ij}^{ab}} \, . \end{equation*}\]

The second-order energy correction follows:

\[\begin{equation*} E_{0}^{(2)} = - \frac{1}{4} \sum_{ijab} \frac{\langle ij \| ab \rangle \langle ab \| ij \rangle}{\varepsilon_{ij}^{ab}} = - \frac{1}{4} \sum_{ijab} t_{ijab} \langle ij || ab \rangle \end{equation*}\]

and the total MP2 energy is given by \(E_{\mathrm{MP2}} = E_{\mathrm{HF}} + E_0^{(2)}\). For a closed-shell, restricted reference using real MOs, \(E_0^{(2)}\) can be written as

\[\begin{equation*} E_0^{(2)} = - \sum_{ij}^{N_{\mathrm{O}}} \sum_{ab}^{N_{\mathrm{V}}} \frac{\langle ij | ab \rangle}{\varepsilon_{ij}^{ab}} [ 2 \langle ij | ab \rangle - \langle ij | ba \rangle ], \end{equation*}\]

where \(N_{\mathrm{O}}\) is the numer of occupied spatial orbitals (corresponding to half the number of electrons), and \(N_{\mathrm{V}}\) the number of virtual spatial orbitals, which we can further rearrange into to two terms, opposite-spin (OS) and same-spin (SS):

\[\begin{equation*} E_0^{(2)} = - \sum_{ij}^{N_{\mathrm{O}}} \sum_{ab}^{N_{\mathrm{V}}} \frac{\langle ij | ab \rangle\langle ij | ab \rangle}{\varepsilon_{ij}^{ab}} - \sum_{ij}^{N_{\mathrm{O}}} \sum_{ab}^{N_{\mathrm{V}}} \frac{\langle ij | ab \rangle[ \langle ij | ab \rangle - \langle ij | ba \rangle ]}{\varepsilon_{ij}^{ab}} = E_{\mathrm{OS}}^{(2)} + E_{\mathrm{SS}}^{(2)}. \end{equation*}\]

Implementation#

To compute \(E_0^{(2)}\) we need to:

  1. Obtain the reference closed-shell determinant from a Hartree-Fock calculation.

  2. Transform the AO basis ERI tensor to MO basis.

  3. Assemble the energy denominators.

  4. Combine the results of steps 2 and 3 to form the perturbative correction.

Obtaining the MP2 energy correction

Obtaining the HF reference#

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.
* Warning * Setting MKL_THREADING_LAYER to "GNU".
                                                                                                                          
                                            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.3436381577 a.u.                                                                      
                                                                                                                          
* Info * Overlap matrix computed in 0.00 sec.                                                                             
                                                                                                                          
* Info * Kinetic energy matrix computed in 0.00 sec.                                                                      
                                                                                                                          
* Info * Nuclear potential matrix computed in 0.00 sec.                                                                   
                                                                                                                          
* Info * Orthogonalization matrix computed in 0.03 sec.                                                                   
                                                                                                                          
* Info * SAD initial guess computed in 0.00 sec.                                                                          
                                                                                                                          
* Info * Starting Reduced Basis SCF calculation...                                                                        
* Info * ...done. SCF energy in reduced basis set: -75.979046359571 a.u. Time: 0.12 sec.                                  
                                                                                                                          
* Info * Overlap matrix computed in 0.00 sec.                                                                             
                                                                                                                          
* Info * Kinetic energy matrix computed in 0.00 sec.                                                                      
                                                                                                                          
* Info * Nuclear potential matrix computed in 0.00 sec.                                                                   
                                                                                                                          
* Info * Orthogonalization matrix computed in 0.00 sec.                                                                   
                                                                                                                          
                                                                                                                          
               Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change               
               --------------------------------------------------------------------------------------------               
                  1       -76.025869744700    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.026983601529   -0.0000025449      0.00247905      0.00043111      0.00166621                
                  5       -76.026984175646   -0.0000005741      0.00021766      0.00003639      0.00052311                
                  6       -76.026984187024   -0.0000000114      0.00003206      0.00000540      0.00012084                
                  7       -76.026984187255   -0.0000000002      0.00000208      0.00000029      0.00001663                
                  8       -76.026984187256   -0.0000000000      0.00000053      0.00000008      0.00000087                
                                                                                                                          
               *** SCF converged in 8 iterations. Time: 0.65 sec.                                                         
                                                                                                                          
               Spin-Restricted Hartree-Fock:                                                                              
               -----------------------------                                                                              
               Total Energy                       :      -76.0269841873 a.u.                                              
               Electronic Energy                  :      -85.3706223449 a.u.                                              
               Nuclear Repulsion Energy           :        9.3436381577 a.u.                                              
               ------------------------------------                                                                       
               Gradient Norm                      :        0.0000005316 a.u.                                              
                                                                                                                          
                                                                                                                          
               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.000 Energy:  -20.54819 a.u.                                                                  
               (   1 O   1s  :    -1.00)                                                                                  
                                                                                                                          
               Molecular Orbital No.   2:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -1.34520 a.u.                                                                  
               (   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.000 Energy:   -0.70585 a.u.                                                                  
               (   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.000 Energy:   -0.57109 a.u.                                                                  
               (   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.000 Energy:   -0.49457 a.u.                                                                  
               (   1 O   1p+1:    -0.63) (   1 O   2p+1:    -0.49)                                                        
                                                                                                                          
               Molecular Orbital No.   6:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.18787 a.u.                                                                  
               (   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.000 Energy:    0.25852 a.u.                                                                  
               (   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.000 Energy:    0.79749 a.u.                                                                  
               (   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.000 Energy:    0.87271 a.u.                                                                  
               (   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.000 Energy:    1.16315 a.u.                                                                  
               (   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"]

Transforming the integrals#

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:

\[\begin{equation*} \langle pq | rs \rangle = \sum_{\mu\nu\kappa\lambda} C_{\mu p}C_{\nu r} (\mu\nu|\kappa\lambda) C_{\kappa q} C_{\lambda s}, \end{equation*}\]

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:

\[\begin{equation*} \langle pq | rs \rangle = \sum_{\mu} C_{\mu p} \left(\sum_{\nu} C_{\nu r} \left (\sum_{\kappa} \left(\sum_{\lambda} (\mu\nu|\kappa\lambda) C_{\lambda s} \right) C_{\kappa q} \right)\right). \end{equation*}\]

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:

\[\begin{equation*} \langle ij | ab \rangle = \sum_{\mu} C_{\mu i} \left(\sum_{\nu} C_{\nu j} \left(\sum_{\kappa} \left(\sum_{\lambda} (\mu\kappa|\nu\lambda) C_{\lambda b} \right) C_{\kappa a}\right)\right). \end{equation*}\]
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

mknb = np.einsum("mknl,lB->mknB", mknl, C[:, N_O:])
print(f"mknb.shape={mknb.shape}")
mnab = np.einsum("mknB,kA->mnAB", mknb, C[:, N_O:])
print(f"mnab.shape={mnab.shape}")
mjab = np.einsum("mnAB,nJ->mJAB", mnab, C[:, :N_O])
print(f"mjab.shape={mjab.shape}")
ijab = np.einsum("mJAB,mI->IJAB", mjab, C[:, :N_O])
print(f"ijab.shape={ijab.shape}")
mknb.shape=(24, 24, 24, 19)
mnab.shape=(24, 24, 19, 19)
mjab.shape=(24, 5, 19, 19)
ijab.shape=(5, 5, 19, 19)

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 ingredients to compute the opposite-spin and same-spin components of the MP2 energy correction:

\[\begin{equation*} E_0^{(2)} = - \sum_{ij}^{N_{\mathrm{O}}} \sum_{ab}^{N_{\mathrm{V}}} \frac{\langle ij | ab \rangle\langle ij | ab \rangle}{\varepsilon_{ij}^{ab}} - \sum_{ij}^{N_{\mathrm{O}}} \sum_{ab}^{N_{\mathrm{V}}} \frac{\langle ij | ab \rangle[ \langle ij | ab \rangle - \langle ij | ba \rangle ]}{\varepsilon_{ij}^{ab}} = E_{\mathrm{OS}}^{(2)} + E_{\mathrm{SS}}^{(2)}. \end{equation*}\]
e_mp2_ss = 0.0
e_mp2_os = 0.0

# extract the occupied subset of the orbital energies
e_ij = epsilon[:N_O]
# extract the virtual subset of the orbital energies
e_ab = epsilon[N_O:]

for i in range(N_O):
    for j in range(N_O):
        for a in range(N_V):
            for b in range(N_V):
                # enegy denominators
                e_ijab = e_ab[a] + e_ab[b] - e_ij[i] - e_ij[j] 
                
                # update opposite-spin component of the energy
                e_mp2_os -= (ijab[i, j, a, b] * ijab[i, j, a, b]) / e_ijab
                
                # update same-spin component of the energy
                e_mp2_ss -= ijab[i, j, a, b] * (ijab[i, j, a, b]  - ijab[i, j, b, a]) / e_ijab
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}")
Opposite-spin MP2 energy:      -0.151630828959
Same-spin MP2 energy:          -0.051381873614
MP2 energy:                    -0.203012702573

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)

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. This can easily be verified numerically on a small case

h2o_2_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
"""

mol = vlx.Molecule.from_xyz_string(h2o_2_xyz)

basis = vlx.MolecularBasis.read(mol, "cc-pvdz")

scfdrv = vlx.ScfRestrictedDriver()
scfdrv.compute(mol, basis)

mp2drv = vlx.Mp2Driver()
mp2drv.compute_conventional(mol, basis, scfdrv.mol_orbs)
                                                                                                                          
                                            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: 19.2164448492 a.u.                                                                     
                                                                                                                          
* Info * Overlap matrix computed in 0.00 sec.                                                                             
                                                                                                                          
* Info * Kinetic energy matrix computed in 0.00 sec.                                                                      
                                                                                                                          
* Info * Nuclear potential matrix computed in 0.00 sec.                                                                   
                                                                                                                          
* Info * Orthogonalization matrix computed in 0.00 sec.                                                                   
                                                                                                                          
* Info * SAD initial guess computed in 0.00 sec.                                                                          
                                                                                                                          
* Info * Starting Reduced Basis SCF calculation...                                                                        
* Info * ...done. SCF energy in reduced basis set: -151.958092575961 a.u. Time: 0.10 sec.                                 
                                                                                                                          
* Info * Overlap matrix computed in 0.00 sec.                                                                             
                                                                                                                          
* Info * Kinetic energy matrix computed in 0.00 sec.                                                                      
                                                                                                                          
* Info * Nuclear potential matrix computed in 0.00 sec.                                                                   
                                                                                                                          
* Info * Orthogonalization matrix computed in 0.00 sec.                                                                   
                                                                                                                          
                                                                                                                          
               Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change               
               --------------------------------------------------------------------------------------------               
                  1      -152.051738140321    0.0000000000      0.13991946      0.01118706      0.00000000                
                  2      -152.053865418939   -0.0021272786      0.02711307      0.00233482      0.03667246                
                  3      -152.053962061620   -0.0000966427      0.00577036      0.00075430      0.00716290                
                  4      -152.053967147969   -0.0000050863      0.00344551      0.00042382      0.00235525                
                  5      -152.053968254474   -0.0000011065      0.00030861      0.00003629      0.00072555                
                  6      -152.053968277233   -0.0000000228      0.00004512      0.00000528      0.00017099                
                  7      -152.053968277688   -0.0000000005      0.00000279      0.00000028      0.00002339                
                  8      -152.053968277689   -0.0000000000      0.00000068      0.00000007      0.00000112                
                                                                                                                          
               *** SCF converged in 8 iterations. Time: 0.72 sec.                                                         
                                                                                                                          
               Spin-Restricted Hartree-Fock:                                                                              
               -----------------------------                                                                              
               Total Energy                       :     -152.0539682777 a.u.                                              
               Electronic Energy                  :     -171.2704131269 a.u.                                              
               Nuclear Repulsion Energy           :       19.2164448492 a.u.                                              
               ------------------------------------                                                                       
               Gradient Norm                      :        0.0000006818 a.u.                                              
                                                                                                                          
                                                                                                                          
               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.   6:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -0.70584 a.u.                                                                  
               (   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.   7:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -0.57109 a.u.                                                                  
               (   4 O   2s  :     0.15) (   4 O   3s  :     0.36) (   4 O   1p0 :    -0.55)                              
               (   4 O   2p0 :    -0.36) (   5 H   1s  :    -0.21) (   6 H   1s  :    -0.21)                              
                                                                                                                          
               Molecular Orbital No.   8:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -0.57109 a.u.                                                                  
               (   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.   9:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -0.49457 a.u.                                                                  
               (   1 O   1p+1:    -0.30) (   1 O   2p+1:    -0.23) (   4 O   1p+1:    -0.56)                              
               (   4 O   2p+1:    -0.44)                                                                                  
                                                                                                                          
               Molecular Orbital No.  10:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -0.49457 a.u.                                                                  
               (   1 O   1p+1:    -0.56) (   1 O   2p+1:    -0.44) (   4 O   1p+1:     0.30)                              
               (   4 O   2p+1:     0.23)                                                                                  
                                                                                                                          
               Molecular Orbital No.  11:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.18787 a.u.                                                                  
               (   1 O   3s  :     1.01) (   1 O   1p0 :     0.19) (   1 O   2p0 :     0.33)                              
               (   2 H   2s  :    -0.83) (   3 H   2s  :    -0.83)                                                        
                                                                                                                          
               Molecular Orbital No.  12:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.18787 a.u.                                                                  
               (   4 O   3s  :     1.01) (   4 O   1p0 :     0.19) (   4 O   2p0 :     0.33)                              
               (   5 H   2s  :    -0.83) (   6 H   2s  :    -0.83)                                                        
                                                                                                                          
               Molecular Orbital No.  13:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.25852 a.u.                                                                  
               (   4 O   1p-1:    -0.28) (   4 O   2p-1:    -0.67) (   5 H   2s  :     1.48)                              
               (   6 H   2s  :    -1.48)                                                                                  
                                                                                                                          
               Molecular Orbital No.  14:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.25852 a.u.                                                                  
               (   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.  15:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.79749 a.u.                                                                  
               (   4 O   1p-1:    -0.26) (   4 O   2p-1:    -0.49) (   5 H   1s  :     0.95)                              
               (   5 H   2s  :    -0.66) (   5 H   1p0 :     0.16) (   6 H   1s  :    -0.95)                              
               (   6 H   2s  :     0.66) (   6 H   1p0 :    -0.16)                                                        
                                                                                                                          
               *** MP2 correlation energy:      -0.203012702573 a.u.                                                      
                                                                                                                          
print("MP2 Energy correction of 2 water molecules", mp2drv.e_mp2)
print("Twice the energy of 1 water molecule", (e_mp2_os + e_mp2_ss) * 2)
MP2 Energy correction of 2 water molecules -0.4060254139814834
Twice the energy of 1 water molecule -0.4060254051461055

We can see that the two energies match. MP2 is size consistent! 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 normalisation, lowering the weight of these determinants as the size of the system increases.

Further reading#

  • 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.