12.1. The linearized 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 linearized coupled cluster correction on top of an RHF (see The Self-Consistent Field Module) or RpCCD (see The pCCD module) reference function. Both modules are explained in greater detail below.

12.1.1. Supported features¶

The perturbation theory module supports spin-restricted orbitals and the DenseLinalgFactory and CholeskyLinalgFactory.

12.1.2. Quick Guide: RLCC¶

This is a short introduction to the RLCC 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.1.2.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 RLCCD calculation with an RHF reference function in PyBEST,

lccd = RHFLCCD(lf, occ_model)
lccd_ = lccd(kin, ne, eri, hf_)

The results are returned as an IOData container, while all results are saved to the pybest-results/checkpoint_RHFLCCD.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_2

The doubles amplitudes

Note

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

To include single excitations as well, you have to use the RHFLCCSD class. The overall input structure is equivalent to RHFLCCD,

lccsd = RHFLCCSD(lf, occ_model)
lccsd_ = lccsd(kin, ne, eri, hf_)

The corresponding IOData container (return value) will also contain the singles amplitudes in addition to the RHFLCCD attributes,

t_1

The singles amplitudes

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

12.1.2.2. RpCCD reference function¶

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

An RLCC correction on top of a pCCD wave function (with and without orbital optimization) is performed in a very similar way as a conventional RHFLCC calculation. We assume that you have performed a restricted pCCD calculation, whose results are stored in the IOData container pccd_ (see The pCCD module). The code snippet below shows how to perform an RLCCD calculation with a pCCD reference function in PyBEST,

lccd = RpCCDLCCD(lf, occ_model)
lccd_ = lccd(kin, ne, eri, pccd_)

Instead of creating an instance of RHFLCCD, you have to use the RpCCDLCCD class. The results are returned as an IOData container, while all results are saved to the pybest-results/checkpoint_RpCCDLCCD.h5 file. Specifically, the IOData container contains similar attributes as the RHFLCCD IOData container. In addition, the pCCD pair amplitudes t_p are stored in the container. Thus, we have (amongst others) the following attributes

e_tot

The total RpCCDLCCD energy (including all external terms)

e_corr

The RpCCDLCCD correlation energy

e_ref

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

t_2

The doubles amplitudes (all pair amplitudes equal zero)

t_p

The pair amplitudes of pCCD

Note

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

To include single excitations as well, you have to use the RpCCDLCCSD class. The overall input structure is equivalent to RpCCDLCCD,

lccsd = RpCCDLCCSD(lf, occ_model)
lccsd_ = lccsd(kin, ne, eri, pccd_)

The corresponding IOData container (return value) will also contain the singles amplitudes in addition to the RpCCDLCCD attributes,

t_1

The singles amplitudes

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

12.1.2.3. 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 LCC correction. The only exception is the RHF-LCC module, where a froze core can be defined even if no frozen core was present in the RHF calculation. To freeze some (occupied) orbitals, the number of frozen cores has to be specified when an instance of some LCC class is created. For instance, if you wish to freeze the first 4 (occupied) orbitals in a RHF-LCCSD calculation, specify the ncore argument during the initialization of RHFLCCSD

# Select 4 frozen core orbital
#-----------------------------
lccsd = RHFLCCSD(lf, occ_model, ncore=4)
lccsd_ = lccsd(kin, ne, eri, hf_)

This syntax is working for all LCC modules mentioned above.

12.1.2.4. Restart options¶

The RLCC module supports restarts from PyBEST’s internal checkpoint files. Furthermore, two different restart options are supported

Restart only CC amplitudes (restart_t keyword)

Restart complete LCC calculation (restart keyword)

Note that only one option can be selected at a time. If you are not sure which restart option to use, it is always safer to restart from a set of CC amplitudes only (option 1). In order to restart a calculation, you have to specify the restart_t/restart keyword, which contains the name of the checkpoint file used for restarts. Since the restart behavior is similar in all RLCC modules, we select RHFLCCD as an example.

If your RLCC calculation did not converge and you wish to restart it and take the amplitudes from a previous calculation as guess, specify the restart_t keyword. Note that in this case, the IOData container of the reference calculation still has to be passed as an argument,

lccd = RHFLCCD(lf, occ_model)
# Restart only amplitudes from checkpoint_RHFLCCD_d.h5
lccd_ = lccd(kin, ne, eri, hf_, restart_t='pybest-results/checkpoint_RHFLCCD.h5')

PyBEST will read the amplitudes stored in pybest-results/checkpoint_RHFLCCD.h5 and provide those as guess amplitudes to the LCC solver. If you wish to read all informations from file (for instance, you want to skip the (RHF or pCCD) reference calculation), you need to substitute the IOData container in the function call (either an RHF or pCCD container) by the keyword restart that contains the file name of the checkpoint file,

lccd = RHFLCCD(lf, occ_model)
# Restart from previous file checkpoint_RHFLCCD_d.h5
lccd_ = lccd(kin, ne, eri, restart='pybest-results/checkpoint_RHFLCCD.h5')

In this case, PyBEST will take all required input arguments (for instance, the molecular orbitals, reference energy, initial guess for all amplitudes, and, if a pCCD reference is chosen, the frozen pair amplitudes).

Note

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

If the arguments of the function call still contain an IOData container of the reference calculation (RHF or pCCD), PyBEST will read the molecular orbitals from this container as well as from the restart file and then project the molecular orbitals read from disk on the molecular orbitals passed by the IOData container. Thus, the final orbitals will change if those molecular orbitals differ. To prevent this projection, use the first option (restart_t keyword), which tells PyBEST to read in only the amplitudes from the file and ignore all other information.

12.1.3. RLCC on top of RHF¶

Please read the Quick Guide documentation in RHF reference function before consulting this part of the documentation. This section on the RLCC module builds upon the Quick Guide above and contains only additional information that might be useful for the user.

12.1.3.1. RHF-LCCD calculations¶

In addition to the IOData container attributes mentioned in the Quick Guide above, the RHFLCCD container includes the following information

e_corr_d

The RHF-LCCD correlation energy wrt double excitations

e_corr_s0

The RHF-LCCD correlation energy wrt the seniority zero sector of all double excitations

e_corr_s2

The RHF-LCCD correlation energy wrt the seniority 2 sector of all double excitations

e_corr_s4

The RHF-LCCD correlation energy wrt the seniority 4 sector of all double excitations

orb_a

A copy of the orbitals used in the RHF-LCCD calculation

olp

The overlap integrals used in the RHF-LCCD calculation (for restart purposes)

12.1.3.2. RHF-LCCSD (CEPA(0)) calculations¶

If single excitations are included in the cluster operator, the corresponding correlation energy is stored as a separate component in the IOData container (in addition to the attributes mentioned in RHF-LCCD calculations)

e_corr_s

The PT2SDd correlation energy wrt single excitations

t_1

The singles amplitudes

12.1.3.3. Summary of keyword arguments¶

The RLCC 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. For convergence difficulties see Troubleshooting in LCC calculations.

restart

(str) if provided, a restart from a given checkpoint file is performed (see Restart options)

restart_t

(str) if provided, only the amplitudes are restarted from a given checkpoint file (see Restart options)

indextrans

(str) 4-index Transformation. Choice between tensordot (default) and einsum. tensordot is faster than einsum, but requires more memory. If DenseLinalgFactory is used, the memory requirement scales as \(2N^4\) for einsum and \(3N^4\) for tensordot, respectively. Due to the storage of the two-electron integrals, the total amount of memory increases to \(3N^4\) for einsum and \(4N^4\) for tensordot, respectively

guesstype

(str) type of initial guess. One of mp2 (default) or random (random numbers)

guess

(1-dim np.array) external guess for cluster amplitudes (default None). If provided, guesstype is ignored. The elements of array have to be indexed in C-like order with T1 amplitudes appearing before T2 amplitudes

Note

The restart and restart_t keywords have precedence over the guess keyword. This keyword will be deprecated in future versions.

solver

(str) type of solver used. Possible values are krylov (scipy solver) or pbqn (perturbation-based quasi-Newton solver) (default krylov)

Note

krylov is more stable but also much slower than pbqn. The convergence for the latter may break down for difficult cases.

maxiter

(int) maximum number of iterations (default 200)

tthreshold

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

Some keyword arguments are only working together with the pbqn solver. If the default solver krylov has been selected, they have no effect.

ethreshold

(float) optimization threshold for energy (default 1e-7)

jacobian

(int) the jacobian approximation used in the pbqn solver. Either 1 (only diagonal elements of Fock matrix) or 2 (contains additional amplitude-independent terms) are allowed (default 1)

diismax

(int) maximum number of DIIS vectors used in the DIIS algorithm of the pbqn solver (default 2)

diisstart

(int) first iteration when DIIS sets in (default 0)

diisreset

(boolean) reset DIIS vectors if DIIS diverges (default True)

12.1.4. RLCC on top of pCCD¶

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

Please read the Quick Guide documentation in RpCCD reference function before consulting this part of the documentation. This section on the RLCC module builds upon the Quick Guide above and contains only additional information that might be useful for the user.

12.1.4.1. RpCCD-LCCD calculations¶

The IOData container of the RpCCDLCCD module has the same attributes as described for the RHFLCCD method (see RHF-LCCD calculations) except that the t_2 amplitudes only contain all broken-pair amplitudes and the pCCD pair amplitudes are also returned, that is,

t_2

The doubles amplitudes (all pair amplitudes equal zero)

t_p

The pair amplitudes of pCCD

12.1.4.2. RpCCD-LCCSD calculations¶

As in the case of RHF-LCCSD, the RpCCDLCCSD module additionally return the correlation energy associated with all single excitations and the corresponding T1 amplitudes,

e_corr_s

The PT2SDd correlation energy wrt single excitations

t_1

The singles amplitudes

All remaining attributes are summarized in RpCCD-LCCD calculations and RHF-LCCD calculations

12.1.5. Troubleshooting in LCC calculations¶

The amplitude solver converges very slowly:

Usually, a calculation should convergence within 10-15 iterations (using the krylov solver). There are, however, cases where the default solver shows poor convergence. Usually this means that the CC model is not appropriate to describe electron correlation in the given system. In case of pCCD-LCC, this can also indicate that the corresponding orbitals are suboptimal and a better solution for the orbitals is to be expected. For instance, the orbital optimizer in pCCD got stuck in some local minimum. You can try to re-optimize the orbitals within pCCD by, for instance, imposing a small perturbation on the present solution.

Note

In any case, the quality of the pCCD solution will strongly affect the performance/convergence of the LCC correction.
The pbqn solver diverges:

The pbqn solver might be unstable for difficult cases. You can switch to the krylov solver in such cases. Note, however, that the latter is rather slow.

12.1.6. \(\Lambda\)-equations¶

The LCC module supports the solution of the corresponding \(\Lambda\)-equations, which are used to, for instance, determine (selected elements of) the 1-, 2-, 3-, 4-particle reduced density matrices (RDM’s). The optimization of the \(\lambda\)-amplitudes is performed by setting the keyword argument l=True in the corresponding LCC function call. The following code snippet shows how to optimize the \(\lambda\)-amplitudes for both pCCD-LCC flavors.

lccd = RpCCDLCCD(lf, occ_model)
lccd_ = lccd(kin, ne, eri, oopccd_, l=True)

lccsd = RpCCDLCCSD(lf, occ_model)
lccsd_ = lccsd(kin, ne, eri, oopccd_, l=True)

The IOData container (return value) contains the corresponding \(\lambda\)-amplitudes as attributes as well as all RDM’s that are required for the orbital entanglement module (see Orbital Entanglement for pCCD-LCC)

Note

We recommend to use the same convergence threshold for the CC amplitude and the \(\Lambda\)-equations. Typically, the default convergence thresholds should be sufficient. If higher accuracy is, however, desired, we recommend to set a tight convergence threshold of 1e-12 (tthreshold=1e-12).

12.1.7. Example Python scripts¶

Several complete examples can be found in the directory data/examples/lcc.

12.1.7.1. RHF-LCC calculation on the water molecule¶

This is a basic example on how to perform a RHF-LCC calculation in PyBEST. This script performs a RHF-LCCD calculation followed by a RHF-LCCSD calculation on the water molecule using the cc-pVDZ basis set.

Listing 12.1 data/examples/lcc/rhflcc_water_cc-pvdz.py¶

from pybest.cc.rlccsd import RHFLCCD, RHFLCCSD
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-LCCD calculation ####################################################
################################################################################
lccd = RHFLCCD(lf, occ_model)
lccd_ = lccd(kin, ne, eri, hf_)

################################################################################
### Do RHF-LCCSD calculation ###################################################
################################################################################
lccsd = RHFLCCSD(lf, occ_model)
lccsd_ = lccsd(kin, ne, eri, hf_)

12.1.7.2. RHF-LCC calculation on the water molecule using a frozen core¶

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

Listing 12.2 data/examples/lcc/rhflcc_frozen_core_water_cc-pvdz.py¶

from pybest.cc.rlccsd import RHFLCCD, RHFLCCSD
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-LCCD calculation and 1 frozen core orbital ##########################
################################################################################
lccd = RHFLCCD(lf, occ_model)
lccd_ = lccd(kin, ne, eri, hf_)

################################################################################
### Do RHF-LCCSD calculation and 1 frozen core orbital #########################
################################################################################
lccsd = RHFLCCSD(lf, occ_model, ncore=1)
lccsd_ = lccsd(kin, ne, eri, hf_)

12.1.7.3. pCCD-LCC calculation on the water molecule¶

This is a basic example on how to perform a pCCD-LCC calculation in PyBEST. This script performs a pCCD-LCCD calculation followed by a pCCD-LCCSD calculation on the water molecule using the cc-pVDZ basis set.

Listing 12.3 data/examples/lcc/rpccdlcc_water_cc-pvdz.py¶

from pybest.cc.rlccsd import RpCCDLCCD, RpCCDLCCSD
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.geminals import RpCCD
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 OO-pCCD optimization ####################################################
###############################################################################
oopccd = RpCCD(lf, occ_model)
oopccd_ = oopccd(kin, ne, eri, hf_)

################################################################################
### Do RpCCD-LCCSD calculation #################################################
################################################################################
lccd = RpCCDLCCD(lf, occ_model)
lccd_ = lccd(kin, ne, eri, oopccd_)

################################################################################
### Do RpCCD-LCCSD calculation #################################################
################################################################################
lccsd = RpCCDLCCSD(lf, occ_model)
lccsd_ = lccsd(kin, ne, eri, oopccd_)

12.1.7.4. pCCD-LCC calculation on the water molecule using a frozen core¶

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

Listing 12.4 data/examples/lcc/rpccdlcc_frozen_core_water_cc-pvdz.py¶

from pybest.cc.rlccsd import RpCCDLCCD, RpCCDLCCSD
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.geminals import RpCCD
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 OO-pCCD optimization and 1 frozen core orbital ##########################
###############################################################################
oopccd = RpCCD(lf, occ_model, ncore=1)
oopccd_ = oopccd(kin, ne, eri, hf_)

################################################################################
### Do RpCCD-LCCSD calculation and 1 frozen core orbital #######################
################################################################################
lccd = RpCCDLCCD(lf, occ_model, ncore=1)
lccd_ = lccd(kin, ne, eri, oopccd_)

################################################################################
### Do RpCCD-LCCSD calculation and 1 frozen core orbital #######################
################################################################################
lccsd = RpCCDLCCSD(lf, occ_model, ncore=1)
lccsd_ = lccsd(kin, ne, eri, oopccd_)