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(a)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(a)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
>
>