Two components in the Fermi Hubbard model

[Rafael Barfknecht]

Hello,

I am trying to run simulations for a two-component Fermi Hubbard model (calculating static properties such as densities and correlations).

Up to now, the two internal components to be different spins and label the oferators accordingly (such as as f_0 and f_1). However, this approach is returning strange results in the single-particle density matrix, even in a non-interacting case of only 2 particles (results are different if I calculate it for individual particles, one at a time).

So I considered calculating the same properties for a system of two fermions with the same total spin s=1/2, but with two different flavors. Then I labeled the operators as f_0_1 and f_1_1, where the first index labels the flavor, as indicated in the documentation.

However, I am getting the following error:

('Quantum number transfers:', [(0, 1), (1, 0), (0, 0), (0, 0)], [(1, 1), (1, 1), (0, 1), (1, 0)])
Traceback (most recent call last):
File "FermiHubbard.py", line 116, in <module>
MainFiles = mps.WriteFiles(parameters, Operators, H)
File "/usr/local/lib/python2.7/dist-packages/MPSPyLib/tools.py", line 914, in WriteFiles
p['unique_ID'], PostProcess)
File "/usr/local/lib/python2.7/dist-packages/MPSPyLib/Ops.py", line 439, in write
self.OperatorstoQOperators()
File "/usr/local/lib/python2.7/dist-packages/MPSPyLib/Ops.py", line 659, in OperatorstoQOperators
" the manual for more detail")
Exception: Operator fsplus is not injective! See the manual for more detail</module>

Can anyone help? Code is attached.

Best regards,

Rafael

import MPSPyLib as mps
import numpy as np
import sys
import time

t0=time.time()

#number of sites
L=128
#size of system
ltrap=5.
indexsite=np.linspace(0,L-1,L)
xgrid=np.linspace(-ltrap,ltrap,L)
#lattice spacing
a=2*ltrap/L
#hopping parameter
t=0.5/(a*a)
#particle number
nup=1
ndw=1
#interaction strength
g=0.0
U=g/a

#offset energy
#s=np.cos(np.pi/(L+1.))/(a*a)
s=1.0/(a*a)
scale=s

#single-well parameters
w=1.0

#double-well parameters
mu=1.0
a1=2.0
a2=2.0
om1=1.0
om2=1.0
C=1.5
D=0.8
a=0.3765
om=1.0644
x1=-0.7376
x2=1.1376

#choice of potential (dw=0 for single-well, dw=1 for double-well)
dw=0

if dw == 0:
   def trap(x):
      return 0.5*w*w*x*x

else:
   def trap(x):
      conds = [x < x1, (x > x1) & (x < x2), x > x2]
      funcs = [lambda x: 0.5*mu*om1*om1*(x+a1)*(x+a1),lambda x: -0.5*mu*om*om*(x-a)*(x-a)+C,lambda x: 0.5*mu*om2*om2*(x-a2)*(x-a2)+D]
      return np.piecewise(x, conds, funcs)

trappot=trap(xgrid)

#print('pot',trappot)

file=open('potential.dat','w')
for j in indexsite:
    l=int(j)
    print>>file, xgrid[l],trappot[l]
file.close()

#Build operators
Operators = mps.BuildFermiOperators(Nmax=2,Nmin=0,spin=0.5)
Operators['interaction']=np.dot(Operators['nf_0'],Operators['nf_1'])

# Define Hamiltonian MPO
H = mps.MPO(Operators)
H.AddMPOTerm('bond', ['fdagger_0', 'f_0'], hparam='t', weight=-1.0, Phase=True)
H.AddMPOTerm('bond', ['fdagger_1', 'f_1'], hparam='t', weight=-1.0, Phase=False)
H.AddMPOTerm('site', 'interaction', hparam='U', weight=1.0)
H.AddMPOTerm('site', 'nf_0', hparam='trap', weight=1.0)
H.AddMPOTerm('site', 'nf_1', hparam='trap', weight=1.0)
H.AddMPOTerm('site', 'nf_0', hparam='s', weight=1.0)
H.AddMPOTerm('site', 'nf_1', hparam='s', weight=1.0)

# Observables
myObservables = mps.Observables(Operators)
myObservables.AddObservable('site', 'nf_0', 'number_up')
myObservables.AddObservable('site', 'nf_1', 'number_dw')
myObservables.AddObservable('corr', ['fdagger_0', 'f_0'], 'spdm_up')
myObservables.AddObservable('corr', ['fdagger_1', 'f_1'], 'spdm_dw')

# Convergence data
myConv = mps.MPSConvParam(variance_tol=1e-5,max_bond_dimension=500,max_num_sweeps=30)
# Define statics
parameters = [{ 
    # Directories
    'simtype':'Finite',
    'job_ID':'FermiHUbbard_',
    'unique_ID':'L_'+str(L)+'_nup_'+str(nup)+'_ndw_'+str(ndw)+'_g_'+str(g),
    'Write_Directory':'TMP/', 
    'Output_Directory':'OUTPUTS/', 
    # System size and Hamiltonian parameters
    'L':L,
    't':t,
    'trap':trappot,
    'U':U,
    's':scale,
    'verbose':2,
    'logfile':True,
    # Specification of symmetries and good quantum numbers
    'Abelian_generators':['nf_0','nf_1'],
    'Abelian_quantum_numbers':[nup,ndw],
    'MPSObservables':myObservables,
    'MPSConvergenceParameters':myConv
}]

# Write Fortran-readable main files
MainFiles = mps.WriteFiles(parameters, Operators, H)

mps.runMPS(MainFiles, RunDir='')

Outputs = mps.ReadStaticObservables(parameters)

upnum = Outputs[0]['number_up']
dwnum = Outputs[0]['number_dw']
energ = Outputs[0]['energy']
convr = Outputs[0]['converged']
spdmu = Outputs[0]['spdm_up']
spdmd = Outputs[0]['spdm_dw']

t1=time.time()

print('total energy',energ)
print('convergence', convr)
print('elapsed time',t1-t0)

#file=open(str(nup)+'+'+str(ndw)+'_density_maj_dw_'+str(dw)+'_g_'+str(g)+'.dat','w')
#for j in indexsite:
#   l=int(j)
#   print>>file, xgrid[l],upnum[l]/(a*nup)
#file.close()

#file=open(str(nup)+'+'+str(ndw)+'_density_min_dw_'+str(dw)+'_g_'+str(g)+'.dat','w')
#for j in indexsite:
#   l=int(j)        
#   print>>file, xgrid[l],dwnum[l]/(a*ndw)
#file.close()

#file=open(str(nup)+'+'+str(ndw)+'_energy_dw_'+str(dw)+'_g_'+str(g)+'.dat','w')
#print>>file, energ
#file.close()

file=open(str(nup)+'+'+str(ndw)+'_spdm_maj_dw_'+str(dw)+'_g_'+str(g)+'.dat','w')
for j in indexsite:
   for k in indexsite:    
      l=int(j)
      m=int(k)     
      print>>file, xgrid[l],xgrid[m],spdmu[l][m]/(a*nup)
file.close()

file=open(str(nup)+'+'+str(ndw)+'_spdm_min_dw_'+str(dw)+'_g_'+str(g)+'.dat','w')
for j in indexsite:
   for k in indexsite:    
      l=int(j)
      m=int(k)     
      print>>file, xgrid[l],xgrid[m],spdmd[l][m]/(a*ndw)
file.close()