# Problem in code for Monte Carlo method for frustrated Ising model

I'm studying Monte Carlo method (it's related to my school project). The code I'm using I got from a paper by Jacques Kotze (https://arxiv.org/abs/0803.0217). In his paper, he uses the Monte Carlo method for a non-frustrated Ising model. I've tried using his code (unchanged), and I can reproduce all the results of the paper (specific heat, magnetization, susceptibility, Binder cumulant). He has an error in the definition of the cumulant, that I fixed.

Now, when I try to use the same code, but for the frustrated Ising model, I am getting bad results. I'll try to explain as best as I can what I'm doing. The frustrated Ising model is a model with a Hamiltonian $$H=J_1\sum_{i,j}S_{i}S_{j}+J_2\sum_{i,j}S_{i}S_{j}$$ (sometimes we use $$p=\frac{J_2}{J_1}$$). Now if I want to add the second part of the Hamiltonian in my code (that is, the code provided by Kotze), I tried changing the part of the code where energy is being calculated: I changed the energy part from:

//energy for specific position
e =-1* lat[pos.x][pos.y] *
(lat[left][pos.y] + lat[right][pos.y] + lat[pos.x][up] + lat[pos.x][down]);


This can be seen in picture: I tried changing this energy to:

e =1 * lat[pos.x][pos.y] *
(lat[left][pos.y] + lat[right][pos.y] + lat[pos.x][up] + lat[pos.x][down])+
0.1*lat[pos.x][pos.y]*(lat[left][down]+lat[left][up]+lat[right][down]+lat[right][up]);


This can be seen in picture: 1 and 0.1 represent $$J_1$$ and p respectively. I understood this part of the code to be: calculate energy using nearest neighbors and next nearest neighbors. lat[pos.x][pos.y] represents a point on a square lattice and up, down, left and right are defined with periodic boundary condition (in the code it's defined above the part with energy e=1*.....). I imagend it to be like: I think that I'm making a mistake here, that is, I did not write next nearest neigbors properly. If I use this code now (with added next nearest neighbors) and for lattice size LxL where L=8 I get for specific heat If I compare this specific heat with specific heat in paper by A K Murtazaev, M K Ramazanov and F A Kassan-Ogly (https://iopscience.iop.org/article/10.1088/1742-6596/510/1/012026) where they got My specific heat has the right shape, has peak at approximately the same critical temperature but my values for specific heat are huge. If I were to use bigger lattice size the values for specific heat will get even bigger. Clearly, I'm making a mistake somewhere. Another mistake that I'm getting is that my values for specific heat for different temperatures are literaly the same for p=0.2,p=0.3,p=0.4. I get the same graph, peak at the same place, it's like I did not change p's.

I know that there are better Monte Carlo methods for frustrated models, but I need to do it with this MC method (that is, with this code). Can someone spot where I'm making a mistake. It's been a long time since I did anything in C++, so maybe I'm missing something obvious. I hope I explained properly my problem. Thank you!!

I don't know if I need to, but I will add code by Kotze here:

#include <iostream>

#include <math.h>

#include <stdlib.h>

#include <fstream>
//random number generator from "Numerical Recipes in C" (Ran1.c)
#include <random>
//file to output data into

using namespace std;

// ran1 random number generator

#define IA 16807
#define IM 2147483647
#define AM (1.0 / IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1 + (IM - 1) / NTAB)
#define EPS 1.2e-7
#define RNMX (1.0 - EPS)

float ran1(long * idum) {
int j;
long k;
static long iy = 0;
static long iv[NTAB];
float temp;

if ( * idum <= 0 || !iy) {
if (-( * idum) < 1) * idum = 1;
else *idum = -( * idum);
for (j = NTAB + 7; j >= 0; j--) {
k = ( * idum) / IQ;
* idum = IA * ( * idum - k * IQ) - IR * k;
if ( * idum < 0) * idum += IM;
if (j < NTAB) iv[j] = * idum;
}
iy = iv;
}
k = ( * idum) / IQ;
* idum = IA * ( * idum - k * IQ) - IR * k;
if ( * idum < 0) * idum += IM;
j = iy / NDIV;
iy = iv[j];
iv[j] = * idum;
if ((temp = AM * iy) > RNMX) return RNMX;
else return temp;
}

//structure for a 2d lattice with coordinates x and y

struct lat_type {
int x;
int y;
};
const int size = 8; //lattice size
const int lsize = size - 1; //array size for lattice
const int n = size * size; //number of spin points on lattice
float T = 5.0; //starting point for temperature
const float minT = 0.5; //minimum temperature
float change = 0.1; //size of steps for temperature loop
int lat[size + 1][size + 1]; //2d lattice for spins
long unsigned int mcs = 10000; //number of Monte Carlo steps
int transient = 1000; //number of transient steps
double norm = (1.0 / float(mcs * n)); //normalization for averaging
long int seed = 436675; //seed for random number generator
//function for random initialization of lattice
initialize(int lat[size + 1][size + 1]) {
for (int y = size; y >= 1; y--) {
for (int x = 1; x <= size; x++) {
if (ran1( & seed) >= 0.5)
lat[x][y] = 1;
else
lat[x][y] = -1;
}
}

for (int i = 0; i <= size; i++) {
for (int j = 0; j <= size; j++) {
cout<<lat[i][j]<< " ";
}
cout<<endl;
}
}
//output of lattice configuration to the screen
output(int lat[size + 1][size + 1]) {
for (int y = size; y >= 1; y--) {
for (int x = 1; x <= size; x++) {
if (lat[x][y] < 0)
cout << " − ";
else
cout << " + ";
}
cout << endl;
}
}
//function for choosing random position on lattice
choose_random_pos_lat(lat_type & pos) {
pos.x = (int) ceil(ran1( & seed) * (size));
pos.y = (int) ceil(ran1( & seed) * (size));
if (pos.x > size || pos.y > size) {
cout << "error in array size allocation for random position on lattice!";
exit;
}
}
//function for calculating energy at a particular position on lattice
int energy_pos(lat_type & pos) {
//periodic boundary conditions
int up, down, left, right, e;
if (pos.y == size)
up = 1;
else
up = pos.y + 1;
if (pos.y == 1)
down = size;
else
down = pos.y - 1;

if (pos.x == 1)
left = size;
else
left = pos.x - 1;
if (pos.x == size)
right = 1;
else
right = pos.x + 1;
//energy for specific position
e =1 * lat[pos.x][pos.y] *
(lat[left][pos.y] + lat[right][pos.y] + lat[pos.x][up] + lat[pos.x][down])+
0.1*lat[pos.x][pos.y]*(lat[left][down]+lat[left][up]+lat[right][down]+lat[right][up]);
return e;
}
//function for testing the validity of flipping a spin at a selected position
bool test_flip(lat_type pos, int & de) {
de = -2 * energy_pos(pos); //change in energy for specific spin
if (de < 0)
return true; //flip due to lower energy
else if (ran1( & seed) < exp(-de / T))
return true; //flip due to heat bath
else
return false; //no flip
}
//flip spin at given position
flip(lat_type pos) {
lat[pos.x][pos.y] = -lat[pos.x][pos.y];
}
//function for disregarding transient results
transient_results() {
lat_type pos;
int de = 0;
for (int a = 1; a <= transient; a++) {
for (int b = 1; b <= n; b++) {
choose_random_pos_lat(pos);
if (test_flip(pos, de)) {
flip(pos);
}
}
}
}
//function for calculating total magnetization of lattice
int total_magnetization() {
int m = 0;
for (int y = size; y >= 1; y--) {
for (int x = 1; x <= size; x++) {
m += lat[x][y];
}
}
return m;
}

//function for calculating total energy of lattice
int total_energy() {
lat_type pos;
int e = 0;
for (int y = size; y >= 1; y--) {
pos.y = y;
for (int x = 1; x <= size; x++) {
pos.x = x;
e += energy_pos(pos);
}
}
return e;
}
//main program
int main() {

ofstream DATA1("Mavg.dat",ios::out);
ofstream DATA2("Mabsavg.dat",ios::out);
ofstream DATA3("Msqavg.dat",ios::out);
ofstream DATA4("X.dat",ios::out);
ofstream DATA5("X'.dat",ios::out);
ofstream DATA6("Eavg.dat",ios::out);
ofstream DATA7("Esqavg.dat",ios::out);
ofstream DATA8("C.dat",ios::out);
ofstream DATA9("U.dat",ios::out);
ofstream DATA10("Mqavg.dat",ios::out);

//declaring variables to be used in calculating the observables
double E = 0, Esq = 0, Esq_avg = 0, E_avg = 0, etot = 0, etotsq = 0;
double M = 0, Msq = 0, Msq_avg = 0, M_avg = 0, mtot = 0, mtotsq = 0;
double Mabs = 0, Mabs_avg = 0, Mq_avg = 0, mabstot = 0, mqtot = 0;
int de = 0;
lat_type pos;
//initialize lattice to random configuration
initialize(lat);
//Temperature loop
for (; T >= minT; T = T - change) {
//transient function
transient_results();
//observables adopt equilibrated lattice configurations values
M = total_magnetization();
Mabs = abs(total_magnetization());
E = total_energy();
//initialize summation variables at each temperature step
etot = 0;
etotsq = 0;
mtot = 0;
mtotsq = 0;
mabstot = 0;
mqtot = 0;
//Monte Carlo loop
for (int a = 1; a <= mcs; a++) {
//Metropolis loop
for (int b = 1; b <= n; b++) {
choose_random_pos_lat(pos);
if (test_flip(pos, de)) {
flip(pos);
E += 2 * de;
M += 2 * lat[pos.x][pos.y];
Mabs += abs(lat[pos.x][pos.y]);
}
}
//keep summation of observables
etot += E / 2.0; //so as not to count the energy for each spin twice
etotsq += E / 2.0 * E / 2.0;
mtot += M;
mtotsq += M * M;
mqtot += M * M * M * M;
mabstot += (sqrt(M * M));
}
//average observables
E_avg = etot * norm;
Esq_avg = etotsq * norm;
M_avg = mtot * norm;
Msq_avg = mtotsq * norm;
Mabs_avg = mabstot * norm;
Mq_avg = mqtot * norm;
//output data to file
DATA1<< T << "\t " << M_avg<<endl;

DATA2 << T << "\t"<< Mabs_avg<<endl;

DATA3<<T<<"\t" << Msq_avg << endl; //<M>;<|M|>;<M^2> per spin

DATA4<<T<<"\t" << (Msq_avg - (M_avg * M_avg * n)) / (T) <<endl; //susceptibility per spin (X)

DATA5<<T<<"\t" << (Msq_avg - (Mabs_avg * Mabs_avg * n)) / (T) <<endl; //susceptibility per spin (X’)

DATA6<<T<<"\t" << E_avg <<endl; //<E>;<E^2> per spin

DATA7<<T<<"\t" << Esq_avg <<endl; //<E>;<E^2> per spin

DATA8<<T<<"\t" << (Esq_avg - (E_avg * E_avg * n)) / (T * T) <<endl; //heat capacity (C) per spin

DATA9<<T<<"\t" << 1 - ((Mq_avg/(n*n*n)) / (3 * (Msq_avg * Msq_avg/(n*n)))) << endl;

DATA10<<T<<"\t" << Mq_avg << endl;

}
return 0;
}

• I would be hesitant in calling that "C++ code" as it's really C code that uses some C++ libraries (and mostly for the print statements). Aug 10, 2022 at 11:35

The Hamiltonian seems reasonable to me, and I did not find significant errors in the variables defined in the code you post (there may still be, but so far I did not find any).

However, there are other quite common mistakes in your Monte-Carlo simulation:

1. Thermalization Steps

The starting configuration of your simulation is randomly generated. Therefore, before you start to samples configurations, you need to perform thermalization steps to ensure the configurations you sample are thermalized at $$T$$.

1. Autocorrelation Lengths

In Monte-Carlo simulation, you have to make sure the sampled configurations $$\sigma_i$$'s are uncorrelated with respect to the observable you want to measure. Take the energy as example. We want $${1 \over m-1}\sum_{i=1}^{m-1}{E(\sigma_i)E(\sigma_{i+1})}-\bigg({1 \over m}\sum_{i=1}^{m}{E(\sigma_i)}\bigg)^2 \sim 0$$ To do this, you have to first measure the autocorrelation length of the observable. Let us assume it is $$l$$, and when sampling configurations, we can only take $$\sigma_i$$ every $$al$$ updates where $$a$$ is sufficiently larger than $$1$$ (in practice some may just average the Markov chain of the length $$mal$$). This concept is critical particularly when the sampling is near the phase transition. And the local update usually has much larger autocorrelation lengths than the global updates although it is much harder to construct some efficient global updates method.

I am not sure whether these are the problems which lead to strangely large heat capacity. My suggestion is you can print everything out and look for anomalies which can lead to this phenomenon.

I have also written the code myself according to the principles I have mentioned about Monte-Carlo simulation, and the results are good to me. These are the graphs for $$L=8$$ and $$(J_1,J_2)=(1,0.4)$$ and $$(J_1,J_2)=(1,0.1)$$. The plots agree with the literature you gave in terms of the location of the peaks and the confined values of heat capacities.

• I did not find significant errors in the variables defined in the code you post This might be the case, but it's really bad code b/c it doesn't follow any reasonably good practices (e.g., #define instead of constexpr, lack of function return types like void, etc). Aug 10, 2022 at 11:33