c c*********************************************************************** c c program : 2dnew c c 2d ising model c c algorithm : swendsen and wang based on an idea by hu c c version : 0.1 as of september c c dieter w. heermann c c*********************************************************************** c integer spin(0:30,0:30) integer lbond(1:30,1:30),ubond(1:30,1:30) integer stack(1:5000) integer nsd(1:1000) integer nsc(1:1000) integer nsa(1:1000) integer nsa1(1:1000) c real jkt real mbond,sbond c c-----set and/or read the input parameter c write(*,*) 'lattice size = ' read(*,*) l write(*,*) 'temp = ' read(*,*) ttc write(*,*) 'mcsmax = ' read(*,*) mcsmax write(*,*) 'iskip = ' read(*,*) iskip write(*,*) 'iseed = ' read(*,*) iseed c write(*,*) write(*,*) 'program 2dnew2' write(*,*) write(9,*) 'simulation parameters' write(9,*) '---------------------' write(9,*) write(9,*) 'lattice dim = ',l write(9,*) 'temperature = ',ttc write(9,*) 'number of samples = ',mcsmax write(9,*) 'iskip = ',iskip write(9,*) 'seed for generator = ',iseed write(9,*) c c-----set the other parameters c arg = 1.0 + sqrt( 2.0 ) jkt = 0.5 * alog( arg ) / ttc pb = 1.0 - exp( - 2.0 * jkt ) write(9,*) 'pb is ',pb lsq = l * l iseed = 4 * iseed + 1 rcum = 0.0 sus = 0.0 c maxcla = 0 suscl = 0.0 prob = 0.0 c mbond = 0.0 sbond = 0.0 c stackm = 5000 c c-----set initial configuration ( random ) c do 10 j=1,l do 10 i=1,l r = sunif(iseed) - 0.5 spin(i,j) = sign(1.,r) 10 continue c c-----set boundray bonds c do 20 j=1,l lbond(j,1) = 0 ubond(1,j) = 0 20 continue c c *************************************************** c * * c * Monte Carlo Part * c * * c *************************************************** c do 1000 mcs=1,mcsmax c c generate the bonds c cb = 0.0 do 300 j=2,l do 300 i=1,l rb = pb - sunif(iseed) ib = sign(1.,rb) ib = ( ib + 1 ) / 2 ileft = (spin(i,j-1) + spin(i,j)) / 2 ileft = abs(ileft) ib = ib * ileft cb = cb + ib lbond(i,j) = ib 300 continue do 310 j=1,l do 310 i=2,l rb = pb - sunif(iseed) ib = sign(1.,rb) ib = ( ib + 1 ) / 2 iup = (spin(i-1,j) + spin(i,j)) / 2 iup = abs(iup) ib = ib * iup cb = cb + ib ubond(i,j) = ib 310 continue mbond = mbond + cb sbond = sbond + cb * cb c c assign new random orientation c to the clusters c icl = 0 do 400 i=1,l do 400 j=1,l c ileft = lbond(i,j) iup = ubond(i,j) if ((ileft .eq. 1).and.(iup .eq. 0)) spin(i,j) = spin(i,j-1) if ((iup .eq. 1).and.(ileft.eq. 0)) spin(i,j) = spin(i-1,j) if ((ileft .eq. 1).and.(iup .eq. 1)) then c search for the up label isrch = spin(i-1,j) 401 inew = stack(isrch) if ( inew .gt. 0 ) then isrch = inew goto 401 else iu = isrch end if c search for the left label isrch = spin(i,j-1) 402 inew = stack(isrch) if ( inew .gt. 0 ) then isrch = inew goto 402 else il = isrch end if mini = min0(iu,il) maxi = max0(iu,il) if ( mini .ne. maxi ) then stack(maxi) = mini spin(i-1,j) = iu spin(i,j-1) = il end if spin(i,j) = mini end if if ((ileft.eq.0).and.(iup.eq.0)) then icl = icl + 1 if ( icl .gt. stackm ) stop 99 r = sunif(iseed) - 0.5 inew = sign(1.,r) stack(icl) = -2 + inew spin(i,j) = icl end if c 400 continue c rm = 0.0 do 410 i=1,l do 410 j=1,l isrch = spin(i,j) 420 inew = stack(isrch) if ( inew .gt. 0 ) then isrch = inew goto 420 else im = inew + 2 spin(i,j) = im rm = rm + im nsd(isrch) = nsd(isrch) + 1 end if 410 continue c itest = (mcs/iskip) * iskip if ( itest .eq. mcs ) then c compute the magnetic properties icount = icount + 1 rm = rm / lsq sus = sus + rm * rm rcum = rcum + abs(rm) c write(*,*) 'mc step ',mcs,' m = ',rm c c compute the cluster properties c maxcl = 0 do 500 i=1,icl ns = nsd(i) nsd(i) = 0 if ( ns .gt. 0 ) then nsc(ns) = nsc(ns) + 1 nsa(ns) = nsa(ns) + 1 if ( ns .gt. maxcl ) maxcl = ns end if 500 continue prob = prob + maxcl if ( maxcl .gt. maxcla ) maxcla = maxcl c cm = 0.0 do 510 i=1,maxcl-1 ns = nsc(i) nsa1(i) = nsa1(i) + ns nsc(i) = 0 if ( ns .ne. 0 ) then c write(9,*) i,ns end if cm = cm + i * i * ns 510 continue cm = cm / lsq suscl = suscl + cm c c write(9,*) 'mean cluster size (sus.) = ',cm c end if c 1000 continue c c *************************************************** c * * c * A n a l y s i s * c * * c *************************************************** c sus = sus - rcum * rcum /icount sus = sus / icount write(9,*) write(9,*) 'results: program 2dnew2' write(9,*) '-----------------------' write(9,*) write(9,*) 'sample size = ',icount write(9,*) 'average magnetization = ',rcum/icount write(9,*) 'susceptibility = ',sus * lsq c prob = (prob / icount)/ lsq write(9,*) 'percolation probability = ',prob c suscl = suscl / icount write(9,*) 'average sus from cluster = ',suscl c mbond = (mbond/lsq) / icount sbond = (sbond/lsq) / icount write(9,*) 'mean number of bonds = ',mbond write(9,*) 'fluctuation of bonds = ',sbond c write(9,*) 'average cluster size distribution' write(9,*) cm1 = 0.0 cm2 = 0.0 do 1510 i=1,maxcla-1 ns = nsa(i) if ( ns .ne. 0 ) then write(9,*) i,ns cm2 = cm2 + i * i * ns end if 1510 continue cm2 = (cm2 / lsq) / icount write(9,*) 'mean cluster size (sus.) = ',cm2 do 1520 i=1,maxcla-1 ns = nsa1(i) if ( ns .ne. 0 ) then write(9,*) i,ns cm2 = cm2 + i * i * ns end if 1520 continue cm2 = (cm2 / lsq) / icount write(9,*) 'mean cluster size (sus.) = ',cm2 c stop end c mz_soft : random number generators c C**********************************************************************CSUN 600 C CSUN 610 C CSUN 620 C 0 , 1 - U N I F O R M DISTRIBUTION CSUN 630 C CSUN 640 C CSUN 650 C**********************************************************************CSUN 660 C**********************************************************************CSUN 670 C CSUN 680 C FOR DETAILS SEE: CSUN 690 C CSUN 700 C AHRENS, J.H., DIETER, U. AND GRUBE, A. CSUN 710 C PSEUDO-RANDOM NUMBERS: A NEW PROPOSAL CSUN 720 C FOR THE CHOICE OF MULTIPLICATORS CSUN 730 C COMPUTING, 6 (1970), 121 - 138 CSUN 740 C CSUN 750 C**********************************************************************CSUN 760 C SUN 770 REAL FUNCTION SUNIF(IR) SUN 780 DOUBLE PRECISION R,FACTOR,TWO28 SUN 790 C SUN 800 C FACTOR - INTEGER OF THE FORM 8*K+5 AS CLOSE AS POSSIBLE SUN 810 C TO 2**26 * (SQRT(5)-1)/2 (GOLDEN SECTION) SUN 820 C TWO28 = 2**28 (I.E. 28 SIGNIFICANT BITS FOR DEVIATES) SUN 830 C SUN 840 DATA FACTOR /41475557.0D0/, TWO28 /268435456.0D0/ SUN 850 C SUN 860 C RETURNS SAMPLE U FROM THE 0,1 -UNIFORM DISTRIBUTION SUN 870 C BY A MULTIPLICATIVE CONGRUENTIAL GENERATOR OF THE FORM SUN 880 C R := R * FACTOR (MOD 1) . SUN 890 C IN THE FIRST CALL R IS INITIALIZED TO SUN 900 C R := IR / 2**28 , SUN 910 C WHERE IR MUST BE OF THE FORM IR = 4*K+1. SUN 920 C THEN R ASSUMES ALL VALUES 0 < (4*K+1)/2**28 < 1 DURING SUN 930 C A FULL PERIOD 2**26 OF SUNIF. SUN 940 C THE PARAMETER IR IS USED ONLY IN THE FIRST CALL FOR SUN 950 C INITIALIZATION OF SUNIF. THEREAFTER (WHEN NEGATIVE) SUN 960 C IR BECOMES A DUMMY VARIABLE. SUN 970 C SUN 980 IF (IR .GE. 0) GO TO 1 SUN 990 C SUN 1000 C STANDARD CASE: SAMPLING SUN 1010 C SUN 1020 R=DMOD(R*FACTOR,1.0D0) SUN 1030 SUNIF=SNGL(R) SUN 1040 RETURN SUN 1050 C SUN 1060 C FIRST CALL: INITIALIZATION SUN 1070 C SUN 1080 1 R=DBLE(FLOAT(IR))/TWO28 SUN 1090 R=DMOD(R*FACTOR,1.0D0) SUN 1100 SUNIF=SNGL(R) SUN 1110 IR=-1 SUN 1120 RETURN SUN 1130 END SUN 1140 C SEX 10 C**********************************************************************CSEX 20