11. The 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.

11.1. Supported features

The perturbation theory module supports spin-restricted orbitals, the DenseLinalgFactory and CholeskyLinalgFactory. PyBEST offers various flavors of perturbation theory with a pCCD reference function (beyond MP2).

If a dynamical correlation is added on top of pCCD, different choices for the zeroth-order Hamiltonian \(\hat{H}_0\), perturbation \(\hat{V}\), dual space \(\langle \tilde{\Psi} \vert\), and projection manifold \(\hat{T}\) are available [boguslawski2017a]. Specifically, you can choose between the following combinations (see also [boguslawski2017a] for more details),

Model	\(\hat{H}_0\)	\(\hat{V}\)	\(\langle \tilde{\Psi} \vert\)	\(\hat{T}\)
PT2SDd	\(\hat{F}_N^{\mathrm{d}}\)	\({F}_N^{\mathrm{o}} + \hat{W}_N^\prime\)	\(\langle 0 \vert\)	\(\hat{T}_2^\prime\), \(\hat{T}_1+\hat{T}_2^\prime\)
PT2MDd	\(\hat{F}_N^{\mathrm{d}}\)	\({F}_N^{\mathrm{o}} + \hat{W}_N^\prime\)	\(\langle \mathrm{pCCD} \vert\)	\(\hat{T}_2^\prime\), \(\hat{T}_1+\hat{T}_2^\prime\)
PT2SDo	\(\hat{F}_N^{\mathrm{d}}\) \(+\hat{F}_N^{\mathrm{o}}\)	\(\hat{W}_N^\prime\)	\(\langle 0 \vert\)	\(\hat{T}_2^\prime\), \(\hat{T}_1+\hat{T}_2^\prime\)
PT2MDo	\(\hat{F}_N^{\mathrm{d}}\) \(+\hat{F}_N^{\mathrm{o}}\)	\(\hat{F}_N - \bar{F}_N\) \(+ \hat{W}_N^\prime\)	\(\langle \mathrm{pCCD} \vert\)	\(\hat{T}_2^\prime\), \(\hat{T}_1+\hat{T}_2^\prime\)
PT2b	\(\hat{F}_N^{\mathrm{d}}\) \(+\hat{F}_N^{\mathrm{o}}\)	\(\hat{H}_N^\prime\)	\(\langle \mathrm{pCCD} \vert\)	\(\hat{T}_2^\prime\), \(\hat{T}_1+\hat{T}_2^\prime\), \(\hat{T}_2\), \(\hat{T}_1+\hat{T}_2\)

For Moller-Plesset Perturbation Theory of second order with a restricted, closed-shell Hartree-Fock reference function, see The Moller-Plesset Perturbation Theory module.

11.2. Quick Guide: General input structure of PT2

We assume that you have performed a restricted pCCD calculation, whose results are stored in the IOData container oopccd_output (see The pCCD module). Furthermore, we will assume the following default 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

11.2.1. Setting up a PT2 calculation

The code snippet below shows how to perform a PT2 correction on top of pCCD in PyBEST,

pt2 = PT2XX(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output)

where PT2XX indicates one of the Perturbation Theory models listed in the above table.

Note

The order of the arguments in the PT2XX function call does not matter. PyBEST recognizes all required arguments automatically.

The results are returned as an IOData container, while all results are saved to the pybest-results/checkpoint_**PT2XX**.h5 file. Specifically, the IOData container contains the following attributes

e_tot:

The total PT2XX energy (including all external terms)

e_corr:

The PT2XX correlation energy

e_ref:

The energy of the reference wave function (here, OOpCCD)

t_1:

The singles amplitudes (if calculated)

t_2:

The doubles amplitudes

All PT2 modules feature additional options, which are discussed in greater detail below.

Note

None of the PT2 modules supports a restart option, yet.

11.2.2. Defining a frozen core

By default, all (occupied and virtual) orbitals are active. If a frozen core has been selected in the pCCD reference calculation, the same orbitals have to be frozen in the chosen flavor of PT2 correction. To freeze some (occupied) orbitals, the number of frozen cores has to be specified when creating an instance of some occupation model. The perturbation theory module extracts then this information,

# Select one frozen core orbital
#-------------------------------
occ_model = AufbauOccModel(gobasis, ncore=1)
# Perform PT2 calculation, ncore is stored in occ_model
#------------------------------------------------------
pt2 = PT2XX(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output)

This syntax is working for all PT2 modules mentioned below.

11.3. The PT2SD module

The PT2SD module features two distinct flavors of perturbative corrections, which differ in the choice of the zeroth-order Hamiltonian \(\hat{H}_0\): (i) PT2SDd and (ii) PT2SDo

11.3.1. The PT2SDd correction

To perform a PT2SDd correction, you have to create an instance of the PT2SDd class. Due to the choice of the zeroth-order Hamiltonian, the singles and doubles correction are uncoupled and are determined simultaneously,

pt2 = PT2SDd(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output)

The function call returns an IOData container. It has the following attributes

e_tot:

The total PT2SDd energy (including all external terms)

e_corr:

The PT2SDd correlation energy (singles+doubles)

e_corr_s:

The PT2SDd correlation energy wrt single excitations

e_corr_d:

The PT2SDd correlation energy wrt double excitations

e_corr_s0:

The PT2SDd correlation energy wrt the seniority zero sector of all double excitations

e_corr_s2:

The PT2SDd correlation energy wrt the seniority 2 sector of alldouble excitations

e_corr_s4:

The PT2SDd correlation energy wrt the seniority 4 sector of alldouble excitations

e_ref:

The energy of the reference wave function (here, OOpCCD)

t_1:

The singles amplitudes

t_2:

The doubles amplitudes (all pair amplitudes equal zero)

t_p:

The pair amplitudes of OOpCCD

Note

The excitation manifold consists of both single and double excitations as they are uncoupled. By default, both are determined simultaneously.

11.3.2. The PT2SDo correction

If an off-diagonal zeroth-order Hamiltonian \(\hat{H}_0\) is chosen, the single and double amplitudes equations are coupled. Thus, the total and correlation energies associated with each excitation manifold will be affected by the presence of the other manifold. In PyBEST, you have the option to consider the doubles and the singles and doubles excitation manifold separately.

11.3.2.1. The PT2SDo correction including only double excitations

If not specified otherwise, PyBEST includes only double excitations in a PT2SDo calculation. Similar to above, an instance of the PT2SDo class has to be created before a function call initiates the optimization of the PT2 amplitudes,

################################################################################
### Do PT2 corrections without single excitations ##############################
################################################################################
pt2 = PT2SDo(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output)

Again, the function call returns an IOData container. In contrast to the PT2SDd flavor, all information related to single excitations is missing and we have the following attributes

e_tot:

The total PT2SDo energy (including all external terms)

e_corr:

The PT2SDo correlation energy (singles+doubles)

e_corr_d:

The PT2SDo correlation energy wrt double excitations

e_corr_s0:

The PT2SDo correlation energy wrt the seniority zero sector of all double excitations

e_corr_s2:

The PT2SDo correlation energy wrt the seniority 2 sector of alldouble excitations

e_corr_s4:

The PT2SDo correlation energy wrt the seniority 4 sector of alldouble excitations

e_ref:

The energy of the reference wave function (here, OOpCCD)

t_2:

The doubles amplitudes (all pair amplitudes equal zero)

t_p:

The pair amplitudes of OOpCCD

11.3.2.2. The PT2SDo correction including single and double excitations

The contribution of single excitations can be switched on using the singles keyword,

################################################################################
### Do PT2 corrections with single excitations #################################
################################################################################
pt2 = PT2SDo(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output, singles=True)

The function call then returns an IOData container, where, in addition, to the attributes above you can access the following internal PyBEST objects

e_corr_s:

The PT2SDo correlation energy wrt single excitations

t_1:

The singles amplitudes

11.3.2.3. Summary of keyword arguments

indextrans:

(str) 4-index transformation. The choice between cupy, tensordot (default) and einsum. tensordot is faster than einsum. If DenseLinalgFactory is used, the memory requirement scales roughly as \(3N^4\). Due to the storage of the two-electron integrals, the total amount of memory increases to \(4N^4\)

Note

If Cupy is not available or unsuccessful, td is selected instead.

eref:

(float) pCCD reference energy (default 0.0). eref gets overwritten if an IOData container is passed

threshold:

(float) optimization threshold for amplitudes (default 1e-6)

maxiter:

(int) maximum number of iterations (default 200)

singles:

(boolean) account for single excitations in projection manifold (default False)

guess:

(1-dim np.array) initial guess (default None). In not provided, an initial guess containing the MP2 amplitudes is generated. If a user-defined guess is provided, the elements of \(t_{ij}^{ab}\) have to be indexed in C-like order

11.4. The PTMD module

Similar to PT2SD, the PTMD module features two distinct flavors of perturbative corrections, which differ in the choice of the zeroth-order Hamiltonian \(\hat{H}_0\): (i) PT2MDd and (ii) PT2MDo

11.4.1. The PT2MDd correction

To perform a PT2MDd correction, you have to create an instance of the PT2MDd class. Similar to the PT2SDd method, the choice of the zeroth-order Hamiltonian decouples the singles and doubles equations and they are thus determined simultaneously by default. Since the PT2MDd method requires (an approximate) wave function overlap of the pCCD reference function as input, this integral has to be calculated prior the PT2MDd function call. It is defined in the RpCCD and ROOpCCD classes. PyBEST automatically calculates and stores this overlap integral. Note that this is only an approximate overlap integral where at most quadratic terms in the pCCD amplitudes are included. The PT2MDd calculation can be started as follows

pt2 = PT2MDd(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output)

The function call returns an IOData container. It has similar attributes as the PT2SDd IOData container,

e_tot:

The total PT2MDd energy (including all external terms)

e_corr:

The PT2MDd correlation energy (singles+doubles)

e_corr_s:

The PT2MDd correlation energy wrt single excitations

e_corr_d:

The PT2MDd correlation energy wrt double excitations

e_corr_s0:

The PT2MDd correlation energy wrt the seniority zero sector of all double excitations

e_corr_s2:

The PT2MDd correlation energy wrt the seniority 2 sector of alldouble excitations

e_corr_s4:

The PT2MDd correlation energy wrt the seniority 4 sector of alldouble excitations

e_ref:

The energy of the reference wave function (here, OOpCCD)

t_1:

The singles amplitudes

t_2:

The doubles amplitudes (all pair amplitudes equal zero)

t_p:

The pair amplitudes of OOpCCD

Note

The excitation manifold consists of both single and double excitations as they are uncoupled. By default, both are determined simultaneously.

11.4.2. The PT2MDo correction

If an off-diagonal zeroth-order Hamiltonian \(\hat{H}_0\) is chosen, the single and double amplitudes equations are coupled. Thus, the total and correlation energies associated with each excitation manifold will be affected by the presence of the other manifold. In PyBEST, you have the option to consider the doubles and the singles and doubles excitation manifold separately.

11.4.2.1. The PT2MDo correction including only double excitations

As in the case of PT2SDo calculations, PyBEST includes only double excitations in a PT2MDo calculation if not specified otherwise,

pt2 = PT2MDo(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output)

Again, the function call returns an IOData container. In contrast to the PT2MDd flavor, we have only the following attributes

e_tot:

The total PT2MDo energy (including all external terms)

e_corr:

The PT2MDo correlation energy (singles+doubles)

e_corr_d:

The PT2MDo correlation energy wrt double excitations

e_corr_s0:

The PT2MDo correlation energy wrt the seniority zero sector of all double excitations

e_corr_s2:

The PT2MDo correlation energy wrt the seniority 2 sector of alldouble excitations

e_corr_s4:

The PT2MDo correlation energy wrt the seniority 4 sector of alldouble excitations

e_ref:

The energy of the reference wave function (here, OOpCCD)

t_2:

The doubles amplitudes (all pair amplitudes equal zero)

t_p:

The pair amplitudes of OOpCCD

where all information related to single excitations is missing.

11.4.2.2. The PT2MDo correction including single and double excitations

The contribution of single excitations can be switched on using the singles keyword,

pt2 = PT2MDo(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output, singles=True)

The function call then returns an IOData container, where we can access the additional internal PyBEST objects

e_corr_s:

The PT2MDo correlation energy wrt single excitations

t_1:

The singles amplitudes

11.4.2.3. Summary of keyword arguments

indextrans:

(str) 4-index transformation. The choice between cupy, tensordot (default) and einsum. tensordot is faster than einsum. If DenseLinalgFactory is used, the memory requirement scales roughly as \(3N^4\). Due to the storage of the two-electron integrals, the total amount of memory increases to \(4N^4\)

Note

If Cupy is not available or unsuccessful, td is selected instead.

eref:

(float) pCCD reference energy (default 0.0). eref gets overwritten if an IOData container is passed

threshold:

(float) optimization threshold for amplitudes (default 1e-6)

maxiter:

(int) maximum number of iterations (default 200)

singles:

(boolean) account for single excitations in projection manifold (default False)

guess:

(1-dim np.array) initial guess (default None). In not provided, an initial guess containing the MP2 amplitudes is generated. If a user-defined guess is provided, the elements of \(t_{ij}^{ab}\) have to be indexed in C-like order

11.5. The PT2b module

In the PT2b module, you have the choice between four different flavors of the PT2b correction, including the origanally formulated PTb model [limacher2014] and three extensions thereof. These Perturbation Theory corrections differ in their choice of the projection manifold.

11.5.1. The original PT2b correction (PTb)

In its original formulation, PT2b (also called PTb) includes all double excitations in the projection manifold, that is both electron-pair and broken-pair doubles are contained in the excitation operator of the PT2b/PTb model. Although these electron-pair contributions do not affect the energy expression (as the corresponding terms cancel), the corresponding pair amplitudes are optimized in the PT2b/PTb amplitude equations.

11.5.1.1. The PT2b correction including all double excitations

If not specified otherwise, the PT2b module optimizes the originally formulated PT2b/PTb wave function (that is, its amplitudes). To invoke a PT2b calculation, only the pCCD output container oopccd_output is required as an input argument,

# Do the actual PT2b calculation without singles/with pairs
# ----------------------------------------------------------
pt2 = PT2b(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output)

Again, the function call returns an IOData container with the following attributes

e_tot:

The total PT2b energy (including all external terms)

e_corr:

The PT2b correlation energy (singles+doubles)

e_corr_d:

The PT2b correlation energy wrt double excitations

e_corr_s0:

The PT2b correlation energy wrt the seniority zero sector of all double excitations

e_corr_s2:

The PT2b correlation energy wrt the seniority 2 sector of alldouble excitations

e_corr_s4:

The PT2b correlation energy wrt the seniority 4 sector of alldouble excitations

e_ref:

The energy of the reference wave function (here, OOpCCD)

t_2:

The doubles amplitudes (including pairs)

t_p:

The pair amplitudes of OOpCCD

Note

In contrast to the PT2SD and PT2MD models, the pair amplitudes of the original PT2b method are optimized and thus can have non-zero values.

11.5.1.2. The PT2b correction including all single and double excitations

Single excitations can be included in the PT2b model using the singles keyword in the function call,

# Do the actual PT2b calculation with singles/with pairs
# -------------------------------------------------------
pt2 = PT2b(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output, singles=True)

The corresponding IOData container (return value) contains, in addition, the single excitation amplitudes and the corresponding correlation energy contribution,

e_corr_s:

The PT2b correlation energy wrt single excitations

t_1:

The singles amplitudes

11.5.2. The modified PT2b correction

Since the original PT2b model does include electron pair excitations in the excitation operator of Perturbation Theory, double-counting effects cannot be excluded. Although these pair amplitudes do not contribute to the correlation energy expression, they do couple to the remaining double and pCCD pair amplitudes. Thus, electron-pair effects are treated twice, once in the pCCD reference and once again in the PT2b model. The modified PT2b correction in PyBEST offers the possibility to exclude those electron-pair excitations from the projection manifold in the Perturbation Theory model. Similar to the original PT2b/PTb model, we can distinguish two special cases.

11.5.2.1. The PT2b correction including only broken-pair double excitations

The electron-pair excitations can be excluded from the PT2b/PTb model using the keyword excludepairs. By default, this keyword is assumed to be False and thus a conventional PT2b/PTb calculation is performed. The code snippet below shows how to modify the PT2b function call to exclude electron-pair excitations from the doubles excitation manifold,

# Do the actual PT2b calculation without singles/without pairs
# -------------------------------------------------------------
pt2 = PT2b(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output, excludepairs=True)

The function call returns an IOData container, which as the same attributes as the original PT2b class as described in The PT2b correction including all double excitations, except for the doubles amplitudes,

t_2:

The doubles amplitudes (all pair amplitudes equal zero)

11.5.2.2. The PT2b correction including all single and only broken-pair double excitations

As mentioned in The PT2b correction including all single and double excitations, single excitations can be included in the PT2b model using the singles keyword in the function call, while all electron-pair contributions are excluded by activating the excludepairs keyword,

# Do the actual PT2b calculation with singles/without pairs
# ----------------------------------------------------------
pt2 = PT2b(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output, singles=True, excludepairs=True)

The corresponding IOData container (return value) contains, in addition, the single excitation amplitudes and the corresponding correlation energy contribution,

e_corr_s:

The PT2b correlation energy wrt single excitations

t_1:

The singles amplitudes

11.5.3. Summary of keyword arguments

indextrans:

(str) 4-index transformation. The choice between cupy, tensordot (default) and einsum. tensordot is faster than einsum. If DenseLinalgFactory is used, the memory requirement scales roughly as \(3N^4\). Due to the storage of the two-electron integrals, the total amount of memory increases to \(4N^4\)

Note

If Cupy is not available or unsuccessful, td is selected instead.

eref:

(float) pCCD reference energy (default 0.0). eref gets overwritten if an IOData container is passed

threshold:

(float) optimization threshold for amplitudes (default 1e-6)

maxiter:

(int) maximum number of iterations (default 200)

singles:

(boolean) account for single excitations in projection manifold (default False)

excludepairs:

(boolean) exclude electron-pair excitations in projection manifold (default False)

guess:

(1-dim np.array) initial guess (default None). In not provided, an initial guess containing the MP2 amplitudes is generated. If a user-defined guess is provided, the elements of \(t_{ij}^{ab}\) have to be indexed in C-like order

11.6. Example Python scripts

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

11.6.1. PT2SDd calculation on the water molecule

This is a basic example of how to perform a PT2SDd calculation in PyBEST. This script performs a PT2SDd calculation on the water molecule using the cc-pVDZ basis set.

Listing 11.1 data/examples/perturbation_theory/pt2sdd_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.geminals import ROOpCCD
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.pt import PT2SDd
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 OO-pCCD optimization ####################################################
###############################################################################
oopccd = ROOpCCD(lf, occ_model)
oopccd_output = oopccd(kin, ne, eri, hf_output)

################################################################################
### Do PT2 corrections #########################################################
################################################################################
pt2 = PT2SDd(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output)

11.6.2. PT2SDo calculation on the water molecule

This is a basic example on how to perform a PT2SDo calculation in PyBEST. This script performs a series of PT2SDo calculations on the water molecule using the cc-pVDZ basis set.

Listing 11.2 data/examples/perturbation_theory/pt2sdo_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.geminals import ROOpCCD
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.pt import PT2SDo
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 OO-pCCD optimization ####################################################
###############################################################################
oopccd = ROOpCCD(lf, occ_model)
oopccd_output = oopccd(kin, ne, eri, hf_output)

################################################################################
### Do PT2 corrections without single excitations ##############################
################################################################################
pt2 = PT2SDo(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output)

################################################################################
### Do PT2 corrections with single excitations #################################
################################################################################
pt2 = PT2SDo(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output, singles=True)

11.6.3. PT2MDd calculation on the water molecule

This is a basic example on how to perform a PT2MDd calculation in PyBEST. This script performs a PT2MDd calculation on the water molecule using the cc-pVDZ basis set.

Listing 11.3 data/examples/perturbation_theory/pt2mdd_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.geminals import ROOpCCD
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.pt import PT2MDd
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 OO-pCCD optimization ####################################################
###############################################################################
oopccd = ROOpCCD(lf, occ_model)
oopccd_output = oopccd(kin, ne, eri, hf_output)

################################################################################
### Do PT2 corrections #########################################################
################################################################################
pt2 = PT2MDd(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output)

11.6.4. PT2MDo calculation on the water molecule

This is a basic example on how to perform a PT2MDo calculation in PyBEST. This script performs a series of PT2MDo calculations on the water molecule using the cc-pVDZ basis set.

Listing 11.4 data/examples/perturbation_theory/pt2mdo_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.geminals import ROOpCCD
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.pt import PT2MDo
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 OO-pCCD optimization ####################################################
###############################################################################
oopccd = ROOpCCD(lf, occ_model)
oopccd_output = oopccd(kin, ne, eri, hf_output)

###############################################################################
### Do PT2 corrections ########################################################
###############################################################################
# Do the actual PT2 calculation without singles ###############################
# -----------------------------------------------------------------------------
pt2 = PT2MDo(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output)

# Do the actual PT2 calculation with singles ##################################
# -----------------------------------------------------------------------------
pt2 = PT2MDo(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output, singles=True)

11.6.5. PT2b calculation on the water molecule

This is a basic example on how to perform a PT2b calculation in PyBEST. This script performs a series of PT2b calculations on the water molecule using the cc-pVDZ basis set.

Listing 11.5 data/examples/perturbation_theory/pt2mdo_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.geminals import ROOpCCD
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.pt import PT2b
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 OO-pCCD optimization ####################################################
###############################################################################
oopccd = ROOpCCD(lf, occ_model)
oopccd_output = oopccd(kin, ne, eri, hf_output)

################################################################################
### Do PT2 corrections #########################################################
################################################################################
# Do the actual PT2b calculation without singles/with pairs
# ----------------------------------------------------------
pt2 = PT2b(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output)

# Do the actual PT2b calculation with singles/with pairs
# -------------------------------------------------------
pt2 = PT2b(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output, singles=True)

# Do the actual PT2b calculation without singles/without pairs
# -------------------------------------------------------------
pt2 = PT2b(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output, excludepairs=True)

# Do the actual PT2b calculation with singles/without pairs
# ----------------------------------------------------------
pt2 = PT2b(lf, occ_model)
pt2_output = pt2(kin, ne, eri, oopccd_output, singles=True, excludepairs=True)