Dear developers,
I have recently discovered this very useful library and I tried to implement the t-J model to get used to the way the library works. I have been trying to reproduce the results in the papers:

- Troyer, M. and Tsunetsugu, H. and Rice, T. M. and Riera, J. and Dagotto, E., Spin gap and superconductivity in the one-dimensional t-J model with Coulomb repulsion, PhysRevB.48.4002 (1993)
- Ogata, Masao and Luchini, M. U. and Sorella, S. and Assaad, F. F., Phase diagram of the one-dimensional t-J model, PhysRevLett.66.2388 (1991)

but I encountered some problems in the calculation of Luttinger parameters that I am not sure where they are coming from.

In the papers I have read that when the (even) number of particles N is a multiple of 4, one takes anti-periodic boundary conditions and periodic in the other case (N is thus in the form N=4m or N=4m+2, where m is an integer). As far as I understand periodic bc are automatically implemented, thus, to perform the antiperiodic bc, I thought that this is equivalent to a flux in the boundary term exp(iπ); I can rewrite the model absorbing this flux into complex hopping coefficients t ---> t * exp(iπ/L) and keep the periodic bc then so that I can use the library with the translational symmetry in action. The python script is

import pyalps import numpy as np import math

Pi = math.pi bc = 'ape' L = 16 J = 2 phi = Pi/L Npt = [2, 4, 6] parms = [] #prepare the input parameters for N in Npt: for K in [0, 1]: if N%4 == 0: bc = 'ape' else: bc = 'pe' if bc == 'ape': tr = 'exp(' + str(phi)+'*I)' tl = 'exp(-' + str(phi)+'*I)' elif bc == 'pe': tr = 1 tl = 1 else: print 'Boundary conditions not defined properly' parms.append({ 'MODEL_LIBRARY' : "my_models.xml", 'LATTICE' : "chain lattice", 'MODEL' : "t-J complex", 'tr' : tr, 'tl' : tl, 'J' : J, 'L' : L, 'N_total' : N, 'Sz_total' : 0, 'TOTAL_MOMENTUM' : 2*Pi*K/L, 'CONSERVED_QUANTUMNUMBERS' : 'N,S,Sz', #'MEASURE_STRUCTURE_FACTOR[Structure Factor S]' : 'Sz', #'MEASURE_CORRELATIONS[Diagonal spin correlations]=' : 'Sz', #'MEASURE_CORRELATIONS[Offdiagonal spin correlations]' : 'Splus:Sminus' })
As you can see I diagonalize both in the sector with K = 0 and K = 2π/L (to calculate charge velocity). I have modified the file my_models.xml for the t-J as follows:

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

<HAMILTONIAN name="t-J complex">
  <PARAMETER name="mu" default="0"/>
  <PARAMETER name="tr" default="1"/>
  <PARAMETER name="tl" default="1"/>
  <PARAMETER name="J" default="1"/>
  <PARAMETER name="V" default="0"/>
  <PARAMETER name="t'" default="0"/>
  <PARAMETER name="J'" default="0"/>
  <PARAMETER name="V'" default="0"/>
  <PARAMETER name="t''" default="0"/>
  <PARAMETER name="J''" default="0"/>
  <PARAMETER name="V''" default="0"/>
  <PARAMETER name="t0" default="t"/>
  <PARAMETER name="t1" default="t'"/>
  <PARAMETER name="t2" default="t''"/>
  <PARAMETER name="V0" default="V"/>
  <PARAMETER name="V1" default="V'"/>
  <PARAMETER name="V2" default="V''"/>
  <PARAMETER name="J0" default="J"/>
  <PARAMETER name="J1" default="J'"/>
  <PARAMETER name="J2" default="J''"/>
  <BASIS ref="t-J"/>
  <SITETERM site="i">
    <PARAMETER name="mu#" default="mu"/>
    -mu#*n(i)
  </SITETERM> 
  <BONDTERM source="i" target="j">
    <PARAMETER name="V#" default="0"/>
    <PARAMETER name="J#" default="0"/>
    -tr*fermion_hop_r(i,j) -tl*fermion_hop_l(i,j) + J#*exchange(i,j)+(V#-J#/4)*n(i)*n(j)
  </BONDTERM>
</HAMILTONIAN>

Here the printed data for N_total = 4

[[x=[0] y=[-7.5655955] props={'observable': 'Energy', 'tl': 'exp( - 0.19634954084899999827 * I)', 'MODEL_LIBRARY': 'my_models.xml', 'Sz_total': 0.0, 'J': 2.0, 'tr': 'exp(0.19634954084899999827 * I)', 'TOTAL_MOMENTUM': 0.0, 'filename': './parm1a.task3.out.h5', 'LATTICE': 'chain lattice', 'hdf5_path': '/spectrum/sectors/0/results/Energy', 'SEED': 1318589528.0, 'CONSERVED_QUANTUMNUMBERS': 'N,S,Sz', 'N_total': 4.0, 'MODEL': 't-J complex', 'L': 16.0}]], [[x=[0] y=[-7.31510355] props={'observable': 'Energy', 'tl': 'exp( - 0.19634954084899999827 * I)', 'MODEL_LIBRARY': 'my_models.xml', 'Sz_total': 0.0, 'J': 2.0, 'tr': 'exp(0.19634954084899999827 * I)', 'TOTAL_MOMENTUM': 0.392699081699, 'filename': './parm1a.task4.out.h5', 'LATTICE': 'chain lattice', 'hdf5_path': '/spectrum/sectors/0/results/Energy', 'SEED': 1587024985.0, 'CONSERVED_QUANTUMNUMBERS': 'N,S,Sz', 'N_total': 4.0, 'MODEL': 't-J complex', 'L': 16.0}]]

What happens is that the curve for v_c is not correct (here below the values)

N v_c

2.0 0.387678357479 4.0 0.637872515084 6.0 0.946138404787 8.0 0.932824458966 10.0 1.07412501501 12.0 0.625143469145 14.0 0.608197205496

and for the case N = 4m (with anti-perdiodic bc) the compressibility is negative. It looks like the values for N = 4m+2 are correct though (by checking the PRL by Ogata et al. cited above). Therefore, probably the case N=4m with anti-periodic bc has some mistake in my coding. I am not sure if what is not working is the hamiltonian definition, the way I implemented the complex coefficients or I don't know what else. Maybe the momentum definition when I take k = 2π/L? In the online documentation I found that this can be a source of mistake, but it's not clear to me how to do it correctly in the python script. From the documentation it looks like one has to pass a string but I get an error, while in this way 'TOTAL_MOMENTUM' : 2*Pi*K/L no error occurs.
What is weird is that in the case of complex coefficients it is given a string in the dictionary of parameters, while for real hopping parameters one gives just a number. I don't get any error when I run, so I thought that the code takes into account these two possibilities. Moreover, I am not sure that the string definition of the parameters 'tr' and 'tl' is correct. I had to invent this fancy way of combining different strings to put the parameters π and L in the phase of the hopping, but I get no error in the run and I also find it printed in the data list correctly as you can check.
As a remark, in my runs I found that the ground state is always in the K = 0 sector when I follow the boundary condition prescription, in agreement with the papers above.

I hope you can give some hint of what is going on. Thank you in advance!

---
Marco Di Liberto
PhD Student

ITP (Institute for Theoretical Physics)
Utrecht University