# 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
]),
}