15.3. Example Python scripts
Several complete examples can be found in the directory data/examples/ip.
15.3.1. IP-pCCD calculations on the NO molecule for 1 and 3 unpaired electrons
This is a basic example of how to perform an IP-pCCD calculation in PyBEST. This script performs an RHF calculation followed by a pCCD calculation and two different IP flavors with a pCCD reference function on the NO 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)
from pybest import context
from pybest.gbasis import (
compute_eri,
compute_kinetic,
compute_nuclear,
compute_nuclear_repulsion,
compute_overlap,
get_gobasis,
)
from pybest.geminals import ROOpCCD
from pybest.ip_eom import RIPpCCD
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")
obasis = get_gobasis("cc-pvdz", fn_xyz)
###############################################################################
## Define Occupation model, expansion coefficients, and overlap ###############
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
orb_a = lf.create_orbital(obasis.nbasis)
olp = compute_overlap(obasis)
# we need to add 1 additional electron
# If we do not specify the number of frozen core orbitals (ncore),
# then ncore will be calculated automatically
occ_model = AufbauOccModel(obasis, charge=-1, ncore=0)
###############################################################################
## 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_output = hf(kin, ne, eri, external, olp, orb_a)
###############################################################################
## Do OO-pCCD optimization ####################################################
###############################################################################
oopccd = ROOpCCD(lf, occ_model)
oopccd_output = oopccd(kin, ne, eri, hf_output)
###############################################################################
### Do RIP-pCCD calculation for 1 unpaired electron ###########################
###############################################################################
ip = RIPpCCD(lf, occ_model, alpha=1)
ip_output = ip(kin, ne, eri, oopccd_output, nroot=4)
###############################################################################
### Do RIP-pCCD calculation for 3 unpaired electrons ##########################
###############################################################################
ip = RIPpCCD(lf, occ_model, alpha=3)
ip_output = ip(kin, ne, eri, oopccd_output, nroot=4)
15.3.2. DIP-pCCD calculations on the \(\textrm{O}_2\) molecule for 0 and 2 unpaired electrons
This is a basic example of how to perform a DIP-pCCD calculation in PyBEST. This script performs an RHF calculation followed by a pCCD calculation and two different IP flavors 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)
from pybest import context
from pybest.gbasis import (
compute_eri,
compute_kinetic,
compute_nuclear,
compute_nuclear_repulsion,
compute_overlap,
get_gobasis,
)
from pybest.geminals import RpCCD
from pybest.ip_eom import RDIPpCCD
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")
obasis = get_gobasis("cc-pvdz", fn_xyz)
###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
orb_a = lf.create_orbital(obasis.nbasis)
olp = compute_overlap(obasis)
# we need to add 2 additional electrons
# If we do not specify the number of frozen core orbitals (ncore),
# then ncore will be calculated automatically
occ_model = AufbauOccModel(obasis, charge=-2, ncore=0)
###############################################################################
## 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_output = hf(kin, ne, eri, external, olp, orb_a)
###############################################################################
## Do OO-pCCD optimization ####################################################
###############################################################################
oopccd = RpCCD(lf, occ_model)
oopccd_output = oopccd(kin, ne, eri, hf_output)
###############################################################################
### Do RDIP-pCCD calculation for 0 unpaired electrons #########################
### The lowest-lying states corresponding to 3Sigma, 1Delta, 1Sigma ###########
###############################################################################
ip = RDIPpCCD(lf, occ_model, alpha=0)
ip_output = ip(kin, ne, eri, oopccd_output, nroot=4)
###############################################################################
### Do RDIP-pCCD calculation for 2 unpaired electrons (high-spin DIP model) ###
### The lowest-lying state corresponding to 3Sigma ############################
###############################################################################
ip = RDIPpCCD(lf, occ_model, alpha=2)
ip_output = ip(kin, ne, eri, oopccd_output, nroot=1)
15.3.3. IP-CCSD calculations on the NO molecule for 1 unpaired electron
This is a basic example of how to perform an IP-CCSD calculation in PyBEST. This script performs an RHF calculation followed by a CCSD calculation and an IP-EOM calculation with a CCSD reference function on the NO molecule using the cc-pVDZ basis set.
data/examples/ip/rip_ccsd_no_cc-pvdz.py
from pybest import context
from pybest.cc import RCCSD
from pybest.gbasis import (
compute_eri,
compute_kinetic,
compute_nuclear,
compute_nuclear_repulsion,
compute_overlap,
get_gobasis,
)
from pybest.ip_eom import RIPCCSD
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")
obasis = get_gobasis("cc-pvdz", fn_xyz)
###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
orb_a = lf.create_orbital(obasis.nbasis)
olp = compute_overlap(obasis)
# we need to add 1 additional electron
# If we do not specify the number of frozen core orbitals (ncore),
# then ncore will be calculated automatically
occ_model = AufbauOccModel(obasis, charge=-1, ncore=0)
###############################################################################
## 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_output = hf(kin, ne, eri, external, olp, orb_a)
###############################################################################
## Do CCSD optimization #######################################################
###############################################################################
ccsd = RCCSD(lf, occ_model)
ccsd_output = ccsd(kin, ne, eri, hf_output)
###############################################################################
### Do RIP-CCSD calculation for 1 unpaired electron ###########################
###############################################################################
ip = RIPCCSD(lf, occ_model, alpha=1)
ip_output = ip(kin, ne, eri, ccsd_output, nroot=8, nhole=2)
15.3.4. IP-fpCCSD calculations on the NO molecule for 1 unpaired electron
This is a basic example of how to perform an IP-fpCCSD calculation in PyBEST. This script performs an RHF calculation followed by a OOpCCD and subsequent fpCCSD calculation and an IP-EOM calculation with a fpCCSD reference function on the NO molecule using the cc-pVDZ basis set.
data/examples/ip/rip_fpccsd_no_cc-pvdz.py
from pybest import context
from pybest.cc import RfpCCSD
from pybest.gbasis import (
compute_eri,
compute_kinetic,
compute_nuclear,
compute_nuclear_repulsion,
compute_overlap,
get_gobasis,
)
from pybest.geminals import ROOpCCD
from pybest.ip_eom import RIPfpCCSD
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")
obasis = get_gobasis("cc-pvdz", fn_xyz)
###############################################################################
## Define Occupation model, expansion coefficients, and overlap ###############
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
orb_a = lf.create_orbital(obasis.nbasis)
olp = compute_overlap(obasis)
# we need to add 1 additional electron
# If we do not specify the number of frozen core orbitals (ncore),
# then ncore will be calculated automatically
occ_model = AufbauOccModel(obasis, charge=-1)
###############################################################################
## 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_output = hf(kin, ne, eri, external, olp, orb_a)
###############################################################################
## Do OO-pCCD optimization ####################################################
###############################################################################
oopccd = ROOpCCD(lf, occ_model)
oopccd_output = oopccd(kin, ne, eri, hf_output)
###############################################################################
## Do fpCCSD optimization #####################################################
###############################################################################
fpccsd = RfpCCSD(lf, occ_model)
fpccsd_output = fpccsd(kin, ne, eri, oopccd_output)
###############################################################################
### Do RIP-fpCCSD calculation for 1 unpaired electron #########################
###############################################################################
ip = RIPfpCCSD(lf, occ_model, alpha=1)
ip_output = ip(kin, ne, eri, fpccsd_output, nroot=8, nhole=2)