分子動力學(xué)的模擬程序設(shè)計

1、分子動力學(xué)的模擬程序設(shè)計一、課題名稱：分子動力學(xué)的模擬程序設(shè)計二、班級和姓名：*三、主要內(nèi)容：1. 研究的內(nèi)容和算法分子動力學(xué)模擬方法確定性的演化過程。按該體系內(nèi)部的內(nèi)稟動力學(xué)規(guī)律計算、確定位形的轉(zhuǎn)變。出發(fā)點：物理系統(tǒng)確定的微觀描述。具體的說就是描述系統(tǒng)的哈密頓、拉格朗日或者牛頓運動方程，用這些方程去驅(qū)動粒子的位置、速度和取向隨時間的演化。模擬過程：(1) 通過對物理體系的微觀數(shù)學(xué)描述建立一組方程組，直接求解每個分子的運動方程。得到每個時刻、每個分子的坐標和動量。(2) 利用統(tǒng)計方法計算系統(tǒng)的靜態(tài)和動態(tài)過程。計算機模擬的問題1 觀測時間是有限的對于某些問題，有限觀測時間可以看成是無限長的。如對

2、一個分子系統(tǒng)的計算，它的觀測時間遠大于分子時間尺度。2 觀測體系是有限的。引入周期性、全反射、漫反射等邊界條件。分子動力學(xué)模擬的計算對體系的分子動力學(xué)方程組進行求解時，需要將運動的連續(xù)性方程離散化變成有限差分方程，常用的方法有：歐拉法，龍格-庫塔法、誤差來源：(1)動力學(xué)模型的近似程度；(2)數(shù)值求解法的近似階數(shù)；(3) 數(shù)值計算的舍入誤差。分子動力學(xué)的適用范圍原則上，分子動力學(xué)方法所適用的微觀物理體系并無任何限制。這個方法適用的體系既可以是少體系統(tǒng)，也可以是多體系統(tǒng)；既可以是點粒子體系，也可以是具有內(nèi)部結(jié)構(gòu)的體系；處理的微觀客體既可以是分子，也可以是其它的微觀粒子。微正則系綜(NVE)的分子

3、動力學(xué)模擬步驟如下：(1)給定初始空間位置(r(0),(r(1)ii計算第n步時粒子所受的力F=Fitni(3) 計算第n+1步粒子的空間位置r(M=2(n)-r(T+F(n)加/miiii(4) 計算第n步的速度v(”)=(r(”+1)-r(”-1)/2h，根據(jù)速度得iii出系統(tǒng)的溫度。(5) 返回第(2)步，開始下一步計算。2.源程序1)簡立方結(jié)構(gòu)/*#include<cmath>#include<cstdlib>#include<fstream>#include<iostream>#include<string>*/#inclu

4、de<iostream>#include<string.h>#include<math.h>#include<iomanip>#include<stdlib.h>#include<fstream>#include<cstdio>usingnamespacestd;constintN=64;doublerN3;doublevN3;doubleaN3;/numberofparticles/positions/velocities/accelerationsdoubleL=10;doublevMax=0.1;/lin

5、earsizeofcubicalvolume，簡單立方/maximuminitialvelocitycomponentvoidinitialize()/initializepositionsintn=int(ceil(pow(N,1.0/3);/numberofatomsineachdirectiondoublea=L/n;/latticespacing，晶格常數(shù)intp=0;/particlesplacedsofarfor(intx=0;x<n;x+)for(inty=0;y<n;y+)for(intz=0;z<n;z+)if(p<N)rp0=(x+0.5)*a;rp

6、1=(y+0.5)*a;rp2=(z+0.5)*a;+p;/initializevelocitiesfor(intp=0;p<N;p+)for(inti=0;i<3;i+)vpi=vMax*(2*rand()/double(RAND_MAX)-1);/setallaccelerationstozerovoidcomputeAccelerations()for(inti=0;i<N;i+)for(intk=0;k<3;k+)aik=0;for(inti=0;i<N-1;i+)for(intj=i+1;j<N;j+)/loopoveralldistinctpai

7、rsi,jdoublerij3;/positionofirelativetojdoublerSqd=0;for(intk=0;k<3;k+)rijk=rik-rjk;rSqd+=rijk*rijk;doublef=24*(2*pow(rSqd,-7)-pow(rSqd,-4);for(intk=0;k<3;k+)aik+=rijk*f;ajk-=rijk*f;voidvelocityVerlet(doubledt)computeAccelerations();for(inti=0;i<N;i+)for(intk=0;k<3;k+)rik+=vik*dt+0.5*aik*

8、dt*dt;vik+=0.5*aik*dt;computeAccelerations();for(inti=0;i<N;i+)for(intk=0;k<3;k+)vik+=0.5*aik*dt;doubleinstantaneousTemperature()doublesum=0;for(inti=0;i<N;i+)for(intk=0;k<3;k+)sum+=vik*vik;returnsum/(3*(N-1);intmain()initialize();doubledt=0.01;ofstreamFileTemp("Temp.dat");for(

9、inti=0;i<5000;i+)cout<<i<<'n'velocityVerlet(dt);FileTemp<<fixed<<setw(15)<<setprecision(5)<<i<<setw(20)<<setprecision(10)<<instantaneousTemperature()<<'n'FileTemp.close();2)面心立方結(jié)構(gòu)/*#include<cmath>#include<cstdlib

10、>#include<fstream>#include<iostream>#include<string>*/#include<iostream>#include<string.h>#include<math.h>#include<iomanip>#include<stdlib.h>#include<fstream>#include<cstdio>usingnamespacestd;/simulationparametersintN=64;doublerho=1.0;do

11、ubleT=1.5;/numberofparticles/density(numberperunitvolume)/temperature/functiondeclarationsvoidinitialize();/allocatesmemory,callsfollowing2functionsvoidinitPositions();/placesparticlesonanfeelattice,面心立方結(jié)構(gòu)voidinitVelocities();/initialMaxwell-BoltzmannvelocitydistributionvoidrescaleVelocities();/adju

12、sttheinstanteoustemperaturetoTdoublegasdev();/Gaussiandistributedrandomnumbersdouble*r;double*v;double*a;/positions/velocities/accelerationsvoidinitialize()r=newdouble*N;v=newdouble*N;a=newdouble*N;for(inti=0;i<N;i+)ri=newdouble3;vi=newdouble3;ai=newdouble3;initPositions();initVelocities();double

13、L;/linearsizeofcubicalvolumevoidinitPositions()/computesideofcubefromnumberofparticlesandnumberdensityL=pow(N/rho,1.0/3);/findMlargeenoughtofitNatomsonanfcclatticeintM=1;while(4*M*M*M<N)+M;doublea=L/M;/latticeconstantofconventionalcell/4atomicpositionsinfccunitcelldoublexCell4=0.25,0.75,0.75,0.25

14、;doubleyCell4=0.25,0.75,0.25,0.75;doublezCell4=0.25,0.25,0.75,0.75;intn=0;/atomsplacedsofarfor(intx=0;x<M;x+)for(inty=0;y<M;y+)for(intz=0;z<M;z+)for(intk=0;k<4;k+)if(n<N)rn0=(x+xCellk)*a;rn1=(y+yCellk)*a;rn2=(z+zCellk)*a;+n;doublegasdev()staticboolavailable=false;staticdoublegset;doub

15、lefac,rsq,v1,v2;if(!available)dov1=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;elseavailable=false;returngset;voidinitVelocities()/Gaussianwithunitvariancefor(intn=0;

16、n<N;n+)for(inti=0;i<3;i+)vni=gasdev();/Adjustvelocitiessocenter-of-massvelocityiszerodoublevCM3=0,0,0;for(intn=0;n<N;n+)for(inti=0;i<3;i+)vCMi+=vni;for(inti=0;i<3;i+)vCMi/=N;for(intn=0;n<N;n+)for(inti=0;i<3;i+)vni-=vCMi;/Rescalevelocitiestogetthedesiredinstantaneoustemperaturere

17、scaleVelocities();voidrescaleVelocities()doublevSqdSum=0;for(intn=0;n<N;n+)for(inti=0;i<3;i+)vSqdSum+=vni*vni;doublelambda=sqrt(3*(N-1)*T/vSqdSum);for(intn=0;n<N;n+)for(inti=0;i<3;i+)vni*=lambda;/setallaccelerationstozerovoidcomputeAccelerations()for(inti=0;i<N;i+)for(intk=0;k<3;k+

18、)aik=0;for(inti=0;i<N-1;i+)for(intj=i+1;j<N;j+)/loopoveralldistinctpairsi,jdoublerij3;/positionofirelativetojdoublerSqd=0;for(intk=0;k<3;k+)rijk=rik-rjk;/closestimageconventionif(abs(rijk)>0.5*L)if(rijk>0)rijk-=L;elserijk+=L;rSqd+=rijk*rijk;doublef=24*(2*pow(rSqd,-7)-pow(rSqd,-4);for(intk=0;k<3;k+)aik+=rijk*f;ajk-=rijk*f;voidvelocityVerlet(doubledt)computeAccelerations();for(inti=0;i<N;i+)for(intk=0;k<3;k+)rik+=vik*dt+0.5*aik*dt*dt;/useperiodicboundaryconditionsif(rik<0)rik+=L;if(rik>=L)rik-=L;vik+=0.5*aik*dt;computeAccelerations();for(in