3.2. Defining a model Hamiltonian¶

3.2.1. Supported features¶

HORTON implements the 1- dimensional Hubbard model Hamiltonian with open and periodic boundary conditions. Note that only DenseLinalgFactory is supported for model Hamiltonians.

3.2.2. The Hubbard model Hamiltonian¶

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

(1)\[\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. The \(t\) and \(U\) are user specified parameters and \(\sigma\) is the electron spin.

3.2.2.1. Preliminaries¶

In contrast to the molecular Hamiltonian, the Hubbard Hamiltonian does not require a molecule and GOBasis instance. Instead, a DenseLinalgFactory instance can be immediately created,

lf = DenseLinalgFactory(n)

Note that CholeskyLinalgFactory is not supported.

3.2.2.2. Defining the Hamiltonian and boundary conditions (open or periodic)¶

Before the hopping and on-site repulsion terms can be calculated, you need to create an instance of Hubbard:

modelham = Hubbard()

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

modelham = Hubbard(pbc=False)

The nearest-neighbor hopping term (\(t\) in eq. (1)) is calculated using the method compute_kinetic():

hopping = modelham.compute_kinetic(lf, t)

The the on-site repulsion \(U\) is calculated using the method compute_er():

onsite = modelham.compute_er(lf, U)

Finally, all terms of the 1-dimensional Hubbard Hamiltonian are combined together and passed to the effective Hamiltonian class, REffHam, which can then be passed to the restricted Hartree-Fock or DFT modules,

terms = [
        RTwoIndexTerm(hopping, 'kin'),
        RDirectTerm(onsite, 'hartree'),
        RExchangeTerm(onsite, 'x_hf'),
        ]
ham = REffHam(terms)

3.2.2.3. Filling the lattice¶

The number of electrons/spin-half particles on the lattice is assigned using the Aufbau occupation number model. An instance of the class AufbauOccModel will be used later for setting the occupations.

occ_model = AufbauOccModel(m)

Note that the number of electrons/spin-half particles must be even and orbitals must be restricted.

3.2.2.4. Generate initial guess orbitals and overlap matrix¶

The instance of class DenseExpansion by using the method create_expansion() is used to initiate the orbitals:

orb = lf.create_expansion(n)

An initial guess can be obtained with guess_core_hamiltonian():

guess_core_hamiltonian(olp, hopping, orb)

The overlap matrix of the Hubbard model can be computed using the method compute_overlap():

olp = modelham.compute_overlap(lf)

3.2.3. Example Python script¶

3.2.3.1. Restricted Hartree-Fock calculations using the 1-dim Hubbard model Hamiltonian with PBC¶

This example shows a restricted Hartree-Fock calculation for the half-filled Hubbard model. Both the number of electron spins and sites is 6. The \(t\) parameter is set to -1, while the \(U\) parameter is equal to 2. Periodic boundary conditions are used.

data/examples/hamiltonian/hubbard.py¶

from horton import *  # pylint: disable=wildcard-import,unused-wildcard-import


# Define Occupation model, expansion coefficients and overlap
# -----------------------------------------------------------
lf = DenseLinalgFactory(6)
occ_model = AufbauOccModel(3)
modelham = Hubbard(pbc=True)
orb = lf.create_expansion(6)
olp = modelham.compute_overlap(lf)


# One and two-body interaction terms
# ----------------------------------

# t-param, t = -1
hopping = modelham.compute_kinetic(lf, -1)
# U-param, U = 2
onsite = modelham.compute_er(lf, 2)


# Perform initial guess
# ---------------------
guess_core_hamiltonian(olp, hopping, orb)
terms = [
    RTwoIndexTerm(hopping, 'kin'),
    RDirectTerm(onsite, 'hartree'),
    RExchangeTerm(onsite, 'x_hf'),
]
ham = REffHam(terms)


# Do a Hartree-Fock calculation
# -----------------------------
scf_solver = PlainSCFSolver()
scf_solver(ham, lf, olp, occ_model, orb)


# CODE BELOW IS FOR horton-regression-test.py ONLY. IT IS NOT PART OF THE EXAMPLE.
rt_results = {
    'energy': ham.cache["energy"],
    'orb_energies': orb.energies,
    'kinetic': ham.cache["energy_kin"],
    'hartree': ham.cache["energy_hartree"],
    'exchange': ham.cache["energy_x_hf"],
}
# BEGIN AUTOGENERATED CODE. DO NOT CHANGE MANUALLY.
import numpy as np  # pylint: disable=wrong-import-position
rt_previous = {
    'energy': -5.0,
    'exchange': -2.9999999999999982,
    'hartree': 5.9999999999999964,
    'kinetic': -7.999999999999998,
    'orb_energies': np.array([
        -1.0000000000000007, -2.3948019172601905e-16, 1.2845788926350238e-16,
        1.9999999999999991, 2.0000000000000004, 2.9999999999999991
    ]),
}