2.3. Model Hamiltonians

2.3.1. Supported Features

PyBEST includes a flexible framework for model Hamiltonians relevant to quantum chemistry and condensed matter physics. These models are implemented through a hierarchy of classes based on the abstract base class PhysModBase. The following Hamiltonians are currently supported:

The Hückel and PPP models are built using adjacency matrices derived from covalent radii, allowing general applicability to π-conjugated systems. The Hubbard model is implemented as an extension of the Hückel model, and in the absence of an XYZ file, it defaults to a 1D lattice-based site model with support for both open and periodic boundary conditions. The PPP model extends the Hubbard model with long-range Coulomb interactions (V-terms). The Hamiltonian for 1D contact interactions can be combined with any arbitrary external potential.

Note

Only DenseLinalgFactory is supported for model Hamiltonians (see also PyBEST objects: the LinalgFactory).

If you use model Hamiltonians, please cite [karimi2025].

2.3.2. Preliminaries

Unlike full molecular Hamiltonians, physical model Hamiltonians do not require an atomic basis set to be defined. Thus, no Basis instance is needed.

To begin, create a DenseLinalgFactory instance for the number of atoms (or sites) in the system:

# The number of sites is represented as a `LinalgFactory` object (indicating the number of supported atoms).
lf = DenseLinalgFactory(n)

The xyz structure can optionally be provided using the .xyz format by passing the xyz_file argument. In the 1D Hubbard model, we default to a site model if an xyz file is not provided (xyz_file defaults to None).

To work with physical model Hamiltonians, you must assign a proper occupation model representing your physical problem, that is, how many sites or states are filled by electrons. Currently, only the AufbauOccModel is supported.

# We assume n sites with one electron each, resulting in n/2 doubly-filled sites
occ_model = AufbauOccModel(lf, nel=n)

This allows later reuse in any model class constructor:

# Calculate overlap matrix for the on-site basis
olp = modelham.compute_overlap()

2.3.3. Computing the Hamiltonian

The general procedure to calculate model Hamiltonians in PyBEST is uniform across models and involves the following steps:

  1. Instantiate the model Hamiltonian class (e.g., Huckel, Hubbard, and PPP), providing the linear algebra factory, occupation model, and optionally the xyz structure.

  2. Specify model-specific parameters as a dictionary, such as hopping strengths, on-site energies, interaction parameters (U, V), boundary conditions, etc.

  3. Call the instance with the parameter dictionary to compute the Hamiltonian, which returns an object containing the one-body and two-body terms.

  4. The one-body terms are calculated using the method compute_one_body() with hopping parameters, and the two-body terms via compute_two_body() with interaction parameters.

  5. Model Hamiltonian classes (e.g., Huckel, Hubbard, PPP) support an optional "rhf" boolean parameter. When set to True (default), the model will internally perform a Restricted Hartree-Fock (RHF) calculation during its evaluation and return the result.

    To skip the RHF calculation (e.g., if you only want the tight-binding data), set "rhf": False.

  6. The resulting Hamiltonian can then be used seamlessly in wave function methods such as RHF, pCCD, or other post-HF calculations.

Examples for each model are provided in their respective sections below.

2.3.3.1. Specifying Model Hamiltonian Parameters

Model parameters are passed as a dictionary when invoking the model instance, allowing flexible tuning of physical quantities:

  • For the Hückel model, typical parameters include:

    • hopping: nearest-neighbor hopping integral (e.g., -1.0)

    • on_site: on-site energy (often 0.0)

    • rhf: boolean to enable RHF calculation.

  • For the Hubbard model, the Hückel parameters dictionary contains the additional keys:

    • u: on-site interaction strength

    • pbc: boolean for periodic boundary conditions (for the 1D-Hubbard model).

  • For the PPP model, additional keys are:

    • k: dielectric constant for screening in the Coulomb interaction

    • hubbard: boolean to enable Hubbard-type interactions

    • u and hopping similarly defined, often in electronvolts

  • For the 1D Contact Interaction, accepted parameter keys are:

    • grid: tuple specifying the spatial domain and discretization, given as (x_min, x_max, dx) (e.g., (-10.0, 10.0, 0.05)), which defines the DVR grid.

    • mass: mass of the particle in atomic units (default is 1.0), representing the fermion mass used in kinetic energy calculations.

    • g_coupling: coupling strength of the contact interaction (e.g., 2.0), scaling the two-body interaction integrals.

Changing these parameters allows modeling different physical regimes or fitting to experimental data.

Note

In Hückel, Hubbard and PPP model Hamiltonians, if you want to run only the RHF calculation with the given Hamiltonian, set "rhf": True in the parameters dictionary and call the model Hamiltonian instance once.

2.3.4. Hückel Hamiltonian

The Hückel model Hamiltonian features the tight-binding and Hartree-Fock methods for electronic structure calculations.

(2.2)\[\hat{H}_{\rm Huckel} = -\sum_{i,j} t_{ij} \left( c_i^\dagger c_j + c_j^\dagger c_i \right)\]

where \(t_{ij}\) is the hopping integral between atoms \(i\) and \(j\). The molecular structure is encoded through the adjacency matrix derived from covalent radii and interatomic distances. The list of currently supported atoms is

  • Carbon

  • Oxygen

  • Nitrogen

  • Fluorine

  • Aluminum

  • Silicon

  • Phosphorus

  • Sulfur

  • Chlorine

  • Gallium

  • Germanium

  • Arsenic

  • Selenium

  • Bromine

2.3.4.1. Usage Example

To construct the Hückel Hamiltonian, create an instance of the Huckel class providing a linear algebra factory and an optional xyz structure file:

# Instantiate the Hückel class providing the xyz structure filename (an xyz file, here ``"mol.xyz"``)
modelham = Huckel(lf, occ_model, xyz_file="mol.xyz")

# Use modelham to compute Hückel Hamiltonian with the following parameters
# "rhf"=True runs an RHF calculation with the defined Hamiltonian
huckel_output = modelham(
    parameters={"on_site": 0.0, "hopping": -1.0, "rhf": True}
)

# Set the parameters internally (without running RHF) so that compute_one_body
# and compute_two_body methods can access these parameters when called
modelham(parameters={"on_site": 0.0, "hopping": -1.0})

# compute the one-body term (tight-binding Hamiltonian)
# Calculate hopping term
one_body = modelham.compute_one_body()
# Compute the two-body term (zero tensor for Hückel, required for RHF/post-HF)
two_body = modelham.compute_two_body()

Note

The hopping parameter hopping is directly used in the method compute_one_body() to calculate the hopping term.

The two-body term in the Hückel model is constructed using compute_two_body(), which returns a zero-initialized four-index tensor of electron repulsion integrals labeled "eri". While the Hückel model neglects electron–electron interactions, this zero-valued tensor is still required to satisfy the interface expected by PyBEST’s Hartree–Fock (RHF) and post-HF methods. It has to be created in order to run RHF (or any post-RHF) calculations, as otherwise PyBEST will raise an error.

By defining both the one-body and two-body parts, the Hückel Hamiltonian is fully defined and can be seamlessly used in any wave function-based method available in PyBEST.

2.3.5. Hubbard Hamiltonian

The Hubbard model extends the Hückel framework by adding an on-site electron-electron interaction (U-term), capturing basic electron correlation effects. It can be understood as a tight-binding model with an additional on-site repulsion (or attraction) term parameterized by \(U\).

(2.3)\[\hat{H}_{\rm Hubbard} = -\sum_{i,j} t_{ij} \left( c_i^\dagger c_j + c_j^\dagger c_i \right) +U\sum_j n_{j\uparrow} n_{j\downarrow}\]

2.3.5.1. Usage Example

# Instantiate the Hubbard class providing the xyz structure filename  (an xyz file, here "mol.xyz")
modelham = Hubbard(lf, occ_model, xyz_file="mol.xyz")

# If the xyz file is not provided, a site model is assumed
# Instantiate the 1D-Hubbard class with periodic boundary condition (parameters are set below)
modelham = Hubbard(lf, occ_model, pbc=True)

# Use modelham to compute the Hubbard Hamiltonian with the following parameters
hubbard_output = modelham(
    parameters={"on_site": 0.0, "hopping": -1.0, "u": 1.0}
)

The nearest-neighbor hopping term (hopping in eq. (2.3)) is calculated using the method compute_one_body() of the Hubbard class, while the on-site repulsion term (\(U\)) is calculated using the compute_two_body() method:

modelham(parameters={"on_site": 0.0, "hopping": -1.0, "u": 1.0})

# Calculate hopping term
one_body = modelham.compute_one_body()
# Calculate on-site electron-electron repulsion term
two_body = modelham.compute_two_body()

Note

The two-body term in the Hubbard model is constructed using compute_two_body(), which defines a minimal form of electron–electron interaction by assigning the on-site repulsion parameter \(U\) to the diagonal elements of the four-index electron repulsion tensor. As in the Huckel models within PyBEST, both the one-body and two-body components must be defined to construct a valid Hamiltonian. The inclusion of the two-body term is necessary to run HF or any post-HF methods, as PyBEST requires an electron repulsion integrals labeled "eri" for all wave function-based calculations.

  • The hopping term is a two-index object of DenseLinalgFactory.

  • The on-site interaction is a four-index object of DenseLinalgFactory.

2.3.6. PPP Hamiltonian

The Pariser-Parr-Pople (PPP) model extends the Hubbard model by including long-range Coulomb interactions in addition to the on-site interaction. It is particularly well-suited for accurate modeling of large π-conjugated systems.

(2.4)\[\hat{H}_{\mathrm{PPP}} = - \sum_{i,j} t_{ij} \left( c_i^\dagger c_j + c_j^\dagger c_i \right) + U \sum_i n_{i\uparrow} n_{i\downarrow} + \sum_{i<j} V_{ij} (n_i - 1)(n_j - 1)\]

with

\[V_{i,j} = \frac{U}{\kappa_{i,j} \sqrt{1 + 0.6117 R_{i,j}^2}}\]

2.3.6.1. Usage Example

# Initialize PPP class with the xyz structure file.
modelham = PPP(lf, occ_model, xyz_file="mol.xyz")

# Use modelham to compute PPP Hamiltonian with the following parameters
# Multiply parameters like hopping and U by `electronvolt` if eV units are needed.
ppp_output = modelham(
    parameters={
        "on_site": 0.0,
        "hopping": -2.7 * electronvolt,
        "u": 8.0 * electronvolt,
        "k": 2.0,
        "hubbard": False,
    }
)

# Overwrite parameters WITHOUT doing an RHF calculation
modelham(parameters={"on_site": 0.0, "hopping": -1.0, "u": 1.0, "k":2.0, "hubbard": False})

# Calculate the hopping (one-body) term, which is the same as the Hubbard model
# Integrals are evaluated for the second set of parameters
one_body = modelham.compute_one_body()

# Calculate electron-electron repulsion (two-body) term,
# Including long-range interactions using the Ohno formula (see code for reference), with an optional on-site Hubbard term.
# Integrals are evaluated for the second set of parameters
two_body = modelham.compute_two_body()

Note

  • The parameters["hopping"] and parameters["u"] values are directly used in the corresponding methods to calculate the hopping and on-site interaction terms.

2.3.7. Contact Interaction Hamiltonian

This is a simplified 1D model using delta-function (contact) interactions, often used as a pedagogical or benchmark model in low-dimensional quantum systems [dvr-1991].

(2.5)\[\hat{H}_{ii'} = \hat{T}_{ii'} + \hat{V}_{ii'}\]

Here, \(\hat{T}_{ii'}\) represents the 1D kinetic energy operator, and \(\hat{V}_{ii'}\) is an arbitrary external potential

\[\hat{V}_{ii'} = \delta_{ii'} V(x_i)\]

with \(x_i = i \Delta x\) being a grid point, where \(i\) can take values of \(0, \pm 1, \pm 2, \dots\). This choice of the Hamiltonian leads to a uniform grid in the coordinate space and corresponds to using a Fourier basis set.

For the 1D Contact Interaction Hamiltonian, a user-specified grid, particle mass, and 1-variable external potential function can be defined. If not provided, the following default values for the function arguments are assumed:

domain

list containing [x_min, x_max, dx] which defines the grid for DVR (discrete variable representation). The default grid is [-10, 10] with a step size of 0.1, which is suitable for the simple harmonic oscillator.

mass

mass of the particle in atomic units (a.u.) to avoid confusion.

potential

external 1-variable function defining the potential. The default external potential is a harmonic oscillator.

2.3.7.1. Usage Example

# Define model parameters
 grid = (-10.0, 10.0, 0.05)
 g_coupling = 2.0
 mass = 1.0

 # Instantiate the ContactInteraction1D class with default values for
 # the grid (domain argument), external potential, and electron mass
 modelham = ContactInteraction1D(lf, occ_model, domain=grid, mass=mass)

 # Compute the one- and two-body interaction terms (integral tensors)
 one = modelham.compute_one_body()
 two = modelham.compute_two_body()
 # scale 2-body integrals by the coupling strength
 two.iscale(g_coupling)

Note

This is a minimal model used to test 1D correlation effects.

2.3.8. Example Python script

2.3.8.1. The Hückel Model

This example builds a Hückel model using some xyz structure. The hopping term \(t\) is set to -1.0 and the on-site energy \(\varepsilon\) to 0.0. RHF is enabled, and hence an RHF optimization is performed.

Listing 2.9 data/examples/hamiltonian/huckel.py
from pybest.linalg import DenseLinalgFactory
from pybest.modelhamiltonians.huckel_model import Huckel
from pybest.occ_model import AufbauOccModel

# get the xyz file from pybest/src/pybest/data/test
fn_xyz = context.get_fn("test/c22h12.xyz")

# Number of sites represented as a `LinalgFactory` object (indicating the number of supported atoms).
lf = DenseLinalgFactory(22)

# Define the occupation model where `nel` is the number of C-H bonding and lone-pair electrons.
occ_model = AufbauOccModel(lf, nel=22)

# t=-1 and epsilon=0 are default hopping and on-site parameters, respectively.
modelham = Huckel(lf, occ_model, xyz_file=fn_xyz)
huckel_output = modelham(
    parameters={"on_site": 0.0, "hopping": -1.0, "rhf": True}
)

2.3.8.2. The Hubbard Model

This example builds a Hubbard model using some xyz structure. The hopping \(t\) is set to -1.0, on-site energy \(\varepsilon\) to 0.0, and Hubbard interaction \(U\) to 1.0. RHF is enabled.

Listing 2.10 data/examples/hamiltonian/hubbard.py
from pybest.linalg import DenseLinalgFactory
from pybest.modelhamiltonians import Hubbard
from pybest.occ_model import AufbauOccModel

# get the xyz file from pybest/src/pybest/data/test
coord = context.get_fn("test/c22h12.xyz")

# Number of sites represented as a `LinalgFactory` object (indicating the number of supported atoms).
lf = DenseLinalgFactory(22)

# Define the occupation model where `nel` is the number of C-H bonding and lone-pair electrons.
occ_model = AufbauOccModel(lf, nel=22)

# t=-1 and epsilon=0 are default hopping and on-site parameters, respectively.
modelham = Hubbard(lf, occ_model, xyz_file=coord)
hubbard_output = modelham(
    parameters={"on_site": 0.0, "hopping": -1.0, "u": 1.0},
)

2.3.8.3. The 1D-Hubbard Model

This example shows how to set up the Hamiltonian, orbitals, and overlap matrix for the half-filled Hubbard model with 6 sites (and hence 3 doubly occupied sites). The hopping term \(t\) is set to -1, while the on-site interaction \(U\) is equal to 2. Periodic boundary conditions are used.

Listing 2.11 data/examples/hamiltonian/hubbard1d.py
from pybest.modelhamiltonians import Hubbard
from pybest.occ_model import AufbauOccModel

# Number of sites represented as a `LinalgFactory` object.
lf = DenseLinalgFactory(6)

# If xyz_file is not provided, this would be a site model.
# Define the occupation model where `nel` is the number of electrons in the 1D Hubbard model.
occ_model = AufbauOccModel(lf, nel=6)

# t=-1 and epsilon=0 are default hopping and on-site parameters, respectively.
modelham = Hubbard(lf, occ_model, pbc=True)
hubbard_output = modelham(
    parameters={"on_site": 0.0, "hopping": -1.0, "u": 1.0},
)

2.3.8.4. The PPP Model

This example shows the setup for a PPP model using molecular structure coordinates. It uses a hopping term \(t = -2.7\,\mathrm{eV}\), on-site repulsion \(U = 8.0\,\mathrm{eV}\), and dielectric constant \(k = 1.0\). Hubbard-type interaction and RHF are enabled.

Listing 2.12 data/examples/hamiltonian/ppp.py
from pybest import context
from pybest.linalg import DenseLinalgFactory
from pybest.modelhamiltonians.ppp_model import PPP
from pybest.occ_model import AufbauOccModel
from pybest.units import electronvolt

# get the xyz file from pybest/src/pybest/data/test
coord = context.get_fn("test/benzeo.xyz")

# Number of sites represented as a `LinalgFactory` object (indicating the number of supported atoms).
lf = DenseLinalgFactory(22)

# Define the occupation model where `nel` is the number of C-H bonding and lone-pair electrons.
occ_model = AufbauOccModel(lf, nel=22)
orb_a = lf.create_orbital()

# t: hopping, u: e-e repulsion, k: dielectric constant, u_p=u/k, hubbard: hubbard term in ppp.
modelham = PPP(lf, occ_model, xyz_file=coord)

ppp_output = modelham(
    parameters={
        "on_site": 0.0,
        "hopping": -2.7 * electronvolt,
        "u": 8.0 * electronvolt,
        "k": 1.0,
        "u_p": 0.8 * electronvolt,
        "hubbard": True,
        "rhf": True,
    }
)

2.3.8.5. The 1D Contact Interaction Model

This example shows the setup for a one-dimensional harmonic trap model. It uses a user-defined uniform grid and a coupling strength to scale the two-body interactions. An RHF calculation is performed outside the model Hamiltonian class, which requires the calculation of the overlap integrals (as implemented in the ContactInteraction1D class) and initial orbitals.

Listing 2.13 data/examples/hamiltonian/contact1d.py
from pybest.modelhamiltonians import ContactInteraction1D
from pybest.occ_model import AufbauOccModel
from pybest.wrappers import RHF

# Define LinalgFactory for no_orbs = 30
# The maximal number of orbitals is determined by the grid size
no_orbs = 30

# numer of fermions  = 2
no_fermions = 2

# Define grid, mass and coupling strength as in
# paper of Grining, 2015,  Eq (1).
# several literature values for 2 paired (up and down) electrons:
# g /   total energy (Eh)
# 2     1.536605
# 0     0.999991
# -4   -1.816517
# note that for g=0 system corresponds to non-interacting oscillators
grid = (-10.0, 10.0, 0.05)
mass = 1.0
g_coupling = 2.0

# Define linear algebra and occupation model
lf = DenseLinalgFactory(no_orbs)
occ_model = AufbauOccModel(lf, nel=no_fermions)

# Calculate a set of 2-body integrals (no_orbs^4)
# by numerical integration, make 1-body integrals
modelham = ContactInteraction1D(lf, occ_model, domain=grid)
olp = modelham.compute_overlap()
one = modelham.compute_one_body()
two = modelham.compute_two_body()
orbs = lf.create_orbital()

# Scale 2-body integrals by the coupling strength
two.iscale(g_coupling)