C ******************************************************************* C ** FICHE F.36. MONTE CARLO SIMULATION OF HARD LINES IN 2D ** C ******************************************************************* C ******************************************************************* C ** TWO SEPARATE PARTS: FORTRAN AND BASIC VERSIONS. ** C ******************************************************************* C ******************************************************************* C ** FICHE F.36 - PART A ** C ** FORTRAN PROGRAM FOR MONTE CARLO SIMULATION OF HARD LINES ** C ******************************************************************* PROGRAM HLINES COMMON / BLOCK1 / RX, RY, EX, EY C ******************************************************************* C ** MONTE CARLO SIMULATION OF HARD LINES IN 2D. ** C ** ** C ** THIS PROGRAM SETS UP AN INITIAL CONFIGURATION OF HARD LINES ** C ** AND CONDUCTS AN ENTIRE SWEEP IN DENSITY FROM LOW TO HIGH ** C ** OR HIGH TO LOW, PRINTING OUT THE ORDER PARAMETER REGULARLY. ** C ** THE LINE LENGTH IS TAKEN TO BE UNITY THROUGHOUT. ** C ** ** C ** PRINCIPAL VARIABLES: ** C ** ** C ** INTEGER N NUMBER OF LINES ** C ** REAL RX(N),RY(N) C-O-M POSITIONS OF LINES ** C ** REAL EX(N),EY(N) UNIT VECTOR ALONG LINE ** C ******************************************************************* INTEGER N PARAMETER ( N = 25 ) REAL RX(N), RY(N), EX(N), EY(N) INTEGER MAXST PARAMETER ( MAXST = 15 ) INTEGER NEQUIL(MAXST), NPROD(MAXST) REAL DENS(MAXST) REAL MAXDEN PARAMETER ( MAXDEN = N / 4.0 ) LOGICAL LOSTRT INTEGER STATE, NSTATE REAL BOX, NEWBOX, MAXDIS, MAXROT, AVORD, RATDIS, RATROT C ******************************************************************* C ** ENTER PARAMETERS ** WRITE(*,'(1H1,'' **** PROGRAM HLINES **** '')') WRITE(*,'(//'' MONTE CARLO SIMULATION OF HARD LINES IN 2D '')') WRITE(*,'('' HOW MANY STATE POINTS WOULD YOU LIKE ? '')') READ (*,*) NSTATE WRITE(*,'('' LOW-DENSITY START CONFIGURATION (.T./.F.) ? '')') READ (*,*) LOSTRT WRITE(*,'('' MAXIMUM ALLOWED DENSITY IS '',F10.6)') MAXDEN WRITE(*,'('' DENS = DENSITY '')') WRITE(*,'('' NEQUIL = NUMBER OF SWEEPS FOR EQUILIBRATION '')') WRITE(*,'('' NPROD = NUMBER OF SWEEPS FOR PRODUCTION '')') DO 10 STATE = 1, NSTATE WRITE(*,'('' DENS, NEQUIL, NPROD FOR STATE '',I5)') STATE READ (*,*) DENS(STATE), NEQUIL(STATE), NPROD(STATE) IF ( DENS(STATE) .GT. MAXDEN ) THEN WRITE(*,'('' SORRY - THAT ONE IS TOO BIG! '')') STOP ENDIF 10 CONTINUE WRITE(*,'('' NUMBER OF STATE POINTS = '',I6)') NSTATE WRITE(*,'('' DESIRED PARAMETERS: '')') WRITE(*,'('' STATE DENSITY EQUILIBR PRODUCT '')') DO 20 STATE = 1, NSTATE WRITE(*,'(1X,I5,F10.6,2I10)') + STATE, DENS(STATE), NEQUIL(STATE), NPROD(STATE) 20 CONTINUE C ** SET UP INITIAL CONFIGURATION ** BOX = SQRT ( REAL ( N ) / DENS(1) ) IF ( LOSTRT ) THEN WRITE(*,'('' SETTING UP LOW-DENSITY START '')') WRITE(*,'('' BOX LENGTH = '',F15.5)') BOX CALL LODEN ( BOX ) ELSE WRITE(*,'('' SETTING UP HIGH-DENSITY START '')') WRITE(*,'('' BOX LENGTH = '',F15.5)') BOX CALL HIDEN ( BOX ) ENDIF MAXDIS = 0.1 MAXROT = 0.1 WRITE(*,10001) C ******************************************************************* C ** MAIN LOOP STARTS ** C ******************************************************************* DO 1000 STATE = 1, NSTATE NEWBOX = SQRT( REAL(N) / DENS(STATE) ) CALL EQUIL ( NEQUIL(STATE), BOX, NEWBOX, : MAXDIS, MAXROT, RATDIS, RATROT ) DENS(STATE) = REAL(N) / BOX ** 2 WRITE(*,10002) : STATE, DENS(STATE), MAXDIS, RATDIS, MAXROT, RATROT CALL PRODUC ( NPROD(STATE), BOX, AVORD, : MAXDIS, MAXROT, RATDIS, RATROT ) WRITE(*,10003) : STATE, DENS(STATE), MAXDIS, RATDIS, MAXROT, RATROT, AVORD 1000 CONTINUE C ******************************************************************* C ** MAIN LOOP ENDS ** C ******************************************************************* STOP 10001 FORMAT (1X,'......STATE ...DENSITY ', : '....MAXDIS RATDIS.... ....MAXROT RATROT.... ', : 'ORDER.....') 10002 FORMAT (1X,'EQUIL ',I5,1X,F10.5,1X, : 2(F10.5,1X,G10.3,1X)) 10003 FORMAT (1X,'PRODU ',I5,1X,F10.5,1X, : 2(F10.5,1X,G10.3,1X),F10.5) END SUBROUTINE HIDEN ( BOX ) COMMON / BLOCK1 / RX, RY, EX, EY C ******************************************************************* C ** ROUTINE TO SET UP A HIGH DENSITY START CONFIGURATION ** C ** ** C ** LINE CENTRES ARE TRANSLATIONALLY DISORDERED, BUT ALL THE ** C ** LINES POINT IN THE SAME DIRECTION. THIS DIRECTION IS CHOSEN ** C ** RANDOMLY. THE CONFIGURATION CANNOT CONTAIN OVERLAPS, BUT WE ** C ** CHECK JUST IN CASE. ** C ******************************************************************* INTEGER N PARAMETER ( N = 25 ) REAL RX(N), RY(N), EX(N), EY(N) REAL BOX REAL TWOPI PARAMETER ( TWOPI = 2.0 * 3.1415927 ) INTEGER I LOGICAL OVRLAP REAL TH, COSTH, SINTH, RANF, DUMMY, COSVAL, SINVAL C ******************************************************************* TH = ( RANF ( DUMMY ) - 0.5 ) * TWOPI COSTH = COS ( TH ) SINTH = SIN ( TH ) DO 100 I = 1, N OVRLAP = .TRUE. 50 IF ( OVRLAP ) THEN RX(I) = ( RANF ( DUMMY ) - 0.5 ) * BOX RY(I) = ( RANF ( DUMMY ) - 0.5 ) * BOX EX(I) = COSTH EY(I) = SINTH OVRLAP = .FALSE. CALL DNCHEK ( BOX, I, OVRLAP ) GOTO 50 ENDIF 100 CONTINUE RETURN END SUBROUTINE LODEN ( BOX ) COMMON / BLOCK1 / RX, RY, EX, EY C ******************************************************************* C ** ROUTINE TO SET UP A LOW DENSITY START CONFIGURATION ** C ** ** C ** LINE CENTRES ARE TRANSLATIONALLY DISORDERED, AND ALL THE ** C ** LINE DIRECTIONS ARE SELECTED RANDOMLY. ** C ** WE ENSURE NO OVERLAPS DURING CONSTRUCTION. ** C ******************************************************************* INTEGER N PARAMETER ( N = 25 ) REAL RX(N), RY(N), EX(N), EY(N) REAL BOX REAL TWOPI PARAMETER ( TWOPI = 2.0 * 3.1415927 ) INTEGER I LOGICAL OVRLAP REAL RANF, DUMMY, TH C ******************************************************************* DO 100 I = 1, N OVRLAP = .TRUE. 50 IF ( OVRLAP ) THEN RX(I) = ( RANF ( DUMMY ) - 0.5 ) * BOX RY(I) = ( RANF ( DUMMY ) - 0.5 ) * BOX TH = ( RANF ( DUMMY ) - 0.5 ) * TWOPI EX(I) = COS ( TH ) EY(I) = SIN ( TH ) OVRLAP = .FALSE. CALL DNCHEK ( BOX, I, OVRLAP ) GOTO 50 ENDIF 100 CONTINUE RETURN END SUBROUTINE EQUIL ( NSWEEP, BOX, NEWBOX, : MAXDIS, MAXROT, RATDIS, RATROT ) COMMON / BLOCK1 / RX, RY, EX, EY C ******************************************************************* C ** CARRIES OUT NSWEEP SWEEPS OF EQUILIBRATION. ** C ** ** C ** ATTEMPTS TO CHANGE BOX SIZE FROM BOX TO NEWBOX THROUGHOUT. ** C ******************************************************************* INTEGER N PARAMETER ( N = 25 ) REAL RX(N), RY(N), EX(N), EY(N) INTEGER NSWEEP REAL BOX, NEWBOX, MAXDIS, MAXROT, RATDIS, RATROT LOGICAL SCALED, OVRLAP INTEGER SWEEP REAL FACBOX, TRYBOX C ******************************************************************* FACBOX = 1.0 / 2.0 SCALED = .FALSE. DO 1000 SWEEP = 1, NSWEEP C ** CONDUCT USUAL MONTE CARLO MOVE SWEEP ** CALL MOVE ( BOX, MAXDIS, MAXROT, RATDIS, RATROT ) C ** NOW ATTEMPT TO SCALE BOX ** IF ( .NOT. SCALED ) THEN C ** FIRST ATTEMPT COMPLETE SCALING ** CALL SCALE ( BOX, NEWBOX ) CALL CHECK ( NEWBOX, OVRLAP ) IF ( OVRLAP ) THEN C ** IT FAILED: UNDO SCALING ** SCALED = .FALSE. CALL SCALE ( NEWBOX, BOX ) C ** MAKE ONE ATTEMPT AT PARTIAL SCALING ** TRYBOX = BOX + FACBOX * ( NEWBOX - BOX ) CALL SCALE ( BOX, TRYBOX ) CALL CHECK ( TRYBOX, OVRLAP ) IF ( OVRLAP ) THEN C ** IT FAILED: UNDO SCALING AND ** C ** DECREASE FACBOX FOR NEXT TIME ** CALL SCALE ( TRYBOX, BOX ) FACBOX = FACBOX / 2.0 ELSE C ** IT WORKED: STICK WITH IT AND ** C ** INCREASE FACBOX FOR NEXT TIME ** BOX = TRYBOX IF ( FACBOX .LT. 0.49 ) FACBOX = FACBOX * 2.0 ENDIF ELSE C ** IT WORKED FIRST TIME ** SCALED = .TRUE. BOX = NEWBOX ENDIF ENDIF 1000 CONTINUE IF ( .NOT. SCALED ) THEN WRITE(*,'('' **** FAILED TO ACHIEVE NEW DENSITY ****'')') ENDIF NEWBOX = BOX RETURN END SUBROUTINE PRODUC ( NSWEEP, BOX, AVORD, : MAXDIS, MAXROT, RATDIS, RATROT ) COMMON / BLOCK1 / RX, RY, EX, EY C ******************************************************************* C ** CARRIES OUT NSWEEP SWEEPS OF PRODUCTION. ** C ******************************************************************* INTEGER N PARAMETER ( N = 25 ) REAL RX(N), RY(N), EX(N), EY(N) INTEGER NSWEEP REAL BOX, AVORD, MAXDIS, MAXROT, RATDIS, RATROT INTEGER SWEEP REAL ORD C ******************************************************************* AVORD = 0.0 DO 1000 SWEEP = 1, NSWEEP CALL MOVE ( BOX, MAXDIS, MAXROT, RATDIS, RATROT ) CALL ORDER ( ORD ) AVORD = AVORD + ORD 1000 CONTINUE AVORD = AVORD / REAL ( NSWEEP ) RETURN END SUBROUTINE MOVE ( BOX, MAXDIS, MAXROT, RATDIS, RATROT ) COMMON / BLOCK1 / RX, RY, EX, EY C ******************************************************************* C ** ATTEMPTS ONE TRANSLATIONAL OR ONE ROTATIONAL MOVE PER LINE ** C ******************************************************************* INTEGER N PARAMETER ( N = 25 ) REAL RX(N), RY(N), EX(N), EY(N) REAL MAXROT, MAXDIS, BOX, RATDIS, RATROT REAL TWOPI PARAMETER ( TWOPI = 2.0 * 3.1415927 ) LOGICAL ROTMOV, OVRLAP INTEGER I, TRYROT, TRYDIS, MAKROT, MAKDIS REAL RANF, DUMMY, TH, COSTH, SINTH, DRX, DRY REAL RXOLD, RYOLD, EXOLD, EYOLD REAL MINDIS, MINROT PARAMETER ( MINDIS = 0.002, MINROT = 0.002 ) C ******************************************************************* TRYROT = 0 TRYDIS = 0 MAKROT = 0 MAKDIS = 0 C ** LOOP OVER ALL LINES ** DO 1000 I = 1, N ROTMOV = ( RANF ( DUMMY ) .LT. 0.5 ) IF ( ROTMOV ) THEN C ** ATTEMPT ROTATIONAL MOVE ** TRYROT = TRYROT + 1 TH = ( RANF ( DUMMY ) - 0.5 ) * MAXROT COSTH = COS ( TH ) SINTH = SIN ( TH ) EXOLD = EX(I) EYOLD = EY(I) EX(I) = EXOLD * COSTH - EYOLD * SINTH EY(I) = EYOLD * COSTH + EXOLD * SINTH OVRLAP = .FALSE. CALL DNCHEK ( BOX, I, OVRLAP ) CALL UPCHEK ( BOX, I, OVRLAP ) IF ( OVRLAP ) THEN EX(I) = EXOLD EY(I) = EYOLD ELSE MAKROT = MAKROT + 1 ENDIF ELSE C ** ATTEMPT DISPLACEMENT ** TRYDIS = TRYDIS + 1 DRX = ( RANF ( DUMMY ) - 0.5 ) * MAXDIS DRY = ( RANF ( DUMMY ) - 0.5 ) * MAXDIS RXOLD = RX(I) RYOLD = RY(I) RX(I) = RX(I) + DRX RY(I) = RY(I) + DRY OVRLAP = .FALSE. CALL DNCHEK ( BOX, I, OVRLAP ) CALL UPCHEK ( BOX, I, OVRLAP ) IF ( OVRLAP ) THEN RX(I) = RXOLD RY(I) = RYOLD ELSE MAKDIS = MAKDIS + 1 ENDIF ENDIF 1000 CONTINUE C ** ADJUST MAXIMUM DISPLACEMENTS FOR NEXT TIME ** IF ( TRYDIS .GT. 0 ) THEN RATDIS = REAL ( MAKDIS ) / REAL ( TRYDIS ) IF ( RATDIS .GT. 0.5 ) THEN MAXDIS = MAXDIS + MINDIS ELSE MAXDIS = MAXDIS - MINDIS ENDIF IF ( MAXDIS .GT. BOX ) MAXDIS = BOX IF ( MAXDIS .LT. MINDIS ) MAXDIS = MINDIS ENDIF IF ( TRYROT .GT. 0 ) THEN RATROT = REAL(MAKROT) / REAL(TRYROT) IF ( RATROT .GT. 0.5 ) THEN MAXROT = MAXROT + MINROT ELSE MAXROT = MAXROT - MINROT ENDIF IF ( MAXROT .GT. TWOPI ) MAXROT = TWOPI IF ( MAXROT .LT. MINROT ) MAXROT = MINROT ENDIF RETURN END SUBROUTINE UPCHEK ( BOX, I, OVRLAP ) COMMON / BLOCK1 / RX, RY, EX, EY C ******************************************************************* C ** ROUTINE TO CHECK FOR OVERLAPS WITH LINES J>I ** C ******************************************************************* INTEGER N PARAMETER ( N = 25 ) REAL RX(N), RY(N), EX(N), EY(N) LOGICAL OVRLAP INTEGER I REAL BOX REAL TOL PARAMETER ( TOL = 1.0E-6 ) INTEGER J REAL RXI, RYI, EXI, EYI REAL RXIJ, RYIJ, RIJSQ, EXJ, EYJ, DET, DI, DJ REAL BOXINV C ******************************************************************* BOXINV = 1.0 / BOX RXI = RX(I) RYI = RY(I) EXI = EX(I) EYI = EY(I) J = I + 1 10 IF ( ( .NOT. OVRLAP ) .AND. ( J .LE. N ) ) THEN RXIJ = RXI - RX(J) RYIJ = RYI - RY(J) RXIJ = RXIJ - ANINT ( RXIJ * BOXINV ) * BOX RYIJ = RYIJ - ANINT ( RYIJ * BOXINV ) * BOX RIJSQ = RXIJ ** 2 + RYIJ ** 2 IF ( RIJSQ .LT. 1.0 ) THEN EXJ = EX(J) EYJ = EY(J) DET = EXI * EYJ - EXJ * EYI IF ( ABS ( DET ) .GT. TOL ) THEN DI = ( EXJ * RYIJ - EYJ * RXIJ ) / DET DJ = ( EXI * RYIJ - EYI * RXIJ ) / DET OVRLAP = ( ABS ( DI ) .LT. 0.5 ) .AND. : ( ABS ( DJ ) .LT. 0.5 ) ENDIF ENDIF J = J + 1 GOTO 10 ENDIF RETURN END SUBROUTINE DNCHEK ( BOX, I, OVRLAP ) COMMON / BLOCK1 / RX, RY, EX, EY C ******************************************************************* C ** ROUTINE TO CHECK FOR OVERLAPS WITH LINES J )" A$ 10DIM X(49),Y(49),T(49),CA(49),SA(49),XS1(49),YS1(49),XS2(49),YS2(49) 15MODE4 16@%=&20208 20N=49 30NR=7 40INPUT "BOX LENGTH (400-700)";LB 50LB2=LB/2 60INPUT "LINE LENGTH ";L 61INPUT "NUMBER OF CYCLES (1000)"; CYCLE 62PRINT'''TAB(5);"BLEEPS FOR TRIAL OVERLAP" 63 PRINT'TAB(5);"OUTPUTS ORDER PARAMETER" 64 PRINTTAB(5);"EVERY FOUR CYCLES" 68 INPUTTAB(5,31)"( press )" A$ 69LSQ=L*L 70L2=L/2 80DP=0 90NA=0 100DMAX=50 110TMAX=PI/4 120KSET=4*N 130MAXKB=CYCLE*N 134MODE1 136COLOUR2 137GCOL0,1 140PROCBOX 150PROCSETUP 160I=0 170REPEAT 180KB=KB+1 190I=I+1 200IF I>N THEN I=1 210DX=(2*RND(1)-1)*DMAX 220DY=(2*RND(1)-1)*DMAX 230DT=(2*RND(1)-1)*TMAX 240XI=X(I)+DX 250YI=Y(I)+DY 260TI=T(I)+DT 270IF XI>LB THEN XI=XI-LB ELSE IF XI<0 THEN XI=XI+LB 280IF YI>LB THEN YI=YI-LB ELSE IF YI<0 THEN YI=YI+LB 290IF TI>PI THEN TI=TI-PI ELSE IF TI<0 THEN TI=TI+PI 300CI=COS(TI) 310SI=SQR(1-CI^2) 320PROCMOVE 325K=0 330K=K+1 340IF K=I THEN 465 350XDIFF=X(K)-XI 360YDIFF=Y(K)-YI 370XDIFF=XDIFF-LB*(XDIFF DIV LB2) 380YDIFF=YDIFF-LB*(YDIFF DIV LB2) 390RSQ=XDIFF^2+YDIFF^2 400IF RSQ>LSQ THEN 465 410S1=(CA(K)*SI-SA(K)*CI)*L 420LM1=(YDIFF*CA(K)-XDIFF*SA(K))/S1 430IF ABS(LM1)>0.5 THEN 465 440LM2=(YDIFF*CI-XDIFF*SI)/S1 450IF ABS(LM2)>0.5 THEN 465 451SOUND 1,-10,53,5 455PROCMOVEBACK 460GOTO 540 465REM 470IF K