分子动力学的模拟程序设计Word下载.docx
- 文档编号:18981629
- 上传时间:2023-01-02
- 格式:DOCX
- 页数:13
- 大小:36.84KB
分子动力学的模拟程序设计Word下载.docx
《分子动力学的模拟程序设计Word下载.docx》由会员分享,可在线阅读,更多相关《分子动力学的模拟程序设计Word下载.docx(13页珍藏版)》请在冰豆网上搜索。
(3)数值计算的舍入误差。
分子动力学的适用范围
原则上,分子动力学方法所适用的微观物理体系并无任何限制。
这个方法适用的体系既可以是少体系统,也可以是多体系统;
既可以是点粒子体系,也可以是具有内部结构的体系;
处理的微观客体既可以是分子,也可以是其它的微观粒子。
微正则系综(NVE)的分子动力学模拟
步骤如下:
(1)给定初始空间位置
(2)计算第n步时粒子所受的力
(3)计算第n+1步粒子的空间位置
(4)计算第n步的速度,根据速度得出系统的温度。
(5)返回第
(2)步,开始下一步计算。
2.源程序
1)简立方结构
/*
#include<
cmath>
cstdlib>
fstream>
iostream>
string>
*/
string.h>
math.h>
iomanip>
stdlib.h>
cstdio>
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<
n;
x++)
for(inty=0;
y<
y++)
for(intz=0;
z<
z++)
{
if(p<
N)
{
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<
N;
p++)
for(inti=0;
i<
3;
i++)
v[p][i]=vMax*(2*rand()/double(RAND_MAX)-1);
}
voidcomputeAccelerations()
i++)//setallaccelerationstozero
for(intk=0;
k<
k++)
a[i][k]=0;
N-1;
i++)//loopoveralldistinctpairsi,j
for(intj=i+1;
j<
j++)
{
doublerij[3];
//positionofirelativetoj
doublerSqd=0;
k++)
rij[k]=r[i][k]-r[j][k];
rSqd+=rij[k]*rij[k];
doublef=24*(2*pow(rSqd,-7)-pow(rSqd,-4));
a[i][k]+=rij[k]*f;
a[j][k]-=rij[k]*f;
voidvelocityVerlet(doubledt)
computeAccelerations();
r[i][k]+=v[i][k]*dt+0.5*a[i][k]*dt*dt;
v[i][k]+=0.5*a[i][k]*dt;
doubleinstantaneousTemperature()
doublesum=0;
sum+=v[i][k]*v[i][k];
returnsum/(3*(N-1));
intmain()
initialize();
doubledt=0.01;
ofstreamFileTemp("
Temp.dat"
);
5000;
{
cout<
<
i<
'
\n'
;
velocityVerlet(dt);
FileTemp<
fixed\
<
setw(15)<
setprecision(5)<
i\
setw(20)<
setprecision(10)<
instantaneousTemperature()<
'
FileTemp.close();
2)面心立方结构
//simulationparameters
intN=64;
doublerho=1.0;
//density(numberperunitvolume)
doubleT=1.5;
//temperature
//functiondeclarations
voidinitialize();
//allocatesmemory,callsfollowing2functions
voidinitPositions();
//placesparticlesonanfcclattice,面心立方结构
voidinitVelocities();
//initialMaxwell-Boltzmannvelocitydistribution
voidrescaleVelocities();
//adjusttheinstanteoustemperaturetoT
doublegasdev();
//Gaussiandistributedrandomnumbers
double**r;
double**v;
double**a;
r=newdouble*[N];
v=newdouble*[N];
a=newdouble*[N];
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<
N)
++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
M;
z++)
4;
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<
n++)
v[n][i]=gasdev();
//Adjustvelocitiessocenter-of-massvelocityiszero
doublevCM[3]={0,0,0};
vCM[i]+=v[n][i];
vCM[i]/=N;
v[n][i]-=vCM[i];
//Rescalevelocitiestogetthedesiredinstantaneoustemperature
rescaleVelocities();
voidrescaleVelocities()
doublevSqdSum=0;
vSqdSum+=v[n][i]*v[n][i];
doublelambda=sqrt(3*(N-1)*T/vSqdSum);
v[n][i]*=lambda;
//closestimageconvention
if(abs(rij[k])>
0.5*L)
{
if(rij[k]>
0)
rij[k]-=L;
else
rij[k]+=L;
//useperiodicboundaryconditions
if(r[i][k]<
r[i][k]+=L;
if(r[i][k]>
=L)
r[i][k]-=L;
Temp2.dat"
4000;
if(i%200==0)
3.计算结果及具体分析讨论
1.两种结构最终结果均收敛到某个值附近,面心立方结构收敛速度更快,并且收敛以后的方差明显更小。
2.该模拟仅符合小角度摆动,摆动角度越大偏离就越大,与实际符合程度也就越差。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 分子 动力学 模拟 程序设计