计算物理实验报告

计算物理实验报告实验内容MD算法及程序专业凝聚态姓名王伟学号SX1408035指导教师丁长庚日期2014-11-25计算物理实验报告-1-一、目的意义(1)了解基本的分子动力学理论。(2)进一步掌握熟悉Verlet算法、Lennard-Jones势以及分子演化过程。(3)编写二维分子动力学演化程序，并计算相应的物理量。二、内容要求编写一个简单的二维MD程序。初始时将给定的粒子放置到正方晶格上，每个粒子速度按随机分布（-0.5,0.5）区间。用Lennard-Jones势模拟，用Verlet算法计算演化，输出粒子动能、势能、总能随时间变化曲线。参数：Npart=256T=1box=16Tmax=0.008deltaT=0.0000008三、问题解决的方法与算法方法：Lennard-Jones势、Verlet算法、周期性边界条件。1、Lennard-Jones势；1261(|rr|)=4*()ijijijijUurr对于不同分子动力学元胞内的粒子，如果相互作用是短程力我们可以在长度=2cboxr处截断，这时为了使势能可导，(r)cU必须要足够小，以使截断不会显著地影响模拟结果，则增加一项截断补偿12614*()cccUrr即是0ijcUUU||||ijcijcrrrrrr，由势能可求得两分子间的作用力，即()xUxUfrxrr2、Verlet算法；已经知道力，用牛顿定律和级数展开32'''4()t(tt)(t)v(t)tt(t)3!ftrrrm32'''4()t(t-t)(t)-v(t)tt-(t)3!ftrrrm243()(tt)2(t)-(t-t)t(t)(tt)-(t-t)(t)(t)2tftrrrmrrv计算物理实验报告-2-3、周期性边界条件；为了将分子动力学在元胞有限矩形内的模拟，扩展到真实大系统的模拟，我们通常采用周期性边界条件。采用这种边界条件，我们可以消除引入元胞后的表面效应，构造出一个准无穷大的面积来精确地代表宏观系统。该边界条件的具体操作：当有一个粒子穿过基本分子动力学元胞的矩形时，就让这个粒子以相同的速度穿过此边对边，重新进入元胞内，以此时的距离代表实际的距离来模拟计算Lennard-Jones势和力。即是X方向，xr1=xr1-box*int(xr1/box)Y方向，xr2=xr2-box*int(xr2/box)4、在计算机上对分子系统的模拟的实际步骤；Step1：设定模拟采用的模型。Step2：给定初始条件。Setp3：趋于平衡的演化计算过程。Setp4：演化的同时计算相应的物理量。Setp5：做出相应的图形。四、计算程序(Fortran)programmddoubleprecision::t=0,tmax=0.008doubleprecision::en,etot,temp,f(2,256)callinitopen(1,file=temp.txt)open(2,file=en.txt)open(3,file=etot.txt)dowhile(t.lt.tmax)callforce(f,en)callintergrate(f,en,etot,temp)write(1,*)t,tempwrite(2,*)t,enwrite(3,*)t,etott=t+0.0000008enddoendmodulevariableimplicitnoneinteger::i,j,npart=256doubleprecision,dimension(2)::sumvdoubleprecision::sumv2doubleprecision::xm(2,256),x(2,256),v(2,256),v2(256)doubleprecision::box=16,delta=0.0000008,rc,rc2,ecutendsubroutineinit计算物理实验报告-3-usevariabledoubleprecision::fsdoi=1,npartx(1,i)=mod(i-1,16)x(2,i)=(i-1)/16callrandom_number(ranvel)v(1,i)=ranvel-0.5callrandom_number(ranvel)v(2,i)=ranvel-0.5v2(i)=v(1,i)**2+v(2,i)**2enddosumv=sum(v)/npartsumv2=sum(v2)/npartfs=sqrt(3*1.0d0/sumv2)doi=1,npartv(:,i)=(v(:,i)-sumv)*fsxm(:,i)=x(:,i)-delta*v(:,i)enddorc=real(box)/2.0d0rc2=rc**2ecut=4*(1.0d0/rc2**3-1.0d0)/rc2**3returnendsubroutineforce(f,en)usevariabledoubleprecision::en,xr1,xr2,r2,r2i,r6idoubleprecision::f(2,256)en=0doi=1,npartf(1,i)=0f(2,i)=0enddodoi=1,npart-1doj=i+1,npartxr1=x(1,i)-x(1,j)xr2=x(2,i)-x(2,j)xr1=xr1-box*int(xr1/box)xr2=xr2-box*int(xr2/box)r2=xr1**2+xr2**2if(r2.lt.rc2)thenr2i=1.0d0/r2r6i=r2i**3计算物理实验报告-4-ff=48*r2i*r6i*(r6i-0.5)f(1,i)=f(1,i)+ff*xr1f(2,i)=f(2,i)+ff*xr2f(1,j)=f(1,j)-ff*xr1f(2,j)=f(2,j)-ff*xr2en=en+4*r6i*(r6i-1)-ecutendifenddoenddoen=en/npartreturnendsubroutineintergrate(f,en,etot,temp)usevariabledoubleprecision::en,etot,temp,f(2,256)doubleprecision,dimension(2)::xxdoi=1,npartxx=2*x(:,i)-xm(:,i)+delta**2*f(:,i)v(:,i)=(xx-xm(:,i))/(2*delta)v2(i)=v(1,i)**2+v(2,i)**2xm(:,i)=x(:,i)x(:,i)=xxenddosumv=sum(v)sumv2=sum(v2)temp=0.5*sumv2/npartetot=en+tempreturnend五、计算结果与分析计算物理实验报告-5-1.平均动能随时间演化图；0.0000.0020.0040.0060.0081.421.431.441.451.461.471.481.491.501.51KineticEnergyTime2.平均势能随时间演化图；0.0000.0020.0040.0060.008-1.03-1.02-1.01-1.00-0.99-0.98-0.97-0.96-0.95EnergyTime计算物理实验报告-6-3.平均总能随时间演化图；0123456789x10-30.47260.47260.47260.47260.47260.47260.47260.47260.47260.47260.4726