      PROGRAM conduct
*     many demon algorithm for Ising chain
*     heat added at spin 1 and subtracted at spin N
      IMPLICIT NONE
      REAL*8 Ed(1000), Edsum(1000), Msum(1000), accept
      INTEGER N, imcs, nmcs, nvisit, spin(1000)
      DATA Ed, Edsum, Msum/3000*0.0D0/, accept/0.0D0/
      CALL initial(N,spin,nmcs,nvisit)
      DO imcs = 1, nmcs
         CALL change(N,spin,Ed,accept)
         IF(mod(imcs,nvisit) .EQ. 0) CALL heat(N,spin,Edsum)
         CALL data(N,spin,Ed,Edsum,Msum)
      END DO
      CALL output(N,Edsum,Msum,nmcs,accept)
      END
C
      SUBROUTINE initial(N,spin,nmcs,nvisit)
      IMPLICIT NONE
      REAL dmy, rnd
      INTEGER i, nmcs, N, spin(1000), nvisit
      dmy = rnd(-1)
      WRITE(*,'(A,$)') 'number of spins = '
      READ(*,*) N
      WRITE(*,'(A,$)') 'number of MC steps per spin = '
      READ(*,*) nmcs
      WRITE(*,'(A,$)') 'MC steps between updates of end spins = '
      READ(*,*) nvisit
*     initial random configuration
      DO i = 1, N
         IF(rnd(0) .GT. 0.5D0) THEN
            spin(i) = 1
         ELSE
            spin(i) = -1
         END IF
      END DO
      END
C
      SUBROUTINE change(N,spin,Ed,accept)
*     do one Monte Carlo step
      IMPLICIT NONE
      REAL*8 de, Ed(1000), accept
      REAL rnd
      INTEGER isite, i, N, spin(1000)
      DO i = 2, N - 1
*        pick spin at random from spins 2 to N - 1
         isite = int(rnd(0)*(N - 2) + 2)
*        trial energy change
         de = 2*spin(isite)*(spin(isite - 1) + spin(isite + 1))
         IF(de .LE. Ed(isite)) THEN
            spin(isite) = - spin(isite)
            accept = accept + 1
            Ed(isite) = Ed(isite) - de
         END IF
      END DO
      END
C
      SUBROUTINE heat(N,spin,Edsum)
*     attempt to add energy at spin 1 and remove it at spin N
*     possible only if spins 1 and 2 are aligned
*     and spins N and N - 1 are not aligned
      IMPLICIT NONE
      REAL*8 Edsum(1000)
      INTEGER N, spin(1000)
      IF((spin(1)*spin(2) .EQ. 1) .AND. (spin(N)*spin(N-1) .EQ. -1))
     + THEN
         Edsum(1) = Edsum(1) + 2
         Edsum(N) = Edsum(N) - 2
         spin(1) = -spin(1)
         spin(N) = -spin(N)
      END IF
      END
C
      SUBROUTINE data(N,spin,Ed,Edsum,Msum)
      IMPLICIT NONE
      REAL*8 Ed(1000), Edsum(1000), Msum(1000)
      INTEGER i, N, spin(1000)
      DO i = 2, N - 1
         Edsum(i) = Edsum(i) + Ed(i)
         Msum(i) = Msum(i) + spin(i)
      END DO
      END
C
      SUBROUTINE output(N,Edsum,Msum,nmcs,accept)
      IMPLICIT NONE
      REAL*8 Edave, M, accept_prob, norm, temperature,
     +       Edsum(1000), Msum(1000), accept
      INTEGER i, N, nmcs
      norm = 1.0D0/nmcs
      accept_prob = accept*norm/(N - 2)
      WRITE(*,*) 'acceptance probability = ', accept_prob
      WRITE(*,*)
      WRITE(*,'(4A13)') 'i', 'Ed(i)', 'T', 'M(i)'
      WRITE(*,*)
      DO i = 2, N-1
         Edave = Edsum(i)*norm
         temperature = 0
         IF(Edave .NE. 0) THEN
            IF((1 + 4/Edave) .GT. 0) THEN
               temperature = 4/log(1 + 4/Edave)
            END IF
         END IF
         M = Msum(i)*norm
         WRITE(*,'(I13,3F13.6)') i, Edave, temperature, M
      END DO
      END
