PROGRAM conduct
! many demon algorithm for Ising chain
! heat added at spin 1 and subtracted at spin N
DIM spin(1000),Ed(1000),Edsum(1000),Msum(1000)
CALL initial(N,spin(),nmcs,nvisit)
FOR imcs = 1 to nmcs
    CALL change(N,spin(),Ed(),accept)
    IF mod(imcs,nvisit) = 0 then CALL heat(N,spin(),Edsum())
    CALL data(N,spin(),Ed(),Edsum(),Msum())
NEXT imcs
CALL output(N,Edsum(),Msum(),nmcs,accept)
END

SUB initial(N,spin(),nmcs,nvisit)
    RANDOMIZE
    INPUT prompt "number of spins = ": N
    INPUT prompt "number of MC steps per spin = ": nmcs
    INPUT prompt "MC steps between updates of end spins = ": nvisit
    ! initial random configuration
    FOR i = 1 to N
        IF rnd > 0.5 then
           LET spin(i) = 1
        ELSE
           LET spin(i) = -1
        END IF
    NEXT i
END SUB

SUB change(N,spin(),Ed(),accept)  ! spin flip dynamics
    ! do one Monte Carlo step
    FOR i = 2 to N - 1
        ! pick spin at random from spins 2 to N - 1
        LET isite = int(rnd*(N - 2) + 2)
        ! trial energy change
        LET de = 2*spin(isite)*(spin(isite - 1) + spin(isite + 1))
        IF de <= Ed(isite) then
           LET spin(isite) = - spin(isite)
           LET accept = accept + 1
           LET Ed(isite) = Ed(isite) - de
        END IF
    NEXT i
END SUB

SUB 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
    IF (spin(1)*spin(2) = 1) and (spin(N)*spin(N-1) = -1) then
       LET Edsum(1) = Edsum(1) + 2
       LET Edsum(N) = Edsum(N) - 2
       LET spin(1) = -spin(1)
       LET spin(N) = -spin(N)
    END IF
END SUB

SUB data(N,spin(),Ed(),Edsum(),Msum())
    FOR i = 2 to N - 1
        LET Edsum(i) = Edsum(i) + Ed(i)
        LET Msum(i) = Msum(i) + spin(i)
    NEXT i
END SUB

SUB output(N,Edsum(),Msum(),nmcs,accept)
    LET norm = 1/nmcs
    LET accept_prob = accept*norm/(N - 2)
    PRINT "acceptance probability = "; accept_prob
    PRINT
    PRINT tab(2);"i";tab(16);"Ed(i)";tab(35);"T";tab(46);"M(i)"
    PRINT
    FOR i = 2 to N-1
        LET edave = Edsum(i)*norm
        LET temperature = 0
        IF Edave <> 0 then
           IF (1 + 4/Edave) > 0 then
              LET temperature = 4/log(1 + 4/Edave)
           END IF
        END IF
        LET M = Msum(i)*norm
        PRINT i,Edave,temperature,M
    NEXT i
END SUB