http://ccl.net/cca/software/SOURCES/FORTRAN/allen-tildesley-book/f.36.shtml
CCL f.36
```********************************************************************************
** FICHE F.36.  MONTE CARLO SIMULATION OF HARD LINES IN 2D                    **
** This FORTRAN code is intended to illustrate points made in the text.       **
** To our knowledge it works correctly.  However it is the responsibility of  **
** the user to test it, if it is to be used in a research application.        **
********************************************************************************

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 ?        '')')
WRITE(*,'('' LOW-DENSITY START CONFIGURATION (.T./.F.) ?   '')')
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

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
```
 Modified: Sat Jul 22 04:41:20 2000 GMT Page accessed 8691 times since Sat Apr 17 21:34:21 1999 GMT