c     =======================================================================
c     PROGRAM: The autocorrelation test for pseudorandom numbers using the
c              2-d Ising model with the Wolff algorithm.
c     -----------------------------------------------------------------------
c     for    testing of pseudorandom numbers; based on calculating some
c            physical quantities of the 2-d Ising model as well as their
c            integrated autocorrelation functions.
c     by     I. Vattulainen, vattulai@convex.csc.fi
c     alg    autocorrelation, the Ising model, Wolff updating
c     ref    Phys. Rev. Lett. 73, 2513 (1994).
c     title  acorrt_iw.f
c     size
c     prec   single/double
c     lang   Fortran77
c     -----------------------------------------------------------------------
c     The author of this software is I. Vattulainen.  Copyright (c) 1993.
c     Permission to use, copy, modify, and distribute this software for
c     any purpose without fee is hereby granted, provided that this entire
c     notice is included in all copies of any software which is or includes
c     a copy or modification of this software and in all copies of the
c     supporting documentation for such software.
c     THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
c     WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKE ANY
c     REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
c     OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
c     -----------------------------------------------------------------------
c     In the autocorrelation test, we calculate averages of some physical
c     quantities such as the energy, the susceptibility, and the updated
c     cluster size in the 2-d Ising model. Additionally, their autocorrelation
c     functions and corresponding integrated autocorrelation values are
c     determined. The chosen temperature corresponds to the critical point,
c     and the Wolff updating method is utilized. The results of the
c     simulation (test) are a function of the pseudorandom number generator
c     used. Since the exact value is known only for the energy, other results
c     must be compared with corresponding results of other generators;
c     i.e. the test must be performed comparatively between several
c     pseudorandom number generators.
c     -----------------------------------------------------------------------
c     Main parameters in the test are as follows:
c     L         linear size of the lattice.
c     DK        interaction / temperature*Boltzmann constant.
c     MCBEG     number of passed initial configurations
c               (initialization time; a minimum of 10000 is suggested).
c     ISCANS    number of independent samples.
c     ICSAMP    maximum detected autocorrelation time.
c     INEWEST   tells the index of the latest (newest) configuration.
c     IZ        contains the 2-d lattice.
c     D*ROO     contains the autocorrelation function.
c     D*TAU     gives the integrated autocorrelation time, where
c     D*SIGMA   is its' error estimate.
c     ------------------------------------------------------
c     ******************
c     PRACTICAL DETAILS:
c     ******************
c
c     Required modification to this code:
c                   Initialization and implementation of the pseudorandom
c                   number generator (PRNG), which will be tested. For
c                   initialization, there is a reserved space in the code
c                   below. The PRNG may be appended to the end of the file,
c                   or may be called separately during compiling. In this
c                   version, the default generator is GGL.
c
c                   During the test, random numbers are called in two cases:
c                   in the formation of the random (initial) lattice and
c                   during its updating. The former is in the main program
c                   and the latter in the subroutine WOLFF. As a consequence,
c                   if you are testing a generator which produces a sequence
c                   of random numbers at a time (this is not the case with
c                   GGL), the code must be modified correspondingly.
c
c     Output files: general.dat         General information of the test.
c                   susc_data.dat       Autocorrelation function and its
c                                       integrated value for susceptibility
c                                       (including the error estimate).
c                   ener_data.dat       Corresponding results for the energy.
c                   clus_data.dat       Corresponding results for the
c                                       flipped cluster sizes.
c                   clus_distr.dat      Normalized probability distribution
c                                       for the flipped cluster size (as a
c                                       function of cluster size).
c
c     Keep in mind that the given results are *NOT* final. This is due to
c     the calculation of errors in a Monte Carlo simulation, in which the
c     integrated autocorrelation time plays an important role. In order to
c     find the correct error estimates, perform the following steps after
c     the test:
c        (1) First of all, determine the estimates for the integrated
c            autocorrelation times. To do that, for each of the quantities
c            energy, susceptibility, and flipped cluster size, you must find
c            the regime (plateau) in which the results have already converged
c            but where the statistical errors are not yet very large. An
c            estimate for the corresponding integrated autocorrelation time
c            (tau) is then an average over this plateau.
c        (2) The given values for the error estimates of energy,
c            susceptibility and the average flipped cluster size (given in
c            a file general.dat) must then be modified: when the
c            corresponding integrated autocorrelation time (tau) is known,
c            the error estimate must then be multiplied by tau. A square
c            root of the product then gives the correct error estimate.
c
c     For further practical details of the method see U. Wolff [Phys. Lett. B
c     228, 379 (1989) or I. Vattulainen [cond-mat@babbage.sissa.it No.
c     9411062].
c     =======================================================================

      IMPLICIT NONE
      DOUBLE PRECISION
     +     DK
      INTEGER
     +     L, L2, MCBEG, ISCANS, ICSAMP
      PARAMETER( L = 16, L2 = L*L )
c     -------------------------------------
c     Coefficients at criticality (some empirical estimates for a finite
c     system are also given):
c     For an infinite lattice (L --> infinity):
c      PARAMETER( DK = 0.44068679350977D0 )
c     For a lattice size L = 64:
c      PARAMETER( DK = 0.43821982133122D0 )
c     For a lattice size L = 16:
c      PARAMETER( DK = 0.43098188945039D0 )
c     -------------------------------------
      PARAMETER( DK = 0.44068679350977D0 )
c     Please choose MCBEG > 9999 (time for equilibration).
      PARAMETER( MCBEG = 10000 )
c     Good statistics require ISCANS > 1e6.
      PARAMETER( ISCANS = 10000000 )
      PARAMETER( ICSAMP = 25 )

c     ===============================================================

      INTEGER
     +     IZ(L,L), IN1(L), IP1(L), ISTACK(L2,2),
     +     IX, IY, ITRSUM, ISUM, IPOS(ICSAMP),
     +     INEWEST, IDT, ILOCAL, ICNT, IFIN(2),
     +     IFAC, ICLSIZ, ICLDST(L2), I, J,
     +     ITRIES, IS, IDIFF, IKS

      DOUBLE PRECISION
     +     PROB,
c     Susceptibility variables:
     +     DSAUTO(ICSAMP), DSUS, DSUSOLD(ICSAMP), DSUSAVE, DSUS2,
     +     DSROO(0:ICSAMP), DSKAPPA(ICSAMP), DSRW(ICSAMP),
     +     DSTAU(ICSAMP), DSXW(ICSAMP), DSYW(ICSAMP),
     +     DSSIGMA(ICSAMP), DSAVERAGE, DSERROR,
c     Energy variables:
     +     DEAUTO(ICSAMP), DENE, DENEOLD(ICSAMP), DENEAVE, DENE2,
     +     DEROO(0:ICSAMP), DEKAPPA(ICSAMP), DERW(ICSAMP),
     +     DETAU(ICSAMP), DEXW(ICSAMP), DEYW(ICSAMP),
     +     DESIGMA(ICSAMP), DEAVERAGE, DEERROR,
c     Cluster variables:
     +     DCLAUTO(ICSAMP), DCL, DCLOLD(ICSAMP), DCLAVE, DCL2,
     +     DCLROO(0:ICSAMP), DCLKAPPA(ICSAMP), DCLRW(ICSAMP),
     +     DCLTAU(ICSAMP), DCLXW(ICSAMP), DCLYW(ICSAMP),
     +     DCLSIGMA(ICSAMP), DCLAVR, DCLERROR, DCLRAT

      REAL*8
     +     DSEED, DISEED

      REAL
     +     RAN, GGL

      OPEN(20,STATUS='NEW',FILE='susc_data.dat',FORM='FORMATTED',
     +     ACCESS='SEQUENTIAL')
      OPEN(21,STATUS='NEW',FILE='ener_data.dat',FORM='FORMATTED',
     +     ACCESS='SEQUENTIAL')
      OPEN(22,STATUS='NEW',FILE='general.dat',FORM='FORMATTED',
     +     ACCESS='SEQUENTIAL')
      OPEN(23,STATUS='NEW',FILE='clus_data.dat',FORM='FORMATTED',
     +     ACCESS='SEQUENTIAL')
      OPEN(24,STATUS='NEW',FILE='clus_distr.dat',FORM='FORMATTED',
     +     ACCESS='SEQUENTIAL')

c     --------------------------
c     Initialize some variables:
c     --------------------------
      ISUM = 0
      ITRSUM = 0
c
      DSUS = 0.D0
      DSUSAVE = 0.D0
      DSUS2 = 0.D0
c
      DENE = 0.D0
      DENEAVE = 0.D0
      DENE2 = 0.D0
c
      DCL = 0.D0
      DCLAVE = 0.D0
      DCLRAT = 0.D0
      DCL2 = 0.D0
c
      DO 801 I=1,ICSAMP
         DSUSOLD(I) = 0.D0
         DENEOLD(I) = 0.D0
         DCLOLD(I) = 0.D0
 801  CONTINUE
      DO 13 I=1,ICSAMP
         DSAUTO(I) = 0.D0
         DEAUTO(I) = 0.D0
         DCLAUTO(I) = 0.D0
 13   CONTINUE
      DO 12 I=1,L
         IN1(I) = I - 1
         IP1(I) = I + 1
 12   CONTINUE
      IN1(1) = L
      IP1(L) = 1
      DO 949 I=1,L2
         ICLDST(I) = 0
 949  CONTINUE

c     ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
c     First of all, initialize the pseudorandom number generator.
c     In the case of GGL, only one number is needed:

      DSEED   = 14159.D0
c     ```````````````````````````````````````````````````````````
      DISEED = DSEED

c     ----------------------------------------------------------
c     Probability PROB to unite equal lattice sites as clusters:
c     ----------------------------------------------------------
      PROB = 1.D0 - DEXP(-2.D0*DK)

c     --------------------------------------
c     Create a random initial configuration:
c     --------------------------------------
      DO 15 IY=1,L
         DO 14 IX=1,L
            RAN = GGL(DSEED)
            IF(RAN.LT.0.5D0) THEN
               IZ(IY,IX) = -1
            ELSE
               IZ(IY,IX) = 1
            ENDIF
 14      CONTINUE
 15   CONTINUE

c     ------------------------------------------------------------------
c     Perform the dynamic part to achieve equilibrium.
c     ------------------------------------------------------------------
c     Pass a total of MCBEG configurations. Determine MCBEG by studying
c     the order parameter of the system (as function of T), for example.
c     ------------------------------------------------------------------

      DO 18 I=1,MCBEG
         ITRIES = 0
 711     CALL WOLFF(L,L2,IZ,ISTACK,IN1,IP1,PROB,
     +        DSEED,ISUM,ICLSIZ)
         ITRIES = ITRIES + ISUM
         IF(ITRIES.LT.L2) GOTO 711
 18   CONTINUE

c     ----------------------------------------------
c     Take a number of ICSAMP samples to start with:
c     ----------------------------------------------
      DO 19 I=1,ICSAMP
         CALL WOLFF(L,L2,IZ,ISTACK,IN1,IP1,PROB,
     +        DSEED,ISUM,ICLSIZ)
         CALL SUSCEP(L,IZ,DSUS)
         DSUSOLD(I) = DSUS
         CALL ENERG(L,L2,IZ,DENE,IN1,IP1)
         DENEOLD(I) = DENE
         DCLOLD(I) = DFLOAT(ICLSIZ)
 19   CONTINUE

c     -------------------------------------------------------------------
c     Set the initial state for the vector IPOS, which gives the order of
c     appearance of the configurations. IPOS(1) gives the newest, IPOS(2)
c     the second newest, ..., and IPOS(IC) the oldest configuration.
c     -------------------------------------------------------------------
      INEWEST = ICSAMP
      DO 166 I=1,ICSAMP
         IPOS(I) = ICSAMP - I + 1
 166  CONTINUE

C     =================================================
c     !!!!!!!!!!!!!! MAIN PROGRAM STARTS !!!!!!!!!!!!!!
c     =================================================
c     Take a number of ISCANS independent measurements.
C     Each measurement corresponds to flipping of a
c     single cluster.
c     =================================================

      DO 20 ICNT=1,ISCANS

         CALL WOLFF(L,L2,IZ,ISTACK,IN1,IP1,PROB,
     +        DSEED,ISUM,ICLSIZ)
c     -----------------------------------------------------
c     Calculate the susceptibility autocorrelation function
c     corresponding to the latest sample.
c     -----------------------------------------------------
         CALL SUSCEP(L,IZ,DSUS)
         DSUSAVE = DSUSAVE + DSUS
         DSUS2 = DSUS2 + DSUS*DSUS
         DO 48 IS=1,ICSAMP
            DSAUTO(IS) = DSAUTO(IS) + DSUSOLD(IPOS(IS))*DSUS
 48      CONTINUE
c     ---------------------------------------------
c     Calculate the energy autocorrelation function
c     corresponding to the latest sample.
c     ---------------------------------------------
         CALL ENERG(L,L2,IZ,DENE,IN1,IP1)
         DENEAVE = DENEAVE + DENE
         DENE2 = DENE2 + DENE*DENE
         DO 78 IS=1,ICSAMP
            DEAUTO(IS) = DEAUTO(IS) + DENEOLD(IPOS(IS))*DENE
 78      CONTINUE
c     ---------------------------------------------------
c     Then calculate the distribution of updated clusters
c     and their autocorrelation function.
c     ---------------------------------------------------
         DCL = DFLOAT(ICLSIZ)
         DCLAVE = DCLAVE + DCL
         DCL2 = DCL2 + DCL*DCL
         DO 98 IS=1,ICSAMP
            DCLAUTO(IS) = DCLAUTO(IS) + DCLOLD(IPOS(IS))*DCL
 98      CONTINUE
         ICLDST(ICLSIZ) = ICLDST(ICLSIZ) + 1
c     -----------------------------------------------------------------
c     Replace the latest lattice with the second latest and so forth...
c     -----------------------------------------------------------------
         INEWEST = INEWEST + 1
         IF(INEWEST.GT.ICSAMP) INEWEST = 1
         DSUSOLD(INEWEST) = DSUS
         DENEOLD(INEWEST) = DENE
         DCLOLD(INEWEST) = DCL
c     -----------------------------------
c     Update the order of configurations.
c     -----------------------------------
         DO 167 I=1,ICSAMP
            IPOS(I) = IPOS(I) + 1
            IF(IPOS(I).GT.ICSAMP) IPOS(I) = 1
 167     CONTINUE

 20   CONTINUE

c     ==================
c     End of sampling...
c     ==================
C     ===============================================
c     !!!!!!!!!!!!!! MAIN PROGRAM ENDS !!!!!!!!!!!!!!
C     ===============================================

c     -----------------------------------------------
c     Susceptibility case:
c     Calculate the estimator for the integrated
c     autocorrelation time. In addition, determine
c     the error estimate. For further details see
c     [U. Wolff, Phys. Lett. B {\bf 228}, 379 (1989).
c     -----------------------------------------------

      DSUSAVE = DSUSAVE/ISCANS
c     Average of samples:
      DSAVERAGE = DSUSAVE
      DSUSAVE = DSUSAVE*DSUSAVE
c     Average of squares of samples:
      DSUS2 = DSUS2/ISCANS
c     Monte Carlo error estimate squared divided by \tau:
      DSERROR = (DSUS2 - DSUSAVE)*2.D0/DFLOAT(ISCANS)

      DO 331 I=1,ICSAMP
         DSAUTO(I) = DSAUTO(I)/ISCANS
 331  CONTINUE
      DSROO(0) = 1.D0
      DO 332 I=1,ICSAMP
         DSROO(I) = (DSAUTO(I) - DSUSAVE)/(DSUS2 - DSUSAVE)
 332  CONTINUE
      DO 333 I=1,ICSAMP
         DSKAPPA(I) = DSROO(I)/(DSROO(I-1))
 333  CONTINUE
      DO 334 I=1,ICSAMP
         DSRW(I) = DSROO(I)/(1.D0 - DSKAPPA(I))
 334  CONTINUE
      DO 336 I=1,ICSAMP
         DSTAU(I) = 0.5D0 + DSRW(I)
         IF(I.GE.2) THEN
            DO 335 J=1,(I-1)
               DSTAU(I) = DSTAU(I) + DSROO(J)
 335        CONTINUE
         ENDIF
 336  CONTINUE
c     Error:
      DO 338 I=1,ICSAMP
         DSXW(I) = 0.5D0
         DSYW(I) = 0.5D0
         IF(I.GE.2) THEN
            DO 337 J=1,(I-1)
               DSXW(I) = DSXW(I) + DSROO(J)*DSROO(J)
               DSYW(I) = DSYW(I) +
     +              (DSROO(J) - DSROO(J-1))*(DSROO(J) - DSROO(J-1))
 337        CONTINUE
         ENDIF
 338  CONTINUE
      DO 339 I=1,ICSAMP
         DSSIGMA(I) =
     +        DSTAU(I)*DSTAU(I)*(DFLOAT(I) - 0.5D0 +
     +        (1.D0 + DSKAPPA(I))/(1.D0 - DSKAPPA(I))) +
     +        DSKAPPA(I)*DSXW(I)/((1.D0 - DSKAPPA(I))**2.D0) +
     +        (DSKAPPA(I)**2.D0)*DSYW(I)/((1.D0 - DSKAPPA(I))**4.D0)
 339  CONTINUE
      DO 340 I=1,ICSAMP
         DSSIGMA(I) = DSSIGMA(I)*4.D0/DFLOAT(ISCANS)
 340  CONTINUE
      DO 341 I=1,ICSAMP
         DSSIGMA(I) = DSQRT(DSSIGMA(I))
 341  CONTINUE

c     -----------------------------------------------
c     Energy case:
c     Calculate the estimator for the integrated
c     autocorrelation time. In addition, determine
c     the error estimate. For further details see
c     [U. Wolff, Phys. Lett. B {\bf 228}, 379 (1989).
c     -----------------------------------------------

      DENEAVE = DENEAVE/ISCANS
      DEAVERAGE = DENEAVE
      DENEAVE = DENEAVE*DENEAVE
      DENE2 = DENE2/ISCANS
      DEERROR = (DENE2 - DENEAVE)*2.D0/DFLOAT(ISCANS)

      DO 731 I=1,ICSAMP
         DEAUTO(I) = DEAUTO(I)/ISCANS
 731  CONTINUE
      DEROO(0) = 1.D0
      DO 732 I=1,ICSAMP
         DEROO(I) = (DEAUTO(I) - DENEAVE)/(DENE2 - DENEAVE)
 732  CONTINUE
      DO 733 I=1,ICSAMP
         DEKAPPA(I) = DEROO(I)/(DEROO(I-1))
 733  CONTINUE
      DO 734 I=1,ICSAMP
         DERW(I) = DEROO(I)/(1.D0 - DEKAPPA(I))
 734  CONTINUE
      DO 736 I=1,ICSAMP
         DETAU(I) = 0.5D0 + DERW(I)
         IF(I.GE.2) THEN
            DO 735 J=1,(I-1)
               DETAU(I) = DETAU(I) + DEROO(J)
 735        CONTINUE
         ENDIF
 736  CONTINUE
c     Error:
      DO 738 I=1,ICSAMP
         DEXW(I) = 0.5D0
         DEYW(I) = 0.5D0
         IF(I.GE.2) THEN
            DO 737 J=1,(I-1)
               DEXW(I) = DEXW(I) + DEROO(J)*DEROO(J)
               DEYW(I) = DEYW(I) +
     +              (DEROO(J) - DEROO(J-1))*(DEROO(J) - DEROO(J-1))
 737        CONTINUE
         ENDIF
 738  CONTINUE
      DO 739 I=1,ICSAMP
         DESIGMA(I) =
     +        DETAU(I)*DETAU(I)*(DFLOAT(I) - 0.5D0 +
     +        (1.D0 + DEKAPPA(I))/(1.D0 - DEKAPPA(I))) +
     +        DEKAPPA(I)*DEXW(I)/((1.D0 - DEKAPPA(I))**2.D0) +
     +        (DEKAPPA(I)**2.D0)*DEYW(I)/((1.D0 - DEKAPPA(I))**4.D0)
 739  CONTINUE
      DO 740 I=1,ICSAMP
         DESIGMA(I) = DESIGMA(I)*4.D0/DFLOAT(ISCANS)
 740  CONTINUE
      DO 741 I=1,ICSAMP
         DESIGMA(I) = DSQRT(DESIGMA(I))
 741  CONTINUE

c     -------------
c     Cluster case:
c     -------------------------------------------------------
      DCLAVE = DCLAVE/ISCANS
      DCLAVR = DCLAVE
      DCLRAT = DCLAVR/L2
      DCLAVE = DCLAVE*DCLAVE
      DCL2 = DCL2/ISCANS
      DCLERROR = (DCL2 - DCLAVE)*2.D0/DFLOAT(ISCANS)

      DO 931 I=1,ICSAMP
         DCLAUTO(I) = DCLAUTO(I)/ISCANS
 931  CONTINUE
      DCLROO(0) = 1.D0
      DO 932 I=1,ICSAMP
         DCLROO(I) = (DCLAUTO(I) - DCLAVE)/(DCL2 - DCLAVE)
 932  CONTINUE

      DO 933 I=1,ICSAMP
         DCLKAPPA(I) = DCLROO(I)/(DCLROO(I-1))
 933  CONTINUE
      DO 934 I=1,ICSAMP
         DCLRW(I) = DCLROO(I)/(1.D0 - DCLKAPPA(I))
 934  CONTINUE
      DO 936 I=1,ICSAMP
         DCLTAU(I) = 0.5D0 + DCLRW(I)
         IF(I.GE.2) THEN
            DO 935 J=1,(I-1)
               DCLTAU(I) = DCLTAU(I) + DCLROO(J)
 935        CONTINUE
         ENDIF
 936  CONTINUE
c     Error:
      DO 938 I=1,ICSAMP
         DCLXW(I) = 0.5D0
         DCLYW(I) = 0.5D0
         IF(I.GE.2) THEN
            DO 937 J=1,(I-1)
               DCLXW(I) = DCLXW(I) + DCLROO(J)*DCLROO(J)
               DCLYW(I) = DCLYW(I) +
     +              (DCLROO(J) - DCLROO(J-1))*
     +              (DCLROO(J) - DCLROO(J-1))
 937        CONTINUE
         ENDIF
 938  CONTINUE
      DO 939 I=1,ICSAMP
         DCLSIGMA(I) =
     +        DCLTAU(I)*DCLTAU(I)*(DFLOAT(I) - 0.5D0 +
     +        (1.D0 + DCLKAPPA(I))/(1.D0 - DCLKAPPA(I))) +
     +        DCLKAPPA(I)*DCLXW(I)/((1.D0 - DCLKAPPA(I))**2.D0) +
     +        (DCLKAPPA(I)**2.D0)*DCLYW(I)/((1.D0 - DCLKAPPA(I))**4.D0)
 939  CONTINUE
      DO 940 I=1,ICSAMP
         DCLSIGMA(I) = DCLSIGMA(I)*4.D0/DFLOAT(ISCANS)
 940  CONTINUE
      DO 941 I=1,ICSAMP
         DCLSIGMA(I) = DSQRT(DCLSIGMA(I))
 941  CONTINUE

c     ------------------
c     Write the results:
c     ------------------
c     General data:
      WRITE(22,*) ' ========================================',
     ,     '=========================='
      WRITE(22,*) ' Parameters in the autocorrelation test: '
      WRITE(22,*) '      Initial seed:                      ',DISEED
      WRITE(22,*) '      Lattice size (L):                  ',L
      WRITE(22,*) '      Number of samples (ISCANS):        ',ISCANS
      WRITE(22,*) '      Coupling constant (DK):            ',DK
      WRITE(22,*) '      Number of omitted confs. (MCBEG):  ',MCBEG
      WRITE(22,*)
      WRITE(22,*) ' Average energy of the system:           ',
     +     DEAVERAGE
      WRITE(22,*) ' and its squared (error estimate / tau): ',
     +     DEERROR
      WRITE(22,*)
      WRITE(22,*) ' Average susceptibility of the system:   ',
     +     DSAVERAGE/DFLOAT(L2)
      WRITE(22,*) ' and its squared (error estimate / tau): ',
     +     DSERROR/DFLOAT(L2)
      WRITE(22,*)
c     The average cluster size is given here in a non-normalized form.
c     (not divided by the system size L2).
      WRITE(22,*) ' Average flipped cluster size:           ',
     +     DCLAVR
      WRITE(22,*) ' and its squared (error estimate / tau): ',
     +     DCLERROR
      WRITE(22,*) ' ----------------------------------------',
     ,     '--------------------------'
      WRITE(22,*) ' The notation <squared (error estimate / tau)>',
     +     ' means the following:'
      WRITE(22,*) ' to get the correct error estimate, multiply',
     +     ' the given value by'
      WRITE(22,*) ' the integrated autocorrelation time tau (see',
     +     ' comments in the code '
      WRITE(22,*) ' for details), and then take a square root.'
      WRITE(22,*) ' ========================================',
     ,     '=========================='
      CLOSE(22)
c     Susceptibility results:
      WRITE(20,*) '      ==================================='
      WRITE(20,*) '      AUTO-CORRELATION OF SUSCEPTIBILITY:'
      WRITE(20,*) '      ==================================='
      IFIN(1) = 0
      IFIN(2) = 1
      WRITE(20,*) IFIN(1), DFLOAT(IFIN(2))
      DO 350 I=1,ICSAMP
         WRITE(20,*) I, DSROO(I)
 350  CONTINUE
      WRITE(20,*) '      ====================================='
      WRITE(20,*) '      INTEGRATED AC-TIME OF SUSCEPTIBILITY:'
      WRITE(20,*) '      ====================================='
      DO 351 I=1,ICSAMP
         WRITE(20,*) I,
     +        DSTAU(I)*DCLRAT,
     +        DSSIGMA(I)*DCLRAT
 351  CONTINUE
      CLOSE(20)
c     Energy results:
      WRITE(21,*) '      ==========================='
      WRITE(21,*) '      AUTO-CORRELATION OF ENERGY:'
      WRITE(21,*) '      ==========================='
      WRITE(21,*) IFIN(1), DFLOAT(IFIN(2))
      DO 750 I=1,ICSAMP
         WRITE(21,*) I, DEROO(I)
 750  CONTINUE
      WRITE(21,*) '      ============================='
      WRITE(21,*) '      INTEGRATED AC-TIME OF ENERGY:'
      WRITE(21,*) '      ============================='
      DO 751 I=1,ICSAMP
         WRITE(21,*) I,
     +        DETAU(I)*DCLRAT,
     +        DESIGMA(I)*DCLRAT
 751  CONTINUE
      CLOSE(21)
c     Cluster results:
      WRITE(23,*) '      ====================================='
      WRITE(23,*) '      AUTO-CORRELATION OF FLIPPED CLUSTERS:'
      WRITE(23,*) '      ====================================='
      WRITE(23,*) IFIN(1), DFLOAT(IFIN(2))
      DO 950 I=1,ICSAMP
         WRITE(23,*) I, DCLROO(I)
 950  CONTINUE
      WRITE(23,*) '      ======================================='
      WRITE(23,*) '      INTEGRATED AC-TIME OF FLIPPED CLUSTERS:'
      WRITE(23,*) '      ======================================='
      DO 951 I=1,ICSAMP
         WRITE(23,*) I,
     +        DCLTAU(I)*DCLRAT,
     +        DCLSIGMA(I)*DCLRAT
 951  CONTINUE
      CLOSE(23)

c     Distribution of flipped clusters:
      DO 432 I=1,L2
         IKS = ICLDST(I)
         IF(IKS.NE.0) THEN
            WRITE(24,*) I, DFLOAT(ICLDST(I))/ISCANS
         ENDIF
 432  CONTINUE
      CLOSE(24)

      STOP
      END

c     ============================================================
c     Wolff Monte Carlo algorithm (discrete case). The idea is by
c     U. Wolff [Phys. Rev. Lett. 60 (1988) 1461] and the program
c     is a modification of Ref. [Physica A 167 (1990) 565] by
c     J. S. Wang and R. H. Swendsen.
c     ------------------------------
c     Ilpo Vattulainen, Tampere univ. of tech., Funland.
c     14.4.1993
c     ============================================================

      SUBROUTINE WOLFF(L,L2,IZ,ISTACK,IN1,IP1,PROB,
     +     DSEED,ISUM,ICLSIZ)

      INTEGER
     +     L, L2, IZ(L,L), ISTACK(L2,2), IN1(L), IP1(L), NN(4,2),
     +     ISIZE, I, J, IPT, KY, KX, ISUM, ISHERE, NNY, NNX,
     +     ICLSIZ
      DOUBLE PRECISION
     +     PROB
      REAL*8
     +     DSEED
      REAL
     +     RAN, RANS, RANX, RANY

      IPT = 0
      ISUM = 0

c     ------------------------------------------------------------------
c     Start a Monte Carlo loop to update a total of MC*L2 lattice sites.
c     ------------------------------------------------------------------

c     ISIZE calculates the number of tried updates.
c     ICLSIZ calculates the size of the flipped cluster.

      ISIZE = 1
      ICLSIZ = 1

c     ------------------------------------------------------------------------
c     Choose a starting position and change its state immediately. Other parts
c     of a cluster will also be changed immediately after recognition.
c     ------------------------------------------------------------------------
      RANY = GGL(DSEED)
      RANX = GGL(DSEED)
      KY = RANY*L + 1
      IF(KY.GT.L) KY = L
      KX = RANX*L + 1
      IF(KX.GT.L) KX = L
      ISHERE = IZ(KY,KX)
      IZ(KY,KX) = -ISHERE
c     -----------------------------
c     Periodic boundary conditions:
c     -----------------------------
 120  NN(1,1) = KY
      NN(1,2) = IN1(KX)
      NN(2,1) = KY
      NN(2,2) = IP1(KX)
      NN(3,1) = IN1(KY)
      NN(3,2) = KX
      NN(4,1) = IP1(KY)
      NN(4,2) = KX
c     ------------------------
c     Check nearest neighbors:
c     ------------------------
      DO 125 J=1,4
         NNY = NN(J,1)
         NNX = NN(J,2)
         IF(ISHERE.NE.IZ(NNY,NNX)) GOTO 125
         ISIZE = ISIZE + 1
         RANS = GGL(DSEED)
         IF(RANS.GT.PROB) GOTO 125
         IZ(NNY,NNX) = -IZ(NNY,NNX)
         ICLSIZ = ICLSIZ + 1
c     -----------------------------------------------------------------
c     Increment the stack of lattice sites which are part of a cluster,
c     but whose nearest neighbors have not yet been checked:
c     -----------------------------------------------------------------
         IPT = IPT + 1
         ISTACK(IPT,1) = NNY
         ISTACK(IPT,2) = NNX
 125  CONTINUE
c     --------------------------------------------------------------
c     When IPT equals zero the cluster has reached its maximum size.
c     Otherwise go back and check the remaining possibilities.
c     --------------------------------------------------------------
      IF(IPT.EQ.0) GOTO 140
      KY = ISTACK(IPT,1)
      KX = ISTACK(IPT,2)
      IPT = IPT - 1
      GOTO 120
c     --------------------------------------------
c     Calculate the total amount of tried updates:
c     --------------------------------------------
 140  ISUM = ISIZE

      RETURN
      END

c     =============================================================
c
c     Calculate the susceptibility based on the definition of Wolff
c     [U. Wolff, Phys. Lett. B {\bf 228}, 379 (1989).
c
c     =============================================================
      SUBROUTINE SUSCEP(L,IZ,DSUS)

      INTEGER L, IZ(L,L), IZNUM
      DOUBLE PRECISION DSUS

      IZNUM = 0
      DO 299 IY=1,L
         DO 298 IX=1,L
            IZNUM = IZNUM + IZ(IY,IX)
 298     CONTINUE
 299  CONTINUE
      DSUS = DFLOAT(IZNUM)
      DSUS = DSUS*DSUS
      DSUS = DSUS/(L*L)

      RETURN
      END

c     ==============================================================
c
c     Calculate the energy of the system (in variables of U/J).
c
c     ==============================================================

      SUBROUTINE ENERG(L,L2,IZ,DENE,IN1,IP1)

      INTEGER L, L2,IZ(L,L), IN1(L), IP1(L), IZNUM
      DOUBLE PRECISION DENE

      IZNUM = 0
      DO 299 IY=1,L
         DO 298 IX=1,L
            IZNUM = IZNUM + IZ(IY,IX)*(IZ(IN1(IY),IX) +
     +           IZ(IY,IN1(IX)) + IZ(IP1(IY),IX) + IZ(IY,IP1(IX)))
 298     CONTINUE
 299  CONTINUE
      DENE = DFLOAT(IZNUM)/(2*L2)

      RETURN
      END

c     ==============================================================
C     A pseudorandom number generator GGL
C     ------------------------------------------------------------

      REAL FUNCTION GGL (DS)

      DOUBLE PRECISION
c     +     D1,
     +     DS, D2
c      DATA D1/2147483648.D0/
      DATA D2/2147483647.D0/
      DS  = DMOD(16807.D0*DS,D2)
c     Generate U(0,1] distributed random numbers:
      GGL = DS/D2

      RETURN
      END

C     ============================================================