National High Magnetic Field Lab, Tallahassee, Florida
Theory Winter school, January 9-13, 2012

Anders W. Sandvik, Boston University

Lecture slides Lecture 1 Lecture 2

Detailed review of numerical methods for quantum spin systems (including SSE)
Computational Studies of Quantum Spin Systems, A. W. Sandvik, AIP Conf. Proc. 1297, 135 (2010)

Material for afternoon tutorial, January 11

A basic variant of the Stochastic Series Exapsnion (SSE) Quantum Monte Carlo program is provided here. It computes some essential properties if the S=1/2 Heisenberg model in one and two dimensions.

The lattice size N=Lx*Ly, with the lengths Lx and Ly given by the user as input. In cases where Lx and Ly are not the same, Lx should be the longer size. For Ly=1, the system is then a one-dimensional chain, with periodic boundaries in the x direction. For Ly=2 it is a 2-leg ladder, with periodic boundary conditions in the x direction but not in the y direction (since in this case that would double the effective interaction strength in the y direction). For all other cases the boundaries are periodic in both directions. The figure below summarizes the lattices and boundary conditions.

Note that in all cases except the single chain, Lx and Ly should be even numbers (Lx should be even also for the Ly=1 chain). Otherwise there is a sign problem in the simulations, which is not taken into account in this program.

Running instructions

Compile the programs (e.g., with g95 or gfortran).

Input parameter for ssebasic are read from a file 'read.in', in the following format:

```    lx, ly, beta
nbins, msteps, isteps

lx,ly  : lattice size in x and y direction
beta   : inverse dimensionless temperature J/T
nbins  : Number of bins (averages written to file 'res.dat' after each bin
msteps : Number of MC sweeps in each bin (measurements after each sweep)
isteps : Number of MC sweeps for equilibration (no measurements)
```

Output is written after each completed bin to 'res.dat'. Thiese raw data should be processed by the program 'res.f90'. It produces two files:

- e.dat: internal energy and specific heat (with error bars)

- m.dat: squared sublattice magnetization and uniform susceptibility (with error bars)

In addition, the simulation program also produces a file 'cut.dat' in which the new truncation of the series expansion is written each time it is increased.

An example

Running the program with these oarameters in 'read.in'

```    8   8      2.
10  10000  10000
```
Produced the following output (results of 10 bins) in 'res.dat' (after running a few seconds on a MacBook Air):
```    -0.60545312499999993       0.45663118750002241       0.12939503187647464       7.49062499999999937E-002
-0.60682343749999990       0.44244269437507455       0.12865359015698782       7.48968750000000016E-002
-0.60736640624999993       0.42738905609377298       0.12857748019236501       7.65843749999999962E-002
-0.60639140624999999       0.47596175609373859       0.13036687638192762       7.63156249999999980E-002
-0.60511015625000009       0.39663114359370866       0.12531685491597871       7.72062500000000040E-002
-0.60642031250000006       0.42985882437494638       0.13028933087043010       7.74468749999999984E-002
-0.60657500000000009       0.52802733499993337       0.13139164458652658       7.64906250000000065E-002
-0.60733203124999990       0.42871590234381074       0.12873853782797756       7.75593749999999998E-002
-0.60921562500000004       0.53947162249994562       0.13185469109301354       7.49156249999999996E-002
-0.60873593750000010       0.53154232437492510       0.12912797510686908       7.72593750000000051E-002
```
Running 'res.f90' gave this output in e.dat (energy and specific heat)
```    -0.60694234     0.00040847     0.46566718     0.01608355
```
and in 'm.dat' (sublattice magnetization squared and susceptibility):
```     0.12937120     0.00057857     0.07635813     0.00034245
```
Note that the specific heat is much noisier (has much larger error bar) than the energy. It is very difficult (requires very long runs to get small error bars) to compute the specific heat at low temperatures (roughly beta>4)

Test of the correctness of the program

It is important to test QMC programs against exact results, whenever possible. It is always possible to check against exact diagonalization for small systems. As an example, for a chain of L=32 the ground state energy per spin is -0.4439539825. Let's see if we can get agreement with this number by choosing beta=1/T large enough (T should be much smaller than the finite-size gap, which scales approximately as 1/L). These are results versus beta using 20 bins of 100000 MC steps each:

```    1.0    -0.20477922     0.00016893
2.0    -0.34135748     0.00016988
4.0    -0.41897368     0.00007303
8.0    -0.43776850     0.00003716
16.0    -0.44240240     0.00002830
32.0    -0.44377160     0.00002360
64.0    -0.44394747     0.00002178

exact   -0.44395398
```
The last two data points agree with the exact results within error bars (here in the fifth digit). Note that beta was successively doubled, which is useful when one wants to check the convergence in calculations aimed at the ground state. Note also that the error bar decreases with increasing beta. This is because the 'imaginary time' dimension of the system increases, and the energy (estimated using the average length of the expansion) is automatically averaged over the whole space-time volume. For the largest beta this calculation took about 5 minutes on a MacBook Air.