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;
}