      program mc_npt
c
c  	--Lennard-Jones fluid
c
      implicit none
      integer iseed, equil, prod, nsamp, ii, icycl, ndispl, attempt, 
     &        nacc, ncycl, nmoves, imove, nvol, accv, attv
      real en, ent, vir, virt, vmax, dr, ran2, ran, succ
 
c     --common blocks declaration:
 
      include 'npt_params.inc'
 
      write (6, *) '**************** mc_npt ***************'

c     ---initialize sysem
 
      call readdat(equil, prod, nsamp, ndispl, dr, nvol, vmax, iseed, 
     &             succ)
      nmoves = ndispl + nvol

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) then
		write (6, *) ' start equilibration '
		write (66, *) ' start equilibration '
	    endif
         else
            ncycl = prod
            if (ncycl.ne.0) then
		write (6, *) ' start production '
		write (66, *) ' start production '
            endif
         end if
         attempt = 0
         nacc = 0
         attv = 0
         accv = 0

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

         call adjust(attempt, nacc, dr, attv, accv, vmax, succ)
         do icycl = 1, ncycl
            do imove = 1, nmoves
               ran = ran2(iseed)*(ndispl+nvol)
               if (ran.lt.ndispl) then
c                 ---attempt to displace a particle
                  call mcmove(en, vir, attempt, nacc, dr, iseed)
               else
c                 ---attempt to change the volume
                  call mcvol(en, vir, attv, accv, vmax, iseed)
               end if
            end do

c           ---sample averages

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

c              ---write intermediate configuration to file

               call store(8, dr, vmax)

c              ---adjust maximum displacements

               call adjust(attempt, nacc, dr, attv, accv, vmax, succ)
            end if
         end do

         if (ncycl.ne.0) then
            if (attempt.ne.0) write (6, 99003) attempt, nacc, 
     &                               100.*float(nacc)/float(attempt)
            if (attv.ne.0) write (6, 99004) attv, accv, 100.*float(accv)
     &                            /float(attv)
c           ---test total energy
            call toterg(ent, virt)
            if (abs(ent-en).gt.1.d-6) then
               write (6, *) 
     &                    ' ######### problems energy ################ '
            end if
            if (abs(virt-vir).gt.1.d-6) then
               write (6, *) 
     &                    ' ######### problems virial ################ '
            end if
            write (6, 99002) ent, en, ent - en, virt, vir, virt - vir
         end if

      end do
      call store(21, dr, vmax)
      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, '%)')
99004 format (' number of att. to chan. volume    : ', i10, /, 
     &        ' success: ', i10, '(= ', f5.2, '%)')
      end



      subroutine adjust(attemp, nacc, dr, attv, accv, vmax, succ)
c
c     sets maximum displacement and maximum volume change
c     such that 50% of the move will be accepted
c
      implicit none
      include 'npt_params.inc'
      integer attemp, nacc, attempp, naccp, accvp, attv, accv, attvp
      real dro, frac, vmaxo, dr, vmax, succ
      save naccp, attempp, attvp, accvp
 
c     ---displacement:

      if (attemp.eq.0.or.attempp.ge.attemp) then
         naccp = nacc
         attempp = attemp
      else
         frac = float(nacc-naccp)/float(attemp-attempp)
         dro = dr
         dr = dr*abs(frac/(succ/100.0))
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

c     ---volume:

      if (attv.eq.0.or.attvp.ge.attv) then
         accvp = accv
         attvp = attv
      else
         frac = float(accv-accvp)/float(attv-attvp)
         vmaxo = vmax
         vmax = vmax*abs(frac/(succ/100.))
c        ---limit the change:
         if (vmax/vmaxo.gt.1.5) vmax = vmaxo*1.5
         if (vmax/vmaxo.lt.0.5) vmax = vmaxo*0.5
         write (6, 99002) vmax, vmaxo, frac, attv - attvp, accv - accvp
c        ---store nacc and attemp for next use
         accvp = accv
         attvp = attv
      end if

      return
99001 format (' max. displ. set to      : ', f6.3, ' (old : ', f6.3, 
     &        ')', /, ' frac. acc.: ', f4.2, ' attempts: ', i7, 
     &        ' succes: ', i7)
99002 format (' max. vol. chan. set to: ', f6.3, ' (old : ', f6.3, ')', 
     &        /, ' frac. acc.: ', f4.2, ' attempts: ', i7, ' succes: ', 
     &        i7)
      end


 
      function corp(r, rho)
      implicit none
      include 'npt_params.inc'
      real sig3, ri3, r, corp, rho
 
      sig3 = sig2*sqrt(sig2)
      ri3 = sig3/(r*r*r)
      corp = 4.*pi*eps4*(rho**2)*sig3*(2.*ri3*ri3*ri3/9.-ri3/3.)
      return
      end



      function coru(r, rho)
      implicit none
      include 'npt_params.inc'
      real sig3, ri3, r, coru, rho
 
      sig3 = sig2*sqrt(sig2)
      ri3 = sig3/(r*r*r)
      coru = 2.*pi*eps4*(rho*sig3)*(ri3*ri3*ri3/9.-ri3/3.)
      return
      end



      subroutine ener(en, vir, r2)
c
c  ---calculates energy (en) and virial (vir) for given
c     distance squared between (r2) two particles
 
      implicit none
      real r2, r2i, r6i, en, vir
      include 'npt_params.inc'
 
      if (r2.le.rc2) then
         r2i = sig2/r2
         r6i = r2i*r2i*r2i
         if (shift) then
            en = eps4*(r6i*r6i-r6i) - ecut
         else
            en = eps4*(r6i*r6i-r6i)
         end if
         vir = eps48*(r6i*r6i-0.5*r6i)
      else
         en = 0.
         vir = 0.
      end if
      return
      end



      subroutine eneri(xi, yi, zi, i, jb, en, vir)
      implicit none
      include 'npt_params.inc'
 
      real xi, yi, zi, en, dx, dy, dz, r2, vir, virij, enij
      integer i, j, jb
c
      en = 0.
      vir = 0.
      do j = jb, npart
         if (j.ne.i) then
            dx = xi - x(j)
            dy = yi - y(j)
            dz = zi - z(j)
            if (dx.gt.hbox) then
               dx = dx - box
            else
               if (dx.lt.-hbox) dx = dx + box
            end if
            if (dy.gt.hbox) then
               dy = dy - box
            else
               if (dy.lt.-hbox) dy = dy + box
            end if
            if (dz.gt.hbox) then
               dz = dz - box
            else
               if (dz.lt.-hbox) dz = dz + box
            end if
            r2 = dx*dx + dy*dy + dz*dz
            call ener(enij, virij, r2)
            en = en + enij
            vir = vir + virij
         end if
      end do
      vir = vir
      return
      end
 


      subroutine lattice(struc)
c
c     ---place `npart' particles on a lattice with density 'rho'
c
      implicit none
      include 'npt_params.inc'
      integer i, j, k, itel, n
      real del
      character struc*3
 
      del = (box**3)**(1./3.)
      if (struc.eq.'fcc') then
c        ---put particles on a cubic fcc-lattice
         n = nint((npart/4.)**(1./3.))
         write (6, *) ' generate fcc-lattice '
         if (n*n*n*4.ne.npart) then
            write (6, *) 
     &                  ' number of particles incomp. with fcc lattice '
     &                  , npart, n
            stop 'lattice'
         end if
         del = del/n
         print *, del, box
         itel = 0
         do i = 0, n - 1
            do j = 0, n - 1
               do k = 0, n - 1
                  x(itel+1) = k*del
                  y(itel+1) = j*del
                  z(itel+1) = i*del
 
                  x(itel+2) = (k+0.5)*del
                  y(itel+2) = (j+0.5)*del
                  z(itel+2) = i*del
 
                  x(itel+3) = k*del
                  y(itel+3) = (j+0.5)*del
                  z(itel+3) = (i+0.5)*del
 
                  x(itel+4) = (k+0.5)*del
                  y(itel+4) = j*del
                  z(itel+4) = (i+0.5)*del
                  itel = itel + 4
               end do
            end do
         end do
      else
c        ---put particels on a simple cubic lattice
         write (6, *) ' generate simple cubic lattice'
         n = int(npart**(1./3.)) + 1
         if (n.eq.0) n = 1
         del = del/float(n)
         itel = 0
         do i = 0, n - 1
            do j = 0, n - 1
               do k = 0, n - 1
                  if (itel.lt.npart) then
                     itel = itel + 1
                     x(itel) = k*del
                     y(itel) = j*del
                     z(itel) = i*del
                  end if
               end do
            end do
         end do
      end if
      write (6, 99001) itel

      return
99001 format (' initialisation on lattice: ', /, i10, 
     &        ' particles placed on a lattice')
      end
 
 

      subroutine mcmove(en, vir, attempt, nacc, dr, iseed)
c
c     attempts to displace a randomly selected particle
c
      implicit none
      include 'npt_params.inc'
      real enn, eno, en, ran2, xn, yn, zn, viro, virn, vir, 
     &                 dr
      integer o, attempt, nacc, jb, iseed
 
      attempt = attempt + 1
      jb = 1
c     ---select a particle at random
      o = int(npart*ran2(iseed)) + 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) + (ran2(iseed)-0.5)*dr
      yn = y(o) + (ran2(iseed)-0.5)*dr
      zn = z(o) + (ran2(iseed)-0.5)*dr
c     ---calculate energy new configuration:
      call eneri(xn, yn, zn, o, jb, enn, virn)
c     ---acceptance test
      if (ran2(iseed).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) xn = xn + box
         if (xn.gt.box) xn = xn - box
         if (yn.lt.0) yn = yn + box
         if (yn.gt.box) yn = yn - box
         if (zn.lt.0) zn = zn + box
         if (zn.gt.box) zn = zn - box
         x(o) = xn
         y(o) = yn
         z(o) = zn
      end if
      return
      end



      subroutine mcvol(en, vir, attempt, acc, vmax, iseed)
c
c     attempts to change the volume
c
      implicit none
      include 'npt_params.inc'
      real enn, en, ran2, vir, virn, boxn, yy, vo, vn, coru, 
     &                 lnvn, vmax, rhon, rhoo, rcn, f, arg, eno, virt
      integer attempt, acc, i, iseed
 
      attempt = attempt + 1
c     ---calulate new volume by making random walk in ln v
      vo = box*box*box
      lnvn = log(vo) + (ran2(iseed)-0.5)*vmax
      vn = exp(lnvn)
      boxn = vn**(1./3.)
      if (.not.shift) then
c        ---calculate new energy using scaling
         if (tailco) then
c           ---substract tail correction
            rhoo = npart/vo
            eno = en - npart*coru(rc, rhoo)
         end if
         yy = (vo/vn)*(vo/vn)
         enn = eno*yy*(2-yy) - vir*yy*(1-yy)/6
         virn = -12*eno*yy*(yy-1) + vir*yy*(2*yy-1)
         if (tailco) then
c            ---add tail correction
            rhon = npart/vn
            rcn = rc*boxn/box
            enn = enn + npart*coru(rcn, rhon)
         end if
      else
c        ---determine new coordinates
         f = boxn/box
         do i = 1, npart
            x(i) = f*x(i)
            y(i) = f*y(i)
            z(i) = f*z(i)
         end do
         box = boxn
         hbox = 0.5*box
         rc = rc*f
         rc2 = rc*rc
c        ---recalculate new energy shift (dirty trick)
         ecut = 0.
         call ener(ecut, virt, rc2)
c        ---calculate new total energy
         call toterg(enn, virn)
      end if
c     ---acceptance:
      arg = exp(-beta*(p*(vn-vo)-(npart+1)*log(vn/vo)/beta+enn-en))
      if (ran2(iseed).lt.arg) then
c        ---accepted
         acc = acc + 1
         en = enn
         vir = virn
         if (shift) then
            return
         end if
c        ---rescale coordinates
         f = boxn/box
         do i = 1, npart
            x(i) = f*x(i)
            y(i) = f*y(i)
            z(i) = f*z(i)
         end do
         box = boxn
         hbox = box/2
         rc = rc*f
         rc2 = rc*rc
      else if (shift) then
c        ---restore the old configuration
         f = 1/f
         do i = 1, npart
            x(i) = f*x(i)
            y(i) = f*y(i)
            z(i) = f*z(i)
         end do
         box = boxn*f
         hbox = 0.5*box
         rc = rc*f
         rc2 = rc*rc
c        ---recalculate energy shift (dirty trick)
         ecut = 0
         call ener(ecut, virt, rc2)
      end if

      return
      end



      subroutine readdat(equil, prod, nsamp, ndispl, dr, nvol, vmax, 
     &                   iseed, succ)
c
c     reads input data and model parameters
c
c     ---input parameters: file: fort.15
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    iseed      : seed random number generator
c    dr         : maximum displacement
c    vmax       : maximum volume change
c    succ       : optimal percentance of accepted attemps
c                 the program adjusts vmax or dr in just a way that
c                 on average succ% of the moves are accepted
c    ndispl     : number of attemps to displace a particle per mc cycle
c    nvol       : number of attemps to change the volume  per mc cycle
c    npart      : total numbero fo particles
c    temp       : temperature
c    rho        : density
c    p          : pressure
c    struc      : structure to initalize the lattice
c               = 'fcc' fcc-lattice (face centered cubic)
c                 otherwise simple cubic bcc-lattice (body centered cubic)
c
c
c     ---input parameters: file: fort.25
c    tailco = .true. : tail corrections are applied
c             .false.: no-tail corrections
c    shift  = .true. : potential is shifted
c             .false.: potential is not-shifted
c    eps    = epsilon lennard-jones potential
c    sig    = sigma 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.15!!)
c    dr     = optimized maximum displacement old configurations
c    vmax   = optimized maximum volume change old configurations
c    x(1),y(1),z(1)            : position first particle 1
c        ...
c    x(npart),y(npart),z(npart): position particle last particle
 
      implicit none
      include 'npt_params.inc'
      integer ibeg, equil, prod, i, ndispl, nsamp, nvol, iseed
      real eps, sig, coru, corp, vir, boxf, hboxf, rhof, 
     &                 rho, dr, vmax, succ
      character struc*3
 
c     ---read simulation data
 
      read (15, *) ibeg, equil, prod, nsamp, iseed
      read (15, *) dr, vmax, succ
      read (15, *) ndispl, nvol
      read (15, *) npart, temp, rho, p, struc
 
      if (npart.gt.npmax) then
         write (6, *) ' error: number of particles too large'
         stop
      end if
c     ---read model parameters
      read (15, *) tailco, shift
      read (15, *) eps, sig, mass, rc

c     ---read/generate configuration

      box = (npart/rho)**(1.0/3.0)
      hbox = 0.50*box
      if (ibeg.eq.0) then
c        ---generate configuration form lattice
         call lattice(struc)
      else
         write (6, *) ' read conf from disk '
         read (11, *) boxf, hboxf
         read (11, *) npart
         read (11, *) dr, vmax
         rhof = npart/boxf**3
         if (ibeg.eq.2) then
c           ---scale density to input value
            if (abs(boxf-box).gt.1.d-6) then
               write (6, 99007) rho, rhof
            end if
         else
c           ---use density from disk
            rho = rhof
            box = boxf
            hbox = hboxf
         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, nvol, vmax
      write (6, 99003) npart, temp, p, rho, box
      write (6, 99004) eps, sig, mass

c     ---calculate parameters:
      beta = 1./temp
      pi = 4.*atan(1.)

c     ---calculate cut-off radius potential

      rc = min(rc, hbox)
      rc2 = rc*rc
      eps4 = 4.*eps
      eps48 = 48.*eps
      sig2 = sig*sig
      if (shift) then

c     ---calculate energy of the shift (dirty trick)

         ecut = 0.
         call ener(ecut, vir, rc2)
         write (6, 99005) rc, ecut
      end if

      if (tailco) then
         write (6, 99006) rc, coru(rc, rho), corp(rc, rho)
      end if

      return
99001 format ('  number of equilibration cycles             :', i10, /, 
     &        '  number of production cycles                :', i10, /, 
     &        '  sample frequency                           :', i10, /)
99002 format ('  number of att. to displ. a part. per cycle :', i10, /, 
     &        '  maximum displacement                       :', f10.3, 
     &        /, '  number of att. to change volume  per cycle :', i10, 
     &        /, '  maximum change volume                      :', 
     &        f10.3, //)
99003 format ('  number of particles                        :', i10, /, 
     &        '  temperature                                :', f10.3, 
     &        /, '  pressure                                   :', 
     &        f10.3, /, '  density                                    :'
     &        , f10.3, /, 
     &        '  box length                                 :', f10.3, 
     &        /)
99004 format ('  model parameters: ', /, '     epsilon: ', f5.3, /, 
     &        '     sigma  : ', f5.3, /, '     mass   : ', f5.3)
99005 format (' simulations with truncated and shifted potential: ', /, 
     &        ' potential truncated at :', f10.3, /, 
     &        ' energy shif            :', f10.6, //)
99006 format (' simulations with tail correction: ', /, 
     &        ' potential truncated at  rc =:', f10.3, /, 
     &        ' tail corrections:   energy = ', f10.3, ' pressure ', 
     &        f10.3, /, /)
99007 format (' requested density: ', f5.2, 
     &        ' different from density on disk: ', f5.2, /, 
     &        ' rescaling of coordinates!')
      end
 


      subroutine sample(i, en, vir)
c
c     writes quantities to file
c
      implicit none
      include 'npt_params.inc'
      integer i
      real en, enp, vir, press, corp, vol, rho
 
      if (npart.ne.0) then
         enp = en/float(npart)
         vol = box**3
         rho = npart/vol
         press = rho/beta + vir/(3.*vol)
         if (tailco) press = press + corp(rc, rho)
      else
         rho = 0.
         enp = 0.
         press = 0.
      end if
      write (66, *) i, enp, press, rho
      return
      end



      subroutine store(iout, dr, vmax)
c
c     writes configuration to disk
c
      implicit none
      include 'npt_params.inc'
      integer iout, i
      real dr, vmax
 
      write (iout, *) box, hbox
      write (iout, *) npart
      write (iout, *) dr, vmax
      do i = 1, npart
         write (iout, *) x(i), y(i), z(i)
      end do
      rewind (iout)
      return
      end



      subroutine toterg(ener, vir)
c
c     ---calculates total energy
c
      implicit none
      include 'npt_params.inc'
 
      real xi, yi, zi, ener, eni, coru, viri, vir, rho
      integer i, jb
 
      ener = 0
      vir = 0
      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
c     ---add tail corrections
      if (tailco) then
         rho = npart/(box**3)
         ener = ener + npart*coru(rc, rho)
      end if
      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
