PROGRAM Square_Percolation REAL P INTEGER SITE(150,150,150),CLUSTER(800),M,L cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c here is a discription of some of the major variables : c p - probability of occupation c site - the lattice array c cluster - an array of proper labels,each cell is the proper label for c the cluster with the same index as the cell c m - m parameter of BP c L - the size of the lattice wanted + 2 c L1 =L+1- the array size transfered between subroutines c max - an upper limit for the number of clusters cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc CALL USER(L1,P,M) CALL INITIAL(P,L1,SITE,CLUSTER) C CALL PLOTS(SITE,L1,CLUSTER) CALL FIX(L1,SITE,CLUSTER) c CALL PLOTS(SITE,L1,CLUSTER) OPEN(UNIT=1,FILE='percol.xyz',IOSTAT=KODE) WRITE(1,FMT='(A)') ' ' CLOSE(1) CALL BOOTS(L1,M,SITE) C CALL PLOTS(SITE,L1,CLUSTER) STOP END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE USER(L1,P,M) C this subroutine gets from the user the lattice size, probability and m c parameter for the calculation c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc INTEGER M write(*,*) 'This Program simulates bootstrap percolation in ' write(*,*)' a 2D square lattice.' write(*,*)'You need to state the lattice size,probability and' write(*,*)'m (minimum number of nearest neighbours)' write(*,*) write(*,*)'The lattice size is (L<145) ?' read(*,*) L IF(L.GT.145 ) THEN L=10 !default value END IF L1=L+2 WRITE(*,*)' The probablity of site occupation ?' READ(*,*) P IF(P.LE.0 .OR.P .GE. 1) THEN P=0.5 !default value END IF WRITE (*,*)'The m parameter(number of nearest neighbours is ?' READ(*,*) M IF(M.LT.1 .OR. M.GT.4)THEN M=2 !default value END IF RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE PLOTS(SITE,L1,CLUSTER) c this subroutine plot the array on screen , and is useful to check the results c when dealing with a small array.It is NOT used otherwise CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER SIT(L1-2),SITE(L1,L1,L1),L1,CLUSTER(2*L1) L=L1-1 MAX=5*L K=2 DO 5,I=2,L DO 1, J=2,L IF (SITE(I,J,K).EQ.MAX) THEN SIT(J-1)=0 ELSE SIT(J-1)=SITE(I,J,K) END IF 1 CONTINUE WRITE(*,*) SIT 5 CONTINUE WRITE(*,*) 7 CONTINUE WRITE(*,*) CLUSTER END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE FIX(L1,SITE,CLUSTER) C CHANGES LABELS IN SITE TO THE PROPER ONES.once this subroutine is used the c array contains zeros for empty cells,and the cluster number at each c occupied site. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER SITE(L1,L1,L1),L1,CLUSTER(5*L1) L=L1-1 MAX=5*L K=2 DO 9,I=1,L1 DO 9, J=1,L1 MS=SITE(I,J,K) IF (MS.NE.MAX)THEN 8 MS=CLUSTER(MS) IF (MS.NE.CLUSTER(MS))THEN GOTO 8 END IF SITE(I,J,K)=MS ELSE SITE(I,J,K)=0 END IF 9 CONTINUE END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE INITIAL(P,L1,SITE,CLUSTER) C This subroutine initiates the lattice, as far as occupation of the sites c and cluster labeling ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc REAL P INTEGER SITE(L1,L1,L1),INDEX,CLUSTER(5*L1),MAX,TOP,LEFT,BACK,NEW REAL SEED c READ(*,*)SEED SEED=487922791 INDEX=0 ! label index L=L1-1 MAX=L*5 DO 10,I=1,L DO 10,J=1,L SITE(I,J,1)=MAX !boundaries are set to max SITE(1,I,J)=MAX SITE(I,1,J)=MAX SITE(I,L1,J)=MAX SITE(I,J,L1)=MAX SITE(L1,I,J)=MAX 10 CONTINUE DO 11,I=1,MAX CLUSTER(I)=MAX 11 CONTINUE C INITIAL CONDITIONS ARE SET K=2 DO 20, I=2,L DO 20,J=2,L IF (RANNOS(SEED) .GT. P) THEN SITE(I,J,K)=MAX GOTO 20 !NEW SITE IS EMPTY cccccccccccccccc Hoshen-Kopelmen algorithm starts here ELSE TOP=SITE(I-1,J,K) LEFT=SITE(I,J-1,K) IF ((TOP+LEFT).EQ. 2*MAX) THEN INDEX=INDEX+1 !NEW CLUSTER CLUSTER(INDEX)=INDEX SITE(I,J,K)=INDEX GOTO 20 END IF IF(TOP.EQ.MAX.OR.TOP.EQ.CLUSTER(TOP))THEN GOTO 15 END IF MS=TOP 13 MS=CLUSTER(MS) ! finding the proper label of the site above IF (MS.NE.CLUSTER(MS))THEN GOTO 13 END IF CLUSTER(TOP)=MS 15 NEW=MIN(CLUSTER(TOP),CLUSTER(LEFT)) ! finding the proper label for SITE(I,J,K)=NEW ! the new site IF(TOP.NE.MAX.AND.CLUSTER(TOP).NE.NEW)THEN CLUSTER(TOP)=NEW END IF IF(LEFT.NE.MAX.AND.CLUSTER(LEFT).NE.NEW)THEN CLUSTER(LEFT)=NEW END IF END IF 20 CONTINUE END cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE BOOTS(L1,B,SITE) C This subroutine performs the culling process associated with BP CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER SITE(L1,L1,L1),B LOGICAL FLAG ! a flag to show whether any change was done L=L1-1 35 FLAG=.FALSE. CALL FILESAV(L1,SITE) !SAVES THE INITIAL ARRAY DO 55,I=2,L DO 55,J=2,L K=2 IF(SITE(I,J,K).NE.0)THEN K1=SITE(I-1,J,K) K2=SITE(I+1,J,K) K3=SITE(I,J-1,K) K4=SITE(I,J+1,K) KS=K1+K2+K3+K4 IF(KS.LT.(B*SITE(I,J,K))) THEN ! are there enough nearest neighbours ? SITE(I,J,K)=0 FLAG=.TRUE. ! a change was made CALL FILESAV(L1,SITE) END IF END IF 55 CONTINUE IF (FLAG)THEN GOTO 35 END IF RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE FILESAV(L1,SITE) C THIS SUBROUTINE SAVES THE ARRAY AS AN XYZ FILE , AND CREATES A TEXT C file OF THE XYZ FILES CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER SITE(L1,L1,L1),N,NUM1,NUM2,SIT(8000,4) CHARACTER FILENAM*11,GENFILE*11,JUNK*11,CLUSS*2,CLUS*2 SAVE N !xyz file index DATA N /1/ GENFILE='percol.xyz'! the name of the file containing the file list FILENAM='file .xyz'! the name of the xyz files FILENAM(5:5)=CHAR(N/100+48) FILENAM(6:6)=CHAR(MOD(N,100)/10+48) FILENAM(7:7)=CHAR(MOD(N,10)+48) !adding the index to the file name L=L1-1 M=0 DO 100,I=2,L DO 100,J=2,L K=2 c IF(SITE(I,J,K).NE.0) THEN M=M+1 !the number of sites in the xyz file SIT(M,1)=SITE(I,J,K) SIT(M,2)=I-1 SIT(M,3)=J-1 SIT(M,4)=1 C END IF 100 CONTINUE c at the moment the program saves the empty sites as well. this can be avoided by removing the c comment marks of the IF expression. AVIZ changes the dimension of the picture according to the c number of sites,so in order to maintain constant size all the sites are being saved. OPEN(UNIT=1,FILE=FILENAM,IOSTAT=KODE) IF (KODE.NE.0)THEN WRITE(*,*)'FILE ERROR1' STOP END IF WRITE(1,FMT='(I5)')M WRITE(1,FMT='(A10)')'##CLUSTER' DO 101,I=1,M CLUSS=CLUS(SIT(I,1)) WRITE(1,102) CLUSS,SIT(I,2)*1.,SIT(I,3)*1.,SIT(I,4)*1. 101 CONTINUE 102 FORMAT (A6,' ',F6.1,' ',F6.1,' ',F6.1) CLOSE(UNIT=1) OPEN(UNIT=1,FILE=GENFILE,IOSTAT=KODE) IF (KODE.NE.0)THEN WRITE(*,*)'FILE ERROR2' STOP END IF 110 READ(UNIT=1,FMT=130,END=120) GOTO 110 120 WRITE(UNIT=1,FMT=130)FILENAM CLOSE(UNIT=1) 130 FORMAT(A) N=N+1 END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CHARACTER*(*) FUNCTION CLUS(N1) C THE FUNCTION RETURNS A LABEL FOR THE CLUSTER - changes the number into c a two letter combination CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER N1,N2,N3 N2=N1/26 N3=MOD(N1,26) CLUS(1:1)=CHAR(65+N2) CLUS(2:2)=CHAR(65+N3) RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL FUNCTION RANNOS(DSEED) C returns a uniformly distributed random number between 0 and 1 c from Koonin & Meredith - Computational physics CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL DSEED REAL D2P31M,D2P31 DATA D2P31M/2147483647.D0/ DATA D2P31 /2147483711.D0/ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DSEED = MOD(16807.D0*DSEED,D2P31M) RANNOS = DSEED / D2P31 RETURN END