wzj1.cpp

根据分子间相互作用——Lennard-Jones模型进行分子间相互作用模拟

//--------------------
//Author:
//Computer:Legend 800
//OS:windows 2000
//Compile:VC++ 6.0
//--------------------

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "iostream.h"
#define  Lx 6
#define  Ly 6

double dt;

void Initial(double *x,double *y, double *vx,double *vy)  //初值子函数
{   
//对位置取初值，假定16个粒子均匀分布在空腔内
	int i,j,k;
	float vmax;

	k=1;
	for (i=0;i<4;i++)
		for(j=0;j<4;j++)
		{
			x[k]=0.75+1.5*j;
			y[k]=0.75+1.5*i;
			k=k+1;	
		}
      
 //对速率取初值，按随机取数
	printf("\n请输入vmax: ");
	scanf("%f",&vmax);
    for (i=1;i<=16;i++)
    {
		vx[i]=vmax*(rand()%100)/100.;
	    vy[i]=vmax*(rand()%100)/100.;
	}
}


void Accel(double *x,double *y ,double *ax,double *ay) //加速度子函数
{   
	int i,j;
	double r,F,dx,dy;
	for (i=1;i<=16;i++)
	{
		ax[i]=0;
		ay[i]=0;
	}

	for(i=1;i<=15;i++)
	{
		for(j=i+1;j<=16;j++)
		{ 
			dx=x[i]-x[j];
			dy=y[i]-y[j];          
			r=sqrt((pow(dx,2)+pow(dy,2)));       //两粒子间距离  
			F=24*(1/(pow(r,8)-1/(pow(r,14))));   //粒子间作用加速度
        
			ax[i]=ax[i]+F*dx;
			ay[i]=ay[i]+F*dy;
			ax[j]=ax[j]-F*dx;
			ay[j]=ay[j]-F*dy;
		}
	}

}

	 
void Verlet( double *x,double *y,double *vx,double *vy,double *ax,double *ay)    //Verlet方程子函数 
{  
	int i;
    for (i=1;i<=16;i++)
	{  
		x[i]=x[i]+vx[i]*dt+0.5*ax[i]*dt*dt;          //考虑周期边界条件   
		if(x[i]>6) x[i]=x[i]-(int)(x[i]/6)*6;  
		if(x[i]<0) x[i]=x[i]+(int)(x[i]/6+1)*6; 
		
		y[i]=y[i]+vy[i]*dt+0.5*ay[i]*dt*dt; 
		if(y[i]>6) y[i]=y[i]-(int)(y[i]/6)*6;
		if(y[i]<0) y[i]=y[i]+(int)(y[i]/6+1)*6; 
		   
        vx[i]=vx[i]+0.5*ax[i]*dt;
		vy[i]=vy[i]+0.5*ay[i]*dt;

		Accel(x,y,ax,ay);

		vx[i]=vx[i]+0.5*ax[i]*dt;
		vy[i]=vy[i]+0.5*ay[i]*dt;
	}
}


void main()
{ 
	int i,j;
	double	x[17],y[17],vx[17],vy[17],ax[17],ay[17];
	double T; 
	int    isum[100],n;
	double vv;
	FILE *fp;
 
	cout<<"请输入步长input dt=";
	cin>>dt;
	printf("\ndt=%f\n",dt);

    printf("\n请输入演算步骤:");
	scanf("%d",&n);
   

	Initial(x,y,vx,vy);
	Accel(x,y,ax,ay);

	for(i=1;i<=n;i++)
	{  

		Verlet(x,y,vx,vy,ax,ay);

		T=0;
		for(j=1;j<17;j++)
		{	
			T=T+0.5*(vx[j]*vx[j]+vy[j]*vy[j])*119.8; 
		}
			printf("平衡温度T=%f\n\n",T/16);
	}


	   for(i=0;i<100;i++)
	   {
		   isum[i]=0;
	   }

	 if((fp=fopen("e:\\作业源程序\\01.txt","w"))==NULL)
	{
		printf("can not open the file!");
		exit(0);
	}

		for(i=1;i<=10000;i++)
		{  
			Verlet(x,y,vx,vy,ax,ay);
	     		
			for(j=1;j<17;j++)
			{	
				vv=sqrt(vx[j]*vx[j]+vy[j]*vy[j]); 
				isum[(int)(vv/0.1)]=isum[(int)(vv/0.1)]+1;    
			}
		}

			for(i=0;i<100;i++)
			  fprintf(fp,"%f       %d\n",0.1*i,isum[i]);     //将速率统计数写入文本
		
        fclose(fp);
}