program microcanonicalmoleculardynamics
c=========================================================================
c
c molecular dynamics
c
c microcanonical ensemble
c
c application to argon. the lennard-jones potential is truncated at
c rcoff and not smoothly continued to zero. initial the npart particles
c are placed on an fcc lattice. the velocities are drawn from a boltzmann
c distribution with temperature tref.
c
c input parameters are as follews
c
c npart numberof particles(must be a multiple of 4)
c side side length of the cubical box in sigma units
c tref reduced temperature
c rcoff cutoff of the potential in sigma units
c h basic time step
c irep velocities scaling every irep'th time step
c istop sto p of scaling of the velocities at istop
c timemx number of integration steps
c iseed seed for the random number generator
c
c version 1.0 as of august 1985
c
c dieter w.heermann
c
c=========================================================================
c
real x(1:768),vh(1:768),f(1:768)
c
real fy(1:256),fz(1:256),y(1:256),z(1:256)
real potential(1:1000000),force(1:1000000)
c
real h,hsq,hsq2,tref
real kx,ky,kz
integer(2) ihour, iminite ,isecond, imsecond
integer(2) ehour, eminite ,esecond, emsecond
real delttime,etime,itime
c
integer(2) fihour, fiminite ,fisecond, fimsecond
integer(2) fehour, feminite ,fesecond, femsecond
c
integer(2) fihourt, fiminitet ,fisecondt, fimsecondt
integer(2) fehourt, feminitet ,fesecondt, femsecondt
real fdelttime,fetime,fitime
c
integer clock,timemx
c
equivalence (fy,f(257)),(fz,f(513)),(y,x(257)),(z,x(513))
c
open(unit=6,file='debug.txt',status='new')
open(unit=7,file='energy.txt',status='new')
open(unit=8,file='temppres.txt',status='new')
open(unit=13,file='lattice.lsp',status='new')!,
c * form='formatted',action='write')
!open(unit=9,file='maxwell.txt',status='new')
!open(unit=10,file='rannum.txt',status='new')
c ,access='sequential',form='formatted')
c
c-------------------------------------------------------------------------
c
c definition of the simulation parameters
c
c-------------------------------------------------------------------------
c
npart =256
den =0.83134
side =6.75284
tref =0.722
rcoff =3.0
h =0.032
irep =200
istop =5000
timemx =10000
iseed= 4711
c
c-------------------------------------------------------------------------
c
c set the other parameters
c
c-------------------------------------------------------------------------
c
write(6,*) 'molecular dynamics simulation program'
write(6,*) '-------------------------------------'
write(6,*)
write(6,*) 'munber of particles is ',npart
write(6,*) 'side length of the box is',side
write(6,*) 'cut off is ',rcoff
write(6,*) 'reduced temperature is ',tref
write(6,*) 'basic time step is
',h
write(6,*) 'density ',den
write(6,*) 'istop ',istop
write(6,*) 'irep ',irep
write(6,*) 'timemx ',timemx
write(6,*) 'iseed ',iseed
write(6,*)
c
a =side/4.0
sideh =side*0.5
hsq =h*h
hsq2 =hsq*0.5
npartm =npart-1
rcoffs =rcoff*rcoff
tscale =16.0/(1.0*npart-1.0)
vaver =1.13*sqrt(tref/24.0)
alat =a/1.2
c
n3 =3*npart
iof1 =npart
iof2 =2*npart
c
write(6,*)
write(6,*) 'a =side/4.0 ',a
write(6,*) 'sideh =side*0.5',sideh
write(6,*) 'hsq =h*h ',hsq
write(6,*) 'hsq2 =hsq*0.5 ',hsq2
write(6,*) 'npartm =npart-1',npartm
write(6,*) 'rcoffs =rcoff*rcoff',rcoffs
write(6,*) 'tscale =16.0/(1.0*npart-1.0)',tscale
write(6,*) 'vaver =1.13*sqrt(tref/24.0)',vaver
write(6,*) 'n3 =3*npart',n3
write(6,*) 'iof1 =npart',iof1
write(6,*) 'iof2 =2*npart',iof2
write(6,*)
c call ranset(iseed)
c
c
c-------------------------------------------------------------------------
c
c this part of the program prepares the initial configuration
c
c-------------------------------------------------------------------------
c
c
c set up fcc lattice for the atoms inside the box
c
c==========================================================================
ihour=0
iminite=0
isecond=0
imsecond=0
ehour=0
eminite=0
esecond=0
emsecond=0
delttime=0.0
etime=0.0
itime=0.0
c
fihourt=0
fiminitet=0
fisecondt=0
fimsecondt=0
fehourt=0
feminitet=0
fesecondt=0
femsecondt=0
fdelttime=0.0
fetime=0.0
fitime=0.0
iw=0
call gettim(ihour, iminite ,isecond,imsecond)
pstep=(rcoffs-0.500)/1000000
do 199 i=1,1000001
rd=0.500+pstep*(i-1)
potential(i)=rd**(-6.0)-rd**(-3.0)
force(i)=rd**(-7.0)-0.5*rd**(-4.0)
!write(6,*)potential(i),force(i)
199 continue
write(13,*) '(defun lattice()'
!write(13,123) boxsize
c123 !format('(command "box" (list 0 0 0)',1x,'(list ',3e23.8,') "")')
ijk=0
do 10 lg=0,1
do 10 i=0,3
do 10 j=0,3
do 10 k=0,3
ijk=ijk+1
x(ijk)=i*a+lg*a*0.5
y(ijk)=j*a+lg*a*0.5
z(ijk)=k*a
iw=iw+1
write(13,122) x(ijk),y(ijk),z(ijk),alat/2.0
!122 format('(command "sphere" (list ',3f23.8,')',f23.8,')')
!write(6,*) 'x(',ijk,')=',x(ijk)
!write(6,*) 'y(',ijk,')=',y(ijk)
!write(6,*) 'z(',ijk,')=',z(ijk)
10 continue
do 15 lg=1,2
do 15 i=0,3
do 15 j=0,3
do 15 k=0,3
ijk=ijk+1
x(ijk)=i*a+(2-lg)*a*0.5
y(ijk)=j*a+(lg-1)*a*0.5
z(ijk)=k*a+a*0.5
iw=iw+1
write(13,122) x(ijk),y(ijk),z(ijk),alat/2.0
122 format('(command "sphere" (list ',3f12.6,')',f12.6,')')
!write(6,*) 'x(',ijk,')=',x(ijk)
!write(6,*) 'y(',ijk,')=',y(ijk)
!write(6,*) 'z(',ijk,')=',z(ijk)
15 continue
write(
*,*) iw
write(13,*) ')'
close(13)
c
c assign velocities distributed normally
c
c==========================================================================
call mxwell(vh,n3,h,tref,iseed)
c
do 50 i=1,n3
f(i)=0.0
50 continue
c
c
c-------------------------------------------------------------------------
c
c start of the actual molecular dynamics program.
c
c the equation of motion are integrated using the 'summed form'
c (g.dahlquist and a.bjoerk, numerical methods, prentice hall
c (1974)). every irep'th step the velocities are rescaled so as
c to give the specfied temperature tref.
c
c version 1.0 august 1985
c dieter w.heermann
c
c
c-------------------------------------------------------------------------
c
do 200 clock=1,timemx
write(*,*) clock
c
c =====advance positions one
c basic time step===========
do 210 i=1,n3
x(i)=x(i)+vh(i)+f(i)
210 continue
c
c =============apply periodic boundary conditions==============
do 215 i=1,n3
if (x(i).lt.0) x(i)=x(i)+side
if (x(i).gt.side) x(i)=x(i)-side
215 continue
c
c =============compute the partial velocities==============
do 220 i=1,n3
vh(i)=vh(i)+f(i)
220 continue
c
c
c =========================================================
c = =
c = this part computes the forces on the particles =
c = =
c =========================================================
c
c
vir=0.0
epot=0.0
ipotential=0
do 225 i=1,n3
f(i)=0.0
225 continue
c
c
do 270 i=1,npart
xi=x(i)
yi=y(i)
zi=z(i)
do 270 j=i+1,npart
xx=xi-x(j)
yy=yi-y(j)
zz=zi-z(j)
c
if (xx.lt.-sideh) xx=xx+side
if (xx.gt. sideh) xx=xx-side
c
if (yy.lt.-sideh) yy=yy+side
if (yy.gt. sideh) yy=yy-side
c
if (zz.lt.-sideh) zz=zz+side
if (zz.gt. sideh) zz=zz-side
rd=xx*xx+yy*yy+zz*zz
!write(6,*) xx,yy,zz,rd
if (rd.gt.rcoffs) goto 270
c
call gettim(fihour, fiminite ,fisecond,fimsecond)
fihourt=fihour+fihourt
fiminitet=fiminite+fiminitet
fisecondt=fisecond+fisecondt
fimsecondt=fimsecond+fimsecondt
c
ipotential=int((rd-0.500)/pstep+0.5)
iw=iw+1
!write(6,*) ipotential
epot=epot+potential(ipotential)
r148=force(ipotential)
c
call gettim(fehour,feminite,fesecond,femsecond)
fehourt=fehour+fehourt
feminitet=feminite+feminitet
fesecondt=fesecond+fesecondt
femsecondt=femsecond+femsecondt
c
vir =vir-rd*r148
kx=xx*r148
f(i)=f(i)+kx
f(j)=f(j)-kx
ky=yy*r148
fy(i)=fy(i)+ky
fy(j)=fy(j)-ky
kz=zz*r148
fz(i)=fz(i)+kz
fz(j)=fz(j)-kz
270 continue
c
do 275 i=1,n3
f(i)=f(i)*hsq2
275 continue
c
c
c =========================================================
c = =
c = end of the force calculation =
c = =
c =========================================================
c
c
c ==================compute the velocities=================
do 300 i=1,n3
vh(i)=vh(i)+f(i)
300 continue
c
c ================compute the kinetic energy===============
ekin=0
do 305 i=1,n3
ekin=ekin+vh(i)*vh(i)
305 continue
ekin=ekin/hsq
c
c ================compute the average velocity=============
vel=0.0
count=0.0
do 306 i=1,npart
vx=vh(i)*vh(i)
vy=vh(i+iof1)*vh(i+iof1)
vz=vh(i+iof2)*vh(i+iof2)
sq=sqrt(vx+vy+vz)
sqt=sq/h
if (sqt.gt.vaver) count=count+1
vel=vel+sq
306 continue
vel=vel/h
c
if (clock.le.istop) then
c ==========normalize the velocities to obtain the========
c =============specified reference tem perature===========
if ((clock/irep)*irep.eq.clock) then
!write(6,*) 'velocity adjustment'
ts=tscale*ekin
!write(6,*) 'temperture before scaling is',ts
sc=tref/ts
sc=sqrt(sc)
!write(6,*) 'scale factor is',sc
do 310 i=1,n3
vh(i)=vh(i)*sc
310 continue
ekin=tref/tscale
endif
endif
c
c
c =========================================================
c = =
c = compute various quantities =
c = =
c =========================================================
c
c
if ((clock.lt.6000)) then !.and.((clock/2)*2.eq.clock)
ek=24.0*ekin
epot=4.0*epot
etot=ek+epot
temp=tscale*ekin
pres=den*16.0*(ekin-vir)/npart
vel =vel/npart
rp =(count/256.0)*100.0
write(7,*) clock,ek,epot,etot,temp
write(8,*) clock,temp,pres,vel,rp
endif
c
200 continue
c
call gettim(ehour, eminite ,esecond,emsecond)
c
etime=ehour*3600+eminite*60+esecond+emsecond*0.01
itime=ihour*3600+iminite*60+isecond+imsecond*0.01
delttime=etime-itime
fetime=fehourt*3600+feminitet*60+fesecondt+femsecondt*0.01
fitime=fihourt*3600+fiminitet*60+fisecondt+fimsecondt*0.01
!111=(fehourt-fihourt)*3600
!222=(feminitet-fiminitet)*60
!333=fesecondt-fisecondt
!444=(femsecondt-fimsecondt)*0.01
!fdelttime=111+222+333+444
fdelttime=fetime-fitime
c
write(6,*) 'total steps:',timemx, ' md steps'
!write(6,*) 'initial time:',ihour, iminite ,isecond,imsecond
!write(6,*) 'end time:',ehour, eminite ,esecond,emsecond
write(6,*) 'total time:',delttime,' seconds'
write(6,*) 'force time:',fdelttime,' seconds'
write(6,*) ' ratio:',fdelttime/delttime*100
write(6,*) ' times:',iw
stop
end
c
c
c===================================================================
c
c mxwell returns upon calling a maxwell distributed deviates
c for the specified temperature tref. all deviates are scaled
c by a factor
c
c calling parameters are as follows
c
c vh vector of length npart (must be a multiple of 2)
c n3 vector
length
c h scaling factor with which all elements of vh are multiplied
c tref temperature
c
c version 1.0 as of august 1985
c dieter w.heermann
c
c===================================================================
c
c
subroutine mxwell(vh,n3,h,tref,iseed)
c
real vh(1:n3)
c
real u1,u2,v1,v2,s,r
c
npart =n3/3
iof1 =npart
iof2 =2*npart
tscale=16.0/(1.0*npart-1.0)
c
do 10 i=1,n3,2
c
1 u1=ranf(iseed)
u2=ranf(iseed)
!write(10,*) u1,u2
c
v1=2.0*u1-1.0
v2=2.0*u2-1.0
s=v1*v1+v2*v2
c
if (s.ge.1.0) goto 1
r =-2.0*alog(s)/s
vh(i) =v1*sqrt(r)
vh(i+1) =v2*sqrt(r)
!write(9,*) vh(i),vh(i+1)
10 continue
c
ekin=0.0
sp=0.0
do 20 i=1,npart
sp=sp+vh(i)
20 continue
sp=sp/npart
do 21 i=1,npart
vh(i)=vh(i)-sp
ekin=ekin+vh(i)*vh(i)
21 continue
!write(6,*)'total linear momentum in x direction ',sp
sp=0.0
do 22 i=iof1+1,iof2
sp=sp+vh(i)
22 continue
sp=sp/npart
do 23 i=iof1+1,iof2
vh(i)=vh(i)-sp
ekin=ekin+vh(i)*vh(i)
23 continue
!write(6,*) 'total linear momentum in y direction ',sp
sp=0.0
do 24 i=iof2+1,n3
sp=sp+vh(i)
24 continue
sp=sp/npart
do 25 i=iof2+1,n3
vh(i)=vh(i)-sp
ekin=ekin+vh(i)*vh(i)
25 continue
!write(6,*) 'total linear momentum in z direction ',sp
!write(6,*) 'velocity adjustment'
ts=tscale*ekin
!write(6,*) 'temperature before scaling is ',ts
sc=tref/ts
sc=sqrt(sc)
!write(6,*) 'scale factor is ',sc
sc=sc*h
do 30 i=1,n3
vh(i)=vh(i)*sc
30 continue
end
c
c===================================================================
c
c random function generate random number between 0 and 1
c
c calling parameters are as follows
c
c vh vector of length npart (must be a multiple of 2)
c n3 vector length
c h scaling factor with which all elements of vh are multiplied
c tref temperature
c
c version 1.0 as of august 1985
c dieter w.heermann
c
c===================================================================
c
real function ranf(iseed)
double precision s,u,v
iseed=dble(iseed)
s=65536.0
u=2053.0
v=13849.0
m=iseed/s
iseed=iseed-m*s
iseed=u*iseed+v
m=iseed/s
iseed=iseed-m*s
ranf=iseed/s
return
end