Sample code to study Z2 gauge theory

julia

function main(Ns = 40, dim = 4)

    latt = ones(Int64, ntuple(i1 -> i1 <= dim ? Ns : dim, dim + 1))
    links = CartesianIndices(latt)

    function mvup(x::CartesianIndex, d)
        # key routine to optimize: KISS
        f1(i1) = x[i1]
        f2(i1) = x[i1] == Ns ? 1 : x[i1] + 1
        f(i1) = i1 != d ? f1(i1) : f2(i1)
        return CartesianIndex(ntuple(f, length(x)))

    end

    function mvdown(x::CartesianIndex, d)
        # key routine to optimize: KISS
        f1(i1) = x[i1]
        f2(i1) = x[i1] == 1 ? Ns : x[i1] - 1
        f(i1) = i1 != d ? f1(i1) : f2(i1)
        return CartesianIndex(ntuple(f, length(x)))

    end

    function coldstart()
        latt .= 1
    end

    function randomstart()
        latt .= rand([-1, 1])
    end

    function staplecal(link::CartesianIndex)

        x = CartesianIndex(link.I[1:end-1])
        d = link.I[end]

        # following M. Creutz
        # staples around link(1->4)
        #    dperp        6--5
        #    ^            |  |
        #    |            1--4
        #    |            |  |
        #    -----> d     2--3

        staplesum = 0.0
        for dperp = 1:length(x)
            if dperp != d

                # plaquette 1234
                x = mvdown(x, dperp)
                link1 = latt[x, dperp]
                link2 = latt[x, d]
                x = mvup(x, d)
                link3 = latt[x, dperp]
                x = mvup(x, dperp)
                # plaquette 4561
                link4 = latt[x, dperp]
                x = mvup(x, dperp)
                x = mvdown(x, d)
                link5 = latt[x, d]
                x = mvdown(x, dperp)
                link6 = latt[x, dperp]
                staplesum += link1 * link2 * link3 + link4 * link5 * link6

            end
        end
        return staplesum
    end

    function update(beta)

        action = 0.0
        for link in links

            staplesum = staplecal(link)
            # calculate the Boltzmann weight
            bplus = exp(beta * staplesum * 1)
            bplus /= (bplus + 1 / bplus)

            # the heatbath algorithm
            if rand() < bplus
                latt[link] = 1
                action += staplesum
            else
                latt[link] = -1
                action -= staplesum
            end
        end

        return 1.0 - action / Ns^dim / dim / 6.0
    end

    beta_arr1 = range(1, stop = 0, length = 20)
    beta_arr2 = range(0, stop = 1, length = 20)

    coldstart()
    # randomstart()
    println("#cold -> hot")
    for beta in beta_arr1
        println(beta, " ", update(beta))
    end

    # already hot
    println("#hot -> cold")
    for beta in beta_arr2
        println(beta, " ", update(beta))
    end

end

@time main()

python


import numpy as np
import itertools

# N = 20 is already slow
N = 20
latt = np.ones(N**4*4).reshape([N, N, N, N, 4])

rng = np.random.default_rng()


# utility
def moveup(xvec, d):
    #  xvec is mutable
    xvec[d] += 1
    xvec[d] = xvec[d] % N   # stays between [0, 1, ..., N-1]
    return None


def movedown(xvec, d):
    #  xvec is mutable
    xvec[d] -= 1
    xvec[d] = xvec[d] % N
    return None


def coldstart():
    latt[:] = 1
    return None


def randomstart():

    sites = itertools.product(range(N), range(N), range(N), range(N), range(4))

    for site in sites:

        spin = rng.integers(2)
        if spin == 0:
            spin = -1

        latt[site] = spin

    return None


def update(beta):

    sites = itertools.product(range(N), range(N), range(N), range(N), range(4))

    action = 0.

    for site in sites:

        *x, d = site

        # following M. Creutz
        # staples around link(1->4)
        #    dperp        6--5
        #    ^            |  |
        #    |            1--4
        #    |            |  |
        #    -----> d     2--3

        staplesum = 0.
        for dperp in range(4):
            if dperp != d:

                # plaquette 1234
                movedown(x, dperp)
                staple = latt[x[0], x[1], x[2], x[3], dperp]
                staple *= latt[x[0], x[1], x[2], x[3], d]
                moveup(x, d)
                staple *= latt[x[0], x[1], x[2], x[3], dperp]
                moveup(x, dperp)
                staplesum += staple

                # plaquette 4561
                staple = latt[x[0], x[1], x[2], x[3], dperp]
                moveup(x, dperp)
                movedown(x, d)
                staple *= latt[x[0], x[1], x[2], x[3], d]
                movedown(x, dperp)
                staple *= latt[x[0], x[1], x[2], x[3], dperp]
                staplesum += staple


        # calculate the Boltzmann weight
        bplus = np.exp(beta*staplesum)
        bminus = np.exp(-beta*staplesum)
        bplus = bplus/(bplus+bminus)

        # the heatbath algorithm
        r = rng.uniform()

        if r < bplus:
            latt[site] = 1
            action += staplesum
        else:
            latt[site] = -1
            action -= staplesum

    return 1. - action/N**4/4./6.


def main():

    # beta_c = 0.44

    beta_arr1 = np.linspace(1., 0., 40)
    beta_arr2 = np.linspace(0., 1., 40)

    coldstart()
    print('# cold -> hot')
    for beta in beta_arr1:
        action = update(beta)
        print(beta, action)

    # already hot
    print('# hot -> down')
    for beta in beta_arr2:
        action = update(beta)
        print(beta, action)

    return None


main()

fortran 90

program main
        use iso_fortran_env
        implicit none
        integer, parameter :: dp = REAL64

        integer, parameter :: N = 40
        integer, dimension(N,N,N,N,4) :: latt

        real(dp) :: beta
        real(dp) :: dbeta = 0.01
        real(dp) :: actionS
        integer :: i1

        ! coldstart
        latt = 1

        ! cold -> hot
        do i1 = 1, 100
                beta = 1. - i1*dbeta
                call sweep(N, beta, latt, actionS)
                print*, beta, actionS
        end do

        ! hot -> cold
        do i1 = 1, 100
                beta = i1*dbeta
                call sweep(N, beta, latt, actionS)
                print*, beta, actionS
        end do

        contains

        subroutine moveup(xvec, d, N)
                implicit none
                integer, intent(in) :: d
                integer, intent(in) :: N
                integer, dimension(4), intent(inout) :: xvec
                xvec(d) = xvec(d) + 1
                if(xvec(d) > N) then
                        xvec(d) = xvec(d) - N
                else
                end if
        end subroutine moveup

        subroutine movedown(xvec, d, N)
                implicit none
                integer, intent(in) :: d
                integer, intent(in) :: N
                integer, dimension(4), intent(inout) :: xvec
                xvec(d) = xvec(d) - 1
                if(xvec(d) < 1) then
                        xvec(d) = xvec(d) + N
                else
                end if
        end subroutine movedown

        subroutine sweep(N, beta, lat, actionS)
                
                use iso_fortran_env
                implicit none
                integer, parameter :: dp = REAL64

                integer, intent(in) :: N
                real(dp), intent(in) :: beta
                integer, dimension(N,N,N,N,4), intent(inout) :: lat
                real(dp), intent(out) :: actionS

                integer :: i1, i2, i3, i4, d, dperp
                real(dp) :: r_num
                real(dp), dimension(N,N,N,N,4):: rand_arr

                integer, dimension(4) :: xvec
                real(dp) :: staple
                real(dp) :: staplesum
                real(dp) :: bplus
                real(dp) :: bminus

                call random_seed()
                call random_number(rand_arr)

                actionS = 0.
                do i1 = 1, N
                do i2 = 1, N
                do i3 = 1, N
                do i4 = 1, N
                do d = 1, 4

                xvec = [ i1, i2, i3, i4 ]

                staplesum = 0.
                do dperp = 1, 4
                if (dperp /= d) then
                call movedown(xvec, dperp, N)
                staple = lat(xvec(1), xvec(2), xvec(3), xvec(4), dperp)
                staple = staple*lat(xvec(1), xvec(2), xvec(3), xvec(4), d)
                call moveup(xvec, d, N)
                staple = staple*lat(xvec(1), xvec(2), xvec(3), xvec(4), dperp)
                call moveup(xvec, dperp, N)
                staplesum = staplesum + staple

                staple = lat(xvec(1), xvec(2), xvec(3), xvec(4), dperp)
                call moveup(xvec, dperp, N)
                call movedown(xvec, d, N)
                staple = staple*lat(xvec(1), xvec(2), xvec(3), xvec(4), d)
                call movedown(xvec, dperp, N)
                staple = staple*lat(xvec(1), xvec(2), xvec(3), xvec(4), dperp)
                staplesum = staplesum + staple
                else
                end if
                end do

                bplus = exp(beta * staplesum)
                bminus = exp(-beta * staplesum)
                bplus = bplus / (bplus + bminus)

                r_num = rand_arr(i1, i2, i3, i4, d)
                !call random_number(r_num)

                if (r_num <= bplus) then
                        lat(i1, i2, i3, i4, d) = 1
                        actionS = actionS + staplesum/N**4/24.
                else
                        lat(i1, i2, i3, i4, d) = -1
                        actionS = actionS - staplesum/N**4/24.
                end if

                end do
                end do
                end do
                end do
                end do

                actionS = 1. - actionS

        end subroutine sweep

end program main

c (from M. Creutz)

/* Z_2 lattice gauge simulation */
/* Michael Creutz <creutz@bnl.gov>     */
/* http://thy.phy.bnl.gov/~creutz/z2.c */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* the lattice is of dimensions SIZE**4  */
#define SIZE 40
int link[SIZE][SIZE][SIZE][SIZE][4]; /* last index gives link direction */

/* utility functions */
void moveup(int x[],int d) {
  x[d]+=1;
  if (x[d]>=SIZE) x[d]-=SIZE;
  return;
}
void movedown(int x[],int d) {
  x[d]-=1;
  if (x[d]<0) x[d]+=SIZE;
  return;
}
void coldstart(){  /* set all links to unity */
  int x[4],d;
  for (x[0]=0;x[0]<SIZE;x[0]++)
    for (x[1]=0;x[1]<SIZE;x[1]++)
      for (x[2]=0;x[2]<SIZE;x[2]++)
        for (x[3]=0;x[3]<SIZE;x[3]++)
          for (d=0;d<4;d++)
	    link[x[0]][x[1]][x[2]][x[3]][d]=1;
  return;
}
/* for a random start: call coldstart() and then update once at beta=0 */

/* do a Monte Carlo sweep; return energy */
double update(double beta){
  int x[4],d,dperp,staple,staplesum;
  double bplus,bminus,action=0.0;
  for (x[0]=0; x[0]<SIZE; x[0]++)
    for (x[1]=0; x[1]<SIZE; x[1]++)
      for (x[2]=0; x[2]<SIZE; x[2]++)
        for (x[3]=0; x[3]<SIZE; x[3]++)
          for (d=0; d<4; d++) {
            staplesum=0;
            for (dperp=0;dperp<4;dperp++){
              if (dperp!=d){
                /*  move around thusly:
                    dperp        6--5
                    ^            |  |
                    |            1--4
                    |            |  |
                    -----> d     2--3  */
                /* plaquette 1234 */
                movedown(x,dperp);
                staple=link[x[0]][x[1]][x[2]][x[3]][dperp]
                  *link[x[0]][x[1]][x[2]][x[3]][d];
                moveup(x,d);
                staple*=link[x[0]][x[1]][x[2]][x[3]][dperp];
                moveup(x,dperp);
                staplesum+=staple;
                /* plaquette 1456 */
                staple=link[x[0]][x[1]][x[2]][x[3]][dperp];
                moveup(x,dperp);
                movedown(x,d);
                staple*=link[x[0]][x[1]][x[2]][x[3]][d];
                movedown(x,dperp);
                staple*=link[x[0]][x[1]][x[2]][x[3]][dperp];
                staplesum+=staple;
              }
	    }
            /* calculate the Boltzmann weight */
            bplus=exp(beta*staplesum);
            bminus=1/bplus;
            bplus=bplus/(bplus+bminus);
            /* the heatbath algorithm */
            if ( drand48() < bplus ){
              link[x[0]][x[1]][x[2]][x[3]][d]=1;
              action+=staplesum;
            }
            else{
              link[x[0]][x[1]][x[2]][x[3]][d]=-1;
              action-=staplesum;
            }
          }
  action /= (SIZE*SIZE*SIZE*SIZE*4*6);
  return 1.-action;
}

/******************************/
int main(){
  double beta, dbeta, action;
  srand48(1234L);  /* initialize random number generator */
  /* do your experiment here; this example is a thermal cycle */
  dbeta=.01;
  coldstart();
  /* heat it up */
  for (beta=1; beta>0.0; beta-=dbeta){
    action=update(beta);
    printf("%g\t%g\n",beta,action);
  }
  printf("\n\n");
  /* cool it down */
  for (beta=0; beta<1.0; beta+=dbeta){
    action=update(beta);
    printf("%g\t%g\n",beta,action);
  }
  printf("\n\n");
  exit(0);
}

cpp

// Z2 4D lattice
// cpp adaptation of M. Creutz C code
//

# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# include <iostream>
# include <random>

using namespace std;

// the lattice is of dimensions SIZE**4
#define SIZE 40

class config {

  private:
  double lattice[SIZE][SIZE][SIZE][SIZE][4];
  public:
  // basic feature of lattice_config
  config ();

  void coldstart();
  void randomstart();
  double staplecal(const int *y, int d);
  double sweep(double beta);

  };

config::config(){
  int x[4];
  int d;
  for (x[0]=0;x[0]<SIZE;x[0]++)
    for (x[1]=0;x[1]<SIZE;x[1]++)
      for (x[2]=0;x[2]<SIZE;x[2]++)
        for (x[3]=0;x[3]<SIZE;x[3]++)
          for (d=0;d<4;d++){
	        lattice[x[0]][x[1]][x[2]][x[3]][d]=1.;}
  return;
}

// utility
// 
// move +1 step along d-th direction
void moveup(int *x, int d) {
  x[d] = x[d] + 1;
  if (x[d]>=SIZE) x[d]=x[d]-SIZE; 
  return;
}

void movedown(int *x, int d) {
  x[d] = x[d] - 1;
  if (x[d]<0) x[d]=x[d]+SIZE; 
  return;
}


void config::coldstart()
  {
  int x[4];
  int d;

  for (x[0]=0;x[0]<SIZE;x[0]++)
    for (x[1]=0;x[1]<SIZE;x[1]++)
      for (x[2]=0;x[2]<SIZE;x[2]++)
        for (x[3]=0;x[3]<SIZE;x[3]++)
          for (d=0;d<4;d++){
	    lattice[x[0]][x[1]][x[2]][x[3]][d]=1.;}
  return;
  }


void config::randomstart()
  {
  int x[4],d;
  // usage: unif(rng) to get a random number
  std::uniform_real_distribution<double> unif(0., 1.);
  std::random_device rand_dev;     // Use random_device to get a random seed.
  std::mt19937_64 rng(rand_dev()); // mt19937 is a good pseudo-random number generator.

  for (x[0]=0; x[0]<SIZE; x[0]++)
    for (x[1]=0; x[1]<SIZE; x[1]++)
      for (x[2]=0; x[2]<SIZE; x[2]++)
        for (x[3]=0; x[3]<SIZE; x[3]++)
          for (d=0; d<4; d++) 
          {
            if ( unif(rng) > 0.5 ){
              lattice[x[0]][x[1]][x[2]][x[3]][d]=1.;
            }
            else{ 
              lattice[x[0]][x[1]][x[2]][x[3]][d]=-1.;
            }
          }
  return;
  }

// staple calculation
double config::staplecal(const int *y, int d)
  {
  double staple,staplesum;    
  int x[4];
  int dperp;

  // copy the address of current site
  x[0] = y[0];
  x[1] = y[1];
  x[2] = y[2];
  x[3] = y[3];

            staplesum=0.;
            for (dperp=0;dperp<4;dperp++){
              if (dperp!=d){
                /*  move around thusly:
                    dperp        6--5
                    ^            |  |
                    |            1--4
                    |            |  |
                    -----> d     2--3  */
                /* plaquette 1234 */
                movedown(x,dperp);
                staple=lattice[x[0]][x[1]][x[2]][x[3]][dperp]
                  *lattice[x[0]][x[1]][x[2]][x[3]][d];
                moveup(x,d);
                staple=staple*lattice[x[0]][x[1]][x[2]][x[3]][dperp];  
                moveup(x,dperp);
                staplesum=staplesum+staple;
                /* plaquette 1456 */
                staple=lattice[x[0]][x[1]][x[2]][x[3]][dperp];
                moveup(x,dperp);
                movedown(x,d);
                staple=staple*lattice[x[0]][x[1]][x[2]][x[3]][d];
                movedown(x,dperp);
                staple=staple*lattice[x[0]][x[1]][x[2]][x[3]][dperp];
                staplesum=staplesum+staple;
              }
	          }

  return staplesum;
  }

  // do one Monte Carlo sweep
  double config::sweep(double beta)
  {

  int x[4],d;
  double staplesum;
  double bplus,bminus;
  double action = 0.0;

  double current_link;

  // usage: unif(rng) to get a random number
  std::uniform_real_distribution<double> unif(0., 1.);
  std::random_device rand_dev;          // Use random_device to get a random seed.
  std::mt19937_64 rng(rand_dev()); // mt19937 is a good pseudo-random number generator.

  for (x[0]=0; x[0]<SIZE; x[0]++)
    for (x[1]=0; x[1]<SIZE; x[1]++)
      for (x[2]=0; x[2]<SIZE; x[2]++)
        for (x[3]=0; x[3]<SIZE; x[3]++)
          for (d=0; d<4; d++) 
          {

            staplesum = staplecal(x, d);

            /* calculate the Boltzmann weight */
            bplus=exp(beta*staplesum);
            bminus=exp(-1.*beta*staplesum);
            bplus=bplus/(bplus+bminus);
            /* the heatbath algorithm */
            if ( unif(rng) < bplus ){
              lattice[x[0]][x[1]][x[2]][x[3]][d]=1.;
	      action +=staplesum;
            }
            else{ 
              lattice[x[0]][x[1]][x[2]][x[3]][d]=-1.;
	      action -=staplesum;
            }
          }
  action /= (SIZE*SIZE*SIZE*SIZE*4*6);
  return 1.-action;
  }


int main(){
  
  // parameters
  double beta_min = 0.;
  double beta_max = 1.;
  const int N_beta = 100;
  double action = 0.;
  double beta, dbeta;

  //static config current;
  static config current;

  dbeta= (beta_max-beta_min)/(N_beta);
  current.coldstart();
  //current.randomstart();
  /* heat it up */
  for (beta=beta_max; beta>beta_min; beta-=dbeta){
    action = current.sweep(beta);
    cout << beta << " " << action  << endl;
  }

  /* cool it down */
  for (beta=beta_min; beta<beta_max; beta+=dbeta){
    action = current.sweep(beta);
    cout << beta << " " << action  << endl;
  }

  return 0;
}