价层电子对互斥理论的模拟验证(C++实现)

价层电子对互斥理论的基础是，分子或离子的几何构型主要决定于与中心原子相关的电子对之间的排斥作用。该电子对既可以是成键的，也可以是没有成键的（叫做孤对电子）。只有中心原子的价层电子才能够对分子的形状产生有意义的影响。分子中电子对间的排斥的三种情况为：孤对电子间的排斥（孤-孤排斥）；孤对电子和成键电子对之间的排斥（孤-成排斥）；成键电子对之间的排斥（成-成排斥）。分子会尽力避免这些排斥来保持稳定。

模型：

中心原子位置固定

中心原子与周围原子间以弹簧力相连：

F ⃗_a=k(r ⃗-r ⃗_0)

周围原子间有静电排斥作用：

F ⃗_ab=  (Z_a Z_b·((r_a ) ⃗-(r_b ) ⃗))/(4πε_0〖((r_a) ⃗-(r_b ) ⃗)〗^3 )

为了得到稳定的静态结构，引入摩擦阻力：

f ⃗=-6π"η" "v"  ⃗"r"

周围粒子初始位置在弹簧力平衡位置所在球壳上随机放置，初速为0。根据粒子的位置计算合力，由v=v0+F/m Δt更新速度，梯形法积分计算位移，循环若干步。

/*价层电子对互斥理论的模拟验证
 *By Aba
 *2015.1.17
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
const int SIZE=sizeof(double);
const int MAX_LINE=255;
const double PI=3.141592653589793;
double dt=1e-3;
double eta=0.01; //粘度系数
double K1=1.00,K2=5.00; //K1：斥力；K2：弹簧系数
double *px,*py,*pz,*vx,*vy,*vz,*fx,*fy,*fz,*M,*L,*Q;
int N,TYPES,*BTYPE,TIMES;
char stmp[MAX_LINE];
double sqr(double x) //平方
{
	return x*x;
}
double tqr(double x) //立方
{
	return x*x*x;
}
int output(const char* s)
{
	FILE *fp,*fv;
	int i;
	char stmp[255];
	sprintf(stmp,"xyz_%s.txt",s);
	fp=fopen(stmp,"w");
	sprintf(stmp,"v_%s.txt",s);
	fv=fopen(stmp,"w");
	fprintf(fp,"%d\n",N);
	fprintf(fv,"%d\n",N);
	for(i=0;i<N;i++)
	{
		fprintf(fp,"%lf %lf %lf\n",px[i],py[i],pz[i]);
		fprintf(fv,"%lf %lf %lf\n",vx[i],vy[i],vz[i]);
	}
	fclose(fp);
	fclose(fv);
	return 0;
}
int initial()
{
	FILE *source=NULL,*fp=NULL,*fv=NULL;
	int i,j,k,num;
	double a,b;
	if((source=fopen("source.txt","r")))
	{
		fgets(stmp,MAX_LINE,source);
		sscanf(stmp,"%d%d%d%lf%lf%lf%lf",&N,&TYPES,&TIMES,&K1,&K2,&eta,&dt);
		M=(double*)malloc(SIZE*N);
		L=(double*)malloc(SIZE*N);
		Q=(double*)malloc(SIZE*N);
		px=(double*)malloc(SIZE*N);
		py=(double*)malloc(SIZE*N);
		pz=(double*)malloc(SIZE*N);
		vx=(double*)malloc(SIZE*N);
		vy=(double*)malloc(SIZE*N);
		vz=(double*)malloc(SIZE*N);
		fx=(double*)malloc(SIZE*N);
		fy=(double*)malloc(SIZE*N);
		fz=(double*)malloc(SIZE*N);
		BTYPE=(int*)malloc(sizeof(int)*N);
		j=0;
		for(i=0;i<TYPES;i++)
		{
		 	fgets(stmp,MAX_LINE,source);
			sscanf(stmp,"%d%lf%lf%lf",&num,M+j,L+j,Q+j);
			BTYPE[j]=i;
			for(k=1;k<num;k++)
			{
				M[j+k]=M[j];
				L[j+k]=L[j];
				Q[j+k]=Q[j];
				BTYPE[j+k]=i;
			}
			j+=k;
		}
		fclose(source);
	}
	else exit(1);
	if((fp=fopen("xyz_continue.txt","r")))
	{
		if((fv=fopen("v_continue.txt","r")))
		{
			fscanf(fp,"%d",&i);
			fscanf(fv,"%d",&i);
			for(i=0;i<N;i++)
			{
				fscanf(fp,"%lf%lf%lf",px+i,py+i,pz+i);
				fscanf(fv,"%lf%lf%lf",vx+i,vy+i,vz+i);
			}
			fclose(fv);
			fclose(fp);
			return 1;
		}
		fclose(fp);
	}
	srand((unsigned)time(NULL));
	for(i=0;i<N;i++)
	{
		a=rand()*1.0/(RAND_MAX)*PI*2-PI;
		b=rand()*1.0/(RAND_MAX)*PI*2-PI;
		px[i]=L[i]*sin(a)*cos(b);
		py[i]=L[i]*sin(a)*sin(b);
		pz[i]=L[i]*cos(a);
		vx[i]=vy[i]=vz[i]=0;
	}
	output("start");
	return 0;
}
int force()
{
	int i,j;
	//弹簧力：F=-K2*(L-L0)
	for(i=0;i<N;i++)
	{
		double L1=sqrt(sqr(px[i])+sqr(py[i])+sqr(pz[i]));
		double F0=K2*(L[i]/L1-1);
		fx[i]=px[i]*F0;
		fy[i]=py[i]*F0;
		fz[i]=pz[i]*F0;
	}
//	getchar();
	//库仑排斥力 F=k1/R^2
	for(i=0;i<N-1;i++)
	{
		for(j=i+1;j<N;j++)
		{
			double dx=px[j]-px[i];
			double dy=py[j]-py[i];
			double dz=pz[j]-pz[i];
			double L2=sqr(dx)+sqr(dy)+sqr(dz);
			double L1=sqrt(L2);
			double k=K1*Q[i]*Q[j]/L2/L1; //力场：f=KQ1Q2/d^2,除L1是为了求分量
			dx*=k;
			dy*=k;
			dz*=k;
			fx[j]+=dx;
			fy[j]+=dy;
			fz[j]+=dz;
			fx[i]-=dx;
			fy[i]-=dy;
			fz[i]-=dz;
		}
	}
	//计算摩擦力 F=6*pi*η*v*r=eta*v
	for(i=0;i<N;i++)
	{
		fx[i]-=eta*vx[i];
		fy[i]-=eta*vy[i];
		fz[i]-=eta*vz[i];
	}
	return 0;
}
int run()
{
	int i;
	double vx_old,vy_old,vz_old;
	for(i=0;i<N;i++)
	{
		vx_old=vx[i];
		vy_old=vy[i];
		vz_old=vz[i];
		vx[i]+=fx[i]*dt/M[i];
		vy[i]+=fy[i]*dt/M[i];
		vz[i]+=fz[i]*dt/M[i];
		px[i]+=(vx[i]+vx_old)*dt/2;
		py[i]+=(vy[i]+vy_old)*dt/2;
		pz[i]+=(vz[i]+vz_old)*dt/2;
	}
	return 0;
}
double energy()
{
    double esum=0,etmp;
    int i,j;
    etmp=0;
    for(i=0;i<N;i++)
    {
        etmp += K2/2*sqr(sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])-L[i]);
    }
    esum+=etmp;
    etmp=0;
    for(i=0;i<N-1;i++) //计算排斥能
        for(j=i+1;j<N;j++)
            etmp += K1*Q[i]*Q[j]/sqrt(sqr(px[i]-px[j])+sqr(py[i]-py[j])+sqr(pz[i]-pz[j]));
    esum+=etmp;
    return esum;
}
void dataout()
{
	int i,*TPS;
	double ene;
	char TABLE[]={'O','C','N','F'},stmp2[MAX_LINE];
	FILE *out;
	TPS=(int*)malloc(TYPES*sizeof(int));
	memset(TPS,0,TYPES*sizeof(int));
	strcpy(stmp,"Al");
	for(i=0;i<N;i++)
	{
		TPS[BTYPE[i]]++;
	}
	for(i=0;i<TYPES;i++)
	{
		sprintf(stmp2,"%c%d",TABLE[i],TPS[i]);
		strcat(stmp,stmp2);
	}
	strcpy(stmp2,stmp);
	strcat(stmp2,".xyz");
	out=fopen(stmp2,"w");
	fprintf(out,"%d\n%s\nAl 0.0 0.0 0.0\n",N+1,stmp);
	for(i=0;i<N;i++)
	{
		fprintf(out,"%c %lf %lf %lf\n",TABLE[BTYPE[i]],px[i],py[i],pz[i]);
	}
	fclose(out);
	ene=energy();
	if(!(out=fopen("energy.txt","a")))out=fopen("energy.txt","w");
	fprintf(out,"%s\t%lf\n",stmp,ene);
	fclose(out);
	return;
}

int main(int argc,char **argv)
{
	int i;
	initial();
	for(i=0;i<TIMES;i++)
	{
		force();
		run();
	}
	run();
	output("end");
	dataout();
}

数据从同目录的source.txt文件输入，格式举例：

7 1 100000 0.1 1.0 0.05 0.01 //粒子总数，粒子种类数，步长，K1-斥力系数，K2-弹簧的劲度系数，eta-阻力系数，dt 

6 1.0 1.0 -1.0 //粒子A个数，质量，据中心平衡键长，电荷量
1 1.0 1.0 -1.2 //粒子B个数，质量，据中心平衡键长，电荷量

将粒子B的平衡键长改为短于粒子A的平衡键长，或将电荷量增大以模拟孤对电子。