PROGRAM Triangle_Percolation REAL P INTEGER SITE(150,150),CLUSTER(500),M,L1 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) CALL FIX(L1,SITE,CLUSTER) OPEN(UNIT=1,FILE='percol.xyz',IOSTAT=KODE) WRITE(1,FMT='(A)') ' ' CLOSE(1) CALL BOOTS(L1,M,SITE) 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 triangluar 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<50) ?' read(*,*) L IF(L.GT.50 ) THEN L=10 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 END IF WRITE (*,*)'The m parameter(number of nearest neighbours is ?' READ(*,*) M IF(M.LT.1 .OR. M.GT.4)THEN M=2 END IF RETURN 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,CLUSTER(10*L1) L=L1-1 MAX=10*L DO 9,I=1,L1 DO 9, J=1,L1 MS=SITE(I,J) IF (MS.NE.MAX)THEN 8 MS=CLUSTER(MS) IF (MS.NE.CLUSTER(MS))THEN GOTO 8 END IF SITE(I,J)=MS ELSE SITE(I,J)=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),INDEX,CLUSTER(10*L1),MAX,TOPR,LEFT,TOPL,NEW REAL SEED c READ(*,*)SEED SEED=487922791 INDEX=0 ! label index L=L1-1 MAX=L*10 DO 10,I=1,L SITE(1,I)=MAX !boundaries are set to max SITE(I,1)=MAX SITE(I,L1)=MAX SITE(L1,I)=MAX 10 CONTINUE DO 11,I=1,MAX CLUSTER(I)=MAX 11 CONTINUE C INITIAL CONDITIONS ARE SET DO 20, I=2,L DO 20,J=2,L IF (RANNOS(SEED) .GT. P) THEN SITE(I,J)=MAX GOTO 20 !NEW SITE IS EMPTY cccccccccccccccc Hoshen-Kopelmen algorithm starts here ELSE I1=MOD(I,2) ! determines whether the shift is left or right TOPL=SITE(I-1,J+I1) LEFT=SITE(I,J-1) TOPR=SITE(I-1,J+(I1-1)) IF ((TOPL+LEFT+TOPR).EQ. 3*MAX) THEN INDEX=INDEX+1 !NEW CLUSTER CLUSTER(INDEX)=INDEX SITE(I,J)=INDEX GOTO 20 END IF IF(TOPL.EQ.MAX.OR.TOPL.EQ.CLUSTER(TOPL))THEN GOTO 14 END IF MS=TOPL 13 MS=CLUSTER(MS) IF (MS.NE.CLUSTER(MS))THEN GOTO 13 END IF CLUSTER(TOPL)=MS 14 IF(TOPR.EQ.MAX.OR.TOPR.EQ.CLUSTER(TOPR)) THEN GOTO 15 END IF MS=TOPR 16 MS=CLUSTER(MS) IF (MS.NE.CLUSTER(MS)) THEN GOTO 16 END IF CLUSTER(TOPR)=MS 15 NEW=MIN(CLUSTER(TOPR),CLUSTER(LEFT),CLUSTER(TOPL)) SITE(I,J)=NEW IF(TOPR.NE.MAX.AND.CLUSTER(TOPR).NE.NEW)THEN CLUSTER(TOPR)=NEW END IF IF(LEFT.NE.MAX.AND.CLUSTER(LEFT).NE.NEW)THEN CLUSTER(LEFT)=NEW END IF IF(TOPL.NE.MAX.AND.CLUSTER(TOPL).NE.NEW)THEN CLUSTER(TOPL)=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),B LOGICAL FLAG L=L1-1 35 FLAG=.FALSE. !RELATIVE POSITION IN THE TRIANGLE LATTICE CALL FILESAV(L1,SITE) DO 55,I=2,L DO 55,J=2,L J1=MOD(I,2) IF(SITE(I,J).NE.0)THEN K1=SITE(I-1,J+J1) K2=SITE(I-1,J+J1-1) K3=SITE(I,J-1) K4=SITE(I,J+1) K5=SITE(I+1,J+J1) K6=SITE(I+1,J+J1-1) KS=K1+K2+K3+K4+K5+K6 IF(KS.LT.(B*SITE(I,J))) THEN SITE(I,J)=0 FLAG=.TRUE. 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 OF THE XYZ FILES CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL SIT(2500,4) INTEGER SITE(L1,L1),N CHARACTER FILENAM*11,GENFILE*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 P2=SQRT(3.)/2. DO 100,I=2,L DO 100,J=2,L c IF (SITE(I,J).NE.0.OR.SITE(I,J).NE.10*L)THEN M=M+1 SIT(M,1)=SITE(I,J) SIT(M,2)=(I-1)*P2 SIT(M,3)=(J+1+MOD(I,2)*.5) 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(INT(SIT(I,1))) ! GIVES A TWO LETTER NAME TO THE CLUSTER WRITE(1,102) CLUSS,SIT(I,2),SIT(I,3),SIT(I,4) 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 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