分子动力学的模拟程序设计.docx
- 文档编号:5988225
- 上传时间:2023-01-02
- 格式:DOCX
- 页数:13
- 大小:36.84KB
分子动力学的模拟程序设计.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;x for(inty=0;y for(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;p for(inti=0;i<3;i++) v[p][i]=vMax*(2*rand()/double(RAND_MAX)-1); } voidcomputeAccelerations() { for(inti=0;i for(intk=0;k<3;k++) a[i][k]=0; for(inti=0;i for(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;i for(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;i for(intk=0;k<3;k++) v[i][k]+=0.5*a[i][k]*dt; } doubleinstantaneousTemperature() { doublesum=0; for(inti=0;i for(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;x for(inty=0;y for(intz=0;z for(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;n for(inti=0;i<3;i++) v[n][i]=gasdev(); //Adjustvelocitiessocenter-of-massvelocityiszero doublevCM[3]={0,0,0}; for(intn=0;n for(inti=0;i<3;i++) vCM[i]+=v[n][i]; for(inti=0;i<3;i++) vCM[i]/=N; for(intn=0;n for(inti=0;i<3;i++) v[n][i]-=vCM[i]; //Rescalevelocitiestogetthedesiredinstantaneoustemperature rescaleVelocities(); } voidrescaleVelocities() { doublevSqdSum=0; for(intn=0;n for(inti=0;i<3;i++) vSqdSum+=v[n][i]*v[n][i]; doublelambda=sqrt(3*(N-1)*T/vSqdSum); for(intn=0;n for(inti=0;i<3;i++) v[n][i]*=lambda; } voidcomputeAccelerations() { for(inti=0;i for(intk=0;k<3;k++) a[i][k]=0; for(inti=0;i for(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;i for(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;i for(intk=0;k<3;k++) v[i][k]+=0.5*a[i][k]*dt; } doubleinstantaneousTemperature() { doublesum=0; for(inti=0;i for(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.该模拟仅符合小角度摆动,摆动角度越大偏离就越大,与实际符合程度也就越差。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 分子 动力学 模拟 程序设计
![提示](https://static.bdocx.com/images/bang_tan.gif)