计算物理实验报告

计算物理实验报告

实验内容MD算法及程序

专业凝聚态

姓名王伟

学号SX1408035

指导教师丁长庚

日期2014-11-25

一、目的意义

(1)了解基本的分子动力学理论。

(2)进一步掌握熟悉Verlet 算法、Lennard-Jones 势以及分子演化过程。 (3)编写二维分子动力学演化程序，并计算相应的物理量。

二、内容要求

编写一个简单的二维MD 程序。初始时将给定的粒子放置到正方晶格上，每个粒子速度按随机分布（-0.5,0.5）区间。用Lennard-Jones 势模拟，用Verlet 算法计算演化，输出粒子动能、势能、总能随时间变化曲线。

参数：Npart=256 T=1 box=16

Tmax=0.008 deltaT=0.0000008

三、问题解决的方法与算法

方法：Lennard-Jones 势、Verlet 算法、周期性边界条件。

1、Lennard-Jones 势；

126

1

(|r r |)=4*(

)ij i j i

j

i j

U u r r =--∑∑ 对于不同分子动力学元胞内的粒子，如果相互作用是短程力我们可以

在长度=2

c box

r 处截断，这时为了使势能可导，(r )c U 必须要足够小，以

使截断不会显著地影响模拟结果，则增加一项截断补偿126

1

4*(

)c c c U r r =-

即是0ij c

U U U -?=?

?

||||i j c i j c r r r r r r -≤->，由势能可求得两分子间的作用力，即 ()x U x U

f r x r r

??=-

=-?? 2、 Verlet 算法；

已经知道力，用牛顿定律和级数展开

32'''

4()t (t t )(t )v (t )t t (

t )

3!f t r r r m ?+?=

+?+?++O ?

32

'''

4

()

t (t -t )(t )-v (t )t

t -(t )

3!

f t r r r m

??

=?+?+O ? 243()(t t)2(t)-(t -t)t (t )(t t)-(t -t)(t)(t )2t f t r r r m

r r v ?

+?=?+?+O ?????+???=+O ????

3、周期性边界条件；

为了将分子动力学在元胞有限矩形内的模拟，扩展到真实大系统的模拟，

我们通常采用周期性边界条件。采用这种边界条件，我们可以消除引入元

胞后的表面效应，构造出一个准无穷大的面积来精确地代表宏观系统。该

边界条件的具体操作：当有一个粒子穿过基本分子动力学元胞的矩形时，

就让这个粒子以相同的速度穿过此边对边，重新进入元胞内，以此时的距

离代表实际的距离来模拟计算Lennard-Jones势和力。即是

X方向，xr1=xr1-box*int(xr1/box)

Y方向，xr2=xr2-box*int(xr2/box)

4、在计算机上对分子系统的模拟的实际步骤；

Step1：设定模拟采用的模型。

Step2：给定初始条件。

Setp3：趋于平衡的演化计算过程。

Setp4：演化的同时计算相应的物理量。

Setp5：做出相应的图形。

四、计算程序(Fortran)

program md

double precision :: t=0,tmax=0.008

double precision :: en,etot,temp,f(2,256)

call init

open(1,file="temp.txt")

open(2,file="en.txt")

open(3,file="etot.txt")

do while(t.lt.tmax)

call force(f,en)

call intergrate(f,en,etot,temp)

write(1,*) t, temp

write(2,*) t, en

write(3,*) t, etot

t=t+0.0000008

end do

end

module variable

implicit none

integer :: i,j,npart=256

double precision, dimension(2) :: sumv

double precision :: sumv2

double precision :: xm(2,256),x(2,256),v(2,256),v2(256)

double precision :: box=16,delta=0.0000008,rc,rc2,ecut

end

subroutine init

use variable

double precision :: fs

do i=1, npart

x(1,i)=mod(i-1,16)

x(2,i)=(i-1)/16

call random_number(ranvel)

v(1,i) = ranvel -0.5

call random_number(ranvel)

v(2,i) = ranvel -0.5

v2(i)=v(1,i)**2+v(2,i)**2

end do

sumv=sum(v)/npart

sumv2=sum(v2)/npart

fs=sqrt(3*1.0d0/sumv2)

do i=1,npart

v(:,i)=(v(:,i)-sumv)*fs

xm(:,i)=x(:,i)-delta*v(:,i)

end do

rc=real(box)/2.0d0

rc2=rc**2

ecut=4*(1.0d0/rc2**3-1.0d0)/rc2**3 return

end

subroutine force(f,en)

use variable

double precision :: en,xr1,xr2,r2,r2i,r6i double precision :: f(2,256)

en=0

do i=1,npart

f(1,i)=0

f(2,i)=0

end do

do i=1,npart-1

do j=i+1,npart

xr1=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**2

if(r2.lt.rc2) then

r2i=1.0d0/r2

r6i=r2i**3

ff=48*r2i*r6i*(r6i-0.5)

f(1,i)=f(1,i)+ff*xr1

f(2,i)=f(2,i)+ff*xr2

f(1,j)=f(1,j)-ff*xr1

f(2,j)=f(2,j)-ff*xr2

en=en+4*r6i*(r6i-1)-ecut

end if

end do

end do

en=en/npart

return

end

subroutine intergrate(f,en,etot,temp)

use variable

double precision :: en,etot,temp,f(2,256)

double precision, dimension(2) :: xx

do i=1,npart

xx=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)**2

xm(:,i)=x(:,i)

x(:,i)=xx

end do

sumv=sum(v)

sumv2=sum(v2)

temp=0.5*sumv2/npart

etot=en+temp

return

end

五、计算结果与分析

1.平均动能随时间演化图；

0.000

0.002

0.004

0.006

0.008

1.42

1.431.441.451.461.471.481.49

1.501.51K i n e t i c E n e r g y

Time

2.平均势能随时间演化图；

0.0000.0020.0040.0060.008

-1.03

-1.02-1.01-1.00-0.99-0.98

-0.97-0.96-0.95E n e r g y

Time

3.平均总能随时间演化图；

0123456789

x 10-3