******************************************************************************** ** 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