13.2. The EOM module with a pCCD reference function¶

Please, first read the Quick Guide on EOM-pCCD available here.

13.2.1. EOM-pCCD and EOM-pCCD+S¶

If you use this module, please cite [boguslawski2016a] and [boguslawski2017c].

In addition to the IOData container attributes mentioned in the Quick Guide above, the REOMpCCD and REOMpCCDS containers include the following information

orb_a

A copy of the orbitals used in the pCCD reference calculation

olp

The overlap integrals used in the pCCD calculation

Due to their simplicity, the effective Hamiltonian of any EOM-pCCD model can be explicitly constructed and exactly diagonalized. If you wish to solve for the eigenvalues and eigenvalues iteratively, you can use the Davidson diagonalizer using the corresponding davidson keyword argument. For instance, if you target the excited states of large systems, the exact diagonalization might be prohibitive. For very large systems, we thus recommend to keep the default setting and exploit the Davidson diagonalizer,

################################################################################
eom = REOMpCCD(lf, occ_model)

If you, however, want to use exact diagonalization, you can do this by setting the davidson keyword argument to False.

Note

The davidson keyword argument works for all REOCC modules in PyBEST.

13.2.2. Summary of keyword arguments¶

The REOCC module supports various keyword arguments that allow us to steer the optimization of excited states (excitation energies and the eigenvectors of the targeted states). 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.

indextrans

(str) 4-index Transformation. The 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

nroot

(int) the number of targeted states, including the ground state (default 2)

davidson

(boolean) type of solver used. If set to True PyBEST’s internal Davidson diagonalizer is used, where only the lowest nroot states are optimized. If set to False, only the lowest nroot are returned and printed (default True)

Note

For very large systems switching davidson to True might save memory. For small- and medium-sized systems, an exact diagonalization is typically of similar cost due to the simplicity of the model.

threshold

(float) printing threshold for the eigenvectors. If a (row) element of civ_ee is larger than the threshold in absolute value, the corresponding excitation contribution is printed (default 0.1)

tco

(str) the flavor of tensor contraction operations used. Either td or einsum (default td)

Note

tco has only effect if multiple operations are supported. Otherwise, PyBEST always refers to einsum.

Some keyword arguments are only working together with the davidson solver. If exact diagonalization is used, they have no effect.

tolerance

(float) optimization threshold for each excitation energy (default 1e-6)

tolerancev

(float) optimization threshold for each excited state eigenvector (default 1e-4)

maxiter

(int) maximum number of total Davidson diagonalization steps (default 200)

nguessv

(int) total number of guess vectors (default (nroots-1)*4+1, that is, 4 vectors per excited state)

maxvectors

(int) maximum number of Davidson vectors. If additional vectors are added, a subspace collapse is performed (default (nroots-1)*10, that is, 10 vectors per excited state)

todisk

(boolean) if set to True all intermediates are stored to disk. This reduces memory requirements. However, due to the intensive I/O operations, the code is slowed downed significantly (default False).

13.2.3. Example Python scripts¶

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

13.2.3.1. EOM-pCCD and EOM-pCCD+S calculations on the water molecule¶

This is a basic example of how to perform a EOM-pCCD and EOM-pCCD+S calculation in PyBEST. This script performs a RHF calculation followed by a pCCD calculation and two different EOM flavors with a pCCD reference function on the water molecule using the cc-pVDZ basis set.

Listing 13.1 data/examples/eom/reompccd_water_cc-pvdz.py¶

from pybest import context
from pybest.ee_eom.eom_pccd import REOMpCCD
from pybest.gbasis import get_gobasis
from pybest.gbasis.gobasis import (
    compute_dipole,
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.geminals.rpccd import RpCCD
from pybest.linalg.dense import DenseLinalgFactory
from pybest.ee_eom.eom_pccd import REOMpCCD
from pybest.ee_eom.eom_pccd_s import REOMpCCDS
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 REOM-pCCD calculation ###################################################
################################################################################
eom = REOMpCCD(lf, occ_model)
eom_ = eom(kin, ne, eri, oopccd_, nroot=3)

################################################################################
### Do REOM-pCCD+S calculation #################################################
################################################################################
eom = REOMpCCDS(lf, occ_model)
eom_ = eom(kin, ne, eri, oopccd_, nroot=3)

################################################################################
### Do REOM-pCCD calculation davidson              #############################
################################################################################
eom = REOMpCCD(lf, occ_model)
eom_ = eom(kin, ne, eri, oopccd_, nroot=3, davidson=True)

13.2.3.2. EOM-pCCD and EOM-pCCD+S calculations on the water molecule using a frozen core¶

This is a basic example on how to perform a EOM-pCCD and EOM-pCCD+S calculation in PyBEST. This script performs a RHF calculation followed by a pCCD calculation and two different EOM flavors with a pCCD reference function on the water molecule using the cc-pVDZ basis set and freezing the lowest orbital.

Listing 13.2 data/examples/eom/reompccd_core_water_cc-pvdz.py¶

from pybest import context
from pybest.ee_eom.eom_pccd import REOMpCCD
from pybest.gbasis import get_gobasis
from pybest.gbasis.gobasis import (
    compute_dipole,
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.geminals.rpccd import RpCCD
from pybest.linalg.dense import DenseLinalgFactory
from pybest.ee_eom.eom_pccd import REOMpCCD
from pybest.ee_eom.eom_pccd_s import REOMpCCDS
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, ncore=1)
oopccd_ = oopccd(kin, ne, eri, hf_)

################################################################################
### Do REOM-pCCD calculation ###################################################
################################################################################
eom = REOMpCCD(lf, occ_model, ncore=1)
eom_ = eom(kin, ne, eri, oopccd_, nroot=3)

################################################################################
### Do REOM-pCCD+S calculation #################################################
################################################################################
eom = REOMpCCDS(lf, occ_model, ncore=1)
eom_ = eom(kin, ne, eri, oopccd_, nroot=3)

################################################################################
### Do REOM-pCCD calculation using exact diagonalization #######################
################################################################################
eom = REOMpCCD(lf, occ_model, ncore=1)
eom_ = eom(kin, ne, eri, oopccd_, nroot=3, davidson=False)