12.1. The Configuration Interaction 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. PyBEST allows you also to perform Restricted Configuration Interaction calculations on top of any single-reference wave function. 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

Currently, PyBEST supports the following variants of CI calculations

restricted Configuration Interaction Singles (RCIS)

restricted Configuration Interaction Doubles (RCID)

restricted Configuration Interaction Singles and Doubles (RCISD)

For all methods, PyBest supports two techniques of calculations which introduce the wave function as a Slater determinant (SD) or configuration state function (CSF). The CSF is a default.

The details are shown below.

12.1.1. Quick Guide: RCIS/RCID/RCISD

This is a short step-by-step introduction to different CI calculations available in PyBEST. Example input file for the water molecule can be found in the section RCI calculations for the water molecule.

12.1.1.1. RCI on top of RHF

We assume that you have performed a restricted Hartree-Fock calculation, whose results are stored in the IOData container hf_out (see The Self-Consistent Field Module). The code snippet below demonstrates how to perform default RCIS/RCID/RCISD calculations with an RHF reference function in PyBEST,

rcis = RCIS(lf, occ_model)
rcis_out = rcis(kin, ne, eri, hf_out)

rcid = RCID(lf, occ_model)
rcid_out = rcid(kin, ne, eri, hf_out)

rcisd = RCISD(lf, occ_model)
rcisd_out = rcisd(kin, ne, eri, hf_out)

Note

By default, only one root is calculated. That is, for RCIS only the lowest-lying excited state is determined, while for RCID and RCISD only the ground-state root is optimized. See Summary of keyword arguments for more details.

The results are returned as an IOData container and can be found in the ./pybest-results directory, where all data is stored in the checkpoint_CIS.h5, checkpoint_CID.h5, and checkpoint_CISD.h5 files, respectively. Specifically, the IOData container contains the following attributes

e_corr:

The RCIS/RCID/RCISD correlation energy

civ:

The eigenvectors of the RCIS/RCID/RCISD Hamiltonian

e_ref:

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

12.1.1.2. Summary of keyword arguments

The RCIS, RCID, and RCISD modules support various keyword arguments that allow the user to steer the process of the RCIS/RCID/RCISD Hamiltonian diagonalization. In the following, all supported keyword arguments are listed together with their default values. Please note that the default values should be sufficient to reach convergence in simple systems.

nroot:

(int) the number of states (default: 1; ground state for RCID/RCISD; first excited state for RCIS)

nguessv:

(int) total number of guess vectors (default: (nroot-1)*4+1)

tolerance:

(float) tolerance for energies (default: 1e-6)

tolerancev:

(float) tolerance for eigenvectors (default: 1e-5)

maxiter:

(int) maximum number of iterations (default: 200).

davidson:

(boolean) if True, a Davidson diagonalization is performed, otherwise an exact diagonalization is performed (for RCIS only, default: True).

threshold:

(float) printing tolerance for contributions of RCI wave function (default: 0.1).

maxvectors:

(int) maximum number of Davidson vectors (default: nroot*10)

scc:

(boolean) determines whether a size-consistency correction is to be calculated (default: True)

threshold_c_0:

a threshold that helps verifying the accuracy of Davidson-type corrections

print_csf:

(boolean) decides in which variant (True: CSF; False: SD) the results will be printed (default: False)

12.1.1.3. Relation between Configuration State Function and Slater Determinant

A Configuration State Function (CSF) is a symmetry-adapted linear combination of Slater determinants (SD). Below, we illustrate the exact relations between CSF and SD for singly- and doubly-excited configurations. To distinguish between the SD and CSF representation, the individual components of a SD will be denoted with normal size letters, while capitalized letters are used for CSFs.

The relation for single excitations is as follows

(12.1)\[\vert ^{\textrm{A}} _{\textrm{I}}\rangle = \frac{1}{\sqrt{2}} (\vert ^a_i \rangle + \vert ^{\bar{a}} _{\bar{i}} \rangle,\]

while the relation for double excitations is more complicated and can be expressed as a set of equations

(12.2)\[\vert ^{\textrm{A} \textrm{A}} _{\textrm{I}\textrm{I}}\rangle = \vert ^{a \bar{a}} _{i\bar{i}}\rangle\]

(12.3)\[\vert ^{\textrm{A} \textrm{B}} _{\textrm{I} \textrm{I}}\rangle = \frac{1}{\sqrt{2}} (\vert ^{a \bar{b}} _{i\bar{i}} \rangle + \vert ^{b \bar{a}} _{i\bar{i}} \rangle),\]

(12.4)\[\vert ^{\textrm{A} \textrm{A}} _{\textrm{I} \textrm{J}}\rangle = \frac{1}{\sqrt{2}} (\vert ^{a \bar{a}} _{i\bar{j}} \rangle + \vert ^{a \bar{a}} _{j\bar{i}} \rangle),\]

(12.5)\[_{A}\vert ^{\textrm{A} \textrm{B}} _{\textrm{I} \textrm{J}}\rangle = \frac{1}{\sqrt{12}} (\vert ^{ab} _{ij} \rangle + \vert ^{\bar{a} \bar{b}} _{\bar{i}\bar{j}} \rangle + \vert ^{a \bar{b}} _{i\bar{j}} \rangle + \vert ^{b \bar{a}} _{j\bar{i}} \rangle - \vert ^{a \bar{b}} _{j\bar{i}} \rangle - \vert ^{b \bar{a}} _{i\bar{j}} \rangle),\]

(12.6)\[_{B}\vert ^{\textrm{A} \textrm{B}} _{\textrm{I} \textrm{J}}\rangle = \frac{1}{2} (\vert ^{a \bar{b}} _{i\bar{j}} \rangle + \vert ^{b \bar{a}} _{j\bar{i}} \rangle + \vert ^{a \bar{b}} _{j\bar{i}} \rangle + \vert ^{b \bar{a}} _{i\bar{j}} \rangle),\]

12.1.1.4. Setting up calculations using CSFs and SDs

By default, all variants of RCI classes perform a calculation using the CSF representation. To change this and use a SD basis instead, the csf argument has to be set to False in the initialization of an instance of the chosen RCI class. The following code snippet shows how to use this option

rcis = RCIS(lf, occ_model, csf=False)

rcid = RCID(lf, occ_model, csf=False)

rcisd = RCISD(lf, occ_model, csf=False)

12.1.1.5. Frozen core RCI

By default, all (occupied and virtual) orbitals are active. To freeze some (occupied) orbitals, the number of frozen cores has to be specified during the initialization of some occupation module class. The code snippet below shows how to freeze the first (occupied) orbital in a RCIS, RCID, RCISD calculation, by specifying the ncore argument during the initialization of the chosen occupation model

# Select one frozen core orbital
#-------------------------------
occ_model = AufbauOccModel(gobasis, ncore=1)
# Perform CI calculation, ncore is stored in occ_model
#-----------------------------------------------------
rcis = RCIS(lf, occ_model)
rcis_out = rcis(kin, ne, er, hf_out)

rcid = RCID(lf, occ_model)
rcid_out = rcid(kin, ne, er, hf_out)

rcisd = RCISD(lf, occ_model)
rcisd_out = rcisd(kin, ne, er, hf_out)

12.1.1.6. Size-consistency Corrections

The RCI module allows you to calculate Davidson-type corrections for the ground state. The following variants of Davidson-type corrections are supported

Davidson: [davidson-corr]

(12.7)\[E_{DC}=(1-{c_{0}}^2)(E_{RCI} - E_{RF})\]
Renormalized Davidson: [scc-overview]

(12.8)\[E_{RDC}=\left(\frac{1-{c_{0}}^2}{{c_{0}}^2}\right)(E_{RCI} - E_{RF})\]
Modified Pople: [scc-overview]

(12.9)\[E_{PC}=E_{RDC}\left(1-\frac{2}{n_e}\right)\]
Meissner: [meissner-overview]

(12.10)\[E_{MC}=E_{RDC}\left( \frac{(n_e-2)(n_e-3)}{n_e(n_e-1)} \right)\]
Duch and Diercksen: [meissner-overview]

(12.11)\[E_{DDC}=E_{RCI}\left(\frac{1-{c_{0}}^2}{2\left(\frac{n_e-1}{n_e-2}\right)c_0^2-1} \right),\]

where \(E_{RCI}\) indicates the total energy of the RCI method, \(E_{RF}\) is the energy of the reference method, \(c_0\) is the contribution of the reference determinant of the reference method, and \(n_e\) denotes the total number of electrons in the system.

Note

Please note that PyBEST supports size-consistency calculations for two variants of the RCI module: RCID and RCISD. The size-consistency corrections are calculated directly by setting the scc keyword argument to True (see also Summary of keyword arguments). By default, all size-consistency corrections are calculated.

12.2. RCI calculations for the water molecule

This is a basic example on how to perform a CI calculation in PyBEST. This script performs an RCIS, RCID, RCISD calculation on the water molecule using the cc-pVDZ basis set.

Listing 12.1 data/examples/rci/rci_water_cc-pvdz.py

from pybest.ci import RCIS, RCID, RCISD
from pybest.context import context
from pybest.gbasis import (
    get_gobasis,
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model 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(obasis)
orb_a = lf.create_orbital(obasis.nbasis)
olp = compute_overlap(obasis)
###############################################################################
## Construct Hamiltonian ######################################################
###############################################################################
kin = compute_kinetic(obasis)
ne = compute_nuclear(obasis)
er = compute_eri(obasis)
external = compute_nuclear_repulsion(obasis)

###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
hf = RHF(lf, occ_model)
hf_out = hf(kin, ne, er, external, olp, orb_a)

###############################################################################
## Do RHF-CIS calculation using Davidson diagonalization
###############################################################################
rcis = RCIS(lf, occ_model)
rcis_out = rcis(kin, ne, er, hf_out)
###############################################################################
## Do RHF-CID calculation using Davidson diagonalization
###############################################################################
rcid = RCID(lf, occ_model)
rcid_out = rcid(kin, ne, er, hf_out)
###############################################################################
## Do RHF-CISD calculation using Davidson diagonalization
###############################################################################
rcisd = RCISD(lf, occ_model)
rcisd_out = rcisd(kin, ne, er, hf_out)
###############################################################################
## Do RHF-CISD calculation for the ground state and 4 excited states
###############################################################################
rcisd = RCISD(lf, occ_model)
rcisd_out = rcisd(kin, ne, er, hf_out, nroot=5)
###############################################################################
## Do RHF-CIS calculation using SD representation
###############################################################################
rcis = RCIS(lf, occ_model, csf=False)
rcis_out = rcis(kin, ne, er, hf_out)
###############################################################################
## Do RHF-CID calculation using SD representation
###############################################################################
rcid = RCID(lf, occ_model, csf=False)
rcid_out = rcid(kin, ne, er, hf_out)
###############################################################################
## Do RHF-CISD calculation using SD representation
###############################################################################
rcisd = RCISD(lf, occ_model, csf=False)
rcisd_out = rcisd(kin, ne, er, hf_out)