0
$\begingroup$

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: Energy for specific position

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: Energy for frustrated case

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: Frustrated model 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 Specific heat for p=0.1

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 Specific heat for various p

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[0];
  }
  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);
          //adjust observables
          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;
}
Improve this question
$\endgroup$
1
  • $\begingroup$ 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). $\endgroup$
    – Kyle Kanos
    Aug 10, 2022 at 11:35

1 Answer 1

Reset to default
1
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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). $\endgroup$
    – Kyle Kanos
    Aug 10, 2022 at 11:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.