Download this notebook (.ipynb + data files)
3. The square Hubbard Model at half filling#
Author: Ryan Levy
import os
# This are in case you have a GPU enabled JAX but are running on CPU.
os.environ['JAX_PLATFORMS'] = 'cpu'
# If you want to use GPU enabled Jax in multiply notebooks with the same
os.environ["XLA_PYTHON_CLIENT_PREALLOCATE"] = "false"
# use x64 numbers
from jax import config
config.update("jax_enable_x64", True)
import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt
import numpy as np
# all the imports used later in the tutorial, but put here for convenience
import afqmctools.systems.lattice as lat
import afqmctools.utils.visualize as vis
import afqmctools.utils.io as io
import afqmctools.hamiltonian.model.director as ham
from afqmctools.wavefunction.converter import read_wavefunction
from afqmctools.wavefunction.model import write_free_electron_wfn, make_free_elec,write_wfn
from afqmctools.inputs.from_hdf import write_json
from afqmctools.analysis.rdm import average_afqmc_rdm
from afqmctools.observables.greens import greens_1body
import afqmctools.observables.spin as spobs
from stats.scalar_dat import analyze_scalar_data
from tutorial_utils import run_afqmc
import autohf
from afqmctools.inputs.from_autohf import autohf_to_afqmc
from pathlib import Path
scratch_dir = Path("data")
scratch_dir.mkdir(parents=True, exist_ok=True)
3.1. Physics Background#
The Hubbard model at Half Filling (\(n=1\)) has a special “symmetry” known as particle-hole symmetry, that allows auxiliary field methods to be sign free. However, care must be taken to make sure that phaseless/constrained path AFQMC is used properly to exploit said symmetry.
[1]: Benchmark study of the two-dimensional Hubbard model with auxiliary-field quantum Monte Carlo method
M Qin, H Shi, S Zhang Phys. Rev. B 94, 085103 4 August, 2016
DOI: https://doi.org/10.1103/PhysRevB.94.085103
3.2. Lattice and Hamiltonian Setup#
# take from 10.1103/PhysRevB.45.10741, U => E_0
E0_n1_lookup = {4 : -13.62185,
8 : -8.46888}
U = 4.0
lattice_dims = (4,4)
lattice_params = {
'L1' : lattice_dims[0],
'L2' : lattice_dims[1],
'boundary1' : 'PBC',
'boundary2' : 'PBC'
}
lattice = lat.get_lattice(params=lattice_params)
Ne = int(lattice.N_sites)//2
nelec = (Ne,Ne)
nelec
vis.plot_lattice(lattice)
plt.show()
The Hamiltonian will be
builder = ham.HamiltonianBuilder(
lattice=lattice,
spin_symm="collinear" # we have no spin-flip terms
)
# add standard Hubbard terms
builder.nth_neighbor_hopping(1.0)
builder.onsite_hubbard(U)
builder.finalize()
builderHF = ham.HamiltonianBuilder(
lattice=lattice,
spin_symm="collinear" # we have no spin-flip terms
)
# add standard Hubbard terms
builderHF.nth_neighbor_hopping(1.0)
builderHF.onsite_hubbard(4) # if you want an effective U, you can change this
builderHF.finalize()
# get the 1 body Hamiltonian
T = builderHF.hamiltonian.get_one_body().toarray().reshape(2,lattice.N_sites,lattice.N_sites)
3.3. Version 1: GHF trial#
Because of the open shell at half filling, we’ll follow Ref. 1 and use a GHF trial w/ AFM constrained to the X-Y plane (eq 14)
where \(M_i= (-1)^i\Delta_i\)
3.3.1. Aside: custom Hartree-Fock solver#
We’ll want to optimize \(M\) directly, and to do this we’ll use AutoHF’s custom ansatz power to build a modification of the DIAG ansatz
To do that we’ll need:
orbitalFuncThis takes astateand produces the orbitalsrdmFuncthis takes astateand produces the rdm \(\langle c^\dagger_{\sigma i}c_{\sigma^\prime j}\rangle\)state0the initial statestate0_refthe initial state that should correspond to e.g. \(U=0\)
N = lattice.N_sites # will be captured by closure
Ne = N//2 # n_up = n_dn = Ne
# get the (-1)^i which is really the checkerboard pattern for a 4x4
signs = np.array([ (-1)**np.sum(s.coord) for s in lattice.sites])
# convert T -> GHF format
# TODO this should be in afqmctools
TGHF_ = np.block([[T[0],np.zeros_like(T[1])],
[np.zeros_like(T[0]),T[1]]])
TGHF_ = jnp.array(TGHF_)
def orbitalFunc(state):
# take the hopping matrix TGHF, and "set" M on the off-diagonal blocks
TGHF = TGHF_+0.0
TGHF = TGHF.at[N+np.arange(N),np.arange(N)].set(signs*state)
TGHF = TGHF.at[np.arange(N),N+np.arange(N)].set(signs*state)
lams,vecs = jnp.linalg.eigh(TGHF)
return jnp.stack([vecs[:,:2*Ne]])
def rdmFunc(state):
o = orbitalFunc(state)
return autohf.solver.sd_to_rdm_noncollinear(o)
def checkPH(state):
'''check the spectrum of the 'faux' HF state is PH symmetric'''
# take the hopping matrix TGHF, and "set" M on the off-diagonal blocks
TGHF = TGHF_+0.0
TGHF = TGHF.at[N+np.arange(N),np.arange(N)].set(signs*state)
TGHF = TGHF.at[np.arange(N),N+np.arange(N)].set(signs*state)
lams,vecs = jnp.linalg.eigh(TGHF)
# verify symmetry
return np.allclose(lams,-lams[::-1])
hf_settings = {
"verbose": False,
"plot": False,
"steps": 200,
"seed": 1479,
"noncollinear": True,
"gpu": False,
"ansatz": "CUSTOM",
"opt_method": "lbfgs",
"batch_size": 16,
"nelec": [Ne, Ne],
}
rng = np.random.default_rng(42)
state0_ref = 0.0 * np.arange(N).reshape(N)
state0 = state0_ref + rng.normal(scale=0.01, size=(hf_settings["batch_size"], N))
dataHFC = autohf.solver.lattice_hf(
hamiltonian=autohf.AutoHFHamiltonian(builderHF.hamiltonian),
lattice=lattice,
settings=hf_settings,
state2orbitals=orbitalFunc,
state2rdm=rdmFunc,
state0=state0,
state0_ref=state0_ref,
)
Let’s check out what M values we got
dataHFC[0]["state"]
plt.imshow(dataHFC[0]["state"].reshape(4,4))
Hmm… This looks like we haven’t fully converged into a translationally invariant solution, but’s close. Let’s increase the precision!
⚠️ Warning This notebook is setup for the CPU, the following cell takes a while (the GPU is a bit faster if you have one)
wfn_fname = "autoHF_wfn.h5"
hf_settings = {
"verbose": False,
"plot": False,
"steps": 2000, # increase
"seed": 1479,
"noncollinear": True,
"gpu": False,
"ansatz": "CUSTOM",
"opt_method": "lbfgs",
"num_batches": 16,
"nelec": [Ne, Ne],
}
rng = np.random.default_rng(42)
state0_ref = 0.0 * np.arange(N).reshape(N)
state0 = state0_ref + rng.normal(scale=0.01, size=(hf_settings["num_batches"], N))
dataHFC2 = autohf.solver.lattice_hf(
hamiltonian=autohf.AutoHFHamiltonian(builderHF.hamiltonian),
lattice=lattice,
settings=hf_settings,
state2orbitals=orbitalFunc,
state2rdm=rdmFunc,
state0=state0,
state0_ref=state0_ref,
jaxoptargs=dict(tol=1e-14), # secret sauce
)
autohf_to_afqmc(
dataHFC2,
output_fname = scratch_dir / wfn_fname
)
dataHFC2[0]["state"]
checkPH(dataHFC2[0]["state"])
plt.imshow(dataHFC2[0]["state"].reshape(4,4))
Let’s check that this is a better energy than the previous energy
dataHFC2[0]["E"] < dataHFC[0]["E"]
Great!
Now let’s check if we can guess a better state, by removing translation invariance
dataHFC2[1]["energy_func"](dataHFC2[0]["state"]) < dataHFC2[1]["energy_func"](dataHFC2[0]["state"]/np.abs(dataHFC2[0]["state"]))
uni_state = np.ones(lattice.N_sites)*np.mean(np.abs(dataHFC2[0]["state"]))
np.isclose(dataHFC2[1]["energy_func"](dataHFC2[0]["state"]),dataHFC2[1]["energy_func"](uni_state))
So we see we’ve reached a minimum, and that a constant \(M\) is the same as our nearly constant optimized \(M_i\)
3.3.2. Constructing the input file#
First we’ll manually create a non-interacting state as our RHF initial state for the walkers
ham_fname = f"hamU{U}_afqmc.h5"
io.write_model_hamiltonian(builder.hamiltonian, scratch_dir / ham_fname,
nelec=nelec,spin_symm="collinear")
afqmc_params = {
"timestep": 0.005,
"steps": 10000,
"population_control_interval": 2,
"walker_ortho_interval": 2,
"n_walkers_per_mpi_task": 40,
"measure_interval_multiplier": 1,
"estimator": {
"name": "energy",
"overwrite": True,
"print_components": True
},
}
write_json(
scratch_dir / "afqmc.json",
fwfn0=scratch_dir / wfn_fname,
fham0=scratch_dir / ham_fname,
exec_opts=afqmc_params,
series=0
)
run_afqmc(
run_dir = scratch_dir,
input_file = "afqmc.json",
np = 16,
output_file = None,
)
!energy_stats {scratch_dir / "qmc.s000.scalar.dat"} -t -e 20.0
settings = dict(
fname = scratch_dir / "qmc.s000.scalar.dat",
nequil = 14,
trace = True,
return_data= True,
verbose = False,
)
data = analyze_scalar_data(settings)
E = data["mean_real"]
err = data["error_real"]
E0_n1_lookup[U]
plt.errorbar([0],[E],yerr=[err],capsize=6,fmt='o-')
plt.axhline(y=E0_n1_lookup[U],ls='--')
plt.ylim(-13.63,-13.6)
plt.show()
3.3.3. For further thinking#
What if we didn’t use the correct signs (\(M_i = (-1)^i\Delta_i\)) for the HF ansatz? Could we construct an even more constrained wave function?
What do the local spin/charge observables of the optimized ansatz look like?
What happens for \(U=8\)? Do we need a different wave function?
What happens if you start from a free electron state instead of the HF state? (
(coeffs,free_wfn),spin_symm = make_free_elec(ham_fname,nelec,spin_symm="nc"))
3.4. Version 2: multi-Slater trial#
To take advantage of particle hole symmetry at half filling with a spin decomposition, besides using a GHF in the X-Y plane, we could use two particle-hole symmetric determinants
# TODO!
3.5. Version 3: Charge decomposition#
The last alternative is to use a particle hole symmetry trial with the charge decomposition. Because we are focusing on positive \(U\) interactions, the charge decomposition will create a complex phase. At half-filling, due to the particle symmetry the walker weights will end up being all real and positive so we don’t have to worry about this
# TODO!