# 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

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.

```
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)
```