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 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 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 be stored in the same order at all stages of the calculations. 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.orbital_utils 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

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

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

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

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

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",
)