Hi Marco,

I am too busy this weekend to help find the problem but I have a suggestion: can you first try your trick for antiperiodic bc on a simple chain of noninteracting and even spinless fermions?

Matthias

On May 17, 2013, at 19:12, Marco Di Liberto m.f.diliberto@uu.nl wrote:

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