c++ - C++中的Ising模型

我正在用C++为2D Ising模型编写代码。这是代码应执行的操作:

生成随机的NxN晶格，每个位点的值均为+1或-1。

随机选择一个站点

如果翻转时的位点(+1到-1或-1到+1)是能量较低的状态，则翻转状态即。如果dE 4.对所有NxN站点执行此操作，而无需重复。这被认为是一次扫除。

做100次扫频。

我在执行步骤1，步骤2和步骤3时遇到问题，不胜感激!对于第1步，我设法创建并显示了一个晶格，但是似乎无法提取位置(x，y)上的站点的值。第2步和第3步，如何根据接受概率使用某种 bool(boolean) 表达式进行翻转？

 #include <cstdlib>
#include <ctime>
using namespace std;
#include <iostream>
int main() //random generation of spin configuration
{
int L;              //Total number of spins L = NxN
int N = 30          //A square lattice of length 30
double B=1;         //magnetic field
double M;           //Total Magnetization = Sum Si
double E;           //Total Energy
int T = 1.0;
int nsweeps = 100;      //number of sweeps
int de;             //change in energy when flipped
double Boltzmann;       //Boltzmann factor
int x,y;            //randomly chosen lattice site
int i,j,a,c;            //counters
  int ROWS = 5;
  int COLS = 5;
  int matrix[ROWS][COLS];
  srand ( static_cast<unsigned> ( time ( 0 ) ) );
  for ( int i = 0; i < ROWS; i++ )
  {
    for ( int j = 0; j < COLS; j++ )
    {
      matrix[i][j] = rand () % 2 *2-1;
    }
  }

 // showing the matrix on the screen
    for(int i=0;i<ROWS;i++)  // loop 3 times for three lines
    {
        for(int j=0;j<COLS;j++)  // loop for the three elements on the line
        {
            cout<<matrix[i][j];  // display the current element out of the array
        }
    cout<<endl;  // when the inner loop is done, go to a new line
    }
    return 0;  // return 0 to the OS.

//boundary conditions and range
if(x<0) x += N;
if(x>=L) x -= N;
if(y<0) y += N;
if(y>=L) y -= N;

//counting total energy of configuration
{  int neighbour = 0;    // nearest neighbour count

   for(int i=0; i<L; i++)
      for(int j=0; j<L; j++)
    {  if(spin(i,j)==spin(i+1, j))     // count from each spin to the right and above
              neighbour++;
           else
              neighbour--;
           if(spin(i, j)==spin(i, j+1))
              neighbour++;
           else
              neighbour--;
    }

    E = -J*neighbour - B*M;

//flipping spin
int x = int(srand48()*L);   //retrieves spin from randomly choosen site
int y = int(srand48()*L);

int delta_M = -2*spin(x, y);    //calculate change in Magnetization M
int delta_neighbour = spin(spinx-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
int delta_neighbour = -2*spin(x,y)* int delta_neighbour;

double delta_E = -J*delta_neighbour -B*delta_M;


//flip or not
if (delta_E<=0)
    {  (x, y) *= -1;     // flip spin and update values
           M += delta_M;
           E += delta_E;

        }



}

最佳答案

跟进我的评论:



为了帮助您入门:

将您的晶格存储为std::vector<int> lattice(N*N)

使用(x,y)访问元素data[x+N*y]。

例:

#include <vector>

struct IsingModel
{
    unsigned size_;

    std::vector<int> lattice_;

    // access element (x,y)
    int& at(int x, int y) {
        return lattice_[x + y*size_];
    }
    int at(int x, int y) const {
        return lattice_[x + y*size_];
    }

    // generate size x size lattice
    IsingModel(unsigned size)
    : size_(size), lattice_(size*size, +1) {
    }

    static int BoolToSpin(bool v) {
        return v ? +1 : -1;
    }

    // initialize spin randomly
    void initializeRandom() {
        for(int y=0; y<size_; y++) {
            for(int x=0; x<size_; x++) {
                at(x,y) = BoolToSpin(rand()%2);
            }
        }
    }

    static int Energy(int a, int b) {
        return (a == b) ? +1 : -1;
    }

    // compute total energy
    unsigned computeTotalEnergy() const {
        unsigned energy = 0;
        for(int y=1; y<size_-1; y++) {
            for(int x=1; x<size_-1; x++) {
                energy += Energy(at(x,y), at(x+1,y));
                energy += Energy(at(x,y), at(x,y+1));
            }
        }
        return energy ;
    }

 };

 #include <iostream>
 #include <cstdlib>
 #include <ctime>

 int main() {
     srand(static_cast<unsigned>(time(0))); // intialize random number generator
     IsingModel im(10);
     im.initializeRandom();
     unsigned energy = im.computeTotalEnergy();
     std::cout << energy << std::endl; // print energy
 }