..
    : PyBEST: Pythonic Black-box Electronic Structure Tool
    : Copyright (C) 2016-- The PyBEST Development Team
    :
    : This file is part of PyBEST.
    :
    : PyBEST is free software; you can redistribute it and/or
    : modify it under the terms of the GNU General Public License
    : as published by the Free Software Foundation; either version 3
    : of the License, or (at your option) any later version.
    :
    : PyBEST is distributed in the hope that it will be useful,
    : but WITHOUT ANY WARRANTY; without even the implied warranty of
    : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    : GNU General Public License for more details.
    :
    : You should have received a copy of the GNU General Public License
    : along with this program; if not, see <http://www.gnu.org/licenses/>
    : --

.. _modphysham:

Model Hamiltonians
##################

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 :py:class:`~pybest.modelhamiltonians.physical_model_base.PhysModBase`. The following Hamiltonians are currently supported:

- Hückel model Hamiltonian (see :ref:`huckelham`)
- Hubbard model Hamiltonian (see :ref:`hubbardham`)
- Pariser-Parr-Pople (PPP) model Hamiltonian (see :ref:`pppham`)
- Contact interaction model in 1D (see :ref:`contactham`)

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 :ref:`user_linalg`).

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

Preliminaries
=============

Unlike full molecular Hamiltonians, physical model Hamiltonians do not require an atomic basis set to be defined. Thus, no :py:class:`~pybest.core.Basis` instance is needed.

To begin, create a :py:class:`~pybest.linalg.dense.DenseLinalgFactory` instance for the number of atoms (or sites) in the system:

.. code-block:: python

   # 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 :py:class:`~pybest.occ_model.AufbauOccModel` is supported.

.. code-block:: python

   # 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:

.. code-block:: python

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

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., :py:class:`Huckel`, :py:class:`Hubbard`, and :py:class:`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 :py:meth:`compute_one_body` with hopping parameters,
   and the two-body terms via :py:meth:`compute_two_body` with interaction parameters.

5. Model Hamiltonian classes (e.g., :py:class:`Huckel`, :py:class:`Hubbard`, :py:class:`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.

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.


.. _huckelham:

Hückel Hamiltonian
==================

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

.. math::
   :label: huckel

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

where :math:`t_{ij}` is the hopping integral between atoms :math:`i` and :math:`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

Usage Example
-------------

To construct the Hückel Hamiltonian, create an instance of the :py:class:`~pybest.modelhamiltonians.huckel_model.Huckel`
class providing a linear algebra factory and an optional xyz structure file:

.. code-block:: python

   # 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
   :py:meth:`~pybest.modelhamiltonians.huckel_model."py:func:`~pybest.modelhamiltonians.huckel_model.compute_one_body`
   to calculate the hopping term.

   The two-body term in the Hückel model is constructed using
   :py:meth:`~pybest.modelhamiltonians.huckel_model."py:func:`~pybest.modelhamiltonians.huckel_model.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.


.. _hubbardham:

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 :math:`U`.

.. math::
   :label: hubbard

   \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}

Usage Example
-------------
.. code-block:: python

   # 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. :eq:`hubbard`) is calculated
using the method :py:meth:`~pybest.modelhamiltonians.Hubbard_model.compute_one_body`
of the :py:class:`~pybest.modelhamiltonians.Hubbard` class,
while the on-site repulsion term (:math:`U`) is calculated using the
:py:meth:`~pybest.modelhamiltonians.Hubbard_model.compute_two_body` method:

.. code-block:: python

   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
   :py:meth:`~pybest.modelhamiltonians.Hubbard_model.compute_two_body`,
   which defines a minimal form of electron–electron interaction by assigning the on-site repulsion
   parameter :math:`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``.

.. _pppham:

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.

.. math::
   :label: ppp

   \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

.. math::

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


Usage Example
-------------

.. code-block:: python

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

.. _contactham:

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]_.


.. math::
   :label: contact-hamiltonian

   \hat{H}_{ii'} = \hat{T}_{ii'} + \hat{V}_{ii'}

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

.. math::
   \hat{V}_{ii'} = \delta_{ii'} V(x_i)

with :math:`x_i = i \Delta x` being a grid point,
where :math:`i` can take values of :math:`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.

Usage Example
-------------

.. code-block:: python

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


Example Python script
=====================

The Hückel Model
----------------

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


.. literalinclude:: ../src/pybest/data/examples/hamiltonian/huckel.py
    :caption: data/examples/hamiltonian/huckel.py
    :lines: 3-

The Hubbard Model
-----------------

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

.. literalinclude:: ../src/pybest/data/examples/hamiltonian/hubbard.py
    :caption: data/examples/hamiltonian/hubbard.py
    :lines: 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 :math:`t` is set to -1, while the on-site interaction :math:`U` is equal to 2.
Periodic boundary conditions are used.

.. literalinclude:: ../src/pybest/data/examples/hamiltonian/hubbard1d.py
    :caption: data/examples/hamiltonian/hubbard1d.py
    :lines: 3-

The PPP Model
-------------

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

.. literalinclude:: ../src/pybest/data/examples/hamiltonian/ppp.py
    :caption: data/examples/hamiltonian/ppp.py
    :lines: 3-

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 :py:class:`~pybest.modelhamiltonians.contact_interaction_1d.ContactInteraction1D` class)
and initial orbitals.

.. literalinclude:: ../src/pybest/data/examples/hamiltonian/contact1d.py
   :caption: data/examples/hamiltonian/contact1d.py
   :lines: 9-45