PROGRAM demon
*     demon algorithm for the d = 1 Ising model in zero magnetic field 
      IMPLICIT NONE
      REAL*8 E, Ecum, Ed, Edcum, M, M2cum, Mcum
      INTEGER IR_, N, spin(1000), accept, imcs, mcs
      DATA accept/0/
      CALL GWopen(IR_, 0)
      CALL initial(N,spin,E,Ed,M,mcs,Ecum,Edcum,Mcum,M2cum)
      CALL setupscreen(N,spin)
      DO imcs = 1, mcs
         CALL change(N,spin,E,Ed,M,accept)
         CALL mydata(E,Ed,M,Ecum,Edcum,Mcum,M2cum)
      END DO
      CALL averages(N,Ecum,Edcum,Mcum,M2cum,mcs,accept)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(N,spin,E,Ed,M,mcs,Ecum,Edcum,Mcum,M2cum)
      IMPLICIT NONE
      REAL*8 Etot, E, Ed, M, Ecum, Edcum, Mcum, M2cum
      REAL dmy, rnd
      INTEGER isite, N, spin(1000), mcs
      dmy = rnd(-1)
      WRITE(*,'(A,$)') 'number of spins = '
      READ(*,*) N
*     choose total energy to be multiple of 4J
*     coupling constant J is unity
      WRITE(*,'(A,$)') 'desired total energy = '
      READ(*,*) Etot
      Etot = 4*int(Etot/4)
      WRITE(*,'(A,$)') 'number of Monte Carlo steps per spin = '
      READ(*,*) mcs
*     initial configuration of spins in minimum energy state
      DO isite = 1, N
         spin(isite) = 1
      END DO
      M = N                                 ! net magnetization
*     compute initial system energy
      E = -N                                ! periodic boundary conditions
      Ed = (Etot - E)
      WRITE(*,*) 'total energy = ', E + Ed
*     initialize sums
      Ecum = 0
      Edcum = 0
      Mcum = 0
      M2cum = 0
      END
C
      SUBROUTINE setupscreen(N,spin)
      IMPLICIT NONE
      REAL pi, r, xwin, ywin
      PARAMETER(pi = 3.1415926)
      INTEGER IR_, i, N, spin(1000)
      r = N/(2*pi)
      CALL compute_aspect_ratio(r+2,xwin,ywin)
      CALL GWindow(IR_, -xwin, -ywin, xwin, ywin)
      DO i = 1, N
         CALL showspin(N,spin(i),i)
      END DO
      END
C
      SUBROUTINE change(N,spin,E,Ed,M,accept)
      IMPLICIT NONE
      REAL*8 de, left, right, E, Ed, M
      REAL rnd
      INTEGER i, isite, N, spin(1000), accept
      DO i = 1, N
         isite = int(rnd(0)*N + 1)          ! random spin
*        determine neighboring spin values
         IF(isite .EQ. 1) THEN
            left = spin(N)
         ELSE
            left = spin(isite - 1)
         END IF
         IF(isite .EQ. N) THEN
            right = spin(1)
         ELSE
            right = spin(isite + 1)
         END IF
*        trial energy change 
         de =  2*spin(isite)*(left + right)
         IF(de .LE. Ed) THEN
*           spin flip dynamics
            spin(isite) = -spin(isite)
            M = M + 2*spin(isite)
            accept = accept + 1             ! number of changes accepted
            Ed = Ed - de
            E = E + de
         END IF
         CALL showspin(N,spin(isite),isite)
      END DO
      END
C
      SUBROUTINE mydata(E,Ed,M,Ecum,Edcum,Mcum,M2cum)
*     accumulate data
      IMPLICIT NONE
      REAL*8 E, Ed, M, Ecum, Edcum, Mcum, M2cum
      Ecum = Ecum + E
      Edcum = Edcum + Ed
      Mcum = Mcum + M
      M2cum = M2cum + M*M
      END
C
      SUBROUTINE averages(N,Ecum,Edcum,Mcum,M2cum,mcs,accept)
      IMPLICIT NONE
      REAL*8 Eave, Edave, M2ave, Mave, T, accept_prob, norm, Ecum,
     +       Edcum, Mcum, M2cum
      INTEGER N, mcs, accept
      norm = 1.0D0/mcs                          ! collected data after every attempt
      Edave = Edcum*norm
      WRITE(*,*) 'mean demon energy =', Edave
      T = 4/log(1 + 4/Edave)
      WRITE(*,*) 'T =', T
      Eave = Ecum*norm
      WRITE(*,*) '<E> = ', Eave
      Mave = Mcum*norm
      WRITE(*,*) '<M> =', Mave
      M2ave = M2cum*norm
      WRITE(*,*) '<M*M> =', M2ave
      accept_prob = accept*norm/N
      WRITE(*,*) 'acceptance probability =', accept_prob
      END
C
      SUBROUTINE showspin(N,dir,isite)
      IMPLICIT NONE
      REAL pi, r, theta, x, y
      PARAMETER(pi = 3.141592)
      INTEGER IR_, N, dir, isite
      r = N/(2*pi)
      theta = isite/r
      x = r*cos(theta)
      y = r*sin(theta)
      CALL GWerase(IR_, isite, 0)
      CALL GWsetogn(IR_, 0, isite)
      IF(dir .EQ. 1) THEN
         CALL GWsrect(IR_, x-0.25, y-0.25, x+0.25, y+0.25, 13)
      ELSE
         CALL GWsrect(IR_, x-0.25, y-0.25, x+0.25, y+0.25, 16)
      END IF
      END