一道分子模拟练习题

最近帮人做的分子模拟练习作业，一个简单的C++程序，题目如下：

1、构筑一个100´100的二维方格，随机挑选一半格点置为0，其余格点置为1。如果相邻的两个格点的赋值均为1，则有相互作用能e/kT=0.3，否则相互作用能为0，如果一个格点处于盒子边缘，则令它与盒子壁面的相互作用能为e/kT=0.15 。系统总的能量为所有相邻格点相互作用能之和。如此产生10000个样本，统计系统总能量的平均值和方差。
2、在上题中，产生初始样本并统计出系统总能量后，我们也可以用下面的方法产生新样本：在盒子中随机挑选两个格点，如果它们的赋值不同，则交换它们的赋值，即将赋值为1的格点置为0，将赋值为0的格点置为1。如果两个格点的赋值相等，则不作任何动作。每做10000次取一个样本，计算系统的总能量。共取10000个样本，统计系统总能量的平均值。比较上述两种方法获得的结果差别和所花时间的差别。

自己写程序的能力还是很弱，原本想将产生系统分布与计算能量写成两个函数的，搞了半天还是写在一起了。计算的速度也不快，暂且先记录下来，看日后能不能再优化一下：

第一题

// 分子模拟课程作业ppt4，第一题
#include<iostream>
#include<cstdio>
#include<cmath>
#include<string>
#include<cstdlib>
#include<ctime>

using namespace std;

const int MAXN = 100;
const double kT1 = 0.3;
const double kT2 = 0.15;


int main(){

int step=0;
double average=0,variance=0;
double sum[10000];
for (int i=0;i<10000;i++)
    sum[i]=0;
srand((unsigned)time(NULL));

while (step<10000){
    int box[MAXN][MAXN];
    int x,y;
    int n=0;
    double energy=0.0;
    
    for (int i=0;i<MAXN;i++){
        for (int j=0;j<MAXN;j++)
            box[i][j]=0;
    }

    while (n<((MAXN*MAXN)/2)){
        x=rand()%MAXN;
        y=rand()%MAXN;
        if (box[x][y] == 0){
            box[x][y] =1;
            n++;
            //cout << box[x][y] << "\t" << n << endl;
        }
    }
    int num = 0;
    for (int i=0;i<MAXN;i++){
        for (int j=0;j<MAXN;j++){
            cout << box[i][j] << " ";
            num++;
            if (num%MAXN==0)
                cout << endl;
        }
    }
    cout << endl << "Martix complete!" << endl;
    
    for (int i=0;i<MAXN;i++){
        for (int j=0;j<MAXN;j++){
            if (i==0 || j==0 || i==(MAXN-1) || j==(MAXN-1)){
                if (box[i][j]==1)
                    energy += kT2;
            } 
            if (i!=(MAXN-1) && j!=(MAXN-1)){
                if (box[i][j]==1 && box[i+1][j]==1){
                    energy += kT1;
                }
                if (box[i][j]==1 && box[i][j+1]==1){
                    energy += kT1;
                }
            }
            if (i==(MAXN-1)){
                if (box[i][j]==1 && box[i][j+1]==1){
                    energy += kT1;
                }
            }
            if (j==(MAXN-1)){
                if (box[i][j]==1 && box[i+1][j]==1){
                    energy += kT1;
                }
            }
        }
    }
    sum[step] = energy;
    energy = 0;
    cout << "Step: " << step+1 << "\t" << "Energy: " << sum[step] << endl << endl;
    step++;
}

for (int i=0;i<step;i++){
    average += sum[i];
}
average = average/step;
for (int i=0; i<step;i++){
    variance += ((sum[i]-average)*(sum[i]-average));
}
variance = sqrt(variance/step);
cout << "The average of the energy is: " << average << endl;
cout << "The variance of the energy is: " << variance << endl;
return 0;
}

第二题

// 分子模拟课程作业ppt4，第二题
#include<iostream>
#include<cstdio>
#include<cmath>
#include<string>
#include<cstdlib>
#include<ctime>

using namespace std;

const int MAXN = 100;
const double kT1 = 0.3;
const double kT2 = 0.15;


int main(){

int step=0;
double average=0,variance=0;
double sum[10000];
for (int i=0;i<10000;i++)
    sum[i]=0;
srand((unsigned)time(NULL));
    
while (step<10000){
    int box[MAXN][MAXN];
    int x,y;
    int n=0,m=0;
    double energy=0.0;
    
    for (int i=0;i<MAXN;i++){
        for (int j=0;j<MAXN;j++)
            box[i][j]=0;
    }
    while (n<((MAXN*MAXN)/2)){
        x=rand()%MAXN;
        y=rand()%MAXN;
        if (box[x][y] == 0){
            box[x][y] =1;
            n++;
            //cout << box[x][y] << "\t" << n << endl;
        }
    }
    while (m<10000){
        int x1,x2,y1,y2;
        x1=rand()%MAXN;
        x2=rand()%MAXN;
        y1=rand()%MAXN;
        y2=rand()%MAXN;
        int temp;
        if (x1!=x2 && y1!=y2){
            if (box[x1][y1] != box[x2][y2]){
                temp = box[x1][y1];
                box[x1][y1] = box[x2][y2];
                box[x2][y2] = temp;
            }
        }
        m++;
    }
    int num = 0;
    for (int i=0;i<MAXN;i++){
        for (int j=0;j<MAXN;j++){
            cout << box[i][j] << " ";
            num++;
            if (num%MAXN==0)
                cout << endl;
        }
    }
    cout << endl << "Martix complete!" << endl;

    for (int i=0;i<MAXN;i++){
        for (int j=0;j<MAXN;j++){
            if (i==0 || j==0 || i==(MAXN-1) || j==(MAXN-1)){
                if (box[i][j]==1)
                    energy += kT2;
            } 
            if (i!=(MAXN-1) && j!=(MAXN-1)){
                if (box[i][j]==1 && box[i+1][j]==1){
                    energy += kT1;
                }
                if (box[i][j]==1 && box[i][j+1]==1){
                    energy += kT1;
                }
            }
            if (i==(MAXN-1)){
                if (box[i][j]==1 && box[i][j+1]==1){
                    energy += kT1;
                }
            }
            if (j==(MAXN-1)){
                if (box[i][j]==1 && box[i+1][j]==1){
                    energy += kT1;
                }
            }
        }
    }
    sum[step] = energy;
    //cout << energy << endl;
    energy = 0;
    cout << "Step: " << step+1 << "\t" << "Energy: " << sum[step] << endl << endl;
    step++;
}

for (int i=0;i<step;i++){
    average += sum[i];
    //cout << average << " " << sum[i] << endl;
}

average = average/step;
for (int i=0; i<step;i++){
    variance += ((sum[i]-average)*(sum[i]-average));
}
variance = sqrt(variance/step);
cout << "The average of the energy is: " << average << endl;
cout << "The variance of the energy is: " << variance << endl;
return 0;
}