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.4 ** C ** TWO ROUTINES THAT TOGETHER IMPLEMENT VELOCITY VERLET METHOD. ** C ** ** C ** REFERENCE: ** C ** ** C ** SWOPE ET AL., J. CHEM. PHYS. 76, 637, 1982. ** C ** ** C ** ROUTINES SUPPLIED: ** C ** ** C ** SUBROUTINE MOVEA ( DT, M ) ** C ** MOVES POSITIONS AND PARTIALLY UPDATES VELOCITIES. ** C ** SUBROUTINE MOVEB ( DT, M, K ) ** C ** COMPLETES VELOCITY MOVE AND CALCULATES KINETIC ENERGY. ** C ** ** C ** PRINCIPAL VARIABLES: ** C ** ** C ** INTEGER N NUMBER OF MOLECULES ** C ** REAL DT TIMESTEP ** C ** REAL M ATOMIC MASS ** C ** REAL K KINETIC ENERGY ** C ** REAL RX(N),RY(N),RZ(N) POSITIONS ** C ** REAL VX(N),VY(N),VZ(N) VELOCITIES ** C ** REAL FX(N),FY(N),FZ(N) FORCES ** C ** ** C ** USAGE: ** C ** ** C ** AT THE START OF A TIMESTEP, MOVEA IS CALLED TO ADVANCE THE ** C ** POSITIONS AND 'HALF-ADVANCE' THE VELOCITIES. THEN THE FORCE ** C ** ROUTINE IS CALLED, AND THIS IS FOLLOWED BY MOVEB WHICH ** C ** COMPLETES THE ADVANCEMENT OF VELOCITIES. ** C ******************************************************************* SUBROUTINE MOVEA ( DT, M ) COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ, FX, FY, FZ C ******************************************************************* C ** FIRST PART OF VELOCITY VERLET ALGORITHM ** C ** ** C ** USAGE: ** C ** ** C ** THE FIRST PART OF THE ALGORITHM IS A TAYLOR SERIES WHICH ** C ** ADVANCES POSITIONS FROM T TO T + DT AND VELOCITIES FROM ** C ** T TO T + DT/2. AFTER THIS, THE FORCE ROUTINE IS CALLED. ** C ******************************************************************* INTEGER N PARAMETER ( N = 108 ) REAL DT, M REAL RX(N), RY(N), RZ(N) REAL VX(N), VY(N), VZ(N) REAL FX(N), FY(N), FZ(N) INTEGER I REAL DT2, DTSQ2 C ******************************************************************* DT2 = DT / 2.0 DTSQ2 = DT * DT2 DO 100 I = 1, N RX(I) = RX(I) + DT * VX(I) + DTSQ2 * FX(I) / M RY(I) = RY(I) + DT * VY(I) + DTSQ2 * FY(I) / M RZ(I) = RZ(I) + DT * VZ(I) + DTSQ2 * FZ(I) / M VX(I) = VX(I) + DT2 * FX(I) / M VY(I) = VY(I) + DT2 * FY(I) / M VZ(I) = VZ(I) + DT2 * FZ(I) / M 100 CONTINUE RETURN END SUBROUTINE MOVEB ( DT, M, K ) COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ, FX, FY, FZ C ******************************************************************* C ** SECOND PART OF VELOCITY VERLET ALGORITHM ** C ** ** C ** USAGE: ** C ** ** C ** THE SECOND PART OF THE ALGORITHM ADVANCES VELOCITIES FROM ** C ** T + DT/2 TO T + DT. THIS ASSUMES THAT FORCES HAVE BEEN ** C ** COMPUTED IN THE FORCE ROUTINE AND STORED IN FX, FY, FZ. ** C ******************************************************************* INTEGER N PARAMETER ( N = 108 ) REAL DT, M, K REAL RX(N), RY(N), RZ(N) REAL VX(N), VY(N), VZ(N) REAL FX(N), FY(N), FZ(N) INTEGER I REAL DT2 C ******************************************************************* DT2 = DT / 2.0 K = 0.0 DO 200 I = 1, N VX(I) = VX(I) + DT2 * FX(I) / M VY(I) = VY(I) + DT2 * FY(I) / M VZ(I) = VZ(I) + DT2 * FZ(I) / M K = K + VX(I) ** 2 + VY(I) ** 2 + VZ(I) ** 2 200 CONTINUE K = 0.5 * M * K RETURN END