分子动力学的模拟程序设计.docx
《分子动力学的模拟程序设计.docx》由会员分享，可在线阅读，更多相关《分子动力学的模拟程序设计.docx（13页珍藏版）》请在冰豆网上搜索。
分子动力学的模拟程序设计
分子动力学的模拟程序设计
一、课题名称：
分子动力学的模拟程序设计
二、班级和姓名：
***
三、主要内容：
1.研究的内容和算法
分子动力学模拟方法–确定性的演化过程。
按该体系内部的内禀动力学规律计算、确定位形的转变。
出发点：
物理系统确定的微观描述。
具体的说就是描述系统的哈密顿、拉格朗日或者牛顿运动方程，用这些方程去驱动粒子的位置、速度和取向随时间的演化。
模拟过程：
（1）通过对物理体系的微观数学描述建立一组方程组，直接求解每个分子的运动方程。
得到每个时刻、每个分子的坐标和动量。
（2）利用统计方法计算系统的静态和动态过程。
计算机模拟的问题
1观测时间是有限的
对于某些问题，有限观测时间可以看成是无限长的。
如对一个分子系统的计算，它的观测时间远大于分子时间尺度。
2观测体系是有限的。
引入周期性、全反射、漫反射等边界条件。
分子动力学模拟的计算
对体系的分子动力学方程组进行求解时，需要将运动的连续性方程离散化变成有限差分方程，常用的方法有：
欧拉法，龙格-库塔法、辛普生法等。
误差来源：
（1）动力学模型的近似程度；
（2）数值求解法的近似阶数；
（3）数值计算的舍入误差。
分子动力学的适用范围
原则上，分子动力学方法所适用的微观物理体系并无任何限制。
这个方法适用的体系既可以是少体系统，也可以是多体系统；既可以是点粒子体系，也可以是具有内部结构的体系；处理的微观客体既可以是分子，也可以是其它的微观粒子。
微正则系综（NVE）的分子动力学模拟
步骤如下：
（1）给定初始空间位置
（2）计算第n步时粒子所受的力
（3）计算第n+1步粒子的空间位置
（4）计算第n步的速度，根据速度得出系统的温度。
（5）返回第
（2）步，开始下一步计算。
2.源程序
1）简立方结构
/*
#include
#include
#include
#include
#include
*/
#include
#include
#include
#include
#include
#include
#include
usingnamespacestd;
constintN=64;//numberofparticles
doubler[N][3];//positions
doublev[N][3];//velocities
doublea[N][3];//accelerations
doubleL=10;//linearsizeofcubicalvolume，简单立方
doublevMax=0.1;//maximuminitialvelocitycomponent
voidinitialize（）
{
//initializepositions
intn=int（ceil（pow（N,1.0/3）））;//numberofatomsineachdirection
doublea=L/n;//latticespacing,晶格常数
intp=0;//particlesplacedsofar
for（intx=0;xfor（inty=0;yfor（intz=0;z{
if（p{
r[p][0]=（x+0.5）*a;
r[p][1]=（y+0.5）*a;
r[p][2]=（z+0.5）*a;
}
++p;
}
//initializevelocities
for（intp=0;pfor（inti=0;i<3;i++）
v[p][i]=vMax*（2*rand（）/double（RAND_MAX）-1）;
}
voidcomputeAccelerations（）
{
for（inti=0;ifor（intk=0;k<3;k++）
a[i][k]=0;
for（inti=0;ifor（intj=i+1;j{
doublerij[3];//positionofirelativetoj
doublerSqd=0;
for（intk=0;k<3;k++）
{
rij[k]=r[i][k]-r[j][k];
rSqd+=rij[k]*rij[k];
}
doublef=24*（2*pow（rSqd,-7）-pow（rSqd,-4））;
for（intk=0;k<3;k++）
{
a[i][k]+=rij[k]*f;
a[j][k]-=rij[k]*f;
}
}
}
voidvelocityVerlet（doubledt）
{
computeAccelerations（）;
for（inti=0;ifor（intk=0;k<3;k++）
{
r[i][k]+=v[i][k]*dt+0.5*a[i][k]*dt*dt;
v[i][k]+=0.5*a[i][k]*dt;
}
computeAccelerations（）;
for（inti=0;ifor（intk=0;k<3;k++）
v[i][k]+=0.5*a[i][k]*dt;
}
doubleinstantaneousTemperature（）
{
doublesum=0;
for（inti=0;ifor（intk=0;k<3;k++）
sum+=v[i][k]*v[i][k];
returnsum/（3*（N-1））;
}
intmain（）
{
initialize（）;
doubledt=0.01;
ofstreamFileTemp（"Temp.dat"）;
for（inti=0;i<5000;i++）
{
cout<velocityVerlet（dt）;
FileTemp<<<}
FileTemp.close（）;
}
2）面心立方结构
/*
#include
#include
#include
#include
#include
*/
#include
#include
#include
#include
#include
#include
#include
usingnamespacestd;
//simulationparameters
intN=64;//numberofparticles
doublerho=1.0;//density（numberperunitvolume）
doubleT=1.5;//temperature
//functiondeclarations
voidinitialize（）;//allocatesmemory,callsfollowing2functions
voidinitPositions（）;//placesparticlesonanfcclattice，面心立方结构
voidinitVelocities（）;//initialMaxwell-Boltzmannvelocitydistribution
voidrescaleVelocities（）;//adjusttheinstanteoustemperaturetoT
doublegasdev（）;//Gaussiandistributedrandomnumbers
double**r;//positions
double**v;//velocities
double**a;//accelerations
voidinitialize（）
{
r=newdouble*[N];
v=newdouble*[N];
a=newdouble*[N];
for（inti=0;i{
r[i]=newdouble[3];
v[i]=newdouble[3];
a[i]=newdouble[3];
}
initPositions（）;
initVelocities（）;
}
doubleL;//linearsizeofcubicalvolume
voidinitPositions（）
{
//computesideofcubefromnumberofparticlesandnumberdensity
L=pow（N/rho,1.0/3）;
//findMlargeenoughtofitNatomsonanfcclattice
intM=1;
while（4*M*M*M++M;
doublea=L/M;//latticeconstantofconventionalcell
//4atomicpositionsinfccunitcell
doublexCell[4]={0.25,0.75,0.75,0.25};
doubleyCell[4]={0.25,0.75,0.25,0.75};
doublezCell[4]={0.25,0.25,0.75,0.75};
intn=0;//atomsplacedsofar
for（intx=0;xfor（inty=0;yfor（intz=0;zfor（intk=0;k<4;k++）
if（n{
r[n][0]=（x+xCell[k]）*a;
r[n][1]=（y+yCell[k]）*a;
r[n][2]=（z+zCell[k]）*a;
++n;
}
}
doublegasdev（）
{
staticboolavailable=false;
staticdoublegset;
doublefac,rsq,v1,v2;
if（!
available）
{
do
{
v1=2.0*rand（）/double（RAND_MAX）-1.0;
v2=2.0*rand（）/double（RAND_MAX）-1.0;
rsq=v1*v1+v2*v2;
}
while（rsq>=1.0||rsq==0.0）;
fac=sqrt（-2.0*log（rsq）/rsq）;
gset=v1*fac;
available=true;
returnv2*fac;
}
else
{
available=false;
returngset;
}
}
voidinitVelocities（）
{
//Gaussianwithunitvariance
for（intn=0;nfor（inti=0;i<3;i++）
v[n][i]=gasdev（）;
//Adjustvelocitiessocenter-of-massvelocityiszero
doublevCM[3]={0,0,0};
for（intn=0;nfor（inti=0;i<3;i++）
vCM[i]+=v[n][i];
for（inti=0;i<3;i++）
vCM[i]/=N;
for（intn=0;nfor（inti=0;i<3;i++）
v[n][i]-=vCM[i];
//Rescalevelocitiestogetthedesiredinstantaneoustemperature
rescaleVelocities（）;
}
voidrescaleVelocities（）
{
doublevSqdSum=0;
for（intn=0;nfor（inti=0;i<3;i++）
vSqdSum+=v[n][i]*v[n][i];
doublelambda=sqrt（3*（N-1）*T/vSqdSum）;
for（intn=0;nfor（inti=0;i<3;i++）
v[n][i]*=lambda;
}
voidcomputeAccelerations（）
{
for（inti=0;ifor（intk=0;k<3;k++）
a[i][k]=0;
for（inti=0;ifor（intj=i+1;j{
doublerij[3];//positionofirelativetoj
doublerSqd=0;
for（intk=0;k<3;k++）
{
rij[k]=r[i][k]-r[j][k];
//closestimageconvention
if（abs（rij[k]）>0.5*L）
{
if（rij[k]>0）
rij[k]-=L;
else
rij[k]+=L;
}
rSqd+=rij[k]*rij[k];
}
doublef=24*（2*pow（rSqd,-7）-pow（rSqd,-4））;
for（intk=0;k<3;k++）
{
a[i][k]+=rij[k]*f;
a[j][k]-=rij[k]*f;
}
}
}
voidvelocityVerlet（doubledt）
{
computeAccelerations（）;
for（inti=0;ifor（intk=0;k<3;k++）
{
r[i][k]+=v[i][k]*dt+0.5*a[i][k]*dt*dt;
//useperiodicboundaryconditions
if（r[i][k]<0）
r[i][k]+=L;
if（r[i][k]>=L）
r[i][k]-=L;
v[i][k]+=0.5*a[i][k]*dt;
}
computeAccelerations（）;
for（inti=0;ifor（intk=0;k<3;k++）
v[i][k]+=0.5*a[i][k]*dt;
}
doubleinstantaneousTemperature（）
{
doublesum=0;
for（inti=0;ifor（intk=0;k<3;k++）
sum+=v[i][k]*v[i][k];
returnsum/（3*（N-1））;
}
intmain（）
{
initialize（）;
doubledt=0.01;
ofstreamFileTemp（"Temp2.dat"）;
for（inti=0;i<4000;i++）
{
cout<velocityVerlet（dt）;
FileTemp<<<if（i%200==0）
rescaleVelocities（）;
}
FileTemp.close（）;
}
3.计算结果及具体分析讨论
1.两种结构最终结果均收敛到某个值附近，面心立方结构收敛速度更快，并且收敛以后的方差明显更小。
2.该模拟仅符合小角度摆动，摆动角度越大偏离就越大，与实际符合程度也就越差。
升级会员 