      program mc_nvt
c                                                    
c   equation of state of the lennard-jones fluid     
c                                                    

c   data files:
c	fort.10 - input  - job parameters
c	fort.11 - input  - initial atomic positions (if available)
c	fort.21 - output - interim atomic positions during run, and 
c                          final positons at the end 
c	fort.22 - output - sequence of atomic positions in pdb format
c	fort.30 - output - instantaneous energy and pressure

      implicit none
      include 'params.inc'
      integer equil, prod, nsamp, ii, icycl, ndispl, attempt, 
     &        nacc, ncycl, nmoves, imove
      real en, ent, vir, virt, dr, av1, av2, press
 
      write (6, *) 
      write (6, *) '**************** mc_nvt ***************'

c     ---initialize system

      call readdat(equil, prod, nsamp, ndispl, dr)

      nmoves = ndispl
      av1    = 0.0
      av2    = 0.0

c     ---total energy of the system

      call toterg(en, vir)
      write (6, 99001) en, vir

c     ---start mc-cycle

      do ii = 1, 2

c        --- ii=1 equilibration
c        --- ii=2 production

         if (ii.eq.1) then
            ncycl = equil
            if (ncycl.ne.0) write (6, *) ' start equilibration '
         else
            if (ncycl.ne.0) write (6, *) ' start production '
            ncycl = prod
         end if
        
         attempt = 0
         nacc    = 0

c        ---intialize the subroutine that adjust the maximum displacement

         call adjust(attempt, nacc, dr)

         do icycl = 1, ncycl

            do imove = 1, nmoves

c              ---attempt to displace a particle

               call mcmove(en, vir, attempt, nacc, dr)
            end do

            if (ii.eq.2) then

c              ---instantaneous sample averages

               if (mod(icycl,nsamp).eq.0) then
                  call sample(icycl, en, vir, press)
                  av1 = av1 + press
                  av2 = av2 + 1.0

               endif
            endif

            if (mod(icycl,ncycl/5).eq.0) then
               write (6, *) '======>> done ', icycl, ' out of ', ncycl

c              ---write intermediate configuration to file

               call store(21, dr)

c              ---adjust maximum displacements

               call adjust(attempt, nacc, dr)
            end if
         end do
         if (ncycl.ne.0) then
            if (attempt.ne.0) write (6, 99003) attempt, nacc, 
     &                               100.0*real(nacc)/real(attempt)

c           ---test total energy

            call toterg(ent, virt)
            if (abs(ent-en).gt.1.e-6) then
               write (6, *) 
     &                    ' ######### problems energy ################ '
            end if
            if (abs(virt-vir).gt.1.e-6) then
               write (6, *) 
     &              ' ######### problems virial ################ '
            end if
            write (6, 99002) ent, en, ent - en, virt, vir, virt - vir
            write (6,*)

c     print chemical potential and pressure   

            if(ii.eq.2) then
               write (6,*) 'average pressure                  : ',
     &              av1/av2
            end if
         end if
      end do
     
c     write final state in pdb format

               call writepdb(22)

      stop
 
99001 format (' total energy initial configuration: ', f12.5, /, 
     &        ' total virial initial configuration: ', f12.5,/)
99002 format (/,' total energy end of simulation    : ', f12.5, /, 
     &        '       running energy              : ', f12.5, /, 
     &        '       difference                  :  ', e12.5, /, 
     &        ' total virial end of simulation    : ', f12.5, /, 
     &        '       running virial              : ', f12.5, /, 
     &        '       difference                  :  ', e12.5)
99003 format (' number of att. to displ. a part.  : ', i10, /, 
     &        ' success: ', i10, '(= ', f5.2, '%)')
      end



      subroutine adjust(attemp, nacc, dr)
c
c     adjusts maximum displacement such that 50% of the
c     movels will be accepted
c
c  attemp (input)  number of attemps that have been performed 
c                  to displace a particle
c  nacc   (input)  number of successful attemps to 
c                  displace a particle
c  dr     (output) new maximum displacement
c
      implicit none
      include 'params.inc'
      integer attemp, nacc, attempp, naccp
      real dro, frac, dr
      save naccp, attempp
      data naccp/0/
      data attempp/0/
 
      if (attemp.eq.0.or.attempp.ge.attemp) then
         naccp = nacc
         attempp = attemp
      else
         frac = real(nacc-naccp)/real(attemp-attempp)
         dro  = dr
         dr   = dr*abs(frac/0.5)

c        ---limit the change:

         if (dr/dro.gt.1.5) dr = dro*1.5
         if (dr/dro.lt.0.5) dr = dro*0.5
         if (dr.gt.hbox/2.) dr = hbox/2.
         write (6, 99001) dr, dro, frac, attemp - attempp, nacc - naccp

c        ---store nacc and attemp for next use

         naccp = nacc
         attempp = attemp
      end if
      return
99001 format ('max. displ. set to : ', f6.3, ' ( old : ', f6.3, ')', /, 
     &        'frac. acc.: ', f4.2, ' attempts: ', i7, ' success: ', i7)
      end



      subroutine ener(en, vir, r2)
c
c     calculate energy between a (single) pair of atoms
c
c     en : (output) energy
c     vir: (output) virial
c     r2 : (input) distance squared between two particles
c
      implicit none
      real r2, r2i, r6i, en, vir
      include 'params.inc'
 
      en = 0.0
      vir = 0.0

      if (r2.lt.rc2) then
         r2i = sig2/r2
         r6i = r2i*r2i*r2i
         en  = eps4*(r6i*r6i-r6i)
         vir = eps48*(r6i*r6i-0.5*r6i)

      end if
      
      return
      end

 
      subroutine eneri(xi, yi, zi, i, jb, en, vir)
c
c    calculates the energy of particle i with particles j=jb,npart
c
c     xi (input)    x coordinate particle i
c     yi (input)    y coordinate particle i
c     zi (input)    z coordinate particle i
c     i  (input)    particle number (excluded !!!)
c     en  (output)  energy particle i
c     vir (output)  virial particle i
c
c
      implicit none
      include 'params.inc'
      real xi, yi, zi, en, dx, dy, dz, 
     &     r2, vir, virij, enij
      integer i, j, jb

      en  = 0.0
      vir = 0.0

      do j = jb, npart

ccccccccccccccccccccccccccc
c     excluse particle i  c
ccccccccccccccccccccccccccc

         if (j.ne.i) then

            dx = xi - x(j)
            dy = yi - y(j)
            dz = zi - z(j)

ccccccccccccccccccccccccccccccccccc
c     nearest image convention    c
ccccccccccccccccccccccccccccccccccc

            if (dx.gt.hbox) then
               dx = dx - box
            elseif (dx.lt.-hbox) then
               dx = dx + box
            end if

            if (dy.gt.hbox) then
               dy = dy - box
            elseif (dy.lt.-hbox) then
               dy = dy + box
            end if

            if (dz.gt.hbox) then
               dz = dz - box
            elseif (dz.lt.-hbox) then
               dz = dz + box
            end if

            r2 = dx*dx + dy*dy + dz*dz

cccccccccccccccccccccccccccccc
c     calculate the energy   c
cccccccccccccccccccccccccccccc

            call ener(enij, virij, r2)

            en  = en + enij
            vir = vir + virij
         endif
      enddo

      return
      end



      subroutine lattice
c
c     place `npart' particles on a simple cubic
c     lattice with density 'rho'
c
      implicit none
      include 'params.inc'
      integer i, j, k, itel, n
      real dx, dy, dz, del
 
      n = int(npart**(1.0/3.0)) + 1
      if (n.eq.0) n = 1
      del = box/real(n)
      itel = 0
      dx = -del
      do i = 1, n
         dx = dx + del
         dy = -del
         do j = 1, n
            dy = dy + del
            dz = -del
            do k = 1, n
               dz = dz + del
               if (itel.lt.npart) then
                  itel = itel + 1
                  x(itel) = dx
                  y(itel) = dy
                  z(itel) = dz
               end if
            end do
         end do
      end do
      write (6, 99001) itel
      return
99001 format (' initialization on lattice: ', /, i10, 
     &        ' particles placed on a lattice')
      end



      subroutine mcmove(en, vir, attempt, nacc, dr)
c
c     attempts to displace a randomly selected particle
c
c
c     ener   (input/output) : total energy
c     vir    (input/output) : total virial
c     attemp (input/output) : number of attemps that have been
c                             performed to displace a particle
c     nacc   (input/output) : number of successful attemps
c                             to displace a particle
c     dr     (input)        : maximum displacement
c
      implicit none
      include 'params.inc'
      real enn, eno, en, ran_uniform, xn, yn, 
     &     zn, viro, virn, vir, dr
      integer o, attempt, nacc, jb
 
      attempt = attempt + 1
      jb = 1

c     ---select a particle at random

      o = int(real(npart)*ran_uniform()) + 1

c     ---calculate energy old configuration

      call eneri(x(o), y(o), z(o), o, jb, eno, viro)

c     ---give particle a random displacement

      xn = x(o) + (ran_uniform()-0.5)*dr 
      yn = y(o) + (ran_uniform()-0.5)*dr
      zn = z(o) + (ran_uniform()-0.5)*dr

c     ---calculate energy new configuration:

      call eneri(xn, yn, zn, o, jb, enn, virn)

c     ---acceptance test

      if (ran_uniform().lt.exp(-beta*(enn-eno))) then

c        --accepted

         nacc = nacc + 1
         en = en + (enn-eno)
         vir = vir + (virn-viro)

c        ---put particle in simulation box

         if (xn.lt.0.0) then
            xn = xn + box
         elseif (xn.gt.box) then
            xn = xn - box
         endif

         if (yn.lt.0.0) then
            yn = yn + box
         elseif (yn.gt.box) then
            yn = yn - box
         endif

         if (zn.lt.0.0) then
            zn = zn + box
         elseif (zn.gt.box) then
            zn = zn - box
         endif

         x(o) = xn
         y(o) = yn
         z(o) = zn
      end if
      return
      end



      subroutine genrand(i)
      implicit none
      
      integer i,j
      real x,ran_uniform

      do j=1,i 
         x = ran_uniform()
      enddo  
      
      return 
      end




      function ran_uniform()
      implicit none

      integer idum
      real ran_uniform,ran2
      save idum
      data idum/-10/

      ran_uniform = ran2(idum)

      return
      end



      function ran2(idum)
      implicit none
      integer idum,im1,im2,imm1,ia1,ia2,
     &     iq1,iq2,ir1,ir2,ntab,ndiv
      real ran2,am,eps,rnmx
      parameter (im1=2147483563,im2=2147483399,
     &     am=1.0/im1,imm1=im1-1,ia1=40014,
     &     ia2=40692,iq1=53668,iq2=52774,ir1=12211,
     &     ir2=3791,ntab=32,ndiv=1+imm1/ntab,
     &     eps=1.2e-7,rnmx=1.0-eps)
      integer idum2,j,k,iv(ntab),iy
      save iv,iy,idum2
      data idum2/123456789/, iv/ntab*0/, iy/0/

      if (idum.le.0) then
         idum=max(-idum,1)
         idum2=idum
         do j=ntab+8,1,-1
            k=idum/iq1
            idum=ia1*(idum-k*iq1)-k*ir1
            if (idum.lt.0) idum=idum+im1
            if (j.le.ntab) iv(j)=idum
         enddo
         iy=iv(1)
      endif

      k=idum/iq1
      idum=ia1*(idum-k*iq1)-k*ir1
      if (idum.lt.0) idum=idum+im1
      k=idum2/iq2
      idum2=ia2*(idum2-k*iq2)-k*ir2
      if (idum2.lt.0) idum2=idum2+im2
      j=1+iy/ndiv
      iy=iv(j)-idum2
      iv(j)=idum
      if(iy.lt.1)iy=iy+imm1
      ran2=min(am*iy,rnmx)

      return
      end



      subroutine readdat(equil, prod, nsamp, ndispl, dr)
      implicit none

c     ---read input data and model parameters
c
c     ---input parameters: file: fort.10
c    ibeg  =  0 : initialize from a lattice
c             1 : read configuration from disk
c    equil      : number of monte carlo cycles during equilibration
c    prod       : number of monte carlo cycles during production
c    nsamp      : number of monte carlo cycles between two sampling periods
c    dr         : maximum displacement
c    ndispl     : number of attemps to displace a particle per mc cycle
c    npart      : total numbero fo particles
c    temp       : temperature
c    rho        : density
c    istart     : number of initial ran2 calls
c    eps        : epsilon in lennard-jones potential
c    sig        : sigma in lennard-jones potential
c    mass       : mass of the particle
c    rc         : cut-off radius of the potential
c
c     ---input parameters: file: fort.11 (restart file
c                to continue a simulation from disk)
c    boxf   = box length old configuration (if this one
c             does not correspond to the requested density, the positions
c             of the particles are rescaled!
c    npart  = number of particles (over rules fort.10!!)
c    dr     = optimized maximum displacement old configurations
c
c
c    x(1),y(1),z(1)            : position first particle 1
c        ...
c    x(npart),y(npart),z(npart): position particle last particle
 
      include 'params.inc'

      integer ibeg, equil, prod, i, ndispl, nsamp, istart
      real eps, sig, boxf, rhof, rho, dr
 
c     ---read simulation data

      read (10, *) ibeg, equil, prod, nsamp
      read (10, *) dr, ndispl
      read (10, *) npart, temp, rho
      read (10, *) istart

c     ---initialize random number generator for istart calls
      
      call genrand(istart)
 
      if (npart.gt.npmax) then
         write (6, *) ' error: number of particles too large'
         stop
      end if

      box  = (real(npart)/rho)**(1./3.)
      hbox = box/2

c     ---read model parameters

      read (10, *) eps, sig, mass, rc

c     ---read/generate configuration

      if (ibeg.eq.0) then

c        ---generate configuration form lattice

         call lattice
      else
         write (6, *) ' read configuration from disk '
         read (11, *) npart, boxf, dr
         rhof = real(npart)/boxf**3
         if (abs(boxf-box).gt.1.e-6) then
            write (6, 99007) rho, rhof
         end if
         do i = 1, npart
            read (11, *) x(i), y(i), z(i)
            x(i) = x(i)*box/boxf
            y(i) = y(i)*box/boxf
            z(i) = z(i)*box/boxf
         end do
         rewind (11)
      end if

c     ---write input data

      write (6, 99001) equil, prod, nsamp
      write (6, 99002) ndispl, dr
      write (6, 99003) npart, temp, rho, box
      write (6, 99004) eps, sig, mass

c     ---calculate parameters:

      beta = 1/temp
      
c     ---calculate cut-off radius potential

      rc    = min(rc, hbox)
      rc2   = rc*rc
      eps4  = 4.0*eps
      eps48 = 48.0*eps
      sig2  = sig*sig
      
      return

99001 format ('  number of equilibration cycles             :', i10, /, 
     &        '  number of production cycles                :', i10, /, 
     &        '  sampling frequency                         :', i10, /)
99002 format ('  number of attempted particle moves / cycle :', i10, /, 
     &        '  maximum displacement                       :', f10.3, 
     &        /)
99003 format ('  number of particles                        :', i10, /, 
     &        '  temperature                                :', f10.3, 
     &        /, '  density                                    :', 
     &        f10.3, /, '  box length                                 :'
     &        , f10.3, /)
99004 format ('  model parameters: ', /, '     epsilon: ', f5.3, /, 
     &        '     sigma  : ', f5.3, /, '     mass   : ', f5.3,/)
99007 format (' requested density: ', f5.2, 
     &        ' different from density on disk: ', f5.2, /, 
     &        ' rescaling of coordinates!')
      end



      subroutine sample(i, en, vir, press)
c
c      write quantities (pressure and energy) to file
c
c      ener (input) : total energy
c      vir  (input) : total virial
c
      implicit none
      include 'params.inc'
      integer i
      real en, enp, vir, press,vol
 
      if (npart.ne.0) then
         enp   = en/real(npart)
         vol   = box**3
         press = (real(npart)/vol)/beta + vir/(3.0*vol)
      else
         enp   = 0.0
         press = 0.0
      end if

      write (30, '(i6,2e16.7)') i, enp, press

      return
      end



      subroutine store(iout, dr)
c
c     writes configuration to disk
c
c     iout (input) file number
c     dr   (input) maximum displacement
c
      implicit none
      include 'params.inc'
      integer iout, i
      real dr
 
      write (iout, *) npart, box, dr

      do i = 1, npart
         write (iout, '(3e22.14)') x(i), y(i), z(i)
      end do

      rewind (iout)

      return
      end



      subroutine toterg(ener, vir)
c
c     calculates total energy of the system
c     only used in the beginning or at the end of the program
c
c     ener (output) : total energy
c     vir  (output) : total virial
c
      implicit none
      include 'params.inc'
      real xi, yi, zi, ener, eni, viri, vir
      integer i, jb
 
      ener = 0.0
      vir  = 0.0

c     --sum over pairs: 1 <= i < jb <= n

      do i = 1, npart - 1
         xi = x(i)
         yi = y(i)
         zi = z(i)
         jb = i + 1

         call eneri(xi, yi, zi, i, jb, eni, viri)

         ener = ener + eni
         vir  = vir + viri
      end do

      return
      end



      subroutine writepdb(iout)
      implicit none
      include 'params.inc'
      integer iout,case,i

      data case/ 0/

      save case

      case = case + 1

      write(iout,'(a,i9)') 'pdb output',case

      do i=1,npart

         write(iout,'(a6,i5,1x,a4,1x,a3,1x,a1,i4,4x,3f8.3)') 
     &		'ATOM  ',i,'   C','XYZ','C',case,x(i),y(i),z(i)

      enddo

      write(iout,'(a)') 'end case'

      return
      end
