******************************************************************************** ** FICHE F.8. SHAKE ALGORITHM FOR CONSTRAINT DYNAMICS OF A CHAIN MOLECULE ** ** 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. ** ******************************************************************************** SUBROUTINE MOVE ( DT, TOL, MAXIT, NB, BOX, K, WC ) COMMON / BLOCK1 / RX, RY, RZ, OX, OY, OZ, FX, FY, FZ COMMON / BLOCK2 / DSQ, M C ******************************************************************* C ** CONSTRAINT DYNAMICS OF A CHAIN OF ATOMS USING SHAKE. ** C ** ** C ** THIS ROUTINE COMPUTES CONSTRAINT EFFECTS IN AN ITERATIVE WAY. ** C ** WE APPLY BOND LENGTH CONSTRAINTS TO ADJACENT ATOMS ONLY IN A ** C ** CHAIN MOLECULE WHICH MAY BE CYCLIC. THE GENERALIZATION TO ** C ** TO MORE COMPLICATED SYSTEMS IS STRAIGHTFORWARD. THE ** C ** CONSTRAINT EQUATIONS ARE LINEARIZED, AND EACH CONSTRAINT IS ** C ** TREATED IN TURN, UNTIL BOND LENGTHS ARE SATISFIED TO WITHIN A ** C ** SPECIFIED TOLERANCE. ** C ** IN THIS EXAMPLE WE TAKE A 6-ATOM CHAIN. ** C ** ** C ** REFERENCE: ** C ** ** C ** RYCKAERT ET AL., J. COMP. PHYS. 23, 327, 1977. ** C ** ** C ** PRINCIPAL VARIABLES: ** C ** ** C ** INTEGER N NUMBER OF MOLECULES ** C ** INTEGER NA NUMBER OF ATOMS PER MOL. ** C ** REAL DT TIMESTEP ** C ** REAL TOL BOND LENGTH TOLERANCE ** C ** INTEGER MAXIT MAXIMUM ALLOWED ITERATIONS ** C ** INTEGER NB NUMBER OF BONDS ** C ** REAL BOX BOX LENGTH ** C ** REAL K KINETIC ENERGY ** C ** REAL WC CONSTRAINT VIRIAL ** C ** REAL RX(N,NA),RY(N,NA),RZ(N,NA) ATOM POSITIONS AT TIME T ** C ** REAL OX(N,NA),OY(N,NA),OZ(N,NA) OLD POSITIONS AT TIME T-DT ** C ** REAL FX(N,NA),FY(N,NA),FZ(N,NA) ATOM FORCES AT TIME T ** C ** REAL DSQ(NA) SQUARED BOND LENGTHS ** C ** REAL M(NA) ATOMIC MASSES ** C ** ** C ** USAGE: ** C ** ** C ** THE ROUTINE IS CALLED AFTER COMPUTATION OF THE FORCES. ** C ** THE VERLET ALGORITHM IS USED TO ADVANCE THE POSITIONS FROM ** C ** RX,RY,RZ TO PX,PY,PZ, WITHOUT ANY CONSTRAINTS APPLIED. ** C ** THE ROUTINE THEN USES THE DESIRED SQUARED BOND LENGTHS IN DSQ ** C ** TO CONSTRAIN THE NEW POSITIONS. ** C ** DSQ(A) IS THE SQUARED BOND LENGTH BETWEEN ATOM A AND A+1. ** C ** NB IS THE NUMBER OF SUCH BONDS. IF NB = NA - 1 THE MOLECULE ** C ** IS NON-CYCLIC WHILE IF NB = NA THE MOLECULE IS CYCLIC. ** C ** THE ROUTINE ALSO REQUIRES THE DESIRED TOLERANCE AND AN UPPER ** C ** LIMIT TO THE NUMBER OF ITERATIONS IN CASE OF NON-CONVERGENCE. ** C ** THE ROUTINE USES TWO LOGICAL ARRAYS TO KEEP TRACK OF WHETHER ** C ** OR NOT WE HAVE MOVED (I.E. CORRECTED) THE ATOM POSITIONS: ** C ** MOVING(A) A=1,NA SAYS WHETHER WE ARE MOVING ATOM A THIS TIME ** C ** MOVED(A) A=1,NA SAYS WHETHER WE MOVED ATOM A LAST TIME. ** C ** THIS IS SO THAT WE CAN STOP CORRECTING THE POSITIONS OF ATOMS ** C ** WHENEVER POSSIBLE, SO AS TO CUT DOWN ON UNNECESSARY WORK. ** C ** THE ROUTINE RETURNS WITH NEW (T+DT) CONSTRAINED POSITIONS IN ** C ** RX,RY,RZ AND OLD (T) POSITIONS IN OX,OY,OZ. ** C ** THE ROUTINE ALSO CALCULATES THE KINETIC ENERGY AT TIME T ** C ** AND THE CONSTRAINT CONTRIBUTION TO THE VIRIAL WC. ** C ******************************************************************* INTEGER N PARAMETER ( N = 108 ) INTEGER NA PARAMETER ( NA = 6 ) INTEGER MAXIT, NB REAL RX(N,NA), RY(N,NA), RZ(N,NA) REAL OX(N,NA), OY(N,NA), OZ(N,NA) REAL FX(N,NA), FY(N,NA), FZ(N,NA) REAL DSQ(NA) REAL M(NA) REAL TOL, DT, K, WC, BOX REAL RXI(NA), RYI(NA), RZI(NA) REAL PXI(NA), PYI(NA), PZI(NA) LOGICAL MOVING(NA) LOGICAL MOVED(NA) LOGICAL DONE INTEGER IT, A, B, I REAL PXAB, PYAB, PZAB, PABSQ REAL RXAB, RYAB, RZAB, RABSQ, DIFFSQ, RPAB REAL GAB, DX, DY, DZ, TOL2, DTSQ, VXIA, VYIA, VZIA REAL BOXINV, RPTOL, RMA, RMB PARAMETER ( RPTOL = 1.0E-6 ) C ******************************************************************* IF ( ( NB .NE. NA ) .AND. ( NB .NE. NA-1 ) ) STOP 'NB IN ERROR' BOXINV = 1.0 / BOX DTSQ = DT ** 2 TOL2 = 2.0 * TOL K = 0.0 WC = 0.0 C ** LOOP OVER ALL MOLECULES ** DO 2000 I = 1, N C ** VERLET ALGORITHM ** DO 100 A = 1, NA RXI(A) = RX(I,A) RYI(A) = RY(I,A) RZI(A) = RZ(I,A) PXI(A) = 2.0 * RX(I,A) - OX(I,A) + DTSQ * FX(I,A) / M(A) PYI(A) = 2.0 * RY(I,A) - OY(I,A) + DTSQ * FY(I,A) / M(A) PZI(A) = 2.0 * RZ(I,A) - OZ(I,A) + DTSQ * FZ(I,A) / M(A) MOVING(A) = .FALSE. MOVED(A) = .TRUE. 100 CONTINUE IT = 0 DONE = .FALSE. C ** BEGIN ITERATIVE LOOP ** 1000 IF ( ( .NOT. DONE ) .AND. ( IT .LE. MAXIT ) ) THEN DONE = .TRUE. DO 300 A = 1, NB B = A + 1 IF ( B .GT. NA ) B = 1 IF ( MOVED(A) .OR. MOVED(B) ) THEN PXAB = PXI(A) - PXI(B) PXAB = PXAB - ANINT ( PXAB * BOXINV ) * BOX PYAB = PYI(A) - PYI(B) PYAB = PYAB - ANINT ( PYAB * BOXINV ) * BOX PZAB = PZI(A) - PZI(B) PZAB = PZAB - ANINT ( PZAB * BOXINV ) * BOX PABSQ = PXAB ** 2 + PYAB ** 2 + PZAB ** 2 RABSQ = DSQ(A) DIFFSQ = RABSQ - PABSQ IF ( ABS(DIFFSQ) .GT. ( RABSQ * TOL2 ) ) THEN RXAB = RXI(A) - RXI(B) RXAB = RXAB - ANINT ( RXAB * BOXINV ) * BOX RYAB = RYI(A) - RYI(B) RYAB = RYAB - ANINT ( RYAB * BOXINV ) * BOX RZAB = RZI(A) - RZI(B) RZAB = RZAB - ANINT ( RZAB * BOXINV ) * BOX RPAB = RXAB * PXAB + RYAB * PYAB + RZAB * PZAB IF ( RPAB .LT. ( RABSQ * RPTOL ) ) THEN STOP 'CONSTRAINT FAILURE ' ENDIF RMA = 1.0 / M(A) RMB = 1.0 / M(B) GAB = DIFFSQ / ( 2.0 * ( RMA + RMB ) * RPAB ) WC = WC + GAB * RABSQ DX = RXAB * GAB DY = RYAB * GAB DZ = RZAB * GAB PXI(A) = PXI(A) + RMA * DX PYI(A) = PYI(A) + RMA * DY PZI(A) = PZI(A) + RMA * DZ PXI(B) = PXI(B) - RMB * DX PYI(B) = PYI(B) - RMB * DY PZI(B) = PZI(B) - RMB * DZ MOVING(A) = .TRUE. MOVING(B) = .TRUE. DONE = .FALSE. ENDIF ENDIF 300 CONTINUE DO 400 A = 1, NA MOVED(A) = MOVING(A) MOVING(A) = .FALSE. 400 CONTINUE IT = IT + 1 GOTO 1000 ENDIF C ** END ITERATIVE LOOP ** IF ( .NOT. DONE ) THEN WRITE(*,'('' TOO MANY CONSTRAINT ITERATIONS '')') WRITE(*,'('' MOLECULE '',I5)') I STOP ENDIF DO 500 A = 1, NA VXIA = 0.5 * ( PXI(A) - OX(I,A) ) / DT VYIA = 0.5 * ( PYI(A) - OY(I,A) ) / DT VZIA = 0.5 * ( PZI(A) - OZ(I,A) ) / DT K = K + ( VXIA ** 2 + VYIA ** 2 + VZIA ** 2 ) * M(A) RX(I,A) = PXI(A) RY(I,A) = PYI(A) RZ(I,A) = PZI(A) OX(I,A) = RXI(A) OY(I,A) = RYI(A) OZ(I,A) = RZI(A) 500 CONTINUE 2000 CONTINUE C ** END LOOP OVER MOLECULES ** WC = WC / DTSQ / 3.0 K = 0.5 * K RETURN END