Dear ALPS Developers and Users,

I am trying to implement Fermi-Hubbard model in saw-tooth lattice and use
ALPS MPS for further calculations. However, I am not able to get the correct
results. I have checked my custom lattice using printgraph, also it seems to
work fine with the Bose-Hubbard model. It will be really helpful if someone can
look into the scripts below and let me know if there is a mistake? Or if there is
some issue with ALPS MPS itself?

parameters file:
--------------------

import pyalps

parms = []
p = dict()
p['SWEEPS'            ] = 10
p['NUMBER_EIGENVALUES'] = 1
p['optimization'      ] = 'twosite'
p['MAXSTATES'         ] = 1200
p['TRUNCATION'        ] = 1.e-10
p['init_state'        ] = 'default'
p['LATTICE_LIBRARY'   ] = 'mylattices.xml'
p['LATTICE'           ] = 'sawtooth open chain lattice'
p['MODEL_LIBRARY'     ] = 'mymodel.xml'
p['MODEL'             ] = 'fermion Hubbard'
p['L'                 ] = 10
p['t0'                ] = -1.4142135623731
p['t1'                ] = -1.0
p['t2'                ] = -1.4142135623731
p['U'                 ] = 0.0
p['Nup_total'         ] = 5
p['Ndown_total'       ] = 5
p['COMPLEX'           ] = 1
p['CONSERVED_QUANTUMNUMBERS'           ] = 'Nup,Ndown'
p['MEASURE[EnergyVariance]'            ] = 1
p['MEASURE_LOCAL[n-up]'                ] = 'n_up'
p['MEASURE_LOCAL[n-down]'              ] = 'n_down'
p['storagedir'] = 'storage'
parms.append(p)

## write the input file and run the simulation
infiles=pyalps.writeInputFiles('sim',parms)
res=pyalps.runApplication('mps_optim',infiles,writexml=True)
---------------------------------------------------------------------------------------------

lattice file
--------------

<LATTICES>
<LATTICE name="chain lattice" dimension="1">
  <PARAMETER name="a" default="1"/>
  <BASIS><VECTOR>a</VECTOR></BASIS>
  <RECIPROCALBASIS><VECTOR>2*pi/a</VECTOR></RECIPROCALBASIS>
</LATTICE>

<UNITCELL name="sawtooth" dimension="1">
  <VERTEX id="1" type="0"></VERTEX>
  <VERTEX id="2" type="1"></VERTEX>
  <EDGE type="0"><SOURCE vertex="1" offset="0"/><TARGET vertex="2" offset="0"/></EDGE>
  <EDGE type="1"><SOURCE vertex="1" offset="0"/><TARGET vertex="1" offset="1"/></EDGE>
  <EDGE type="0"><SOURCE vertex="2" offset="0"/><TARGET vertex="1" offset="1"/></EDGE>
</UNITCELL>

<LATTICEGRAPH name = "sawtooth open chain lattice" vt_skip="true">
 <FINITELATTICE>
   <LATTICE ref="chain lattice"/>
   <EXTENT dimension="1" size ="L"/>
   <BOUNDARY type="open"/>
 </FINITELATTICE>
 <UNITCELL ref="sawtooth"/>
</LATTICEGRAPH>

</LATTICES>

------------------------------------------------------------------------------------------------------------
model file
-------------

<MODELS>
<SITEBASIS name="fermion">
  <QUANTUMNUMBER name="Nup" min="0" max="1" type="fermionic"/>
  <QUANTUMNUMBER name="Ndown" min="0" max="1" type="fermionic"/>
  <OPERATOR name="Splus" matrixelement="1">
    <CHANGE quantumnumber="Ndown" change="-1"/>
    <CHANGE quantumnumber="Nup" change="1"/>
  </OPERATOR>
  <OPERATOR name="Sminus" matrixelement="1">
    <CHANGE quantumnumber="Nup" change="-1"/>
    <CHANGE quantumnumber="Ndown" change="+1"/>
  </OPERATOR>
  <OPERATOR name="Sz" matrixelement="(Nup-Ndown)/2"/>
  <OPERATOR name="Nup" matrixelement="Nup"/>
  <OPERATOR name="Ndown" matrixelement="Ndown"/>
  <OPERATOR name="c_down" matrixelement="1">
    <CHANGE quantumnumber="Ndown" change="-1"/>
  </OPERATOR>
  <OPERATOR name="cdag_down" matrixelement="1">
    <CHANGE quantumnumber="Ndown" change="1"/>
  </OPERATOR>
  <OPERATOR name="c_up" matrixelement="1">
    <CHANGE quantumnumber="Nup" change="-1"/>
  </OPERATOR>
  <OPERATOR name="cdag_up" matrixelement="1">
    <CHANGE quantumnumber="Nup" change="1"/>
  </OPERATOR>
  <OPERATOR name="n" matrixelement="Nup+Ndown"/>
  <OPERATOR name="n_up" matrixelement="Nup"/>
  <OPERATOR name="n_down" matrixelement="Ndown"/>
  <OPERATOR name="double_occupancy" matrixelement="Nup*Ndown"/>
</SITEBASIS>

<BASIS name="fermion">
  <SITEBASIS ref="fermion"/>
  <CONSTRAINT quantumnumber="Nup" value="Nup_total"/>
  <CONSTRAINT quantumnumber="Ndown" value="Ndown_total"/>
</BASIS>

<SITEOPERATOR name="fielddag_ud" site="x">
  cdag_up(x)*cdag_down(x)
</SITEOPERATOR>
<SITEOPERATOR name="field_du" site="x">
  c_down(x)*c_up(x)
</SITEOPERATOR>

<SITEOPERATOR name="nsqu" site="x">
  n(x)*n(x)
</SITEOPERATOR>
<SITEOPERATOR name="Szsqu" site="x">
  Sz(x)*Sz(x)
</SITEOPERATOR>

<BONDOPERATOR name="fermion_hop" source="x" target="y">
  cdag_up(x)*c_up(y)+cdag_up(y)*c_up(x)+cdag_down(x)*c_down(y)+cdag_down(y)*c_down(x)
</BONDOPERATOR>

<HAMILTONIAN name="harm fermion Hubbard">
  <PARAMETER name="mu" default="0"/>
  <PARAMETER name="t" default="1"/>
  <PARAMETER name="V" default="0"/>
  <PARAMETER name="t'" default="0"/>
  <PARAMETER name="V'" default="0"/>
  <PARAMETER name="U" default="0"/>
  <PARAMETER name="t0" default="t"/>
  <PARAMETER name="t1" default="t'"/>
  <PARAMETER name="V0" default="V"/>
  <PARAMETER name="V1" default="V'"/>
  <PARAMETER name="K" default="0"/>
  <BASIS ref="fermion"/>
  <SITETERM site="i">
    <PARAMETER name="mu#" default="mu"/>
    <PARAMETER name="U#" default="U"/>
    <PARAMETER name="K#" default="K"/>
    -mu#*n(i)+U#*n_up(i)*n_down(i)+K#*n(i)*(x-0.5*L)^2
  </SITETERM>
  <BONDTERM source="i" target="j">
    <PARAMETER name="t#" default="0"/>
    <PARAMETER name="V#" default="0"/>
    -t#*fermion_hop(i,j) + V#*n_up(i)*n_up(j)
  </BONDTERM>
</HAMILTONIAN>

<HAMILTONIAN name="fermion Hubbard">
  <PARAMETER name="mu" default="0"/>
  <PARAMETER name="t" default="1"/>
  <PARAMETER name="V" default="0"/>
  <PARAMETER name="t'" default="0"/>
  <PARAMETER name="V'" default="0"/>
  <PARAMETER name="U" default="0"/>
  <PARAMETER name="t0" default="t"/>
  <PARAMETER name="t1" default="t'"/>
  <PARAMETER name="V0" default="V"/>
  <PARAMETER name="V1" default="V'"/>
  <BASIS ref="fermion"/>
  <SITETERM site="i">
    <PARAMETER name="mu#" default="mu"/>
    <PARAMETER name="U#" default="U"/>
    -mu#*n(i)+U#*n_up(i)*n_down(i)
  </SITETERM>
  <BONDTERM source="i" target="j">
    <PARAMETER name="t#" default="0"/>
    <PARAMETER name="V#" default="0"/>
    -t#*fermion_hop(i,j) + V#*n(i)*n(j)
  </BONDTERM>
</HAMILTONIAN>
</MODELS>
-----------------------------------------------------------------------------------------------------------

Best regards-
Manpreet