C ----------------------------------------------------------------
C       PROGRAM ICE (FLEXIBLE MODEL) JUL.12,1995 LIJU.                       
C       1. USE EWALD SUMMATION.                                  
C       2. ABLE TO CONTINUE FROM LAST CONFIGURATION
C       3. WITH HYDROGEN ATOM VELOCITY AUTOCORRELATION FUNCTIONAL
C       4. START WITH ZERO-DIPOLE SELF-CONSISTENT CONFIGURATION 
C       5. NEW AHARMONIC INTRAMOLECULAR POTENTIAL
C ----------------------------------------------------------------
C       USE THE FOLLOWING UNITS:                                 
C       ENERGY UNIT = THE ELECTROSTATIC ENERGY BETWEEN TWO        
C                     ELECTRONS WITH 1 ANGSTROMS APART           
C                   = 2.3070795E-18 JOULE                           
C       MASS   UNIT = 1 AMU = 1.66054E-27 KG                      
C       LENGTH UNIT = 1 ANGSTROM                                 
C       TIME   UNIT = 2.6828315E-15 SECOND = 2.6828315E-3 PS          
C ----------------------------------------------------------------
      
      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NPD = 648)
      PARAMETER (ID  = 50000)
      PARAMETER (IV  = 300)
      COMMON /EWALD/ ALPHA,ERR(ID),CLI(ID),DERR,A(9,9,9)
      COMMON /POS/ X(NPD),Y(NPD),Z(NPD),X1(NPD),Y1(NPD),Z1(NPD)
      DOUBLE PRECISION X1,Y1,Z1
      COMMON /FORCE/ FX(NPD),FY(NPD),FZ(NPD),A12,A6
      COMMON /PARM/ RDEL,IGOFR(3,1500),NM,NPO,NPH,NP,PMASS(NPD),
     1              ITYPE(NPD),Q(NPD),VOLUME,CUBE,FROG
      COMMON /ENERGY/ TOTE,POTE,KINE,PLJ12,PLJ6,ELEC,VINTRA,VIRIAL
      COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2
      COMMON /CONTROL/ KB,KSTART,SCALER,SEED,DELTA,BON,ANG,PM2
      COMMON /ADDENDUM/ PMX,PMY,PMZ

      REAL KINE,MO,MH
      
      REAL MSD,KINESUM
      DIMENSION V2(IV,216),VX(IV,216),VY(IV,216),VZ(IV,216),CORR(5000),
     1 	        X2(NPD),Y2(NPD),Z2(NPD),DX(NPD),DY(NPD),DZ(NPD)
C     2          VMX(216),VMY(216),VMZ(216)

      DATA LP,LP_GR,LP_CORR,LP_POSI,LP_CONFIG,LP_BAK /25,26,27,28,29,30/
      
C     SELECT 300 TIME ORIGINS, EACH WITH 10 STEPS APART
C     TOTAL CORRELATION LENGTH = 5000 STEPS = 2 ps
 
      KVSTEP = 3000/IV
      AVO=6.0221367E+23
      BOLZ=1.380658E-23
      AMU=1.6605E-24
      
C     ONE TIME STEP = 0.0004ps
      DELTA  = 0.149096206
      NM = 216
      NPO = 1
      NPH = 2
      NP = NM*(NPO+NPH)

C     BERENDSEN T,P COUPLING CONSTANT
      TAUT = 0.05/2.6828315E-3
      TAUP = 0.5/2.6828315E-3
      BETA = 4.9E-5
      FACTORT = 1.
      FACTORP = 1.

C     A12   = 1895.294007
C     A6    = -1.883540782
C     QO    = -0.8476

      PI = 3.141592654
      P2 = 2.*PI
      READ *,A12,A6,QO
      QH    = -QO/2.
      QOO   = QO*QO
      QOH   = QO*QH
      QHH   = QH*QH
      MO = 15.9994
      MH = 1.00797
      AMASS = MO+2.*MH
      HPART = MH/AMASS
      OPART = MO/AMASS
      BOND  = 0.9572
      ANGLE = (104+52./100.)/180.*PI
      HHDIST= 2.*SIN(ANGLE/2.)*BOND
      
C     RANDOM NUMBER GENERATOR SEED
      SEED  = 0.2573
      
C     WATER OXYGEN AND HYDROGEN
      DO 1100 I = 1,NP-2,3
      PMASS(I) = MO
      PMASS(I+1)=MH
      PMASS(I+2)=MH
      ITYPE(I) = 0
      ITYPE(I+1)=1
      ITYPE(I+2)=1
      Q(I) = QO
      Q(I+1) = QH
      Q(I+2) = QH
1100  CONTINUE
      
      READ *,T,P,MAXKB,KSTART
      READ *,LAST,K_CORR,K_POSI
      READ *,DENSITY

      TOTALMASS = NM*AMASS
      VOLUME = TOTALMASS*1.66054/DENSITY
      CUBE = VOLUME**(1./3.)
      FROG = P2/CUBE
C     SCALER IS THE IDEAL TOTAL KINETIC ENERGY OF CM
      SCALER = 1.5*BOLZ*T*NM/2.3070795E-18     

      IF (LAST.EQ.0) THEN
       CALL DIAMOND
       CALL INTVEL
      ELSE
       DO 377 I=1,NP
       READ *,X(I),Y(I),Z(I),X1(I),Y1(I),Z1(I)
377    CONTINUE
      ENDIF
      
      RDEL = CUBE/2./500.
      DO 153 I=1,1500
      IGOFR(1,I) = 0
      IGOFR(2,I) = 0
      IGOFR(3,I) = 0
 153  CONTINUE
      DO 154 I = 1,5000
      CORR(I) = 0.
 154  CONTINUE

      OPEN (UNIT=LP,        FORM='FORMATTED',FILE="temp.out")
      OPEN (UNIT=LP_CONFIG, FORM='FORMATTED',FILE="config")
      OPEN (UNIT=LP_BAK,    FORM='FORMATTED',FILE="config.bak")
      OPEN (UNIT=LP_GR,     FORM='FORMATTED',FILE="gr.out")
      IF (K_CORR.EQ.1) 
     1   OPEN (UNIT=LP_CORR, FORM='FORMATTED',FILE="corr.out")
      IF (K_POSI.EQ.1) 
     1   OPEN (UNIT=LP_POSI, FORM='FORMATTED',FILE="posi.out")
      
C     -------------------------------------------------------------
C     KWRITE IS THE FREQUENCY OF WRITING ON TEMP.OUT
C     KWBAK IS THE FREQUENCY OF WRITING ON CONFIG.BAK
C     KWGR IS THE FREQUENCY OF WRITING ON GR.OUT
C     --------------------------------------------------------------
      KWRITE = 20
      KWGR   = 20000
      KWBAK  = 10160
     
C     SET CONSTS AND TABLES FOR EWALD SUMMATION
      ALPHA = 5.714/CUBE
      DERR = ALPHA*CUBE/2.*SQRT(3.)*1.2
      DERR = DERR/ID
      C = DERR
      WW = 2./SQRT(PI)
      
      DO 500 I=1,ID
      ERR(I) = ERFC(C)
      CLI(I) = WW*EXP(-C*C)
      C = C+DERR
500   CONTINUE
      
C     ----------------------------------------------------------
C     PRINT PARAMETERS.
      
 5    WRITE(LP,10)
10    FORMAT(///)
      WRITE(LP,20)
20    FORMAT(7X,50('*'))
      WRITE(LP,30)
30    FORMAT(7X,'*',T57,'*')
      WRITE(LP,40)
40    FORMAT(7X,'*',3X,'MOLECULAR DYNAMICS SIMULATION USING SPC/F',
     A T57,'*')
      WRITE(LP,50) NM
50    FORMAT(7X,'*',3X,'FLEXIBLE MODEL FOR',I4,
     A ' WATER MOLECULES,',T57,'*')
      WRITE(LP,70)
70    FORMAT(7X,'*',3X,'WITH DIAMOND CUBIC STRUCTURE.',
     A       T57,'*')
      WRITE(LP,30)
      WRITE(LP,20)
      WRITE(LP,30)
      WRITE(LP,80) T,P
80    FORMAT(7X,'*',2X,'T(K)   =',F8.3,T30,'P(BAR) = ',
     1       F7.4,T57,'*')
      WRITE(LP,90) MAXKB,DELTA*2.68279E-3
90    FORMAT(7X,'*',2X,'MAXKB  =',I6,T30,  'TIME STEP(ps) = ',
     1       F7.4,T57,'*')
      WRITE(LP,30)
      WRITE(LP,111) QO,QH
111   FORMAT(7X,'*',2X,'QO  = ',F7.4,T30,'QH  = ',F7.4,T57,'*')
      WRITE(LP,112) A12,A6
112   FORMAT(7X,'*',2X,'A12 = ',F11.6,T30,'A6  = ',F11.7,T57,'*')
      WRITE(LP,30)
      WRITE(LP,110) 9*9*9-1
110   FORMAT(7X,'*',2X,'EWALD SUMMATION OVER ',I3,
     1       ' RECIPROCAL LATTICES',T57,'*')
      WRITE(LP,30)
      WRITE(LP,20)
      WRITE(LP,190)
190   FORMAT(///,4X,'KB',5X,'TIME',4X,'TEMP',7X,'PRES',4X,'DENSITY',
     1           4X,'POTE',5X,'TOTE',7X,'MSD')

      DO 55 I=1,NP
      DX(I) = 0.
      DY(I) = 0.
      DZ(I) = 0.
      X2(I) = 0.
      Y2(I) = 0.
      Z2(I) = 0.
55    CONTINUE

      TSUM    = 0.
      VINSUM  = 0.
      KINESUM = 0.
      POTESUM = 0.
      PRESUM  = 0.
      BONSUM  = 0.
      ANGSUM  = 0.
      PM2SUM  = 0.
      PMXSUM  = 0.
      PMYSUM  = 0.
      PMZSUM  = 0.
      DENSUM  = 0.
   
C     -----------------------------------------------------------
C     ENTER MAIN LOOP OF SIMULATION  

      DO 310 KB=1,MAXKB

      CUBE = CUBE*FACTORP
      FROG = FROG/FACTORP
      ALPHA = ALPHA/FACTORP
      VOLUME = VOLUME*FACTORP**3
      DENSITY = DENSITY/FACTORP**3

      DO 320 I=1,NP
      DXX = DELTA*(X1(I)+0.5*DELTA*X2(I))
      DYY = DELTA*(Y1(I)+0.5*DELTA*Y2(I))
      DZZ = DELTA*(Z1(I)+0.5*DELTA*Z2(I))
        
      DX(I) = DX(I) + DXX
      DY(I) = DY(I) + DYY
      DZ(I) = DZ(I) + DZZ

      X(I) = X(I)+DXX
      Y(I) = Y(I)+DYY
      Z(I) = Z(I)+DZZ

      IF (X(I).GT.CUBE) X(I)=X(I)-CUBE
      IF (X(I).LT.0)    X(I)=X(I)+CUBE
      IF (Y(I).GT.CUBE) Y(I)=Y(I)-CUBE
      IF (Y(I).LT.0)    Y(I)=Y(I)+CUBE
      IF (Z(I).GT.CUBE) Z(I)=Z(I)-CUBE
      IF (Z(I).LT.0)    Z(I)=Z(I)+CUBE

      X(I) = X(I)*FACTORP
      Y(I) = Y(I)*FACTORP
      Z(I) = Z(I)*FACTORP
 
320   CONTINUE

      CALL INTERACT 
  
      KINE = 0.
      DO 330 I=1,NP
      XX2 = FX(I)/PMASS(I) 
      YY2 = FY(I)/PMASS(I) 
      ZZ2 = FZ(I)/PMASS(I) 
      X1(I) = FACTORT*X1(I)+DELTA*(X2(I)+XX2)/2.
      Y1(I) = FACTORT*Y1(I)+DELTA*(Y2(I)+YY2)/2.
      Z1(I) = FACTORT*Z1(I)+DELTA*(Z2(I)+ZZ2)/2.
      X2(I) = XX2
      Y2(I) = YY2
      Z2(I) = ZZ2
      KINE = KINE+PMASS(I)*(X1(I)*X1(I)+Y1(I)*Y1(I)+Z1(I)*Z1(I))
330   CONTINUE
      KINE = KINE/2.

      MSD   = 0.
      IP    = 1
      DO 66 I=1,NM
C      VMX(I) = X1(IP)*OPART+(X1(IP+1)+X1(IP+2))*HPART 
C      VMY(I) = Y1(IP)*OPART+(Y1(IP+1)+Y1(IP+2))*HPART 
C      VMZ(I) = Z1(IP)*OPART+(Z1(IP+1)+Z1(IP+2))*HPART
      XX     = DX(IP)*OPART+(DX(IP+1)+DX(IP+2))*HPART 
      YY     = DY(IP)*OPART+(DY(IP+1)+DY(IP+2))*HPART 
      ZZ     = DZ(IP)*OPART+(DZ(IP+1)+DZ(IP+2))*HPART 
      MSD    = MSD+XX*XX+YY*YY+ZZ*ZZ
      IP     = IP+3
66    CONTINUE
      MSD   = MSD/NM

      IF (MOD(KB,KWGR).EQ.0.OR.KB.EQ.MAXKB) THEN
C     ------------- PAIR DISTRIBUTION FUCTIONAL -------------
      IF (KB.EQ.MAXKB) KWGR = MOD(MAXKB-1,KWGR)+1
C     CALC PAIR DISTRIBUTION FUNCTIONS G(R)s 
      DO 475 I = 1,500
      DV = 4.*PI*RDEL**3*I*I
      OP = FLOAT(NM*NPO)
      HP = FLOAT(NM*NPH)
      GROO = 2.*VOLUME*FLOAT(IGOFR(1,I))/(DV*KWGR*OP*OP)
      GRHH = 2.*VOLUME*FLOAT(IGOFR(3,I))/(DV*KWGR*HP*HP)
      GROH =    VOLUME*FLOAT(IGOFR(2,I))/(DV*KWGR*OP*HP)
      WRITE(LP_GR,4100) I,RDEL*I,GROO,GRHH,GROH
4100  FORMAT(2X,I4,4(2X,F7.4))
475   CONTINUE

      IF (KB.NE.MAXKB) THEN
      WRITE(LP_GR,4200)
4200  FORMAT('B')
      DO 476 I = 1,1500
      IGOFR(1,I) = 0
      IGOFR(2,I) = 0
      IGOFR(3,I) = 0
476   CONTINUE
      ENDIF

      ENDIF
C     ------------- PAIR DISTRIBUTION FUCTIONAL -------------

C     PRESSURE WILL BE IN BARS & ENERGIES WILL BE IN KJOULE/MOL.
      PRES = (VIRIAL+2.*KINE)/3./VOLUME*2.3070795E+7
      TNOW = KINE *T/(SCALER*3)
      KINE = KINE *2.3070795E-18*AVO/NM/1000.
      POTE = POTE *2.3070795E-18*AVO/NM/1000.
      VINTRA = VINTRA*2.3070795E-18*AVO/NM/1000.
    
C     ACCUMMULATE SUMS FOR PROPERTY AVERAGES.
      IF (KB.GT.KSTART) THEN 
      BONSUM  = BONSUM  + BON
      ANGSUM  = ANGSUM  + ANG
      PRESUM  = PRESUM  + PRES
      TSUM    = TSUM    + TNOW
      POTESUM = POTESUM + POTE
      KINESUM = KINESUM + KINE
      VINSUM  = VINSUM  + VINTRA
      PMXSUM  = PMXSUM  + PMX
      PMYSUM  = PMYSUM  + PMY
      PMZSUM  = PMZSUM  + PMZ
      PM2SUM  = PM2SUM  + PM2
      DENSUM  = DENSUM  + DENSITY
      ENDIF

      IF(MOD(KB,KWRITE).EQ.0) THEN
C     --------------------------------------------------------
C     PROPERTY CALCULATION & PRINT-OUT AT INTERVALS OF KWRITE 
C     REAL TIME IS IN ps
      REALTIME = DELTA*KB*2.6828315E-3
44    WRITE(LP,270) KB,REALTIME,TNOW,PRES,DENSITY,POTE,KINE+POTE,MSD
270   FORMAT(I6,2(F9.3),F11.3,F9.3,2(F10.3),F9.2)
      ENDIF
C     --------------------------------------------------------

      FACTORT = 1+DELTA/2./TAUT*(T/TNOW-1)
      FACTORP = 1+BETA*DELTA/3./TAUP*(PRES-P)

      IF (K_CORR.EQ.1) THEN
C     --------- VELOCITY CORRELATION FUNCTIONAL --------------
C     CALCULATE VELOCITY AUTOCORRELATION FUNCTION 
      I = KB+8000-MAXKB
      KPT = I/KVSTEP

      IF (MOD(I,KVSTEP).EQ.0.AND.KPT.GE.1.AND.KPT.LE.IV) THEN
      IP = 2
      DO 1300 IM = 1,NM
      VX(KPT,IM) = X1(IP)
      VY(KPT,IM) = Y1(IP)
      VZ(KPT,IM) = Z1(IP)
      V2(KPT,IM) = X1(IP)*X1(IP)+Y1(IP)*Y1(IP)+Z1(IP)*Z1(IP)
      IP = IP+3
1300  CONTINUE
      ENDIF

      IF (KPT.GT.IV) KPT = IV

      DO 1600 KP = 1,KPT
      J = I-KP*KVSTEP
      IF (J.GE.1.AND.J.LE.5000) THEN 
      IP = 2
      DO 1400 IM = 1,NM
      CORR(J) = CORR(J)+(X1(IP)*VX(KP,IM)+Y1(IP)*VY(KP,IM)+
     1          Z1(IP)*VZ(KP,IM))/V2(KP,IM)
      IP = IP+3
1400  CONTINUE
      ENDIF

1600  CONTINUE
      ENDIF
C     --------- VELOCITY CORRELATION FUNCTIONAL --------------


      IF (MOD(KB,KWBAK).EQ.0) THEN
      WRITE (LP_BAK,221) DENSITY
      DO 284 I=1,NP
      WRITE (LP_BAK,222) X(I),Y(I),Z(I),X1(I),Y1(I),Z1(I)
 284  CONTINUE
      WRITE (LP_BAK,4200)
      ENDIF

      print *,'kb = ',kb
      print *,'pote =',pote
      print *,'tote =',pote+kine
      print *,'T  =  ',tnow
      print *,'pres = ',pres
      print *,'density = ',density
      print *,' '

310   CONTINUE 

C               END OF MAIN LOOP    
C     -----------------------------------------------------------

      IF (K_CORR.EQ.1) THEN
      DO 1500 I = 1,5000
      CORR(I) = CORR(I)/NM/IV
      WRITE(LP_CORR,113) CORR(I)
113   FORMAT(F9.6)
1500  CONTINUE
      ENDIF
      
      IF (K_POSI.EQ.1) THEN
      DO 1700 I = 1,NP-2,3
      WRITE(LP_POSI,233) X(I),Y(I),Z(I),X(I+1),Y(I+1),Z(I+1),
     1                   X(I+2),Y(I+2),Z(I+2) 
233   FORMAT(9(F11.5))
1700  CONTINUE
      ENDIF
    
      TSUM    = TSUM/(MAXKB-KSTART)
      ANGSUM  = ANGSUM/(MAXKB-KSTART)
      BONSUM  = BONSUM/(MAXKB-KSTART)
      PRESUM  = PRESUM/(MAXKB-KSTART)
      VINSUM  = VINSUM/(MAXKB-KSTART)
      KINESUM = KINESUM/(MAXKB-KSTART)
      POTESUM = POTESUM/(MAXKB-KSTART)
      DENSUM  = DENSUM/(MAXKB-KSTART)
      PMXSUM  = PMXSUM/(MAXKB-KSTART)
      PMYSUM  = PMYSUM/(MAXKB-KSTART)
      PMZSUM  = PMZSUM/(MAXKB-KSTART)
      PM2SUM  = PM2SUM/(MAXKB-KSTART)

      DIE = (PM2SUM-PMXSUM*PMXSUM-PMYSUM*PMYSUM-PMZSUM*PMZSUM)*QH*QH
      DIE = 1.+DIE*4.*PI*2.3070795E-18/(3*BOLZ*TSUM)/VOLUME

      WRITE (LP,501) KSTART
 501  FORMAT(//,5X,'  THE RUNNING AVERAGE BEGINS FROM STEP',I6)
      WRITE (LP,601) MAXKB
 601  FORMAT(5X,   '                             TO   STEP',I6,/)
      WRITE (LP,701) TSUM
 701  FORMAT(5X,'  THE AVERAGED TEMPERATURE IS ',F7.2,' KELVIN')
      WRITE (LP,1001) PRESUM
 1001 FORMAT(5X,'  THE AVERAGED PRESSURE IS',F10.2,' BAR')
      WRITE (LP,1011) DENSUM
 1011 FORMAT(5X,'  THE AVERAGED DENSITY IS ',F7.5,' G/CM^3')
      WRITE (LP,801) BONSUM
 801  FORMAT(5X,'  THE AVERAGED BOND LENGTH IS ',F7.5,' ANGSTROM')
      WRITE (LP,901) ANGSUM*180/PI
 901  FORMAT(5X,'  THE AVERAGED BOND ANGLE  IS ',F7.3,' DEGREE',/)
      WRITE (LP,1101) KINESUM
 1101 FORMAT(5X,'  THE AVERAGED KINETIC ENERGY IS',F7.2,' KJOULE/MOL')
      WRITE (LP,1301) VINSUM
 1301 FORMAT(5X,'  THE AVERAGED VINTRA IS',F7.4,' KJOULE/MOL')
      WRITE (LP,1201) POTESUM
 1201 FORMAT(5X,'  THE AVERAGED POTENTIAL ENERGY IS',F7.2,' KJOULE/MOL')
      WRITE (LP,1501) DIE
 1501 FORMAT(5X,'  THE STATIC DIELECTRIC CONSTANT IS',F7.2)

 221  FORMAT(1X,E13.7)
      WRITE (LP_CONFIG,221) DENSITY
 222  FORMAT(6(1X,E13.7))
      DO 244 I=1,NP
      WRITE (LP_CONFIG,222) X(I),Y(I),Z(I),X1(I),Y1(I),Z1(I)
 244  CONTINUE

      STOP
      END
C     ------------------- PROGRAM END HERE -------------------

C     --------------------------------------------------------
      SUBROUTINE RANDOM
      IMPLICIT REAL(A-H,O-Z)
      COMMON /CONTROL/ KB,KSTART,SCALER,SEED,DELTA,BON,ANG,PM2
      DATA K,J,M,RM / 5701,3612,566927,566927.0 /
      IX=INT(SEED*RM)
      KR = J*IX + K
      IRAND = MOD(KR,M)
      SEED = (FLOAT(IRAND) + 0.5) / RM
      RETURN
      END
C     ---------------------------------------------------------
     
C     ---------------------------------------------------------
      SUBROUTINE INTVEL
C     ASSIGN INITIAL VELOCITIES TO MOLECULES

      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NPD = 648)
      COMMON /POS/ X(NPD),Y(NPD),Z(NPD),X1(NPD),Y1(NPD),Z1(NPD)
      DOUBLE PRECISION X1,Y1,Z1
      COMMON /FORCE/ FX(NPD),FY(NPD),FZ(NPD),A12,A6
      COMMON /PARM/ RDEL,IGOFR(3,1500),NM,NPO,NPH,NP,PMASS(NPD),
     1              ITYPE(NPD),Q(NPD),VOLUME,CUBE,FROG
      COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2
      COMMON /CONTROL/ KB,KSTART,SCALER,SEED,DELTA,BON,ANG,PM2
      REAL MO,MH

      PX=0.
      PY=0.
      PZ=0.
       
      DO 360 I=1,NP
      IF (MOD(I,3).EQ.1) THEN
      CALL RANDOM
      SEED1 = SEED
      CALL RANDOM
      SEED2 = SEED
      CALL RANDOM
      SEED3 = SEED
      ENDIF
      X1(I) = SEED1 
      Y1(I) = SEED2 
      Z1(I) = SEED3 
      PX = PX + PMASS(I)*X1(I)
      PY = PY + PMASS(I)*Y1(I)
      PZ = PZ + PMASS(I)*Z1(I)
360   CONTINUE

C     SCALE VELOCITIES SO THAT TOTAL MOMEMTUM=0.
      SUMKINE = 0.
      C = 1./FLOAT(NM)/(NPO*MO+NPH*MH)
      DO 370 I=1,NP
      X1(I)=X1(I)-PX*C
      Y1(I)=Y1(I)-PY*C
      Z1(I)=Z1(I)-PZ*C
      SUMKINE = SUMKINE+PMASS(I)*(X1(I)*X1(I)+Y1(I)*Y1(I)+Z1(I)*Z1(I))
370   CONTINUE
      SUMKINE = SUMKINE/2.

C     SCALE VELOCITIES TO DESIRED TEMPERATURE
C     MULTIPLY BY 6, FOR WATER ONLY. JUN.14
      FACTOR=SQRT(6*SCALER/SUMKINE)
      DO 380 I=1,NP
      X1(I)=X1(I)*FACTOR
      Y1(I)=Y1(I)*FACTOR
      Z1(I)=Z1(I)*FACTOR
380   CONTINUE

      RETURN
      END
C     --------------------------------------------------------------

C     --------------------------------------------------------------
      SUBROUTINE INTERACT

      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NPD = 648)
      PARAMETER (ID  = 50000)
      COMMON /EWALD/ ALPHA,ERR(ID),CLI(ID),DERR,A(9,9,9)
      COMMON /POS/ X(NPD),Y(NPD),Z(NPD),X1(NPD),Y1(NPD),Z1(NPD)
      DOUBLE PRECISION X1,Y1,Z1
      COMMON /FORCE/ FX(NPD),FY(NPD),FZ(NPD),A12,A6
      COMMON /PARM/ RDEL,IGOFR(3,1500),NM,NPO,NPH,NP,PMASS(NPD),
     1              ITYPE(NPD),Q(NPD),VOLUME,CUBE,FROG
      COMMON /ENERGY/ TOTE,POTE,KINE,PLJ12,PLJ6,ELEC,VINTRA,VIRIAL
      COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2
      COMMON /CONTROL/ KB,KSTART,SCALER,SEED,DELTA,BON,ANG,PM2
      COMMON /ADDENDUM/ PMX,PMY,PMZ

      REAL KINE,MO,MH

      DIMENSION SINN(NPD,9,9,9),COSN(NPD,9,9,9)
      DIMENSION XR(3,3),YR(3,3),ZR(3,3),DR(3,3)
      DIMENSION CMX(216),CMY(216),CMZ(216)

      QUOTAH = MH/(MO+2*MH)
      QUOTAO = MO/(MO+2*MH)

      BON   = 0.
      ANG   = 0.
      PMX   = 0.
      PMY   = 0.
      PMZ   = 0.      
      ELEC  = 0.
      PLJ6  = 0.
      PLJ12 = 0.
      VINTRA = 0.
      VIRIAL = 0.
      W = 2.*ALPHA/SQRT(PI)
      RCUT = CUBE/2.
      TUCR = -RCUT
      
C     RECIPROCAL SPACE SUMMATION

      DO 805 I = 1,NP
C     THE INERT V2 TERM
      ELEC  = ELEC - Q(I)*Q(I)*W/2.
      FX(I) = 0.
      FY(I) = 0.
      FZ(I) = 0.
805   CONTINUE

      DO 1000 I=1,NP

      SX = X(I)*FROG
      SY = Y(I)*FROG
      SZ = Z(I)*FROG

      SINX = SIN(SX)
      COSX = COS(SX)
      SINY = SIN(SY)
      COSY = COS(SY)
      SINZ = SIN(SZ)
      COSZ = COS(SZ)
 
      SINN(I,1,1,1) = SIN(-4.*SX-4.*SY-4.*SZ)
      COSN(I,1,1,1) = COS(-4.*SX-4.*SY-4.*SZ)

      DO 1100 J=2,9
      SINN(I,J,1,1) = SINN(I,J-1,1,1)*COSX+COSN(I,J-1,1,1)*SINX
      COSN(I,J,1,1) = COSN(I,J-1,1,1)*COSX-SINN(I,J-1,1,1)*SINX
1100  CONTINUE

      DO 1200 J=1,9
      DO 1200 K=2,9
      SINN(I,J,K,1) = SINN(I,J,K-1,1)*COSY+COSN(I,J,K-1,1)*SINY
      COSN(I,J,K,1) = COSN(I,J,K-1,1)*COSY-SINN(I,J,K-1,1)*SINY
1200  CONTINUE

      DO 1300 J=1,9
      DO 1300 K=1,9
      DO 1300 L=2,9
      SINN(I,J,K,L) = SINN(I,J,K,L-1)*COSZ+COSN(I,J,K,L-1)*SINZ
      COSN(I,J,K,L) = COSN(I,J,K,L-1)*COSZ-SINN(I,J,K,L-1)*SINZ
1300  CONTINUE
1000  CONTINUE

      DO 1400 J=1,9
      DO 1400 K=1,9
      DO 1400 L=1,9

      QX = (J-5)*FROG
      QY = (K-5)*FROG 
      QZ = (L-5)*FROG
      Q2 = QX*QX+QY*QY+QZ*QZ
      IF (Q2.EQ.0.) GOTO 1400
      C1 = EXP(-Q2/ALPHA/ALPHA/4.)
      A(J,K,L) = 4.*PI*C1/Q2/VOLUME
      QK = 2./Q2+1./2./ALPHA/ALPHA

      COSSUM = 0.
      SINSUM = 0.
      
      DO 1500 I=1,NP
      COSSUM = COSSUM + Q(I)*COSN(I,J,K,L)
      SINSUM = SINSUM + Q(I)*SINN(I,J,K,L)
1500  CONTINUE

      DO 1600 I=1,NP
      CC = A(J,K,L)*Q(I)*(COSSUM*SINN(I,J,K,L)-SINSUM*COSN(I,J,K,L))
200   FX(I) = FX(I)+QX*CC
      FY(I) = FY(I)+QY*CC
      FZ(I) = FZ(I)+QZ*CC
1600  CONTINUE
      
      ELEC = ELEC+A(J,K,L)*(COSSUM*COSSUM+SINSUM*SINSUM)/2.
1400  CONTINUE

      DO 522 I=1,NP
      IMI = (I-1)/3
      II  = I-IMI*3

      DO 5520 J=I+1,NP
      IMJ = (J-1)/3
      IJ  = J-IMJ*3

      RX = X(I)-X(J)
      RY = Y(I)-Y(J)
      RZ = Z(I)-Z(J)
 
      IF (RX.GT.RCUT) RX = RX - CUBE
      IF (RX.LT.TUCR) RX = RX + CUBE
      IF (RY.GT.RCUT) RY = RY - CUBE
      IF (RY.LT.TUCR) RY = RY + CUBE
      IF (RZ.GT.RCUT) RZ = RZ - CUBE
      IF (RZ.LT.TUCR) RZ = RZ + CUBE

      R2 = RX*RX+RY*RY+RZ*RZ
      R1 = SQRT(R2)
      R3 = R2*R1
      IJT = ITYPE(I)+ITYPE(J)+1

      IF (IMI.NE.IMJ) THEN
      IR = R1/RDEL+1
      IF (IR.LE.500) IGOFR(IJT,IR) = IGOFR(IJT,IR)+1
      ENDIF

C     CALCULATE EWALD FACTOR IN REAL SPACE SUMMATION
      PRODUCT = ALPHA*R1
      I0 = PRODUCT/DERR
      XX = I0*DERR
      ARGIN = PRODUCT-XX

      IF (I0.GE.ID) GOTO 5520
      GACTOR = ERR(I0)-CLI(I0)*ARGIN*(1-XX*ARGIN)

C     IF THEY ARE THE SAME MOLECULE
      IF (IMI.EQ.IMJ) THEN
       GACTOR = GACTOR-1.
       XR(II,IJ) = RX
       YR(II,IJ) = RY
       ZR(II,IJ) = RZ     
       DR(II,IJ) = R1
      ENDIF

      ARGIN = PRODUCT*PRODUCT-XX*XX
      E1 = 1.-ARGIN*(1-ARGIN*0.5)
      FACTOR = GACTOR + E1*PRODUCT*CLI(I0)

      IF ( IJT .EQ. 1 ) THEN
       R6  = R3*R3
       R12 = R6*R6
       E1  = A12/R12
       E2  = A6/R6
       PLJ12 = PLJ12 + E1
       PLJ6  = PLJ6  + E2
       ELEC  = ELEC  + GACTOR*QOO/R1
       FF   = 12.*E1 + 6.*E2
       FF   = FF/R2 + QOO*FACTOR/R3 

      ELSEIF ( IJT .EQ. 2 ) THEN
       IF (R1.LT.0.3) GOTO 5520
       ELEC = ELEC + GACTOR*QOH/R1
       FF = FACTOR*QOH/R3

      ELSEIF ( IJT .EQ. 3 ) THEN
       ELEC = ELEC + GACTOR*QHH/R1
       FF = FACTOR*QHH/R3

      ENDIF

      FFX= FF*RX
      FFY= FF*RY
      FFZ= FF*RZ
      FX(I) = FX(I) + FFX
      FY(I) = FY(I) + FFY
      FZ(I) = FZ(I) + FFZ 
      FX(J) = FX(J) - FFX
      FY(J) = FY(J) - FFY
      FZ(J) = FZ(J) - FFZ

5520  CONTINUE
      
      IF (II.EQ.3) THEN
       PMX = PMX+XR(1,2)+XR(1,3)
       PMY = PMY+YR(1,2)+YR(1,3)
       PMZ = PMZ+ZR(1,2)+ZR(1,3)
       BON = BON+DR(1,2)+DR(1,3)
       ANG = ANG+ACOS((DR(1,2)**2+DR(1,3)**2-DR(2,3)**2)/
     1       2./DR(1,2)/DR(1,3))
       DR1 = DR(1,2)-BOND
       DR2 = DR(1,3)-BOND
       DR3 = DR(2,3)-HHDIST

       DR11 = DR1*DR1
       DR22 = DR2*DR2
       DR33 = DR3*DR3
       
       DR12 = DR1*DR2
       DR13 = DR1*DR3
       DR23 = DR2*DR3

       DR111 = DR11*DR1
       DR222 = DR22*DR2
       DR333 = DR33*DR3

       DR112 = DR11*DR2
       DR113 = DR11*DR3
       DR122 = DR22*DR1
       DR123 = DR12*DR3
       DR133 = DR33*DR1
       DR223 = DR22*DR3
       DR233 = DR33*DR2

       F1 = 3.959608852*DR1+0.2514577891*DR2-0.5348000472*DR3
       G1 = -13.52143988*DR11-0.768247786*DR12
     1      +1.233829205*DR13-0.3841238925*DR22
     2      +0.944571334*DR23-0.6852521196*DR33

       F2 = 3.959608852*DR2+0.2514577891*DR1-0.5348000472*DR3
       G2 = -13.52143988*DR22-0.768247786*DR12
     1      +1.233829205*DR23-0.3841238925*DR11
     2      +0.944571334*DR13-0.6852521196*DR33

       F3 = -0.5348000472*(DR1+DR2)+0.8804540624*DR3
       G3 = 0.6169146026*(DR11+DR22)+0.6267843126*DR33+
     1      0.944571334*DR12-1.370504239*(DR13+DR23)


       H1 =  24.83326248*DR111-1.017614014*DR333
     1      +6.876821916*DR113-6.767337447*DR112
     2      -1.100316147*DR133-5.500673460*DR122
     3      +5.872225848*DR123-2.255779147*DR222
     4      +2.936112915*DR223+1.223369784*DR233


       H2 = -1.017614014*DR333+24.83326249*DR222
     1      +6.876821910*DR223-2.255779149*DR111
     2      -5.500673460*DR112+2.936112924*DR113
     3      -6.767337441*DR122+5.872225830*DR123
     4      +1.223369784*DR133-1.100316147*DR233

       H3 = -3.052842042*(DR133+DR233)+2.292273970*(DR111+DR222)
     1      -1.100316147*(DR113+DR223)+2.936112924*(DR112+DR122)
     2      +2.446739568*DR123+1.274486629*DR333


C      SET CUTOFF SO THAT WON'T COLLAPSE DUE TO 3RD AND 4TH ORDER TERMS.

       IF (DR1.GT.0.07.OR.DR1.LT.-0.07.OR.DR2.GT.0.07.OR.
     1    DR2.LT.-0.07.OR.DR3.GT.0.14.OR.DR3.LT.-0.14) THEN
       G1 = 0.
       H1 = 0.
       G2 = 0.
       H2 = 0.
       G3 = 0.
       H3 = 0.
       ENDIF
  
       VINTRA = VINTRA+(F1/2.+G1/3.+H1/4.)*DR1+(F2/2.+G2/3.+H2/4.)*DR2
     1          +(F3/2.+G3/3.+H3/4.)*DR3
       F1 = F1+G1+H1
       F2 = F2+G2+H2
       F3 = F3+G3+H3

       FX(I-2) = FX(I-2)-F1*XR(1,2)/DR(1,2)-F2*XR(1,3)/DR(1,3) 
       FY(I-2) = FY(I-2)-F1*YR(1,2)/DR(1,2)-F2*YR(1,3)/DR(1,3) 
       FZ(I-2) = FZ(I-2)-F1*ZR(1,2)/DR(1,2)-F2*ZR(1,3)/DR(1,3) 

       FX(I-1) = FX(I-1)+F1*XR(1,2)/DR(1,2)-F3*XR(2,3)/DR(2,3) 
       FY(I-1) = FY(I-1)+F1*YR(1,2)/DR(1,2)-F3*YR(2,3)/DR(2,3) 
       FZ(I-1) = FZ(I-1)+F1*ZR(1,2)/DR(1,2)-F3*ZR(2,3)/DR(2,3) 

       FX(I) = FX(I)+F2*XR(1,3)/DR(1,3)+F3*XR(2,3)/DR(2,3) 
       FY(I) = FY(I)+F2*YR(1,3)/DR(1,3)+F3*YR(2,3)/DR(2,3) 
       FZ(I) = FZ(I)+F2*ZR(1,3)/DR(1,3)+F3*ZR(2,3)/DR(2,3) 

       VIRIAL = VIRIAL-F1*DR(1,2)-F2*DR(1,3)-F3*DR(2,3)
      ENDIF
522   CONTINUE

C     CALCULATE TOTAL POTENTIAL ENERGIES
      BON  = BON/2./NM
      ANG  = ANG/NM
      PM2  = PMX*PMX+PMY*PMY+PMZ*PMZ

C     THE CUTOFF REMEDY
      PLJ6 = PLJ6+A6*NM/VOLUME*PI*2./3./RCUT/RCUT/RCUT*NM
      POTE = PLJ12+PLJ6+ELEC+VINTRA
      VIRIAL = VIRIAL+12.*PLJ12+6.*PLJ6+ELEC

      RETURN
      END
C     ----------------------------------------------------------------



C     ----------------------------------------------------------------     
      SUBROUTINE DIAMOND

      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NPD = 648)
      COMMON /POS/ X(NPD),Y(NPD),Z(NPD),X1(NPD),Y1(NPD),Z1(NPD)
      DOUBLE PRECISION X1,Y1,Z1
      COMMON /PARM/ RDEL,IGOFR(3,1500),NM,NPO,NPH,NP,PMASS(NPD),
     1              ITYPE(NPD),Q(NPD),VOLUME,CUBE,FROG
      COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2

      DIMENSION SX(24),SY(24),SZ(24)
      DIMENSION VECTOR(4,3),MTYPE(16)
      DATA VECTOR/ 1,1,-1,-1,
     Y             1,-1,-1,1,
     Z             1,-1,1,-1/
      DATA MTYPE /4,2,4,2,3,1,3,1,-2,-1,-2,-1,-3,-4,-3,-4/

      N = NINT((NM/8)**(1./3.))
      A = CUBE/N

      SX(1) = 0.
      SY(1) = 0.
      SZ(1) = 0.

      SX(7) = 1./2.
      SY(7) = 1./2.
      SZ(7) = 0. 

      SX(10) = 0.
      SY(10) = 1./2.
      SZ(10) = 1./2. 

      SX(4) = 1./2.
      SY(4) = 0.
      SZ(4) = 1./2.

      DO 100 I = 13,22,3
      SX(I) = SX(I-12)+0.25
      SY(I) = SY(I-12)+0.25
      SZ(I) = SZ(I-12)+0.25
100   CONTINUE

      IO = -2
      IB = 1
      DO 500 I = 1,24

      IF (MOD(I,3).EQ.1) THEN
      SX(I) = SX(I)*A
      SY(I) = SY(I)*A
      SZ(I) = SZ(I)*A
      IO = IO+3

      ELSE
      IC = MTYPE(IB)
      IB = IB+1
      IF (IC.LT.0) THEN
      IC = -IC
      SX(I) = SX(IO)-VECTOR(IC,1)/SQRT(3.)*BOND
      SY(I) = SY(IO)-VECTOR(IC,2)/SQRT(3.)*BOND
      SZ(I) = SZ(IO)-VECTOR(IC,3)/SQRT(3.)*BOND
      ELSE
      SX(I) = SX(IO)+VECTOR(IC,1)/SQRT(3.)*BOND
      SY(I) = SY(IO)+VECTOR(IC,2)/SQRT(3.)*BOND
      SZ(I) = SZ(IO)+VECTOR(IC,3)/SQRT(3.)*BOND
      ENDIF

      ENDIF
 500  CONTINUE

      M  = 0
      DO 200 I = 0,N-1
      DO 200 J = 0,N-1
      DO 200 K = 0,N-1
      DO 300 L = 1,24
      X(L+M) = SX(L)+I*A
      Y(L+M) = SY(L)+J*A
      Z(L+M) = SZ(L)+K*A
300   CONTINUE
      M = M + 24
200   CONTINUE

      DO 700 I= 1,NP
      IF(X(I).LT.0)     X(I) = X(I)+CUBE
      IF(X(I).GT.CUBE)  X(I) = X(I)-CUBE
      IF(Y(I).LT.0)     Y(I) = Y(I)+CUBE
      IF(Y(I).GT.CUBE)  Y(I) = Y(I)-CUBE
      IF(Z(I).LT.0)     Z(I) = Z(I)+CUBE
      IF(Z(I).GT.CUBE)  Z(I) = Z(I)-CUBE
700   CONTINUE
      
      RETURN
      END
C     ---------------------------------------------------