CCL Home Page
Up Directory CCL f.10
********************************************************************************
** FICHE F.10.  HARD SPHERE MOLECULAR DYNAMICS PROGRAM                        **
** 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.        **
********************************************************************************



        PROGRAM SPHERE

        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ
        COMMON / BLOCK2 / COLTIM, PARTNR

C    *******************************************************************
C    ** MOLECULAR DYNAMICS OF HARD SPHERE ATOMS.                      **
C    **                                                               **
C    ** THIS PROGRAM TAKES IN A HARD-SPHERE CONFIGURATION (POSITIONS  **
C    ** AND VELOCITIES), CHECKS FOR OVERLAPS, AND THEN CONDUCTS A     **
C    ** MOLECULAR DYNAMICS SIMULATION RUN FOR A SPECIFIED NUMBER OF   **
C    ** COLLISIONS.  THE PROGRAM IS FAIRLY EFFICIENT, BUT USES NO     **
C    ** SPECIAL NEIGHBOUR LISTS, SO IS RESTRICTED TO A SMALL NUMBER   **
C    ** OF PARTICLES (<500).  IT IS ALWAYS ASSUMED THAT COLLISIONS    **
C    ** CAN BE PREDICTED BY LOOKING AT NEAREST NEIGHBOUR PARTICLES IN **
C    ** THE MINIMUM IMAGE CONVENTION OF PERIODIC BOUNDARIES.          **
C    ** THE BOX IS TAKEN TO BE OF UNIT LENGTH.                        **
C    ** HOWEVER, RESULTS ARE GIVEN IN UNITS WHERE SIGMA=1, KT=1.      **
C    **                                                               **
C    ** PRINCIPAL VARIABLES:                                          **
C    **                                                               **
C    ** INTEGER N                    NUMBER OF ATOMS                  **
C    ** REAL    RX(N),RY(N),RZ(N)    ATOM POSITIONS                   **
C    ** REAL    VX(N),VY(N),VZ(N)    ATOM VELOCITIES                  **
C    ** REAL    COLTIM(N)            TIME TO NEXT COLLISION           **
C    ** INTEGER PARTNR(N)            COLLISION PARTNER                **
C    ** REAL    SIGMA                ATOM DIAMETER                    **
C    **                                                               **
C    ** ROUTINES REFERENCED:                                          **
C    **                                                               **
C    ** SUBROUTINE READCN ( CNFILE )                                  **
C    **    READS IN CONFIGURATION                                     **
C    ** SUBROUTINE CHECK ( SIGMA, OVRLAP, E )                         **
C    **    CHECKS CONFIGURATION AND CALCULATES ENERGY                 **
C    ** SUBROUTINE UPLIST ( SIGMA, I )                                **
C    **    SEEKS COLLISIONS WITH J>I                                  **
C    ** SUBROUTINE DNLIST ( SIGMA, I )                                **
C    **    SEEKS COLLISIONS WITH J I                         **
C    *******************************************************************

        INTEGER     N
        PARAMETER ( N = 108 )

        REAL        TIMBIG
        PARAMETER ( TIMBIG = 1.0E10 )

        INTEGER     I
        REAL        SIGMA
        REAL        RX(N), RY(N), RZ(N), VX(N), VY(N), VZ(N)
        REAL        COLTIM(N)
        INTEGER     PARTNR(N)

        INTEGER     J
        REAL        RXI, RYI, RZI, RXIJ, RYIJ, RZIJ
        REAL        VXI, VYI, VZI, VXIJ, VYIJ, VZIJ
        REAL        RIJSQ, VIJSQ, BIJ, TIJ, DISCR, SIGSQ

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

        IF ( I .EQ. N ) RETURN

        SIGSQ = SIGMA ** 2
        COLTIM(I) = TIMBIG
        RXI = RX(I)
        RYI = RY(I)
        RZI = RZ(I)
        VXI = VX(I)
        VYI = VY(I)
        VZI = VZ(I)

        DO 100 J = I + 1, N

           RXIJ = RXI - RX(J)
           RYIJ = RYI - RY(J)
           RZIJ = RZI - RZ(J)
           RXIJ = RXIJ - ANINT ( RXIJ )
           RYIJ = RYIJ - ANINT ( RYIJ )
           RZIJ = RZIJ - ANINT ( RZIJ )
           VXIJ = VXI - VX(J)
           VYIJ = VYI - VY(J)
           VZIJ = VZI - VZ(J)
           BIJ  = RXIJ * VXIJ + RYIJ * VYIJ + RZIJ * VZIJ

           IF ( BIJ .LT. 0.0 ) THEN

              RIJSQ = RXIJ ** 2 + RYIJ ** 2 + RZIJ ** 2
              VIJSQ = VXIJ ** 2 + VYIJ ** 2 + VZIJ ** 2
              DISCR = BIJ ** 2 - VIJSQ * ( RIJSQ - SIGSQ )

              IF ( DISCR .GT. 0.0 ) THEN

                 TIJ = ( -BIJ - SQRT ( DISCR ) ) / VIJSQ

                 IF ( TIJ .LT. COLTIM(I) ) THEN

                    COLTIM(I) = TIJ
                    PARTNR(I) = J

                 ENDIF

              ENDIF

           ENDIF

100     CONTINUE

        RETURN
        END



        SUBROUTINE DNLIST ( SIGMA, J )

        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ
        COMMON / BLOCK2 / COLTIM, PARTNR

C    *******************************************************************
C    ** LOOKS FOR COLLISIONS WITH ATOMS I < J                         **
C    *******************************************************************

        INTEGER     N
        PARAMETER ( N = 108 )

        REAL        TIMBIG
        PARAMETER ( TIMBIG = 1.E10 )

        INTEGER     J
        REAL        SIGMA
        REAL        RX(N), RY(N), RZ(N), VX(N), VY(N), VZ(N)
        REAL        COLTIM(N)
        INTEGER     PARTNR(N)

        INTEGER     I
        REAL        RXJ, RYJ, RZJ, RXIJ, RYIJ, RZIJ
        REAL        VXJ, VYJ, VZJ, VXIJ, VYIJ, VZIJ
        REAL        RIJSQ, VIJSQ, BIJ, TIJ, DISCR, SIGSQ

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

        IF ( J .EQ. 1 ) RETURN

        SIGSQ = SIGMA ** 2
        RXJ = RX(J)
        RYJ = RY(J)
        RZJ = RZ(J)
        VXJ = VX(J)
        VYJ = VY(J)
        VZJ = VZ(J)

        DO 100 I = 1, J - 1

           RXIJ = RX(I) - RXJ
           RYIJ = RY(I) - RYJ
           RZIJ = RZ(I) - RZJ
           RXIJ = RXIJ - ANINT ( RXIJ )
           RYIJ = RYIJ - ANINT ( RYIJ )
           RZIJ = RZIJ - ANINT ( RZIJ )
           VXIJ = VX(I) - VXJ
           VYIJ = VY(I) - VYJ
           VZIJ = VZ(I) - VZJ
           BIJ  = RXIJ * VXIJ + RYIJ * VYIJ + RZIJ * VZIJ

           IF ( BIJ .LT. 0.0 ) THEN

              RIJSQ = RXIJ ** 2 + RYIJ ** 2 + RZIJ ** 2
              VIJSQ = VXIJ ** 2 + VYIJ ** 2 + VZIJ ** 2
              DISCR = BIJ ** 2 - VIJSQ * ( RIJSQ - SIGSQ )

              IF ( DISCR .GT. 0.0 ) THEN

                 TIJ = ( - BIJ - SQRT ( DISCR ) ) / VIJSQ

                 IF ( TIJ .LT. COLTIM(I) ) THEN

                    COLTIM(I) = TIJ
                    PARTNR(I) = J

                 ENDIF

              ENDIF

           ENDIF

100     CONTINUE

        RETURN
        END




        SUBROUTINE BUMP ( SIGMA, I, J, W )

        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ

C    *******************************************************************
C    ** COMPUTES COLLISION DYNAMICS FOR PARTICLES I AND J.            **
C    **                                                               **
C    ** IT IS ASSUMED THAT I AND J ARE IN CONTACT.                    **
C    ** THE ROUTINE ALSO COMPUTES COLLISIONAL VIRIAL W.               **
C    *******************************************************************

        INTEGER     N
        PARAMETER ( N = 108 )

        INTEGER     I, J
        REAL        SIGMA, W
        REAL        RX(N), RY(N), RZ(N), VX(N), VY(N), VZ(N)

        REAL        RXIJ, RYIJ, RZIJ, FACTOR
        REAL        DELVX, DELVY, DELVZ, SIGSQ

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

        SIGSQ = SIGMA ** 2
        RXIJ = RX(I) - RX(J)
        RYIJ = RY(I) - RY(J)
        RZIJ = RZ(I) - RZ(J)
        RXIJ = RXIJ - ANINT ( RXIJ )
        RYIJ = RYIJ - ANINT ( RYIJ )
        RZIJ = RZIJ - ANINT ( RZIJ )

        FACTOR = ( RXIJ * ( VX(I) - VX(J) ) +
     :             RYIJ * ( VY(I) - VY(J) ) +
     :             RZIJ * ( VZ(I) - VZ(J) ) ) / SIGSQ
        DELVX = - FACTOR * RXIJ
        DELVY = - FACTOR * RYIJ
        DELVZ = - FACTOR * RZIJ

        VX(I) = VX(I) + DELVX
        VX(J) = VX(J) - DELVX
        VY(I) = VY(I) + DELVY
        VY(J) = VY(J) - DELVY
        VZ(I) = VZ(I) + DELVZ
        VZ(J) = VZ(J) - DELVZ

        W = DELVX * RXIJ + DELVY * RYIJ + DELVZ * RZIJ

        RETURN
        END



        SUBROUTINE READCN ( CNFILE )

        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ

C    *******************************************************************
C    ** READS CONFIGURATION FROM CHANNEL 10                           **
C    *******************************************************************

        INTEGER     N
        PARAMETER ( N = 108 )

        INTEGER     CNUNIT
        PARAMETER ( CNUNIT = 10 )

        INTEGER     NN
        CHARACTER   CNFILE*(*)
        REAL        RX(N), RY(N), RZ(N), VX(N), VY(N), VZ(N)

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

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

        READ ( CNUNIT ) NN
        IF ( NN .NE. N ) STOP 'N ERROR IN READCN'
        READ ( CNUNIT ) RX, RY, RZ
        READ ( CNUNIT ) VX, VY, VZ

        CLOSE ( UNIT = CNUNIT )

        RETURN
        END



        SUBROUTINE WRITCN ( CNFILE )

        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ

C    *******************************************************************
C    ** WRITES CONFIGURATION TO CHANNEL 10                            **
C    *******************************************************************

        INTEGER     N
        PARAMETER ( N = 108 )

        INTEGER     CNUNIT
        PARAMETER ( CNUNIT = 10 )

        CHARACTER   CNFILE*(*)
        REAL        RX(N), RY(N), RZ(N), VX(N), VY(N), VZ(N)

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

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

        WRITE ( CNUNIT ) N
        WRITE ( CNUNIT ) RX, RY, RZ
        WRITE ( CNUNIT ) VX, VY, VZ

        CLOSE ( UNIT = CNUNIT )

        RETURN
        END



Modified: Sat Jul 22 04:40:52 2000 GMT
Page accessed 10434 times since Sat Apr 17 21:34:08 1999 GMT