C----------------------------------------------------------------------C
C                                                                      C
C       Metropolis Monte Carlo update of the Ising model               C
C      (Metropolis et al., J. Chem. Phys. 21 (1953) 1087).             C
C                                                                      C
C    Based on the code in the book "Statistical Field Theory",         C
C    G. Parisi, (Addison-Wesley, 1988).                                C
C                                                                      C
C       Paul Coddington, August 1992, June 1993                        C
C                                                                      C
C----------------------------------------------------------------------C

      PROGRAM METROPOLIS
      INCLUDE "common.h"
      INTEGER i,seed,start

C
C  Initialize.
C
C --- Read in parameters.
      CALL input(seed,start)
 
C --- Initialize the random number generator.
      CALL randinit(seed)
 
C --- Initialize the spins.
      CALL init_spins(start)
 
C --- Initialize the measurement data.
      CALL init_measure
 
C --- Set up periodic boundary conditions.
      CALL neighbors

C
C  Thermalize the system.
C
      CALL metupdate(therm_sweeps)

C
C  Metropolis update, taking measurements every 'sweeps' iterations.
C
      DO i=1,measurements
        CALL metupdate(sweeps)
        CALL measure
      ENDDO

C
C  Print out the results.
C
      CALL output

      STOP
      END




C----------------------------------------------------------------------C
C                                                                      C
C    The Metropolis Monte Carlo update.                                C
C                                                                      C
C    For the Ising model new_spin = - old_spin, so                     C
C    the old energy   oldE = - old_spin * (sum of neighboring spins)   C
C    the new energy   newE = - new_spin * (sum of neighboring spins)   C
C                          = - oldE                                    C
C                                                                      C
C----------------------------------------------------------------------C

      SUBROUTINE metupdate(iterations)
      INCLUDE "common.h"
      INTEGER iterations
      INTEGER n, i,j, new_spin, newE,oldE,deltaE
      REAL rand
	
      DO n=1,iterations

        accept = 0

        DO i=1,side
        DO j=1,side
C
C  New trial spin (= - old spin)
C
          new_spin = - spin(i,j)
C
C  Calculate old and new energy (new energy = - old energy).
C
          oldE = - spin(i,j) * 
     &             ( spin(i,pos(j)) + spin(i,neg(j)) +
     &               spin(pos(i),j) + spin(neg(i),j) )
          newE = - oldE
          deltaE = newE - oldE
C
C  Metropolis accept/reject 
C 
          IF ( ( deltaE.LE.0 ) .OR.
     &         ( EXP(-beta*deltaE).GT.rand() ) ) THEN
            spin(i,j) = new_spin
            accept = accept + 1
          ENDIF

        ENDDO
        ENDDO

      ENDDO
      RETURN 
      END


