Hubbard calculations

Hubbard calculations#

# Force the local gqcpy to be imported
import sys
sys.path.insert(0, '../../build/gqcpy/')

import gqcpy
import numpy as np

Specifying the Hubbard model for a four-site chain#

If we would like to investigate a four-site chain, we should specify the adjacency matrix and the Hubbard model parameters t and U. Then, we can make a hopping matrix.

A = gqcpy.AdjacencyMatrix.Linear(4)
print(A.matrix())

[[0 1 0 0]
 [1 0 1 0]
 [0 1 0 1]
 [0 0 1 0]]

t = 1.0
U = 3.5
H = gqcpy.HoppingMatrix.Homogeneous(A, t)
print(H.matrix())

[[-0. -1. -0. -0.]
 [-1. -0. -1. -0.]
 [-0. -1. -0. -1.]
 [-0. -0. -1. -0.]]

This hopping matrix can be used to construct the Hubbard model Hamiltonian.

hubbard_hamiltonian = gqcpy.HubbardHamiltonian(H, [U]*4, [0.0]*4)
print(hubbard_hamiltonian.core().parameters())

[[ 0. -1.  0.  0.]
 [-1.  0. -1.  0.]
 [ 0. -1.  0. -1.]
 [ 0.  0. -1.  0.]]

Dense Hubbard calculations#

The first thing we should take care of is the use of the correct ONV basis.

K = 4  # number of sites
N_P = 2  # number of electron pairs

onv_basis = gqcpy.SpinResolvedONVBasis(K, N_P, N_P)

In GQCP, the Hubbard model is treated as a model Hamiltonian. What this means is that, behind the scenes, it is not a regular SQHamiltonian, but combined with the associated SpinResolvedONVBasis, it uses specialized matrix construction and matrix-vector product routines.

Even though it isn’t a SQHamiltonian, we could say it behaves as such, so doing Hubbard (CI) calculations resembles doing FCI calculations. We’ll have to specify:

a solver;
an associated environment,

and then combine these arguments in the CI QCMethod.

solver = gqcpy.EigenproblemSolver.Dense_d()
environment = gqcpy.CIEnvironment.Dense(hubbard_hamiltonian, onv_basis)

energy = gqcpy.CI(onv_basis).optimize(solver, environment).groundStateEnergy()
print(energy)

-2.1358716082319393

This is exactly the energy we find from a reference implementation by Ward Poelmans.