Dear Prof. Troyer,

thanks for you answer. Here what I have done: I was able to check that the t-J results for L=16 with periodic bc are numerically correct and what is wrong are the ones with anti-periodic bc (I extrapolated v_c from the graphs in the paper by Ogata and checked also Bethe-ansatz results for the ground state energy).

EXTRAPOLATION

N N/L v_c

2 0.125891 0.389187

4 0.252045 0.703023

6 0.376964 0.941431

MY NUMERICS

N N/L v_c

2 0.125891 0.387678

4 0.252045 0.637872

6 0.376964 0.946138

I tried running for non-interacting spinless fermions. I did not impose any constraint on the first place and I looked at all momentum sectors. Number of particles is always odd to achieve total momentum zero. I have found for this system a curious behavior. I put L=12. For N=5 I use antiperiodic bc and find that the zero momentum sector has lowest energy -7.20976982 but there is another sector with 'TOTAL_MOMENTUM': 3.66519142918 that has energy -7.20976947. Slightly larger than the former one. If I run for the same filling and periodic bc, the former feature disappears and the zero momentum state has energy -7.46410162. Same thing happens with N=3. I am not sure of what kind of things you wanted me to look for the spinless case.

I also ran a simulation for the Hubbard model with L=12 sites and U=0 and some simulation with L=16. I used the same choice as I was doing for the tJ: N=4m+2 periodic and N=4m anti-periodic. I extrapolated the expected value from the theoretical curves in the thermodynamic limit (Schulz).

CASE L=12

N v_c v_c(expected)

2 0.511745 0.64

4 0.988616 1.12

6 1.398114 1.53

As you can see, the trend in v_c is below the real curve but, contrary to the t-J, for different bcs there is the same systematic error more or less. However, only for antiperiodic I find negative compressibility. I saw in the paper by Schulz that a simulation on L=16 sites fits well the v_c curve, so I tried to increase the size.

CASE L=16

N v_c v_c(expected)

2 0.38767 0.48

4 0.38767 0.896

SAME v_c AS IN THE t-J?? IS IT A COINCIDENCE?

I can also perform the calculation for N=6 for L=16 but the diagonalization of a single sector lasts 15min on my laptop.. which is unexpectedly long. How come? Is it really so long? From the paper by Schulz I have seen that 20 years ago it was doable to diagonalize the Hubbard model for L=16 up to half-filling.. and now it looks like it's very time consuming even at low filling. Is this lowness library related?

Anyway from the calculation at L=16 that I have done it looks like that the points are not sitting on the theoretical curve, contrary to the paper by Schulz for the calculation of K_rho with 16 sites or to my former numerical calculation for the t-J of v_c (only for periodic bc).

---

Marco Di Liberto

PhD Student

ITP (Institute for Theoretical Physics)

Utrecht University

On 17 May 2013 19:35, Matthias Troyer <troyer@phys.ethz.ch> wrote:

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