12.2. The 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 an RHF (see The Self-Consistent Field Module) reference function. Both modules are explained in greater detail below.

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

12.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_ (see The Self-Consistent Field Module). The code snippet below shows how to perform an RCCSD calculation with an RHF reference function in PyBEST,

ccsd = RCCSD(lf, occ_model)
ccsd_ = ccsd(kin, ne, eri, hf_)

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 RHFLCCD energy (including all external terms)

e_corr

The RHFLCCD correlation energy

e_ref

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

t_1

The single excitation cluster amplitudes

t_2

The double excitation cluster amplitudes

Note

The order of arguments does not matter because PyBEST assigns them in an automatic fashion.

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_ = ccd(kin, ne, eri, hf_)

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

12.2.1.2. 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 RCC 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 RCCSD

# Select 4 frozen core orbital
#-----------------------------
ccsd = RCCSD(lf, occ_model, ncore=4)
ccsd_ = ccsd(kin, ne, eri, hf_)

This syntax is working for all CC modules mentioned above.

12.2.1.3. Restart options¶

The RCC module supports restarts from PyBEST’s internal checkpoint files. If your RCC calculation did not converge and you wish to restart it and take the amplitudes from a previous calculation as guess, specify the path for checkpoint file using initguess keyword.

ccd = RCCD(lf, occ_model)
# Restart only amplitudes from checkpoint_RCCD_d.h5
ccd_ = ccd(
    kin, ne, eri, hf_, initguess='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 RCC 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.

12.2.1.4. Summary of keyword arguments¶

The RCC module 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.

checkpoint

(int) frequency of checkpointing

checkpoint_path

(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 amplitudes, one of: “random” - random numbers, “mp2” - performs MP2 calculations before, “const” - uses constant number, used for debugging, path to file - path to .h5 file containing amplitudes (saved as checkpoint_method.h5 by default), numpy.ndarray - 1d vector with amplitudes.

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

12.2.1.5. RCC 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 12.5 data/examples/rccsd/rccsd_frozen_core_water_cc-pvdz.py¶

from pybest.cc.rccsd import RCCD, RCCSD
from pybest.context import context
from pybest.gbasis import get_gobasis
from pybest.gbasis.gobasis import (
    compute_kinetic,
    compute_nuclear,
    compute_eri,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.linalg.dense import DenseLinalgFactory
from pybest.scf import AufbauOccModel
from pybest.wrappers.hf 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(5)
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_ = hf(kin, ne, eri, external, olp, orb_a)

###############################################################################
### Do RHF-CCD calculation and 1 frozen core orbital ##########################
###############################################################################
ccd = RCCD(lf, occ_model, ncore=1)
ccd_ = ccd(kin, ne, eri, hf_)

###############################################################################
### Do RHF-CCSD calculation and 1 frozen core orbital #########################
###############################################################################
ccsd = RCCSD(lf, occ_model, ncore=1)
ccsd_ = ccsd(kin, ne, eri, hf_)

12.2.1.6. RCC calculation on the water molecule¶

This is a basic example on how to perform a CCSD 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 12.6 data/examples/rccsd/rccsd_water_cc-pvdz.py¶

from pybest.cc.rccsd import RCCD, RCCSD
from pybest.context import context
from pybest.gbasis import get_gobasis
from pybest.gbasis.gobasis import (
    compute_kinetic,
    compute_nuclear,
    compute_eri,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.linalg.dense import DenseLinalgFactory
from pybest.scf import AufbauOccModel
from pybest.wrappers.hf 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(5)
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_ = hf(kin, ne, eri, external, olp, orb_a)

###############################################################################
### Do RHF-CCD calculation ####################################################
###############################################################################
ccd = RCCD(lf, occ_model)
ccd_ = ccd(kin, ne, eri, hf_)

###############################################################################
### Do RHF-CCSD calculation ###################################################
###############################################################################
ccsd = RCCSD(lf, occ_model)
ccsd_ = ccsd(kin, ne, eri, hf_)

12.2.1.7. 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 12.7 data/examples/rccsd/rccsd_frozen_core_water_cc-pvdz.py¶

from pybest.cc.rccsd import RCCD, RCCSD
from pybest.context import context
from pybest.gbasis import get_gobasis
from pybest.gbasis.gobasis import (
    compute_kinetic,
    compute_nuclear,
    compute_eri,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.linalg.dense import DenseLinalgFactory
from pybest.scf import AufbauOccModel
from pybest.wrappers.hf 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(5)
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_ = hf(kin, ne, eri, external, olp, orb_a)

###############################################################################
### Do RHF-CCSD calculation that do not converge  #############################
###############################################################################
ccsd = RCCSD(lf, occ_model)
ccsd_ = ccsd(kin, ne, eri, hf_, maxiter=1)

###############################################################################
### Do RHF-CCSD calculation with amplitides from file #########################
###############################################################################
ccsd = RCCSD(lf, occ_model)
ccsd_ = ccsd(kin, ne, eri, hf_, initguess="pybest-results/checkpoint_RCCSD.h5")