Ising model simulation offset critical temperature and interal ernergy

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.

#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> 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++) {
}
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);

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 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 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:

magnetization is below: ($$Q$$=magnetizatiom)

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:

But it shows me :

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

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