13.3. The restricted frozen-pair 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 frozen-pair coupled-cluster correction on top of a RpCCD (restricted pair coupled cluster doubles) 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.3.1. Quick Guide: RfpCCD and RfpCCSD¶

This is a short introduction to the RfpCC 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.3.1.1. ROOpCCD 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 code snippet below shows how to perform a variational orbital-optimized pCCD calculations followed by an RfpCCSD calculation with the ROOpCCD reference function in PyBEST.

pccd = ROOpCCD(lf, occ_model)
pccd_output = pccd(hf_output, kin, ne, eri)
fpccsd = RfpCCSD(lf, occ_model)
fpccsd_output = fpccsd(kin, ne, eri, pccd_output)

Note

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

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

e_tot

the total RfpCCSD energy (sum of ROOpCCD energy and RfpCC energy)

e_corr

the RfpCCSD correlation energy

e_corr_s

the single excitation contribution in RfpCCSD correlation energy

e_corr_d

the double excitation contribution in RfpCCSD correlation energy

e_ref

the energy of the reference determinant (here, the pCCD determinant energy)

t_1

the single excitation cluster amplitudes

t_2

the double excitation cluster amplitudes

method

the name of the model (here: RfpCCSD)

nocc

the number of all occupied orbitals (including core orbitals)

nvirt

the number of all virtual orbitals

ncore

the number of (frozen) core orbitals (doubly 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 RfpCCD class. The overall input structure is equivalent to RfpCCSD.

fpccd = RfpCCD(lf, occ_model)
fpccd_output = fpccd(kin, ne, eri, pccd_output)

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

13.3.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 an fpCCSD calculation, specify the ncore argument during the initialization of some occupation model class,

# Select 4 frozen core orbital
#-----------------------------
occ_model = AufbauOccModel(gobasis, ncore=4)
# Perform CC calculation, ncore is stored in occ_model
#-----------------------------------------------------
fpccsd = RfpCCSD(lf, occ_model)
fpccsd_output = fpccsd(kin, ne, eri, pccd_output)

This syntax is working for all fpCC modules mentioned above.

13.3.1.3. Restart options¶

The RfpCCD and RfpCCSD 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.

fpccd = RfpCCD(lf, occ_model)
# Restart only amplitudes from checkpoint_RCCD_d.h5
fpccd_output = fpccd(
    kin, ne, eri, pccd_output, restart='pybest-results/checkpoint_RfpCCD.h5'
)

PyBEST will read the amplitudes stored in pybest-results/checkpoint_RfpCCD.h5 and provide those as guess amplitudes to the RfpCCD 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.3.1.4. Summary of keyword arguments¶

The RfpCCD and RfpCCSD classes support 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.

e_ref

(float) reference energy - if you want to specify it manually instead of using the pCCD output IOData container

filename

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

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

restart

(str) path to file containing some (un)converged CC result.

solver

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

t_p

(float) pair-excitation doubles amplitudes - if you want to specify it manually instead of using the pCCD output container

threshold

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

13.3.1.5. RfpCCD and RfpCCSD calculation on the water molecule¶

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

Listing 13.8 data/examples/rccsd/fpccsd_water.py¶

from pybest.cc import RfpCCD, RfpCCSD
from pybest.context import context
from pybest.gbasis import (
    compute_cholesky_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
from pybest.geminals import ROOpCCD
from pybest.linalg import CholeskyLinalgFactory
from pybest.localization import PipekMezey
from pybest.occ_model import AufbauOccModel
from pybest.part import get_mulliken_operators
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 = CholeskyLinalgFactory(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_cholesky_eri(obasis, threshold=1e-8)
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)
###############################################################################
# Localize orbitals to improve pCCD convergence ###############################
###############################################################################
mulliken = get_mulliken_operators(obasis)
loc = PipekMezey(lf, occ_model, mulliken)
loc(orb_a, "occ")
loc(orb_a, "virt")
###############################################################################
# Do pCCD #####################################################################
###############################################################################
pccd = ROOpCCD(lf, occ_model)
pccd_output = pccd(hf_output, kin, ne, eri, e_core=nr)
###############################################################################
### Do pCCD-fpCCD calculation #################################################
###############################################################################
ccd = RfpCCD(lf, occ_model)
ccd_output = ccd(pccd_output, kin, ne, eri)
###############################################################################
### Do pCCD-fpCCSD calculation ################################################
###############################################################################
ccsd = RfpCCSD(lf, occ_model)
ccsd_output = ccsd(pccd_output, kin, ne, eri)