Forcing S=0 and total number of electrons in the Fermi-Hubbard model (Julia)

Question

I'm trying to adapt the code here ( http://itensor.org/docs.cgi?page=classes/autompo&vers=cppv3 )from C++ to Julia, specifically for the 1D Hubbard model. I'm a little confused on how I can enforce the number of electrons, number of sites, and number of spin-up and spin-down electrons.

Here is the code:

N = 2 # number of sites
t = 1.0 # hopping term
u = 0.0 # potential term (onsite)

sites = siteinds("Electron", N) 

ampo = AutoMPO()

for i=1:N
    ampo += u ,"Nupdn", i
end

for i=1:N-1
    ampo += -t, "Cdagup", i, "Cup", i+1
    ampo += -t, "Cdagup", i+1, "Cup", i
    ampo += -t, "Cdagdn", i, "Cdn", i+1
    ampo += -t, "Cdagdn", i+1, "Cdn", i
end

H = MPO(ampo, sites)

N = 2 above enforces that we have a two site Hubbard model. I'm assuming that I've also enforced two electrons here in the line:

sites = siteinds("Electron", N)

since I'm creating N=2 electrons (though I'm not sure if this is entirely true).

How might I enforce total S=0 (half of the electrons spin up, the other half spin down)? And how do I control how many electrons I have?

Apologies for the basic question!

Answer 1

No problem! This is a common question since the way the total quantum numbers (QNs), such as particle number, is fixed is kind of subtle in ITensor. It is fixed by your choice of initial state. So if you prepare the MPS you are using as an initial guess for DMRG (I assume you are planning to do a DMRG calculation) to be a state with a certain particle number, then that particle number will remain the same as long as you construct your MPS using site indices that carry QN information.

The line of code

sites = siteinds("Electron",N)

is actually creating an array of N objects of type Index, so N electron sites, not N electrons.

To make these sites carry QN information, you need to change this line to:

sites = siteinds("Electron",N,conserve_qns=true)

Then, to prepare an initial state with two electrons, say one being an up-spin electron and the other being a down-spin electron, you can use code such as:

states = ["0" for n=1:N]
states[1] = "Up"
states[2] = "Dn"
psi = productMPS(sites,states)

You can also call psi = randomMPS(sites,states,10) to create a randomized MPS which has the site indices sites, the same total QN as given by the product state states, and a bond dimension of 10.

Then if you plug this psi MPS into DMRG as an initial state together with a Hamiltonian that conserves particle number, the way ITensor works is that the particle number is guaranteed not to change. This is because when you use Index objects that have QNs in them, only operations which change these QNs by a definite amount are allowed, and your Hamiltonian will change them by an amount which is zero if it conserves them.

energy,psi = dmrg(H,psi,sweeps)

Please post a comment if you have any follow-up questions.

There is a full example of doing DMRG for the extended Hubbard model with QN conservation here:
https://github.com/ITensor/ITensors.jl/blob/master/examples/dmrg/exthubbard.jl

Best,
Miles

Comments

Ah I didn't even think to look at QNs! Thanks a bunch!