c
c     Metropolis evolution of a Lennard Jones 13 atom cluster
c      interacting with a Java 3D environment (JIMD)
c
c     
c
c

      program mc_lj  
      include 'lj.inc'
      include "jimd_f.h"        ! for interative molecular dynamics 
      real *8 xx(n),xy(n),xz(n)
      idum=2005
 1    call init(x,y,z,v,h)
      
      idum=2002                 ! seed for random number generator
      do i=1,nstep
         ist=ist+1              ! to compute cv
         call propagate(x,y,z,v,xo,yo,zo,vo,h,idum)
         
         if(mod(i,ipasinf).eq.0) then ! display system
            bx=radius2/2  
            do ion=1,n
               xx(ion)=x(ion)+bx
               xy(ion)=y(ion)+bx
               xz(ion)=z(ion)+bx
               call JIMD_pack_particle(ion,xx(ion),xy(ion),xz(ion))
            enddo
            call JIMD_send_data
            
         endif
         
      enddo
      
      goto 1
      call JIMD_finalize
      end
c     
C     ************************************************************
c
      subroutine temperature_changed(new_temperature)
c     
C     ************************************************************

      include 'lj.inc'
      real*8 new_temperature
      temp=new_temperature
    
      end
c  
c     
C     ************************************************************
   
      subroutine register_temperature_observer
c
C     ************************************************************

      include "jimd_f.h"
	 include 'lj.inc'
      external temperature_changed


      real*8 min_value, max_value, initial_value, step_size
      integer jimdstatus

      min_value     = 0.001
      max_value     = 3.0
      initial_value = temp
      step_size     = 0.01

      call JIMD_register_parameter_observer(temperature_changed,
     :     "Temperature",min_value,max_value,initial_value,step_size,
     :     jimdstatus)
      if(jimdstatus .ne. JIMD_SUCCESS) then
         write(*,*) 'Error while registering temperature observer'
      endif
	 return
	 end

C
C------------------init----------------------------------
Cc     
C     ************************************************************

      subroutine init(x,y,z,v,h)
c     
C     ************************************************************

      include 'lj.inc'
      real *8 radiu
      radiu=radius
      ib=radius2
      call JIMD_simple_init(n,ib,'TestParticle',radiu,4729,
     :     jimdstatus)
      call REGISTER_TEMPERATURE_OBSERVER
      open(10,file='xyz.in',status='old')
      do i=1,n
         read(10,*) x(i),y(i),z(i)
      enddo
      close(10)
      open(10,file='temp.in',status='old')
         read(10,*) temp
         h=(16./(1.-sqrt(temp)))**(1./6.)
         h=h-(16./(1.+sqrt(Temp)))**(1./6.)
         h=0.45*h   ! optimal step size
      close(10)

      open(unit=41,file='cv.out',status='UNKNOWN')
      
      
      do i=1,n
         v(i)=0.
      enddo
      
      do i=1,n
         do j=i+1,n
            w=vlj(x,y,z,i,j)
            v(i)=v(i)+w
            v(j)=v(j)+w
         enddo
      enddo
      
      return
      end
c
c--------------LJ potential-------------------
c
      function vlj(x,y,z,i,j)
      include 'lj.inc'
      r2=(x(i)-x(j))**2+(y(i)-y(j))**2+(z(i)-z(j))**2
      r6=1./r2**3
      vlj=4.*r6*(r6-1.)
      return
      end
c
c------------------propagate--------------------
c

      subroutine propagate(x,y,z,v,xo,yo,zo,vo,h,idum)
      include 'lj.inc'
      
      do k=1,n
         xo(k)=x(k)
         yo(k)=y(k)
         zo(k)=z(k)
         vo(k)=v(k)
      enddo

      i=1+n*ran0(idum)

      do k=1,n
         if (i.ne.k) v(k)=v(k)-vlj(x,y,z,i,k)
      enddo
      
      ddx=h*(ran0(idum)-0.5)    !random step
      ddy=h*(ran0(idum)-0.5)
      ddz=h*(ran0(idum)-0.5)
      ddxn=ddx/n 
      ddyn=ddy/n
      ddzn=ddz/n

      rmax=0.
      do k=1,n
         if (k.ne.i) then
            x(k)=x(k)-ddxn      !to conserve center of mass
            y(k)=y(k)-ddyn
            z(k)=z(k)-ddzn
         else
            x(k)=x(k)-ddxn+ddx
            y(k)=y(k)-ddyn+ddy
            z(k)=z(k)-ddzn+ddz
         endif
         r2=x(k)**2+y(k)**2+z(k)**2
         rmax=max(r2,rmax)      !to ensure atoms stay around
      enddo

      v(i)=0.
      do k=1,n
         if (k.ne.i) then
            w=vlj(x,y,z,k,i)
            v(k)=v(k)+w
            v(i)=v(i)+w
         endif
      enddo

      de=v(i)-vo(i)
      w=min(1.,exp(-de/temp))
      if (ran0(idum).gt.w.or.rmax.gt.radius2) then   ! reject move
         do k=1,n
            x(k)=xo(k)
            y(k)=yo(k)
            z(k)=zo(k)
            v(k)=vo(k)
         enddo
      endif

      return
      end
      
c
c----------ran0-------------
c
c completely basic (and wrong) random number generator
c replace with your own !
c
      function ran0(idum)
      idum=123456789*idum+987654321
      idum=mod(idum,2**30)
      ran0=idum/(2**30*1.0)
      if(ran0.lt.0.0) ran0=ran0+1.0
      return
      end

