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 ', + ' 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 ============================================================