Ising Modelのプログラム

http://qube.phys.kindai.ac.jp/mc.pdfを参考にしてメトロポリス法でイジングモデルのシミュレーションを行ないました。
機会があれば解説記事を書きたいとお思います。とりあえずはコードだけ。

#include <iostream>
#include <vector>
#include <cmath>
#include <ctime>
using namespace std;

typedef vector<int> Sv;
typedef vector<Sv> Sm;
const int L = 20;
const int M = 200;
const double J = 1.0;
const double B = 0.0;

double beta = 0.01;
Sm spin(L, Sv(L, 1));

double ran() {
  return (double)rand()/((double)RAND_MAX+1);
}

double E(Sm spin) {
  double sum = 0.0;
  for(int i = 0; i < L-1 ; i++) {
    for( int j = 0; j < L-1 ; j++) {
      sum += -J*(spin[i][j]*spin[i][j+1] + spin[i][j]*spin[i+1][j]);
    }
  }
  for(int i = 0; i < L; i++) {
    for(int j = 0; j < L ; j++) {
      sum += -B*spin[i][j];
    }
  }
  return sum;
}
 
double dE(int i, int j, int Pre) {
  return -J *( Pre - spin[i][j] )* ( (i == 0 ? 0 : spin[i-1][j])
				     + (i == L-1 ? 0 : spin[i+1][j])
				     + (j == L-1 ? 0 : spin[i][j+1])
				     + (j == 0 ? 0 : spin[i][j-1]))
    -B*(Pre - spin[i][j]);
}    
void sweep_1(int i, int j) {
  double r = ran();
  int Pre = -spin[i][j];
  double difE = dE(i, j, Pre);
  if(difE < 0) {
    spin[i][j] = Pre;
  } else if(r < exp(-beta*difE) ) {
    spin[i][j] = Pre;
  }
}
void sweep_all() {
  for(int i = 0; i < L ; i++) {
    for(int j = 0; j < L ; j++) {
      sweep_1(i, j);
    }
  }
}
void output() {
  for(beta = 0; beta < 1; beta += 0.01) {
    int sum = 0;
    for(int i = 0; i < L; i++ ) {
      for(int j = 0; j < L; j++) {
	spin[i][j] = pow(-1.0, rand());
      }
    }
    for(int i = 1; i < M ; i++) {
      for(int j = 0; j< L; j++) {
	for(int k = 0; k < L; k++) {
	  sum += spin[j][k];
	}
      }
      sweep_all();
    }
    cout << sum / M << endl;
  }
}
int main() {
  srand( (unsigned)time(NULL));
  output();
  return 0;
}