      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                                          
      DIMENSION SUM(3)                                                  
                                                                        
      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                                                                       
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/                     
CX00001 CALL MAINIX (1)                                                 
CY00001 CALL MAINIY (1)                                                 
CX60000 CALL CNTINX ()                                                  
CY60000 CALL CNTINY ()                                                  
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 INSTRUMENTED VERSION'/
     *       1X,'$Revision: 1.A.1 $ $Author: psharma $'/)     
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)                                                    
      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                                                       
CX10501 CALL DOINX( 501 )                                               
CY10501 CALL DOINY( 501 )                                               
      DO 100 I=1,NVAR                                                   
      IX(I)=II                                                          
  100 II=II+NATMO                                                       
CX20502 CALL DOOUTX( 502 )                                              
CY20502 CALL DOOUTY( 502 )                                              
      IXF=IX(NVAR-2)                                                    
      IYF=IX(NVAR-1)                                                    
      IZF=IX(NVAR)                                                      
      IMY=NMOL+1                                                        
      IMZ=NMOL+IMY                                                      
CX10503 CALL DOINX( 503 )                                               
CY10503 CALL DOINY( 503 )                                               
      DO 101 I=1,NDVAR                                                  
  101 VAR(I)=0.0D0                                                      
CX20504 CALL DOOUTX( 504 )                                              
CY20504 CALL DOOUTY( 504 )                                              
C                                                                       
C.....CHECK WHAT TO BE DONE (NEW RUN, CONTINUATION, ETC.)               
C                                                                       
      IF(NSAVE.GT.0) GOTO 200                                           
C                                                                       
C.....INITIALIZATION OR EQUILIBRIATION RUN                              
C                                                                       
      IF(NSAVE.EQ.0) GOTO 200                                           
      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,TVIR,TKIN,TTMV,ELPST)     
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                                                          
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.EQ.TSTEP) GO TO 300                                      
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                                                       
CX10505 CALL DOINX( 505 )                                               
CY10505 CALL DOINY( 505 )                                               
      DO 210 I=7,NVAR,3                                                 
      FAC=FAC*VF                                                        
      DO 211 J=1,NATMO3                                                 
      VAR(IJ)=VAR(IJ)*FAC                                               
  211 IJ=IJ+1                                                           
  210 CONTINUE                                                          
CX20506 CALL DOOUTX( 506 )                                              
CY20506 CALL DOOUTY( 506 )                                              
C                                                                       
C.....START MOLECULAR DYNAMIC LOOP                                      
C                                                                       
  300 IF(NFMC.LT.0) THEN                                                
       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,D)       
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)                    
CX70000 CALL CNTOTX ()                                                  
CY70000 CALL CNTOTY ()                                                  
CX50000 CALL MAINOX ()                                                  
CY50000 CALL MAINOY ()                                                  
      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')     
CX70000 CALL CNTOTX ()                                                  
CY70000 CALL CNTOTY ()                                                  
CX50000 CALL MAINOX ()                                                  
CY50000 CALL MAINOY ()                                                  
C                                                                       
      END                                                               
C***********************************************************************
      SUBROUTINE MDMAIN(NFSV,NFRST,NSTEP,NRST,NPRINT,NSAVE,LKT,NORD1,D) 
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)                                                  
CX30003 CALL SBINX (3)                                                  
CY30003 CALL SBINY (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)                                                    
CX10507 CALL DOINX( 507 )                                               
CY10507 CALL DOINY( 507 )                                               
      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     
      RI = I                                                            
      RNSAVE = NSAVE                                                    
      CALL SAVEDT(NFRST,VAR,NXYZ,TVIR,RI,RNSAVE,D,D,D)                  
      REWIND NFRST                                                      
      WRITE (6,1005) IRST+1,NFRST                                       
      ENDIF                                                             
 2000 CONTINUE                                                          
CX20508 CALL DOOUTX( 508 )                                              
CY20508 CALL DOOUTY( 508 )                                              
C                                                                       
    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)                                       
CX40004 CALL SBOUTX (4)                                                 
CY40004 CALL SBOUTY (4)                                                 
      RETURN                                                            
      END                                                               
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            
CX30005 CALL SBINX (5)                                                  
CY30005 CALL SBINY (5)                                                  
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                                                      
CX40006 CALL SBOUTX (6)                                                 
CY40006 CALL SBOUTY (6)                                                 
      RETURN                                                            
      END                                                               
C***********************************************************************
      SUBROUTINE KINETI(VAR,NMOL,SUM,HMAS,OMAS)                         
C***********************************************************************
C                                                                       
C.....THIS ROUTINE EVALUATES KINETIC ENERGY                             
C                                                                       
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      DIMENSION VAR(1),SUM(1)                                           
CX30007 CALL SBINX (7)                                                  
CY30007 CALL SBINY (7)                                                  
C                                                                       
      JJ=1                                                              
CX10509 CALL DOINX( 509 )                                               
CY10509 CALL DOINY( 509 )                                               
      DO 100 K=1,3                                                      
      S=0.0D0                                                           
      DO 110 J=1,NMOL                                                   
      S=S+(VAR(JJ)*VAR(JJ)+VAR(JJ+2)*VAR(JJ+2))*HMAS                    
     *   +VAR(JJ+1)*VAR(JJ+1)*OMAS                                      
  110 JJ=JJ+3                                                           
      SUM(K)=S                                                          
  100 CONTINUE                                                          
CX20510 CALL DOOUTX( 510 )                                              
CY20510 CALL DOOUTY( 510 )                                              
CX40008 CALL SBOUTX (8)                                                 
CY40008 CALL SBOUTY (8)                                                 
      RETURN                                                            
      END                                                               
C***********************************************************************
      SUBROUTINE NRMLKT(VAR,SUM,FKIN,NATMO,TVIR,TKIN,TTMV,ELPST)        
C***********************************************************************
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      DIMENSION VAR(1),SUM(1)                                           
CX30009 CALL SBINX (9)                                                  
CY30009 CALL SBINY (9)                                                  
C                                                                       
C.....NORMALIZE EACH COMPOENT OF KINETIC ENERGY TO KT/2                 
C                                                                       
      KMAX=0                                                            
CX10511 CALL DOINX( 511 )                                               
CY10511 CALL DOINY( 511 )                                               
      DO 300 J=1,3                                                      
      XT=SQRT(0.5D0/(SUM(J)*FKIN))                                      
      KMIN=KMAX+1                                                       
      KMAX=KMAX+NATMO                                                   
      DO 300 K=KMIN,KMAX                                                
      VAR(K)=VAR(K)*XT                                                  
  300 CONTINUE                                                          
CX20512 CALL DOOUTX( 512 )                                              
CY20512 CALL DOOUTY( 512 )                                              
      TVIR=0.0D0                                                        
      TKIN=0.0D0                                                        
      TTMV=0.0D0                                                        
      ELPST=0.0D0                                                       
CX40010 CALL SBOUTX (10)                                                
CY40010 CALL SBOUTY (10)                                                
      RETURN                                                            
      END                                                               
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)                                                   
CX30011 CALL SBINX (11)                                                 
CY30011 CALL SBINY (11)                                                 
C                                                                       
      WRITE (NF) TENG,POTA,POTR,ATEMP,NV,V,POTRF,VIR                    
CX40012 CALL SBOUTX (12)                                                
CY40012 CALL SBOUTY (12)                                                
      RETURN                                                            
      END                                                               
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)                                                    
CX30013 CALL SBINX (13)                                                 
CY30013 CALL SBINY (13)                                                 
CX10513 CALL DOINX( 513 )                                               
CY10513 CALL DOINY( 513 )                                               
C                                                                       
      DO 100 I=2,NV,NS                                                  
      IF(V(I).LE.BOXL) GO TO 10                                         
      V(I-1)=V(I-1)-BOXL                                                
      V(I  )=V(I  )-BOXL                                                
      V(I+1)=V(I+1)-BOXL                                                
      GO TO 100                                                         
   10 IF(V(I).GE.0.0D0) GO TO 100                                       
      V(I-1)=V(I-1)+BOXL                                                
      V(I  )=V(I  )+BOXL                                                
      V(I+1)=V(I+1)+BOXL                                                
  100 CONTINUE                                                          
CX20514 CALL DOOUTX( 514 )                                              
CY20514 CALL DOOUTY( 514 )                                              
CX40014 CALL SBOUTX (14)                                                
CY40014 CALL SBOUTY (14)                                                
      RETURN                                                            
      END                                                               
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)                                                  
CX30015 CALL SBINX (15)                                                 
CY30015 CALL SBINY (15)                                                 
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                                                      
CX10515 CALL DOINX( 515 )                                               
CY10515 CALL DOINY( 515 )                                               
      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                                                          
CX20516 CALL DOOUTX( 516 )                                              
CY20516 CALL DOOUTY( 516 )                                              
C                                                                       
C.....PREDICTOR-CORRECTOR CONSTANTS FOR 2ND ORDER DIFFERENTIAL EQUATION 
C                                                                       
      PCC(3)=1.0D0                                                      
      N1=N-1                                                            
      GO TO (1,2,3,4,5,6,7), N1                                         
    1 WRITE (6,10)                                                      
CX40016 CALL SBOUTX (16)                                                
CY40016 CALL SBOUTY (16)                                                
      RETURN                                                            
    2 WRITE (6,10)                                                      
CX40016 CALL SBOUTX (16)                                                
CY40016 CALL SBOUTY (16)                                                
      RETURN                                                            
    3 PCC(1)=1.0D0/6.0D0                                                
      PCC(2)=5.0D0/6.0D0                                                
      PCC(4)=1.0D0/3.0D0                                                
CX40016 CALL SBOUTX (16)                                                
CY40016 CALL SBOUTY (16)                                                
      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                                               
CX40016 CALL SBOUTX (16)                                                
CY40016 CALL SBOUTY (16)                                                
      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                                               
CX40016 CALL SBOUTX (16)                                                
CY40016 CALL SBOUTY (16)                                                
      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                                              
CX40016 CALL SBOUTX (16)                                                
CY40016 CALL SBOUTY (16)                                                
      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                                             
CX40016 CALL SBOUTX (16)                                                
CY40016 CALL SBOUTY (16)                                                
      RETURN                                                            
   10 FORMAT('0***** ERROR: THE ORDER HAS TO BE GREATER THAN 2 ****')   
CX40016 CALL SBOUTX (16)                                                
CY40016 CALL SBOUTY (16)                                                
      END                                                               
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                                             
      DATA XMIN,YMIN,ZMIN/3*0.0D0/                                      
C                                                                       
      RANDOM(L)=2.0D0*RAND(0)-1.0D0                                     
CX30017 CALL SBINX (17)                                                 
CY30017 CALL SBINY (17)                                                 
      XMAS(2)=SQRT(OMAS*HMAS)                                           
      XMAS(1)=HMAS                                                      
      XMAS(3)=HMAS                                                      
C                                                                       
C.....ASSIGN POSITIONS                                                  
C                                                                       
      IF(NFMC.GT.0) GO TO 10                                            
      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                                                        
CX10517 CALL DOINX( 517 )                                               
CY10517 CALL DOINY( 517 )                                               
      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                                                    
      DO 1111 L=1,NATOMS                                                
      IJK=IJK+1                                                         
      XC(IJK)=XT(L)                                                     
      YC(IJK)=YT(L)                                                     
      ZC(IJK)=Z                                                         
 1111 CONTINUE                                                          
      Z=Z+XS                                                            
 1110 CONTINUE                                                          
      YT(2)=YT(2)+XS                                                    
 1100 CONTINUE                                                          
      XT(2)=XT(2)+XS                                                    
 1000 CONTINUE                                                          
CX20518 CALL DOOUTX( 518 )                                              
CY20518 CALL DOOUTY( 518 )                                              
      GO TO 12                                                          
   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                                                       
CX10519 CALL DOINX( 519 )                                               
CY10519 CALL DOINY( 519 )                                               
      DO 300 I=1,NATMO                                                  
      XMIN=MIN(XMIN,XC(I))                                              
      YMIN=MIN(YMIN,YC(I))                                              
      ZMIN=MIN(ZMIN,ZC(I))                                              
  300 CONTINUE                                                          
CX20520 CALL DOOUTX( 520 )                                              
CY20520 CALL DOOUTY( 520 )                                              
CX10521 CALL DOINX( 521 )                                               
CY10521 CALL DOINY( 521 )                                               
      DO 400 I=1,NATMO                                                  
      XC(I)=XC(I)-XMIN                                                  
      YC(I)=YC(I)-YMIN                                                  
      ZC(I)=ZC(I)-ZMIN                                                  
  400 CONTINUE                                                          
CX20522 CALL DOOUTX( 522 )                                              
CY20522 CALL DOOUTY( 522 )                                              
C                                                                       
C.....ASSIGN RANDOM MOMENTA                                             
C                                                                       
   12 SUX=RAND(1)                                                       
      SUMX=0.0D0                                                        
      SUMY=0.0D0                                                        
      SUMZ=0.0D0                                                        
      IJ=1                                                              
CX10523 CALL DOINX( 523 )                                               
CY10523 CALL DOINY( 523 )                                               
      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                                                          
CX20524 CALL DOOUTX( 524 )                                              
CY20524 CALL DOOUTY( 524 )                                              
C                                                                       
C.....FIND AVERAGE MOMENTA PER ATOM                                     
C                                                                       
      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                                                               
CX10525 CALL DOINX( 525 )                                               
CY10525 CALL DOINY( 525 )                                               
      DO 100 I=1,NMOL                                                   
      SUX=SUX+((VX(J)-SUMX)**2+(VX(J+2)-SUMX)**2)/HMAS                  
     *       +(VX(J+1)-SUMX)**2/OMAS                                    
      SUY=SUY+((VY(J)-SUMY)**2+(VY(J+2)-SUMY)**2)/HMAS                  
     *       +(VY(J+1)-SUMY)**2/OMAS                                    
      SUZ=SUZ+((VZ(J)-SUMZ)**2+(VZ(J+2)-SUMZ)**2)/HMAS                  
     *       +(VZ(J+1)-SUMZ)**2/OMAS                                    
      J=J+NATOMS                                                        
  100 CONTINUE                                                          
CX20526 CALL DOOUTX( 526 )                                              
CY20526 CALL DOOUTY( 526 )                                              
      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                                                               
CX10527 CALL DOINX( 527 )                                               
CY10527 CALL DOINY( 527 )                                               
      DO 200 I=1,NMOL                                                   
      DO 200 L=1,NATOMS                                                 
      VX(J)=(VX(J)-SUMX)*SUX/XMAS(L)                                    
      VY(J)=(VY(J)-SUMY)*SUY/XMAS(L)                                    
      VZ(J)=(VZ(J)-SUMZ)*SUZ/XMAS(L)                                    
      J=J+1                                                             
  200 CONTINUE                                                          
CX20528 CALL DOOUTX( 528 )                                              
CY20528 CALL DOOUTY( 528 )                                              
CX40018 CALL SBOUTX (18)                                                
CY40018 CALL SBOUTY (18)                                                
      RETURN                                                            
      END                                                               
C***********************************************************************
      FUNCTION RAND(N)                                                  
C***********************************************************************
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      COMMON /RANDNO/ R3(127),R1,I2                                     
      SAVE S, T, RMC                                                    
      SAVE IW                                                           
      DATA S, T, RMC /0.0D0, 1.0D0, 1.0D0/                              
      DATA IW /-1/                                                      
CX30019 CALL SBINX (19)                                                 
CY30019 CALL SBINY (19)                                                 
C                                                                       
      IF(R1.LT.1.0D0.AND.N.EQ.0) GO TO 60                               
      IF(IW.GT.0) GO TO 30                                              
   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)                                                     
CX10529 CALL DOINX( 529 )                                               
CY10529 CALL DOINY( 529 )                                               
      DO 20 I=1,IC                                                      
   20 RMC=0.5D0*RMC                                                     
CX20530 CALL DOOUTX( 530 )                                              
CY20530 CALL DOOUTY( 530 )                                              
      RM=0.015625D0*0.015625D0                                          
   30 I2=127                                                            
      IR=MOD(ABS(N),8190)+1                                             
   40 R1=0.0D0                                                          
CX10531 CALL DOINX( 531 )                                               
CY10531 CALL DOINY( 531 )                                               
      DO 50 I=1,IKT                                                     
      IR=MOD(17*IR,8191)                                                
   50 R1=(R1+REAL(IR/2))*RM                                             
CX20532 CALL DOOUTX( 532 )                                              
CY20532 CALL DOOUTY( 532 )                                              
      IR=MOD(17*IR,8191)                                                
      R1=(R1+REAL(IR/ID))*RMC                                           
      R3(I2)=R1                                                         
      I2=I2-1                                                           
      IF(I2.GT.0) GO TO 40                                              
   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                                                           
CX40020 CALL SBOUTX (20)                                                
CY40020 CALL SBOUTY (20)                                                
      RETURN                                                            
      END                                                               
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,NMOL,NORDER,  
     *               NATMO,NATMO3,NMOL1                                 
      DIMENSION V(1),C(1)                                               
CX30021 CALL SBINX (21)                                                 
CY30021 CALL SBINY (21)                                                 
C                                                                       
      IK=1                                                              
      JIZ=2                                                             
C.....LOOP OVER F(X), F'(X), F''(X), .....                              
CX10533 CALL DOINX( 533 )                                               
CY10533 CALL DOINY( 533 )                                               
      DO 1000 I=1,NORDER                                                
C.....LOOP OVER X, Y, AND Z OF NATOMS*NMOL                              
      DO 1100 K=1,NATMO3                                                
      JI=JIZ                                                            
      IKL=IK+NATMO3                                                     
C.....SUM OVER TAYLOR SERIES                                            
      S=0.0D0                                                           
      DO 1110 L=I,NORDER                                                
      S=S+C(JI)*V(IKL)                                                  
      IKL=IKL+NATMO3                                                    
      JI=JI+1                                                           
 1110 CONTINUE                                                          
      V(IK)=V(IK)+S                                                     
      IK=IK+1                                                           
 1100 CONTINUE                                                          
      JIZ=JIZ+NOR1                                                      
 1000 CONTINUE                                                          
CX20534 CALL DOOUTX( 534 )                                              
CY20534 CALL DOOUTY( 534 )                                              
CX40022 CALL SBOUTX (22)                                                
CY40022 CALL SBOUTY (22)                                                
      RETURN                                                            
      END                                                               
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)                               
CX30023 CALL SBINX (23)                                                 
CY30023 CALL SBINY (23)                                                 
CX10535 CALL DOINX( 535 )                                               
CY10535 CALL DOINY( 535 )                                               
C                                                                       
      DO 1000 I=1,NT                                                    
      Y=FR(I)-SD(I)                                                     
      JJ=I                                                              
      DO 1100 J=1,NOR1                                                  
      VAR(JJ)=VAR(JJ)+PCC(J)*Y                                          
      JJ=JJ+NT                                                          
 1100 CONTINUE                                                          
 1000 CONTINUE                                                          
CX20536 CALL DOOUTX( 536 )                                              
CY20536 CALL DOOUTY( 536 )                                              
CX40024 CALL SBOUTX (24)                                                
CY40024 CALL SBOUTY (24)                                                
      RETURN                                                            
      END                                                               
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,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 /FRCNST/ FC11,FC12,FC13,FC33, FC111,FC333,FC112,FC113,     
     *                FC123,FC133, FC1111,FC3333,FC1112,                
     *                FC1122,FC1113,FC1123,FC1133,FC1233,FC1333         
      DIMENSION VR1(3),VR2(3),DR11(3),DR23(3),DT1(3),DT3(3)             
      DIMENSION V(1),FR(1),VM(1)                                        
CX30025 CALL SBINX (25)                                                 
CY30025 CALL SBINY (25)                                                 
C                                                                       
      ICZ=1                                                             
CX10537 CALL DOINX( 537 )                                               
CY10537 CALL DOINY( 537 )                                               
      DO 1000 I=1,NMOL                                                  
      IC=ICZ                                                            
      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                           
      IC=ICZ                                                            
      DO 1300 J=1,3                                                     
      T1=F1*DR11(J)+F2*DT1(J)                                           
      FR(IC  )=T1                                                       
      T2=F3*DR23(J)+F2*DT3(J)                                           
      FR(IC+2)=T2                                                       
      FR(IC+1)=-(T1+T2)                                                 
      IC=IC+NATMO                                                       
 1300 CONTINUE                                                          
      ICZ=ICZ+NATOMS                                                    
 1000 CONTINUE                                                          
CX20538 CALL DOOUTX( 538 )                                              
CY20538 CALL DOOUTY( 538 )                                              
      VIR=0.0D0                                                         
CX10539 CALL DOINX( 539 )                                               
CY10539 CALL DOINY( 539 )                                               
      DO 2000 I=1,NATMO3                                                
 2000 VIR=VIR+V(I)*FR(I)                                                
CX20540 CALL DOOUTX( 540 )                                              
CY20540 CALL DOOUTY( 540 )                                              
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
      RETURN                                                            
      END                                                               
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)                                
      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)              
      DIMENSION X(1),Y(1),Z(1),FX(1),FY(1),FZ(1),XM(1),YM(1),ZM(1)      
      DIMENSION XL(14),YL(14),ZL(14),RL(14),RS(14),FF(14),GG(14)        
CX30027 CALL SBINX (27)                                                 
CY30027 CALL SBINY (27)                                                 
C                                                                       
      IW1=1                                                             
      IWO=2                                                             
      IW2=3                                                             
CX10541 CALL DOINX( 541 )                                               
CY10541 CALL DOINY( 541 )                                               
      DO 1000 I=1,NMOL1                                                 
      JW1=IW1                                                           
      JWO=IWO                                                           
      JW2=IW2                                                           
      DO 1100 J=I+1,NMOL                                                
      JW1=JW1+NATOMS                                                    
      JWO=JWO+NATOMS                                                    
      JW2=JW2+NATOMS                                                    
      CALL CSHIFT(X(IW1),X(JW1),XM(I),XM(J),XL,BOXH,BOXL)               
      CALL CSHIFT(Y(IW1),Y(JW1),YM(I),YM(J),YL,BOXH,BOXL)               
      CALL CSHIFT(Z(IW1),Z(JW1),ZM(I),ZM(J),ZL,BOXH,BOXL)               
      KC=0                                                              
      DO 1110 K=1,9                                                     
      RS(K)=XL(K)*XL(K)+YL(K)*YL(K)+ZL(K)*ZL(K)                         
 1110 IF(RS(K).GT.CUT2) KC=KC+1                                         
      IF(KC.EQ.9) GO TO 1100                                            
      DO 1120 K=1,14                                                    
      FF(K)=0.0D0                                                       
 1120 CONTINUE                                                          
      IF(RS(1).GE.CUT2) GO TO 10                                        
      FF(1)=QQ4/(RS(1)*SQRT(RS(1)))+REF4                                
      VIR=VIR+FF(1)*RS(1)                                               
   10 CONTINUE                                                          
      DO      1130      K=2,5                                           
      IF(RS(K).GE.CUT2) GO TO 11                                        
      FF(K)=-QQ2/(RS(K)*SQRT(RS(K)))-REF2                               
      VIR=VIR+FF(K)*RS(K)                                               
   11 IF(RS(K+4).GT.CUT2) GO TO 1130                                    
      RL(K+4)=SQRT(RS(K+4))                                             
      FF(K+4)=QQ/(RS(K+4)*RL(K+4))+REF1                                 
      VIR=VIR+FF(K+4)*RS(K+4)                                           
 1130 CONTINUE                                                          
      IF(KC.NE.0) GO TO 20                                              
      RS(10)=XL(10)*XL(10)+YL(10)*YL(10)+ZL(10)*ZL(10)                  
      RL(10)=SQRT(RS(10))                                               
      FF(10)=AB1*EXP(-B1*RL(10))/RL(10)                                 
      VIR=VIR+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                                             
      VIR=VIR+FTEMP*RS(K-5)                                             
      RS(K)=XL(K)*XL(K)+YL(K)*YL(K)+ZL(K)*ZL(K)                         
      RL(K)=SQRT(RS(K))                                                 
      FF(K)=(AB3*EXP(-B3*RL(K))-AB4*EXP(-B4*RL(K)))/RL(K)               
      VIR=VIR+FF(K)*RS(K)                                               
 1140 CONTINUE                                                          
C                                                                       
C.....CALCULATE X-COMPONENT FORCES                                      
C                                                                       
   20 CONTINUE                                                          
      DO      1150      K=1,14                                          
 1150 GG(K)=FF(K)*XL(K)                                                 
      G110=GG(10)+GG(1)*C1                                              
      G23=GG(2)+GG(3)                                                   
      G45=GG(4)+GG(5)                                                   
      FX(IWO)=FX(IWO)+G110+GG(11)+GG(12)+C1*G23                         
      FX(JWO)=FX(JWO)-G110-GG(13)-GG(14)-C1*G45                         
      TT1=GG(1)*C2                                                      
      TT=G23*C2+TT1                                                     
      FX(IW1)=FX(IW1)+GG(6)+GG(7)+GG(13)+TT+GG(4)                       
      FX(IW2)=FX(IW2)+GG(8)+GG(9)+GG(14)+TT+GG(5)                       
      TT=G45*C2+TT1                                                     
      FX(JW1)=FX(JW1)-GG(6)-GG(8)-GG(11)-TT-GG(2)                       
      FX(JW2)=FX(JW2)-GG(7)-GG(9)-GG(12)-TT-GG(3)                       
C                                                                       
C.....CALCULATE Y-COMPONENT FORCES                                      
C                                                                       
      DO 1160 K=1,14                                                    
 1160 GG(K)=FF(K)*YL(K)                                                 
      G110=GG(10)+GG(1)*C1                                              
      G23=GG(2)+GG(3)                                                   
      G45=GG(4)+GG(5)                                                   
      FY(IWO)=FY(IWO)+G110+GG(11)+GG(12)+C1*G23                         
      FY(JWO)=FY(JWO)-G110-GG(13)-GG(14)-C1*G45                         
      TT1=GG(1)*C2                                                      
      TT=G23*C2+TT1                                                     
      FY(IW1)=FY(IW1)+GG(6)+GG(7)+GG(13)+TT+GG(4)                       
      FY(IW2)=FY(IW2)+GG(8)+GG(9)+GG(14)+TT+GG(5)                       
      TT=G45*C2+TT1                                                     
      FY(JW1)=FY(JW1)-GG(6)-GG(8)-GG(11)-TT-GG(2)                       
      FY(JW2)=FY(JW2)-GG(7)-GG(9)-GG(12)-TT-GG(3)                       
C                                                                       
C.....CALCULATE Z-COMPONENT FORCES                                      
C                                                                       
      DO 1170 K=1,14                                                    
 1170 GG(K)=FF(K)*ZL(K)                                                 
      G110=GG(10)+GG(1)*C1                                              
      G23=GG(2)+GG(3)                                                   
      G45=GG(4)+GG(5)                                                   
      FZ(IWO)=FZ(IWO)+G110+GG(11)+GG(12)+C1*G23                         
      FZ(JWO)=FZ(JWO)-G110-GG(13)-GG(14)-C1*G45                         
      TT1=GG(1)*C2                                                      
      TT=G23*C2+TT1                                                     
      FZ(IW1)=FZ(IW1)+GG(6)+GG(7)+GG(13)+TT+GG(4)                       
      FZ(IW2)=FZ(IW2)+GG(8)+GG(9)+GG(14)+TT+GG(5)                       
      TT=G45*C2+TT1                                                     
      FZ(JW1)=FZ(JW1)-GG(6)-GG(8)-GG(11)-TT-GG(2)                       
      FZ(JW2)=FZ(JW2)-GG(7)-GG(9)-GG(12)-TT-GG(3)                       
 1100 CONTINUE                                                          
      IW1=IW1+NATOMS                                                    
      IWO=IWO+NATOMS                                                    
      IW2=IW2+NATOMS                                                    
 1000 CONTINUE                                                          
CX20542 CALL DOOUTX( 542 )                                              
CY20542 CALL DOOUTY( 542 )                                              
CX10543 CALL DOINX( 543 )                                               
CY10543 CALL DOINY( 543 )                                               
C                                                                       
C.....DIVIDE FINAL FORCES BY MASSES                                     
C                                                                       
      DO 2000 I=1,NATMO,NATOMS                                          
      FX(  I)=FX(  I)*FHM                                               
      FY(  I)=FY(  I)*FHM                                               
      FZ(  I)=FZ(  I)*FHM                                               
      FX(I+2)=FX(I+2)*FHM                                               
      FY(I+2)=FY(I+2)*FHM                                               
      FZ(I+2)=FZ(I+2)*FHM                                               
      FX(I+1)=FX(I+1)*FOM                                               
      FY(I+1)=FY(I+1)*FOM                                               
      FZ(I+1)=FZ(I+1)*FOM                                               
 2000 CONTINUE                                                          
CX20544 CALL DOOUTX( 544 )                                              
CY20544 CALL DOOUTY( 544 )                                              
CX40028 CALL SBOUTX (28)                                                
CY40028 CALL SBOUTY (28)                                                
      RETURN                                                            
      END                                                               
C***********************************************************************
      SUBROUTINE CSHIFT(XA,XB,XMA,XMB,XL,BOXH,BOXL)                     
C***********************************************************************
C                                                                       
C.....THIS ROUTINE CHECKS AND SHIFTS IF NECESSARY THE TWO INTERACTING   
C     WATER MOLECULES TO WITHIN HALF THE BOX LENGTH                     
C                                                                       
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      DIMENSION XA(3),XB(3),XL(14)                                      
CX30029 CALL SBINX (29)                                                 
CY30029 CALL SBINY (29)                                                 
C                                                                       
      XL(1)=XMA-XMB                                                     
      XL(2)=XMA-XB(1)                                                   
      XL(3)=XMA-XB(3)                                                   
      XL(4)=XA(1)-XMB                                                   
      XL(5)=XA(3)-XMB                                                   
      XL(6)=XA(1)-XB(1)                                                 
      XL(7)=XA(1)-XB(3)                                                 
      XL(8)=XA(3)-XB(1)                                                 
      XL(9)=XA(3)-XB(3)                                                 
      XL(10)=XA(2)-XB(2)                                                
      XL(11)=XA(2)-XB(1)                                                
      XL(12)=XA(2)-XB(3)                                                
      XL(13)=XA(1)-XB(2)                                                
      XL(14)=XA(3)-XB(2)                                                
CX10545 CALL DOINX( 545 )                                               
CY10545 CALL DOINY( 545 )                                               
      DO 100 I=1,14                                                     
      IF(ABS(XL(I)).GT.BOXH) XL(I)=XL(I)-SIGN(BOXL,XL(I))               
  100 CONTINUE                                                          
CX20546 CALL DOOUTX( 546 )                                              
CY20546 CALL DOOUTY( 546 )                                              
CX40030 CALL SBOUTX (30)                                                
CY40030 CALL SBOUTY (30)                                                
      RETURN                                                            
      END                                                               
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,NATOMS
      COMMON /WWPOT/ QQ,A1,B1,A2,B2,A3,B3,A4,B4,AB1,AB2,AB3,AB4,C1,C2,  
     *               QQ2,QQ4,REF1,REF2,REF4                             
      COMMON /FRCNST/ FC11,FC12,FC13,FC33, FC111,FC333,FC112,FC113,     
     *                FC123,FC133, FC1111,FC3333,FC1112,                
     *                FC1122,FC1113,FC1123,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)                                       
CX30031 CALL SBINX (31)                                                 
CY30031 CALL SBINY (31)                                                 
C                                                                       
C.....INTRA-MOLECULAR POTENTIAL ENERGY                                  
C                                                                       
      POTA=0.0D0                                                        
      IW1=1                                                             
      IWO=2                                                             
      IW2=3                                                             
CX10547 CALL DOINX( 547 )                                               
CY10547 CALL DOINY( 547 )                                               
      DO 1000 I=1,NMOL                                                  
      XM(I)=C1*X(IWO)+C2*(X(IW1)+X(IW2))                                
      YM(I)=C1*Y(IWO)+C2*(Y(IW1)+Y(IW2))                                
      ZM(I)=C1*Z(IWO)+C2*(Z(IW1)+Z(IW2))                                
      R1=(X(IWO)-X(IW1))**2+(Y(IWO)-Y(IW1))**2+(Z(IWO)-Z(IW1))**2       
      R2=(X(IWO)-X(IW2))**2+(Y(IWO)-Y(IW2))**2+(Z(IWO)-Z(IW2))**2       
      RX=(X(IWO)-X(IW1))*(X(IWO)-X(IW2))+(Y(IWO)-Y(IW1))                
     *  *(Y(IWO)-Y(IW2))+(Z(IWO)-Z(IW1))*(Z(IWO)-Z(IW2))                
      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                                                         
      POTA=POTA+(FC11*(DR1S+DR2S)+FC33*DTS)*0.5D0                       
     *         +FC12*DR1*DR2+FC13*DRP*DT                                
     *   +(FC111*(DR1S*DR1+DR2S*DR2)+FC333*DTS*DT+FC112*DRP*DR1*DR2+    
     *     FC113*(DR1S+DR2S)*DT+FC123*DR1*DR2*DT+FC133*DRP*DTS)*ROHI    
      POTA=POTA                                                         
     *   +(FC1111*(DR1S*DR1S+DR2S*DR2S)+FC3333*DTS*DTS+                 
     *     FC1112*(DR1S+DR2S)*DR1*DR2+FC1122*DR1S*DR2S+                 
     *     FC1113*(DR1S*DR1+DR2S*DR2)*DT+FC1123*DRP*DR1*DR2*DT+         
     *     FC1133*(DR1S+DR2S)*DTS+FC1233*DR1*DR2*DTS+                   
     *     FC1333*DRP*DTS*DT)*ROHI2                                     
      IW1=IW1+NATOMS                                                    
      IWO=IWO+NATOMS                                                    
      IW2=IW2+NATOMS                                                    
 1000 CONTINUE                                                          
CX20548 CALL DOOUTX( 548 )                                              
CY20548 CALL DOOUTY( 548 )                                              
C                                                                       
C.....INTER-MOLECULAR POTENTIAL                                         
C                                                                       
      IW1=1                                                             
      IWO=2                                                             
      IW2=3                                                             
      POTR=0.0D0                                                        
      PTRF=0.0D0                                                        
CX10549 CALL DOINX( 549 )                                               
CY10549 CALL DOINY( 549 )                                               
      DO 2000 I=1,NMOL1                                                 
      JW1=IW1+NATOMS                                                    
      JWO=JW1+1                                                         
      JW2=JW1+2                                                         
      JMIN=I+1                                                          
      DO 2100 J=JMIN,NMOL                                               
      CALL CSHIFT(X(IW1),X(JW1),XM(I),XM(J),XL,BOXH,BOXL)               
      CALL CSHIFT(Y(IW1),Y(JW1),YM(I),YM(J),YL,BOXH,BOXL)               
      CALL CSHIFT(Z(IW1),Z(JW1),ZM(I),ZM(J),ZL,BOXH,BOXL)               
      KC=0                                                              
      DO 2110 K=1,9                                                     
      RS(K)=XL(K)*XL(K)+YL(K)*YL(K)+ZL(K)*ZL(K)                         
 2110 IF(RS(K).GT.CUT2) KC=KC+1                                         
      IF(KC.EQ.9) GO TO 100                                             
      DO 2120 K=1,9                                                     
      IF(RS(K).GT.CUT2) GOTO 10                                         
      RL(K)=SQRT(RS(K))                                                 
      GO TO 2120                                                        
   10 RL(K)=CUTOFF                                                      
      RS(K)=CUT2                                                        
 2120 CONTINUE                                                          
      POTR=POTR-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)                                               
      PTRF=PTRF-REF2*RS(1)-REF1*((RS(6)+RS(7)+RS(8)+RS(9))*0.5D0        
     *                           -RS(2)-RS(3)-RS(4)-RS(5))              
      IF(KC.GT.0) GOTO 100                                              
      DO 2130 K=10,14                                                   
 2130 RL(K)=SQRT(XL(K)*XL(K)+YL(K)*YL(K)+ZL(K)*ZL(K))                   
      POTR=POTR+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)))    
  100 JW1=JW1+NATOMS                                                    
      JWO=JWO+NATOMS                                                    
      JW2=JW2+NATOMS                                                    
 2100 CONTINUE                                                          
      IW1=IW1+NATOMS                                                    
      IWO=IWO+NATOMS                                                    
      IW2=IW2+NATOMS                                                    
 2000 CONTINUE                                                          
CX20550 CALL DOOUTX( 550 )                                              
CY20550 CALL DOOUTY( 550 )                                              
CX40032 CALL SBOUTX (32)                                                
CY40032 CALL SBOUTY (32)                                                
      RETURN                                                            
      END                                                               
