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 ============================================================
.
