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



