Transverse Field Ising model with QMC

[ShHsLin-quantinuum]

Hi all,

Thanks for the effort of keeping the package running.
I was trying to get some benchmark data of TFI energy in 3D. But to test whether my setup is correct, I first try TFI in 2D of size (4,4) with the code below modified from the tutorial.

import numpy as np
import pyalps


params = []
h = 3.1
L = 4
for T in [0.08, 0.04,]:
    params.append({
        # lattice parameters
        'LATTICE'        : "square lattice",
        'L'              : L,

        # model parameters
        'MODEL'          : "spin",
        'local_S'        : 0.5,
        'J0'              : 4,
        'J1'              : 0,
        'J2'              : 0,
        'h'              :  2*h,
        # ensemble parameter
        'T'              : T,

        # QMC parameters
        'THERMALIZATION' : 2000,
        'SWEEPS'         : 5000,
    })
prefix = 'qmc_square_lattice'


# Write the input files for the ALPS codes.
input_file = pyalps.writeInputFiles(prefix, params)

# The following command runs the applications.
res = pyalps.runApplication('dirloop_sse',input_file,Tmin=5)

However, the energy density I got -3.154 is significantly different from what I got from ED -3.303 on PBC lattice, also it doesn't match the OBC lattice energy -3.231.

I tried varying T, THERMALIZATION, SWEEPS, and other parameters. But it seems like I am not getting the correct values. I guess there is some setup that I did wrong here. Any suggestions is highly appreciated.

[Ooolab]

One possibility is that you haven't simulate the system at low enough temperature. Your ED result is lower in energy compared with QMC simulations, which is what we should have. So, could you try several lower temperatures and see if the result approaches your ED calculation?

[ShHsLin-quantinuum]

No, I would guess that is not the culprit. I did already try actually 0.08, 0.04, 0.02 with higher number of

THERMALIZATION=20000
SWEEPS=10000

and I got

Energy:  [-3.153102496606304 +/- 0.0014179427292937712
 -3.156207199734865 +/- 0.0015464035856224017
 -3.1560331920582474 +/- 0.0017122601435828344]

[LodePollet]

I am confused by what you are asking. Do you want to study the transverse field Ising model? In that case, what is the value of Gamma? And you should not use dirloop_sse for that but loop. It seems to me that you are effectively simulating another model, probably an XY-type of model in a magnetic field.
Assuming you want to study the 2d TFIM, I ran the following input file and obtained a ground state energy close to -33.2 (extensive). But as I said, I am not sure what you want to simulate exactly. Hope this helps.

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
import numpy as np
Ls = [4]
betas = np.linspace (1.0 , 20.0 , 21)
parms = []
for beta in betas :
  for L in Ls :
    parms.append (
      {
        'LATTICE'               : "square lattice",
        'MODEL'                 : "spin",
        'local_S'               : 0.5,
        'L'                     : L,
        'Jxy'                   : 0,
        'Jz'                    : 4,
        'h'                     : 6.2,
        'MEASURE[Correlations]' : 1,
        'Gamma'                 : 1.0,
        'BETA'                  : beta,
        'THERMALIZATION'        : 10000,
        'SWEEPS'                : 10000,
        'ALGORITHM'             : 'loop'
      }
    )
input_file = pyalps.writeInputFiles('2D-TFIM', parms)
res = pyalps.runApplication('loop', input_file, Tmin =1)

data=pyalps.loadMeasurements(pyalps.getResultFiles( prefix='2D-TFIM'), 'Energy')
energy = pyalps.collectXY (data, x='BETA', y ='Energy' ,foreach =['L'])

plt.figure ()
pyalps.plot.plot(energy)
plt.xlabel('Beta')
plt.ylabel('Energy')
plt.legend()
plt.show ()

[ShHsLin-quantinuum]

Yes, I want to do 2D-TFIM. Thanks a lot! I leave the details comment below in case in might be helpful for others.

Following your script, I could now produce the correct energy (ED -3.303):

Energy:  [-3.3039098750000013 +/- 0.0016651735712997654
 -3.3032742499999985 +/- 0.0021627242144303988
 -3.2954909999999997 +/- 0.0028834497341627507]
 Ts: [0.02, 0.04, 0.08]

under the setup

        'Jxy'             : 0,
        'Jz'              : 4,
        'Gamma'          : 2*h,
        'ALGORITHM'      : "loop",

The factor 2 and 4 are just the typical convention conversion between of spin-1/2 and paulis, where the above implementation convention gives H = ZZ + h X

Previously, I setup the script differently because I follow the Quantum Monte Carlo tutorial MC-09. Since the example is Heisenberg with field so effectively H ~ XX +YY +ZZ + h Z. I thought I could directly modify the script and get result of TFIM by setting interaction of YY and ZZ to zero. (J1=J2=0). But I guess it is the wrong way to do what I want to do.