c------------------------------------------------------------------------------
      PROGRAM AHB1
c------------------------------------------------------------------------------
c     SSE simulation of the isotropic S=1/2 Heisenberg antiferromagnet on a 
c     simple cubic lattice in 1,2, or 3 dimensions. Uses deterministic loop 
c     updates introduced in Phys. Rev. B59, R14157 (1999). Calculates the
c     internal energy, specific heat, uniform susceptibility, spin stiffness,
c     staggered structure factor, and staggered susceptibility, using SSE 
c     estimators discussed in Phys. Rev. B56, 11678 (1997).
c------------------------------------------------------------------------------
c     INPUT DATA: Simulation parameters from file 'read.in' in the format
c     d,l,beta,nr,mstp,in,istp
c     d    = dimensionality (1,2, or 3) (must be same as dd in ahb1.h)
c     l    = linear system size (must be same as nx in ahb1.h)
c     beta = inverse temperature J/T (maxmimum approximately bm in ahb1.h)
c     nr   = number of bin averages to calculate
c     mstp = number of Monte Carlo steps per bin
c     in   = 0 to start from scratch, 1 to start from saved configuration 
c     istp = number of Monte Carlo steps for equilibration
c------------------------------------------------------------------------------
c     RANDOM NUMBER SEEDS: 4 integers from file 'rand.in, in the format
c     seed1
c     seed2
c     seed3
c     seed4
c------------------------------------------------------------------------------
c     OUTPUT DATA: Bin averages written to the files 'len.dat' and 'mag.dat' 
c     after completion of each bin (see comments in subroutine RESULTS). The 
c     current configuration is written to the file 'conf' after equilibration 
c     and after completion of each bin. Progress reports and possible error 
c     messages are written to the file 'log.txt'.
c------------------------------------------------------------------------------
c     Written by A. W. Sandvik (asandvik@abo.fi), Åbo Akademi University.
c     Last modifeid: Sept. 2, 2002.
c------------------------------------------------------------------------------

      INCLUDE 'ahb1.h'

      INTEGER i,j,nr,in,istp,mstp
      REAL*8  beta

      OPEN(UNIT=10,FILE='read.in',STATUS='old')
      READ(10,*)i,j,beta,nr,mstp,in,istp
      CLOSE(10)
      IF (i.NE.dd.OR.j.NE.nx) THEN
         CALL OPENLOG
         WRITE(12,*)'Wrong system size '
         WRITE(12,*)'D ',dd,i
         WRITE(12,*)'L ',nx,j
         CALL CLOSELOG
         STOP
      ENDIF

      CALL INITRAN
      IF (in.EQ.0) THEN
         CALL INITCONF
      ELSE
         CALL READCONF
      ENDIF
      ap1=0.5d0*beta*DFLOAT(nb)
      dp1=1.d0/(0.5d0*beta*DFLOAT(nb))
      CALL LATTICE

      DO i=1,istp
         CALL DUPDATE
         CALL LUPDATE
         CALL ADJLGTH
         IF (MOD(i,istp/10).EQ.0) THEN
            CALL OPENLOG
            WRITE(12,10)i,l
 10         FORMAT('Equilibrated. ',I7,'   L = : ',I8) 
            CALL CLOSELOG
         ENDIF
      ENDDO
      IF (istp.NE.0) CALL WRITECONF
      DO j=1,nr
         CALL CLEARDAT
         DO i=1,mstp
            CALL DUPDATE
            CALL LUPDATE
            CALL MEASURE
            IF (MOD(i,mstp/10).EQ.0) THEN
               CALL OPENLOG
               WRITE(12,20)j,i
 20            FORMAT(' Done Bin ',I5,'  Step ',I8) 
               CALL CLOSELOG
            ENDIF
         ENDDO
         CALL RESULTS(beta,mstp)
         CALL WRITECONF
      ENDDO

      END

c------------------------------------------------------------------------------
      SUBROUTINE DUPDATE
c------------------------------------------------------------------------------
c     Carries out a sweep of diagonal updates.
c------------------------------------------------------------------------------
c     MAIN DATA STRUCTURES AND VARIABLES:
c     str(i), i=1,...,l, contains the operators coded in the following way:
c     str(i)=0 corresponds to a fill-in unit operator, str(i)=2*b, b=1,...,nb,
c     is a diagonal operator on bond b, str(i)=2*b+1, b=1,...,nb, is an
c     off-diagonal operator on bond b. nh is the current number of non-0
c     elements in str() (the current expansion power). spn(s)=-1,1, are
c     the spin variables at sites s=1,...,N (-1 = down, 1 = up). bst(1,b)
c     and bst(2,b) are the sies (spins) connected by bond b.
c------------------------------------------------------------------------------

      INCLUDE 'ahb1.h'

      INTEGER i,b,op
      REAL*8  RAN,r,p

      DO i=1,l
         op=str(i)
         IF (op.EQ.0) THEN
            b=MIN(INT(RAN()*DFLOAT(nb))+1,nb)
            IF (spn(bst(1,b)).NE.spn(bst(2,b))) THEN
               IF (ap1.GE.DFLOAT(l-nh)) THEN
                  str(i)=2*b
                  nh=nh+1
               ELSEIF (RAN()*DFLOAT(l-nh).LE.ap1) THEN                  
                  str(i)=2*b
                  nh=nh+1
               ENDIF
            ENDIF
         ELSEIF (MOD(op,2).EQ.0) THEN        
            p=dp1*DFLOAT(l-nh+1)
            IF (p.GE.1.d0) THEN
               str(i)=0
               nh=nh-1
            ELSEIF (RAN().LE.p) THEN               
               str(i)=0
               nh=nh-1
            ENDIF
         ELSE
            b=op/2
            spn(bst(1,b))=-spn(bst(1,b))
            spn(bst(2,b))=-spn(bst(2,b))
         ENDIF
      ENDDO

      RETURN
      END

c-----------------------------------------------------------------------------
      SUBROUTINE LUPDATE
c-----------------------------------------------------------------------------
c     Constructs the linked vertex list and carries out deterministic operator
c     loop updates. Every loop is constructed once and each loop is flipped 
c     with probability 1/2. Maps the updated vertex list back to an operator 
c     string and makes the corresponding changes in the spin state. Flips free
c     spins with probability 1/2.
c------------------------------------------------------------------------------
c     DATA STRUCTURES AND CONSTRUCTION OF LINKED VERTEX LIST:
c     The array vrtx() is the linked list. The element vrtx(P), with P=p+i, 
c     p=0,...,nh-1, i=0,1,2,3, corresponds to leg i of vertex p. It contains
c     a link Q=q+j to another element, i.e., vrtx(P)=Q means that leg i of 
c     vertex p is linked to leg j of vertex q. In the construction of vrtx(),
c     two arrays, frst() and last(), are used. All elements frst(s),last(s),
c     s=1,...,N, are first initialized to -1. frst(s) will contain the first 
c     position P in vrtx() corresponding to an operation on site s, i.e., 
c     frst(s)=P=p+i means that leg i of vertex p is the first vertex-leg that 
c     acts on site s. Similarly, last(s) contains the last operation on s. 
c     Whereas frst(s) is set at most once during the construction of vrtx()
c     (never if no operator acts on site s), last(s) is updated during the
c     construction every time an operator acting on s is encountered, and is
c     used to determine what elements in vrtx() should be linked.
c------------------------------------------------------------------------------
c     LOOP CONSTRUCTION AND UPDATE:
c     Bits 30 and 31 in the elements of vrtx() are used as flags to
c     indicate that the corresponding vertex legs have been visited and/or
c     flipped. If bit 30 of vrtx(P), P=p+i, is set, then leg i of vertex
c     p has been visited but the spin has not been flipped. If instead
c     bit 31 is set, then the spin has also been flipped. The next element
c     to visit from an element p1=p+i (vertex p, entrance leg i) is obtaind 
c     using the XOR operation: p2=IEOR(p1,i) is the element corresponding
c     to the correct exit j, i.e., p2=p+j where (i,j) are entrance and
c     exit legs corresponding to the switch-and-reverse process through
c     the vertex. All even elements in vrtx() are tested for the status
c     of bits 30 and 31. A new loop is started whenever none of these bits
c     is set, and the decison of whether or not to flip (bit 30 or 31 to
c     be set in the elements visited) is made at random. When all loops
c     have been constructed, bit 31 in the even elements of vrtx() are used 
c     to update the operators in str(). The array frst() is used with
c     the 31-bit test to update the spins spn(), and those spins s that
c     are not acted on by any operator (frst(s)=-1) are flipped with
c     probability 1/2.
c------------------------------------------------------------------------------

      INCLUDE 'ahb1.h'

      INTEGER j,b,i0,i1,s0,s1,p0,p1,p2
      INTEGER vrtx(0:4*mm),frst(nn+1),last(nn+1)
      REAL*8  RAN

c**** Construction of the linked vertex list vrtx

      DO j=1,nn
         frst(j)=-1      
         last(j)=-1
      ENDDO
      i0=0
      i1=1
      DO j=1,l
         IF (str(j).NE.0) THEN
            b=str(j)/2
            s0=bst(1,b)
            s1=bst(2,b)
            p0=last(s0)
            p1=last(s1)
            IF (p0.NE.-1) THEN
               vrtx(p0)=i0
               vrtx(i0)=p0
            ELSE
               frst(s0)=i0
            ENDIF
            IF (p1.NE.-1) THEN
               vrtx(p1)=i1
               vrtx(i1)=p1
            ELSE
               frst(s1)=i1
            ENDIF
            last(s0)=i0+2
            last(s1)=i1+2
            i0=i0+4
            i1=i1+4
         ENDIF
      ENDDO
      DO s0=1,nn
         i0=frst(s0)
         IF (i0.NE.-1) THEN
            p0=last(s0)
            vrtx(p0)=i0
            vrtx(i0)=p0
         ENDIF
      ENDDO

c**** Loop update; every loop constructed, flipped with probability 1/2

      DO 20 p0=0,4*nh-1,2       
         IF (BTEST(vrtx(p0),30).OR.BTEST(vrtx(p0),31)) GOTO 20
         p1=p0
         b=30+MIN(INT(RAN()*2.d0),1)
 10      p2=IEOR(p1,1)
         vrtx(p1)=IBSET(vrtx(p1),b)
         p1=vrtx(p2)
         vrtx(p2)=IBSET(vrtx(p2),b)
         IF (p1.EQ.p0) GOTO 20
         GOTO 10
 20   CONTINUE

c**** Mapping of updated vertices to operator sequence

      j=0
      DO i0=1,l
         IF (str(i0).NE.0) THEN
            IF (BTEST(vrtx(j),31).XOR.BTEST(vrtx(j+2),31)) THEN
               str(i0)=IEOR(str(i0),1)
            ENDIF
            j=j+4
         ENDIF
      ENDDO

c**** Flipping of spins affected by loop updates, and random free spins 

      DO j=1,nn
         IF (frst(j).NE.-1) THEN
            IF (BTEST(vrtx(frst(j)),31)) spn(j)=-spn(j)
         ELSE
            IF (RAN().LT.0.5d0) spn(j)=-spn(j)
         ENDIF
      ENDDO

      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE MEASURE
c------------------------------------------------------------------------------
c     Accumulates measurements for the average of the expansion power nh
c     (len1) and its square (len2), the uniform susceptibility (usus), the
c     spin stiffness averaged over all directions (rhos), the staggered
c     structure factor (astr) and the staggered susceptibility (asus). In
c     some instances it is assumed that the total number of spins is even.
c------------------------------------------------------------------------------

      INCLUDE 'ahb1.h'

      REAL*8  len1,len2,usus,rhos,astr,asus
      COMMON/MBLOCK/len1,len2,usus,rhos,astr,asus

      INTEGER i,b,op,s1,s2,um,am,jj(0:2)
      REAL*8  astr1,asus1,rhos1

      am=0
      um=0
      DO i=1,nn
         um=um+spn(i)
         am=am+spn(i)*phs(i)
      ENDDO      
      um=um/2
      am=am/2
      astr1=DFLOAT(am)**2
      asus1=0.d0
      jj(0)=0
      jj(1)=0
      jj(2)=0
      DO i=1,l
         op=str(i)
         IF (MOD(op,2).EQ.1) THEN        
            b=op/2
            s1=bst(1,b)
            s2=bst(2,b)
            spn(s1)=-spn(s1)
            spn(s2)=-spn(s2)
            jj((b-1)/nn)=jj((b-1)/nn)+spn(s2)
            am=am+2*spn(s1)*phs(s1)
            asus1=asus1+DFLOAT(am)
            astr1=astr1+DFLOAT(am)**2
         ELSEIF (op.NE.0) THEN
            asus1=asus1+DFLOAT(am)
            astr1=astr1+DFLOAT(am)**2
         ENDIF
      ENDDO
      astr1=astr1/DFLOAT(nh+1)
      IF (nh.NE.0) THEN
         asus1=(asus1**2+DFLOAT(nh)*astr1)/(DFLOAT(nh)*DFLOAT(nh+1))
      ELSE
         asus1=astr1
      ENDIF

      len1=len1+DFLOAT(nh)
      len2=len2+DFLOAT(nh)**2
      usus=usus+(DFLOAT(um)**2)
      rhos=rhos+DFLOAT(jj(0))**2+DFLOAT(jj(1))**2+DFLOAT(jj(2))**2
      astr=astr+astr1
      asus=asus+asus1

      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE RESULTS(beta,mstp)
c------------------------------------------------------------------------------
c     Writes the accumulated bin averages to files len.dat and mag.dat.
c------------------------------------------------------------------------------

      INCLUDE 'ahb1.h'

      REAL*8  len1,len2,usus,rhos,astr,asus
      COMMON/MBLOCK/len1,len2,usus,rhos,astr,asus

      INTEGER mstp
      REAL*8  beta

      len1=len1/(DFLOAT(mstp)*beta*DFLOAT(nn))
      len2=len2/(DFLOAT(mstp)*(beta*DFLOAT(nn))**2)
      usus=beta*usus/(DFLOAT(mstp)*DFLOAT(nn))
      rhos=rhos/(DFLOAT(mstp)*beta*DFLOAT(nb))
      astr=astr/(DFLOAT(mstp)*DFLOAT(nn))
      asus=beta*asus/(DFLOAT(mstp)*DFLOAT(nn))

      OPEN(UNIT=10,FILE='len.dat',STATUS='unknown',ACCESS='append')
      WRITE(10,10)len1,len2
 10   FORMAT(F16.12,' ',F16.12)
      CLOSE(10)

      OPEN(UNIT=10,FILE='mag.dat',STATUS='unknown',ACCESS='append')
      WRITE(10,20)usus,rhos,astr,asus
 20   FORMAT(F12.8,' ',F12.8,' ',F14.8,' ',F18.8)
      CLOSE(10)

      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE CLEARDAT
c------------------------------------------------------------------------------
c     Resets the stored data variables after each bin.
c------------------------------------------------------------------------------

      REAL*8  len1,len2,usus,rhos,astr,asus
      COMMON/MBLOCK/len1,len2,usus,rhos,astr,asus

      len1=0.d0
      len2=0.d0
      rhos=0.d0
      astr=0.d0
      asus=0.d0
      usus=0.d0

      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE ADJLGTH
c------------------------------------------------------------------------------
c     Increases the cut-off l to nh+nh/4 if l is currently less than this.
c     Distributes the nh non-0 operators randomly among the new l positions.
c------------------------------------------------------------------------------

      INCLUDE 'ahb1.h'

      INTEGER i,j,l1
      REAL*8  RAN,r

      l1=nh+nh/4
      IF (l1.LE.l) THEN
         RETURN
      ELSEIF (l1.GT.ll) THEN
         CALL OPENLOG
         WRITE(12,*)'L exceeded maximum ',l1,ll
         CALL CLOSELOG
         STOP
      ENDIF
      j=0
      DO i=1,l
         IF (str(i).NE.0) THEN
            j=j+1
            str(j)=str(i)
         ENDIF
      ENDDO
      DO i=nh+1,l1
         str(i)=0
      ENDDO
      l=l1
      j=nh
      r=DFLOAT(nh)/DFLOAT(l)
      DO i=l,1,-1
         IF (j.GT.0.AND.j.LT.i.AND.RAN().LT.r) THEN
            str(i)=str(j)
            str(j)=0
            j=j-1
         ENDIF
      ENDDO
            
      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE INITCONF
c------------------------------------------------------------------------------
c     Initializes a configuration with random spin state and empty string.
c------------------------------------------------------------------------------

      INCLUDE 'ahb1.h'

      INTEGER i
      REAL*8  RAN

      l=20
      nh=0
      DO i=1,nn
         spn(i)=2*MIN(INT(RAN()*2.d0),1)-1
      ENDDO
      DO i=1,ll
         str(i)=0
      ENDDO

      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE LATTICE
c------------------------------------------------------------------------------
c     Constructs the list of sites bst(1,b) and bst(2,b) connected by bond b.
c     The x-bonds correspond to b=1,...,N, y-bonds b=N+1,...,2*N, and
c     z-bonds b=2*N+1,...,3*N, where N=nx**dd is the number of sites. Also 
c     constructs the staggered phase array phs().
c------------------------------------------------------------------------------

      INCLUDE 'ahb1.h'

      INTEGER is,x1,x2,y1,y2,z1,z2

      is=0
      IF (dd.EQ.1) THEN
         DO x1=1,nx
            is=is+1
            x2=1+MOD(x1,nx)
            bst(1,is)=is
            bst(2,is)=x2
         ENDDO
      ELSEIF (dd.EQ.2) THEN
        DO y1=1,nx
        DO x1=1,nx
            is=is+1
            x2=1+MOD(x1,nx)
            y2=y1
            bst(1,is)=is
            bst(2,is)=x2+(y2-1)*nx
            x2=x1
            y2=1+MOD(y1,nx)
            bst(1,is+nn)=is
            bst(2,is+nn)=x2+(y2-1)*nx
         ENDDO
         ENDDO
      ELSEIF (dd.EQ.3) THEN
        DO z1=1,nx
        DO y1=1,nx
        DO x1=1,nx
            is=is+1
            x2=1+MOD(x1,nx)
            y2=y1
            z2=z1
            bst(1,is)=is
            bst(2,is)=x2+(y2-1)*nx+(z2-1)*nx*nx
            x2=x1
            y2=1+MOD(y1,nx)
            z2=z1
            bst(1,is+nn)=is
            bst(2,is+nn)=x2+(y2-1)*nx+(z2-1)*nx*nx
            x2=x1
            y2=y1
            z2=1+MOD(z1,nx)
            bst(1,is+2*nn)=is
            bst(2,is+2*nn)=x2+(y2-1)*nx+(z2-1)*nx*nx
         ENDDO
         ENDDO
         ENDDO
      ENDIF
      DO is=1,nn
         x1=MOD(is-1,nx)
         y1=MOD(is-1,nx*nx)/nx
         z1=(is-1)/(nx*nx)
         phs(is)=(-1)**(x1+y1+z1)
      ENDDO

      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE READCONF
c------------------------------------------------------------------------------
c     Reads on old configuration stored in the file "conf".
c------------------------------------------------------------------------------

      INCLUDE 'ahb1.h'

      INTEGER i

      OPEN(UNIT=10,FILE='conf',STATUS='old')
      READ(20,*)l,nh
      DO i=1,nn
         READ(10,*)spn(i)
      ENDDO
      DO i=1,l
         READ(10,*)str(i)
      ENDDO
      CLOSE(10) 

      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE WRITECONF
c------------------------------------------------------------------------------
c     Writes the current configuration to the file "conf". 
c------------------------------------------------------------------------------

      INCLUDE 'ahb1.h'

      INTEGER i

      OPEN(UNIT=10,FILE='conf',STATUS='unknown')
      WRITE(10,*)l,nh
      DO i=1,nn
         WRITE(10,*)spn(i)
      ENDDO
      DO i=1,l
         WRITE(10,*)str(i)
      ENDDO
      CLOSE(10)

      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE OPENLOG
c------------------------------------------------------------------------------
c     Opens the log-file.
c------------------------------------------------------------------------------

      INTEGER i

      OPEN(UNIT=12,FILE='log.txt',STATUS='unknown',ACCESS='append')
      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE CLOSELOG
c------------------------------------------------------------------------------
c     Closes the log-file.
c------------------------------------------------------------------------------

      CLOSE(12)
      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE INITRAN
c------------------------------------------------------------------------------

      IMPLICIT NONE

      INTEGER iir,jjr,kkr,nnr,mzran
      COMMON/BRANDOM/iir,jjr,kkr,nnr,mzran

      INTEGER seed(4)

      OPEN(UNIT=10,FILE='rand.in',STATUS='old')
      READ(10,*)seed(1)
      READ(10,*)seed(2)
      READ(10,*)seed(3)
      READ(10,*)seed(4)
      CLOSE(10)      
      iir=1+iabs(seed(1))
      jjr=1+iabs(seed(2))
      kkr=1+iabs(seed(3))
      nnr=seed(4)

      RETURN
      END

c------------------------------------------------------------------------------
      REAL*8 FUNCTION RAN()
c------------------------------------------------------------------------------
c     Mixed generator MZRAN by G. A. Marsaglia.
c------------------------------------------------------------------------------

      IMPLICIT NONE

      INTEGER iir,jjr,kkr,nnr,mzran
      COMMON/BRANDOM/iir,jjr,kkr,nnr,mzran

      mzran=iir-kkr
      IF (mzran.LT.0) mzran=mzran+2147483579
      iir=jjr
      jjr=kkr
      kkr=mzran
      nnr=69069*nnr+1013904243
      mzran=mzran+nnr
      RAN=0.5d0+0.23283064d-9*mzran

      RETURN
      END

c------------------------------------------------------------------------------