17.3. Example Python scripts

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

17.3.1. EA-pCCD calculations on the \(\textrm{O}_2\) molecule for 1 and 3 unpaired electrons

This is a basic example of how to perform an EA-pCCD calculation in PyBEST. This script performs an RHF calculation, followed by a pCCD calculation, and finally an EA calculation with a pCCD reference function on the oxygen molecule using the cc-pVDZ basis set. First, only one unpaired electron is considered (doublet and quartet states). In the second calculation, only the quartet states are targeted (3 unpaired electrons)

Listing 17.1 data/examples/ea_eom/reapccd_o2_cc-pvdz.py.py

import numpy as np

from pybest import context
from pybest.ea_eom import REApCCD
from pybest.gbasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
from pybest.geminals import ROOpCCD
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/o2_ccdz.xyz")
basis = get_gobasis("cc-pvdz", fn_xyz)

#
# Define Occupation model, orbitals, and overlap
#
lf = DenseLinalgFactory(basis.nbasis)
orb_a = lf.create_orbital(basis.nbasis)
olp = compute_overlap(basis)
# We need to start with 2 electron less to obtain O_2^+
occ_model = AufbauOccModel(basis, charge=2)

#
# Construct Hamiltonian
#
kin = compute_kinetic(basis)
ne = compute_nuclear(basis)
eri = compute_eri(basis)
external = compute_nuclear_repulsion(basis)

#
# Do a Hartree-Fock calculation
#
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, external, olp, orb_a)
# swap sigma* and pi* orbitals as HF solution has wrong occupation
hf_output.orb_a.swap_orbitals(np.array([[6, 7]]))

#
# Do OO-pCCD optimization
#
oopccd = ROOpCCD(lf, occ_model)
oopccd_output = oopccd(kin, ne, eri, hf_output)

#
# Do REA-pCCD calculation
#
# 1 unpaired electron; O_2^+
ea_1 = REApCCD(lf, occ_model, alpha=1)
ea_output_1 = ea_1(kin, ne, eri, oopccd_output, nroot=4)

# 3 unpaired electrons; O_2^+
ea_3 = REApCCD(lf, occ_model, alpha=3)
ea_output_3 = ea_3(kin, ne, eri, oopccd_output, nroot=3)

17.3.2. DEA-pCCD calculations on the \(\textrm{O}_2\) molecule for 0 and 2 unpaired electrons

This is a basic example of how to perform a DEA-pCCD calculation in PyBEST. This script performs an RHF calculation followed by a pCCD calculation, and finally a DEA calculation with a pCCD reference function on the \(\textrm{O}_2\) molecule using the cc-pVDZ basis set. First, no unpaired electrons are considered (singlet, triplet, and quintet states). In the second calculation, only the triplet and quintet states are targeted (2 unpaired electrons). In the third only quintet states are targeted (3 unpaired electrons).

Listing 17.2 data/examples/ea/rdeapccd_o2_cc-pvdz.py

import numpy as np

from pybest import context
from pybest.ea_eom import RDEApCCD
from pybest.gbasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
from pybest.geminals import ROOpCCD
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/o2_ccdz.xyz")
basis = get_gobasis("cc-pvdz", fn_xyz)

#
# Define Occupation model, orbitals, and overlap
#
lf = DenseLinalgFactory(basis.nbasis)
orb_a = lf.create_orbital(basis.nbasis)
olp = compute_overlap(basis)
# We need to start with 2 electron less
occ_model = AufbauOccModel(basis, charge=2)

#
# Construct Hamiltonian
#
kin = compute_kinetic(basis)
ne = compute_nuclear(basis)
eri = compute_eri(basis)
external = compute_nuclear_repulsion(basis)

#
# Do a Hartree-Fock calculation
#
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, external, olp, orb_a)
# swap sigma* and pi* orbitals as HF solution has wrong occupation
hf_output.orb_a.swap_orbitals(np.array([[6, 7]]))

#
# Do OO-pCCD optimization
#
oopccd = ROOpCCD(lf, occ_model)
oopccd_output = oopccd(kin, ne, eri, hf_output)

#
# Do REA-pCCD calculation
#
# 0 unpaired electron; O_2
ea_1 = RDEApCCD(lf, occ_model, alpha=0)
ea_output_1 = ea_1(kin, ne, eri, oopccd_output, nroot=4)

# 2 unpaired electrons; O_2
ea_2 = RDEApCCD(lf, occ_model, alpha=2)
ea_output_2 = ea_2(kin, ne, eri, oopccd_output, nroot=3)

# 4 unpaired electrons; O_2
ea_3 = RDEApCCD(lf, occ_model, alpha=4)
ea_output_3 = ea_3(kin, ne, eri, oopccd_output, nroot=3)

17.3.3. EA-CCSD calculations on the NO molecule and one unpaired electron

This is a basic example of how to perform an EA-CCSD calculation in PyBEST. This script performs an RHF calculation, followed by a CCSD calculation, and one EA-CCSD calculation on the NO molecule using the cc-pVDZ basis set and targeting one unpaired electron (doublet and quartet states).

Listing 17.3 data/examples/ea_eom/reaccsd_no_cc-pvdz.py

from pybest import context
from pybest.cc import RCCSD
from pybest.ea_eom import REACCSD
from pybest.gbasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
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/no.xyz")
basis = get_gobasis("cc-pvdz", fn_xyz)
#
# Define Occupation model, orbitals, and overlap
#
lf = DenseLinalgFactory(basis.nbasis)
orb_a = lf.create_orbital(basis.nbasis)
olp = compute_overlap(basis)
# we need to remove 1 electron
# If we do not specify the number of frozen core orbitals (ncore),
# then ncore will be calculated automatically
occ_model = AufbauOccModel(basis, charge=1, ncore=0)
#
# Construct Hamiltonian
#
kin = compute_kinetic(basis)
ne = compute_nuclear(basis)
eri = compute_eri(basis)
external = compute_nuclear_repulsion(basis)
#
# Do a Hartree-Fock calculation
#
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, external, olp, orb_a)
#
# Do CCSD calculation
#
ccsd = RCCSD(lf, occ_model)
ccsd_output = ccsd(kin, ne, eri, hf_output)
#
# Do SEACCSD calculation for 1 unpaired electron
#
ea = REACCSD(lf, occ_model, alpha=1)
ea_output = ea(kin, ne, eri, ccsd_output, nroot=8)

17.3.4. EA-LCCD calculations on the NO molecule and one unpaired electron

This is a basic example of how to perform an EA-LCCD calculation in PyBEST. This script performs an RHF calculation, followed by an LCCD calculation, and one EA-LCCD calculation on the NO molecule using the cc-pVDZ basis set and targeting one unpaired electron (doublet and quartet states).

Listing 17.4 data/examples/ea_eom/realccd_no_cc-pvdz.py

from pybest import context
from pybest.cc import RLCCD
from pybest.ea_eom import REALCCD
from pybest.gbasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
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/no.xyz")
basis = get_gobasis("cc-pvdz", fn_xyz)
#
# Define Occupation model, orbitals, and overlap
#
lf = DenseLinalgFactory(basis.nbasis)
orb_a = lf.create_orbital(basis.nbasis)
olp = compute_overlap(basis)
# we need to remove 1 electron
# If we do not specify the number of frozen core orbitals (ncore),
# then ncore will be calculated automatically
occ_model = AufbauOccModel(basis, charge=1, ncore=0)
#
# Construct Hamiltonian
#
kin = compute_kinetic(basis)
ne = compute_nuclear(basis)
eri = compute_eri(basis)
external = compute_nuclear_repulsion(basis)
#
# Do a Hartree-Fock calculation
#
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, external, olp, orb_a)
#
# Do LCCD calculation
#
lccd = RLCCD(lf, occ_model)
lccd_output = lccd(kin, ne, eri, hf_output)
#
# Do REALCCD calculation
#
ea = REALCCD(lf, occ_model, alpha=1)
ea_output = ea(kin, ne, eri, lccd_output, nroot=8)

17.3.5. EA-fpCCSD calculation on the NO molecule and one unpaired electron

This is a basic example of how to perform an EA-fpCCSD calculation in PyBEST. This script performs an RHF calculation, followed by pCCD and fpCCSD calculations, and an EA-fpCCSD calculation on the NO molecule using the cc-pVDZ basis set, targeting one unpaired electron (doublet and quartet states). A similar example for EA-fpCCD is provided in the data/examples/ea_eom/ directory.

Listing 17.5 data/examples/ea_eom/reafpccsd_no_cc-pvdz.py

from pybest import context
from pybest.cc import RfpCCSD
from pybest.ea_eom import REAfpCCSD
from pybest.gbasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
from pybest.geminals import ROOpCCD
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/no.xyz")
basis = get_gobasis("cc-pvdz", fn_xyz)
#
# Define Occupation model, orbitals, and overlap
#
lf = DenseLinalgFactory(basis.nbasis)
orb_a = lf.create_orbital(basis.nbasis)
olp = compute_overlap(basis)
# we need to remove 1 electron
# If we do not specify the number of frozen core orbitals (ncore),
# then ncore will be calculated automatically
occ_model = AufbauOccModel(basis, charge=1, ncore=0)
#
# Construct Hamiltonian
#
kin = compute_kinetic(basis)
ne = compute_nuclear(basis)
eri = compute_eri(basis)
external = compute_nuclear_repulsion(basis)
#
# Do a Hartree-Fock calculation
#
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, external, olp, orb_a)
#
# Do OOpCCD calculation
#
OOpCCD = ROOpCCD(lf, occ_model)
OOpCCD_output = OOpCCD(kin, ne, eri, hf_output)
#
# Do fpCCSD calculation
#
fpccsd = RfpCCSD(lf, occ_model)
fpccsd_output = fpccsd(kin, ne, eri, OOpCCD_output)
#
# Do REAfpCCSD calculation for 1 unpaired electron
#
ea = REAfpCCSD(lf, occ_model, alpha=1)
ea_output = ea(kin, ne, eri, fpccsd_output, nroot=8)

17.3.6. EA-fpLCCSD calculation on the NO molecule and one unpaired electron

This is a basic example of how to perform an EA-fpLCCSD calculation in PyBEST. This script performs an RHF calculation, followed by pCCD and fpLCCSD calculations, and an EA-fpLCCSD calculation on the NO molecule using the cc-pVDZ basis set, targeting one unpaired electron (doublet and quartet states). A similar example for EA-fpCCD is provided in the data/examples/ea_eom/ directory.

Listing 17.6 data/examples/ea_eom/reafplccsd_no_cc-pvdz.py

from pybest import context
from pybest.cc import RpCCDLCCSD
from pybest.ea_eom import REAfpLCCSD
from pybest.gbasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
from pybest.geminals import ROOpCCD
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/no.xyz")
basis = get_gobasis("cc-pvdz", fn_xyz)
#
# Define Occupation model, orbitals, and overlap
#
lf = DenseLinalgFactory(basis.nbasis)
orb_a = lf.create_orbital(basis.nbasis)
olp = compute_overlap(basis)
# we need to remove 1 electron
# If we do not specify the number of frozen core orbitals (ncore),
# then ncore will be calculated automatically
occ_model = AufbauOccModel(basis, charge=1, ncore=0)
#
# Construct Hamiltonian
#
kin = compute_kinetic(basis)
ne = compute_nuclear(basis)
eri = compute_eri(basis)
external = compute_nuclear_repulsion(basis)
#
# Do a Hartree-Fock calculation
#
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, external, olp, orb_a)
#
# Do OOpCCD calculation
#
OOpCCD = ROOpCCD(lf, occ_model)
OOpCCD_output = OOpCCD(kin, ne, eri, hf_output)
#
# Do fpLCCSD calculation
#
fplccsd = RpCCDLCCSD(lf, occ_model)
fplccsd_output = fplccsd(kin, ne, eri, OOpCCD_output)
#
# Do REAfpLCCSD calculation for 1 unpaired electron
#
ea = REAfpLCCSD(lf, occ_model, alpha=1)
ea_output = ea(kin, ne, eri, fplccsd_output, nroot=8)

17.3.7. Cholesky-decomposed EA examples

For computational efficiency, all EA models are also available using Cholesky-decomposed integrals. Below, we show an example for EA-fpCCSD only. Other examples are collected in the example directory.

Listing 17.7 data/examples/ea_eom/reafpccsd_no_cc-pvdz_cholesky.py

from pybest import context
from pybest.cc import RfpCCSD
from pybest.ea_eom import REAfpCCSD
from pybest.gbasis import (
    compute_cholesky_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
from pybest.geminals import ROOpCCD
from pybest.linalg import CholeskyLinalgFactory
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/no.xyz")
basis = get_gobasis("cc-pvdz", fn_xyz)
#
# Define Occupation model, orbitals, and overlap
#
lf = CholeskyLinalgFactory(basis.nbasis)
orb_a = lf.create_orbital(basis.nbasis)
olp = compute_overlap(basis)
# we need to remove 1 electron
# If we do not specify the number of frozen core orbitals (ncore),
# then ncore will be calculated automatically
occ_model = AufbauOccModel(basis, charge=1, ncore=0)
#
# Construct Hamiltonian
#
kin = compute_kinetic(basis)
ne = compute_nuclear(basis)
eri = compute_cholesky_eri(basis)
external = compute_nuclear_repulsion(basis)
#
# Do a Hartree-Fock calculation
#
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, eri, external, olp, orb_a)
#
# Do OOpCCD calculation
#
OOpCCD = ROOpCCD(lf, occ_model)
OOpCCD_output = OOpCCD(kin, ne, eri, hf_output)
#
# Do fpCCSD calculation
#
fpccsd = RfpCCSD(lf, occ_model)
fpccsd_output = fpccsd(kin, ne, eri, OOpCCD_output)
#
# Do REAfpCCSD calculation for 1 unpaired electron
#
ea = REAfpCCSD(lf, occ_model, alpha=1)
ea_output = ea(kin, ne, eri, fpccsd_output, nroot=8)