13.4. The restricted tailored 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 single-reference tailored coupled-cluster correction on top of a CASCI-type wave function. The module is explained in greater detail below.

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

13.4.1. Quick Guide: RtCC

This is a short introduction to the RtCC module. More information regarding input and output structures 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 electron repulsion integrals

13.4.1.1. \(e^{T_{CAS}}\) reference function in the canonical RHF orbital basis

Currently, CASCI-type calculations are not available in PyBEST and must be performed with external software, e.g., the QC-DMRG-Budapest program. To ensure the correctness of calculations, the one- and two-body Hamiltonian integrals must maintain the same order throughout all stages of the computations. Thus, we suggest either generating integrals in the FCIDUMP format using some external software or generating integrals using PyBEST and read them in the external software. The following snippet allows us to save one- and two-body integrals to file. The integrals from the cas.FCIDUMP file can be used for CASCI-type calculations with some external software.

from pybest.utility import transform_integrals, split_core_active
# Transform Hamiltonian to MO basis
mo_ints = transform_integrals(kin, ne, eri, hf_output.orb_a)
one = mo_ints.one[0]
two = mo_ints.two[0]
# Write all integrals to a FCIDUMP file
data = IOData(one=one, two=two, e_core=e_core, ms2=0, nelec=14)
data.to_file("all.FCIDUMP")
# Get CAS Hamiltonian
mo_ints_cas = split_core_active(one, two, ncore=4, nactive=6, e_core=e_core)
one_cas = mo_ints_cas.one
two_cas = mo_ints_cas.two
# Write CAS integrals to a FCIDUMP file
data = IOData(one=one_cas, two=two_cas, e_core=e_core, ms2=0, nelec=6)
data.to_file("cas.FCIDUMP")

See FCIDUMP format for more details about dumping integrals.

We assume that you have performed some CASCI calculations using cas.FCIDUMP with some external software, and you obtained a file with the CC amplitudes reconstructed from CASCI result called T_DUMP. One example of such a file (obtained with the QC-DMRG-Budapest program) is available in data/examples/dmrg_tccsd/data/T_DUMP.

The following code snippet below shows how to perform an RtCCSD calculation with a DMRG reference function in PyBEST.

# Read integrals
data = IOData.from_file("all.FCIDUMP")
# Do tailored coupled cluster calculations
tcc = RtCCSD(lf, occ_model)
tcc_out = tcc(
    data.one, data.two, data.orb_a,
    external_file="T_DUMP",
    e_ref=energy_of_reference_determinant,
)

Note

The order of arguments does not matter 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_RtCCSD.h5 file. Specifically, the IOData container contains (amongst others) the following attributes

e_tot:

the total RtCCSD energy (sum of reference energy and RCC energy)

e_corr:

the RtCCSD correlation energy

e_corr_s:

the single excitation contribution in RtCCSD correlation energy

e_corr_d:

the double excitation contribution in RtCCSD correlation energy

e_ref:

the energy of the reference determinant

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 RHF model) orbitals

ncore:

the number of (frozen) core orbitals (occupied but not correlated)

nact:

the number of correlated occupied and unoccupied orbitals

converged:

a boolean value indicating if CC converged

13.4.1.2. \(e^{T_{CAS}}\) reference function in ROOpCCD orbital basis

To perform analogous calculations in the ROOpCCD basis, follow the steps from e^{T_{CAS}} reference function in the canonical RHF orbital basis but use the ROOpCCD output instead of the RHF output while transforming the integrals.

13.4.1.3. Summary of keyword arguments

The RtCCSD and RtCCSD 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.

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

external_file:

(str) path to file containing external amplitudes.

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_r:

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

threshold_e:

(float) tolerance for energy (default: 1e-6, works only in combination with “pbqn” solver).

13.4.1.4. RtCCSD calculation on the carbon nitride cation

This is a basic example on how to perform a RtCCSD calculation in PyBEST. These scripts perform a RtCCSD calculation on the carbon nitride action using the CC amplitudes from DMRG calculations.

In the 0th step, we generate one- and two-body integrals in ROOpCCD orbital basis.

Listing 13.9 data/examples/dmrg_tccsd/step0.py

"""Getting one- and two-body integrals in the ROOpCCD orbital basis for
the CN+ molecule.
"""

from pybest.gbasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
from pybest.geminals import ROOpCCD
from pybest.io import IOData
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.utility import split_core_active, transform_integrals
from pybest.wrappers import RHF

###############################################################################
## Set up molecule, define basis set ##########################################
###############################################################################
# get the XYZ file from PyBEST's test data directory
obasis = get_gobasis("cc-pvdz", "data/cn+.xyz")
lf = DenseLinalgFactory(obasis.nbasis)
occ_model = AufbauOccModel(obasis, charge=1)
###############################################################################
## Construct Hamiltonian ######################################################
###############################################################################
kin = compute_kinetic(obasis)
ne = compute_nuclear(obasis)
eri = compute_eri(obasis)
e_core = compute_nuclear_repulsion(obasis)
olp = compute_overlap(obasis)
orb_a = lf.create_orbital()
###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, e_core, olp, orb_a)
###############################################################################
# Do pCCD using RHF orbitals as an inital guess for solver ####################
###############################################################################
pccd = ROOpCCD(lf, occ_model)
pccd_out = pccd(kin, ne, eri, olp, hf_output)
###############################################################################
# Transform Hamiltonian to MO basis
###############################################################################
mo_ints = transform_integrals(kin, ne, eri, hf_output.orb_a)
one = mo_ints.one[0]
two = mo_ints.two[0]
###############################################################################
# Write all integrals to a FCIDUMP file
###############################################################################
data = IOData(one=one, two=two, e_core=e_core, ms2=0, nelec=12)
data.to_file("all.FCIDUMP")
###############################################################################
# Get CAS Hamiltonian
###############################################################################
mo_ints_cas = split_core_active(one, two, ncore=2, nactive=8, e_core=e_core)
one_cas = mo_ints_cas.one
two_cas = mo_ints_cas.two
###############################################################################
# Write CAS integrals to a FCIDUMP file
###############################################################################
data = IOData(one=one_cas, two=two_cas, e_core=e_core, ms2=0, nelec=8)
data.to_file("cas.FCIDUMP")

The generated complete active space integrals (saved as cas.FCIDUMP) can be utilized to compute a matrix product state (MPS) wave function. We extract CC amplitudes from the MPS and save them as data/T_DUMP. In the second step, we perform DMRG-tCCSD calculations based on the data obtained from the previous steps.

Listing 13.10 data/examples/dmrg_tccsd/step2.py

"""DMRG-tailored CCSD calculations performed in the pCCD orbital basis for the
CN+ molecule. DMRG amplitudes have been precomputed using the Budapest DMRG
program.
"""

from pybest.cc import RtCCSD
from pybest.gbasis import get_gobasis
from pybest.io import IOData
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel

###############################################################################
## Set up molecule, define basis set ##########################################
###############################################################################
# get the XYZ file from PyBEST's test data directory
obasis = get_gobasis("cc-pvdz", "data/cn+.xyz")
lf = DenseLinalgFactory(obasis.nbasis)
occ_model = AufbauOccModel(obasis, charge=1)
###############################################################################
# Load all integrals to a FCIDUMP file
###############################################################################
ints = IOData.from_file("all.FCIDUMP")
pccd = IOData.from_file("pybest-results/checkpoint_pccd.h5")
###############################################################################
### Do DMRG-tailored CCSD calculation #########################################
###############################################################################
tcc = RtCCSD(lf, occ_model)
tcc_out = tcc(
    ints.one,
    ints.two,
    ints.orb_a,
    e_ref=pccd.e_ref,
    external_file="data/T_DUMP",
)