13.2. The restricted coupled-cluster 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. With PyBEST, you can perform a coupled-cluster correction on top of any single-reference wave function, eg., the restricted Hartree-Fock (RHF) (see The Self-Consistent Field Module) or ROOpCCD (restricted orbital-optimized pair coupled cluster doubles; see The pCCD module) reference wave function. The module is explained in greater detail below.

If you use this module, please cite [leszczyk2022].

13.2.1. Quick Guide: RCC

This is a short introduction to the RCC module. More information on the input and output structure can be found in the following sections. Similar to the previous modules, 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

13.2.1.1. RHF reference function

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). The following code snippet below shows how to perform an RCCSD calculation with an RHF reference function in PyBEST.

ccsd = RCCSD(lf, occ_model)
ccsd_output = ccsd(kin, ne, eri, hf_output)

The results are returned as an IOData container, while all results are saved to the pybest-results/checkpoint_RCCSD.h5 file. Specifically, the IOData container contains (amongst others) the following attributes

e_tot:

the total RCCSD energy (sum of RHF energy and RCC energy)

e_corr:

the RCCSD correlation energy

e_corr_s:

the single excitation contribution in RCCSD correlation energy

e_corr_d:

the double excitation contribution in RCCSD correlation energy

e_ref:

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

t_1:

the single excitation cluster amplitudes

t_2:

the double excitation cluster amplitudes

method:

the name of the model (here: RCCSD)

nocc:

the number of all occupied orbitals (including core orbitals)

nvirt:

the number of all virtual (unoccupied in the RHF model) orbitals

ncore:

the number of (frozen) core orbitals (occupied but not correlated)

nact:

the number of correlated occupied and unoccupied orbitals

converged:

boolean value indicating if CC converged

To exclude single excitations, you have to use the RCCD class. The overall input structure is equivalent to RCCSD.

ccd = RCCD(lf, occ_model)
ccd_output = ccd(kin, ne, eri, hf_output)

Note

The order of arguments does not matter while calling RCCD instance because PyBEST recognizes them in an automatic fashion.

All results are saved to the pybest-results/checkpoint_RCCD.h5 checkpoint file.

13.2.1.2. ROOpCCD reference function

The following code snippet shows how to perform an RCCSD calculation with an ROOpCCD reference function (that is, ROOpCCD orbitals and its reference determinant) in PyBEST.

pccd = ROOpCCD(lf, occ_model)
pccd_output = pccd(hf_output, kin, ne, eri)
ccsd = RCCSD(lf, occ_model)
ccsd_output = ccsd(kin, ne, eri, pccd_output)

13.2.1.3. Defining a frozen core

By default, all (occupied and virtual) orbitals are active. To freeze some (occupied) orbitals, the number of frozen cores has to be specified when an instance of some occupation model class is created. For instance, if you wish to freeze the first 4 (occupied) orbitals in a RHF-CCSD calculation, specify the ncore argument during the initialization of the occupation model,

# Select 4 frozen core orbital
#-----------------------------
occ_model = AufbauOccModel(gobasis, ncore=4)
# Perform CC calculation, ncore is stored in occ_model
#-----------------------------------------------------
ccsd = RCCSD(lf, occ_model)
ccsd_output = ccsd(kin, ne, eri, hf_output)

This syntax is working for all CC modules mentioned above.

13.2.1.4. Restart options

The RCCD and RCCSD modules support restarts from PyBEST’s internal checkpoint files. If your coupled-cluster calculation did not converge before the number of iterations exceeded some maximum number allowed and you wish to restart it and take the amplitudes from a previous calculation as an initial guess, specify the path to the corresponding checkpoint file using the restart keyword.

ccd = RCCD(lf, occ_model)
# Restart only amplitudes from checkpoint_RCCD_d.h5
ccd_output = ccd(
    kin, ne, eri, hf_output, restart='pybest-results/checkpoint_RCCD.h5'
)

PyBEST will read the amplitudes stored in pybest-results/checkpoint_RCCD.h5 and provide those as guess amplitudes to the RCCD solver.

Note

Currently, PyBEST does not dump the Hamiltonian to disk that is used in the CC calculation. Thus, if the user modifies the Hamiltonian, the calculation might be corrupted and special care is advised.

13.2.1.5. Summary of keyword arguments

The RCCD and RCCSD classes supports various keyword arguments that allow us to steer the optimization of the cluster amplitudes. In the following, all supported keyword arguments are listed together with their default values. Please note, that for most cases the default values should be sufficient to reach convergence.

filename:

(str) path to file where current solution is dumped (default pybest-results/checkpoint_method.h5).

e_ref:

(float) reference energy (default value 0.0).

initguess:

(str or numpy.ndarray) initial guess for the amplitudes, one of: “random”, “mp2”, “const”, path to .h5 file, or IOData instance.

maxiter:

(int) maximum number of iterations (default: 100).

solver:

(str) one of scipy.optimize.root solvers (default: “krylov”).

threshold:

(float) tolerance for amplitudes (default: 1e-8).

13.2.1.6. RCCD and RCCSD calculation on the water molecule

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

Listing 13.5 data/examples/rccsd/rccsd_water.py

from pybest.cc import RCCD, RCCSD
from pybest.context 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.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)
nr = compute_nuclear_repulsion(obasis)
###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, nr, olp, orb_a)
###############################################################################
### Do RHF-CCD calculation ####################################################
###############################################################################
ccd = RCCD(lf, occ_model)
ccd_output = ccd(kin, ne, eri, hf_output)
###############################################################################
### Do RHF-CCSD calculation ###################################################
###############################################################################
ccsd = RCCSD(lf, occ_model)
ccsd_output = ccsd(kin, ne, eri, hf_output)

13.2.1.7. RCCD and RCCSD calculation on the water molecule using a frozen core

This is a basic example on how to perform a CCD and CCSD calculation in PyBEST. This script performs a CCD calculation followed by a CCSD calculation on the water molecule using the cc-pVDZ basis set and freezing the lowest orbital.

Listing 13.6 data/examples/rccsd/rccsd_frozen_core_water.py

from pybest.cc import RCCD, RCCSD
from pybest.context 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.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, ncore=1)
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)
nr = compute_nuclear_repulsion(obasis)
###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, nr, olp, orb_a)
###############################################################################
### Do RHF-CCD calculation and 1 frozen core orbital (stored in occ_model) ####
###############################################################################
ccd = RCCD(lf, occ_model)
ccd_output = ccd(kin, ne, eri, hf_output)
###############################################################################
### Do RHF-CCSD calculation and 1 frozen core orbital (stored in occ_model) ###
###############################################################################
ccsd = RCCSD(lf, occ_model)
ccsd_output = ccsd(kin, ne, eri, hf_output)

13.2.1.8. RCCSD calculation using data from file as an initial guess for amplitudes

This is a basic example on how to perform a CCSD calculation in PyBEST. This script performs a CCSD calculation where the maximum number of iterations is too small to converge. This calculation is followed up by the CCSD calculations that use the previous results that were dumped to file.

Listing 13.7 data/examples/rccsd/rccsd_frozen_core_water.py

from pybest.cc import RCCSD
from pybest.context 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.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)
nr = compute_nuclear_repulsion(obasis)
###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, nr, olp, orb_a)
###############################################################################
### Do RHF-CCSD calculation that do not converge  #############################
###############################################################################
ccsd = RCCSD(lf, occ_model)
ccsd_output = ccsd(kin, ne, eri, hf_output, maxiter=1)
###############################################################################
### Do RHF-CCSD calculation with amplitides from file #########################
###############################################################################
ccsd = RCCSD(lf, occ_model)
ccsd_output = ccsd(
    kin, ne, eri, hf_output, initguess="pybest-results/checkpoint_RCCSD.h5"
)