0
$\begingroup$

I'm writing a code for the Ising model using WHAM (the weighted histogram analysis method),But it seems to produce critical temperature and internal energy wrong.

(newest rewritten code is below)

#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <fstream>
#include <iostream>
#include <map>
#include <vector>

using namespace std;

ofstream out("data.dat");
const long  L = 8;       //模拟尺寸
const long  MCS = 1.1e6; //总的模拟步数
const long  MS = 1e5;     //除去预热统计的步数
const long double J = -1;    // J=1AFM,J=-1FM

//随机初始化自旋分布
long double InitialSpin(long  s[][L]) {
  for (long  i = 0; i < L; i++) {
    for (long  j = 0; j < L; j++) {
      // s[i][j] = 1;
      s[i][j] = 2 * (rand() % 2) - 1;
    }
  }
  return 0;
}

// 计算体系的总能量
long double SumEnergy(long  Spin[][L]) {
  long double e = 0;
  for (long  i = 0; i < L; i++) {
    for (long  j = 0; j < L; j++) {
      e += J * Spin[i][j] *
           (Spin[(i + 1) % L][j] + Spin[(i - 1 + L) % L][j] +
            Spin[i][(j + 1) % L] + Spin[i][(j - 1 + L) % L]) /
           2; //+L为了防止i-1负数
    }
  }
  return e;
}
//计算格点的局域能量
long  DeltaPartEnergy(long  Spin[][L], long  i, long  j) {
  long  e1 = 0, e2 = 0, deltae;

  e1 = J * Spin[i][j] *
       (Spin[(i + 1) % L][j] + Spin[(i - 1 + L) % L][j] + Spin[i][(j + 1) % L] +
        Spin[i][(j - 1 + L) % L]); //+L为了防止i-1负数
  e2 = -J * Spin[i][j] *
       (Spin[(i + 1) % L][j] + Spin[(i - 1 + L) % L][j] + Spin[i][(j + 1) % L] +
        Spin[i][(j - 1 + L) % L]);
  deltae = e2 - e1;

  return deltae;
}
//输出体系的总磁矩
long  Magnetization(long  s[][L]) {
  long  m = 0;
  for (long  i = 0; i < L; i++) {
    for (long  j = 0; j < L; j++) {
      m += s[i][j];
    }
  }
  return m;
}

//生成选中的格点坐标(经过确认,ubu上可以2^32,vscode坑爹的只有16位,远远大于所需位数)
long  Myrand(long  L) {
  long  x;
  x = rand() % L;
  return x;
}
// monte交换
void Montechange(long  Spin[][L], long double T, long  *result) {
  long  x, y, DeltaE, trans = 0;
  long double p, r;
  x = Myrand(L);
  y = Myrand(L);
  DeltaE = DeltaPartEnergy(Spin, x, y);
  if (DeltaE < 0)
    trans = 1;
  else {
    p = 1.0 * rand() / (long double)RAND_MAX;
    r = exp(1.0 * (-DeltaE) / T);
    if (p < r)
      trans = 1;
  }
  if (trans == 1) {
    result[0] = x;
    result[1] = y;
    result[2] = trans;
  }
}

void MonteCarlo(long  Spin[][L], long double T, map<long , long > &ener,
                map<long , long > &mag, map<vector<long >, long > &mage) {

  long  m, trans[3] = {0};
  long double SE;
  std::vector<long > mage1;
  for (long  mcs = 0; mcs < MCS; mcs++) {

    Montechange(Spin, T, trans);
    if (trans[2] == 1) {
      Spin[trans[0]][trans[1]] = -Spin[trans[0]][trans[1]];
    }

    if (mcs == MS - 1) {
      SE = SumEnergy(Spin);
      m = Magnetization(Spin);
    }
    if (mcs >= MS) {
      if (trans[2] == 1) {
        SE = SE - DeltaPartEnergy(Spin, trans[0], trans[1]);
        m = m + 2 * Spin[trans[0]][trans[1]]; //输出总磁矩
      }
      mage1.clear();
      mage1.push_back(SE);
      mage1.push_back(abs(m));
      if (mag.find(abs(m)) == mag.end())
        mag[abs(m)] = 1;
      else
        mag[abs(m)]++;
      // SE=SumEnergy(Spin);
      if (ener.find(SE) == ener.end())
        ener[SE] = 1;
      else
        ener[SE]++;

      if (mage.find(mage1) == mage.end())
        mage[mage1] = 1;
      else
        mage[mage1]++;
    }

    // out << mcs << '\t' << m <<'\t' << SE << endl;
  }
  // for(map<long , long >::iterator it = ener.begin(); it != ener.end(); it++) {
  // out << it->first << '\t' << it->second << endl;
  //}
}
/*
void single_histogram(map<long , long > ener, map<vector<long >, long > mage,
                      vector<long double> T, vector<long double> Tadd, vector<long double> &U,
                      vector<long double> &magT) {
  long double up, down, a;

  for (long  i = 0; i < Tadd.size(); i++) {
    up = 0;
    down = 0;
    for (map<long , long >::iterator it = ener.begin(); it != ener.end(); it++) {
      a = 1.0 * it->second *
          exp(-1.0 * ((1.0 / Tadd[i]) - (1.0 / T[0])) * it->first);
      down += a;
      up += a * it->first;
    }
    U.push_back(up / down);
    up = 0;
    down = 0;
    for (map<vector<long >, long >::iterator it = mage.begin(); it != mage.end();
         it++) {
      a=1.0 * exp(-1.0 * ((1.0 / Tadd[i]) - (1.0 / T[0])) * it->first[0]);
      down+=a*it->second;
      up+=a*it->first[1]*it->second;
    }
    magT.push_back(up/down);
  }
}*/
void normalize(vector<long double> &Z) {
  auto maxPosition = max_element(Z.begin(), Z.end());
  auto minPosition = min_element(Z.begin(), Z.end());
  long i;
  long double max, min, A;
  max = *maxPosition;
  min = *minPosition;
  A = 1.0 / sqrt(max * min);
  for (i = 0; i < Z.size(); i++) {
    Z[i] = Z[i] * A;
  }
}
void iterZ(vector<long double> T, vector<long double> &Z, map<long , long> ener) {
  long  i, j;
  long double delta, epsilon = 1e-10;

    for (i = 0; i < T.size(); i++) {
      Z.push_back(1.0);
    }
    std::vector<long double> Y(Z.size());
    do {
        delta=0.0;
    for (j = 0; j < T.size(); j++) {
      long double up = 0, down = 0,sum=0;
      for (map<long , long >::iterator it = ener.begin(); it != ener.end(); it++) {
        up = 1.0 * it->second;
        for (i = 0; i < T.size(); i++) {
          down +=
              1.0 * (MCS - MS) / Z[i] * exp((1.0 / T[j] - 1.0 / T[i]) * it->first);
        }
        sum+=up/down;
      }
      Y[j] = sum;
    }
    normalize(Y);
    for (j = 0; j < T.size(); j++) {
      delta += ((Y[j] - Z[j]) / Y[j]) * ((Y[j] - Z[j]) / Y[j]);
      Z[j] = Y[j];
    }
    cout<<Z[0]<<endl;
    cout<<Z[1]<<endl;
    cout<<delta<<endl;
  } while (delta > epsilon * epsilon);
}
void multi(vector<long double> T, vector<long double> Tadd, vector<long double> Z,
           vector<long double> &Y, map<long , long > ener) {
  long  i, j;
  for (i = 0; i < Tadd.size(); i++) {
    long double up = 0.0, down = 0.0,sum=0.0;
    for (map<long , long >::iterator it = ener.begin(); it != ener.end(); it++) {
      up = 1.0 * it->second;
      for (j = 0; j < T.size(); j++) {
        down +=
            1.0 * (MCS - MS) / Z[j] * exp((1.0 / Tadd[i] - 1.0 / T[j]) * it->first);
      }
      sum+=up/down;
    }
    Y.push_back(sum);
  }
}
int  main() {
  srand((unsigned)time(NULL)); //初始化时间
  long  s[L][L] = {{0}}, i, j, k;
  //long double enerp, magp;
  std::map<long , long > ener;
  std::map<long , long > mag;
  std::vector<long double> T;
  std::vector<long double> Tadd;
  //std::vector<long double> Tm;
  //std::vector<long double> Te;
  std::map<vector<long >, long > mage;
  std::vector<long double> U;
  std::vector<long double> Z;
  std::vector<long double> Y;
  std::vector<long double> magT;
  for (i = 0; i < 2; i++) {
    Tadd.push_back(2.1+0.1*i);
  }
  for(i=0;i<2;i++){
  T.push_back(2.1+0.1*i);
  }
  for (k = 0; k < T.size(); k++) {
    for (i = 0; i < L; i++) {
      for (j = 0; j < L; j++) {
        // s[i][j] = 1;
        s[i][j] = 2 * (rand() % 2) - 1;
      }
    }
    MonteCarlo(s, T[k], ener, mag, mage);
    //magp = 0.0;
    /*for (map<long , long >::iterator it = mag.begin(); it != mag.end(); it++) {
      magp += (it->first) * (it->second);
    }
    Tm.push_back(magp / (MCS - MS) / L / L);
    enerp = 0.0;
    for (map<long , long >::iterator it = ener.begin(); it != ener.end(); it++) {
      enerp += (it->first) * (it->second);
    }
    Te.push_back(enerp / (MCS - MS) / L / L);*/
  }
  iterZ(T, Z, ener);
  multi(T,Tadd,Z,Y,ener);

  // single_histogram(ener,mage,T,Tadd,U,magT);

  ofstream out3("magA.dat");
  for (i = 0; i < T.size(); i++) {
    out3 << T[i] << '\t' << Z[i]  << endl;
  }
  for (i = 0; i < Tadd.size(); i++) {
    out3 << Tadd[i] << '\t' << Y[i]<< endl;
  }
  out3.close();
  /*
  ofstream out1("mag.dat");
  for (map<long , long >::iterator it = mag.begin(); it != mag.end(); it++) {
    out1 << it->first << '\t' << it->second << endl;
  }
  out1.close();

  ofstream out2("ener.dat");
  for (map<long , long >::iterator it = ener.begin(); it != ener.end(); it++) {
    out2 << it->first << '\t' << it->second << endl;
  }
  out2.close();
  */
  return 0;
}

This is the result of my code:

magnetization-kT magnetization-kT interal ernergy-kT interal ernergy-kT And it's my single-histogram result.

The blue line is simulation one by one,the orange line is single-histogram.

magnetization-kT enter image description here interal ernergy-kT ![interal ernergy-kT

As you see, they are different. I don't know why.

Because my single-histogram code is same to this code until MonteCarlo(){}, so I think there is error in normalize(), iteration(), or main()

And the formula is: enter image description here

magnetization is below: ($Q$=magnetizatiom)

enter image description here

ok, I test this code again. The error is that after multi-histogram templating, the partition function of 1.2 is not a simulation partition simulation of 1.2. What should I do?

After I tried again and again, I find the error is in iterZ() and multi()

According to Newman's Monte Carlo Methods in Statistical Physics, it should converge naturally:

enter image description here

But it shows me :

enter image description here

The first formula is different from the second. It's awful. I don't know why.

$\endgroup$
  • 2
    $\begingroup$ It would help if you remove unnecessary commented out portions of your code and translate the comments to English. $\endgroup$ – Anton Menshov Dec 19 '19 at 16:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.