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