8. The Moller-Plesset Perturbation Theory module

Before reading this section, please read the general introduction mentioned in General Remarks concerning Post-Hartree-Fock Calculations. This part of the Documentation builds upon it.

8.1. Supported features

The MP2 module supports spin-restricted orbitals and the DenseLinalgFactory and CholeskyLinalgFactory. PyBEST offers additional flavors of perturbation theory with an pCCD reference function (see The perturbation theory module).

8.2. Quick Guide: conventional RMP2

We assume that you have performed a restricted Hartree-Fock calculation, whose results are stored in the IOData container hf_output (see The Self-Consistent Field Module). Furthermore, we will assume the following names for all PyBEST objects

lf:

A LinalgFactory instance (see Preliminaries).

occ_model:

An Aufbau occupation model of the AufbauOccModel class

kin:

the kinetic energy integrals

ne:

the nucleus-electron attraction integrals

eri:

the two-electron repulsion integrals

The code snippet below shows how to perform an MP2 calculation in PyBEST,

mp2 = RMP2(lf, occ_model)
mp2_output = mp2(kin, ne, eri, hf_output)

The results are returned as an IOData container, while all results are saved to the checkpoint_MP2_d.h5 file. Specifically, the IOData container contains the following attributes

e_tot:

The total MP2 energy (including all external terms)

e_corr:

The MP2 correlation energy

e_ref:

The energy of the reference determinant (here, the Hartree-Fock determinant)

t_2:

The doubles amplitudes

The MP2 module features additional options, which are discussed in greater detail below.

8.2.1. Defining a frozen core

All post-RHF modules in PyBEST support frozen core orbitals. These orbitals are not included in the amplitude equations. The frozen (core) orbitals are thus by construction doubly occupied. The number of frozen orbitals can be defined when creating an instance of the OccupationModel class,

orb_a = lf.create_orbital(obasis.nbasis)
## Do RMP2 calculation with 1 frozen core orbital #############################
###############################################################################
mp2 = RMP2(lf, occ_model)
mp2_output = mp2(kin, ne, eri, hf_output)

Note

If ncore is set to some nonzero number, the orbitals are frozen with respect to their order in the orb_a attribute passed to RMP2. In general, these are the energetically lowest-lying orbitals or those with the largest occupation numbers.

The number of virtual orbitals cannot be change and hence PyBEST assumes all virtual orbitals to be active. All features of the RMP2 module discussed below work with frozen core orbitals.

8.3. RMP2 with single excitations

The conventional MP2 implementation only accounts for double excitations. The contributions of single excitations can be included by setting the keyword singles to True,

## Do RMP2 calculation with single excitations included #######################
###############################################################################
mp2 = RMP2(lf, occ_model)
mp2_output = mp2(kin, ne, eri, hf_output, singles=True)

Note

If canonical Hartree-Fock orbitals are used in MP2 calculations, the singles correction does not have any affect and the energy contribution is thus zero.

8.4. RMP2 (relaxed) density matrices and RMP2 natural orbitals

PyBEST also supports the calculation of (relaxed and unrelaxed) 1-particle reduced density matrices (1RDM) and the corresponding natural orbitals. By default, those are not calculated. To determine the MP2 natural occupation numbers, activate the natorb keyword,

## Do RMP2 calculation AND determine natural orbitals #########################
###############################################################################
mp2 = RMP2(lf, occ_model)
mp2_output = mp2(kin, ne, eri, hf_output, natorb=True)

The orbital relaxation part is by default neglected. Orbital relaxation effects to the MP2 1RDM are accounted for by activating the relaxation keyword,

## Do RMP2 calculation AND determine natural orbitals INCLUDING ###############
## orbital relaxation contributions ###########################################
###############################################################################
mp2 = RMP2(lf, occ_model)
mp2_output = mp2(kin, ne, eri, hf_output, natorb=True, relaxation=True)

The MP2 natural orbitals are saved to the checkpoint file and copied to the IOData container mp2_output. Thus, the natural MP2 orbitals can be accesses as follows

mp2_output.orb_a

Note

The relaxed 1RDM can only be determined for the conventional MP2 method, that is, for double excitations only.

8.5. Spin component scaled (SCS) RMP2

All above features can also be combined with spin component scaling, where the contributions (to the energy or the relaxed and unrelaxed 1RDM) of the same spin and opposite spin component are calculated. Each scaling factor is passed as an additional keyword argument (fos or fss) to the RMP2 function call, for instance,

## Do SCS RMP2 calculation with fos=0.3 and fss=1.2 ###########################
###############################################################################
mp2 = RMP2(lf, occ_model)
mp2_output = mp2(kin, ne, eri, hf_output, fos=0.3, fss=1.2)

Similarly, spin-component-scaling can be combined with the calculation of (relaxed and unrelaxed) 1RDMs and natural occupation numbers,

## Do SCS RMP2 calculation with fos=0.3 and fss=1.2 ###########################
## AND determine natural orbitals INCLUDING orbital relaxation contributions ##
###############################################################################
mp2 = RMP2(lf, occ_model)
mp2_output = mp2(
    kin, ne, eri, hf_output, natorb=True, relaxation=True, fos=0.3, fss=1.2
)

8.6. Example Python scripts

Several complete examples can be found in the directory data/examples/perturbation_theory.

8.6.1. MP2 calculation on the water molecule

This is a basic example on how to perform a RMP2 calculation in PyBEST. This script performs several RMP2 calculation on the water molecule using the cc-pVDZ basis set and (i) only double excitations, (ii) double and single excitations, (iii) double excitations and the calculation of unrelaxed natural orbitals and occupations numbers, and (iv) double excitations and relaxed natural orbitals and occupations numbers.

Listing 8.1 data/examples/perturbation_theory/mp2_water_cc-pvdz.py

from pybest import context
from pybest.gbasis.gobasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.pt import RMP2
from pybest.wrappers import RHF

###############################################################################
## Set up molecule, define basis set ##########################################
###############################################################################
# get the XYZ file from PyBEST's test data directory
fn_xyz = context.get_fn("test/water.xyz")
obasis = get_gobasis("cc-pvdz", fn_xyz)
###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
occ_model = AufbauOccModel(obasis)
orb_a = lf.create_orbital(obasis.nbasis)
olp = compute_overlap(obasis)
###############################################################################
## Construct Hamiltonian ######################################################
###############################################################################
kin = compute_kinetic(obasis)
ne = compute_nuclear(obasis)
eri = compute_eri(obasis)
external = compute_nuclear_repulsion(obasis)

###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, external, olp, orb_a)

###############################################################################
## Do RMP2 calculation ########################################################
###############################################################################
mp2 = RMP2(lf, occ_model)
mp2_output = mp2(kin, ne, eri, hf_output)

###############################################################################
## Do RMP2 calculation with single excitations included #######################
###############################################################################
mp2 = RMP2(lf, occ_model)
mp2_output = mp2(kin, ne, eri, hf_output, singles=True)

###############################################################################
## Do RMP2 calculation AND determine natural orbitals #########################
###############################################################################
mp2 = RMP2(lf, occ_model)
mp2_output = mp2(kin, ne, eri, hf_output, natorb=True)

###############################################################################
## Do RMP2 calculation AND determine natural orbitals INCLUDING ###############
## orbital relaxation contributions ###########################################
###############################################################################
mp2 = RMP2(lf, occ_model)
mp2_output = mp2(kin, ne, eri, hf_output, natorb=True, relaxation=True)

8.6.2. SCS-MP2 calculation on the water molecule

This is a basic example on how to perform a SCS-RMP2 calculation in PyBEST. This script performs several RMP2 calculation on the water molecule using the cc-pVDZ basis set and (i) only double excitations and (ii) double excitations and the calculation of unrelaxed natural orbitals and occupations numbers, and (iv) double excitations and relaxed natural orbitals and occupations numbers.

Listing 8.2 data/examples/perturbation_theory/scsmp2_water_cc-pvdz.py

from pybest import context
from pybest.gbasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.pt import RMP2
from pybest.wrappers import RHF

###############################################################################
## Set up molecule, define basis set ##########################################
###############################################################################
# get the XYZ file from PyBEST's test data directory
fn_xyz = context.get_fn("test/water.xyz")
obasis = get_gobasis("cc-pvdz", fn_xyz)
###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
occ_model = AufbauOccModel(obasis)
orb_a = lf.create_orbital(obasis.nbasis)
olp = compute_overlap(obasis)
###############################################################################
## Construct Hamiltonian ######################################################
###############################################################################
kin = compute_kinetic(obasis)
ne = compute_nuclear(obasis)
eri = compute_eri(obasis)
external = compute_nuclear_repulsion(obasis)

###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, external, olp, orb_a)

###############################################################################
## Do SCS RMP2 calculation with fos=0.3 and fss=1.2 ###########################
###############################################################################
mp2 = RMP2(lf, occ_model)
mp2_output = mp2(kin, ne, eri, hf_output, fos=0.3, fss=1.2)

###############################################################################
## Do SCS RMP2 calculation with fos=0.3 and fss=1.2 ###########################
## AND determine natural orbitals INCLUDING orbital relaxation contributions ##
###############################################################################
mp2 = RMP2(lf, occ_model)
mp2_output = mp2(
    kin, ne, eri, hf_output, natorb=True, relaxation=True, fos=0.3, fss=1.2
)