15.1. Quick Guide

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. The current version of PyBEST offers ionization potential (IP) calculations with a RpCCD (see The pCCD module) reference function using the Equation of Motion (EOM) formalism. The IP module is explained in greater detail below.

15.1.1. Supported features

The IP module supports spin-restricted orbitals and the DenseLinalgFactory and CholeskyLinalgFactory.

15.1.2. How to: RIP

This is a short introduction to the RIP 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 two-electron repulsion integrals

Note

In all IP modules, only the spin projection is defined, that is, only the number of excess \(\alpha\) electrons over \(\beta\) electrons is defined.

15.1.2.1. RpCCD reference function

This version of PyBEST supports single IP-pCCD and double (D)IP-pCCD calculations with two different excitation operators

  • IP-pCCD (RIPpCCD):

    1. 1 hole (1h) operators

    2. 2 hole, 1 particle (2h1p) operators

  • DIP-pCCD (RDIPpCCD):

    1. 2 hole (2h) operators

    2. 3 hole, 1 particle (3h1p) operators

Complete examples can be found in the following subsections (Example Python scripts).

15.1.2.1.1. IP-pCCD: Doublet and quartet states

If you use this module, please cite [boguslawski2021].

We assume that you have performed a restricted pCCD calculation (either with or without orbital optimization), whose results are stored in the IOData container pccd_output (see The pCCD module). The code snippet below shows how to perform a RIPpCCD calculation. The current version supports the optimization of doublet and quartet states and just quartet states. The former can be optimized by setting the keyword argument alpha to 1 when creating an instance of RIPpCCD, while the latter are accessible by choosing alpha=3. The number of targeted states is passed through the nroot keyword argument in the function call.

# Calculate 3 lowest-lying roots for 1 unpaired electron (S_z=0.5)
# Both doublet and quartet states can be targeted
ip = RIPpCCD(lf, occ_model, alpha=1)
ip_output = ip(kin, ne, eri, pccd_output, nroot=3)

The results are returned as a IOData container, while all results are saved to the pybest-results/checkpoint_IP-EOM-pCCD.h5 file. Specifically, the IOData container contains (amongst others) the following attributes

e_ip:

The ionization energies in \(E_h\)

civ_ip_alpha:

The CI vectors (that is, the eigenvectors) for each state (column) for a given alpha value

The eigenvalues and eigenvectors are stored as numpy arrays, not as instances of the LinalgFactory. Since for one spin projection \(S_z\), several spin multiplets are accessible, PyBEST suggests the multiplicity for each targeted state in the output file.

By default, RIPpCCD includes up to 2h1p terms. For alpha=1, the 2h1p terms can be neglected and only the 1h terms are considered during the diagonalization. The number of hole/particle operators is specified by the keyword argument nhole. The following code snippet shows how to perform a RIPpCCD calculation for 1 unpaired electron including only 1h terms,

# Calculate 3 lowest-lying roots for 1 unpaired electron (S_z=0.5) and 1h terms
# Both doublet and quartet states can be targeted
ip = RIPpCCD(lf, occ_model, alpha=1)
ip_output = ip(kin, ne, eri, pccd_output, nroot=3, nhole=1)

Note

For three unpaired electrons, only 2h1p terms are supported. Specifying nhole=1 will raise an error.

15.1.2.1.2. DIP-pCCD: Singlet, triplet, and quintet states

If you use this module, please cite [boguslawski2021].

We assume that you have performed a restricted pCCD calculation (either with or without orbital optimization), whose results are stored in the IOData container pccd_output (see The pCCD module). The code snippet below shows how to perform a RDIPpCCD calculation. The current version supports the optimization of singlet, triplet, and quintet states, triplet and quintet states, or just quintet states. The first case can be optimized by setting the keyword argument alpha to 0 when creating an instance of RDIPpCCD, the second case is accessible by choosing alpha=2, while the last spin multiplicity can be chosen by setting alpha=4. The number of targeted states is passed through the nroot keyword argument in the function call.

# Calculate 3 lowest-lying roots for 0 unpaired electron (S_z=0.0)
# Singlet, triplet, and quintet states can be targeted
ip = RDIPpCCD(lf, occ_model, alpha=0)
ip_output = ip(kin, ne, eri, pccd_output, nroot=3)

The results are returned as a IOData container, while all results are saved to the pybest-results/checkpoint_DIP-EOM-pCCD.h5 file. Specifically, the IOData container contains (amongst others) the following attributes

e_ip:

The ionization energies in \(E_h\)

civ_ip_alpha:

The CI vectors (that is, the eigenvectors) for each state (column) for a given alpha value

The eigenvalues and eigenvectors are stored as numpy arrays, not as instances of the LinalgFactory. Since for one spin projection \(S_z\), several spin multiplets are accessible, PyBEST suggests the multiplicity for each targeted state in the output file.

By default, RDIPpCCD includes up to 3h1p terms. For alpha=0 and alpha=2, the 3h1p terms can be neglected so that only the 2h terms are considered during the diagonalization. The number of hole/particle operators is specified by the keyword argument nhole. The following code snippet shows how to perform a RDIPpCCD calculation for 0 unpaired electron including only 2h terms,

# Calculate 3 lowest-lying roots for 1 unpaired electron (S_z=0.0) and 2h terms
# Singlet and triplet states can be targeted
ip = RDIPpCCD(lf, occ_model, alpha=0)
ip_output = ip(kin, ne, eri, pccd_output, nroot=3, nhole=2)

Note

For four unpaired electrons (qunitet states), only 3h1p terms are supported. Specifying nhole=2 will raise an error.

15.1.2.2. Defining a frozen core

By default, all (occupied and virtual) orbitals are active. If a frozen core has been selected in the CC reference calculation, the same orbitals have to be frozen in the chosen IP flavor. To freeze some (occupied) orbitals, the number of frozen cores has to be specified when an instance of some OccupationModel class is created. For instance, if you wish to freeze the first 4 (occupied) orbitals in an IP-pCCD calculation, specify the ncore argument during the initialization of some occupation model class,

# Select 4 frozen core orbital
#-----------------------------
occ_model = AufbauOccModel(gobasis, ncore=4)
# Perform IP-CC calculation for 1 unpaired electron (S_z=0.5),
# ncore is stored in occ_model
#-------------------------------------------------------------
ip = RIPpCCD(lf, occ_model, alpha=1)
ip_output = ip(kin, ne, eri, pccd_output)

This syntax is working for all IP modules mentioned above.

15.1.2.3. Restart options

Restart options are not supported yet.