c *************************************************************************** c Copyright(C) M.Kotelyanskii, November, 1997 c This code is just an example, it can be freely distributed and modified c as soon as this notice remains in the file. c It has been tested on the IBM Risc6000 under AIX. Other systems may require c changes, particularly the functions calls for the random numbers generation. c Even though it works properly to the best of the author's knowledge, c he does not assume any responsibility for any damage caused by it c *************************************************************************** c *************************************************************************** c ***** This is a code example to analyze the MD trajectory c ***** produced by the MD program prog7.f c ***** It reads the coordinates saved every Nprint steps from the file c ***** named "xyz.file" c ***** and calculates the mean sqd. displacement (MSD) vs. time c ***** It uses THE SAME parameter statments as prog7.f specifying c ***** Npart, dt, Nprint, etc. c ***** If you change them there DO NOT FORGET TO change them here too. c ***** This code is not optimized for performance, rather to illustrate c ***** the algorithm. c *************************************************************************** implicit real*8(a-h,o-z) c **** Specify parameters: c **** Number of trials parameter(Ncycle=10000, Nshots = 100, Nprint=Ncycle/Nshots) c **** Number of particles parameter(Npart=100) c **** Time step parameter(dt = 0.01, dt2=dt/2, dtsq2 = dt2*dt) c **** Set up arrays c **** Coordinates dimension x(2,Nshots,Npart) c **** MSD and counters to average dimension amsd (Nshots), amsd2 (Nshots), counter(Nshots) do i=1,Nshots amsd(i)=0.0d0 amsd2(i)=0.0d0 counter(i)=0.0d0 enddo c *** open trajectory file and read the snapshots open(unit=7,file="xyz.file",status="old") do ishot = 1, Nshots c *** read the snapshot (each ends with the empty record) read(7,11) time 11 format(e12.5) do ipart = 1, Npart read(7,12) x(1,ishot,ipart), x(2,ishot,ipart) 12 format(2(1x,e12.5)) enddo read(7,*) enddo close(7) c *** Calculate the MSD do ipart = 1,Npart do ishot = 1,Nshots-1 do jshot = ishot+1, Nshots idelta = jshot - ishot dx2= > (x(1,jshot,ipart)-x(1,ishot,ipart))**2 + > (x(2,jshot,ipart)-x(2,ishot,ipart))**2 amsd(idelta) = amsd(idelta) + dx2 amsd2(idelta) = amsd2(idelta) + dx2**2 counter(idelta) = counter(idelta) + 1 enddo enddo enddo c *** printout the result do ishot = 1, Nshots amean = amsd(ishot)/counter(ishot) amean2 = amsd2(ishot)/counter(ishot) error = sqrt( ( amean2 - amean**2)/counter(ishot) ) write(*,*) dt*Nprint*ishot, amean, error , counter(ishot) enddo c *** That's all folks end