C     Cedar Fortran KAP   (Fri, 16 Nov 1990)  o3,r3,d3 11 Dec 1990 11:22:58
C***********************************************************************
      SUBROUTINE BNDRY(BOXL,V,NV,NS)
C***********************************************************************
C
C.....THIS ROUTINE PUTS THE MOLECULES BACK INSIDE THE BOX IF OUT
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION V(1)
C
D        call trace_enter (2001) 
      DO 100 I=2,NV,NS
      IF (V(I) .GT. BOXL) THEN
      V(I-1) = V(I-1) - BOXL
      V(I) = V(I) - BOXL
      V(I+1) = V(I+1) - BOXL
      ELSE
   10 IF (V(I) .LT. 0.0D0) THEN
      V(I-1) = V(I-1) + BOXL
      V(I) = V(I) + BOXL
      V(I+1) = V(I+1) + BOXL
      END IF
      END IF
  100 CONTINUE
D        call trace_exit (2002) 
      RETURN
      END
C     Cedar Fortran KAP   (Fri, 16 Nov 1990)  o3,r3,d3 11 Dec 1990 11:22:58
C***********************************************************************
      SUBROUTINE CNSTNT(N,C)
C***********************************************************************
C
C.....THIS ROUTINE DEFINES CONSTANTS TO BE USED
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /WATER/ OMAS,HMAS,WTMOL,ROH,ANGLE,FHM,FOM,ROHI,ROHI2,NATOMS
      COMMON /WWPOT/ QQ,A1,B1,A2,B2,A3,B3,A4,B4,AB1,AB2,AB3,AB4,C1,C2,
     *               QQ2,QQ4,REF1,REF2,REF4
      COMMON /CNST / UNITT,UNITL,UNITM,BOLTZ,AVGNO,PCC(10)
      COMMON /FRCNST/ FC11,FC12,FC13,FC33, FC111,FC333,FC112,FC113,
     *                FC123,FC133, FC1111,FC3333,FC1112,
     *                FC1122,FC1113,FC1123,FC1133,FC1233,FC1333
      DIMENSION C(N,N)
C
C.....MOLECULAR CONSTANTS FOR WATER IN ANGSTROM, RADIAN, AND A.M.U.
C
      NATOMS=3
      ROH=0.9572
      ROHI=1.0D0/ROH
      ROHI2=ROHI*ROHI
      ANGLE=1.824218
      OMAS=15.99945
      HMAS=1.007825
      WTMOL=OMAS+2.0D0*HMAS
C
C.....UNITS USED TO SCALE VARIABLES (IN C.G.S.)
C
      UNITT=1.0E-15
      UNITL=1.0E-8
      UNITM=1.6605655E-24
      BOLTZ=1.380662E-16
      AVGNO=6.022045E23
C
C.....FORCE CONSTANTS SCALED(DIVIDED) BY (UNITM/UNITT**2)
C
      FC11= 0.512596
      FC33= 0.048098
      FC12=-0.005823
      FC13= 0.016452
      FC111=-0.57191
      FC333=-0.007636
      FC112=-0.001867
      FC113=-0.002047
      FC123=-0.03083
      FC133=-0.0094245
      FC1111= 0.8431
      FC3333=-0.00193
      FC1112=-0.0030
      FC1122= 0.0036
      FC1113=-0.012
      FC1123= 0.0060
      FC1133=-0.0048
      FC1233= 0.0211
      FC1333= 0.006263
C
C.....WATER-WATER INTERACTION PARAMETERS
C
      QQ=0.07152158
      A1=455.313100
      B1=5.15271070
      A2=0.27879839
      B2=2.76084370
      A3=0.60895706
      B3=2.96189550
      A4=0.11447336
      B4=2.23326410
      CM=0.45682590
      AB1=A1*B1
      AB2=A2*B2
      AB3=A3*B3
      AB4=A4*B4
      C1=1.0D0-CM
      C2=0.5D0*CM
      QQ2=2.0D0*QQ
      QQ4=2.0D0*QQ2
C
C.....CALCULATE THE COEFFICIENTS OF TAYLOR SERIES EXPANSION
C     FOR F(X), F'(X), F''(X), ...... (WITH DELTAT**N/N] INCLUDED)
C     IN C(1,1),..... C(1,2),..... C(1,3),....... VECTORS
C
      C(1,1)=1.0D0
D        call trace_enter (2003) 
      DO 1000 N1=2,N
      NN=N1-1
      TN=NN
      C(1,N1) = 1.0D0
      TK=1.0D0
      DO 1100 K1=2,N1
      C(K1,NN) = C(K1-1,NN+1) * TN / TK
      NN = NN - 1
      TN = TN - 1.0D0
      TK = TK + 1.0D0
 1100 CONTINUE
 1000 CONTINUE
C
C.....PREDICTOR-CORRECTOR CONSTANTS FOR 2ND ORDER DIFFERENTIAL EQUATION
C
D        call trace_exit (2004) 
      PCC(3)=1.0D0
      N1=N-1
      GO TO (1,2,3,4,5,6,7), N1
    1 WRITE (6,10)
      RETURN
    2 WRITE (6,10)
      RETURN
    3 PCC(1)=1.0D0/6.0D0
      PCC(2)=5.0D0/6.0D0
      PCC(4)=1.0D0/3.0D0
      RETURN
    4 PCC(1)=19.0D0/120.0D0
      PCC(2)=3.0D0/4.0D0
      PCC(4)=1.0D0/2.0D0
      PCC(5)=1.0D0/12.0D0
      RETURN
    5 PCC(1)=3.0D0/20.0D0
      PCC(2)=251.0D0/360.0D0
      PCC(4)=11.0D0/18.0D0
      PCC(5)=1.0D0/6.0D0
      PCC(6)=1.0D0/60.0D0
      RETURN
    6 PCC(1)=863.0D0/6048.0D0
      PCC(2)=665.0D0/1008.0D0
      PCC(4)=25.0D0/36.0D0
      PCC(5)=35.0D0/144.0D0
      PCC(6)=1.0D0/24.0D0
      PCC(7)=1.0D0/360.0D0
      RETURN
    7 PCC(1)=275.0D0/2016.0D0
      PCC(2)=19087.0D0/30240.0D0
      PCC(4)=137.0D0/180.0D0
      PCC(5)=5.0D0/16.0D0
      PCC(6)=17.0D0/240.0D0
      PCC(7)=1.0D0/120.0D0
      PCC(8)=1.0D0/2520.0D0
      RETURN
   10 FORMAT('0***** ERROR: THE ORDER HAS TO BE GREATER THAN 2 ****')
      END
C     Cedar Fortran KAP   (Fri, 16 Nov 1990)  o3,r3,d3 11 Dec 1990 11:22:58
C***********************************************************************
      SUBROUTINE CORREC(VAR,SD,FR,NT,PCC,NOR1)
C***********************************************************************
C
C.....THIS ROUTINE CALCULATES CORRECTED F(X), F'(X), F"(X), ....
C     FROM CORRECTED F(X) = PREDICTED F(X) + PCC(1)*(FR-SD)
C     WHERE SD IS PREDICTED F"(X) AND FR IS FORCE AT PREDICTED POSITION
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION VAR(1),SD(1),FR(1),PCC(1)
      INTEGER II1, J1, J2, J3, J4, J5
C
D        call trace_enter (2005) 
        CDOALL 1000 I=1,NT

        double precision y
        Y=FR(I)-SD(I)
C      JJ=I
C      II1 = I
C      J2 = (NOR1 + 31) / 32
C      CDOALL 1100 J1=0,32*J2-32,32
C      INTEGER J6, J7
C      J7 = J1 + 1
C      J6 = MIN0 (J1 + 1 + 32 - 1, NOR1) - J1
C      VAR(II1+((-1)+J7)*NT:$J6:NT) = VAR(II1+((-1)+J7)*NT:$J6:NT) + PCC(
C     X  J7:$J6) * Y
C 1100 CONTINUE
       VAR(I:$NOR1:NT) = VAR(I:$NOR1:NT) + PCC(1:$NOR1) * Y
 1000 CONTINUE

D        call trace_exit (2006) 
      RETURN
      END
C     Cedar Fortran KAP   (Fri, 16 Nov 1990)  o3,r3,d3 11 Dec 1990 11:22:58
C***********************************************************************
      SUBROUTINE INITIA(NFMC,XC,YC,ZC,VX,VY,VZ,MD)
C***********************************************************************
C
C.....THIS ROUTINE INITIALIZES POSITIONS IN A CUBE AND RANDOMIZES
C     VELOCITY OF EACH ATOM
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /MDVAR/ TEMP,RHO,TSTEP,BOXL,BOXH,CUTOFF,CUT2,NMOL,NORDER,
     *               NATMO,NATMO3,NMOL1
      COMMON /WATER/ OMAS,HMAS,WTMOL,ROH,ANGLE,FHM,FOM,ROHI,ROHI2,NATOMS
      COMMON /CNST / UNITT,UNITL,UNITM,BOLTZ,AVGNO,PCC(10)
      DIMENSION XC(MD),YC(MD),ZC(MD),VX(2),VY(2),VZ(2)
      DIMENSION XT(3),YT(3),XMAS(3)
      SAVE XMIN, YMIN, ZMIN
      INTEGER I1, I2, I3, I4, I5
      DATA XMIN,YMIN,ZMIN/3*0.0D0/
C
      RANDOM(L)=2.0D0*RAND(0)-1.0D0
      XMAS(2)=SQRT(OMAS*HMAS)
      XMAS(1)=HMAS
      XMAS(3)=HMAS
C
C.....ASSIGN POSITIONS
C
      IF (NFMC .LE. 0) THEN
      WRITE (6,1)
    1 FORMAT('0***** NEW RUN STARTING FROM REGULAR LATTICE *****')
      NS=REAL(NMOL)**(1.0D0/3.0D0)+0.00001D0
      XS=BOXL/NS
      ZERO=XS*0.5D0
      WCOS=ROH*COS(ANGLE*0.5D0)
      WSIN=ROH*SIN(ANGLE*0.5D0)
      IJK=0
      XT(2)=ZERO
D        call trace_enter (2007) 
      DO 1000 I=1,NS
      XT(1) = XT(2) + WCOS
      XT(3) = XT(1)
      YT(2) = ZERO
      DO 1100 J=1,NS
      YT(1) = YT(2) + WSIN
      YT(3) = YT(2) - WSIN
      Z=ZERO
      DO 1110 K=1,NS
      CDOALL 1111 L=1,NATOMS
      XC(L+MAX0 (NATOMS, 0)*((-1)+K)+MAX0 (NATOMS, 0)*MAX0 (NS, 0)*((-1)
     X  +J)+MAX0 (NATOMS, 0)*MAX0 (NS, 0)*MAX0 (NS, 0)*((-1)+I)) = XT(L)
      YC(L+MAX0 (NATOMS, 0)*((-1)+K)+MAX0 (NATOMS, 0)*MAX0 (NS, 0)*((-1)
     X  +J)+MAX0 (NATOMS, 0)*MAX0 (NS, 0)*MAX0 (NS, 0)*((-1)+I)) = YT(L)
      ZC(L+MAX0 (NATOMS, 0)*((-1)+K)+MAX0 (NATOMS, 0)*MAX0 (NS, 0)*((-1)
     X  +J)+MAX0 (NATOMS, 0)*MAX0 (NS, 0)*MAX0 (NS, 0)*((-1)+I)) = Z
 1111 CONTINUE
      IF (MAX0 (NATOMS, 0) .GT. 0) IJK = MAX0 (NATOMS, 0) + MAX0 (NATOMS
     X  , 0) * ((-1) + K) + MAX0 (NATOMS, 0) * MAX0 (NS, 0) * ((-1) + J)
     X   + MAX0 (NATOMS, 0) * MAX0 (NS, 0) * MAX0 (NS, 0) * ((-1) + I)
      Z = Z + XS
 1110 CONTINUE
      YT(2) = YT(2) + XS
 1100 CONTINUE
      XT(2) = XT(2) + XS
 1000 CONTINUE
D        call trace_exit (2008) 
        ELSE
   10 REWIND NFMC
      WRITE (6,2) NFMC
    2 FORMAT('0***** NEW RUN STARTING FROM COORDINATES OF WATERS',
     *       ' IN FT',I2,' *****')
      READ (NFMC,11) XC
      READ (NFMC,11) YC
      READ (NFMC,11) ZC
   11 FORMAT(5E16.8)
      REWIND NFMC
      XMIN = MIN (XMIN, MINVAL (XC(1:$NATMO)))
      YMIN = MIN (YMIN, MINVAL (YC(1:$NATMO)))
      ZMIN = MIN (ZMIN, MINVAL (ZC(1:$NATMO)))
      I2 = (NATMO + 31) / 32
D        call trace_enter (2009) 
      CDOALL 400 I1=0,32*I2-32,32
      INTEGER I6, I7
      I7 = I1 + 1
      I6 = MIN0 (I1 + 1 + 32 - 1, NATMO) - I1
      XC(I7:$I6) = XC(I7:$I6) - XMIN
      YC(I7:$I6) = YC(I7:$I6) - YMIN
      ZC(I7:$I6) = ZC(I7:$I6) - ZMIN
  400 CONTINUE
C
C.....ASSIGN RANDOM MOMENTA
C
D        call trace_exit (2010) 
        END IF
   12 SUX=RAND(1)
      SUMX=0.0D0
      SUMY=0.0D0
      SUMZ=0.0D0
      IJ=1
D        call trace_enter (2011) 
      DO 2000 I=1,NMOL
      DO 2000 J=1,NATOMS
      VX(IJ)=RANDOM(0)*XMAS(J)
      VY(IJ)=RANDOM(0)*XMAS(J)
      VZ(IJ)=RANDOM(0)*XMAS(J)
      SUMX=SUMX+VX(IJ)
      SUMY=SUMY+VY(IJ)
      SUMZ=SUMZ+VZ(IJ)
      IJ=IJ+1
 2000 CONTINUE
C
C.....FIND AVERAGE MOMENTA PER ATOM
C
D        call trace_exit (2012) 
      SUMX=SUMX/(NATOMS*NMOL)
      SUMY=SUMY/(NATOMS*NMOL)
      SUMZ=SUMZ/(NATOMS*NMOL)
C
C.....FIND NORMALIZATION FACTOR SO THAT <K.E.>=KT/2
C
      SUX=0.0D0
      SUY=0.0D0
      SUZ=0.0D0
      J=1
      SUX = SUX + SUM (((VX(1:$NMOL:NATOMS) - SUMX) ** 2 + (VX(3:$NMOL:
     X       NATOMS) - SUMX) ** 2) / HMAS + (VX(2:$NMOL:NATOMS) - SUMX)
     X       ** 2 / OMAS)
      SUY = SUY + SUM (((VY(1:$NMOL:NATOMS) - SUMY) ** 2 + (VY(3:$NMOL:
     X       NATOMS) - SUMY) ** 2) / HMAS + (VY(2:$NMOL:NATOMS) - SUMY)
     X       ** 2 / OMAS)
      SUZ = SUZ + SUM (((VZ(1:$NMOL:NATOMS) - SUMZ) ** 2 + (VZ(3:$NMOL:
     X       NATOMS) - SUMZ) ** 2) / HMAS + (VZ(2:$NMOL:NATOMS) - SUMZ)
     X       ** 2 / OMAS)
      FAC=BOLTZ*TEMP*NATMO/UNITM * (UNITT*TSTEP/UNITL)**2
      SUX=SQRT(FAC/SUX)
      SUY=SQRT(FAC/SUY)
      SUZ=SQRT(FAC/SUZ)
C
C.....NORMALIZE INDIVIDUAL VELOCITIES SO THAT THERE IS NO BULK MOMENTA
C
      XMAS(2)=OMAS
      J=1
D        call trace_enter (2013) 
      CDOALL 200 I=1,NMOL
      VX(NATOMS*((-1)+I)+1:$NATOMS) = (VX(NATOMS*((-1)+I)+1:$NATOMS) -
     X  SUMX) * SUX / XMAS(1:$NATOMS)
      VY(NATOMS*((-1)+I)+1:$NATOMS) = (VY(NATOMS*((-1)+I)+1:$NATOMS) -
     X  SUMY) * SUY / XMAS(1:$NATOMS)
      VZ(NATOMS*((-1)+I)+1:$NATOMS) = (VZ(NATOMS*((-1)+I)+1:$NATOMS) -
     X  SUMZ) * SUZ / XMAS(1:$NATOMS)
  200 CONTINUE
D        call trace_exit (2014) 
      RETURN
      END
C     Cedar Fortran KAP   (9 Mar, 1990)       o3,r3,d3 26 Mar 1990 14:03:56
c   this is the manually transformed version of this subroutine.
c   
c   manual chage diary:
c   1) localization of arrays XL..GG
c   2) natoms=3 assignment
c   3) manual concurrentization for one cluster:
c     - induction variable substitution (JW1, JWO, JW2, IW1, IWO, IW2)
c     - temporary variables for forces. Updates are at the end of the do_1100
c       iteration under a qlock
c     - also temp. var for vir. Update is at the end of each do_1000 iteration
c   4) local declaration of all loop local scalars
c   5) improved update of accumulation variables: fxi/fyi/fzi updated in outer loop
c
C***********************************************************************
      SUBROUTINE INTERF(X,Y,Z,FX,FY,FZ,XM,YM,ZM,VIR)
C***********************************************************************
C
C.....THIS ROUTINE CALCULATES INTER-MOLECULAR INTERACTION FORCES
C     THE DISTANCES ARE ARRANGED IN THE ORDER OF M-M, M-H1, M-H3, H1-M,
C     H3-M, H1-H3, H1-H1, H3-H1, H3-H3, O-O, O-H1, O-H3, H1-O, H3-O
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
c      implicit none
      COMMON /MDVAR/TEMP, RHO, TSTEP, BOXL, BOXH, CUTOFF, CUT2,
     X  NMOL, NORDER, NATMO, NATMO3, NMOL1
      COMMON /WATER/OMAS, HMAS, WTMOL, ROH, ANGLE, FHM, FOM,
     X  ROHI, ROHI2, NATOMS
      COMMON /WWPOT/QQ, A1, B1, A2, B2, A3, B3, A4, B4, AB1, AB2
     X  , AB3, AB4, C1, C2, QQ2, QQ4, REF1, REF2, REF4
      COMMON /CNST/UNITT, UNITL, UNITM, BOLTZ, AVGNO, PCC(10)
      DIMENSION X(1), Y(1), Z(1), FX(1), FY(1), FZ(1), XM(1), YM(1), ZM(
     X  1)
c      DIMENSION XL(NMOL,NMOL1,14), YL(NMOL,NMOL1,14), ZL(NMOL,NMOL1,14)
c     x, RL(NMOL,NMOL1,14), RS(NMOL,NMOL1,14), FF(NMOL,NMOL1,14)
c     x, GG(NMOL,NMOL1,14)
      INTEGER I1
      integer lock
D     real tt, dtime, t1(2)
C
c following line is exp line
      natoms = 3
D     tt = dtime(t1)
      lock = 0
D        call trace_enter (2015) 
      CDOALL 1000 I=1,NMOL1
      integer iw1,iwo,iw2,jw1,jwo,jw2,kc
      INTEGER II1, II2, II3, j, k
      double precision ftemp, tt, tt1, g110, g23, g45

      DOUBLE PRECISION XL(14), YL(14), ZL(14), RL(14), RS(14), FF(14), GG(14)
      double precision virl, fxi1, fxio, fxi2, fxj1, fxjo, fxj2
      double precision fyi1, fyio, fyi2, fyj1, fyjo, fyj2
      double precision fzi1, fzio, fzi2, fzj1, fzjo, fzj2
       REAL   XMI , XIW1 ,XIW2 , XIW11, XJW1 ,XMJ ,XJW2 , XJW11
      virl = 0
      fxi1 = 0
      fxio = 0
      fxi2 = 0
      fyi1 = 0
      fyio = 0
      fyi2 = 0
      fzi1 = 0
      fzio = 0
      fzi2 = 0


      IW1=(I-1)*NATOMS +1
      IWO=IW1 + 1
      IW2=IW1 + 2
      DO 1100 J=I+1,NMOL
      JW1 = (J-1) * NATOMS + 1
      JWO = JW1 +1
      JW2 = JW1 +2
      XMI = XM(I)
      XIW1 = X(IW1)
      XIW2 =  X(IW2)
      XIW11 = X(IWO)
      XJW1  =  X(JW1)
      XMJ   = XM(J)
      XJW2  =  X(JW2)
      XJW11 = X(JWO)
      XL(1) = XMI - XMJ 
      XL(2) = XMI - XJW1
      XL(3) = XMI - XJW2
      XL(4) =XIW1 - XMJ 
      XL(5) = XIW2 - XMJ 
      XL(6) =XIW1 - XJW1
      XL(7) =XIW1 - XJW2
      XL(8) = XIW2 - XJW1
      XL(9) = XIW2 - XJW2
      XL(10) = XIW11 - XJW11
      XL(11) = XIW11 - XJW1
      XL(12) = XIW11 - XJW2
      XL(13) =XIW1 - XJW11
      XL(14) = XIW2 - XJW11
C      XL(1) = XM(I) - XM(J)
C      XL(2) = XM(I) - X(JW1)
C      XL(3) = XM(I) - X(JW1+2)
C      XL(4) = X(IW1) - XM(J)
C      XL(5) = X(IW1+2) - XM(J)
C      XL(6) = X(IW1) - X(JW1)
C      XL(7) = X(IW1) - X(JW1+2)
C      XL(8) = X(IW1+2) - X(JW1)
C      XL(9) = X(IW1+2) - X(JW1+2)
C      XL(10) = X(IW1+1) - X(JW1+1)
C      XL(11) = X(IW1+1) - X(JW1)
C      XL(12) = X(IW1+1) - X(JW1+2)
C      XL(13) = X(IW1) - X(JW1+1)
C      XL(14) = X(IW1+2) - X(JW1+1)
      WHERE (ABS (XL(1:$14)) .GT. BOXH) XL(1:$14) = XL(1:$14
     X  ) - SIGN (BOXL, XL(1:$14))
      XMI = YM(I)
      XIW1 = Y(IW1)
      XIW2 =  Y(IW1+2)
      XIW11 = Y(IW1+1)
      XJW1  =  Y(JW1)
      XMJ   = YM(J)
      XJW2  =  Y(JW1+2)
      XJW11 = Y(JW1+1)
      YL(1) = XMI - XMJ 
      YL(2) = XMI - XJW1
      YL(3) = XMI - XJW2
      YL(4) =XIW1 - XMJ 
      YL(5) = XIW2 - XMJ 
      YL(6) =XIW1 - XJW1
      YL(7) =XIW1 - XJW2
      YL(8) = XIW2 - XJW1
      YL(9) = XIW2 - XJW2
      YL(10) = XIW11 - XJW11
      YL(11) = XIW11 - XJW1
      YL(12) = XIW11 - XJW2
      YL(13) =XIW1 - XJW11
      YL(14) = XIW2 - XJW11

C      YL(1) = YM(I) - YM(J)
C      YL(2) = YM(I) - Y(JW1)
C      YL(3) = YM(I) - Y(JW1+2)
C      YL(4) = Y(IW1) - YM(J)
C      YL(5) = Y(IW1+2) - YM(J)
C      YL(6) = Y(IW1) - Y(JW1)
C      YL(7) = Y(IW1) - Y(JW1+2)
C      YL(8) = Y(IW1+2) - Y(JW1)
C      YL(9) = Y(IW1+2) - Y(JW1+2)
C      YL(10) = Y(IW1+1) - Y(JW1+1)
C      YL(11) = Y(IW1+1) - Y(JW1)
C      YL(12) = Y(IW1+1) - Y(JW1+2)
C      YL(13) = Y(IW1) - Y(JW1+1)
C      YL(14) = Y(IW1+2) - Y(JW1+1)
      WHERE (ABS (YL(1:$14)) .GT. BOXH) YL(1:$14) = YL(1:$14
     X  ) - SIGN (BOXL, YL(1:$14))
 


      XMI = ZM(I)
      XIW1 = Z(IW1)
      XIW2 =  Z(IW1+2)
      XIW11 = Z(IW1+1)
      XJW1  =  Z(JW1)
      XMJ   = ZM(J)
      XJW2  =  Z(JW1+2)
      XJW11 = Z(JW1+1)
      ZL(1) = XMI - XMJ 
      ZL(2) = XMI - XJW1
      ZL(3) = XMI - XJW2
      ZL(4) =XIW1 - XMJ 
      ZL(5) = XIW2 - XMJ 
      ZL(6) =XIW1 - XJW1
      ZL(7) =XIW1 - XJW2
      ZL(8) = XIW2 - XJW1
      ZL(9) = XIW2 - XJW2
      ZL(10) = XIW11 - XJW11
      ZL(11) = XIW11 - XJW1
      ZL(12) = XIW11 - XJW2
      ZL(13) =XIW1 - XJW11
      ZL(14) = XIW2 - XJW11
C     ZL(1) = ZM(I) - ZM(J)
C      ZL(2) = ZM(I) - Z(JW1)
C      ZL(3) = ZM(I) - Z(JW1+2)
C      ZL(4) = Z(IW1) - ZM(J)
C      ZL(5) = Z(IW1+2) - ZM(J)
C      ZL(6) = Z(IW1) - Z(JW1)
C      ZL(7) = Z(IW1) - Z(JW1+2)
C      ZL(8) = Z(IW1+2) - Z(JW1)
C      ZL(9) = Z(IW1+2) - Z(JW1+2)
C      ZL(10) = Z(IW1+1) - Z(JW1+1)
C      ZL(11) = Z(IW1+1) - Z(JW1)
C      ZL(12) = Z(IW1+1) - Z(JW1+2)
C      ZL(13) = Z(IW1) - Z(JW1+1)
C      ZL(14) = Z(IW1+2) - Z(JW1+1)
      WHERE (ABS (ZL(1:$14)) .GT. BOXH) ZL(1:$14) = ZL(1:$14
     X  ) - SIGN (BOXL, ZL(1:$14))
      KC = 0
      RS(1:$9) = XL(1:$9) * XL(1:$9) + YL(1:$9) * YL(1:$9) + ZL(1:$9) * ZL(1:$9)
 1110 KC = KC + COUNT (MASK=RS(1:$9) .GT. CUT2)
      IF (KC .NE. 9) THEN
      FF(1:$14) = 0.0D0
      IF (RS(1) .LT. CUT2) THEN
      FF(1)=QQ4/(RS(1)*SQRT(RS(1)))+REF4
      VIRL=VIRL+FF(1)*RS(1)
      END IF
   10 DO 1130 K=2,5
      IF (RS(K) .LT. CUT2) THEN
      FF(K)=-QQ2/(RS(K)*SQRT(RS(K)))-REF2
      VIRL=VIRL+FF(K)*RS(K)
      END IF
   11 IF (RS(K+4) .LE. CUT2) THEN
      RL(K+4)=SQRT(RS(K+4))
      FF(K+4)=QQ/(RS(K+4)*RL(K+4))+REF1
      VIRL=VIRL+FF(K+4)*RS(K+4)
      END IF
 1130 CONTINUE
      IF (KC .EQ. 0) THEN
      RS(10)=XL(10)*XL(10)+YL(10)*YL(10)+ZL(10)*
     x ZL(10)
      RL(10)=SQRT(RS(10))
      FF(10)=AB1*EXP(-B1*RL(10))/RL(10)
      VIRL=VIRL+FF(10)*RS(10)
      DO 1140 K=11,14
      FTEMP=AB2*EXP(-B2*RL(K-5))/RL(K-5)
      FF(K-5)=FF(K-5)+FTEMP
      VIRL=VIRL+FTEMP*RS(K-5)
      RS(K)=XL(K)*XL(K)+YL(K)*YL(K)+ZL(K)*
     x ZL(K)
      RL(K)=SQRT(RS(K))
      FF(K)=(AB3*EXP(-B3*RL(K))-AB4*EXP(-B4*RL(K)))/
     x RL(K)
      VIRL=VIRL+FF(K)*RS(K)
 1140 CONTINUE
      END IF
C
C.....CALCULATE X-COMPONENT FORCES
C
 1150 GG(1:$14) = FF(1:$14) * XL(1:$14)
      G110=GG(10)+GG(1)*C1
      G23=GG(2)+GG(3)
      G45=GG(4)+GG(5)
      FXIO=FXIO+G110+GG(11)+GG(12)+C1*G23
      FXJO=-G110-GG(13)-GG(14)-C1*G45
      TT1=GG(1)*C2
      TT=G23*C2+TT1
      FXI1=FXI1+GG(6)+GG(7)+GG(13)+TT+GG(4)
      FXI2=FXI2+GG(8)+GG(9)+GG(14)+TT+GG(5)
      TT=G45*C2+TT1
      FXJ1=-GG(6)-GG(8)-GG(11)-TT-GG(2)
      FXJ2=-GG(7)-GG(9)-GG(12)-TT-GG(3)
C
C.....CALCULATE Y-COMPONENT FORCES
C
 1160 GG(1:$14) = FF(1:$14) * YL(1:$14)
      G110=GG(10)+GG(1)*C1
      G23=GG(2)+GG(3)
      G45=GG(4)+GG(5)
      FYIO=FYIO+G110+GG(11)+GG(12)+C1*G23
      FYJO=-G110-GG(13)-GG(14)-C1*G45
      TT1=GG(1)*C2
      TT=G23*C2+TT1
      FYI1=FYI1+GG(6)+GG(7)+GG(13)+TT+GG(4)
      FYI2=FYI2+GG(8)+GG(9)+GG(14)+TT+GG(5)
      TT=G45*C2+TT1
      FYJ1=-GG(6)-GG(8)-GG(11)-TT-GG(2)
      FYJ2=-GG(7)-GG(9)-GG(12)-TT-GG(3)
C
C.....CALCULATE Z-COMPONENT FORCES
C
 1170 GG(1:$14) = FF(1:$14) * ZL(1:$14)
      G110=GG(10)+GG(1)*C1
      G23=GG(2)+GG(3)
      G45=GG(4)+GG(5)
      FZIO=FZIO+G110+GG(11)+GG(12)+C1*G23
      FZJO=-G110-GG(13)-GG(14)-C1*G45
      TT1=GG(1)*C2
      TT=G23*C2+TT1
      FZI1=FZI1+GG(6)+GG(7)+GG(13)+TT+GG(4)
      FZI2=FZI2+GG(8)+GG(9)+GG(14)+TT+GG(5)
      TT=G45*C2+TT1
      FZJ1=-GG(6)-GG(8)-GG(11)-TT-GG(2)
      FZJ2=-GG(7)-GG(9)-GG(12)-TT-GG(3)
      END IF
      call set_qlock(lock)
      fx(jw1) = fx(jw1) + fxj1
      fx(jwo) = fx(jwo) + fxjo
      fx(jw2) = fx(jw2) + fxj2
      fy(jw1) = fy(jw1) + fyj1
      fy(jwo) = fy(jwo) + fyjo
      fy(jw2) = fy(jw2) + fyj2
      fz(jw1) = fz(jw1) + fzj1
      fz(jwo) = fz(jwo) + fzjo
      fz(jw2) = fz(jw2) + fzj2
      call clear_qlock(lock)

 1100 CONTINUE
      call set_qlock(lock)
      fx(iw1) = fx(iw1) + fxi1
      fx(iwo) = fx(iwo) + fxio
      fx(iw2) = fx(iw2) + fxi2
      fy(iw1) = fy(iw1) + fyi1
      fy(iwo) = fy(iwo) + fyio
      fy(iw2) = fy(iw2) + fyi2
      fz(iw1) = fz(iw1) + fzi1
      fz(iwo) = fz(iwo) + fzio
      fz(iw2) = fz(iw2) + fzi2
      vir = vir + virl
      call clear_qlock(lock)
 1000 CONTINUE
D        call trace_exit (2016) 
D     print *,dtime(t1)
C
C.....DIVIDE FINAL FORCES BY MASSES
C
      I1 = (NATMO + 2) / 3
      FX(1:$I1:3) = FX(1:$I1:3) * FHM
      FY(1:$I1:3) = FY(1:$I1:3) * FHM
      FZ(1:$I1:3) = FZ(1:$I1:3) * FHM
      FX(3:$I1:3) = FX(3:$I1:3) * FHM
      FY(3:$I1:3) = FY(3:$I1:3) * FHM
      FZ(3:$I1:3) = FZ(3:$I1:3) * FHM
      FX(2:$I1:3) = FX(2:$I1:3) * FOM
      FY(2:$I1:3) = FY(2:$I1:3) * FOM
      FZ(2:$I1:3) = FZ(2:$I1:3) * FOM
      RETURN
      END
C  last modification: Fri Jan 10 17:11:18 1992
C     Cedar Fortran KAP   (Fri, 15 Feb 1991)  o3,r1,d3  8 Mar 1991 16:31:38
C***********************************************************************
      SUBROUTINE INTRAF(V,FR,VM,VIR)
C***********************************************************************
C
C.....THIS ROUTINE CALCULATES THE INTRA-MOLECULAR FORCE/MASS ACTING ON
C     EACH ATOM IN FR(I).
C     FC11, FC12, FC13, AND FC33 ARE THE QUARDRATIC FORCE CONSTANTS
C     FC111, FC112, ....... ETC. ARE THE CUBIC      FORCE CONSTANTS
C     FC1111, FC1112 ...... ETC. ARE THE QUARTIC    FORCE CONSTANTS
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
       COMMON /MDVAR/TEMP, RHO, TSTEP, BOXL, BOXH, CUTOFF, CUT2,
     X  NMOL, NORDER, NATMO, NATMO3, NMOL1
       COMMON /WATER/OMAS, HMAS, WTMOL, ROH, ANGLE, FHM, FOM,
     X  ROHI, ROHI2, NATOMS
       COMMON /WWPOT/QQ, A1, B1, A2, B2, A3, B3, A4, B4, AB1, AB2
     X  , AB3, AB4, C1, C2, QQ2, QQ4, REF1, REF2, REF4
       COMMON /FRCNST/FC11, FC12, FC13, FC33, FC111, FC333, FC112
     X  , FC113, FC123, FC133, FC1111, FC3333, FC1112, FC1122, FC1113,
     X  FC1123, FC1133, FC1233, FC1333
c      DIMENSION VR1(3),VR2(3),DR11(3),DR23(3),DT1(3),DT3(3)
       DIMENSION V(1),FR(1),VM(1)
       double precision      VIR, VM, FR, V

C
      ICZ=1
D        call trace_enter (2017)

      CDOALL 1000 I=1,NMOL
      integer ic, j, im
      double precision sum, r1, r2, sin, cos, dt, dts, dr1, dr1s, dr2, dr2s
     x, r1s, r2s, f1, f2, f3, t1, t2
      double precision VR1(3),VR2(3),DR11(3),DR23(3),DT1(3),DT3(3)

      IC = ((-1) + I) * NATOMS + 1
      SUM=0.0D0
      R1=0.0D0
      R2=0.0D0
      IM=I
      DO 1100 J=1,3
      VM(IM)=C1*V(IC+1)+C2*(V(IC)+V(IC+2))
      IM=IM+NMOL
      VR1(J)=V(IC+1)-V(IC)
      R1=R1+VR1(J)*VR1(J)
      VR2(J)=V(IC+1)-V(IC+2)
      R2=R2+VR2(J)*VR2(J)
      SUM=SUM+VR1(J)*VR2(J)
      IC=IC+NATMO
 1100 CONTINUE
      R1=SQRT(R1)
      R2=SQRT(R2)
C
C.....CALCULATE COS(THETA), SIN(THETA), DELTA(R1), DELTA(R2),
C.....AND DELTA(THETA)
C
      COS=SUM/(R1*R2)
      SIN=SQRT(1.0D0-COS*COS)
      DT=(ACOS(COS)-ANGLE)*ROH
      DTS=DT*DT
      DR1=R1-ROH
      DR1S=DR1*DR1
      DR2=R2-ROH
      DR2S=DR2*DR2
C
C.....CALCULATE DERIVATIVES OF R1/X1, R2/X3, THETA/X1, AND THETA/X3
C
      R1S=ROH/(R1*SIN)
      R2S=ROH/(R2*SIN)
      DO 1200 J=1,3
      DR11(J) = VR1(J) / R1
      DR23(J) = VR2(J) / R2
      DT1(J)=(-DR23(J)+DR11(J)*COS)*R1S
      DT3(J)=(-DR11(J)+DR23(J)*COS)*R2S
 1200 CONTINUE
C
C.....CALCULATE FORCES
C
      F1=FC11*DR1+FC12*DR2+FC13*DT
      F2=FC33*DT +FC13*(DR1+DR2)
      F3=FC11*DR2+FC12*DR1+FC13*DT
      F1=F1+(3.0D0*FC111*DR1S+FC112*(2.0D0*DR1+DR2)*DR2
     *      +2.0D0*FC113*DR1*DT+FC123*DR2*DT+FC133*DTS)*ROHI
      F2=F2+(3.0D0*FC333*DTS+FC113*(DR1S+DR2S)
     *      +FC123*DR1*DR2+2.0D0*FC133*(DR1+DR2)*DT)*ROHI
      F3=F3+(3.0D0*FC111*DR2S+FC112*(2.0D0*DR2+DR1)*DR1
     *      +2.0D0*FC113*DR2*DT+FC123*DR1*DT+FC133*DTS)*ROHI
      F1=F1+(4.0D0*FC1111*DR1S*DR1+FC1112*(3.0D0*DR1S+DR2S)
     *      *DR2+2.0D0*FC1122*DR1*DR2S+3.0D0*FC1113*DR1S*DT
     *      +FC1123*(2.0D0*DR1+DR2)*DR2*DT+(2.0D0*FC1133*DR1
     *      +FC1233*DR2+FC1333*DT)*DTS)*ROHI2
      F2=F2+(4.0D0*FC3333*DTS*DT+FC1113*(DR1S*DR1+DR2S*DR2)
     *      +FC1123*(DR1+DR2)*DR1*DR2+2.0D0*FC1133*(DR1S+DR2S)
     *      *DT+2.0D0*FC1233*DR1*DR2*DT+3.0D0*FC1333*(DR1+DR2)*DTS)
     *      *ROHI2
      F3=F3+(4.0D0*FC1111*DR2S*DR2+FC1112*(3.0D0*DR2S+DR1S)
     *      *DR1+2.0D0*FC1122*DR1S*DR2+3.0D0*FC1113*DR2S*DT
     *      +FC1123*(2.0D0*DR2+DR1)*DR1*DT+(2.0D0*FC1133*DR2
     *      +FC1233*DR1+FC1333*DT)*DTS)*ROHI2
      DO 1300 J=1,3
      T1=F1*DR11(J)+F2*DT1(J)
      FR(((-1)+I)*NATOMS+((-1)+J)*NATMO+1) = T1
      T2=F3*DR23(J)+F2*DT3(J)
      FR(((-1)+I)*NATOMS+((-1)+J)*NATMO+3) = T2
      FR(((-1)+I)*NATOMS+((-1)+J)*NATMO+2) = -(T1 + T2)
 1300 CONTINUE
 1000 CONTINUE
D        call trace_exit (2018) 
C        call end_interval( 19 )
      VIR=0.0D0
D        call trace_enter (2039) 
 2000 VIR = VIR + DOTPRODUCT (V(1:$NATMO3), FR(1:$NATMO3))
D        call trace_exit (2040) 
      RETURN
      END
C     Cedar Fortran KAP   (Fri, 16 Nov 1990)  o3,r3,d3 11 Dec 1990 11:22:58
C***********************************************************************
      SUBROUTINE KINETI (VAR, NMOL, SUM0, HMAS, OMAS)
C***********************************************************************
C
C.....THIS ROUTINE EVALUATES KINETIC ENERGY
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION VAR(1), SUM0(1)
C
      JJ=1
D        call trace_enter (2019) 
      CDOALL 100 K=1,3
      SUM0(K) = 0.0D0
      SUM0(K) = SUM0(K) + SUM ((VAR(3*(NMOL*((-1)+K))+1:$NMOL:3) * VAR(3
     X   *(NMOL*((-1)+K))+1:$NMOL:3) + VAR(3*(NMOL*((-1)+K))+3:$NMOL:3)
     X   * VAR(3*(NMOL*((-1)+K))+3:$NMOL:3)) * HMAS + VAR(3*(NMOL*((-1)+
     X   K))+2:$NMOL:3) * VAR(3*(NMOL*((-1)+K))+2:$NMOL:3) * OMAS)
  110 JJ = 3 * (NMOL + NMOL * ((-1) + K)) + 1
  100 CONTINUE
D        call trace_exit (2020) 
      RETURN
      END
c manually optimized subroutines: interf, poteng
c subroutine cshift not needed any more (it is inlined in interf and poteng)
c
C     Cedar Fortran KAP   (Fri, 16 Nov 1990)  o3,r3,d3 11 Dec 1990 11:22:58
      PROGRAM MDG
C
C.....DRIVER FOR MOLECULAR DYNAMIC SIMULATION OF FLEXIBLE WATER MOLECULE
C     WRITTEN BY GEORGE C. LIE, IBM, KINGSTON, N.Y.
C     APRIL 14, 1987 VERSION
C
C ******** INPUT TO THE PROGRAM ********
C
C   &MDINP
C
C     TEMP  : TEMPERATURE IN DEGREE K (DEFAULT=298.0).
C     RHO   : DENSITY IN G/C.C. (DEFAULT=0.998).
C     NORDER: ORDER USED TO SOLVE NEWTONIAN EQUATIONS (DEFAULT=5).
C     TSTEP : TIME STEP IN SECONDS (DEFAULT=1.0E-15).
C     NSTEP : NO. OF TIME STEPS FOR THIS RUN.
C     NPRINT: FREQUENCY OF PRINTING INTERMEDIATE DATA, SUCH AS KINETIC
C             ENERGY, POTENTIAL ENERGY, AVERAGE TEMPERATURE, ETC.
C             (ONE LINE PER PRINTING, DEFAULT=100).
C     LKT   : 1 IF RENORMALIZATION OF KINETIC ENERGY IS TO BE DONE JUST
C               BEFORE SAVING DATA FOR RESTART (NEEDED ONLY AT THE
C               BEGINNING WHERE THE ENERGY OR TEMPERATURE IS TOO HIGH);
C             0 OTHERWISE (DEFAULT=0).
C     NSAVE :-1 FOR THE VERY FIRST RUN;
C             0 DURING EQUILIBRATION STAGE;
C             N FREQUENCY OF DATA (X AND V) SAVING DURING DATA
C               COLLECTING STAGE (DEFAULT=10).
C     NRST  : FREQUENCY OF SAVING INTERMEDIATE DATA FOR RESTART.
C     CUTOFF: CUTOFF RADIUS FOR NEGLECTING FORCE AND POTENTIAL.
C             SET TO 0.0 IF HALF THE SIZE OF THE BOX IS TO BE USED
C             (DEFAULT=0.0D0)
C     NFMC  : 0 IF INITIALIZATION IS TO BE STARTED FROM REGULAR LATTICE;
C            11 IF INITIALIZATION IS TO BE STARTED FROM A RESTART FILE;
C             N IF INITIALIZATION IS TO BE STARTED FROM FT-N
C               WHICH CONTAINS THE COORDINATES OF THE WATER MOLECULES
C               IN THE FORMAT OF 5E16.8.  THE ORDERS ARE X OF H, O, H,
C               OF THE 1-ST WATER, X OF H, O, H OF THE 2-ND WATER, ....
C               FOLLOWED (STARTING FROM A NEW LINE) BY Y'S, THEN Z'S
C               (DEFAULT=12);
C            <0 TO RESET THE STATISTICAL COUNTERS.
C     NFSV  : FORTRAN FILE NO. FOR SAVED DATA (DEFAULT=10).
C     NFRST : FORTRAN FILE NO. FOR RESTART DATA (DEFAULT=11).
C
C   &END
C
C ******* FORTRAN FILES NEEDED ********
C
C     FT05F001  INPUT FILE CONTAINING NAMELIST &MDINP.
C     FT10F001  OUTPUT SEQUENTIAL FILE FOR SAVING X AND V DATA FOR
C               ANALYSIS; LOGICAL RECORD LENGTH IS 18*NMOL+10 DOUBLE
C               WORDS.
C     FT11F001  INPUT AND OUTPUT SEQUENTIAL RESTART FILE;
C               LOGICAL RECORD LENGTH IS 18*NMOL*3*(NORDER+1)+10
C               DOUBLE WORDS.
C     FT12F001  NEEDED ONLY IF NFMC=12 AND NSAVE=0;
C               CONTAINING COORDINATES OF ALL THE WATER MOLECULES.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 
      REAL WAL0,WAL1,CPU0,CPU1
 
      PARAMETER (NOMOL=343,MAXODR=7,NATOM=3,MXOD2=MAXODR+2,
     *           NDVAR=NATOM*NOMOL*3*MXOD2)
      COMMON /MDVAR/ TEMP,RHO,TSTEP,BOXL,BOXH,CUTOFF,CUT2,NMOL,NORDER,
     *               NATMO,NATMO3,NMOL1
      COMMON /WATER/ OMAS,HMAS,WTMOL,ROH,ANGLE,FHM,FOM,ROHI,ROHI2,NATOMS
      COMMON /WWPOT/ QQ,A1,B1,A2,B2,A3,B3,A4,B4,AB1,AB2,AB3,AB4,C1,C2,
     *               QQ2,QQ4,REF1,REF2,REF4
      COMMON /CNST / UNITT,UNITL,UNITM,BOLTZ,AVGNO,PCC(10)
      COMMON /MDDATA/ VAR(NDVAR),VM(3*NOMOL),TLC(100),
     *                ELPST,TKIN,TVIR,TTMV,FPOT,FKIN,IX(3*MXOD2),
     *                IRST,NVAR,NXYZ,NXV,IXF,IYF,IZF,IMY,IMZ
      INTEGER II1, I1, I2, I3, I4, I5, I6, II2
C
C.....DEFAULT VALUES FOR THE CONTROL PARAMETERS OF THE DRIVER
C
      DATA  NRST,NPRINT,NSAVE,NFRST,NFSV,NFMC, LKT
     *    / 1000,   100,   10,   11,  10,  12,   0/
C
C.....INPUT VARIABLES
C REMOVED NAMLIST 2/17/89 L. POINTER CSRD
C      NAMELIST /MDINP/ TEMP,RHO,NMOL,NORDER,TSTEP,NSTEP,NPRINT,
C     *                 NFRST,NFSV,NSAVE,CUTOFF,NRST,NFMC,LKT
C
C ******** SYSTEM DEPENDENT ROUTINE CALLS ********
C
C.....TIMER ROUTINE (THERE IS ANOTHER ONE TO CHANGE AT THE END
C     THERE ARE TWO WRITE STATEMENTS USING DATA FROM TIMER CALLS:
C     WRITE (6,1000) AND WRITE (6,1005)
C     CALL FCLOCK(DATE,TIME,ELPS,CPU)
C.....UNFORMATTED VARIABLE SPANNED RECORD IS LONGER THAN LRECL
C     CALL ERRSET(201,256,-1,1,1)
C.....FLOATING POINT UNDERFLOW SUPPRESSION
C     CALL XUFLOW
C     CALL TIMER (CPU,WALL,CPUTOT,WALTOT)
      CALL CPUTIM (CPU0)
      CALL ELAPSE (WAL0)
C
      OPEN(1,FILE='LWV',STATUS='UNKNOWN')
      OPEN(5,FILE='LWI5',STATUS='OLD')
      OPEN(6,FILE='LWO6',STATUS='UNKNOWN')
      OPEN(NFMC,FILE='LWI12',STATUS='OLD')
 
      WRITE(1,1)
1     FORMAT(1X,'OUTPUT FOR PERFECT CLUB BENCHMARK: MDG'/
     *       1X,'$Revision: 1.8 $ $Author: kipp $'/)
C
C.....READ INPUT
C
      NORDER=5
      TEMP  =298.0D0
      RHO   =0.998D0
      TSTEP =1.0D-15
      CUTOFF=0.0D0
C     READ (5,MDINP)
      READ (5,*) TSTEP, NMOL, NSTEP, NORDER, NSAVE, NRST, NPRINT, NFMC
C
C**** OPEN ALL FILES NEEDED HERE. THE LENGTH OF EACH FILE IS GIVEN ABOVE
C
C
C.....SET UP SCALING FACTORS AND CONSTANTS
C
      NORD1=NORDER+1
      CALL CNSTNT(NORD1,TLC)
      WRITE (6,1000) TIME,DATE
      WRITE (6,1001) TEMP,RHO,NMOL,TSTEP,NORDER,NSTEP,NSAVE,NRST
      CALL SYSCNS
      WRITE (6,1002) CUTOFF
C
C.....CALCULATE ARRAY INDEX FOR X,Y,Z,VX,VY,VZ,AX,AY,AZ,....,FX,FY,FZ
C
      II=1
      NVAR=3*(NORDER+2)
      NXYZ=(NVAR-3)*NATMO
      NXV=NATMO*6
      I3 = (NVAR + 31) / 32
D        call trace_enter (2021) 
      CDOALL 100 I1=0,32*I3-32,32
      INTEGER I7, I8
      I8 = I1 + 1
      I7 = MIN0 (I1 + 1 + 32 - 1, NVAR) - I1
      IX(I8:$I7) = [((-1) + I8) * NATMO + 1 :$ I7 : NATMO]
  100 END CDOALL
D        call trace_exit (2022) 
      IXF=IX(NVAR-2)
      IYF=IX(NVAR-1)
      IZF=IX(NVAR)
      IMY=NMOL+1
      IMZ=NMOL+IMY
D        call trace_enter (2023) 
      CDOALL 101 I2=0,27776,32
      INTEGER I9, I10
      I10 = I2 + 1
      I9 = MIN0 (I2 + 1 + 32 - 1, NDVAR) - I2
  101 VAR(I10:$I9) = 0.0D0
D        call trace_exit (2024) 
      I = MAX (1, NDVAR + 1)
C
C.....CHECK WHAT TO BE DONE (NEW RUN, CONTINUATION, ETC.)
C
      IF (NSAVE .LE. 0) THEN
C
C.....INITIALIZATION OR EQUILIBRIATION RUN
C
      IF (NSAVE .NE. 0) THEN
      IRST=0
C
C.....IF NEW RUN INITIALIZE POSITION, VELECITY, AND ACCELERATION
C
      IF(NFMC.EQ.NFRST) THEN
C
C.......INITIALIZATION FROM A RESTART FILE
C
        READ (NFRST)
        READ (NFRST) D,D,D,D,N,(VAR(I),I=1,N)
        REWIND NFRST
        WRITE (6,1003) NFRST
        CALL KINETI(VAR(IX(4)),NMOL,SUM,HMAS,OMAS)
        CALL NRMLKT(VAR(IX(4)),SUM,FKIN,NATMO)
C
C.......INITIALIZATION FROM A REGULAR LATTICE OR GEOMETRY
C
      ELSE
        CALL INITIA(NFMC,VAR,VAR(IX(2)),VAR(IX(3)),VAR(IX(4)),
     *              VAR(IX(5)),VAR(IX(6)),NATMO)
C
C.......ESTIMATE ACCELERATION FROM F/M
C
        CALL INTRAF(VAR,VAR(IXF),VM,VIR)
        CALL INTERF(VAR,VAR(IX(2)),VAR(IX(3)),VAR(IX(7)),VAR(IX(8)),
     *              VAR(IX(9)),VM,VM(IMY),VM(IMZ),VIR)
      ENDIF
      NFMC=-1
      GOTO 300
      END IF
      END IF
C
C.....IF RESTART, GET X,V,A,.... FROM RESTART FILE
C
  200 REWIND NFRST
      READ (NFRST) DTEMP,DRHO,NMOL,OLDTS,NORD,ELPST,TKIN,TTMV,IRST
      IF(DTEMP.NE.TEMP) WRITE (6,1006) TEMP,DTEMP
      IF(DRHO.NE.RHO) WRITE (6,1007) RHO,DRHO
      IF(NORD.NE.NORDER) WRITE (6,1008) NORDER,NORD
      IRST=IRST+1
      READ (NFRST) TVIR,D,D,D,N,(VAR(I),I=1,N)
      REWIND NFRST
      IF (OLDTS .NE. TSTEP) THEN
C
C.....IF TIME STEP CHANGED
C
      WRITE (6,1009) TSTEP,OLDTS
      VF=TSTEP/OLDTS
      TKIN=TKIN*VF
      TVIR=TVIR/VF
      TTMV=TTMV/VF
      FAC=1.0D0
      IJ=NATMO3+1
      II1 = NATMO3 + 1
      II2 = (-1)
D        call trace_enter (2025) 
      DO 210 I=7,NVAR,3
      II2 = II2 + 1
      FAC = FAC * VF
      CDOALL 211 J=1,NATMO3
      VAR(II1+(J+MAX0 (NATMO3, 0)*II2)-1) = VAR(II1+(J+MAX0 (NATMO3, 0)*
     X  II2)-1) * FAC
  211 END CDOALL
      IF (MAX0 (NATMO3, 0) .GT. 0) IJ = II1 + (MAX0 (NATMO3, 0) + MAX0 (
     X  NATMO3, 0) * II2)
  210 CONTINUE
C
C.....START MOLECULAR DYNAMIC LOOP
C
D        call trace_exit (2026) 
        END IF
  300 IF(NFMC.LT.0) THEN
      IF(NDATA.GT.0) WRITE (6,1010)
        ELPST=0.0D0
        TKIN=0.0D0
        TVIR=0.0D0
        TTMV=0.0D0
      ENDIF
      IF(NSAVE.GT.0) NRST=NRST/NSAVE*NSAVE
      CALL MDMAIN(NFSV,NFRST,NSTEP,NRST,NPRINT,NSAVE,LKT,NORD1)
C
C.....FINISHED, PRINT TOTAL TIME USED
C
      CALL CPUTIM (CPU1)
      CALL ELAPSE (WAL1)
C     CALL TIMER (CPU,WALL,CPUTOT,WALTOT)
      WRITE (6,1005) CPU1-CPU0,MOD(WAL1-WAL0+1000000.0,1000000.0)
C
      WRITE(1,*)
      WRITE(1,*) 'ELAPSED CPU TIME IN SECONDS: ',CPU1-CPU0
      WRITE(1,*) 'MFLOP RATE: ',3432.550/(CPU1-CPU0)
      WRITE(1,*) 'ELAPSED WALLCLOCK TIME IN SECONDS: ',
     *            MOD(WAL1-WAL0+1000000.0,1000000.0)
      STOP
C
C.....PRINTING FORMAT
C
 1000 FORMAT('1MOLECULAR DYNAMICS SIMULATION FOR VIBRATING WATER',
     *       ' STARTS AT ',A8,' ON ',A8)
 1001 FORMAT('0TEMPERATURE                = ',F8.2,' K',/,
     1       ' DENSITY                    = ',F8.5,' G/C.C.',/,
     2       ' NUMBER OF MOLECULES        = ',I8,/,
     3       ' TIME STEP                  = ',D8.2,' SEC',/,
     4       ' ORDER USED TO SOLVE F=MA   = ',I8,/,
     5       ' NO. OF TIME STEPS          = ',I8,/,
     6       ' FREQUENCY OF DATA SAVING   = ',I8,/,
     7       ' FREQUENCY TO WRITE RST FILE= ',I8  )
 1002 FORMAT('0SPHERICAL CUTOFF RADIUS    = ',F8.4,' ANGSTROM')
 1003 FORMAT('0***** NEW RUN STARTING FROM RE-START FILE FT',I2,
     *       ' *****')
 1005 FORMAT('0      TOTAL CPU USED =',F10.2,' SECONDS')
C1005 FORMAT('0      TOTAL CPU USED =',F10.2,' SECONDS',/,
C    *       '       TOTAL WALL TIME USED =',F10.2,' SECONDS')
 1006 FORMAT('0***** WARNING: INPUT TEMPERATURE',F8.2,' IS DIFFERENT',
     *       ' FROM RESTART FILE TEMPERATURE',F8.2,' *****')
 1007 FORMAT('0***** WARNING: INPUT DENSITY',F8.2,' IS DIFFERENT',
     *       ' FROM RESTART FILE DENSITY',F8.2,' *****')
 1008 FORMAT('0***** WARNING: INPUT ORDER',F8.2,' IS DIFFERENT',
     *       ' FROM RESTART FILE ORDER',F8.2,' *****')
 1009 FORMAT('0***** WARNING: TIME STEP HAS BEEN CHANGED FROM',E9.2,
     *       ' IN RESTART FILE TO ',E9.2,' IN INPUT')
 1010 FORMAT('0***** STATISTICAL COUNTERS ARE BEEN RESET TO ZEROS')
C
      END
C     Cedar Fortran KAP   (Fri, 16 Nov 1990)  o3,r3,d3 11 Dec 1990 11:22:58
C***********************************************************************
      SUBROUTINE MDMAIN(NFSV,NFRST,NSTEP,NRST,NPRINT,NSAVE,LKT,NORD1)
C***********************************************************************
C
C.....MOLECULAR DYNAMICS LOOP
C
C     IMPLICIT DOUBLE PRECISION(A-H,O-Z),LOGICAL*1($)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (NOMOL=343,MAXODR=7,NATOM=3,MXOD2=MAXODR+2,
     *           NDVAR=NATOM*NOMOL*3*MXOD2)
      COMMON /MDVAR/ TEMP,RHO,TSTEP,BOXL,BOXH,CUTOFF,CUT2,NMOL,NORDER,
     *               NATMO,NATMO3,NMOL1
      COMMON /WATER/ OMAS,HMAS,WTMOL,ROH,ANGLE,FHM,FOM,ROHI,ROHI2,NATOMS
      COMMON /WWPOT/ QQ,A1,B1,A2,B2,A3,B3,A4,B4,AB1,AB2,AB3,AB4,C1,C2,
     *               QQ2,QQ4,REF1,REF2,REF4
      COMMON /CNST / UNITT,UNITL,UNITM,BOLTZ,AVGNO,PCC(10)
      COMMON /MDDATA/ VAR(NDVAR),VM(3*NOMOL),TLC(100),
     *                ELPST,TKIN,TVIR,TTMV,FPOT,FKIN,IX(3*MXOD2),
     *                IRST,NVAR,NXYZ,NXV,IXF,IYF,IZF,IMY,IMZ
      DIMENSION SUM(3)
C
      WRITE (6,1000) IRST,ELPST
      DT=UNITT*TSTEP*NRST
      IF(NSAVE.GT.0) THEN
      WRITE (6,1001) NSAVE
      REWIND NFSV
      WRITE (NFSV) TEMP,RHO,NMOL,NATOMS,TSTEP,NORDER,IRST,NSTEP,NSAVE,
     *             HMAS,OMAS,UNITL,UNITT,UNITM,BOLTZ,AVGNO
      ENDIF
      WRITE (6,1002)
D        call trace_enter (2027) 
      DO 2000 I=1,NSTEP
      TTMV=TTMV+1.0D0
      CALL PREDIC(VAR,TLC,NORD1)
C.....CALL BNDRY(BOXL,VAR,NATMO3,NATOMS)
      CALL INTRAF(VAR,VAR(IXF),VM,VIR)
      CALL INTERF(VAR,VAR(IX(2)),VAR(IX(3)),VAR(IXF),VAR(IYF),VAR(IZF),
     *            VM,VM(IMY),VM(IMZ),VIR)
      CALL CORREC(VAR,VAR(IX(7)),VAR(IXF),NATMO3,PCC,NORD1)
      CALL BNDRY(BOXL,VAR,NATMO3,NATOMS)
      CALL KINETI(VAR(IX(4)),NMOL,SUM,HMAS,OMAS)
      TKIN=TKIN+SUM(1)+SUM(2)+SUM(3)
      TVIR=TVIR-VIR
C.....CHECK IF PRINTING AND/OR SAVING IS TO BE DONE
      IF(MOD(I,NPRINT).EQ.0.OR.NSAVE.GT.0.AND.MOD(I,NSAVE).EQ.0) THEN
      CALL POTENG(VAR,VAR(IX(2)),VAR(IX(3)),VM,VM(IMY),VM(IMZ),
     *            POTA,POTR,POTRF)
      POTA=POTA*FPOT
      POTR=POTR*FPOT
      POTRF=POTRF*FPOT
      XVIR=TVIR*FPOT*0.5D0/TTMV
      AVGT=TKIN*FKIN*TEMP*2.0D0/(3.0D0*TTMV)
      TEN=(SUM(1)+SUM(2)+SUM(3))*FKIN
      XTT=POTA+POTR+POTRF+TEN
      IF(MOD(I,NPRINT).EQ.0)
     *   WRITE (6,1003) I,TEN,POTA,POTR,POTRF,XTT,AVGT,XVIR
      IF (I.EQ.NSTEP) THEN
         WRITE(1,2)
         WRITE(1,5) XTT
         IF (XTT.GT.-3.05 .AND. XTT.LT.-3.) THEN
            WRITE(1,3)
         ELSE
            WRITE(1,4)
         ENDIF
      ENDIF
C
C.....ONLY COORDINATES AND VELOCITIES DATA ARE SAVED
C
      IF(NSAVE.GT.0.AND.MOD(I,NSAVE).EQ.0)
     *        CALL SAVEDT(NFSV,VAR,NXV,TEN,POTA,POTR,AVGT,POTRF,VIR)
C
C.....SAVE INTERMEDIATE DATA FOR RESTART
C
      IF(MOD(I,NRST).NE.0) GO TO 2000
      ELPST=ELPST+DT
C
C.....NORMALIZE KINETIC TO 1.5KT IF LKT=1
C
      IF(LKT.EQ.1) THEN
      CALL NRMLKT(VAR(IX(4)),SUM,FKIN,NATMO,TVIR,TKIN,TTMV,ELPST)
      WRITE (6,1004) TEMP
      ENDIF
C
C.....WRITE RESTART AND SAVED FILES
C
      REWIND NFRST
      WRITE (NFRST) TEMP,RHO,NMOL,TSTEP,NORDER,ELPST,TKIN,TTMV,IRST
      CALL SAVEDT(NFRST,VAR,NXYZ,TVIR,I,NSAVE,D,D,D)
      REWIND NFRST
      WRITE (6,1005) IRST+1,NFRST
      ENDIF
 2000 CONTINUE
C
D        call trace_exit (2028) 
    2 FORMAT(1X,'VALIDATION PARAMETERS:'/)
    3 FORMAT(//1X,'RESULTS FOR THIS RUN ARE:   VALID')
    4 FORMAT(//1X,'RESULTS FOR THIS RUN ARE: INVALID')
    5 FORMAT(1X,F10.6)
 1000 FORMAT('0RESTART',I3,' AFTER',D12.5,' SECONDS')
 1001 FORMAT('0COLLECTING X AND V DATA AT EVERY',
     *       I4,' TIME STEPS ')
 1002 FORMAT('0INTERMEDIATE RESULTS (ENERGY EXPRESSED IN UNITS OF KT',
     *   '/ATOM)',
     *   /,'0       TIME       KINETIC   INTRA POT   INTER POT',
     *   '   REACT POT       TOTAL      <TEMPERATURE>        <VIR>',/)
 1003 FORMAT(5X,I7,F14.5,F12.5,F12.5,F12.5,F12.5,F16.3,F16.5)
 1004 FORMAT('0***** KINETIC ENERGY IS NOW RE-NORMALIZED TO',F8.2,
     *       'K *****',/)
 1005 FORMAT('0***** FILE FOR NEXT RESTART (',I2,') IS',
     *       ' WRITTEN ON FT',I2)
      RETURN
      END
C     Cedar Fortran KAP   (Fri, 16 Nov 1990)  o3,r3,d3 11 Dec 1990 11:22:58
C***********************************************************************
      SUBROUTINE NRMLKT(VAR,SUM,FKIN,NATMO,TVIR,TKIN,TTMV,ELPST)
C***********************************************************************
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION VAR(1),SUM(1)
      INTEGER K1, K2, K3, K4, K5
C
C.....NORMALIZE EACH COMPOENT OF KINETIC ENERGY TO KT/2
C
      KMAX=0
D        call trace_enter (2029) 
      DO 300 J=1,3
      XT=SQRT(0.5D0/(SUM(J)*FKIN))
      KMIN=KMAX+1
      KMAX=KMAX+NATMO
      K2 = (J * NATMO - ((-1) + J) * NATMO + 31) / 32
      CDOALL 300 K1=0,32*K2-32,32
      INTEGER K6, K7
      K7 = ((-1) + J) * NATMO + K1 + 1
      K6 = MIN0 (((-1) + J) * NATMO + K1 + 32, J * NATMO) - (((-1) + J)
     X  * NATMO + K1)
      VAR(K7:$K6) = VAR(K7:$K6) * XT
  300 CONTINUE
D        call trace_exit (2030) 
      TVIR=0.0D0
      TKIN=0.0D0
      TTMV=0.0D0
      ELPST=0.0D0
      RETURN
      END
C     Cedar Fortran KAP   (Fri, 16 Nov 1990)  o3,r3,d3 11 Dec 1990 13:33:58
C***********************************************************************
      SUBROUTINE POTENG(X,Y,Z,XM,YM,ZM,POTA,POTR,PTRF)
C***********************************************************************
C
C.....THIS ROUTINE CALCULATES THE POTENTIAL ENERGY
C     FC11 ,FC12, FC13, AND FC33 ARE THE QUARDRATIC FORCE CONSTANT
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /MDVAR/ TEMP,RHO,TSTEP,BOXL,BOXH,CUTOFF,CUT2,NMOL,NORDER,
     *               NATMO,NATMO3,NMOL1
      COMMON /WATER/OMAS, HMAS, WTMOL, ROH, ANGLE, FHM, FOM, ROHI, ROHI2
     X  , NATOMS
      COMMON /WWPOT/QQ, A1, B1, A2, B2, A3, B3, A4, B4, AB1, AB2, AB3,
     X  AB4, C1, C2, QQ2, QQ4, REF1, REF2, REF4
      COMMON /FRCNST/FC11, FC12, FC13, FC33, FC111, FC333, FC112, FC113
     X  , FC123, FC133, FC1111, FC3333, FC1112, FC1122, FC1113, FC1123,
     X  FC1133, FC1233, FC1333
      DIMENSION X(1), Y(1), Z(1), XL(14), YL(14), ZL(14), RL(14), RS(14)
      DIMENSION XM(1), YM(1), ZM(1)
      INTEGER II1, II2, II3
      INTEGER II5, II6, II11, II21, II31, K1, K2, K3, K4, K5, K6, II32,
     X  II33, II22, II23, II12, II13
C
C.....INTRA-MOLECULAR POTENTIAL ENERGY
C
      POTA = 0.0D0
      IW1 = 1
      IWO = 2
      IW2 = 3
D        call trace_enter (2031) 
      CDOACROSS 1000 I=1,NMOL
      DOUBLE PRECISION DTS, DRP, DR2S, DR1S, DR2, DR1, DT, COS, RX, R2,
     X  R1, potal
      XM(I) = C1 * X(((-1)+I)*NATOMS+2) + C2 * (X(((-1)+I)*NATOMS+1) + X
     X  (((-1)+I)*NATOMS+3))
      YM(I) = C1 * Y(((-1)+I)*NATOMS+2) + C2 * (Y(((-1)+I)*NATOMS+1) + Y
     X  (((-1)+I)*NATOMS+3))
      ZM(I) = C1 * Z(((-1)+I)*NATOMS+2) + C2 * (Z(((-1)+I)*NATOMS+1) + Z
     X  (((-1)+I)*NATOMS+3))
      R1 = (X(((-1)+I)*NATOMS+2) - X(((-1)+I)*NATOMS+1)) ** 2 + (Y(((-1)
     X  +I)*NATOMS+2) - Y(((-1)+I)*NATOMS+1)) ** 2 + (Z(((-1)+I)*NATOMS+
     X  2) - Z(((-1)+I)*NATOMS+1)) ** 2
      R2 = (X(((-1)+I)*NATOMS+2) - X(((-1)+I)*NATOMS+3)) ** 2 + (Y(((-1)
     X  +I)*NATOMS+2) - Y(((-1)+I)*NATOMS+3)) ** 2 + (Z(((-1)+I)*NATOMS+
     X  2) - Z(((-1)+I)*NATOMS+3)) ** 2
      RX = (X(((-1)+I)*NATOMS+2) - X(((-1)+I)*NATOMS+1)) * (X(((-1)+I)*
     X  NATOMS+2) - X(((-1)+I)*NATOMS+3)) + (Y(((-1)+I)*NATOMS+2) - Y((
     X  (-1)+I)*NATOMS+1)) * (Y(((-1)+I)*NATOMS+2) - Y(((-1)+I)*NATOMS+3
     X  )) + (Z(((-1)+I)*NATOMS+2) - Z(((-1)+I)*NATOMS+1)) * (Z(((-1)+I)
     X  *NATOMS+2) - Z(((-1)+I)*NATOMS+3))
      R1 = SQRT (R1)
      R2 = SQRT (R2)
      COS = RX / (R1 * R2)
      DT = (ACOS (COS) - ANGLE) * ROH
      DR1 = R1 - ROH
      DR2 = R2 - ROH
      DR1S = DR1 * DR1
      DR2S = DR2 * DR2
      DRP = DR1 + DR2
      DTS = DT * DT
      potal = (FC11 * (DR1S + DR2S) + FC33 * DTS) * 0.5D0 + FC12 *
     X      DR1 * DR2 + FC13 * DRP * DT + (FC111 * (DR1S * DR1 + DR2S *
     X     DR2) + FC333 * DTS * DT + FC112 * DRP * DR1 * DR2 + FC113 * (
     X     DR1S + DR2S) * DT + FC123 * DR1 * DR2 * DT + FC133 * DRP *
     X     DTS) * ROHI
      potal = potal + (FC1111 * (DR1S * DR1S + DR2S * DR2S) + FC3333 * DTS
     X      * DTS + FC1112 * (DR1S + DR2S) * DR1 * DR2 + FC1122 * DR1S *
     X      DR2S + FC1113 * (DR1S * DR1 + DR2S * DR2) * DT + FC1123 *
     X     DRP * DR1 * DR2 * DT + FC1133 * (DR1S + DR2S) * DTS + FC1233
     X     * DR1 * DR2 * DTS + FC1333 * DRP * DTS * DT) * ROHI2
      CALL AWAIT (1, 1, POTA)
      POTA = POTA + potal
      CALL ADVANCE (1, POTA)
 1000 CONTINUE
C
C.....INTER-MOLECULAR POTENTIAL
C
D        call trace_exit (2032) 
      IW1=1
      IWO=2
      IW2=3
      POTR=0.0D0
      PTRF=0.0D0
D        call trace_enter (2033) 
      CDOACROSS 2000 I=1,NMOL1
      integer iw1,iwo,iw2,jw1,jwo,jw2,kc
      INTEGER II1, II2, II3, j, k
      DOUBLE PRECISION XL(14), YL(14), ZL(14), RL(14), RS(14), potrl, ptrfl

      potrl = 0.0D0
      ptrfl = 0.0D0
      IW1 = 1+ (I-1)*NATOMS
      IWO = IW1 + 1
      IW2 = IW1 + 2
      JW1=IW1+NATOMS
      JWO=JW1+1
      JW2=JW1+2
      JMIN=I+1
      DO 2100 J=I+1,NMOL
      XL(1) = XM(I) - XM(J)
      XL(2) = XM(I) - X(JW1)
      XL(3) = XM(I) - X(JW1+2)
      XL(4) = X(IW1) - XM(J)
      XL(5) = X(IW1+2) - XM(J)
      XL(6) = X(IW1) - X(JW1)
      XL(7) = X(IW1) - X(JW1+2)
      XL(8) = X(IW1+2) - X(JW1)
      XL(9) = X(IW1+2) - X(JW1+2)
      XL(10) = X(IW1+1) - X(JW1+1)
      XL(11) = X(IW1+1) - X(JW1)
      XL(12) = X(IW1+1) - X(JW1+2)
      XL(13) = X(IW1) - X(JW1+1)
      XL(14) = X(IW1+2) - X(JW1+1)
      WHERE (ABS (XL(1:$14)) .GT. BOXH) XL(1:$14) = XL(1:$14) - SIGN (
     X  BOXL, XL(1:$14))
    2 CONTINUE
      YL(1) = YM(I) - YM(J)
      YL(2) = YM(I) - Y(JW1)
      YL(3) = YM(I) - Y(JW1+2)
      YL(4) = Y(IW1) - YM(J)
      YL(5) = Y(IW1+2) - YM(J)
      YL(6) = Y(IW1) - Y(JW1)
      YL(7) = Y(IW1) - Y(JW1+2)
      YL(8) = Y(IW1+2) - Y(JW1)
      YL(9) = Y(IW1+2) - Y(JW1+2)
      YL(10) = Y(IW1+1) - Y(JW1+1)
      YL(11) = Y(IW1+1) - Y(JW1)
      YL(12) = Y(IW1+1) - Y(JW1+2)
      YL(13) = Y(IW1) - Y(JW1+1)
      YL(14) = Y(IW1+2) - Y(JW1+1)
      WHERE (ABS (YL(1:$14)) .GT. BOXH) YL(1:$14) = YL(1:$14) - SIGN (
     X  BOXL, YL(1:$14))
    3 CONTINUE
      ZL(1) = ZM(I) - ZM(J)
      ZL(2) = ZM(I) - Z(JW1)
      ZL(3) = ZM(I) - Z(JW1+2)
      ZL(4) = Z(IW1) - ZM(J)
      ZL(5) = Z(IW1+2) - ZM(J)
      ZL(6) = Z(IW1) - Z(JW1)
      ZL(7) = Z(IW1) - Z(JW1+2)
      ZL(8) = Z(IW1+2) - Z(JW1)
      ZL(9) = Z(IW1+2) - Z(JW1+2)
      ZL(10) = Z(IW1+1) - Z(JW1+1)
      ZL(11) = Z(IW1+1) - Z(JW1)
      ZL(12) = Z(IW1+1) - Z(JW1+2)
      ZL(13) = Z(IW1) - Z(JW1+1)
      ZL(14) = Z(IW1+2) - Z(JW1+1)
      WHERE (ABS (ZL(1:$14)) .GT. BOXH) ZL(1:$14) = ZL(1:$14) - SIGN (
     X  BOXL, ZL(1:$14))
    4 CONTINUE
      KC = 0
      RS(1:$9) = XL(1:$9) * XL(1:$9) + YL(1:$9) * YL(1:$9) + ZL(1:$9) *
     X  ZL(1:$9)
 2110 KC = KC + COUNT (MASK=RS(1:$9) .GT. CUT2)
      IF (KC .NE. 9) THEN
      WHERE (RS(1:$9) .LE. CUT2)
      RL(1:$9) = SQRT (RS(1:$9))
   10 OTHERWISE
      RL(1:$9) = CUTOFF
      RS(1:$9) = CUT2
      END WHERE
      POTRl=POTRl-QQ2/RL(2)-QQ2/RL(3)-QQ2/RL(4)-QQ2/RL(5)
     *         +QQ /RL(6)+QQ /RL(7)+QQ /RL(8)+QQ /RL(9)
     *         +QQ4/RL(1)
      PTRFl=PTRFl-REF2*RS(1)-REF1*((RS(6)+RS(7)+RS(8)+RS(9))*0.5D0
     *                           -RS(2)-RS(3)-RS(4)-RS(5))
      IF (KC .LE. 0) THEN
 2130 RL(10:$5) = SQRT (XL(10:$5) * XL(10:$5) + YL(10:$5) * YL(10:$5) +
     X  ZL(10:$5) * ZL(10:$5))
      POTRl=POTRl+A1* EXP(-B1*RL(10))
     *         +A2*(EXP(-B2*RL( 6))+EXP(-B2*RL( 7))
     *                             +EXP(-B2*RL( 8))+EXP(-B2*RL( 9)))
     *         +A3*(EXP(-B3*RL(11))+EXP(-B3*RL(12))
     *                             +EXP(-B3*RL(13))+EXP(-B3*RL(14)))
     *         -A4*(EXP(-B4*RL(11))+EXP(-B4*RL(12))
     *                             +EXP(-B4*RL(13))+EXP(-B4*RL(14)))
      END IF
      END IF
  100 JW1 = JW1 + NATOMS
      JWO = JWO + NATOMS
      JW2 = JW2 + NATOMS
 2100 CONTINUE
      CALL AWAIT (1, 1, POTR,PTRF)
      POTR = POTR + potrl
      PTRF = PTRF + ptrfl
      CALL ADVANCE (1, POTR,PTRF)
      
 2000 CONTINUE
D        call trace_exit (2034) 
      RETURN
      END
C     Cedar Fortran KAP   (Fri, 15 Feb 1991)  o3,r1,d3 11 Mar 1991 08:28:20
C***********************************************************************
      SUBROUTINE PREDIC(V,C,NOR1)
C***********************************************************************
C
C.....THIS ROUTINE CALCULATES PREDICTED F(X), F'(X), F''(X), ......
C     (WITH DELTAT**N/N] INCLUDED)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
       COMMON /MDVAR/TEMP, RHO, TSTEP, BOXL, BOXH, CUTOFF, CUT2,
     X  NMOL, NORDER, NATMO, NATMO3, NMOL1
      DIMENSION V(1),C(1) 
      INTEGER II1, II2, L1
C      GLOBAL NOR1, C, V, s, length, lenV
C
       DIMENSION s(NATMO3*NORDER)
D        call trace_enter (2035) 
      CDOALL 1000 I=1,NORDER
      integer l1,k,ii1
      L1 = NORDER - I + 1

      DO 1100 K=1,NATMO3
      II1 = ((-1) + I) * NOR1 + 2

      S(natmo3*(i-1)+k)=  DOTPRODUCT (C(II1:$L1), 
     x V(K+MAX0 (NATMO3, 0)*((-1)+I)+ NATMO3:$L1:NATMO3))
 1100 CONTINUE

 1000 CONTINUE

      length = NORDER*natmo3
      lenV = (length+31)/32 
      CDOALL 2000 I=1,length,lenV
      integer high
      high = min(i+lenV-1,length)
      v(i:high)=v(i:high)+s(i:high)
 2000 continue
 2100 continue
D        call trace_exit (2036) 
      RETURN
      END
C     Cedar Fortran KAP   (Fri, 16 Nov 1990)  o3,r3,d3 11 Dec 1990 11:22:58
C***********************************************************************
      FUNCTION RAND(N)
C***********************************************************************
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /RANDNO/ R3(127),R1,I2
      SAVE S, T, RMC
      SAVE IW
      INTEGER I1
      DATA S, T, RMC /0.0D0, 1.0D0, 1.0D0/
      DATA IW /-1/
C
      IF (R1 .GE. 1.0D0 .OR. N .NE. 0) THEN
      IF (IW .LE. 0) THEN
   10 IW=IW+1
      T=T*0.5D0
      R1=S
      S=S+T
      IF(S.GT.R1.AND.S.LT.1.0D0) GO TO 10
      IKT=(IW-1)/12
      IC=IW-12*IKT
      ID=2**(13-IC)
      I1 = IW - 12 * ((IW - 1) / 12)
   20 RMC = RMC * PRODUCT (SPREAD (0.5D0, 1, I1))
      RM=0.015625D0*0.015625D0
      END IF
   30 I2=127
      IR=MOD(ABS(N),8190)+1
   40 R1=0.0D0
D        call trace_enter (2037) 
      DO 50 I=1,IKT
      IR = MOD (17 * IR, 8191)
   50 R1=(R1+REAL(IR/2))*RM
D        call trace_exit (2038) 
      IR=MOD(17*IR,8191)
      R1=(R1+REAL(IR/ID))*RMC
      R3(I2)=R1
      I2=I2-1
      IF(I2.GT.0) GO TO 40
      END IF
   60 IF(I2.EQ.0) I2=127
      T=R1+R3(I2)
      IF(T.GE.1.0D0) T=(R1-0.5D0)+(R3(I2)-0.5D0)
      R1=T
      R3(I2)=R1
      I2=I2-1
      RAND=R1
      RETURN
      END
C     Cedar Fortran KAP   (Fri, 16 Nov 1990)  o3,r3,d3 11 Dec 1990 11:22:58
C***********************************************************************
      SUBROUTINE SAVEDT(NF,V,NV,TENG,POTA,POTR,ATEMP,POTRF,VIR)
C***********************************************************************
C
C.....THIS ROUTINE SAVE VARIABLES V(NV),TEMP,POT,ATEMP,POTRF,VIR
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION V(NV)
C
      WRITE (NF) TENG,POTA,POTR,ATEMP,NV,V,POTRF,VIR
      RETURN
      END
C     Cedar Fortran KAP   (Fri, 16 Nov 1990)  o3,r3,d3 11 Dec 1990 11:22:58
C***********************************************************************
      SUBROUTINE SYSCNS
C***********************************************************************
C
C.....CONSTANTS FOR THE SIMULATION SYSTEM
C
C     IMPLICIT DOUBLE PRECISION(A-H,O-Z),LOGICAL*1($)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (NOMOL=343,MAXODR=7,NATOM=3,MXOD2=MAXODR+2,
     *           NDVAR=NATOM*NOMOL*3*MXOD2)
      COMMON /MDVAR/ TEMP,RHO,TSTEP,BOXL,BOXH,CUTOFF,CUT2,NMOL,NORDER,
     *               NATMO,NATMO3,NMOL1
      COMMON /WATER/ OMAS,HMAS,WTMOL,ROH,ANGLE,FHM,FOM,ROHI,ROHI2,NATOMS
      COMMON /WWPOT/ QQ,A1,B1,A2,B2,A3,B3,A4,B4,AB1,AB2,AB3,AB4,C1,C2,
     *               QQ2,QQ4,REF1,REF2,REF4
      COMMON /CNST / UNITT,UNITL,UNITM,BOLTZ,AVGNO,PCC(10)
      COMMON /MDDATA/ VAR(NDVAR),VM(3*NOMOL),TLC(100),
     *                ELPST,TKIN,TVIR,TTMV,FPOT,FKIN,IX(3*MXOD2),
     *                IRST,NVAR,NXYZ,NXV,IXF,IYF,IZF,IMY,IMZ
C
      TSTEP=TSTEP/UNITT
      NATMO=NATOMS*NMOL
      NATMO3=NATMO*3
      FPOT=UNITM*(UNITL/UNITT)**2/(BOLTZ*TEMP*NATMO)
      FKIN=FPOT*0.5D0/(TSTEP*TSTEP)
      BOXL=(NMOL*WTMOL*UNITM/RHO)**(1.0D0/3.0D0)
      BOXL=BOXL/UNITL
      BOXH=BOXL*0.5D0
      CUTOFF=MAX(BOXH,CUTOFF)
      REF1=-QQ/(CUTOFF*CUTOFF*CUTOFF)
      REF2=2.0D0*REF1
      REF4=2.0D0*REF2
      CUT2=CUTOFF*CUTOFF
      FHM=(TSTEP*TSTEP*0.5D0)/HMAS
      FOM=(TSTEP*TSTEP*0.5D0)/OMAS
      NMOL1=NMOL-1
      RETURN
      END