|
allen-tildesley-book
|
|
README,
f.00,
f.01,
f.02,
f.03,
f.04,
f.05,
f.06,
f.07,
f.08,
f.09,
f.10,
f.11,
f.12,
f.13,
f.14,
f.15,
f.16,
f.17,
f.18,
f.19,
f.20,
f.21,
f.22,
f.23,
f.24,
f.25,
f.26,
f.27,
f.28,
f.29,
f.30,
f.31,
f.32,
f.33,
f.34,
f.35,
f.36,
f.37,
older-version
|
|
|
********************************************************************************
** 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
|
