20.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).

20.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 (Extended Koopmans’ Theorem): Calculate and store ionization potentials (IPs) from Extended Koopmans’ theorem.

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.

20.5.2. Theoretical Framework

20.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.

20.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.

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

Table 20.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}\)

20.5.6. Extended Koopmans’ Theorem

Solves the generalized eigenvalue problem \(\mathbf{F C} = \boldsymbol{\varepsilon}\,\boldsymbol{\gamma C}\) using the one-particle reduced density matrix and the generalized Fock matrix of the reference wavefunction. This method incorporates electron correlation and provides only ionization potentials (no electron affinities).

20.5.7. Quick Guide: Koopmans-based orbital energies

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

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

20.5.7.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.

20.5.7.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.

20.5.7.3. How to: Extended Koopmans’ orbital energies

The code snippet below shows how to perform a Extended Koopmans’ to calculate IPs in PyBEST,

orbital_energyek = ExtendedKoopmans(lf, occ_model)
orben_outputek = orbital_energyek(kin, ne, eri, pccd_output)

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

20.5.8. 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_eks:

single (s) orbital energies based on Extended Koopmans’ (ek) theorem. Note that currently only the occupied orbital energies (or their corresponding IPs) are computed.

20.5.9. Summary of Keyword Arguments

The Koopmans, ModifiedKoopmans, and ExtendedKoopmans 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.

Note

While Koopmans’ and Modified Koopmans’ provide both IPs and EAs, Extended Koopmans’ supports only ionization potentials.

20.5.9.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.

20.5.10. Example Usage

20.5.10.1. Example Python scripts for orbital energies

Several examples are provided in data/examples/properties.

20.5.10.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 20.6 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)

# Compute orbital energies 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},
)

# Compute orbital energies 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"},
)

20.5.10.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 20.7 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,
)

# Compute orbital energies 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},
)

# Compute orbital energies 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"},
)

20.5.10.1.3. EKT example for water using RHF orbitals (cc-pVDZ)

This example calculates IPs based on Extended Koopmans’ method and RHF orbitals using the pCCD method.

Listing 20.8 data/examples/properties/water_extended_koopmans_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.extended_koopmans import ExtendedKoopmans
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)

# Compute orbital energies from extended Koopmans' theorem
# orb_range_o: (int) range of occupied orbitals
# all : (str) all orbiatl energies are printed.
# occupation_number_thresh` (float, optional): Minimum density-matrix element to consider as significant.
# Elements of the pCCD 1-RDM with absolute values below this threshold are ignored.

orbital_energyek = ExtendedKoopmans(lf, occ_model)
orben_outputek = orbital_energyek(
    kin,
    ne,
    eri,
    pccd_output,
    printoptions={"orb_range_o": "all"},
    occupation_number_thresh=5e-3,
)

20.5.10.1.4. EKT example for water using pCCD-optimized orbitals (cc-pVDZ)

This example calculates IPs based on Extended Koopmans’ method using an OOpCCD wavefunction.

Listing 20.9 data/examples/properties/water_extended_koopmans_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.extended_koopmans import ExtendedKoopmans
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,
)

# Compute orbital energies from extended Koopmans' theorem
# orb_range_o: (int) range of occupied orbitals
# all : (str) all orbiatl energies are printed.
# occupation_number_thresh` (float, optional): Minimum density-matrix element to consider as significant.
# Elements of the pCCD 1-RDM with absolute values below this threshold are ignored.

orbital_energyek = ExtendedKoopmans(lf, occ_model)
orben_outputek = orbital_energyek(
    kin,
    ne,
    eri,
    oopccd_output,
    printoptions={"orb_range_o": "all"},
    occupation_number_thresh=5e-5,
)