PROGRAM md IMPLICIT NONE REAL*8 E, area, ke, kecum, pe, pecum, t, vcum, + virial, Lx, Ly, R2cum, ax, ay, dr, dt, dt2, gcum, + nbin, vx, vy, x, xsave, y, ysave c '"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"' REAL*8 vxi, vyi, xi, yi INTEGER ooxx ! limit cycle to 3.5 loops c ._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ INTEGER IC1_, IC2_, IR_, N, ncum, ii, switch CHARACTER*80 flag LOGICAL Ltmp1_ COMMON /GLBLD/Lx, Ly, R2cum(100), ax(36), ay(36), dr, dt, dt2, + gcum(0:1000), nbin, vx(36), vy(36), x(36), + xsave(100), y(36), ysave(100) COMMON /GLBLI/N c '"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"' COMMON /xxoo/ vxi(16), vyi(16), xi(16), yi(16) c ._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ DATA IC1_/2/ DATA IC2_/1/ CALL initial(t,ke,kecum,pecum,vcum,area) CALL set_up_windows(IC1_,IC2_) c111 continue c CALL initial(t,ke,kecum,pecum,vcum,area) c CALL set_up_windows(IC1_,IC2_) c '"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"' 321 continue c ._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ CALL accel(pe,virial) E = ke + pe ! total energy ncum = 0 ! number of times data accumulated flag = '' Ltmp1_ = .TRUE. c DO WHILE(Ltmp1_) c do ii = 1, 8000 ii = 0 switch = -1 c '"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"' ooxx = 0 c ._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ 222 continue c if (ii.lt.4000) then c switch = -1*sign(sin(3.14159268*ii/4000),1.0) c else c switch = 1 c switch = -1 if (ii.gt.40000) then ii = ii - 40000 switch = switch*-1 end if c end if CALL show_positions(flag,IC2_) CALL Verlet(t,ke,pe,virial,switch) CALL show_output(t,ke,pe,virial,kecum,vcum,ncum,area,IC1_) c call check_stop (ncum,flag) Ltmp1_ = .NOT.(flag .EQ. 'stop') ii = ii + 1 c '"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"' ooxx = ooxx + 1 if (ooxx.LT.200000) goto 222 do ii = 1, N x(ii) = xi(ii) y(ii) = yi(ii) vx(ii) = vxi(ii) vy(ii) = vyi(ii) enddo t = 0.D0 goto 321 c ._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ c END DO c goto 111 CALL save_config c CALL GWquit(IR_) call pgclos END subroutine check_stop (nucm, flag) integer ncum character*80 flag if ( mod(ncum,1000) .eq. 0 ) then open (unit=21, file='md_stop.txt', status='old',err= 201) read (*,*) flag close (21) endif 201 continue end C mhl done SUBROUTINE initial(t,ke,kecum,pecum,vcum,area) IMPLICIT NONE REAL*8 t, ke, kecum, pecum, vcum, area, Lx, Ly, + R2cum, ax, ay, dr, dt, dt2, gcum, nbin, vx, vy, x, xsave, + y, ysave, DATA_(64) c '"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"' REAL*8 vxi, vyi, xi, yi c ._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ INTEGER IO1_, IR_, N, i CHARACTER*80 file, heading, response, start LOGICAL Ltmp1_ COMMON /GLBLD/Lx, Ly, R2cum(100), ax(36), ay(36), dr, dt, dt2, + gcum(0:1000), nbin, vx(36), vy(36), x(36), + xsave(100), y(36), ysave(100) c '"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"' COMMON /xxoo/ vxi(16), vyi(16), xi(16), yi(16) c ._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ COMMON /GLBLI/N INTEGER DP_ DATA DATA_/ + 1.09D0,0.98D0,-0.33D0,0.78D0,3.12D0,5.25D0,0.12D0,-1.19D0, + 0.08D0,2.38D0,-0.08D0,-0.10D0,0.54D0,4.08D0,-1.94D0,-0.56D0, + 2.52D0,4.39D0,0.75D0,0.34D0,3.03D0,2.94D0,1.70D0,-1.08D0, + 4.25D0,3.01D0,0.84D0,0.47D0,0.89D0,3.11D0,-1.04D0,0.06D0, + 2.76D0,0.31D0,1.64D0,1.36D0,3.14D0,1.91D0,0.38D0,-1.24D0, + 0.23D0,5.71D0,-1.58D0,0.55D0,1.91D0,2.46D0,-1.55D0,-0.16D0, + 4.77D0,0.96D0,-0.23D0,-0.83D0,5.10D0,4.63D0,-0.31D0,0.65D0, + 4.97D0,5.88D0,1.18D0,1.48D0,3.90D0,0.20D0,0.46D0,-0.51D0/ DATA DP_/1/, IO1_/11/ dt = 0.01D0 dt2 = dt*dt response = '' Ltmp1_ = .TRUE. DO WHILE(Ltmp1_) c WRITE(*,'(A,$)') 'read data statements (d) or file (f)? ' c READ(*,'(A)') start start = 'd' IF((start .EQ. 'd') .OR. (start .EQ. 'D')) THEN response = 'ok' N = 16 c mhl Lx = 6 c mhl Ly = 6 Lx = 12 Ly = 12 DO i = 1, N x(i) = DATA_(DP_) DP_ = DP_ + 1 y(i) = DATA_(DP_) DP_ = DP_ + 1 vx(i) = DATA_(DP_) c mhl lower the initial temperature c vx(i) = vx(i) DP_ = DP_ + 1 vy(i) = DATA_(DP_) c mhl lower the initial temperature c vy(i) = vy(i) DP_ = DP_ + 1 c '"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"'"' xi(i) = x(i) yi(i) = y(i) vxi(i) = vx(i) vyi(i) = vy(i) c ._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ END DO ELSE IF((start .EQ. 'f') .OR. (start .EQ. 'F')) THEN response = 'ok' WRITE(*,'(A,$)') 'file name = ' READ(*,'(A)') file OPEN(IO1_, file = file) READ(IO1_,*) N READ(IO1_,*) Lx READ(IO1_,*) Ly READ(IO1_,'(A)') heading DO i = 1, N READ(IO1_,*) x(i),y(i) END DO READ(IO1_,'(A)') heading DO i = 1, N READ(IO1_,*) vx(i),vy(i) c mhl lower the initial temperature c vx(i) = 0.001*vx(i) c vy(i) = 0.001*vy(i) END DO CLOSE(IO1_) ELSE WRITE(*,*) WRITE(*,*) 'd or f are the only acceptable responses.' END IF Ltmp1_ = response .NE. 'ok' END DO c mhl added do i=1,n xsave(i) = x(i) ysave(i) = y(i) end do ke = 0 ! kinetic energy DO i = 1, N ke = ke + vx(i)*vx(i) + vy(i)*vy(i) END DO ke = 0.5D0*ke area = Lx*Ly t = 0 ! time c initialize sums kecum = 0 pecum = 0 vcum = 0 END C mhl done SUBROUTINE set_up_windows(IC1_,IC2_) IMPLICIT NONE REAL*8 Lx, Ly, R2cum, ax, ay, dr, dt, dt2, gcum, + nbin, vx, vy, x, xsave, y, ysave REAL xwin, ywin, GWaspect INTEGER IR_, N, IC1_, IC2_, pgopen COMMON /GLBLD/Lx, Ly, R2cum(100), ax(36), ay(36), dr, dt, dt2, + gcum(0:1000), nbin, vx(36), vy(36), x(36), + xsave(100), y(36), ysave(100) COMMON /GLBLI/N if( pgopen('/xwin') .le. 0) stop 'Open device fail.' call pgpap (10.0, 0.5) call pgiden call pgstbg(0) call pgsubp (2,1) call pgenv (0.0, real(Lx), 0.0, real(Ly), 0, -1) call pgsch (1.5) call pglab ('', '', 'Program for MD siulation.') call pgpanl (2,1) call headings END C mhl done SUBROUTINE headings IMPLICIT NONE CHARACTER*80 Stmp1_ INTEGER IR_ c call pgsch (2.0) call pgsch (1.5) WRITE(Stmp1_, '(A)') 'time steps ' call pgptxt(2.5, 4.5, 0.0, 1.0, 'time steps = ') WRITE(Stmp1_, '(A)') 'time ' call pgptxt(2.5, 3.5, 0.0, 1.0, 'time = ') WRITE(Stmp1_, '(A)') 'energy ' call pgptxt(2.5, 2.5, 0.0, 1.0, 'energy = ') WRITE(Stmp1_, '(A)') ' ' call pgptxt(2.5, 1.5, 0.0, 1.0, ' = ') WRITE(Stmp1_, '(A)') '

' call pgptxt(2.5, 0.5, 0.0, 1.0, '

= ') END C SUBROUTINE Verlet(t,ke,pe,virial,switch) IMPLICIT NONE REAL*8 xnew, ynew, t, ke, pe, virial, Lx, Ly, R2cum, ax, ay, + dr, dt, dt2, gcum, nbin, vx, vy, x, xsave, y, ysave, pbc, + scale_factor INTEGER N, i ,switch COMMON /GLBLD/Lx, Ly, R2cum(100), ax(36), ay(36), dr, dt, dt2, + gcum(0:1000), nbin, vx(36), vy(36), x(36), + xsave(100), y(36), ysave(100) COMMON /GLBLI/N scale_factor = 1.0 + switch*0.0001 DO i = 1, N xnew = x(i) + vx(i)*dt + 0.5D0*ax(i)*dt2 ynew = y(i) + vy(i)*dt + 0.5D0*ay(i)*dt2 c partially update velocity using old acceleration vx(i) = vx(i) + 0.5D0*ax(i)*dt vy(i) = vy(i) + 0.5D0*ay(i)*dt c mhl gradually reduce or increase temperature vx(i) = scale_factor*vx(i) vy(i) = scale_factor*vy(i) c mhl x(i) = pbc(xnew,Lx) y(i) = pbc(ynew,Ly) END DO CALL accel(pe,virial) ! new acceleration ke = 0 DO i = 1, N c complete the update of the velocity using new acceleration vx(i) = vx(i) + 0.5D0*ax(i)*dt vy(i) = vy(i) + 0.5D0*ay(i)*dt ke = ke + vx(i)*vx(i) + vy(i)*vy(i) END DO ke = 0.5D0*ke t = t + dt END C FUNCTION pbc(pos,L) IMPLICIT NONE REAL*8 pbc, pos, L IF(pos .LT. 0) THEN pbc = pos + L ELSE IF(pos .GT. L) THEN pbc = pos - L ELSE pbc = pos END IF END C SUBROUTINE accel(pe,virial) IMPLICIT NONE REAL*8 dx, dy, fxij, fyij, pot, pe, virial, Lx, Ly, R2cum, ax, + ay, dr, dt, dt2, gcum, nbin, vx, vy, x, xsave, y, ysave, + separation INTEGER N, i, j COMMON /GLBLD/Lx, Ly, R2cum(100), ax(36), ay(36), dr, dt, dt2, + gcum(0:1000), nbin, vx(36), vy(36), x(36), + xsave(100), y(36), ysave(100) COMMON /GLBLI/N DO i = 1, N ax(i) = 0 ay(i) = 0 END DO pe = 0 virial = 0 DO i = 1, N - 1 DO j = i + 1, N dx = separation(x(i) - x(j),Lx) dy = separation(y(i) - y(j),Ly) c acceleration = force because mass = 1 in reduced units CALL force(dx,dy,fxij,fyij,pot) ax(i) = ax(i) + fxij ay(i) = ay(i) + fyij ax(j) = ax(j) - fxij ! Newton's third law ay(j) = ay(j) - fyij pe = pe + pot virial = virial + dx*fxij + dy*fyij END DO END DO END C SUBROUTINE force(dx,dy,fx,fy,pot) IMPLICIT NONE REAL*8 f_over_r, r2, rm2, rm6, dx, dy, fx, fy, pot r2 = dx*dx + dy*dy rm2 = 1/r2 rm6 = rm2*rm2*rm2 f_over_r = 24*rm6*(2*rm6 - 1)*rm2 fx = f_over_r*dx fy = f_over_r*dy pot = 4*(rm6*rm6 - rm6) END C FUNCTION separation(ds,L) IMPLICIT NONE REAL*8 separation, ds, L IF(ds .GT. 0.5D0*L) THEN separation = ds - L ELSE IF(ds .LT. -0.5D0*L) THEN separation = ds + L ELSE separation = ds END IF END C SUBROUTINE show_positions(flag,IC2_) IMPLICIT NONE REAL*8 Lx, Ly, R2cum, ax, ay, dr, dt, dt2, + gcum, nbin, vx, vy, x, xsave, y, ysave real radius parameter (radius = 0.10) INTEGER IR_, N, i, k, IC2_ LOGICAL TBkeyinput CHARACTER*80 flag COMMON /GLBLD/Lx, Ly, R2cum(100), ax(36), ay(36), dr, dt, dt2, + gcum(0:1000), nbin, vx(36), vy(36), x(36), + xsave(100), y(36), ysave(100) COMMON /GLBLI/N call pgpanl (1,1) call pgsci (1) c call pgbbuf do i = 1, n call pgsci (0) call pgcirc (real(xsave(i)), real(ysave(i)), radius) call pgsci (1) call pgbox ('bc',0.0, 0,'bc',0.0 ,0) call pgsci (mod(i,4)+1) call pgcirc ( real(x(i)), real(y(i)), radius ) end do c call pgebuf do i = 1, n xsave(i) = x(i) ysave(i) = y(i) end do END C mhl done SUBROUTINE show_output(t,ke,pe,virial,kecum,vcum,ncum,area,IC1_) IMPLICIT NONE CHARACTER*80 Stmp1_ REAL*8 E, mean_ke, p, t, ke, pe, virial, kecum, vcum, area, Lx, + Ly, R2cum, ax, ay, dr, dt, dt2, gcum, nbin, vx, vy, x, + xsave, y, ysave INTEGER IR_, N, ncum, IC1_ COMMON /GLBLD/Lx, Ly, R2cum(100), ax(36), ay(36), dr, dt, dt2, + gcum(0:1000), nbin, vx(36), vy(36), x(36), + xsave(100), y(36), ysave(100) COMMON /GLBLI/N call pgpanl (2,1) c call pgsch (2.0) call pgsch (1.5) ncum = ncum + 1 c call pgbbuf WRITE(Stmp1_, '(I10)') ncum c call pgtext (1.5,4.5,'time steps = '//Stmp1_) call pgptxt (2.5,4.5, 0.0, 0.0, Stmp1_) WRITE(Stmp1_, '(F6.2)') t ! time c call pgtext (1.5,3.5,'time = '//Stmp1_) call pgptxt (2.5,3.5, 0.0, 0.0, Stmp1_) E = ke + pe ! total energy WRITE(Stmp1_, '(F8.3)') E c call pgtext (1.5,2.5,' ='//Stmp1_) call pgptxt (2.5,2.5, 0.0, 0.0, Stmp1_) kecum = kecum + ke vcum = vcum + virial mean_ke = kecum/ncum ! still need to divide by N p = mean_ke + (0.5D0*vcum)/ncum ! mean pressure * area p = p/area WRITE(Stmp1_, '(F5.2)') mean_ke/N ! mean kinetic temperature c call pgtext (1.5,1.5,' ='//Stmp1_) call pgptxt (2.5,1.5, 0.0, 0.0, Stmp1_) WRITE(Stmp1_, '(F5.2,$)') p c call pgtext (1.5,0.5,'Pressure ='//Stmp1_) call pgptxt (2.5,0.5, 0.0, 0.0, Stmp1_) c call pgebuf END C SUBROUTINE save_config IMPLICIT NONE REAL*8 Lx, Ly, R2cum, ax, ay, dr, dt, dt2, gcum, nbin, vx, vy, + x, xsave, y, ysave INTEGER IR_, IO1_, N, i CHARACTER*80 filename COMMON /GLBLD/Lx, Ly, R2cum(100), ax(36), ay(36), dr, dt, dt2, + gcum(0:1000), nbin, vx(36), vy(36), x(36), + xsave(100), y(36), ysave(100) COMMON /GLBLI/N DATA IO1_/11/ filename = 'md_save.dat' OPEN(IO1_, file = filename) WRITE(IO1_,*) N WRITE(IO1_,*) Lx WRITE(IO1_,*) Ly WRITE(IO1_,'(2A12)') 'x(i)', 'y(i)' * comma added between outputs on the same line so that file * can be read by True BASIC DO i = 1, N WRITE(IO1_, '(3X,F9.4,'', '',F9.5)') x(i), y(i) END DO WRITE(IO1_, '(3X,2A9)') 'vx(i)', 'vy(i)' DO i = 1, N WRITE(IO1_, '(3X,F9.4,'', '',F9.5)') vx(i), vy(i) END DO CLOSE(IO1_) END