19.5. Orbital Energies

This section documents how to compute orbital energies for closed-shell molecules using Koopmans-based methods. The Koopmans’ theorem approximates ionization energies and electron affinities by the orbital energies of a system.

Before proceeding, please ensure you’re familiar with the basics of pCCD models in PyBEST (The pCCD module).

19.5.1. Supported Models

Koopmans’ theorem (Koopmans’ Theorem): Calculate and store single and pair orbital energies from Koopmans’ theorem.
Modified Koopmans’ theorem (Modified Koopmans’ Theorem): Calculate and store single and pair orbital energies from Modified Koopmans’ theorem.
Extended Koopmans’ theorem: Calculate and store orbital energies from Extended Koopmans’ theorem (unreleased).

Note

The Koopmans theorem is implemented for a restricted Hartree–Fock or pCCD reference (The pCCD module) to calculate orbital energies. The modified version only supports a pCCD reference.

19.5.2. Theoretical Framework

19.5.3. Koopmans’ Theorem

Koopmans’ theorem provides an approximate link between orbital energies and electron addition/removal processes within the RHF framework. The theorem states that the negative of an occupied orbital’s energy approximates the ionization potential (IP), while that of a virtual orbital estimates the electron affinity (EA). These expressions derive from differences in total energies when electrons are added or removed, and are evaluated from the diagonal elements of the Fock matrix. Extensions to double ionization potentials (DIPs) and double electron affinities (DEAs) are also possible, including electron–electron repulsion terms.

19.5.4. Modified Koopmans’ Theorem

The modified Koopmans’ theorem extends the original theorem by incorporating correlation effects through a multi-determinant pCCD reference. It replaces the bare Hamiltonian with a similarity-transformed version, \(\hat{H}^{(\mathrm{pCCD})} = e^{-\hat{T}_p} \hat{H} e^{\hat{T}_p}\), where \(\hat{T}_p\) is the pCCD cluster operator. IPs and EAs are then evaluated as commutators and expectation values using this transformed Hamiltonian, capturing additional correlation contributions localized to specific orbitals. These corrections refine the orbital energies beyond the Hartree–Fock approximation.

For example, the modified IP expression becomes \(-\!f_{ii} - \sum_c t^{c\bar{c}}_{i\bar{i}} \langle i\bar{i} | c\bar{c} \rangle\), explicitly accounting for the correlation energy lost when removing an electron from orbital \(i\). The corresponding EA is defined analogously with respect to a virtual orbital.

19.5.5. Summary of Koopmans’ and Modified Koopmans’ equations for IPs, EAs, DIPs, and DEAs

Table 19.1 Summary of Koopmans’ and Modified Koopmans’ Equations for IPs, EAs, DIPs, and DEAs
Orbital Energy	Koopmans’	Modified Koopmans’
IP	\(f_{ii}\)	\(-f_{ii} - \sum\limits_{c} t^{c\bar{c}}_{i\bar{i}} \langle i\bar{i} \vert c\bar{c} \rangle\)
EA	\(-f_{aa}\)	\(-f_{aa} + \sum\limits_{k} t^{a\bar{a}}_{k\bar{k}} \langle k\bar{k} \vert a\bar{a} \rangle\)
DIP (\(m_s=1\))	\(f_{ii} + f_{jj} - \langle ij \Vert ij \rangle\)	\(-f_{ii} - f_{jj} + V_{ijij}\) \(- \sum\limits_{c} t^{c\bar{c}}_{i\bar{i}} \langle i\bar{i} \vert c\bar{c} \rangle\) \(- \sum\limits_{c} t^{c\bar{c}}_{j\bar{j}} \langle j\bar{j} \vert c\bar{c} \rangle\)
DIP (\(m_s=0\))	\(f_{ii} + f_{\bar{j}\bar{j}}\) \(- \langle i\bar{j} \vert i\bar{j} \rangle\)	\(-f_{ii} - f_{\bar{j}\bar{j}} + V_{i\bar{j}i\bar{j}}\) \(- \sum\limits_{c} t^{c\bar{c}}_{i\bar{i}} \langle i\bar{i} \vert c\bar{c} \rangle\) \(- \sum\limits_{c} t^{c\bar{c}}_{j\bar{j}} \langle j\bar{j} \vert c\bar{c} \rangle\) \(+ \sum\limits_{c} t^{c\bar{c}}_{i\bar{i}} \langle i\bar{i} \vert c\bar{c} \rangle \delta_{ij}\)
DEA (\(m_s=1\))	\(f_{aa} + f_{bb} + \langle ab \Vert ab \rangle\)	\(-f_{aa} - f_{bb} - V_{abab}\) \(+ \sum\limits_{k} t^{a\bar{a}}_{k\bar{k}} \langle k\bar{k} \vert a\bar{a} \rangle\) \(+ \sum\limits_{k} t^{b\bar{b}}_{k\bar{k}} \langle k\bar{k} \vert b\bar{b} \rangle\)
DEA (\(m_s=0\))	\(f_{aa} + f_{\bar{b}\bar{b}}\) \(+ \langle a\bar{b} \vert a\bar{b} \rangle\)	\(-f_{aa} - f_{\bar{b}\bar{b}} - V_{a\bar{b}a\bar{b}}\) \(+ \sum\limits_{k} t^{a\bar{a}}_{k\bar{k}} \langle k\bar{k} \vert a\bar{a} \rangle\) \(+ \sum\limits_{k} t^{b\bar{b}}_{k\bar{k}} \langle k\bar{k} \vert b\bar{b} \rangle\) \(- \sum\limits_{k} t^{a\bar{a}}_{k\bar{k}} \langle k\bar{k} \vert a\bar{a} \rangle \delta_{ab}\)

19.5.6. Quick Guide: Koopmans-based orbital energies

If you use this part of the code, please cite [jahani2025a] and [jahani2025b].

PyBEST allows the user to calculate orbital energies using various flavors of Koopmans theorem. The current version of PyBEST supports the Koopmans-based methods for a closed-shell RHF (see The Self-Consistent Field Module) or pCCD wavefunction (see The pCCD module). We assume you have performed an RHF or pCCD calculation whose results are stored in some output container. Furthermore, 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

19.5.6.1. How to: Koopmans’ orbital energies

The code snippet below shows how to perform a Koopmans’ orbital energy calculation in PyBEST,

orbital_energy = Koopmans(lf, occ_model)
e_orb_output = orbital_energy(kin, ne, eri, pccd_output)

where pccd_output is some output container containing the results of a pCCD calculations. Replacing it with an RHF output container will give you the orbital energies for an RHF reference function.

19.5.6.2. How to: Modified Koopmans’ orbital energies

The code snippet below shows how to perform a Modified Koopmans’ orbital energy calculation in PyBEST,

orbital_energy = ModifiedKoopmans(lf, occ_model)
e_orb_output = orbital_energy(kin, ne, eri, pccd_output)

Note

The modified Koopmans’ theorem requires a pCCD reference function and is not supported for RHF.

19.5.7. Output Attributes

In addition to calculating the orbital energies from Koopmans-based method, the following attributes are provided:

e_orb_ks

single (s) orbital energies based on Koopmans’ (k) theorem

e_orb_kp_0

diagonal singlet and off-diagonal open-shell singlet pair (p) orbital energies based on Koopmans’ (k) theorem

e_orb_kp_1

diagonal triplet and off-diagonal open-shell triplet pair (p) orbital energies based on Koopmans’ (k) theorem

e_orb_mks

modification (m) ‘e_orb_ks’ to single (s) orbital energies based on Koopmans’ (k) theorem

e_orb_mkp_0

modification (m) ‘e_orb_kp_0’ to diagonal singlet and off-diagonal open-shell singlet pair (p) orbital energies based on Koopmans’ (k) theorem

e_orb_mkp_1

modification (m) ‘e_orb_kp_1’ to diagonal triplet and off-diagonal open-shell triplet (p) orbital energies based on Koopmans’ (k) theorem

e_orb_ek

orbital energies based on Extended Koopmans’ (ek) theorem (unreleased).

19.5.8. Summary of Keyword Arguments

The Koopmans and ModifiedKoopmans modules support various keyword arguments that allow us to set the ranges for printing the ionization potentials and electron affinities. In the following, all supported keyword arguments are listed together with their default values.

19.5.8.1. Orbital Energies Property Options

printoptions

(dictionary) print level:

orb_range_o

(int): Specifies the range of occupied orbital energies to display. Default is 10.

orb_range_v

(int): Specifies the range of virtual orbital energies to display. Default is 20.

all

(str): Controls the range of all orbital energies. Specific values may vary depending on the calculation requirements.

warning

(boolean) if True, (scipy) solver-specific warnings are printed (default False).

indextrans

(str) 4-index transformation. The choice between cupy, tensordot (default) and einsum. tensordot is faster than einsum.

If DenseLinalgFactory is used, the memory requirement scales roughly as \(3N^4\). Due to the storage of the two-electron integrals, the total amount of memory increases to \(4N^4\)

Note

If Cupy is not available or unsuccessful, td is selected instead.

19.5.9. Example Usage

19.5.9.1. Example Python scripts for orbital energies

Several examples are provided in data/examples/properties.

19.5.9.1.1. Orbital energies for water using RHF orbitals (cc-pVDZ)

This example calculates single and pair orbital energies based on Koopmans-type methods and RHF orbitals using the pCCD method.

Listing 19.9 data/examples/properties/water_orbital_energy_pccd.py

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.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.properties.koopmans import Koopmans
from pybest.properties.modified_koopmans import ModifiedKoopmans
from pybest.wrappers import RHF

#
# Set up molecule, define basis set
#
fn_xyz = context.get_fn("test/water.xyz")
basis = get_gobasis("cc-pvdz", fn_xyz)
#
# Define the Occupation model, orbitals, and overlap
#
lf = DenseLinalgFactory(basis.nbasis)
occ_model = AufbauOccModel(basis)
orb_a = lf.create_orbital(basis.nbasis)
olp = compute_overlap(basis)
#
# 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 pCCD calculation for the water molecule with a frozen 1s orbital
#
pccd = RpCCD(lf, occ_model)
pccd_output = pccd(kin, ne, eri, hf_output)

# Do orbital enrgy from Koopmans' theorem
# orb_range_o: (int) range of occupied orbitals
# orb_range_v: (int) range of virtual orbitals
# all : (str) all orbiatl energies are printed.
#
orbital_energyk = Koopmans(lf, occ_model)
orben_outputk = orbital_energyk(
    kin,
    ne,
    eri,
    pccd_output,
    printoptions={"orb_range_o": 3, "orb_range_v": 15},
)

# Do orbital enrgy from Modified Koopmans' theorem
# orb_range_o: (int) range of occupied orbitals
# orb_range_v: (int) range of virtual orbitals
# all : (str) all orbiatl energies are printed.
#
orbital_energymk = ModifiedKoopmans(lf, occ_model)
orben_outputmk = orbital_energymk(
    kin,
    ne,
    eri,
    pccd_output,
    printoptions={"orb_range_o": "all", "orb_range_v": "all"},
)

19.5.9.1.2. Orbital energies for water using pCCD-optimized orbitals (cc-pVDZ)

This example calculates single and pair orbital energies based on Koopmans-type methods using an OOpCCD wavefunction.

Listing 19.10 data/examples/properties/water_orbital_energy_oopccd.py

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.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.properties.koopmans import Koopmans
from pybest.properties.modified_koopmans import ModifiedKoopmans
from pybest.wrappers import RHF

#
# Set up molecule, define basis set
#
fn_xyz = context.get_fn("test/water.xyz")
basis = get_gobasis("cc-pvdz", fn_xyz)
#
# Define the Occupation model, orbitals, and overlap
#
lf = DenseLinalgFactory(basis.nbasis)
occ_model = AufbauOccModel(basis)
orb_a = lf.create_orbital(basis.nbasis)
olp = compute_overlap(basis)
#
# 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 oo-pCCD optimization for the water molecule with a frozen 1s orbital
#
oopccd = ROOpCCD(lf, occ_model)
oopccd_output = oopccd(
    kin,
    ne,
    eri,
    hf_output,
)

# Do orbital enrgy from Koopmans' theorem
# orb_range_o: (int) range of occupied orbitals
# orb_range_v: (int) range of virtual orbitals
# all : (str) all orbiatl energies are printed.
#
orbital_energyk = Koopmans(lf, occ_model)
orben_outputk = orbital_energyk(
    kin,
    ne,
    eri,
    oopccd_output,
    printoptions={"orb_range_o": 3, "orb_range_v": 15},
)

# Do orbital enrgy from Modified Koopmans' theorem
# orb_range_o: (int) range of occupied orbitals
# orb_range_v: (int) range of virtual orbitals
# all : (str) all orbiatl energies are printed.
#
orbital_energymk = ModifiedKoopmans(lf, occ_model)
orben_outputmk = orbital_energymk(
    kin,
    ne,
    eri,
    oopccd_output,
    printoptions={"orb_range_o": "all", "orb_range_v": "all"},
)