c---------------------------------------------------------------------- ccc ISING2D c---------------------------------------------------------------------- ccc This program is a simulation of the 2-D ising model using ccc the Metropolis Monte Carlo sampling technique. Energy and ccc temperatures are measured in units of J, the interaction ccc between nearest neighbors. Periodic boundary conditions ccc are used throughout. The code is written in an interactive ccc mode. The results are output to the screen and also written ccc to a log file which the user names. ccc The external magnetic field is set to zero. ccc You may want to experiment with other updating techniques, ccc eg. 2-spin flips or cluster moves [cf. Swendsen, R. H., ccc Comp. Phys. Comm. 65, 281 (1991) for review]. ccc You may wish to modify the code to simplify the operation, ccc improve the performace, and help ensure you get good results. I ccc would appreciate it if you would provide me with a copy ccc of your source if you make any significant enhancements. ccc To compile this code on a Sun workstation, do the following: ccc In the directory in which the source code resides, type: ccc f77 -fast ising2d.f ccc This will create an executable called a.out. If you would ccc prefer an executable with another name, such as mc.exe, type ccc f77 -fast -o mc.exe ising2d.f ccc Note that -fast turns optimization on, and makes the code run a ccc bit faster. The adventurous student may wish to fool around ccc with compilation options (man f77) to improve performance. ccc The code has not been tested on other architectures (i.e. PC, ccc SGI, etc. Parts may not work, especially the wtime subroutine. ccc Appropriate modifications will be necessary. ccc You can run the program interactively by typing a.out and ccc answering the questions as you're prompted. Runs can take ccc anywhere from a few minutes to significantly longer, depending ccc on the parameters you give. You may want to automate things ccc and run "batch" jobs by editing the source code. Feel free ccc to do this. ccc This program is for instructional use in ChEg 698 D ccc "Molecular Modeling and Theory" at the University of Notre Dame. ccc This code may be copied but please give reference to it's ccc sources. Large parts of the code come from Prof. David ccc Chandler's group at Berkeley, courtesy of Kevin Yeung and the ccc rest was added by Ed Maginn. I think it works OK, but please ccc report any bugs you find to me at ed@nd.edu. c---------------------------------------------------------------------- program ising2d ccc note to my students and anyone else who does programming: ccc implicit declarations are a bad thing - do *not* use them in ccc research codes! Mea culpa. implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter(maxl=40) common/lattice/latt(maxl,maxl),nsize common/stat/embe(4),etot,emean,efluc,rmmean,rmfluc,rj,rj2,mtot common/misc/iseed character*30 logfile character*70 message c---------------------------------------------------------------------- ccc initialize random number generator. ccc note: function uni changes iseed every time it is called. ccc if you'd like to run a simulation from a different starting ccc point but under the same conditions, read iseed in as a variable ccc and change it from run to run. This will cause a different ccc sequence of random numbers to be generated. You may also wish ccc to use your own favorite random number generator... iseed=305 useed=ustart(iseed) ccc read in user inputs write(*,*) 'input name of log file results written to' read(*,'(a30)') logfile open(1,file=logfile,status='new') message = 'beginning Ising simulation' CALL wtime(message) write(*,*) 'input lattice edge size ( <= 40 ):' read(*,*) nsize if ((nsize.le.1).or.(nsize.gt.maxl)) go to 2000 if (mod(nsize,2).ne.0) nsize=nsize-1 write(1,998) 'size:',nsize ccc How do you want to start out write(*,*) 'initialization of spins:' write(*,*) '1- random' write(*,*) '2- interface' write(*,*) '3- checker board' write(*,*) '4- read from config.dat' write(*,*) 'option:' read(*,*) inin write(1,998) 'initialization:',inin ccc set the initialization if (inin.eq.1) then CALL inran else if (inin.eq.2) then CALL inter else if (inin.eq.3) then CALL inckr else CALL inext end if 1000 write(*,*) write(1,*) write(*,*) 'enter number of MC steps, 0=stop simulation:' read(*,*) nstep if (nstep.le.0) go to 2000 write(*,*) 'Total number of MC moves: ' write(*,8) real(nsize)**2 * real(nstep) 8 format(E15.7) write(1,998) 'number of steps:',nstep write(*,*) 'enter reduced temperature, kT/J:' read(*,*) t write(1,999) 'kT/J:',t if (t.le.1.d-5) t=1d-5 rj=1.d0/t rj2=2.d0*rj ccc get initial energy of lattice CALL echeck write(*,*)'initial energy of lattice ', etot write(*,*)'initial magnetization of lattice ', mtot cccc sets up 'look-up' table for energies embe(2)=dexp(-2.d0*rj2) embe(4)=dexp(-4.d0*rj2) ccc zero the counters emean=0.d0 efluc=0.d0 rmmean=0.d0 rmfluc=0.d0 ccc run 'nstep' number of 'cycles,' each consisting of ccc flipping nsize*nsize spins. the spins are chosen ccc at random. do istep=1,nstep CALL update ccc return with energy and magnetization ccc report values periodically: comment out if you don't like ccc this... if (nstep .gt. 10) then if (mod(istep,int(nstep/10)) .eq. 0) then write(*,*) 'At step: ',istep write(*,*) 'energy, magnetization ' write(*,*) etot,mtot endif endif ccc sum energies, magnetization, and get fluctuations emean=emean+etot efluc=efluc+etot*etot rmmean=rmmean+dfloat(mtot) rmfluc=rmfluc+dfloat(mtot)*dfloat(mtot) enddo ! istep = 1,nstep message = 'Finished with Ising simulation' CALL wtime(message) ccc get means and fluctuations on a per-spin basis emean=emean/dfloat(nstep) efluc=(efluc/dfloat(nstep)-emean*emean)**0.5d0 emean=emean/dfloat(nsize*nsize) efluc=efluc/dfloat(nsize*nsize) rmmean=rmmean/dfloat(nstep) rmfluc=(rmfluc/dfloat(nstep)-rmmean*rmmean)**0.5d0 rmmean=rmmean/dfloat(nsize*nsize) rmfluc=rmfluc/dfloat(nsize*nsize) write(*,999) 'mean energy per spin:',emean write(*,999) 'rms energy fluctuation per spin:',efluc write(*,999) 'mean magnetization per spin:',rmmean write(*,999) 'rms magnetization fluctuation per spin:',rmfluc write(1,999) 'mean energy per spin:',emean write(1,999) 'rms energy fluctuation per spin:',efluc write(1,999) 'mean magnetization per spin:',rmmean write(1,999) 'rms magnetization fluctuation per spin:',rmfluc write(*,*) write(*,*) 'dump spin configuration to file:' write(*,*) 'enter 0 to skip dump' write(*,*) 'enter # to dump in fort.# (11<#<20)' write(*,*) 'option:' read(*,*) idump if ((idump.ge.11).and.(idump.le.20)) then do 50 iy=1,nsize write(idump,997) ((latt(ix,iy)+1)/2,ix=1,nsize) 50 continue close(idump,status='keep') end if go to 1000 2000 write(*,*) 'Finished with ising2d' ccc see if you want to save the last configuration' write(*,*) 'Save last config in config.dat?' write(*,*) '1 = yes' read(*,*) isave if (isave .eq. 1) then open(unit=20,file='config.dat') do iy = 1,nsize write(20,997) ((latt(ix,iy)+1)/2,ix=1,nsize) enddo close(unit=20) endif 997 format(40i2) 998 format(A40,i10) 999 format(A40,f15.8) message = 'Leaving Ising2d' CALL wtime(message) close(1,status='keep') stop 'bye now' end c---------------------------------------------------------------------- ccc UPDATE c---------------------------------------------------------------------- ccc update nsize*nsize spins chosen at random. employ the ccc Metropolis scheme. c---------------------------------------------------------------------- subroutine update implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter(maxl=40) common/lattice/latt(maxl,maxl),nsize common/stat/embe(4),etot,emean,efluc,rmmean,rmfluc,rj,rj2,mtot common/misc/iseed dimension nhxl(4),nhyl(4) data nhxl/-1,1,0,0/ data nhyl/0,0,-1,1/ c---------------------------------------------------------------------- size=dfloat(nsize) nsize1=nsize-1 ccc Pick a spin, any spin... do iflip=1,nsize*nsize ix=int(size*uni())+1 iy=int(size*uni())+1 ccc Do the Metropolis thing... nhsum=0 do nhbr=1,4 nhx=mod(ix+nhxl(nhbr)+nsize1,nsize)+1 nhy=mod(iy+nhyl(nhbr)+nsize1,nsize)+1 nhsum=nhsum+latt(nhx,nhy) enddo itest=nhsum*latt(ix,iy) if (itest.le.0) then ccc This is a move accept etot=etot+dfloat(itest)*rj2 mtot=mtot-2*latt(ix,iy) latt(ix,iy)=-latt(ix,iy) else test=uni() ccc Conditional accept if (test.le.embe(itest)) then etot=etot+dfloat(itest)*rj2 mtot=mtot-2*latt(ix,iy) latt(ix,iy)=-latt(ix,iy) end if end if enddo return end c---------------------------------------------------------------------- subroutine inran c c sets up random initial configurations c appropriate for infinite temperature. c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter(maxl=40) common/lattice/latt(maxl,maxl),nsize common/misc/iseed do 10 iy=1,nsize do 20 ix=1,nsize latt(ix,iy)=int(uni()+0.5d0)*2-1 20 continue 10 continue return end c---------------------------------------------------------------------- subroutine inter c c sets up interfacial initial configuration c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter(maxl=40) common/lattice/latt(maxl,maxl),nsize nsize2=nsize/2 do 10 iy=1,nsize do 20 ix=1,nsize2 latt(ix,iy)=1 latt(ix+nsize2,iy)=-1 20 continue 10 continue return end c---------------------------------------------------------------------- subroutine inckr c c sets up checkered initial configurations c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter(maxl=40) common/lattice/latt(maxl,maxl),nsize do 10 iy=1,nsize do 20 ix=1,nsize latt(ix,iy)=(-1)**(ix+iy) 20 continue 10 continue return end c---------------------------------------------------------------------- subroutine inext c c read in initial configuration c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter(maxl=40) common/lattice/latt(maxl,maxl),nsize open(20,file='config.dat') do 10 iy=1,nsize read(20,*) (latt(ix,iy),ix=1,nsize) do 20 ix=1,nsize latt(ix,iy)=latt(ix,iy)*2-1 20 continue 10 continue close(20,status='keep') return end c---------------------------------------------------------------------- subroutine echeck c c calculates initial energy c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter(maxl=40) common/lattice/latt(maxl,maxl),nsize common/stat/embe(4),etot,emean,efluc,rmmean,rmfluc,rj,rj2,mtot dimension nhxl(4),nhyl(4) data nhxl/-1,1,0,0/ data nhyl/0,0,-1,1/ nsize1=nsize-1 ie=0 im=0 do 10 iy=1,nsize do 20 ix=1,nsize nhsum=0 do 30 nhbr=1,4 nhx=mod(ix+nhxl(nhbr)+nsize1,nsize)+1 nhy=mod(iy+nhyl(nhbr)+nsize1,nsize)+1 nhsum=nhsum+latt(nhx,nhy) 30 continue ie=ie-nhsum*latt(ix,iy) im=im+latt(ix,iy) 20 continue 10 continue etot=dfloat(ie)*rj/2.d0 mtot=im return end c---------------------------------------------------------------------- REAL*8 FUNCTION UNI() c c (Courtesy of Raymond Yee who suggested that I use this c public domain code.) c C***BEGIN PROLOGUE UNI C***DATE WRITTEN 810915 (YYMMDD) C***REVISION DATE 871210 (YYMMDD) C***CATEGORY NO. L6A21 C***KEYWORDS RANDOM NUMBERS, UNIFORM RANDOM NUMBERS C***AUTHOR KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS C MARSAGLIA, GEORGE, SUPERCOMPUTER RES. INST., FLORIDA ST. U. C C***PURPOSE THIS ROUTINE GENERATES REAL (SINGLE PRECISION) UNIFORM C RANDOM NUMBERS ON [0,1) C***DESCRIPTION C Computes real (single precision) uniform numbers on [0,1). C From the book, "Numerical Methods and Software" by C D. Kahaner, C. Moler, S. Nash C Prentice Hall, 1988 C C***REFERENCES MARSAGLIA G., "COMMENTS ON THE PERFECT UNIFORM RANDOM C NUMBER GENERATOR", UNPUBLISHED NOTES, WASH S. U. C***ROUTINES CALLED (NONE) C***END PROLOGUE UNI C PARAMETER( * CSAVE=362436./16777216. , * CD=7654321./16777216., * CM=16777213./16777216. ) C 2**24=16777216 REAL*8 U(17),S,T,USTART,C,UNIB INTEGER*4 I,J,II,JJ,K,KK,I1,J1,K1,L1,M1,ISEED C SAVE U,I,J,K,C C Load data array in case user forgets to initialize. C This array is the result of calling UNI 100000 times C with ISEED=305 and K=64. DATA U/ *0.8668672834288, 0.3697986366357, 0.8008968294805, *0.4173889774680, 0.8254561579836, 0.9640965269077, *0.4508667414265, 0.6451309529668, 0.1645456024730, *0.2787901807898, 0.06761531340295, 0.9663226330820, *0.01963343943798, 0.02947398211399, 0.1636231515294, *0.3976343250467, 0.2631008574685/ DATA I,J,K,C/17,5,24,CSAVE/ C C Basic generator is Fibonacci C UNI = U(I)-U(J) IF(UNI.LT.0.0)UNI = UNI+1.0 U(I) = UNI I = I-1 IF(I.EQ.0)I = 17 J = J-1 IF(J.EQ.0)J = 17 C C Second generator is congruential C C = C-CD IF(C.LT.0.0) C=C+CM C C Combination generator C UNI = UNI-C IF(UNI.LT.0.0)UNI = UNI+1.0 RETURN C ENTRY USTART(ISEED) C C Set up ... C Convert ISEED to four smallish positive integers. C I1 = MOD(ABS(ISEED),177)+1 J1 = MOD(ABS(ISEED),167)+1 K1 = MOD(ABS(ISEED),157)+1 L1 = MOD(ABS(ISEED),147)+1 C C Generate random bit pattern in array based on given seed. C DO 2 II = 1,17 S = 0.0 T = 0.5 C Do for each of the bits of mantissa of word C Loop over K bits, where K is defaulted to 24 but can C be changed by user call to UNIB(K) DO 3 JJ = 1,K M1 = MOD(MOD(I1*J1,179)*K1,179) I1 = J1 J1 = K1 K1 = M1 L1 = MOD(53*L1+1,169) IF(MOD(L1*M1,64).GE.32)S=S+T 3 T = .5*T 2 U(II) = S USTART = FLOAT(ISEED) RETURN C ENTRY UNIB(KK) IF(KK.LE.24)THEN K=24 ELSE K=KK ENDIF UNIB=FLOAT(K) END c---------------------------------------------------------------------- ccc WTIME c---------------------------------------------------------------------- ccc Reports cpui timing to log file c---------------------------------------------------------------------- subroutine wtime(message) character*70 version character*70 message integer istart common /verinfo/version data istart/1/ save istart character fdate*24,name*20 real*4 tarray(2) c---------------------------------------------------------------------- if (istart.eq.1) then istart=0 write(1,100) fdate() 100 format(' *** run_date: ',a24,/) i=hostnm(name) write(1,'(a)') ' *** machine: ',name(1:lentrm(name)),' ***' return endif r=etime(tarray) write(1,120) message(1:lentrm(message)) 120 format(1x,'>>> ',a) write(1,110) r 110 format(' total elapsed time =',f12.1,' seconds',/) return end FUNCTION LENTRM (STRING) CHARACTER *(*) STRING DO 10 LENTRM = LEN (STRING), 1, -1 IF (STRING (LENTRM:LENTRM) .NE. ' ') RETURN 10 CONTINUE LENTRM = 0 END