..
    : 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
==================

The current version of PyBEST supports the :ref:`1-dimensional Hubbard model Hamiltonian
<hubbardham>` with open and periodic boundary conditions. Note that only
``DenseLinalgFactory`` is supported for model Hamiltonians (see also
:ref:`user_linalg`).

.. _hubbardham:

The Hubbard model Hamiltonian
=============================

The Hubbard model Hamiltonian is the simplest model of interacting particles on a lattice and reads

.. math::
    :label: hubbard

    \hat{H}_{\rm Hub} = -t\sum_{j,\sigma} \left( a_{(j+1)\sigma}^{\dagger}a_{j\sigma}
    + a_{j\sigma}^{\dagger}a_{(j+1)\sigma} \right )
    +U\sum_j n_{j\uparrow} n_{j\downarrow},

where the first term is a one-electron term and accounts for the nearest-neighbor
hopping, while the second term is the repulsive on-site interaction. :math:`t`
and :math:`U` are user specified parameters and :math:`\sigma` is the electron spin.


Preliminaries
-------------

In contrast to the molecular Hamiltonian, the Hubbard model Hamiltonian does not
require a molecular geometry and atomic basis set to be defined. Thus, no
:py:class:`~pybest.gbasis.cext.GOBasis` instance is needed. Instead, the
matrix representation of the Hubbard model Hamiltonian can be immediately
evaluated.
To do so, an instance of :py:class:`~pybest.linalg.dense.DenseLinalgFactory`
has to be created,

.. code-block:: python

    # create LinalgFactory instance for n sites
    lf = DenseLinalgFactory(n)

where ``n`` is the number of sites in the Hubbard model.

.. note::

    ``CholeskyLinalgFactory`` is not supported for the Hubbard model Hamiltonian.


Defining the Hamiltonian and boundary conditions (open or periodic)
-------------------------------------------------------------------

In order to calculate the hopping and on-site repulsion terms, an instance of
the :py:class:`~pybest.modelhamiltonians.physmodham.Hubbard` class has to be
created, where an instance of the ``DenseLinalgFactory`` has to be passed,

.. code-block:: python

    # create an instance of the Hubbard class with PBC
    modelham = Hubbard(lf)

By default, the Hubbard Hamiltonian is constructed with periodic boundary conditions.
Open boundary conditions can be enforced during the initialization using

.. code-block:: python

    # create an instance of the Hubbard class with open boundary conditions
    modelham = Hubbard(lf, pbc=False)

The nearest-neighbor hopping term (:math:`t` in eq. :eq:`hubbard`) is calculated
using the method :py:meth:`~pybest.modelhamiltonians.physmodham.Hubbard.compute_one_body`
of the :py:class:`~pybest.modelhamiltonians.physmodham.Hubbard` class,

.. code-block:: python

    # strength of hopping term
    t = -1.0
    # calculate hopping term
    hopping = modelham.compute_one_body(t)

The method requires the
strength of the hopping term :math:`t` as argument.

.. note::

    The negative sign of the hopping term :math:`t` as shown in eq. :eq:`hubbard`
    is not taken into account. Thus, if the hopping term is negative,
    :math:`t` has to be chosen negative as well.

The on-site repulsion term :math:`U` is calculated using the method
:py:meth:`~pybest.modelhamiltonians.physmodham.Hubbard.compute_one_body`
of the :py:class:`~pybest.modelhamiltonians.physmodham.Hubbard` class,

.. code-block:: python

    # strength of the on-site interactions
    U = 2.0
    # calculate on-site interaction
    onsite = modelham.compute_one_body(U)

Similar to the hopping term, the method requires the strength of the on-site interaction
:math:`U` as argument. If :math:`U` is positive, the on-site interaction is
repulsive.

.. note::

    The hopping term is a two-index object of ``DenseLinalgFactory``, while the
    on-site interaction is a four-index object of ``DenseLinalgFactory``


Filling the lattice
-------------------

The number of electrons/spin-half particles on the lattice is assigned in the
Aufbau occupation number model and specifying the total number of electrons
using the ``nel`` keyword argument. An instance of the
:py:class:`~pybest.scf.occ.AufbauOccModel` class is used for setting
the occupations in all electronic structure modules.

.. code-block:: python

    # define the number of doubly filled lattice sites for m electrons
    occ_model = AufbauOccModel(lf, nel=m)

See also :ref:`user_occupation_model` for more details on the occupation model
used in PyBEST.

.. note::

    The number of electrons/spin-half particles must be even and the orbital
    basis (or the on-site basis) must be restricted.


Generating the orbital basis and the overlap matrix
---------------------------------------------------

In all electronic structure calculations with the Hubbard model Hamiltonian,
the model Hamiltonian is evaluated for the local on-site basis and the orbitals
represent the transformation from the on-site basis to the molecular orbital
basis. Note that the on-site basis is orthonormal.
As for the molecular Hamiltonian, the orbitals (or, equivalently, the on-site basis)
are initialized using the :py:meth:`~pybest.linalg.dense.DenseLinalgFactory.create_orbital`
method of the :py:class:`~pybest.linalg.dense.DenseOrbital` class.

.. code-block:: python

    # initialize orbitals for n sites (assuming lf = DenseLinalgFactory(n))
    orb = lf.create_orbital()

See alse :ref:`user_wroking_with_orbitals` for more details on
orbitals and how to work with orbitals in PyBEST.
The overlap matrix of the Hubbard model can be computed using the
:py:meth:`~pybest.modelhamiltonians.physmodham.Hubbard.compute_overlap`
method of the :py:class:`~pybest.modelhamiltonians.physmodham.Hubbard` class,

.. code-block:: python

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


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

The 1-dim Hubbard model Hamiltonian with PBC
--------------------------------------------

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/hubbard.py
    :caption: data/examples/hamiltonian/hubbard.py
    :lines: 3-5,7-29