一道分子模拟练习题

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

1、构筑一个100x100的二维方格，随机挑选一半格点置为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;
}

Email