CCL Home Page
Up Directory CCL f.31
C    *******************************************************************
C    ** THIS FORTRAN CODE IS INTENDED TO ILLUSTRATE POINTS MADE IN    **
C    ** THE TEXT. TO OUR KNOWLEDGE IT WORKS CORRECTLY. HOWEVER IT IS  **
C    ** THE RESPONSIBILITY OF THE USER TO TEST IT, IF IT IS USED IN A **
C    ** RESEARCH APPLICATION.                                         **
C    *******************************************************************

C    *******************************************************************
C    ** FICHE F.31                                                    **
C    ** CONSTANT-NPT MOLECULAR DYNAMICS USING CONSTRAINT ALGORITHM.   **
C    *******************************************************************

        PROGRAM EVAMOR

        COMMON / BLOCK1 / RX, RY, RZ, RX1, RY1, RZ1,
     :                    RX2, RY2, RZ2, RX3, RY3, RZ3,
     :                    PX, PY, PZ, PX1, PY1, PZ1,
     :                    PX2, PY2, PZ2, PX3, PY3, PZ3,
     :                    FX, FY, FZ
        COMMON / BLOCK2 / BOX, BOX1, BOX2, BOX3

C    *******************************************************************
C    ** CONSTANT-NPT MOLECULAR DYNAMICS USING CONSTRAINT ALGORITHM.   **
C    **                                                               **
C    ** THE MODIFIED EQUATIONS OF MOTION ARE AS FOLLOWS:              **
C    **       R1 = P/M + CHI*R                                        **
C    **       P1 = F - CHI*P - XI*P                                   **
C    **       V1 = 3*V*CHI                                            **
C    ** WHERE R1 IS THE TIME DERIVATIVE OF POSITION R                 **
C    ** P1 IS THE TIME DERIVATIVE OF MOMENTUM P                       **
C    ** V1 IS THE TIME DERIVATIVE OF VOLUME V                         **
C    ** AND CHI AND XI ARE LAGRANGE MULTIPLIERS                       **
C    **       CHI = ( - SUM( (R.P) X / (R.R) ) ) / ( SUM(X) + 9PV )   **
C    **       XI  = SUM(P.F) / SUM(P.P) - CHI                         **
C    ** HERE TERMS LIKE (P.F) ARE SCALAR PRODUCTS OVER PAIR TERMS AND **
C    ** PV STANDS FOR PRESSURE * VOLUME.                              **
C    ** WE SOLVE THESE EQUATIONS BY A GEAR 4-VALUE METHOD FOR FIRST   **
C    ** ORDER DIFFERENTIAL EQUATIONS                                  **
C    **                                                               **
C    ** REFERENCES:                                                   **
C    **                                                               **
C    ** EVANS AND MORRISS, CHEM PHYS 77, 63, 1983.                    **
C    ** EVANS AND MORRISS, COMPUT PHYS REP 1, 297, 1984.              **
C    **                                                               **
C    ** PRINCIPAL VARIABLES:                                          **
C    **                                                               **
C    ** INTEGER N                             NUMBER OF MOLECULES     **
C    ** REAL    DT                            TIMESTEP                **
C    ** REAL    RX(N),RY(N),RZ(N)             POSITIONS               **
C    ** REAL    RX1(N),RY1(N),RZ1(N)          FIRST DERIVATIVES       **
C    ** REAL    RX2(N),RY2(N),RZ2(N)          SECOND DERIVATIVES      **
C    ** REAL    RX3(N),RY3(N),RZ3(N)          THIRD DERIVATIVES       **
C    ** REAL    PX(N),PY(N),PZ(N)             MOMENTA                 **
C    ** REAL    PX1(N),PY1(N),PZ1(N)          FIRST DERIVATIVES       **
C    ** REAL    PX2(N),PY2(N),PZ2(N)          SECOND DERIVATIVES      **
C    ** REAL    PX3(N),PY3(N),PZ3(N)          THIRD DERIVATIVES       **
C    ** REAL    BOX,BOX1,BOX2,BOX3            BOX LENGTH AND DERIVS   **
C    ** REAL    FX(N),FY(N),FZ(N)             TOTAL FORCES            **
C    **                                                               **
C    ** ROUTINES REFERENCED                                           **
C    **                                                               **
C    ** SUBROUTINE READCN ( CNFILE )                                  **
C    **    READS IN CONFIGURATION AND BOX VARIABLES                   **
C    ** SUBROUTINE FORCE ( RCUT, V, W, X, RPX, PF )                   **
C    **    CALCULATES FORCES, POTENTIAL, VIRIAL, HYPERVIRIAL ETC.     **
C    ** SUBROUTINE KINET ( K )                                        **
C    **    CALCULATES KINETIC ENERGY                                  **
C    ** SUBROUTINE PREDIC ( DT )                                      **
C    **    PREDICTOR ROUTINE FOR CONFIGURATION AND BOX VARIABLES      **
C    ** SUBROUTINE CORREC ( DT, CHI, CHIPXI )                         **
C    **    CORRECTOR ROUTINE FOR CONFIGURATION AND BOX VARIABLES      **
C    ** SUBROUTINE WRITCN ( CNFILE )                                  **
C    **    WRITES OUT CONFIGURATION AND BOX VARIABLES                 **
C    ** SUBROUTINE SCALE ( TEMPER, PRESUR, WLRC0, XLRC0, RCUT )       **
C    **    SCALES MOMENTA AND BOX SIZE TO GIVE DESIRED T AND P        **
C    *******************************************************************

        INTEGER     N
        PARAMETER ( N = 108 )

        REAL        RX(N), RY(N), RZ(N)
        REAL        RX1(N), RY1(N), RZ1(N)
        REAL        RX2(N), RY2(N), RZ2(N)
        REAL        RX3(N), RY3(N), RZ3(N)
        REAL        PX(N), PY(N), PZ(N)
        REAL        PX1(N), PY1(N), PZ1(N)
        REAL        PX2(N), PY2(N), PZ2(N)
        REAL        PX3(N), PY3(N), PZ3(N)
        REAL        FX(N), FY(N), FZ(N)
        REAL        BOX, BOX1, BOX2, BOX3

        INTEGER     STEP, NSTEP, IPRINT, ISCALE
        REAL        ACV, ACK, ACE, ACP, ACT, ACH, ACD
        REAL        ACVSQ, ACKSQ, ACESQ, ACPSQ, ACTSQ, ACHSQ, ACDSQ
        REAL        AVV, AVK, AVE, AVP, AVT, AVH, AVD
        REAL        FLV, FLK, FLE, FLP, FLT, FLH, FLD
        REAL        DT, DENS, TEMP, RCUT, VOL, PRES, H, NORM
        REAL        TEMPER, PRESUR
        REAL        K, V, W, E
        REAL        KN, VN, EN, HN
        REAL        X, RPX, PF, PP, PV, CHI, CHIPXI
        REAL        SR3, SR9, VLRC, WLRC, XLRC, VLRC0, WLRC0, XLRC0, PI
        CHARACTER   TITLE*80, CNFILE*30
        REAL        FREE

        PARAMETER ( FREE = 3.0 )
        PARAMETER ( PI = 3.1415927 )

C    *******************************************************************

        WRITE(*,'(1H1,'' **** PROGRAM EVAMOR ****                  '')')
        WRITE(*,'(//1X,''MOLECULAR DYNAMICS OF LENNARD-JONES ATOMS '')')
        WRITE(*,'(1X,''CONSTANT-NPT ALGORITHM OF EVANS AND MORRISS '')')

C    ** BASIC SIMULATION PARAMETERS **

        WRITE(*,'('' ENTER RUN TITLE                               '')')
        READ (*,'(A)') TITLE
        WRITE(*,'('' ENTER CONFIGURATION FILENAME                  '')')
        READ (*,'(A)') CNFILE
        WRITE(*,'('' ENTER NUMBER OF STEPS                         '')')
        READ (*,*) NSTEP
        WRITE(*,'('' ENTER INTERVAL BETWEEN PRINTS                 '')')
        READ (*,*) IPRINT
        WRITE(*,'('' ENTER INTERVAL BETWEEN BOX AND MOMENTUM SCALES'')')
        READ (*,*) ISCALE
        WRITE(*,'('' ENTER THE FOLLOWING IN L-J REDUCED UNITS      '')')
        WRITE(*,'('' ENTER TIMESTEP                                '')')
        READ (*,*) DT
        WRITE(*,'('' ENTER POTENTIAL CUTOFF                        '')')
        READ (*,*) RCUT
        WRITE(*,'('' ENTER DESIRED TEMPERATURE                     '')')
        READ (*,*) TEMPER
        WRITE(*,'('' ENTER DESIRED PRESSURE                        '')')
        READ (*,*) PRESUR

        WRITE(*,'(//1X,A)') TITLE
        WRITE(*,'('' CONFIGURATION FILENAME '',A)') CNFILE
        WRITE(*,'('' NUMBER OF STEPS  = '',I6   )') NSTEP
        WRITE(*,'('' PRINT INTERVAL   = '',I6   )') IPRINT
        WRITE(*,'('' SCALE INTERVAL   = '',I6   )') ISCALE
        WRITE(*,'('' TIMESTEP         = '',F10.5)') DT
        WRITE(*,'('' POTENTIAL CUTOFF = '',F10.5)') RCUT
        WRITE(*,'('' DESIRED TEMP.    = '',F10.5)') TEMPER
        WRITE(*,'('' DESIRED PRES.    = '',F10.5)') PRESUR

C    ** READCN MUST READ IN INITIAL CONFIGURATION    **
C    ** AND ASSIGN VALUES TO BOX AND ITS DERIVATIVES **

        CALL READCN ( CNFILE )
        VOL = BOX ** 3
        DENS = REAL ( N ) / VOL

        WRITE(*,'('' INITIAL DENS.    = '',F10.5)') DENS

        IF ( IPRINT .LE. 0 ) IPRINT = NSTEP + 1
        IF ( ISCALE .LE. 0 ) ISCALE = NSTEP + 1

C    ** PREPARE FACTORS FOR LONG-RANGE CORRECTIONS **
C    ** NB: SPECIFIC TO LENNARD-JONES POTENTIAL    **

        SR3 = ( 1.0 / RCUT ) ** 3
        SR9 = SR3 ** 3
        VLRC0 =   8.0 * PI * REAL ( N ) * SR9 / 9.0
     :          - 8.0 * PI * REAL ( N ) * SR3 / 3.0
        WLRC0 =  32.0 * PI * REAL ( N ) * SR9 / 9.0
     :         - 16.0 * PI * REAL ( N ) * SR3 / 3.0
        XLRC0 = 128.0 * PI * REAL ( N ) * SR9 / 9.0
     :         - 32.0 * PI * REAL ( N ) * SR3 / 3.0

C    ** ZERO ACCUMULATORS **

        ACV = 0.0
        ACK = 0.0
        ACE = 0.0
        ACP = 0.0
        ACT = 0.0
        ACH = 0.0
        ACD = 0.0

        ACVSQ = 0.0
        ACKSQ = 0.0
        ACESQ = 0.0
        ACPSQ = 0.0
        ACTSQ = 0.0
        ACHSQ = 0.0
        ACDSQ = 0.0

        FLV = 0.0
        FLK = 0.0
        FLE = 0.0
        FLP = 0.0
        FLT = 0.0
        FLH = 0.0
        FLD = 0.0

        WRITE(*,'(//1X,''**** START OF DYNAMICS ****'')')
        WRITE(*,10001)

C    *******************************************************************
C    ** MAIN LOOP STARTS                                              **
C    *******************************************************************

        DO 1000 STEP = 1, NSTEP

C       ** IMPLEMENT ALGORITHM **

           CALL PREDIC ( DT )

           CALL FORCE ( RCUT, V, W, X, RPX, PF )
           CALL KINET ( K )

C       ** INCLUDE LONG-RANGE CORRECTIONS IN ALGORITHM **

           KN   = K / REAL ( N )
           VOL  = BOX ** 3
           DENS = REAL ( N ) / VOL
           WLRC = WLRC0 * DENS
           XLRC = XLRC0 * DENS
           TEMP = 2.0 * KN / FREE
           PRES = DENS * TEMP + ( W + WLRC ) / VOL
           PV   = PRES * VOL
           PP   = 2.0 * K

           CHI    = - RPX / 9.0 / ( X + XLRC + PV )
           CHIPXI = PF / PP

           CALL CORREC ( DT, CHI, CHIPXI )

C       ** CALCULATE INSTANTANEOUS VALUES   **
C       ** INCLUDING LONG-RANGE CORRECTIONS **

           VOL  = BOX ** 3
           DENS = REAL ( N ) / VOL
           VLRC = VLRC0 * DENS
           WLRC = WLRC0 * DENS
           XLRC = XLRC0 * DENS
           V    = V + VLRC
           W    = W + WLRC
           X    = X + XLRC
           E    = K + V
           VN   = V / REAL ( N )
           EN   = E / REAL ( N )
           TEMP = 2.0 * KN / FREE
           PRES = DENS * TEMP + W / VOL
           H    = E + PRES * VOL
           HN   = H / REAL ( N )

C       ** INCREMENT ACCUMULATORS **

           ACE = ACE + EN
           ACK = ACK + KN
           ACV = ACV + VN
           ACP = ACP + PRES
           ACH = ACH + HN
           ACD = ACD + DENS

           ACESQ = ACESQ + EN ** 2
           ACKSQ = ACKSQ + KN ** 2
           ACVSQ = ACVSQ + VN ** 2
           ACPSQ = ACPSQ + PRES ** 2
           ACHSQ = ACHSQ + HN ** 2
           ACDSQ = ACDSQ + DENS ** 2

C       ** OPTIONALLY PRINT INFORMATION **

           IF ( MOD ( STEP, IPRINT ) .EQ. 0 ) THEN

              WRITE(*,'(1X,I8,7(2X,F10.5))')
     :                 STEP, EN, HN, KN, VN, PRES, TEMP, DENS

           ENDIF

C       ** OPTIONALLY SCALE MOMENTA AND BOX SIZE **

           IF ( MOD ( STEP, ISCALE ) .EQ. 0 ) THEN

              CALL SCALE ( TEMPER, PRESUR, WLRC0, XLRC0, RCUT )

           ENDIF

1000    CONTINUE

C    *******************************************************************
C    ** MAIN LOOP ENDS                                                **
C    *******************************************************************

        WRITE(*,'(/1X,''**** END OF DYNAMICS **** ''//)')

C    ** WRITE OUT FINAL CONFIGURATION **
C    ** INCLUDING BOX,BOX1,BOX2,BOX3  **

        CALL WRITCN ( CNFILE )

C    ** WRITE OUT FINAL AVERAGES **

        NORM = REAL ( NSTEP )
        AVE  = ACE / NORM
        AVK  = ACK / NORM
        AVV  = ACV / NORM
        AVP  = ACP / NORM
        AVH  = ACH / NORM
        AVD  = ACD / NORM

        ACESQ = ( ACESQ / NORM ) - AVE ** 2
        ACKSQ = ( ACKSQ / NORM ) - AVK ** 2
        ACVSQ = ( ACVSQ / NORM ) - AVV ** 2
        ACPSQ = ( ACPSQ / NORM ) - AVP ** 2
        ACHSQ = ( ACHSQ / NORM ) - AVH ** 2
        ACDSQ = ( ACDSQ / NORM ) - AVD ** 2

        IF ( ACESQ .GT. 0.0 ) FLE = SQRT ( ACESQ )
        IF ( ACKSQ .GT. 0.0 ) FLK = SQRT ( ACKSQ )
        IF ( ACVSQ .GT. 0.0 ) FLV = SQRT ( ACVSQ )
        IF ( ACPSQ .GT. 0.0 ) FLP = SQRT ( ACPSQ )
        IF ( ACHSQ .GT. 0.0 ) FLH = SQRT ( ACHSQ )
        IF ( ACDSQ .GT. 0.0 ) FLD = SQRT ( ACDSQ )

        AVT = AVK / 1.5
        FLT = FLK / 1.5

        WRITE(*,'('' AVERAGES'',7(2X,F10.5))')
     :             AVE, AVH, AVK, AVV, AVP, AVT, AVD
        WRITE(*,'('' FLUCTS  '',7(2X,F10.5))')
     :             FLE, FLH, FLK, FLV, FLP, FLT, FLD

        STOP

10001   FORMAT(//1X,'TIMESTEP  ..ENERGY..  .ENTHALPY.  ..KINETIC.',
     :          '  ..POTENT..  .PRESSURE.  ..TEMPER..  ..DENSITY.')
        END



        SUBROUTINE FORCE ( RCUT, V, W, X, RPX, PF )

        COMMON / BLOCK1 / RX, RY, RZ, RX1, RY1, RZ1,
     :                    RX2, RY2, RZ2, RX3, RY3, RZ3,
     :                    PX, PY, PZ, PX1, PY1, PZ1,
     :                    PX2, PY2, PZ2, PX3, PY3, PZ3,
     :                    FX, FY, FZ
        COMMON / BLOCK2 / BOX, BOX1, BOX2, BOX3

C    *******************************************************************
C    ** LENNARD-JONES FORCE ROUTINE IN REDUCED UNITS                  **
C    **                                                               **
C    ** THE POTENTIAL IS V(R) = 4*((1/R)**12-(1/R)**6)                **
C    ** WE INCLUDE SPHERICAL CUTOFF AND MINIMUM IMAGING IN CUBIC BOX. **
C    **                                                               **
C    ** PRINCIPAL VARIABLES:                                          **
C    **                                                               **
C    ** INTEGER N                 NUMBER OF MOLECULES                 **
C    ** REAL    RX(N),RY(N),RZ(N) MOLECULAR POSITIONS                 **
C    ** REAL    PX(N),PY(N),PZ(N) MOLECULAR MOMENTA                   **
C    ** REAL    FX(N),FY(N),FZ(N) MOLECULAR FORCES                    **
C    ** REAL    BOX               SIMULATION BOX LENGTH               **
C    ** REAL    RCUT              PAIR POTENTIAL CUTOFF               **
C    ** REAL    V                 POTENTIAL ENERGY                    **
C    ** REAL    W                 VIRIAL FUNCTION                     **
C    ** REAL    X                 HYPERVIRIAL FUNCTION                **
C    ** REAL    RPX               ADDITIONAL PAIR SUMMATION           **
C    ** REAL    PF                ANOTHER PAIR SUMMATION              **
C    *******************************************************************

        INTEGER     N
        PARAMETER ( N = 108 )

        REAL        RX(N), RY(N), RZ(N)
        REAL        RX1(N), RY1(N), RZ1(N)
        REAL        RX2(N), RY2(N), RZ2(N)
        REAL        RX3(N), RY3(N), RZ3(N)
        REAL        PX(N), PY(N), PZ(N)
        REAL        PX1(N), PY1(N), PZ1(N)
        REAL        PX2(N), PY2(N), PZ2(N)
        REAL        PX3(N), PY3(N), PZ3(N)
        REAL        FX(N), FY(N), FZ(N)
        REAL        BOX, BOX1, BOX2, BOX3

        INTEGER     I, J
        REAL        RCUT, V, W, X, RPX, PF
        REAL        BOXINV, RCUTSQ
        REAL        RXI, RYI, RZI, RXIJ, RYIJ, RZIJ, RIJSQ
        REAL        PXI, PYI, PZI, PXIJ, PYIJ, PZIJ
        REAL        FXI, FYI, FZI, FXIJ, FYIJ, FZIJ
        REAL        SR2, SR6, SR12, VIJ, WIJ, FIJ, XIJ, RP

C    *******************************************************************

C    ** USEFUL QUANTITIES **

        BOXINV = 1.0 / BOX
        RCUTSQ = RCUT ** 2

C    ** ZERO FORCES, POTENTIAL, VIRIAL **

        DO 100 I = 1, N

           FX(I) = 0.0
           FY(I) = 0.0
           FZ(I) = 0.0

100     CONTINUE

        V   = 0.0
        W   = 0.0
        X   = 0.0
        RPX = 0.0
        PF  = 0.0

C    ** OUTER LOOP BEGINS **

        DO 200 I = 1, N - 1

           RXI = RX(I)
           RYI = RY(I)
           RZI = RZ(I)
           FXI = FX(I)
           FYI = FY(I)
           FZI = FZ(I)
           PXI = PX(I)
           PYI = PY(I)
           PZI = PZ(I)

C       ** INNER LOOP BEGINS **

           DO 199 J = I + 1, N

              RXIJ = RXI - RX(J)
              RYIJ = RYI - RY(J)
              RZIJ = RZI - RZ(J)
              RXIJ = RXIJ - ANINT ( RXIJ * BOXINV ) * BOX
              RYIJ = RYIJ - ANINT ( RYIJ * BOXINV ) * BOX
              RZIJ = RZIJ - ANINT ( RZIJ * BOXINV ) * BOX
              RIJSQ  = RXIJ ** 2 + RYIJ ** 2 + RZIJ ** 2

              IF ( RIJSQ .LT. RCUTSQ ) THEN

                 SR2   = 1.0 / RIJSQ
                 SR6   = SR2 * SR2 * SR2
                 SR12  = SR6 ** 2
                 VIJ   = SR12 - SR6
                 V     = V + VIJ
                 WIJ   = VIJ + SR12
                 W     = W + WIJ
                 FIJ   = WIJ * SR2
                 FXIJ  = FIJ * RXIJ
                 FYIJ  = FIJ * RYIJ
                 FZIJ  = FIJ * RZIJ
                 FXI   = FXI + FXIJ
                 FYI   = FYI + FYIJ
                 FZI   = FZI + FZIJ
                 FX(J) = FX(J) - FXIJ
                 FY(J) = FY(J) - FYIJ
                 FZ(J) = FZ(J) - FZIJ
                 XIJ   = WIJ + 2.0 * SR12
                 X     = X + XIJ

                 PXIJ  = PXI - PX(J)
                 PYIJ  = PYI - PY(J)
                 PZIJ  = PZI - PZ(J)
                 RP    = RXIJ * PXIJ + RYIJ * PYIJ + RZIJ * PZIJ
                 RPX   = RPX + RP * XIJ * SR2
                 PF    = PF + PXIJ * FXIJ + PYIJ * FYIJ + PZIJ * FZIJ

              ENDIF

199        CONTINUE

C       ** INNER LOOP ENDS **

           FX(I) = FXI
           FY(I) = FYI
           FZ(I) = FZI

200     CONTINUE

C    ** OUTER LOOP ENDS **

C    ** MULTIPLY RESULTS BY NUMERICAL FACTORS **

        DO 300 I = 1, N

           FX(I) = FX(I) * 24.0
           FY(I) = FY(I) * 24.0
           FZ(I) = FZ(I) * 24.0

300     CONTINUE

        V   = V * 4.0
        W   = W * 24.0 / 3.0
        X   = X * 144.0 / 9.0
        RPX = RPX * 144.0
        PF  = PF * 24.0

        RETURN
        END



        SUBROUTINE KINET ( K )

        COMMON / BLOCK1 / RX, RY, RZ, RX1, RY1, RZ1,
     :                    RX2, RY2, RZ2, RX3, RY3, RZ3,
     :                    PX, PY, PZ, PX1, PY1, PZ1,
     :                    PX2, PY2, PZ2, PX3, PY3, PZ3,
     :                    FX, FY, FZ
        COMMON / BLOCK2 / BOX, BOX1, BOX2, BOX3

C    *******************************************************************
C    ** ROUTINE TO COMPUTE KINETIC ENERGY                             **
C    *******************************************************************

        INTEGER     N
        PARAMETER ( N = 108 )

        REAL        RX(N), RY(N), RZ(N)
        REAL        RX1(N), RY1(N), RZ1(N)
        REAL        RX2(N), RY2(N), RZ2(N)
        REAL        RX3(N), RY3(N), RZ3(N)
        REAL        PX(N), PY(N), PZ(N)
        REAL        PX1(N), PY1(N), PZ1(N)
        REAL        PX2(N), PY2(N), PZ2(N)
        REAL        PX3(N), PY3(N), PZ3(N)
        REAL        FX(N), FY(N), FZ(N)
        REAL        BOX, BOX1, BOX2, BOX3

        REAL        K
        INTEGER     I

C    *******************************************************************

        K = 0.0

        DO 1000 I = 1, N

           K = K + PX(I) ** 2 + PY(I) ** 2 + PZ(I) ** 2

1000    CONTINUE

        K = 0.5 * K

        RETURN
        END



        SUBROUTINE READCN ( CNFILE )

        COMMON / BLOCK1 / RX, RY, RZ, RX1, RY1, RZ1,
     :                    RX2, RY2, RZ2, RX3, RY3, RZ3,
     :                    PX, PY, PZ, PX1, PY1, PZ1,
     :                    PX2, PY2, PZ2, PX3, PY3, PZ3,
     :                    FX, FY, FZ
        COMMON / BLOCK2 / BOX, BOX1, BOX2, BOX3

C    *******************************************************************
C    ** SUBROUTINE TO READ IN INITIAL CONFIGURATION FROM UNIT 10      **
C    *******************************************************************

        INTEGER     N
        PARAMETER ( N = 108 )

        REAL        RX(N), RY(N), RZ(N)
        REAL        RX1(N), RY1(N), RZ1(N)
        REAL        RX2(N), RY2(N), RZ2(N)
        REAL        RX3(N), RY3(N), RZ3(N)
        REAL        PX(N), PY(N), PZ(N)
        REAL        PX1(N), PY1(N), PZ1(N)
        REAL        PX2(N), PY2(N), PZ2(N)
        REAL        PX3(N), PY3(N), PZ3(N)
        REAL        FX(N), FY(N), FZ(N)
        REAL        BOX, BOX1, BOX2, BOX3

        CHARACTER   CNFILE*(*)

        INTEGER     CNUNIT, NN
        PARAMETER ( CNUNIT = 10 )

C     ******************************************************************

        OPEN ( UNIT = CNUNIT, FILE = CNFILE,
     :         STATUS = 'OLD', FORM = 'UNFORMATTED' )

        READ ( CNUNIT ) NN, BOX, BOX1, BOX2, BOX3
        IF ( NN .NE. N ) STOP 'INCORRECT VALUE OF N'
        READ ( CNUNIT ) RX, RY, RZ
        READ ( CNUNIT ) RX1, RY1, RZ1
        READ ( CNUNIT ) RX2, RY2, RZ2
        READ ( CNUNIT ) RX3, RY3, RZ3
        READ ( CNUNIT ) PX, PY, PZ
        READ ( CNUNIT ) PX1, PY1, PZ1
        READ ( CNUNIT ) PX2, PY2, PZ2
        READ ( CNUNIT ) PX3, PY3, PZ3

        CLOSE ( UNIT = CNUNIT )

        RETURN
        END



        SUBROUTINE WRITCN ( CNFILE )

        COMMON / BLOCK1 / RX, RY, RZ, RX1, RY1, RZ1,
     :                    RX2, RY2, RZ2, RX3, RY3, RZ3,
     :                    PX, PY, PZ, PX1, PY1, PZ1,
     :                    PX2, PY2, PZ2, PX3, PY3, PZ3,
     :                    FX, FY, FZ
        COMMON / BLOCK2 / BOX, BOX1, BOX2, BOX3

C    *******************************************************************
C    ** ROUTINE TO WRITE OUT FINAL CONFIGURATION TO UNIT 10           **
C    *******************************************************************

        INTEGER     N
        PARAMETER ( N = 108 )

        REAL        RX(N), RY(N), RZ(N)
        REAL        RX1(N), RY1(N), RZ1(N)
        REAL        RX2(N), RY2(N), RZ2(N)
        REAL        RX3(N), RY3(N), RZ3(N)
        REAL        PX(N), PY(N), PZ(N)
        REAL        PX1(N), PY1(N), PZ1(N)
        REAL        PX2(N), PY2(N), PZ2(N)
        REAL        PX3(N), PY3(N), PZ3(N)
        REAL        FX(N), FY(N), FZ(N)
        REAL        BOX, BOX1, BOX2, BOX3

        CHARACTER   CNFILE*(*)

        INTEGER     CNUNIT
        PARAMETER ( CNUNIT = 10 )

C    *******************************************************************

        OPEN ( UNIT = CNUNIT, FILE = CNFILE,
     :         STATUS = 'OLD', FORM = 'UNFORMATTED' )

        WRITE ( CNUNIT ) N, BOX, BOX1, BOX2, BOX3
        WRITE ( CNUNIT ) RX, RY, RZ
        WRITE ( CNUNIT ) RX1, RY1, RZ1
        WRITE ( CNUNIT ) RX2, RY2, RZ2
        WRITE ( CNUNIT ) RX3, RY3, RZ3
        WRITE ( CNUNIT ) PX, PY, PZ
        WRITE ( CNUNIT ) PX1, PY1, PZ1
        WRITE ( CNUNIT ) PX2, PY2, PZ2
        WRITE ( CNUNIT ) PX3, PY3, PZ3

        CLOSE ( UNIT = CNUNIT )

        RETURN
        END



        SUBROUTINE PREDIC ( DT )

        COMMON / BLOCK1 / RX, RY, RZ, RX1, RY1, RZ1,
     :                    RX2, RY2, RZ2, RX3, RY3, RZ3,
     :                    PX, PY, PZ, PX1, PY1, PZ1,
     :                    PX2, PY2, PZ2, PX3, PY3, PZ3,
     :                    FX, FY, FZ
        COMMON / BLOCK2 / BOX, BOX1, BOX2, BOX3

C    *******************************************************************
C    ** STANDARD TAYLOR SERIES PREDICTORS                             **
C    *******************************************************************

        INTEGER     N
        PARAMETER ( N = 108 )

        REAL        RX(N), RY(N), RZ(N)
        REAL        RX1(N), RY1(N), RZ1(N)
        REAL        RX2(N), RY2(N), RZ2(N)
        REAL        RX3(N), RY3(N), RZ3(N)
        REAL        PX(N), PY(N), PZ(N)
        REAL        PX1(N), PY1(N), PZ1(N)
        REAL        PX2(N), PY2(N), PZ2(N)
        REAL        PX3(N), PY3(N), PZ3(N)
        REAL        FX(N), FY(N), FZ(N)
        REAL        BOX, BOX1, BOX2, BOX3
        REAL        DT

        INTEGER     I
        REAL        C1, C2, C3

C    *******************************************************************

        C1 = DT
        C2 = C1 * DT / 2.0
        C3 = C2 * DT / 3.0

        DO 100 I = 1, N

           RX(I)  = RX(I)  + C1*RX1(I) + C2*RX2(I) + C3*RX3(I)
           RY(I)  = RY(I)  + C1*RY1(I) + C2*RY2(I) + C3*RY3(I)
           RZ(I)  = RZ(I)  + C1*RZ1(I) + C2*RZ2(I) + C3*RZ3(I)
           RX1(I) = RX1(I) + C1*RX2(I) + C2*RX3(I)
           RY1(I) = RY1(I) + C1*RY2(I) + C2*RY3(I)
           RZ1(I) = RZ1(I) + C1*RZ2(I) + C2*RZ3(I)
           RX2(I) = RX2(I) + C1*RX3(I)
           RY2(I) = RY2(I) + C1*RY3(I)
           RZ2(I) = RZ2(I) + C1*RZ3(I)

           PX(I)  = PX(I)  + C1*PX1(I) + C2*PX2(I) + C3*PX3(I)
           PY(I)  = PY(I)  + C1*PY1(I) + C2*PY2(I) + C3*PY3(I)
           PZ(I)  = PZ(I)  + C1*PZ1(I) + C2*PZ2(I) + C3*PZ3(I)
           PX1(I) = PX1(I) + C1*PX2(I) + C2*PX3(I)
           PY1(I) = PY1(I) + C1*PY2(I) + C2*PY3(I)
           PZ1(I) = PZ1(I) + C1*PZ2(I) + C2*PZ3(I)
           PX2(I) = PX2(I) + C1*PX3(I)
           PY2(I) = PY2(I) + C1*PY3(I)
           PZ2(I) = PZ2(I) + C1*PZ3(I)

100     CONTINUE

        BOX  = BOX  + C1*BOX1 + C2*BOX2 + C3*BOX3
        BOX1 = BOX1 + C1*BOX2 + C2*BOX3
        BOX2 = BOX2 + C1*BOX3

        RETURN
        END



        SUBROUTINE CORREC ( DT, CHI, CHIPXI )

        COMMON / BLOCK1 / RX, RY, RZ, RX1, RY1, RZ1,
     :                    RX2, RY2, RZ2, RX3, RY3, RZ3,
     :                    PX, PY, PZ, PX1, PY1, PZ1,
     :                    PX2, PY2, PZ2, PX3, PY3, PZ3,
     :                    FX, FY, FZ
        COMMON / BLOCK2 / BOX, BOX1, BOX2, BOX3

C    *******************************************************************
C    ** GEAR CORRECTOR ALGORITHM.                                     **
C    **                                                               **
C    ** FOR TIMESTEP-SCALED VARIABLES GEAR COEFFICIENTS WOULD BE AS   **
C    ** FOLLOWS (4-VALUE METHOD, 1ST-ORDER D.E.): 3/8, 1, 3/4, 1/6.   **
C    ** CHIPXI IS SHORT FOR ( CHI + XI ) (SEE CHAPTER 7).             **
C    *******************************************************************

        INTEGER     N
        PARAMETER ( N = 108 )

        REAL        RX(N), RY(N), RZ(N)
        REAL        RX1(N), RY1(N), RZ1(N)
        REAL        RX2(N), RY2(N), RZ2(N)
        REAL        RX3(N), RY3(N), RZ3(N)
        REAL        PX(N), PY(N), PZ(N)
        REAL        PX1(N), PY1(N), PZ1(N)
        REAL        PX2(N), PY2(N), PZ2(N)
        REAL        PX3(N), PY3(N), PZ3(N)
        REAL        FX(N), FY(N), FZ(N)
        REAL        BOX, BOX1, BOX2, BOX3
        REAL        DT, CHI, CHIPXI

        INTEGER     I
        REAL        C1, C2, C3, COEFF0, COEFF2, COEFF3
        REAL        CORR, CORRX, CORRY, CORRZ, CORPX, CORPY, CORPZ
        REAL        RX1I, RY1I, RZ1I, PX1I, PY1I, PZ1I

        REAL        GEAR0, GEAR2, GEAR3
        PARAMETER ( GEAR0 = 3.0 / 8.0,
     :              GEAR2 = 3.0 / 4.0,
     :              GEAR3 = 1.0 / 6.0 )

C    *******************************************************************

        C1 = DT
        C2 = C1 * DT / 2.0
        C3 = C2 * DT / 3.0

        COEFF0 = GEAR0 * C1
        COEFF2 = GEAR2 * C1 / C2
        COEFF3 = GEAR3 * C1 / C3

        DO 400 I = 1, N

           RX1I = PX(I) + CHI * RX(I)
           RY1I = PY(I) + CHI * RY(I)
           RZ1I = PZ(I) + CHI * RZ(I)
           CORRX = RX1I - RX1(I)
           CORRY = RY1I - RY1(I)
           CORRZ = RZ1I - RZ1(I)

           RX(I)  = RX(I)  + COEFF0 * CORRX
           RY(I)  = RY(I)  + COEFF0 * CORRY
           RZ(I)  = RZ(I)  + COEFF0 * CORRZ
           RX1(I) = RX1I
           RY1(I) = RY1I
           RZ1(I) = RZ1I
           RX2(I) = RX2(I) + COEFF2 * CORRX
           RY2(I) = RY2(I) + COEFF2 * CORRY
           RZ2(I) = RZ2(I) + COEFF2 * CORRZ
           RX3(I) = RX3(I) + COEFF3 * CORRX
           RY3(I) = RY3(I) + COEFF3 * CORRY
           RZ3(I) = RZ3(I) + COEFF3 * CORRZ

           PX1I = FX(I) - CHIPXI * PX(I)
           PY1I = FY(I) - CHIPXI * PY(I)
           PZ1I = FZ(I) - CHIPXI * PZ(I)
           CORPX = PX1I - PX1(I)
           CORPY = PY1I - PY1(I)
           CORPZ = PZ1I - PZ1(I)

           PX(I)  = PX(I)  + COEFF0 * CORPX
           PY(I)  = PY(I)  + COEFF0 * CORPY
           PZ(I)  = PZ(I)  + COEFF0 * CORPZ
           PX1(I) = PX1I
           PY1(I) = PY1I
           PZ1(I) = PZ1I
           PX2(I) = PX2(I) + COEFF2 * CORPX
           PY2(I) = PY2(I) + COEFF2 * CORPY
           PZ2(I) = PZ2(I) + COEFF2 * CORPZ
           PX3(I) = PX3(I) + COEFF3 * CORPX
           PY3(I) = PY3(I) + COEFF3 * CORPY
           PZ3(I) = PZ3(I) + COEFF3 * CORPZ

400     CONTINUE

        CORR = CHI * BOX - BOX1
        BOX  = BOX  + COEFF0 * CORR
        BOX1 = CHI * BOX
        BOX2 = BOX2 + COEFF2 * CORR
        BOX3 = BOX3 + COEFF3 * CORR

        RETURN
        END



        SUBROUTINE SCALE ( TEMPER, PRESUR, WLRC0, XLRC0, RCUT )

        COMMON / BLOCK1 / RX, RY, RZ, RX1, RY1, RZ1,
     :                    RX2, RY2, RZ2, RX3, RY3, RZ3,
     :                    PX, PY, PZ, PX1, PY1, PZ1,
     :                    PX2, PY2, PZ2, PX3, PY3, PZ3,
     :                    FX, FY, FZ
        COMMON / BLOCK2 / BOX, BOX1, BOX2, BOX3

C    *******************************************************************
C    ** SCALES MOMENTA AND BOX SIZE TO GIVE DESIRED VALUES OF T AND P.**
C    **                                                               **
C    ** WE USE DIRECT MOMENTUM SCALING TO GIVE THE TEMPERATURE AND A  **
C    ** NEWTON-RAPHSON METHOD FOR THE BOX SIZE.                       **
C    ** WE USE THE PRESSURE DERIVATIVE FUNCTION                       **
C    ** BOX * DP/D(BOX) = 3 * VOL * DP/D(VOL) = 3 * ( - P - X / VOL ) **
C    ** WHERE P IS THE PRESSURE AND X THE HYPERVIRIAL FUNCTION.       **
C    ** LONG-RANGE CORRECTIONS ARE TAKEN INTO ACCOUNT.                **
C    *******************************************************************

        INTEGER     N
        PARAMETER ( N = 108 )

        REAL        RX(N), RY(N), RZ(N)
        REAL        RX1(N), RY1(N), RZ1(N)
        REAL        RX2(N), RY2(N), RZ2(N)
        REAL        RX3(N), RY3(N), RZ3(N)
        REAL        PX(N), PY(N), PZ(N)
        REAL        PX1(N), PY1(N), PZ1(N)
        REAL        PX2(N), PY2(N), PZ2(N)
        REAL        PX3(N), PY3(N), PZ3(N)
        REAL        FX(N), FY(N), FZ(N)
        REAL        BOX, BOX1, BOX2, BOX3

        INTEGER     I
        REAL        TEMPER, PRESUR, WLRC0, XLRC0, RCUT
        REAL        K, KN, V, W, X, RPX, PF, FACTOR
        REAL        VOL, DENS, WLRC, XLRC, TEMP, PRES
        REAL        FREE, TOL
        PARAMETER ( FREE = 3.0, TOL = 0.01 )

C    *******************************************************************

C    ** SCALE MOMENTA **

        CALL KINET ( K )
        KN = K / REAL ( N )
        TEMP = 2.0 * KN / FREE
        FACTOR = SQRT ( TEMPER / TEMP )

        DO 100 I = 1, N

           PX(I) = PX(I) * FACTOR
           PY(I) = PY(I) * FACTOR
           PZ(I) = PZ(I) * FACTOR

100     CONTINUE

        TEMP = TEMPER

C    ** SCALE PRESSURE USING NEWTON-RAPHSON PROCEDURE **

        CALL FORCE ( RCUT, V, W, X, RPX, PF )
        VOL  = BOX ** 3
        DENS = REAL ( N ) / VOL
        WLRC = WLRC0 * DENS
        XLRC = XLRC0 * DENS
        W    = W + WLRC
        X    = X + XLRC
        PRES = DENS * TEMP + W / VOL

1000    IF ( ABS ( PRES - PRESUR ) .GT. TOL ) THEN

           FACTOR = 1.0 + ( PRES - PRESUR ) / ( PRES + X / VOL ) / 3.0

           DO 200 I = 1, N

              RX(I) = RX(I) * FACTOR
              RY(I) = RY(I) * FACTOR
              RZ(I) = RZ(I) * FACTOR

200        CONTINUE

           BOX = BOX * FACTOR

           CALL FORCE ( RCUT, V, W, X, RPX, PF )
           VOL  = BOX ** 3
           DENS = REAL ( N ) / VOL
           WLRC = WLRC0 * DENS
           XLRC = XLRC0 * DENS
           W    = W + WLRC
           X    = X + XLRC
           PRES = DENS * TEMP + W / VOL

           GOTO 1000

        ENDIF

        RETURN
        END


Modified: Fri Oct 15 16:00:00 1993 GMT
Page accessed 5796 times since Sat Aug 26 22:58:33 2000 GMT