13.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 a single Slater determinant reference function like RHF (see The Self-Consistent Field Module) or an RpCCD (see The pCCD module) reference function. Both modules are explained in greater detail below.

13.1.1. Supported features

The linearized coupled cluster module supports spin-restricted orbitals and the DenseLinalgFactory and CholeskyLinalgFactory.

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

13.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_output (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 = RLCCD(lf, occ_model)
lccd_output = lccd(kin, ne, eri, hf_output)

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

e_tot:

The total RLCCD energy (including all external terms)

e_corr:

The RLCCD 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 RLCCSD class. The overall input structure is equivalent to RLCCD,

lccsd = RLCCSD(lf, occ_model)
lccsd_output = lccsd(kin, ne, eri, hf_output)

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

t_1:

The singles amplitudes

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

13.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 RLCC calculation. We assume that you have performed a restricted pCCD calculation, whose results are stored in the IOData container pccd_output (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_output = lccd(kin, ne, eri, pccd_output)

Instead of creating an instance of RLCCD, 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 RLCCD 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 wave function (here, the total RpCCD energy)

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_output = lccsd(kin, ne, eri, pccd_output)

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.

13.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 RLCC module, where a frozen 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 RLCCSD calculation, specify the ncore argument during the initialization of RLCCSD

# Select 4 frozen core orbital
#-----------------------------
lccsd = RLCCSD(lf, occ_model, ncore=4)
lccsd_output = lccsd(kin, ne, eri, hf_output)

This syntax is working for all LCC modules mentioned above.

13.1.2.4. Restart options

The RLCC module supports restarts from PyBEST’s internal checkpoint files. Furthermore, only one restart option is currently supported

  1. Restart only CC amplitudes (restart keyword)

In order to restart a calculation, you have to specify the restart keyword, which contains the name of the checkpoint file used for restarts. Since the restart behavior is similar in all RCC modules, we select RLCCD as an example.

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

lccd = RLCCD(lf, occ_model)
# Restart only amplitudes from checkpoint_RLCCD.h5
lccd_output = lccd(kin, ne, eri, hf_output, restart='pybest-results/checkpoint_RLCCD.h5')

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

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.

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

13.1.3.1. RLCCD calculations

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

e_corr_d:

The RLCCD correlation energy wrt double excitations

e_corr_s0:

The RLCCD correlation energy wrt the seniority zero sector of all double excitations

e_corr_s2:

The RLCCD correlation energy wrt the seniority 2 sector of all double excitations

e_corr_s4:

The RLCCD correlation energy wrt the seniority 4 sector of all double excitations

orb_a:

A copy of the orbitals used in the RLCCD calculation

olp:

The overlap integrals used in the RLCCD calculation (for restart purposes)

13.1.3.2. RLCCSD (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 RLCCD calculations)

e_corr_s:

The RLCCSD correlation energy wrt single excitations

t_1:

The singles amplitudes

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

restart:

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

initguess:

(str or numpy.ndarray) initial guess for the amplitudes, one of:

random:

random numbers

mp2:

performs MP2 calculations before (default)

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

{“t_1”: DenseTwoIndex instance, …}:

dictionary of amplitudes IOData(t_1=DenseTwoIndex instace), IOData instance

solver:

(str) type of solver used. Possible values are one of scipy.optimize.root solvers 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 100)

threshold:

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

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

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)

diis:

(dictionary) contains optimization parameters for DIIS used in the PBQN solver, one of:

diismax:

(int) maximum number of DIIS vectors (default 2)

diisreset:

(bolean) if True deletes selected DIIS vectors (default True)

diisstart:

(int) first iteration of DIIS procedureue (default 0)

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

13.1.4.1. RpCCD-LCCD calculations

The IOData container of the RpCCDLCCD module has the same attributes as described for the RLCCD method (see RLCCD 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

13.1.4.2. RpCCD-LCCSD calculations

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

e_corr_s:

The RpCCD-LCCSD correlation energy wrt single excitations

t_1:

The singles amplitudes

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

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

13.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_output = lccd(kin, ne, eri, oopccd_output, l=True)

lccsd = RpCCDLCCSD(lf, occ_model)
lccsd_output = lccsd(kin, ne, eri, oopccd_output, 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 (threshold=1e-12).

13.1.7. Example Python scripts

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

13.1.7.1. RLCC calculation on the water molecule

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

Listing 13.1 data/examples/lcc/rhflcc_water_cc-pvdz.py
from pybest.cc import RLCCD, RLCCSD
from pybest.context import context
from pybest.gbasis import (
    get_gobasis,
    compute_kinetic,
    compute_nuclear,
    compute_eri,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.linalg import DenseLinalgFactory
from pybest.scf 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(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_output = hf(kin, ne, eri, external, olp, orb_a)

################################################################################
### Do RHF-LCCD calculation ####################################################
################################################################################
lccd = RLCCD(lf, occ_model)
lccd_output = lccd(kin, ne, eri, hf_output)

################################################################################
### Do RHF-LCCSD calculation ###################################################
################################################################################
lccsd = RLCCSD(lf, occ_model)
lccsd_output = lccsd(kin, ne, eri, hf_output)

13.1.7.2. RLCC calculation on the water molecule using a frozen core

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

Listing 13.2 data/examples/lcc/rhflcc_frozen_core_water_cc-pvdz.py
from pybest.cc import RLCCD, RLCCSD
from pybest.context import context
from pybest.gbasis import (
    get_gobasis,
    compute_kinetic,
    compute_nuclear,
    compute_eri,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.linalg import DenseLinalgFactory
from pybest.scf 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(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_output = hf(kin, ne, eri, external, olp, orb_a)

################################################################################
### Do RHF-LCCD calculation and 1 frozen core orbital ##########################
################################################################################
lccd = RLCCD(lf, occ_model, ncore=1)
lccd_output = lccd(kin, ne, eri, hf_output)

################################################################################
### Do RHF-LCCSD calculation and 1 frozen core orbital #########################
################################################################################
lccsd = RLCCSD(lf, occ_model, ncore=1)
lccsd_output = lccsd(kin, ne, eri, hf_output)

13.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 13.3 data/examples/lcc/rpccdlcc_water_cc-pvdz.py
from pybest.cc import RpCCDLCCD, RpCCDLCCSD
from pybest.context import context
from pybest.gbasis import (
    get_gobasis,
    compute_kinetic,
    compute_nuclear,
    compute_eri,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.geminals import RpCCD
from pybest.linalg import DenseLinalgFactory
from pybest.scf 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(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_output = hf(kin, ne, eri, external, olp, orb_a)

###############################################################################
## Do OO-pCCD optimization ####################################################
###############################################################################
pccd = RpCCD(lf, occ_model)
pccd_output = pccd(kin, ne, eri, hf_output)

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

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

13.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 13.4 data/examples/lcc/rpccdlcc_frozen_core_water_cc-pvdz.py
from pybest.cc import RpCCDLCCD, RpCCDLCCSD
from pybest.context import context
from pybest.gbasis import (
    get_gobasis,
    compute_kinetic,
    compute_nuclear,
    compute_eri,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.geminals import RpCCD
from pybest.linalg import DenseLinalgFactory
from pybest.scf 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(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_output = hf(kin, ne, eri, external, olp, orb_a)

###############################################################################
## Do OO-pCCD optimization and 1 frozen core orbital ##########################
###############################################################################
pccd = RpCCD(lf, occ_model, ncore=1)
pccd_output = pccd(kin, ne, eri, hf_output)

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

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