!////////////////////////////////////////////////////////////// ! ME 8253 -- COMPUTATIONAL NANOMECHANICS ! Spring 06 ! This program computes MD simulations upto 64 Argon atoms ! using the velocity Verlet algorithm !////////////////////////////////////////////////////////////// Program hw3 Implicit None Integer*4 L, Leq,M,N,Q,I,J,Nbin,Nop,Noiter,True,NbinMax Integer*4, Dimension(:),Allocatable ::B Real*8 Deltat, Energk, EnergPot, Energ, EnergIni, Lx, Ly, Lz Real*8 T, pi, Temp, EneFactor, DenBulk, distij, dr, CorrFactor Real*8 Deltax, Deltay, Deltaz, aDeltax, aDeltay, aDeltaz Real*8, Dimension(128)::Acx,Acy,Acz Real*8, Dimension(:), Allocatable::KinTot,Xc,Yc,Zc,Vex,Vey,Vez Real*8, Dimension(:), Allocatable::g, KinAve Integer :: ipart logical :: restart !C Number of particles Nop !C Periodic box dimensions Lx Ly Lz !C Time step deltat !C Number of iterations Noiter Nop=64 Lx=8.0d0 Ly=8.0d0 Lz=8.0d0 Deltat=1.0d-3 Noiter=15000 restart = .false. Q=5 !C Array dimensions Allocate (Xc(Nop)) Allocate (Yc(Nop)) Allocate (Zc(Nop)) Allocate(Vex(Nop)) Allocate(Vey(Nop)) Allocate(Vez(Nop)) Allocate(KinTot(Noiter+1)) Allocate(KinAve(Noiter+1)) !C Initial conditions If(Nop==2) Then Xc(1)=7.0 Yc(1)=6.0 Zc(1)=0.0 Vex(1)=0.0 Vey(1)=0.0 Vez(1)=0.0 Xc(2)=3.0 Yc(2)=6.0 Zc(2)=0.0 Vex(2)=0.4 Vey(2)=0.0 Vez(2)=0.0 Else If(Nop==64) Then if(restart) then Open(1,File='restart.dat') else !!$ Open(1,File='64IniData') Open(1,File='coords.dat') endif Do L=1, Nop if(restart) then Read(1,*) Xc(L), Yc(L), Zc(L), Vex(L), Vey(L), Vez(L) else Read(1,*) Xc(L), Yc(L), Zc(L) !!$ Cold Start Vex(L)=0.0 Vey(L)=0.0 Vez(L)=0.0 endif End Do !!$ Vex(3)=0.01 !!$ Vey(3)=0.02 !!$ Vez(3)=0.03 End If Call accel(Nop,Xc,Yc,Zc,Acx,Acy,Acz,Lx,Ly,Lz,EnergPot) Print *, 'Initial Potential=', EnergPot Energk=0.0d0 Do I=1,Nop Energk=Energk+ (Vex(I)**2 + Vey(I)**2 + Vez(I)**2) End Do Energk=0.5*Energk Print*, 'Initial Kinetic=',Energk EnergIni = EnergPot + Energk !C Output Files Open(2,File='64Coords1') Open(3,File='64Coords2') Open(4,File='64Kinetic') Open(5,File='64Potential') Open(6,File='64Total') Open(61,File='64EnergyAccuracy') Open(7,File='64KineticAve') Open(8,File='64Correlation') Open(9,File='64Equilibrium') Open(10,file='trajectory.pdb') Open(11,file='cm.dat') !C Modelling M=50 N=Q T=0.0 KinTot(1)=0. KinAve=0. !C Correlation Function Parameters dr=0.05 NbinMax=Int(Sqrt( (0.5*Lx)**2 + (0.5*Ly)**2 + (0.5*Lz)**2 ) / dr ) + 1 Allocate(B(NbinMax)) Allocate(g(NbinMax)) Do I=1, NbinMax B(I)=0. g(I)=0. End Do DenBulk=Nop/(Lx*Ly*Lz) pi=4.*Atan(1.) Do L=1, Noiter If(Nop==64 .and. M==50) Then Open(11,file='restart.dat') write(10,'(A4,I9)') 'MODEL',1 Do ipart=1,64 write(10,'(A4,I7,1x,A4,I4,9x,3f8.3)') 'ATOM',1,'C', & 0,Xc(ipart), Yc(ipart), Zc(ipart) write(11,'(6f10.4)') Xc(ipart), Yc(ipart), Zc(ipart) , Vex(ipart), Vey(ipart), Vez(ipart) Enddo Close(11) write(10,'(a6)') 'ENDMDL' M=0 End If !C Coordinates for 2 particle case every M steps If(Nop==2 .and. M==50) Then write(11,*) abs (( 0.5d0*(Xc(1)+Xc(2)) - 5.0d0 )/ 5.0d0) * 100.0d0 Write(2,*) T,Xc(1) Write(3,*) T,Xc(2) write(10,'(A4,I9)') 'MODEL',1 write(10,'(A4,I7,1x,A4,I4,9x,3f8.3)') 'ATOM',1,'C', & 0,Xc(1), Yc(1), Zc(1) write(10,'(A4,I7,1x,A4,I4,9x,3f8.3)') 'ATOM',2,'C', & 0,Xc(2), Yc(2), Zc(2) write(10,'(a6)') 'ENDMDL' M=0 End If !C Computing the pair-wise distance Do I=1, (Nop-1) Do J=I+1, Nop Deltax=Xc(I)-Xc(J) Deltay=Yc(I)-Yc(J) Deltaz=Zc(I)-Zc(J) !C Minimum image distance aDeltax=Abs(Deltax) aDeltay=Abs(Deltay) aDeltaz=Abs(Deltaz) If(aDeltax>0.5*Lx) Deltax=Deltax-Sign(Lx,Deltax) If(aDeltay>0.5*Ly) Deltay=Deltay-Sign(Ly,Deltay) If(aDeltaz>0.5*Lz) Deltaz=Deltaz-Sign(Lz,Deltaz) Distij=sqrt(Deltax**2 + Deltay**2 + Deltaz**2) !C Bin number Nbin=int(Distij/dr)+1 If(Nbin.Lt.NbinMax) Then B(Nbin)=B(Nbin)+1 endif End Do End Do !C Kin, Pot, Tot Energy and Ave KE every N steps Energ=Energk+EnergPot KinAve(L)=KinTot(L)/L !C If(N==Q) Then if( mod(L,10) ==0 ) then Write(4,*) T, Energk Write(5,*) T, EnergPot Write(6,*) T, Energ endif Write(61,*) T, 100*(Energ/EnergIni-1.0) !C N=0 !C End If Write(7,*) T, KinTot(L)/L M=M+1 N=N+1 Call update(Nop,Xc,Yc,Zc,Vex,Vey,Vez,Acx,Acy,Acz,& Lx,Ly,Lz,Deltat,EnergPot,Energk) KinTot(L+1)=KinTot(L)+ Energk T=T+Deltat End Do !C Equilibrium Time and Temperature If(Nop /= 2) Then True=1 Do L=5000, Noiter-4 If(True==1 .and. & Abs(KinAve(L+5)-KinAve(L))< ((1.e-8)*KinAve(L))) Then Write(9,*) 'Iterations to reach equilibrium=',L Leq=L Write(9,*) 'Equilibrium time=', Leq*Deltat True=0 End If End Do If(True==1) Then Write(9,*)'Equilibrium was not reached' Leq=Noiter End If !C Dimensionless temperature for argon with E/KB= 119.8 EneFactor=119.8 Temp=KinAve(Noiter)/Nop Temp=2.*Temp/3. Write(9,*)'Temperature=',Temp Temp=EneFactor*Temp Write(9,*)'TemperatureEne=',Temp !C Pair correlation function CorrFactor=1./(4.*pi*Nop*DenBulk*(dr)**3) print*,NbinMax Do I=1, NbinMax B(I)=B(I)/Noiter g(I)=CorrFactor*B(I)/I**2 Write(8,*) I*dr, 2*g(I) End Do End If Contains !C Acceleration subroutine Subroutine accel(np,x,y,z,ax,ay,az,sx,sy,sz,pot) Integer*4 i,j,np Real*8 x(128),y(128),z(128),ax(128),ay(128),az(128) Real*8 dx,dy,dz,adx,ady,adz Real*8 ri,ri2,ri4,ri6,ri8,b,sx,sy,sz,pot Do i=1, np ax(i)=0.0 ay(i)=0.0 az(i)=0.0 End Do pot =0.0 !C Compute pair-wise force Do i=1,(np-1) Do j=i+1, np dx=x(i)-x(j) dy=y(i)-y(j) dz=z(i)-z(j) !C minimum image distance adx=Abs(dx) ady=Abs(dy) adz=Abs(dz) !!$ If(adx>0.5*sx) adx=adx-Sign(sx,dx) !!$ If(ady>0.5*sy) ady=ady-Sign(sy,dy) !!$ If(adz>0.5*sz) adz=adz-Sign(sz,dz) If(adx>0.5*sx) dx=dx-Sign(sx,dx) If(ady>0.5*sy) dy=dy-Sign(sy,dy) If(adz>0.5*sz) dz=dz-Sign(sz,dz) !!$ print*,sx,dx,Sign(sx,dx) !!$ If(dx.Ge.0.5*sx) Then !!$ dx = dx - sx !!$ Elseif(dx.Lt.-0.5*sx) Then !!$ dx = dx + sx !!$ Endif !!$ !!$ If(dy.Ge.0.5*sy) Then !!$ dy = dy - sy !!$ Elseif(dy.Lt.-0.5*sy) Then !!$ dy = dy + sy !!$ Endif !!$ !!$ If(dz.Ge.0.5*sz) Then !!$ dz = dz - sz !!$ Elseif(dz.Lt.-0.5*sz) Then !!$ dz = dz + sz !!$ Endif ri=1./sqrt(dx*dx+dy*dy+dz*dz) ri2=ri*ri ri4=ri2*ri2 ri6=ri4*ri2 ri8=ri6*ri2 pot=pot+ (ri6 -1.)*ri6 b=24.*(2.*ri6-1.)*ri8 ax(i)=ax(i) + b*dx ay(i)=ay(i) + b*dy az(i)=az(i) + b*dz !!$ ax(j)=ax(j)+b*dx !!$ ay(j)=ay(j)+b*dy !!$ az(j)=az(j)+b*dz ax(j)=ax(j) - b*dx ay(j)=ay(j) - b*dy az(j)=az(j) - b*dz End Do End Do pot=4.*pot Return End Subroutine accel !C Update subroutine Subroutine update(np,x,y,z,vx,vy,vz,ax,ay,az,sx,sy,sz,dt,pot,ek) Integer*4 i,np Real*8 x(128),y(128),z(128),vx(128),vy(128),vz(128) Real*8 ax(128),ay(128),az(128),dt,hdt2,hdt,sx,sy,sz,pot,ek ek=0. hdt=0.5*dt hdt2=hdt*dt Do i=1,np x(i)=x(i)+dt*vx(i) + hdt2*ax(i) y(i)=y(i)+dt*vy(i) + hdt2*ay(i) z(i)=z(i)+dt*vz(i) + hdt2*az(i) vx(i)=vx(i)+hdt*ax(i) vy(i)=vy(i)+hdt*ay(i) vz(i)=vz(i)+hdt*az(i) End Do Call accel(np,x,y,z,ax,ay,az,sx,sy,sz,pot) !C enforce PBC Do i=1,np x(i)=Mod(x(i)+sx,sx) y(i)=Mod(y(i)+sy,sy) z(i)=Mod(z(i)+sz,sz) End Do Do i=1,np vx(i)=vx(i)+hdt*ax(i) vy(i)=vy(i)+hdt*ay(i) vz(i)=vz(i)+hdt*az(i) ek=ek+(vx(i)**2+vy(i)**2+vz(i)**2) End Do ek=0.5*ek Return End subroutine update End Program hw3