      PROGRAM BDNA                                                      
C================= MAIN  K.N.S. 3/17 84 === REV.001  S.K. 2/04 86 ======
C******THIS IS A SEQUENTIAL VERSION TO RUN ON THE IBM HOST  ************
C***********************************************************************
                   IMPLICIT                                             
     *             DOUBLE PRECISION                 (A-H,O-Z)           
C                                AS                                     
                   DOUBLE PRECISION        A, MATTER,OF,FACT,IN,LIQUID  
                   INTEGER       THIS,CODE,CAN,BE,USED                  
                   LOGICAL       MC,PLUS,MD,WORKS                       
                   REAL CPU, CPU0, WALL, WAL0                           
C                                                                       
C "MCMOLDYN"                                                            
                   DIMENSION   COM                          (116400)    
     *                        ,SITE                          (37965)    
     B                        ,BUFFR                         (20000)    
     Y                                                                  
                   COMMON      /WINDEX/         THIS           (110)    
     *                         /WJNDEX/         CODE            (10)    
     *                         /WKNDEX/         CAN             (10)    
     A                         /CONTRL/         BE              (25)    
     A                         /INDEXQ/         USED           (155)    
     T                         /TPARAM/         IN               (9)    
     T                         /LCONTR/         MC,MD,WORKS     (20)    
     O                         /CORRCC/         OF              (21)    
     *                         /SCRTZI/         LIQUID         (200)    
     L                         /SCRTZL/         AND             (20)    
     A                         /SCRTZR/         SOLID         (7812)    
     A                         /AVERIZ/         MATTER          (26)    
     K                         /LRCCEV/         STOP           (125)    
     S                         /PHYSIC/         FOR             (67)    
     O                         /MOUNTE/         MORE             (4)    
     N                         /CHARLO/         FACT             (3)    
     E                         /PARAMS/         GET              (6)    
     N                         /TFACTS/         A                (5)    
     *                         /LRCVAL/         USER            (26)    
                   COMMON      /IOUNIT/         MANUAL           (4)    
     1                         /TABELS/         FROM           (131)    
     9                         /CONSTS/         THE              (6)    
     8                         /CONVRT/         AUTHOR         (139)    
     4                         /RANDNO/         ST0P1  (128), IDUMMY    
CX00001 CALL MAINIX (1)                                                 
CY00001 CALL MAINIY (1)                                                 
CX60000 CALL CNTINX ()                                                  
CY60000 CALL CNTINY ()                                                  
C                                                       --------------- 
C                                                       --------------- 
C..... DEFINE FILES.....                                                
C            INPUT DATA         : 5 : MDSJE DATAIN                      
C            LINE_PRINTER OUTPUT: 6 : MDSJE DATAOUT                     
C            INPUT  FOR RESTART : 7 : MDSJE NRTIN                       
C            OUTPUT FOR RESTART : 8 : MDSJE NRTOUT                      
C     OPEN(5,FORM='FORMATTED',RECL=80)                                  
C     OPEN(66,FORM='FORMATTED',RECL=132)                                
C     OPEN(6,FORM='FORMATTED',RECL=132,FILE='FTN006,SIZE= 5000')        
C     OPEN(5, FILE = ':HOSTCHAR: MDSJE DATAIN  X' )                     
C     OPEN(6, FILE = ':HOSTCHAR: MDSJE DATAOUT B' )                     
C     OPEN(7,FORM='UNFORMATTED',RECL=32768)                             
C     OPEN(8,FORM='UNFORMATTED',RECL=32768)                             
C     OPEN(9, FILE = ':AATTO:HISTORY',FORM='UNFORMATTED',RECL=32768 )   
C     OPEN(9,FORM='FORMATTED',RECL=80)                                  
C     OPEN(10,FORM='FORMATTED',RECL=80)                                 
C     OPEN(11,FORM='FORMATTED',RECL=80)                                 
C     OPEN(12,FORM='FORMATTED',RECL=80)                                 
C     OPEN(13,FORM='FORMATTED',RECL=80)                                 
C     OPEN(14,FORM='FORMATTED',RECL=80)                                 
C     OPEN(50,FORM='UNFORMATTED',RECL=32768)                            
C     OPEN(51,FORM='UNFORMATTED',RECL=32768)                            
C     OPEN(61,FORM='UNFORMATTED',RECL=32768)                            
C                CPU1=0.D0                                              
C                CPU2=0.D0                                              
C                CPU3=0.D0                                              
C                CPU4=0.D0                                              
             CALL CPUTIM (CPU0)                                         
             CALL ELAPSE (WAL0)                                         
C            CALL XUFLOW                                                
      OPEN(1,FILE='NAV',STATUS='UNKNOWN')                               
      OPEN(5,FILE='NAI5',STATUS='OLD')                                  
      OPEN(6,FILE='NAO6',STATUS='UNKNOWN')                              
      OPEN(8,FILE='NAO8',FORM='UNFORMATTED',STATUS='UNKNOWN')           
      OPEN(9,FILE='NAI9',STATUS='OLD')                                  
      OPEN(10,FILE='NAI10',STATUS='OLD')                                
      OPEN(11,FILE='NAI11',STATUS='OLD')                                
      OPEN(12,FILE='NAI12',STATUS='OLD')                                
      OPEN(13,FILE='NAI13',STATUS='OLD')                                
      OPEN(14,FILE='NAI14',STATUS='OLD')                                
C     OPEN(15,FILE='NAI15',STATUS='UNKNOWN')                            
      OPEN(20,FILE='NAO20',FORM='UNFORMATTED',STATUS='UNKNOWN')         
      OPEN(66,FILE='NAO66',STATUS='UNKNOWN')                            
      WRITE(1,1)                                                        
1     FORMAT(1X,
     * 'OUTPUT FOR PERFECT CLUB BENCHMARK: BDNA INSTRUMENTED VERSION'/
     *       1X,'$Revision: 1.A.1 $ $Author: psharma $'/)                   
                   CALL INPUT                                           
C                                                                       
                   CALL INDEX1                                          
C                                                                       
                   CALL WATER                                           
C                                                                       
                   CALL START  (COM,SITE)                               
C                                                                       
                   CALL CYCLES (COM,SITE,BUFFR)                         
C                                                                       
C                                                                       
              CALL CPUTIM (CPU)                                         
              CALL ELAPSE(WALL)                                         
              WRITE(66,555) CPU-CPU0,MOD(WALL-WAL0+1000000.0,1000000.0) 
  555         FORMAT(4X,'TOTAL CPU TIME USED  =',F15.5,' SECONDS',/,    
     1               4X,'TOTAL WALL TIME USED =',F15.5,' SECONDS',/)    
              WRITE(1,*)                                                
              WRITE(1,*) 'ELAPSED CPU TIME IN SECONDS: ',CPU-CPU0       
C MFLOPS = MFLOP 1 CPU CRAY X-MP/CPU TIME                               
              WRITE(1,*) 'MFLOPS RATE: ', 905.534/(CPU-CPU0)            
              WRITE(1,*) 'ELAPSED WALLCLOCK TIME IN SECONDS: ',         
     *                    MOD(WALL-WAL0+1000000.0,1000000.0)            
CX70000 CALL CNTOTX ()                                                  
CY70000 CALL CNTOTY ()                                                  
CX50000 CALL MAINOX ()                                                  
CY50000 CALL MAINOY ()                                                  
                   STOP                                                 
C                                                                       
                   END                                                  
C============= BLOCK DATAK.N.S. 3/17 84 === REV.001  S.K. 2/04 86 ======
C                                                                       
      BLOCK DATA                                                        
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      LOGICAL         UNTRUE                                            
      COMMON /LCONTR/ UNTRUE(22)                                        
     X       /CONVRT/ ZEROEM(139)                                       
     X       /AVERIZ/ VALUES(26)                                        
     X       /IOUNIT/ NUMBER(4)                                         
     X       /PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK           
     X               ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)         
     X               ,R2(3),EWW,EWI,EWA,EII,EIA                         
     X       /CONTRL/ NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP   
     X               ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION       
     X               ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS               
     X       /CONSTS/ AVSNO,BOLTZ,PLANCK,PI,ELECHG,EPS0                 
     X       /INDEXQ/ NCOM,NSITE,NBUFFR,NREST(152)                      
      COMMON/TABELS/AD(131)                                             
     X      /WINDEX/IE(110)                                             
     X      /WJNDEX/IT(10)                                              
     X      /WKNDEX/IG(10)                                              
     X      /TPARAM/AH(9)                                               
      COMMON/CORRCC/AI(21)                                              
     X      /SCRTZI/AJ(200)                                             
     X      /SCRTZL/AK(20)                                              
     X      /SCRTZR/AL2(7812)                                           
     X      /TFACTS/AM(5)                                               
     X      /PARAMS/AN(6)                                               
     X      /LRCCEV/AP(125)                                             
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *     
      DATA   NCOM      /116400/                                         
     X       NSITE     /37965/                                          
     X       NBUFFR    /20000/                                          
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *     
      DATA NUMBER /5,6,7,8/                                             
     X     UNTRUE /22*.FALSE./                                          
     X     VALUES /26*0.D0/                                             
     X     ZEROEM /139*0.D0/                                            
     X     NREST  /152*0/                                               
     X     NATOMS,NION  /635,20/                                        
     X     NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP,NSITES       
     X    ,NSPECI,NSITET,NENT,NSS,NAVRS,NOP3,NOP4,NSS3                  
     X     /1,0,4*1,2,1520,4,1,0,1500,20,0,0,2000,5*0/                  
     X     PI,BOLTZ,AVSNO/3.141592653589793D0,1.3806D-23,6.0222D+23/    
     X     ELECHG,EPS0,PLANCK/1.602D-19,8.854D-12,6.6256D-34/           
     X     TTRAN,TROT,RHO,VOLM,DT,FNOP,BREAK,DM,QM/2*300.D0,0.D0        
     X     ,300.D0,2*0.D0,1.D0,0.D0,1.D-15,1520.D0,1.0D0,36*0.D0/       
      DATA AD/131*0.D0/                                                 
     X     IE/110*0/                                                    
     X     IT/10*0/                                                     
     X     IG/10*0/                                                     
     X     AH/9*0.D0/                                                   
      DATA AI/21*0.D0/                                                  
     X     AJ/200*0.D0/                                                 
     X     AK/20*0.D0/                                                  
     X     AL2/7812*0.D0/                                               
     X     AM/5*0.D0/                                                   
     X     AN/6*0.D0/                                                   
     X     AP/125*0.D0/                                                 
      END                                                               
C================= START K.N.S. 3/24 84 === REV.001  S.K. 2/04 86 ======
C                                                                       
      SUBROUTINE START(T,S)                                             
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION T(*),S(*)                                               
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/INDEXQ/ NWI,NWJ,NWK                                        
      COMMON/WINDEX/ NOR,IQ, IW, IA,IM,IF,IT,IP,IO,IX(99)               
      COMMON/WJNDEX/ JF1,JX(9)                                          
CX30003 CALL SBINX (3)                                                  
CY30003 CALL SBINY (3)                                                  
      CALL ZERO(NWI,T(1))                                               
      CALL ZERO(NWJ,S(1))                                               
      IF(INMC) THEN                                                     
      CALL MCIN (T(IX(1)) ,T(IX(2))   ,T(IX(3))                         
     X          ,T(IX(IQ)),T(IX(IQ+1)),T(IX(IQ+2)),T(IX(IQ+3))          
     X          ,T(IX(IW)),T(IX(IW+1)),T(IX(IW+2))                      
     X          ,T(IX(IA)),T(IX(IA+1)),T(IX(IA+2))                      
     X          ,S(JX(1)) ,S(JX(2))   ,S(JX(3))                         
     X          ,S(JX(4)) ,S(JX(5))   ,S(JX(6))   )                     
      ENDIF                                                             
      IF(LRES) CALL RESTAR(T,.TRUE.,9)                                  
CX40004 CALL SBOUTX (4)                                                 
CY40004 CALL SBOUTY (4)                                                 
      RETURN                                                            
      END                                                               
C================= CYCLESK.N.S. 3/25 84 === REV.005  S.K. 2/04 86 ======
C                                                                       
      SUBROUTINE CYCLES(T,S,W)                                          
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION T(*),S(*),W(*)                                          
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      COMMON/WINDEX/ NOR,IQ, IW, IA,IM,IF,IT,IP,IO,IX(99)               
      COMMON/WJNDEX/ JF1,JX(9)                                          
      COMMON/WKNDEX/ KX(10)                                             
CX30005 CALL SBINX (5)                                                  
CY30005 CALL SBINY (5)                                                  
      CALL TYPE (3)                                                     
      IF(.NOT.LRES.AND.LPC) THEN                                        
      NSP=NSPECI(1)                                                     
      CALL UPDATE(T(IX(IQ))  ,T(IX(IQ+1)) ,T(IX(IQ+2)) ,T(IX(IQ+3))     
     X           ,T(IX(IF))  ,T(IX(IF+1)) ,T(IX(IF+2))                  
     X           ,S(JX(1))   ,S(JX(2))    ,S(JX(3)))                    
      CALL ACTFOR(T(IX(1))   ,T(IX(2)),T(IX(3)),                        
     X            T(IX(IF))  ,T(IX(IF+1)) ,T(IX(IF+2)),                 
     X            T(IX(IT))  ,T(IX(IT+1)) ,T(IX(IT+2)),                 
     X            S(JX(1))   ,S(JX(2))    ,S(JX(3)),                    
     X            S(JX(4))   ,S(JX(5))    ,S(JX(6)))                    
      CALL PCINIT(T(IX(7))   ,T(IX(8))    ,T(IX(9))                     
     X           ,T(IX(IQ+0)),T(IX(IQ+1)) ,T(IX(IQ+2)) ,T(IX(IQ+3))     
     X           ,T(IX(IQ+4)),T(IX(IQ+5)) ,T(IX(IQ+6)) ,T(IX(IQ+7))     
     X           ,T(IX(IA))  ,T(IX(IA+1)) ,T(IX(IA+2))                  
     X           ,T(IX(IA+3)),T(IX(IA+4)) ,T(IX(IA+5))                  
     X           ,T(IX(IF))  ,T(IX(IF+1)) ,T(IX(IF+2))                  
     X           ,T(IX(IT))  ,T(IX(IT+1)) ,T(IX(IT+2)))                 
      CALL ZERO(20*NSP,T(IX(IQ+8)))                                     
      CALL ZERO(15*NSP,T(IX(IA+6)))                                     
      CALL ZERO( 9*NOP,T(IX(9)))                                        
      ENDIF                                                             
    1 NSTEP=NSTEP+1                                                     
      IF(LPC)                                                           
     XCALL PREDIC(T)                                                    
      CALL UPDATE(T(IX(IQ)) ,T(IX(IQ+1)) ,T(IX(IQ+2)) ,T(IX(IQ+3))      
     X           ,T(IX(IF)) ,T(IX(IF+1)) ,T(IX(IF+2))                   
     X           ,S(JX(1))  ,S(JX(2))    ,S(JX(3)))                     
C                CPU1=0.D0                                              
C                CPU2=0.D0                                              
C                CPU3=0.D0                                              
C                CPU4=0.D0                                              
      CALL ACTFOR(T(IX(1))   , T(IX(2))    ,  T(IX(3)),                 
     X            T(IX(IF))  , T(IX(IF+1)) ,  T(IX(IF+2)),              
     X            T(IX(IT))  , T(IX(IT+1)) ,  T(IX(IT+2)),              
     X            S(JX(1))   , S(JX(2))    ,  S(JX(3)),                 
     X            S(JX(4))   , S(JX(5))    ,  S(JX(6)))                 
      IF(LLF)                                                           
     XCALL MOTION( T(IX(1)) ,T(IX(2)) ,T(IX(3))                         
     X ,T(IX(4)) ,T(IX(5)) ,T(IX(6)) ,T(IX(7)) ,T(IX(8)) ,T(IX(9))      
     X ,T(IX(10)),T(IX(11)),T(IX(12)),T(IX(13)),T(IX(14)),T(IX(15))     
     X ,T(IX(16)),T(IX(17)),T(IX(18)),T(IX(19)),T(IX(20)),T(IX(21))     
     X ,T(IX(22)),S(JX(1)) ,S(JX(2)) ,S(JX(3)) ,S(JX(4)))               
      IF(LPC)                                                           
     XCALL CORREC(T,S)                                                  
      CALL VALUES                                                       
      IF(LTS) THEN                                                      
      IF(LLF)THEN                                                       
      CALL RNMTRN(T(IX(IW)),T(IX(IW+1)),T(IX(IW+2)))                    
      CALL RNMROT(T(IX(IA)),T(IX(IA+1)),T(IX(IA+2)))                    
      ENDIF                                                             
        IF(LPC) THEN                                                    
          II=6                                                          
CX10501 CALL DOINX( 501 )                                               
CY10501 CALL DOINY( 501 )                                               
          DO 3 I=1,4                                                    
          CALL TSCALE(T(IX(1+II)),T(IX(2+II)),T(IX(3+II)),1,.TRUE.)     
          II=II+3                                                       
    3     CONTINUE                                                      
CX20502 CALL DOOUTX( 502 )                                              
CY20502 CALL DOOUTY( 502 )                                              
      CALL RNMTRN(T(IX(IW)),T(IX(IW+1)),T(IX(IW+2)))                    
        IF(.NOT.LQ2PC) THEN                                             
            II=3                                                        
CX10503 CALL DOINX( 503 )                                               
CY10503 CALL DOINY( 503 )                                               
            DO 4 I=1,6                                                  
            CALL TSCALE(T(IX(IA+II)),T(IX(IA+II+1)),T(IX(IA+II+2))      
     X                 ,1,.FALSE.)                                      
            II=II+3                                                     
    4       CONTINUE                                                    
CX20504 CALL DOOUTX( 504 )                                              
CY20504 CALL DOOUTY( 504 )                                              
      CALL RNMROT(T(IX(IA)),T(IX(IA+1)),T(IX(IA+2)))                    
          ELSE                                                          
      CALL RNMROT(T(IX(IA)),T(IX(IA+1)),T(IX(IA+2)))                    
      CALL SCALER(T(IX(IQ+0)),T(IX(IQ+1)),T(IX(IQ+2)),T(IX(IQ+3)),      
     X            T(IX(IQ+4)),T(IX(IQ+5)),T(IX(IQ+6)),T(IX(IQ+7)),      
     X            T(IX(IA+0)),T(IX(IA+1)),T(IX(IA+2)))                  
            ENDIF                                                       
            ENDIF                                                       
            ENDIF                                                       
      IF(LSAVE) THEN                                                    
      IF(MOD(NSTEP,2).EQ.0) CALL RESTAR(T,.FALSE.,8)                    
      ENDIF                                                             
      IF(MOD(NSTEP,2).EQ.0) THEN                                        
      CALL UPDATE(T(IX(IQ)) ,T(IX(IQ+1)) ,T(IX(IQ+2)) ,T(IX(IQ+3))      
     X           ,T(IX(IF)) ,T(IX(IF+1)) ,T(IX(IF+2))                   
     X           ,S(JX(1))  ,S(JX(2))    ,S(JX(3)))                     
      CALL DUMP(T(IX(1)),T(IX(2)),T(IX(3)),T(IX(IW)),T(IX(IW+1))        
     X          ,T(IX(IW+2)),S(JX(1)),S(JX(2)),S(JX(3))                 
     X          ,T(IX(IA+0)),T(IX(IA+1)),T(IX(IA+2))                    
     X          ,T(IX(IQ+0)),T(IX(IQ+1)),T(IX(IQ+2)),T(IX(IQ+3)))       
      ENDIF                                                             
      IF(NSTEP.LT.NSTEPS) GO TO 1                                       
      IF(NAVRS.NE.0) CALL TYPE(4)                                       
CX40006 CALL SBOUTX (6)                                                 
CY40006 CALL SBOUTY (6)                                                 
      RETURN                                                            
      END                                                               
C================= ROTATEK.N.S. 3/26 84 === REV.000  X.X.  X/X XX ======
C                                                                       
      SUBROUTINE ROTATE(QA,QB,QC,QD,TX,TY,TZ,NOP,NBTOP)                 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      LOGICAL NBTOP                                                     
      DIMENSION QA(*),QB(*),QC(*),QD(*),TX(*),TY(*),TZ(*)               
CX30007 CALL SBINX (7)                                                  
CY30007 CALL SBINY (7)                                                  
      SIGN=1.D0                                                         
      IF(.NOT.NBTOP) SIGN=-1.D0                                         
CX10505 CALL DOINX( 505 )                                               
CY10505 CALL DOINY( 505 )                                               
      DO 100 I=1,NOP                                                    
      AA=QA(I)*QA(I)                                                    
      BB=QB(I)*QB(I)                                                    
      CC=QC(I)*QC(I)                                                    
      DD=QD(I)*QD(I)                                                    
      AB=QA(I)*QB(I)                                                    
      AC=QA(I)*QC(I)                                                    
      AD=QA(I)*QD(I)*SIGN                                               
      BC=QB(I)*QC(I)                                                    
      BD=QB(I)*QD(I)*SIGN                                               
      CD=QC(I)*QD(I)*SIGN                                               
      RX=TX(I)*(-AA+BB-CC+DD)                                           
     X  +TY(I)*(+CD-AB)*2.D0                                            
     X  +TZ(I)*(+BC+AD)*2.D0                                            
      RY=TX(I)*(-AB-CD)*2.D0                                            
     X  +TY(I)*(+AA-BB-CC+DD)                                           
     X  +TZ(I)*(+BD-AC)*2.D0                                            
      RZ=TX(I)*(+BC-AD)*2.D0                                            
     X  +TY(I)*(-AC-BD)*2.D0                                            
     X  +TZ(I)*(-AA-BB+CC+DD)                                           
      TX(I)=RX                                                          
      TY(I)=RY                                                          
      TZ(I)=RZ                                                          
  100 CONTINUE                                                          
CX20506 CALL DOOUTX( 506 )                                              
CY20506 CALL DOOUTY( 506 )                                              
CX40008 CALL SBOUTX (8)                                                 
CY40008 CALL SBOUTY (8)                                                 
      RETURN                                                            
      END                                                               
C================= TSCALEK.N.S. 4/2  84 === REV.000  X.X.  X/X XX ======
C                                                                       
      SUBROUTINE TSCALE(X,Y,Z,NNN,LSCALE)                               
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      LOGICAL LFLUCT,LSCALE                                             
      DIMENSION X(*),Y(*),Z(*)                                          
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/AVERIZ/ V(13),VV(13)                                       
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
CX30009 CALL SBINX (9)                                                  
CY30009 CALL SBINY (9)                                                  
      IF(.NOT.LTS) THEN                                                 
CX40010 CALL SBOUTX (10)                                                
CY40010 CALL SBOUTY (10)                                                
      RETURN                                                            
      ENDIF                                                             
      IF(NSTEP.GE.NSSTS) WRITE(66,10) NSTEP                             
   10 FORMAT(//1X,'*** TEMP. SCALING SWITHED OFF AT TIME STEP ',I6//)   
      IF(NSTEP.GE.NSSTS) LTS=.FALSE.                                    
      IF(MOD(NSTEP,NSBTS).NE.0) THEN                                    
CX40010 CALL SBOUTX (10)                                                
CY40010 CALL SBOUTY (10)                                                
      RETURN                                                            
      ENDIF                                                             
      UPTO=REAL(NNN)                                                    
      I=0                                                               
CX10507 CALL DOINX( 507 )                                               
CY10507 CALL DOINY( 507 )                                               
      DO 31 N=1,NTYPES                                                  
      NSP=NSPECI(N)                                                     
      NS=NSITES(N)                                                      
      TSCAL=1.D0                                                        
      RSCAL=1.D0                                                        
      IF(TTS(N).NE.0.D0) TSCAL=SQRT(TTRAN(N)/TTS(N))                    
      IF(RTS(N).NE.0.D0) RSCAL=SQRT(TROT(N)/RTS(N))                     
      IF(LSCALE) SC1=TSCAL                                              
      IF(.NOT.LSCALE) SC1=RSCAL                                         
       IF(SC1.GT.1.5D0.OR.SC1.LT..5D0) THEN                             
        SC1=0.D0                                                        
        ELSE                                                            
        SC1=SC1**UPTO                                                   
       ENDIF                                                            
      IF(.NOT.LSCALE.AND.NS.EQ.1) GO TO 31                              
      DO 30 J=1,NSP                                                     
      I=I+1                                                             
      X(I)=X(I)*SC1                                                     
      Y(I)=Y(I)*SC1                                                     
      Z(I)=Z(I)*SC1                                                     
   30 CONTINUE                                                          
   31 CONTINUE                                                          
CX20508 CALL DOOUTX( 508 )                                              
CY20508 CALL DOOUTY( 508 )                                              
      NAVRS=0                                                           
CX10509 CALL DOINX( 509 )                                               
CY10509 CALL DOINY( 509 )                                               
      DO 40 I=1,13                                                      
      V(I)=0.D0                                                         
      VV(I)=0.D0                                                        
   40 CONTINUE                                                          
CX20510 CALL DOOUTX( 510 )                                              
CY20510 CALL DOOUTY( 510 )                                              
CX40010 CALL SBOUTX (10)                                                
CY40010 CALL SBOUTY (10)                                                
      RETURN                                                            
      END                                                               
C================= UPDATEK.N.S. 3/27 84 === REV.000  X.X.  X/X XX ======
C                                                                       
      SUBROUTINE UPDATE(QA0,QB0,QC0,QD0,TX,TY,TZ,SX,SY,SZ)              
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION QA0(*),QB0(*),QC0(*),QD0(*),TX(*),TY(*),TZ(*)           
      DIMENSION  SX(*) ,SY(*) ,SZ(*)                                    
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
      COMMON/WINDEX/ NOR,IQ, IW, IA,IM,IF,IT,IP,IO,IX(99)               
CX30011 CALL SBINX (11)                                                 
CY30011 CALL SBINY (11)                                                 
      NSP=NSPECI(1)                                                     
      NS =NSITES(1)                                                     
CX10511 CALL DOINX( 511 )                                               
CY10511 CALL DOINY( 511 )                                               
C     CALL ZERO(NSS,SX(1))                                              
C     CALL ZERO(NSS,SY(1))                                              
C     CALL ZERO(NSS,SZ(1))                                              
C     CALL ZERO(NSP,TX(1))                                              
C     CALL ZERO(NSP,TY(1))                                              
C     CALL ZERO(NSP,TZ(1))                                              
      DO 10 I=1,NSS                                                     
      SX(I)=0.0D0                                                       
      SY(I)=0.0D0                                                       
      SZ(I)=0.0D0                                                       
 10   CONTINUE                                                          
CX20512 CALL DOOUTX( 512 )                                              
CY20512 CALL DOOUTY( 512 )                                              
CX10513 CALL DOINX( 513 )                                               
CY10513 CALL DOINY( 513 )                                               
      DO 11 I=1,NSP                                                     
      TX(I)=0.0D0                                                       
      TY(I)=0.0D0                                                       
      TZ(I)=0.0D0                                                       
 11   CONTINUE                                                          
CX20514 CALL DOOUTX( 514 )                                              
CY20514 CALL DOOUTY( 514 )                                              
CX10515 CALL DOINX( 515 )                                               
CY10515 CALL DOINY( 515 )                                               
      DO 100 J=1,NS                                                     
      EX=DSITE(1,J)                                                     
      EY=DSITE(2,J)                                                     
      EZ=DSITE(3,J)                                                     
      DO 50 I=1,NSP                                                     
      TX(I)=EX                                                          
      TY(I)=EY                                                          
      TZ(I)=EZ                                                          
   50 CONTINUE                                                          
      CALL ROTATE(QA0,QB0,QC0,QD0,TX,TY,TZ,NSP,.FALSE.)                 
      DO 100 I=1,NSP                                                    
      ISP=(I-1)*NS+J                                                    
      SX(ISP)=TX(I)                                                     
      SY(ISP)=TY(I)                                                     
      SZ(ISP)=TZ(I)                                                     
  100 CONTINUE                                                          
CX20516 CALL DOOUTX( 516 )                                              
CY20516 CALL DOOUTY( 516 )                                              
CX10517 CALL DOINX( 517 )                                               
CY10517 CALL DOINY( 517 )                                               
      DO 12 I=1,NSP                                                     
      TX(I)=0.0D0                                                       
      TY(I)=0.0D0                                                       
      TZ(I)=0.0D0                                                       
 12   CONTINUE                                                          
CX20518 CALL DOOUTX( 518 )                                              
CY20518 CALL DOOUTY( 518 )                                              
CX40012 CALL SBOUTX (12)                                                
CY40012 CALL SBOUTY (12)                                                
C     CALL ZERO(NSP,TX(1))                                              
C     CALL ZERO(NSP,TY(1))                                              
C     CALL ZERO(NSP,TZ(1))                                              
      RETURN                                                            
      END                                                               
C================= MCINK.N.S.  6/4 85 === REV.001  S.K.  2/04 86 ====== 
C                                                                       
      SUBROUTINE MCIN(X0,Y0,Z0,QA0,QB0,QC0,QD0,WX0,WY0,WZ0              
     X               ,AX0,AY0,AZ0,SX,SY,SZ,FS1,FS2,FS3)                 
C********                                                               
C******** THIS SUBROUTINE IS WRITTEN FOR 500 WATER MOLECULES            
C******** AND 650 ATOMS OF THE BIOMOLECULE. FOR LARGER                  
C******** SYSTEMS YOU HAVE TO INCREASE THE DIMENSIONS.                  
C********                                                               
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION  X0(*), Y0(*), Z0(*)                                    
      DIMENSION QA0(*),QB0(*),QC0(*),QD0(*)                             
      DIMENSION WX0(*),WY0(*),WZ0(*)                                    
      DIMENSION AX0(*),AY0(*),AZ0(*)                                    
      DIMENSION  SX(*), SY(*), SZ(*),FS1(*),FS2(*),FS3(*)               
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONTRL/ NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP    
     X               ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION       
     X               ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS               
      COMMON/SCRTZR/ A(3),B(3),C(3),AA(3,650),BB(3,650),CC(3,650)       
     X               ,AAA(650),BBB(650),CCC(650),AAAA,BBBB,CCCC         
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
      DIMENSION EX(3),EY(3),EZ(3),QP(4)                                 
      DIMENSION PQDATA(11)                                              
      DIMENSION CHAR(650)                                               
      DIMENSION AXION(50),AYION(50),AZION(50)                           
      DIMENSION XAT(650),YAT(650),ZAT(650)                              
      DIMENSION  CX(650)                                                
      DIMENSION FX(500),FY(500),FZ(500),GX(500),GY(500),GZ(500)         
      DIMENSION X6(500),Y6(500),Z6(500)                                 
      CHARACTER*8 NAME                                                  
      DIMENSION NAME(650),ICLAS(650),JCLAS(650),ITITLE(3)               
      DIMENSION FX06(50),FH06(50),FI06(50),                             
     &FX12(50),FH12(50),FI12(50),FX01(50),FH01(50),FI01(50)             
CX30013 CALL SBINX (13)                                                 
CY30013 CALL SBINY (13)                                                 
      WRITE(66,2402)                                                    
 2402 FORMAT(///1X,'  OUTPUT FROM MCIN : '////)                         
C                                                                       
C*********** READ IN MC INPUT DECK *************************************
C                                                                       
C               ION-ION                                                 
C                                                                       
      READ(10,152) NAT,N,NION,MCLAS,NCLAS                               
      READ(10,429) VI06,VI12,VI01                                       
C-----------------------------------------------------------------      
C                                                                       
C               WATER-ION                                               
C                                                                       
      READ(10,129)PQDATA(5),PQDATA(6),PQDATA(4)                         
      READ(10,129)PQDATA(2),PQDATA(3),PQDATA(1)                         
C-----------------------------------------------------------------      
C-----------------------------------------------------------------      
      WRITE(66,923) NAT,N,NION,MCLAS,NCLAS                              
  923 FORMAT(/1X,' NAT,N,NION,MCLAS,NCLAS :',5I6)                       
C-----------------------------------------------------------------      
C                                                                       
C*********** READ IN WATER CONFIGURATION *******************************
C                                                                       
      READ(12,130) (I,X6(L),Y6(L),Z6(L),FX(L),FY(L),FZ(L),              
     *                                 GX(L),GY(L),GZ(L),L=1,N)         
C-----------------------------------------------------------------      
      WRITE(66,924) X6(1),Y6(1),Z6(1),FX(1),FY(1),FZ(1),GX(1),GY(1)     
     X ,GZ(1)                                                           
  924 FORMAT(/1X,' 1ST WATER: ',3F7.3,6F9.5)                            
C-----------------------------------------------------------------      
C                                                                       
C*********** READ IN ION CONFIGURATION *********************************
C                                                                       
      READ(12,131)(I,I,AXION(K),AYION(K),AZION(K),K =1,NION)            
C-----------------------------------------------------------------      
      WRITE(66,925) AXION(1),AYION(1),AZION(1)                          
  925 FORMAT(/1X,' 1ST ION  : ',3F18.14)                                
C-----------------------------------------------------------------      
C                                                                       
C*********** READ IN DNA COORDINATES ***********************************
C                                                                       
      READ (11,128)(NAME(K),XAT(K),YAT(K),ZAT(K),CHAR(K),ICLAS(K),      
     &                                         JCLAS(K),K=1,NAT)        
C-----------------------------------------------------------------      
      WRITE(66,926) NAME(1),XAT(1),YAT(1),ZAT(1),CHAR(1),ICLAS(1)       
     X,JCLAS(1)                                                         
  926 FORMAT(/1X,' 1ST ATOM : ',A8,3F15.8,F10.6,2I3)                    
C-----------------------------------------------------------------      
C                                                                       
C*********** READ IN POTENTIAL INFORMATION: ****************************
C                                                                       
C               ION-ION                                                 
C                                                                       
C-----------------------------------------------------------------      
      WRITE(66,927) VI06,VI12,VI01                                      
  927 FORMAT(/1X,' ION-ION PARAMETERS: ',3(1X,D12.4))                   
      WRITE(66,928) (PQDATA(I),I=1,6)                                   
  928 FORMAT(/1X,' ION-H2O PARAMETERS: ',6(1X,D12.4))                   
C-----------------------------------------------------------------      
CX10519 CALL DOINX( 519 )                                               
CY10519 CALL DOINY( 519 )                                               
C                                                                       
C                ION-DNA                                                
C                                                                       
      DO 3331 K=1,MCLAS                                                 
      READ(13,140) FI06(K)                                              
      READ(13,140) FI12(K)                                              
 3331 CONTINUE                                                          
CX20520 CALL DOOUTX( 520 )                                              
CY20520 CALL DOOUTY( 520 )                                              
CX10521 CALL DOINX( 521 )                                               
CY10521 CALL DOINY( 521 )                                               
      DO 3332 K=1,MCLAS                                                 
      READ(13,140) FI01(K)                                              
      READ(13,140) DUMMY                                                
 3332 CONTINUE                                                          
C-----------------------------------------------------------------      
CX20522 CALL DOOUTX( 522 )                                              
CY20522 CALL DOOUTY( 522 )                                              
      WRITE(66,929) FI06(1),FI12(1),FI01(1)                             
  929 FORMAT(/1X,' 1ST ION-DNA PARAMETER: ',3(1X,D12.4))                
C-----------------------------------------------------------------      
CX10523 CALL DOINX( 523 )                                               
CY10523 CALL DOINY( 523 )                                               
C                                                                       
C                WATER-DNA                                              
C                                                                       
      DO 3333 K=1,NCLAS                                                 
      READ(14,140) FX06(K)                                              
      READ(14,140) FX12(K)                                              
 3333 CONTINUE                                                          
CX20524 CALL DOOUTX( 524 )                                              
CY20524 CALL DOOUTY( 524 )                                              
CX10525 CALL DOINX( 525 )                                               
CY10525 CALL DOINY( 525 )                                               
      DO 3335 K=1,NCLAS                                                 
      READ(14,140) FH06(K)                                              
      READ(14,140) FH12(K)                                              
 3335 CONTINUE                                                          
CX20526 CALL DOOUTX( 526 )                                              
CY20526 CALL DOOUTY( 526 )                                              
CX10527 CALL DOINX( 527 )                                               
CY10527 CALL DOINY( 527 )                                               
      DO 3334 K=1,NCLAS                                                 
      READ(14,140) FX01(K)                                              
      READ(14,140) DUMMY                                                
 3334 CONTINUE                                                          
CX20528 CALL DOOUTX( 528 )                                              
CY20528 CALL DOOUTY( 528 )                                              
CX10529 CALL DOINX( 529 )                                               
CY10529 CALL DOINY( 529 )                                               
      DO 3336 K=1,NCLAS                                                 
      READ(14,140) FH01(K)                                              
      READ(14,140) DUMMY                                                
 3336 CONTINUE                                                          
C-----------------------------------------------------------------      
CX20530 CALL DOOUTX( 530 )                                              
CY20530 CALL DOOUTY( 530 )                                              
      WRITE(66,930) FX06(1),FX12(1),FX01(1)                             
  930 FORMAT(/1X,' 1ST   O-DNA PARAMETER: ',3(1X,D12.4))                
      WRITE(66,931) FH06(1),FH12(1),FH01(1)                             
  931 FORMAT(/1X,' 1ST   H-DNA PARAMETER: ',3(1X,D12.4))                
C-----------------------------------------------------------------      
C*********************************************************************  
  128 FORMAT(A8,3F15.8,F10.6,2I3)                                       
  129 FORMAT(5X,E20.10)                                                 
  130 FORMAT(I5,3F7.3,6F9.5)                                            
  131 FORMAT(2I4,3F18.14)                                               
  140 FORMAT(F20.10,I3)                                                 
  152 FORMAT(   8I10)                                                   
  429 FORMAT(    3E20.10 )                                              
C*********************************************************************  
      QOXY=-0.6331D0                                                    
      QHYD= 0.3165D0                                                    
      PI = 4.D0*ATAN(1.D0)                                              
      AVNO = 6.0222D23                                                  
      CALKJ = 4.184D0                                                   
      HKS=1.D3*CALKJ/AVNO                                               
      QION=1.D0                                                         
      RF=1.D-10/UNITL                                                   
CX10531 CALL DOINX( 531 )                                               
CY10531 CALL DOINY( 531 )                                               
      DO 374 K=1,NAT                                                    
C**   CHAR(K)=REAL(NUCL(K))-CHAR(K)                                     
      J=ICLAS(K)                                                        
      L=JCLAS(K)                                                        
      BB(1,K)=-HKS*FX06(L)*RF**6 /UNITE                                 
      BB(2,K)=-HKS*FH06(L)*RF**6 /UNITE                                 
      BB(3,K)= BB(2,K)                                                  
      BBB(K) =-HKS*FI06(J)*RF**6 /UNITE                                 
      CC(1,K)= HKS*FX12(L)*RF**12/UNITE                                 
      CC(2,K)= HKS*FH12(L)*RF**12/UNITE                                 
      CC(3,K)=CC(2,K)                                                   
      CCC(K) = HKS*FI12(J)*RF**12/UNITE                                 
      AA(1,K)= CHAR(K)*FX01(L)*QOXY*COULF                               
      AA(2,K)= CHAR(K)*FH01(L)*QHYD*COULF                               
      AA(3,K)=AA(2,K)                                                   
      AAA(K) = CHAR(K)*FI01(J)*QION*COULF                               
 374  CONTINUE                                                          
CX20532 CALL DOOUTX( 532 )                                              
CY20532 CALL DOOUTY( 532 )                                              
      QOXY=-0.6802D0                                                    
      QHYD= 0.3401D0                                                    
      B(1)=-HKS*PQDATA(5)*RF**6/UNITE                                   
      B(2)=-HKS*PQDATA(2)*RF**6/UNITE                                   
      B(3)=B(2)                                                         
      C(1)= HKS*PQDATA(6)*RF**12/UNITE                                  
      C(2)= HKS*PQDATA(3)*RF**12/UNITE                                  
      C(3)=C(2)                                                         
      A(1)= QION*PQDATA(4)*QOXY*COULF                                   
      A(2)= QION*PQDATA(1)*QHYD*COULF                                   
      A(3)=A(2)                                                         
C**************************                                             
C*       THIS IS TO TAKE ACCOUNT THE COORDINATES W.R.T.                 
C*                THE CENTER OF THE BOX.                                
      ZSHIFT=67.96D0                                                    
C**************************                                             
      IF(NION.LE.0) GO TO 525                                           
C------------------------------------                                   
      VI12=159889.2497685301                                            
      VI06=-2548.3840453430                                             
      VI01=1.0                                                          
C------------------------------------                                   
      BBBB=    -HKS*VI06*RF**6/ UNITE                                   
      CCCC=     HKS*VI12*RF**12/UNITE                                   
      AAAA= QION*QION*VI01*COULF                                        
      LSPA=.FALSE.                                                      
      IF(LRES) THEN                                                     
CX40014 CALL SBOUTX (14)                                                
CY40014 CALL SBOUTY (14)                                                
      RETURN                                                            
      ENDIF                                                             
      I=NSPECI(1)                                                       
CX10533 CALL DOINX( 533 )                                               
CY10533 CALL DOINY( 533 )                                               
      DO 520 J=1,NION                                                   
      I=I+1                                                             
      X0(I)=AXION(J)*RF                                                 
      Y0(I)=AYION(J)*RF                                                 
      Z0(I)=(AZION(J)-ZSHIFT)*RF                                        
  520 CONTINUE                                                          
CX20534 CALL DOOUTX( 534 )                                              
CY20534 CALL DOOUTY( 534 )                                              
  525 CONTINUE                                                          
C..............................................                         
C..............................................                         
CX10535 CALL DOINX( 535 )                                               
CY10535 CALL DOINY( 535 )                                               
      DO 530 J=1,NAT                                                    
      I=I+1                                                             
      X0(I)=XAT(J)*RF                                                   
      Y0(I)=YAT(J)*RF                                                   
      Z0(I)=ZAT(J)*RF                                                   
C     Z0(I)=(ZAT(J)-ZSHIFT)*RF                                          
  530 CONTINUE                                                          
CX20536 CALL DOOUTX( 536 )                                              
CY20536 CALL DOOUTY( 536 )                                              
      CALL SITES(X0,Y0,Z0,FX,FY,FZ,GX,GY,GZ,QA0,QB0,QC0,QD0             
     X          ,X6,Y6,Z6,SX,SY,SZ,FS1,FS2,FS3,N)                       
      WRITE(66,1110)                                                    
 1110 FORMAT(//1X,' **** ION   COORDINATES: ****'//)                    
      K=NSPECI(1)                                                       
CX10537 CALL DOINX( 537 )                                               
CY10537 CALL DOINY( 537 )                                               
      DO 1210 I=1,NION                                                  
      K=K+1                                                             
      X=X0(K)/RF                                                        
      Y=Y0(K)/RF                                                        
      Z=Z0(K)/RF                                                        
      WRITE(66,1300) K,X,Y,Z                                            
 1210 CONTINUE                                                          
CX20538 CALL DOOUTX( 538 )                                              
CY20538 CALL DOOUTY( 538 )                                              
      WRITE(66,1120)                                                    
 1120 FORMAT(//1X,' **** DNA   COORDINATES: ****'//)                    
CX10539 CALL DOINX( 539 )                                               
CY10539 CALL DOINY( 539 )                                               
      DO 1220 I=1,NAT                                                   
      K=K+1                                                             
      X=X0(K)/RF                                                        
      Y=Y0(K)/RF                                                        
      Z=Z0(K)/RF                                                        
      WRITE(66,1310) K,X,Y,Z,CHAR(I),ICLAS(I),JCLAS(I)                  
 1220 CONTINUE                                                          
CX20540 CALL DOOUTX( 540 )                                              
CY20540 CALL DOOUTY( 540 )                                              
 1300 FORMAT(1X,I6,3(1X,F12.8))                                         
 1310 FORMAT(1X,I6,3(1X,F12.8),1X,F12.6,2I5)                            
      NSP=NSPECI(1)                                                     
      NOP=NSP+NION                                                      
CX10541 CALL DOINX( 541 )                                               
CY10541 CALL DOINY( 541 )                                               
      DO 1500 I=1,NOP                                                   
      WX0(I)=0.D0                                                       
      WY0(I)=0.D0                                                       
      WZ0(I)=0.D0                                                       
 1500 CONTINUE                                                          
CX20542 CALL DOOUTX( 542 )                                              
CY20542 CALL DOOUTY( 542 )                                              
CX10543 CALL DOINX( 543 )                                               
CY10543 CALL DOINY( 543 )                                               
      DO 1600 I=1,NSP                                                   
      AX0(I)=0.D0                                                       
      AY0(I)=0.D0                                                       
      AZ0(I)=0.D0                                                       
 1600 CONTINUE                                                          
CX20544 CALL DOOUTX( 544 )                                              
CY20544 CALL DOOUTY( 544 )                                              
      LSPA=.TRUE.                                                       
      NSTEP=1                                                           
CX40014 CALL SBOUTX (14)                                                
CY40014 CALL SBOUTY (14)                                                
      RETURN                                                            
      END                                                               
C================= SITES K.N.S. 4/22 84 === REV.000  X.X.  X/X XX ======
C                                                                       
      SUBROUTINE SITES(X0,Y0,Z0,FX,FY,FZ,GX,GY,GZ,QA0,QB0,QC0,QD0       
     X                ,X6,Y6,Z6,SX,SY,SZ,FS1,FS2,FS3,NNN)               
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION FX(500),FY(500),FZ(500)                                 
      DIMENSION GX(500),GY(500),GZ(500)                                 
      DIMENSION X0(*),Y0(*),Z0(*)                                       
      DIMENSION X6(500),Y6(500),Z6(500)                                 
      DIMENSION QA0(*),QB0(*),QC0(*),QD0(*)                             
      DIMENSION SX(*),SY(*),SZ(*),FS1(*),FS2(*),FS3(*)                  
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON /CONTRL/ NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP   
     X               ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION       
     X               ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS               
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
      DIMENSION X3(500),Y3(500),Z3(500)                                 
      DIMENSION X4(500),Y4(500),Z4(500)                                 
      DIMENSION X8(500),Y8(500),Z8(500)                                 
      DIMENSION AX(500),AY(500),AZ(500)                                 
CX30015 CALL SBINX (15)                                                 
CY30015 CALL SBINY (15)                                                 
      RF=1.D-10/UNITL                                                   
      NSTEP=1                                                           
C      PI=4.D0*ATAN(1.D0)                                               
      ROH=0.957193D0                                                    
      D2Q=0.2677D0                                                      
      ELBOW=1.82422D0/2.D0                                              
      SHOH=SIN(ELBOW)                                                   
      CHOH=COS(ELBOW)                                                   
      RSOH = ROH*SHOH                                                   
      RCOH = ROH*CHOH                                                   
CX10545 CALL DOINX( 545 )                                               
CY10545 CALL DOINY( 545 )                                               
      DO 800 K=1,NNN                                                    
      Z6(K)=Z6(K)-67.96D0                                               
      XXX = X6(K)                                                       
      YYY = Y6(K)                                                       
      ZZZ = Z6(K)                                                       
      CX = RCOH*FX(K)                                                   
      CY = RCOH*FY(K)                                                   
      CZ = RCOH*FZ(K)                                                   
      BX = RSOH*GX(K)                                                   
      BY = RSOH*GY(K)                                                   
      BZ = RSOH*GZ(K)                                                   
      X3(K) = XXX + CX - BX                                             
      Y3(K) = YYY + CY - BY                                             
      Z3(K) = ZZZ + CZ - BZ                                             
      X4(K) = XXX + CX + BX                                             
      Y4(K) = YYY + CY + BY                                             
      Z4(K) = ZZZ + CZ + BZ                                             
      X8(K) = XXX + D2Q * FX(K)                                         
      Y8(K) = YYY + D2Q * FY(K)                                         
      Z8(K) = ZZZ + D2Q * FZ(K)                                         
  800 CONTINUE                                                          
CX20546 CALL DOOUTX( 546 )                                              
CY20546 CALL DOOUTY( 546 )                                              
CX10547 CALL DOINX( 547 )                                               
CY10547 CALL DOINY( 547 )                                               
      DO 100 I=1,NNN                                                    
      X3(I)=X3(I)-X6(I)                                                 
      Y3(I)=Y3(I)-Y6(I)                                                 
      Z3(I)=Z3(I)-Z6(I)                                                 
      X4(I)=X4(I)-X6(I)                                                 
      Y4(I)=Y4(I)-Y6(I)                                                 
      Z4(I)=Z4(I)-Z6(I)                                                 
      X=X8(I)-X6(I)                                                     
      Y=Y8(I)-Y6(I)                                                     
      Z=Z8(I)-Z6(I)                                                     
      PCNORM=1.D0/SQRT(X**2+Y**2+Z**2)                                  
      XX=X*PCNORM                                                       
      YY=Y*PCNORM                                                       
      ZZ=Z*PCNORM                                                       
      AX(I)=-0.06556D0*XX                                               
      AY(I)=-0.06556D0*YY                                               
      AZ(I)=-0.06556D0*ZZ                                               
      X6(I)=X6(I)-AX(I)                                                 
      Y6(I)=Y6(I)-AY(I)                                                 
      Z6(I)=Z6(I)-AZ(I)                                                 
      X3(I)=X3(I)+AX(I)                                                 
      Y3(I)=Y3(I)+AY(I)                                                 
      Z3(I)=Z3(I)+AZ(I)                                                 
      X4(I)=X4(I)+AX(I)                                                 
      Y4(I)=Y4(I)+AY(I)                                                 
      Z4(I)=Z4(I)+AZ(I)                                                 
      X8(I)=X+AX(I)                                                     
      Y8(I)=Y+AY(I)                                                     
      Z8(I)=Z+AZ(I)                                                     
  100 CONTINUE                                                          
CX20548 CALL DOOUTX( 548 )                                              
CY20548 CALL DOOUTY( 548 )                                              
      NS=NSITES(1)                                                      
CX10549 CALL DOINX( 549 )                                               
CY10549 CALL DOINY( 549 )                                               
      DO 950 I=1,NNN                                                    
      X0(I)=X6(I)*RF                                                    
      Y0(I)=Y6(I)*RF                                                    
      Z0(I)=Z6(I)*RF                                                    
  950 CONTINUE                                                          
CX20550 CALL DOOUTX( 550 )                                              
CY20550 CALL DOOUTY( 550 )                                              
      NSP=NSPECI(1)                                                     
      CALL QUATER(FX,FY,FZ,GX,GY,GZ,QA0,QB0,QC0,QD0,NNN)                
      CALL UPDATE(QA0,QB0,QC0,QD0,FS1,FS2,FS3,SX,SY,SZ)                 
CX10551 CALL DOINX( 551 )                                               
CY10551 CALL DOINY( 551 )                                               
      DO 960 I=1,NNN                                                    
      ISP=(I-1)*NS                                                      
      SX(ISP+1)=AX(I)*RF                                                
      SY(ISP+1)=AY(I)*RF                                                
      SZ(ISP+1)=AZ(I)*RF                                                
      SX(ISP+2)=X3(I)*RF                                                
      SY(ISP+2)=Y3(I)*RF                                                
      SZ(ISP+2)=Z3(I)*RF                                                
      SX(ISP+3)=X4(I)*RF                                                
      SY(ISP+3)=Y4(I)*RF                                                
      SZ(ISP+3)=Z4(I)*RF                                                
      SX(ISP+4)=X8(I)*RF                                                
      SY(ISP+4)=Y8(I)*RF                                                
      SZ(ISP+4)=Z8(I)*RF                                                
  960 CONTINUE                                                          
CX20552 CALL DOOUTX( 552 )                                              
CY20552 CALL DOOUTY( 552 )                                              
CX40016 CALL SBOUTX (16)                                                
CY40016 CALL SBOUTY (16)                                                
      RETURN                                                            
      END                                                               
C================= FCHECKK.N.S. 7/14 85 === REV.000  X.X.  X/X XX ======
C                                                                       
      SUBROUTINE FCHECK (FW,FI,FD)                                      
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION FW(*),FI(*),FD(*)                                       
      COMMON /CONTRL/ NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP   
     X               ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION       
     X               ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS               
CX30017 CALL SBINX (17)                                                 
CY30017 CALL SBINY (17)                                                 
      IF(MOD(NSTEP,10).NE.0) THEN                                       
CX40018 CALL SBOUTX (18)                                                
CY40018 CALL SBOUTY (18)                                                
      RETURN                                                            
      ENDIF                                                             
      FSUW=0.0                                                          
CX10553 CALL DOINX( 553 )                                               
CY10553 CALL DOINY( 553 )                                               
      DO 100 I=1,NSS                                                    
  100 FSUW=FSUW+FW(I)+FI(I)+FD(I)                                       
CX20554 CALL DOOUTX( 554 )                                              
CY20554 CALL DOOUTY( 554 )                                              
      FSUI=0.0                                                          
      N=NSS+NION                                                        
CX10555 CALL DOINX( 555 )                                               
CY10555 CALL DOINY( 555 )                                               
      DO 200 I=NSS+1,N                                                  
  200 FSUI=FSUI+FW(I)+FI(I)+FD(I)                                       
CX20556 CALL DOOUTX( 556 )                                              
CY20556 CALL DOOUTY( 556 )                                              
      FSUD=0.0                                                          
      NSST=NSS+NION+NATOMS                                              
CX10557 CALL DOINX( 557 )                                               
CY10557 CALL DOINY( 557 )                                               
      DO 300 I=N+1,NSST                                                 
  300 FSUD=FSUD+FW(I)+FI(I)+FD(I)                                       
CX20558 CALL DOOUTX( 558 )                                              
CY20558 CALL DOOUTY( 558 )                                              
      FSUT=FSUW+FSUI+FSUD                                               
      WRITE(66,400) FSUT,FSUW,FSUI,FSUD                                 
  400 FORMAT(1X,'*** FORCE SUMS:  TOTAL= ',D12.4,'  WATER=',D12.4,'     
     X ION=  ',D12.4,'   DNA= ',D12.4)                                  
CX40018 CALL SBOUTX (18)                                                
CY40018 CALL SBOUTY (18)                                                
      RETURN                                                            
      END                                                               
C================= QUATERK.N.S. 7/7 85 === REV.000  X.X.  X/X XX ====== 
C                                                                       
      SUBROUTINE QUATER (FX,FY,FZ,GX,GY,GZ,QA0,QB0,QC0,QD0,NNN)         
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION FX(*),FY(*),FZ(*),GX(*),GY(*),GZ(*)                     
      DIMENSION QA0(*),QB0(*),QC0(*),QD0(*)                             
      DIMENSION A1(500),A2(500),A3(500)                                 
      COMMON /CONTRL/ NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP   
     X               ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION       
     X               ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS               
CX30019 CALL SBINX (19)                                                 
CY30019 CALL SBINY (19)                                                 
      PI=4.D0*ATAN(1.D0)                                                
CX10559 CALL DOINX( 559 )                                               
CY10559 CALL DOINY( 559 )                                               
      DO 100 I=1,NNN                                                    
      A2(I)=ACOS(-FZ(I))                                                
      SIN2=SIN(A2(I))                                                   
      A3I=-FY(I)/SIN2                                                   
      IF(ABS(A3I).GT.1) A3I=A3I/ABS(A3I)                                
      A3(I)=ASIN(A3I)                                                   
      IF(COS(A3(I))*SIN2*FX(I).LT.0.D0) A3(I)=PI-A3(I)                  
      A1I=GZ(I)/SIN2                                                    
      IF(ABS(A1I).GT.1) A1I=A1I/ABS(A1I)                                
      A1(I)=ASIN(A1I)                                                   
      ABC=A1(I)                                                         
      DEF=A3(I)                                                         
      A3(I)=ABC                                                         
      A1(I)=DEF                                                         
      QA0(I)=-SIN((A3(I)-A1(I))/2.D0)*SIN(A2(I)/2.D0)                   
      QB0(I)= COS((A3(I)-A1(I))/2.D0)*SIN(A2(I)/2.D0)                   
      QC0(I)= SIN((A3(I)+A1(I))/2.D0)*COS(A2(I)/2.D0)                   
      QD0(I)= COS((A3(I)+A1(I))/2.D0)*COS(A2(I)/2.D0)                   
      QNORM=SQRT(QA0(I)**2+QB0(I)**2+QC0(I)**2+QD0(I)**2)               
      QA0(I)=QA0(I)/QNORM                                               
      QB0(I)=QB0(I)/QNORM                                               
      QC0(I)=QC0(I)/QNORM                                               
      QD0(I)=QD0(I)/QNORM                                               
  100 CONTINUE                                                          
CX20560 CALL DOOUTX( 560 )                                              
CY20560 CALL DOOUTY( 560 )                                              
CX40020 CALL SBOUTX (20)                                                
CY40020 CALL SBOUTY (20)                                                
      RETURN                                                            
      END                                                               
C================= INDEX K.N.S. 3/21 84 === REV.001K.N.S. 4/27 84 ======
C                                                                       
      SUBROUTINE INDEX1                                                 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/INDEXQ/ NWI,NWJ,NWK                                        
      COMMON/WINDEX/ NOR,IQ, IW, IA,IM,IF,IT,IP,IO,IX(99)               
      COMMON/WJNDEX/ JF1,JX(9)                                          
      COMMON/WKNDEX/ KX(9)                                              
CX30021 CALL SBINX (21)                                                 
CY30021 CALL SBINY (21)                                                 
      NSP=NSPECI(1)                                                     
      NS =NSITES(1)                                                     
      NSS3=NS*NSP*3                                                     
      NR=0                                                              
      NDIM=3                                                            
      NOP=NSP+NION                                                      
      IF(LLF) NR=22                                                     
      IF(LPC) NR=76                                                     
      IF(LLF.EQV.LPC) GO TO 940                                         
      IF(NR.EQ.0) GO TO 950                                             
      NTOT=NOP+NATOMS                                                   
      IX(1)=1                                                           
      IX(2)=IX(1)+NTOT                                                  
      IX(3)=IX(2)+NTOT                                                  
      IX(4)=IX(3)+NTOT                                                  
      IF(.NOT.LLF) GO TO 200                                            
CX10561 CALL DOINX( 561 )                                               
CY10561 CALL DOINY( 561 )                                               
      DO 100 I=5,8                                                      
  100 IX(I)=IX(I-1)+NSP                                                 
CX20562 CALL DOOUTX( 562 )                                              
CY20562 CALL DOOUTY( 562 )                                              
CX10563 CALL DOINX( 563 )                                               
CY10563 CALL DOINY( 563 )                                               
      DO 110 I=9,11                                                     
  110 IX(I)=IX(I-1)+NOP                                                 
CX20564 CALL DOOUTX( 564 )                                              
CY20564 CALL DOOUTY( 564 )                                              
CX10565 CALL DOINX( 565 )                                               
CY10565 CALL DOINY( 565 )                                               
      DO 120 I=12,14                                                    
  120 IX(I)=IX(I-1)+NSP                                                 
CX20566 CALL DOOUTX( 566 )                                              
CY20566 CALL DOOUTY( 566 )                                              
CX10567 CALL DOINX( 567 )                                               
CY10567 CALL DOINY( 567 )                                               
      DO 130 I=15,17                                                    
  130 IX(I)=IX(I-1)+NOP                                                 
CX20568 CALL DOOUTX( 568 )                                              
CY20568 CALL DOOUTY( 568 )                                              
CX10569 CALL DOINX( 569 )                                               
CY10569 CALL DOINY( 569 )                                               
      DO 140 I=18,20                                                    
  140 IX(I)=IX(I-1)+NOP                                                 
CX20570 CALL DOOUTX( 570 )                                              
CY20570 CALL DOOUTY( 570 )                                              
CX10571 CALL DOINX( 571 )                                               
CY10571 CALL DOINY( 571 )                                               
      DO 150 I=21,22                                                    
  150 IX(I)=IX(I-1)+NSP                                                 
CX20572 CALL DOOUTX( 572 )                                              
CY20572 CALL DOOUTY( 572 )                                              
      LIMIT=IX(NR)+NSP-1                                                
      IF(LIMIT.GT.NWI) GO TO 900                                        
      IQ=4                                                              
      IW=8                                                              
      IA=11                                                             
      IM=14                                                             
      IF=17                                                             
      IT=20                                                             
  200 IF(.NOT.LPC) GO TO 300                                            
CX10573 CALL DOINX( 573 )                                               
CY10573 CALL DOINY( 573 )                                               
      DO 210 I=5,19                                                     
  210 IX(I)=IX(I-1)+NOP                                                 
CX20574 CALL DOOUTX( 574 )                                              
CY20574 CALL DOOUTY( 574 )                                              
CX10575 CALL DOINX( 575 )                                               
CY10575 CALL DOINY( 575 )                                               
      DO 220 I=20,47                                                    
  220 IX(I)=IX(I-1)+NSP                                                 
CX20576 CALL DOOUTX( 576 )                                              
CY20576 CALL DOOUTY( 576 )                                              
CX10577 CALL DOINX( 577 )                                               
CY10577 CALL DOINY( 577 )                                               
      DO 230 I=48,68                                                    
  230 IX(I)=IX(I-1)+NSP                                                 
CX20578 CALL DOOUTX( 578 )                                              
CY20578 CALL DOOUTY( 578 )                                              
CX10579 CALL DOINX( 579 )                                               
CY10579 CALL DOINY( 579 )                                               
      DO 240 I=69,71                                                    
  240 IX(I)=IX(I-1)+NOP                                                 
CX20580 CALL DOOUTX( 580 )                                              
CY20580 CALL DOOUTY( 580 )                                              
CX10581 CALL DOINX( 581 )                                               
CY10581 CALL DOINY( 581 )                                               
      DO 250 I=72,74                                                    
  250 IX(I)=IX(I-1)+NOP                                                 
CX20582 CALL DOOUTX( 582 )                                              
CY20582 CALL DOOUTY( 582 )                                              
CX10583 CALL DOINX( 583 )                                               
CY10583 CALL DOINY( 583 )                                               
      DO 260 I=75,76                                                    
  260 IX(I)=IX(I-1)+NSP                                                 
CX20584 CALL DOOUTX( 584 )                                              
CY20584 CALL DOOUTY( 584 )                                              
      LIMIT=IX(NR)+NSP-1                                                
      IF(LIMIT.GT.NWI) GO TO 900                                        
      IQ=19                                                             
      IW=4                                                              
      IA=47                                                             
      IM=68                                                             
      IF=71                                                             
      IT=74                                                             
  300 CONTINUE                                                          
      NSP=NSPECI(1)                                                     
      NS=NSITES(1)                                                      
      NSS=NS*NSP                                                        
      NSST=NSS+NION+NATOMS                                              
      NOP3=3*NOP                                                        
      NOP4=4*NOP                                                        
      NSS3=3*NSS                                                        
      MR=4                                                              
      NR=6                                                              
      JX(1)=1                                                           
CX10585 CALL DOINX( 585 )                                               
CY10585 CALL DOINY( 585 )                                               
      DO 500 J=2,MR                                                     
  500 JX(J)=JX(J-1)+NSS                                                 
CX20586 CALL DOOUTX( 586 )                                              
CY20586 CALL DOOUTY( 586 )                                              
CX10587 CALL DOINX( 587 )                                               
CY10587 CALL DOINY( 587 )                                               
      DO 510 J=5,NR                                                     
  510 JX(J)=JX(J-1)+NSST                                                
CX20588 CALL DOOUTX( 588 )                                              
CY20588 CALL DOOUTY( 588 )                                              
      LJMJT=JX(NR)+NSST-1                                               
      IF(LJMJT.GT.NWJ) GO TO 910                                        
      JF1=4                                                             
      KX(1)=1                                                           
      NR=3                                                              
CX10589 CALL DOINX( 589 )                                               
CY10589 CALL DOINY( 589 )                                               
      DO 550 J=2,NR                                                     
  550 KX(J)=KX(J-1)+NENT                                                
CX20590 CALL DOOUTX( 590 )                                              
CY20590 CALL DOOUTY( 590 )                                              
      LKMKT=KX(NR)+NENT-1                                               
      IF(LKMKT.GT.NWK) GO TO 920                                        
      WRITE(66,600) LIMIT,NWI                                           
  600 FORMAT(/1X,'*** NR OF ELEMENTS IN C.O.M. WORK ARRAY = ',I8,'      
     XCURRENT DIMENSION = ',I8)                                         
      WRITE(66,610) LJMJT,NWJ                                           
  610 FORMAT(/1X,'*** NR OF ELEMENTS IN  SITE  WORK ARRAY = ',I8,'      
     XCURRENT DIMENSION = ',I8)                                         
      WRITE(66,620) LKMKT,NWK                                           
  620 FORMAT(/1X,'*** NR OF ELEMENTS IN BUFFER WORK ARRAY = ',I8,'      
     XCURRENT DIMENSION = ',I8)                                         
      WRITE(66,630)                                                     
  630 FORMAT(/1X,'??? JUST IN CASE  -  CHECK THAT THE WORK ARRAY DIMENSI
     XONS IN BLOCK DATA ARE THE SAME AS IN THE MAIN PROGRAM'//)         
CX40022 CALL SBOUTX (22)                                                
CY40022 CALL SBOUTY (22)                                                
      RETURN                                                            
  900 WRITE(66,975) LIMIT,NWI                                           
CX40022 CALL SBOUTX (22)                                                
CY40022 CALL SBOUTY (22)                                                
      STOP                                                              
  910 WRITE(66,980) LJMJT,NWJ                                           
CX40022 CALL SBOUTX (22)                                                
CY40022 CALL SBOUTY (22)                                                
      STOP                                                              
  920 WRITE(66,985) LKMKT,NWK                                           
CX40022 CALL SBOUTX (22)                                                
CY40022 CALL SBOUTY (22)                                                
      STOP                                                              
  940 WRITE(66,990)                                                     
CX40022 CALL SBOUTX (22)                                                
CY40022 CALL SBOUTY (22)                                                
      STOP                                                              
  950 WRITE(66,995)                                                     
CX40022 CALL SBOUTX (22)                                                
CY40022 CALL SBOUTY (22)                                                
      STOP                                                              
  975 FORMAT(1X,']]] NR OF ELEMENTS IN C.O.M. WORK ARRAY ',I8,' EXCEEDS 
     X THE CURRENT DIMENSION',I8,'  - INCREASE THE LIMIT')              
  980 FORMAT(1X,']]] NR OF ELEMENTS IN  SITE  WORK ARRAY ',I8,' EXCEEDS 
     XTHE CURRENT DIMENSION',I8,'  - INCREASE THE LIMIT')               
  985 FORMAT(1X,']]] NR OF ELEMENTS IN BUFFER WORK ARRAY ',I8,' EXCEEDS 
     XTHE CURRENT DIMENSION',I8,'  - INCREASE THE LIMIT')               
  990 FORMAT(1X,']]] INTEGRATOR TYPE NOT SPECIFIED')                    
  995 FORMAT(1X,']]] MC/MD MODE NOT SPECIFIED')                         
CX40022 CALL SBOUTX (22)                                                
CY40022 CALL SBOUTY (22)                                                
      END                                                               
C================= MOTIONK.N.S. 3/31 84 === REV.001K.N.S. 4/27 84 ======
      SUBROUTINE MOTION(X0,Y0,Z0,QA0,QB0,QC0,QD0,WX0,WY0,WZ0            
     X ,AX0,AY0,AZ0,XM,YM,ZM,FX,FY,FZ,TX,TY,TZ,QA,QB,QC,QD)             
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION  X0(*), Y0(*), Z0(*),QA0(*),QB0(*),QC0(*),QD0(*)        
      DIMENSION WX0(*),WY0(*),WZ0(*),AX0(*),AY0(*),AZ0(*)               
      DIMENSION  XM(*), YM(*), ZM(*), FX(*), FY(*), FZ(*)               
      DIMENSION  TX(*), TY(*), TZ(*), QA(*), QB(*), QC(*), QD(*)        
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONSTS/ AVSNO,BOLTZ,PLANCK,PI,ELECHG,EPS0                  
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
CX30023 CALL SBINX (23)                                                 
CY30023 CALL SBINY (23)                                                 
      CALL ZERO(NSS,QA(1))                                              
      CALL ZERO(NSS,QB(1))                                              
      CALL ZERO(NSS,QC(1))                                              
      CALL ZERO(NSS,QD(1))                                              
      ALENGT=22.D-10/UNITL                                              
CX10591 CALL DOINX( 591 )                                               
CY10591 CALL DOINY( 591 )                                               
      DO 10 I=1,NTYPES                                                  
      TE(I)=0.D0                                                        
   10 CONTINUE                                                          
CX20592 CALL DOOUTX( 592 )                                              
CY20592 CALL DOOUTY( 592 )                                              
      I=0                                                               
CX10593 CALL DOINX( 593 )                                               
CY10593 CALL DOINY( 593 )                                               
      DO 500 N=1,NTYPES                                                 
      R2(N)=0.D0                                                        
      VSQ=0.D0                                                          
      NSP=NSPECI(N)                                                     
      FNSP=DBLE(NSP)                                                    
C      RFACT=FNSP*UNITL**2                                              
      R=1.D0/DSUMM(N)                                                   
      DO 100 J=1,NSP                                                    
      I=I+1                                                             
      VX=WX0(I)                                                         
      VY=WY0(I)                                                         
      VZ=WZ0(I)                                                         
      WX0(I)=WX0(I)+FX(I)                                               
      WY0(I)=WY0(I)+FY(I)                                               
      WZ0(I)=WZ0(I)+FZ(I)                                               
      VX=(VX+WX0(I))*0.5D0                                              
      VY=(VY+WY0(I))*0.5D0                                              
      VZ=(VZ+WZ0(I))*0.5D0                                              
      VSQ=VSQ+(VX**2+VY**2+VZ**2)                                       
      X0(I)=X0(I)+WX0(I)*R                                              
      Y0(I)=Y0(I)+WY0(I)*R                                              
      Z0(I)=Z0(I)+WZ0(I)*R                                              
      XM(I)=XM(I)+WX0(I)*R                                              
      YM(I)=YM(I)+WY0(I)*R                                              
      ZM(I)=ZM(I)+WZ0(I)*R                                              
      RM2=XM(I)**2+YM(I)**2+ZM(I)**2                                    
      R2(N)=R2(N)+RM2*UNITL**2/FNSP                                     
  100 CONTINUE                                                          
      TE(N)=0.5D0*R*VSQ                                                 
  500 CONTINUE                                                          
CX20594 CALL DOOUTX( 594 )                                              
CY20594 CALL DOOUTY( 594 )                                              
CX10595 CALL DOINX( 595 )                                               
CY10595 CALL DOINY( 595 )                                               
      DO 501 I=1,NOP                                                    
      X0(I)=X0(I)-2.D0*ALENGT*DBLE(INT(X0(I)/ALENGT))                   
      Y0(I)=Y0(I)-2.D0*ALENGT*DBLE(INT(Y0(I)/ALENGT))                   
      Z0(I)=Z0(I)-2.D0*DBLE(INT(Z0(I)))                                 
  501 CONTINUE                                                          
CX20596 CALL DOOUTX( 596 )                                              
CY20596 CALL DOOUTY( 596 )                                              
CX10597 CALL DOINX( 597 )                                               
CY10597 CALL DOINY( 597 )                                               
      DO 510 I=1,NTYPES                                                 
      RE(I)=0.D0                                                        
  510 CONTINUE                                                          
CX20598 CALL DOOUTX( 598 )                                              
CY20598 CALL DOOUTY( 598 )                                              
      I=0                                                               
CX10599 CALL DOINX( 599 )                                               
CY10599 CALL DOINY( 599 )                                               
      DO 1500 N=1,NTYPES                                                
      RKET=0.D0                                                         
      NSIT=NSITES(N)                                                    
      IF(NSIT.EQ.1) GO TO 1450                                          
      NSP=NSPECI(N)                                                     
      DX=DRI(1,N)                                                       
      DY=DRI(2,N)                                                       
      DZ=DRI(3,N)                                                       
      DO 1200 J=1,NSP                                                   
      I=I+1                                                             
      FX(I)=AX0(I)+TX(I)*0.5D0                                          
      FY(I)=AY0(I)+TY(I)*0.5D0                                          
      FZ(I)=AZ0(I)+TZ(I)*0.5D0                                          
 1200 CONTINUE                                                          
      CALL ROTATE(QA0,QB0,QC0,QD0,FX,FY,FZ,NSP,.TRUE.)                  
      I=I-NSP                                                           
      DO 1300 J=1,NSP                                                   
      I=I+1                                                             
      XS=FX(I)/DX                                                       
      YS=FY(I)/DY                                                       
      ZS=FZ(I)/DZ                                                       
      RKET=RKET+DX*XS*XS+DY*YS*YS+DZ*ZS*ZS                              
      QA(I)=QA0(I)+(-QC0(I)*XS-QD0(I)*YS+QB0(I)*ZS)*0.25D0              
      QB(I)=QB0(I)+(+QD0(I)*XS-QC0(I)*YS-QA0(I)*ZS)*0.25D0              
      QC(I)=QC0(I)+(+QA0(I)*XS+QB0(I)*YS+QD0(I)*ZS)*0.25D0              
      QD(I)=QD0(I)+(-QB0(I)*XS+QA0(I)*YS-QC0(I)*ZS)*0.25D0              
      AX0(I)=AX0(I)+TX(I)                                               
      AY0(I)=AY0(I)+TY(I)                                               
      AZ0(I)=AZ0(I)+TZ(I)                                               
      FX(I)=AX0(I)                                                      
      FY(I)=AY0(I)                                                      
      FZ(I)=AZ0(I)                                                      
 1300 CONTINUE                                                          
      CALL ROTATE(QA,QB,QC,QD,FX,FY,FZ,NSP,.TRUE.)                      
      I=I-NSP                                                           
      DO 1400 J=1,NSP                                                   
      I=I+1                                                             
      XS=FX(I)/DX                                                       
      YS=FY(I)/DY                                                       
      ZS=FZ(I)/DZ                                                       
      QA0(I)=QA0(I)+(-QC(I)*XS-QD(I)*YS+QB(I)*ZS)*0.50D0                
      QB0(I)=QB0(I)+(+QD(I)*XS-QC(I)*YS-QA(I)*ZS)*0.50D0                
      QC0(I)=QC0(I)+(+QA(I)*XS+QB(I)*YS+QD(I)*ZS)*0.50D0                
      QD0(I)=QD0(I)+(-QB(I)*XS+QA(I)*YS-QC(I)*ZS)*0.50D0                
      QNORM=SQRT(QA0(I)**2+QB0(I)**2+QC0(I)**2+QD0(I)**2)               
      QNORM=1.D0/QNORM                                                  
      QA0(I)=QNORM*QA0(I)                                               
      QB0(I)=QNORM*QB0(I)                                               
      QC0(I)=QNORM*QC0(I)                                               
      QD0(I)=QNORM*QD0(I)                                               
 1400 CONTINUE                                                          
 1450 CONTINUE                                                          
      RE(N)=RKET*0.5D0                                                  
 1500 CONTINUE                                                          
CX20600 CALL DOOUTX( 600 )                                              
CY20600 CALL DOOUTY( 600 )                                              
CX40024 CALL SBOUTX (24)                                                
CY40024 CALL SBOUTY (24)                                                
      RETURN                                                            
      END                                                               
C                                                                       
C================= FIXCOMK.N.S. 3/17 84 === REV.000  X.X.  X/X XX ======
C                                                                       
      SUBROUTINE FIXCOM(K,L,M,R,MOMENT,SUM,RMOI)                        
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DOUBLE PRECISION M,MOMENT                                         
      DIMENSION M(*),R(3,*),MOMENT(3,3),RMOI(3,3)                       
      DIMENSION      COM(3),T(3)                                        
CX30025 CALL SBINX (25)                                                 
CY30025 CALL SBINY (25)                                                 
      NDIM=3                                                            
CX10601 CALL DOINX( 601 )                                               
CY10601 CALL DOINY( 601 )                                               
      DO 10 J=1,3                                                       
   10 COM(J)=0.D0                                                       
CX20602 CALL DOOUTX( 602 )                                              
CY20602 CALL DOOUTY( 602 )                                              
      SUM=0.D0                                                          
CX10603 CALL DOINX( 603 )                                               
CY10603 CALL DOINY( 603 )                                               
      DO 20 I=L,K                                                       
      SUM=SUM+M(I)                                                      
      DO 20 J=1,3                                                       
   20 COM(J)=COM(J)+M(I)*R(J,I)                                         
CX20604 CALL DOOUTX( 604 )                                              
CY20604 CALL DOOUTY( 604 )                                              
CX10605 CALL DOINX( 605 )                                               
CY10605 CALL DOINY( 605 )                                               
      DO 30 J=1,3                                                       
   30 COM(J)=COM(J)/SUM                                                 
CX20606 CALL DOOUTX( 606 )                                              
CY20606 CALL DOOUTY( 606 )                                              
CX10607 CALL DOINX( 607 )                                               
CY10607 CALL DOINY( 607 )                                               
      DO 40 I=L,K                                                       
      DO 40 J=1,3                                                       
   40 R(J,I)=R(J,I)-COM(J)                                              
CX20608 CALL DOOUTX( 608 )                                              
CY20608 CALL DOOUTY( 608 )                                              
CX10609 CALL DOINX( 609 )                                               
CY10609 CALL DOINY( 609 )                                               
      DO 60 INERT=1,3                                                   
      DO 60 IA=1,3                                                      
      TT=0.D0                                                           
      MOMENT(INERT,IA)=0.D0                                             
      IF(INERT.EQ.IA) TT=1.D0                                           
      DO 60 I=L,K                                                       
      SQS=0.D0                                                          
      DO 50 J=1,3                                                       
   50 SQS=SQS+R(J,I)**2                                                 
   60 MOMENT(INERT,IA)=MOMENT(INERT,IA)+                                
     X M(I)*(TT*SQS-R(INERT,I)*R(IA,I))                                 
CX20610 CALL DOOUTX( 610 )                                              
CY20610 CALL DOOUTY( 610 )                                              
      CALL JACOBI(MOMENT,RMOI,NDIM)                                     
CX10611 CALL DOINX( 611 )                                               
CY10611 CALL DOINY( 611 )                                               
      DO 120 I=L,K                                                      
      DO 100 J=1,3                                                      
  100 T(J)=R(J,I)                                                       
      DO 120 J=1,3                                                      
  120 R(J,I)=T(1)*RMOI(1,J)+T(2)*RMOI(2,J)+T(3)*RMOI(3,J)               
CX20612 CALL DOOUTX( 612 )                                              
CY20612 CALL DOOUTY( 612 )                                              
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
      RETURN                                                            
      END                                                               
C================= COORDSK.N.S. 3/20 84 === REV.000  X.X.  X/X XX ======
C                                                                       
      SUBROUTINE COORDS(NA,R)                                           
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION R(3,*)                                                  
      COMMON/CONSTS/ AVSNO,BOLTZ,PLANCK,PI,ELECHG,EPS0                  
CX30027 CALL SBINX (27)                                                 
CY30027 CALL SBINY (27)                                                 
      IF(NA.EQ.1) GO TO 2000                                            
      ELBOW=1.82422D0                                                   
C      ELBG=(ELBOW*180.D0)/PI                                           
      ROH=0.0957193D0                                                   
      CHOH=COS(0.5D0*ELBOW)                                             
      SHOH=SIN(0.5D0*ELBOW)                                             
      OG=CHOH*ROH                                                       
      GH=SHOH*ROH                                                       
      R(1,1) = 0.D0                                                     
      R(2,1) = 0.D0                                                     
      R(3,1) = 0.D0                                                     
      R(1,2) = 0.D0                                                     
      R(2,2) = GH                                                       
      R(3,2) = -OG                                                      
      R(1,3) = 0.D0                                                     
      R(2,3) = -GH                                                      
      R(3,3) = -OG                                                      
      RON=0.02677D0                                                     
      R(1,4)=0.D0                                                       
      R(2,4)=0.D0                                                       
      R(3,4)=-RON                                                       
      ROM=0.436D0*0.052917D0                                            
      OGM=CHOH*ROM                                                      
      GHM=SHOH*ROM                                                      
      R(1,5)=0.D0                                                       
      R(2,5)=0.D0                                                       
      R(3,5)=-ROM                                                       
      R(1,6)=0.D0                                                       
      R(2,6)=GHM                                                        
      R(3,6)=-OGM                                                       
      R(1,7)=0.D0                                                       
      R(2,7)=-GHM                                                       
      R(3,7)=-OGM                                                       
CX40028 CALL SBOUTX (28)                                                
CY40028 CALL SBOUTY (28)                                                
      RETURN                                                            
 2000 CONTINUE                                                          
      R(1,1)=0.D0                                                       
      R(2,1)=0.D0                                                       
      R(3,1)=0.D0                                                       
CX40028 CALL SBOUTX (28)                                                
CY40028 CALL SBOUTY (28)                                                
      RETURN                                                            
      END                                                               
C================= UNITS  K.N.S. 3/17 84 === REV.001K.N.S. 4/26 84 ==== 
C                                                                       
      SUBROUTINE UNITS(SITE,Q,AMASS,SUMM,RRI)                           
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION SITE(3,*),Q(8),AMASS(*),SUMM(*),RRI(3,3)                
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
      COMMON/CONSTS/ AVSNO,BOLTZ,PLANCK,PI,ELECHG,EPS0                  
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/TFACTS/ FACTE(3),FACTL(2)                                  
      COMMON/CORRCC/F02,F12,F22,F32,F42,F52,F62                         
     X             ,F01,F11,F21,F31,F41,F51,F61                         
     X             ,F03,F13,F23,F33,F43,F53,F63                         
CX30029 CALL SBINX (29)                                                 
CY30029 CALL SBINY (29)                                                 
C      DGR=20.D-12                                                      
      F02=3.D0/16.D0                                                    
      F12=251.D0/360.D0                                                 
      F22=1.D0                                                          
      F32=11.D0/18.D0                                                   
      F42=1.D0/6.D0                                                     
      F52=1.D0/60.D0                                                    
      F62=1.D0                                                          
      F01=19087.D0/60480.D0                                             
      F11=1.D0                                                          
      F21=137.D0/120.D0                                                 
      F31=5.D0/8.D0                                                     
      F41=17.D0/96.D0                                                   
      F51=1.D0/40.D0                                                    
      F61=1.D0/720.D0                                                   
      F03=863.D0/6048.D0                                                
      F13=665.D0/1008.D0                                                
      F23=1.D0                                                          
      F33=25.D0/36.D0                                                   
      F43=35.D0/144.D0                                                  
      F53=1.D0/24.D0                                                    
      F63=1.D0/360.D0                                                   
C********************************************************************   
      FACTM=AVSNO*1000.D0                                               
CX10613 CALL DOINX( 613 )                                               
CY10613 CALL DOINY( 613 )                                               
      DO 100 I=1,NTYPES                                                 
      IF(SUMM(I).GT.1.D-20) SUMM(I)=SUMM(I)/FACTM                       
      T=SQRT(RRI(1,I)**2+RRI(2,I)**2+RRI(3,I)**2)                       
      IF(T.LT.1.D-30) GO TO 70                                          
      DO 50 J=1,3                                                       
   50 RRI(J,I)=RRI(J,I)/FACTM                                           
   70 CONTINUE                                                          
  100 CONTINUE                                                          
CX20614 CALL DOOUTX( 614 )                                              
CY20614 CALL DOOUTY( 614 )                                              
CX10615 CALL DOINX( 615 )                                               
CY10615 CALL DOINY( 615 )                                               
      DO 200 I=1,NSITET                                                 
      IF(AMASS(I).GT.1.D-20) AMASS(I)=AMASS(I)/FACTM                    
      IF(ABS(Q(I)).GT.1.D-15) Q(I)=Q(I)*ELECHG                          
      T=SQRT(SITE(1,I)**2+SITE(2,I)**2+SITE(3,I)**2)                    
      IF(T.LT.1.D-6) GO TO 170                                          
      DO 150 J=1,3                                                      
  150 SITE(J,I)=SITE(J,I)*1.D-9                                         
  170 CONTINUE                                                          
  200 CONTINUE                                                          
CX20616 CALL DOOUTX( 616 )                                              
CY20616 CALL DOOUTY( 616 )                                              
      CALL MPOLES(SITE,Q)                                               
      TOTMAS=0.D0                                                       
CX10617 CALL DOINX( 617 )                                               
CY10617 CALL DOINY( 617 )                                               
      DO 300 I=1,NTYPES                                                 
  300 TOTMAS=TOTMAS+SUMM(I)*DBLE(NSPECI(I))                             
CX20618 CALL DOOUTX( 618 )                                              
CY20618 CALL DOOUTY( 618 )                                              
      IF(RHO.NE.0.D0) VOLM=(AVSNO/FNOP)*                                
     X (TOTMAS/(RHO*1000.D0))                                           
      VOLUME=VOLM*FNOP/AVSNO                                            
      UNITM=TOTMAS/FNOP                                                 
C     UNITL=0.5D0*VOLUME**(1.D0/3.D0)                                   
      UNITL=16.90D-10                                                   
      UNITT=DT                                                          
      UNITE=UNITM*(UNITL/UNITT)**2                                      
      IF(LPC) THEN                                                      
      UNITE=BOLTZ*TTRAN(1)                                              
      TSTEP=DT/SQRT((UNITM/UNITE)*UNITL**2)                             
      ENDIF                                                             
      COULF=(ELECHG/EPS0)*(ELECHG/UNITE)/(4.D0*PI*UNITL)                
      ENERF=(AVSNO/FNOP)*UNITE                                          
      TEMPF=2.D0/(3.D0*AVSNO*BOLTZ)                                     
CX10619 CALL DOINX( 619 )                                               
CY10619 CALL DOINY( 619 )                                               
      DO 400 I=1,NTYPES                                                 
      DSUMM(I)=SUMM(I)/UNITM                                            
      DO 350 J=1,3                                                      
  350 DRI(J,I)=(RRI(J,I)/UNITM)/((UNITL*1.D+9)*(UNITL*1.D+9))           
  400 CONTINUE                                                          
CX20620 CALL DOOUTX( 620 )                                              
CY20620 CALL DOOUTY( 620 )                                              
CX10621 CALL DOINX( 621 )                                               
CY10621 CALL DOINY( 621 )                                               
      DO 500 I=1,NSITET                                                 
      DQ(I)=Q(I)/ELECHG                                                 
      DO 450 J=1,3                                                      
  450 DSITE(J,I)=SITE(J,I)/UNITL                                        
  500 CONTINUE                                                          
CX20622 CALL DOOUTX( 622 )                                              
CY20622 CALL DOOUTY( 622 )                                              
      AUA=0.52917D0                                                     
      HKCAL=627.503D0                                                   
      FACTE(1)=4.184D0/6.022D20                                         
      FACTE(2)=HKCAL*1.D-3*FACTE(1)                                     
      FACTE(3)=HKCAL*FACTE(1)                                           
      FACTL(1)=1.D10                                                    
      FACTL(2)=FACTL(1)/AUA                                             
CX40030 CALL SBOUTX (30)                                                
CY40030 CALL SBOUTY (30)                                                
      RETURN                                                            
      END                                                               
C================= MPOLESK.N.S. 3/17 84 === REV.000  X.X.  X/X XX ======
C                                                                       
      SUBROUTINE MPOLES(SITE,Q)                                         
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION Q(*),SITE(3,*)                                          
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
CX30031 CALL SBINX (31)                                                 
CY30031 CALL SBINY (31)                                                 
      K=0                                                               
CX10623 CALL DOINX( 623 )                                               
CY10623 CALL DOINY( 623 )                                               
      DO 100 I=1,NTYPES                                                 
      NS=NSITES(I)                                                      
      IF(NS.EQ.1) GO TO 100                                             
      DO 70 IS=1,NS                                                     
      KS=0                                                              
      DO 50 J=1,3                                                       
      DM(J,I)=DM(J,I)+Q(K+IS)*SITE(J,K+IS)                              
      DO 50 L=1,3                                                       
      KS=KS+1                                                           
   50 QM(KS,I)=QM(KS,I)+Q(K+IS)*SITE(J,K+IS)*SITE(L,K+IS)               
   70 CONTINUE                                                          
      K=K+NS                                                            
  100 CONTINUE                                                          
CX20624 CALL DOOUTX( 624 )                                              
CY20624 CALL DOOUTY( 624 )                                              
CX40032 CALL SBOUTX (32)                                                
CY40032 CALL SBOUTY (32)                                                
      RETURN                                                            
      END                                                               
C================= WATER K.N.S. 5/30 85 === REV.000  X.X. X/XX XX ======
C                                                                       
      SUBROUTINE WATER                                                  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      CHARACTER*4  NEU,NLU,NT1,NT2                                      
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/INDEXQ/ NWI,NWJ,NWK,INDEX(11,11),INDS(3,10),NONE           
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/TABELS/ QQ(11,11),A(5),B(5)                                
      COMMON/TPARAM/ NCPBF(3),NPBF,NQ                                   
CX30033 CALL SBINX (33)                                                 
CY30033 CALL SBINY (33)                                                 
      N=0                                                               
CX10625 CALL DOINX( 625 )                                               
CY10625 CALL DOINY( 625 )                                               
      DO 10 I=1,NTYPES                                                  
      NS=NSITES(I)                                                      
      DO 10 J=1,NS                                                      
      N=N+1                                                             
      INDS(I,J)=N                                                       
   10 CONTINUE                                                          
CX20626 CALL DOOUTX( 626 )                                              
CY20626 CALL DOOUTY( 626 )                                              
      CALL PARPOT(NEU,NLU,NT1,NT2)                                      
      CALL CNVERT(NEU,NLU,NT1,NT2)                                      
CX40034 CALL SBOUTX (34)                                                
CY40034 CALL SBOUTY (34)                                                
      RETURN                                                            
      END                                                               
C================= PARPOTK.N.S. 3/24 84 === REV.001K.N.S. 4/26 84 ======
C                                                                       
      SUBROUTINE PARPOT(NEU,NLU,NT1,NT2)                                
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      CHARACTER*4 ZITES,TYPE,ESORT,LSORT,NT1,NT2,NEU,NLU                
      DIMENSION ZITES(9),TYPE(3),ESORT(3),LSORT(2)                      
      SAVE ZITES, TYPE, ESORT, LSORT                                    
      DIMENSION CHRG(11)                                                
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/INDEXQ/ NWI,NWJ,NWK,INDEX(11,11),INDS(3,10),NONE           
      COMMON/TPARAM/ NCPBF(3),NPBF,NQ                                   
      COMMON/TABELS/ QQ(11,11),A(5),B(5)                                
      COMMON/PARAMS/ A1,A2,A3,A4,A5,B1                                  
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/CONSTS/ AVSNO,BOLTZ,PLANCK,PI,ELECHG,EPS0                  
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
      DATA ZITES/'   O','  H1','  H2','   P','   M','  M1','  M2',      
     X'+ION','-ION'/                                                    
      DATA TYPE /' H2O','+ION','-ION'/                                  
      DATA ESORT /'KCAL','MHAR','HART'/                                 
      DATA LSORT /'ANGS','A.U.'/                                        
CX30035 CALL SBINX (35)                                                 
CY30035 CALL SBINY (35)                                                 
      A(1)=1088213.2D0                                                  
      A(2)=1455.427D0                                                   
      A(3)=-273.5954D0                                                  
      A(4)=666.3373D0                                                   
      A(5)=0.D0                                                         
      B(1)=5.152712D0                                                   
      B(2)=2.961895D0                                                   
      B(3)=2.233264D0                                                   
      B(4)=2.760844D0                                                   
      B(5)=0.D0                                                         
      NT1= TYPE(1)                                                      
      NT2=TYPE(1)                                                       
      NPBF=3                                                            
      NCPBF(1)=1                                                        
      NCPBF(2)=2                                                        
      NCPBF(3)=1                                                        
      NEU=ESORT(1)                                                      
      NLU=LSORT(1)                                                      
      CHRG(1)=0.D0                                                      
      CHRG(2)=0.7175D0                                                  
      CHRG(3)=0.7175D0                                                  
      CHRG(4)=-1.435D0                                                  
CX10627 CALL DOINX( 627 )                                               
CY10627 CALL DOINY( 627 )                                               
      DO 100 I=1,11                                                     
      DO 100 J=1,11                                                     
  100 INDEX(I,J)=4                                                      
CX20628 CALL DOOUTX( 628 )                                              
CY20628 CALL DOOUTY( 628 )                                              
      INDEX(1,1)=1                                                      
      INDEX(1,2)=2                                                      
      INDEX(1,3)=2                                                      
      INDEX(2,2)=3                                                      
      INDEX(2,3)=3                                                      
      INDEX(3,3)=3                                                      
CX10629 CALL DOINX( 629 )                                               
CY10629 CALL DOINY( 629 )                                               
      DO 110 I=1,3                                                      
      DO 110 J=I,3                                                      
  110 INDEX(J,I)=INDEX(I,J)                                             
CX20630 CALL DOOUTX( 630 )                                              
CY20630 CALL DOOUTY( 630 )                                              
CX10631 CALL DOINX( 631 )                                               
CY10631 CALL DOINY( 631 )                                               
      DO 120 I=1,4                                                      
      DO 120 J=1,4                                                      
      QQ(I,J)=CHRG(I)*CHRG(J)                                           
  120 CONTINUE                                                          
CX20632 CALL DOOUTX( 632 )                                              
CY20632 CALL DOOUTY( 632 )                                              
      M=4                                                               
      WRITE(66,1505)                                                    
      WRITE(66,1510)                                                    
 1505 FORMAT(////20X,' P O T E N T I A L   M O D E L       : '//)       
 1510 FORMAT(/20X,' NON-COULOMB INTERACTION INDEX TABLE: '//)           
      WRITE(66,1520) (I,I=1,M)                                          
 1520 FORMAT(/20X,3X,10I5//)                                            
CX10633 CALL DOINX( 633 )                                               
CY10633 CALL DOINY( 633 )                                               
      DO 1530 I=1,M                                                     
 1530 WRITE(66,1540)I,(INDEX(I,J),J=1,M)                                
CX20634 CALL DOOUTX( 634 )                                              
CY20634 CALL DOOUTY( 634 )                                              
 1540 FORMAT(20X,I2,1X,10I5)                                            
      WRITE(66,1550)                                                    
 1550 FORMAT(//20X,' CHARGES ON COULOMB TERMS:'//)                      
      WRITE(66,1560) (ZITES(I),I=1,M)                                   
 1560 FORMAT(20X,10(5X,A4)//)                                           
      QSUM=0.D0                                                         
CX10635 CALL DOINX( 635 )                                               
CY10635 CALL DOINY( 635 )                                               
      DO 1570 I=1,M                                                     
      WRITE(66,1580) ZITES(I),(QQ(I,J),J=1,M)                           
      DO 1570 J=1,M                                                     
      QSUM=QSUM+QQ(I,J)                                                 
 1570 CONTINUE                                                          
CX20636 CALL DOOUTX( 636 )                                              
CY20636 CALL DOOUTY( 636 )                                              
 1580 FORMAT(20X,A4,10(1X,F8.5))                                        
      WRITE(66,1585) QSUM                                               
 1585 FORMAT(/20X,'..SUM OVER POINT CHARGES =',F5.2,' (SHOULD BE ZERO)')
CX10637 CALL DOINX( 637 )                                               
CY10637 CALL DOINY( 637 )                                               
      DO 150 I=1,4                                                      
      DO 150 J=1,4                                                      
      QQ(I,J)=QQ(I,J)*COULF                                             
  150 CONTINUE                                                          
CX20638 CALL DOOUTX( 638 )                                              
CY20638 CALL DOOUTY( 638 )                                              
CX40036 CALL SBOUTX (36)                                                
CY40036 CALL SBOUTY (36)                                                
      RETURN                                                            
      END                                                               
C================= CNVERTK.N.S. 5/30 85 === REV.000  X.X. X/XX XX ======
C                                                                       
      SUBROUTINE CNVERT(NEU,NLU,NT1,NT2)                                
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      CHARACTER*4  EU(3),LU(2),NEU,NLU,NT1,NT2                          
      SAVE EU, LU                                                       
      LOGICAL KOLLL,KOLLE                                               
      DIMENSION AC(5),BC(5),AD(5),BD(5)                                 
      LOGICAL MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW               
     X      ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC          
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
      COMMON/TFACTS/ FACTE(3),FACTL(2)                                  
      COMMON/TPARAM/ NCPBF(3),NPBF,NQ                                   
      COMMON/TABELS/ QQ(11,11),A(5),B(5)                                
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/CONSTS/ AVSNO,BOLTZ,PLANCK,PI,ELECHG,EPS0                  
      COMMON/LRCCEV/ ELRCNC,VLRCNC,ELRCEW,VLRCEW,PELRC(11),VILRC(11)    
      DATA EU/'KCAL','HART','MHAR'/                                     
      DATA LU/'ANGS','A.U.'/                                            
CX30037 CALL SBINX (37)                                                 
CY30037 CALL SBINY (37)                                                 
      WRITE(66,15)                                                      
   15 FORMAT(////20X,' P O T E N T I A L  I N F O R M A T I O N '//)    
      KK=0                                                              
      WRITE(66,20) NT1,NT2                                              
   20 FORMAT(//20X,'*** INTERACTION TYPE: ',1X,A4,1X,'-',1X,A4,' ***'//)
      NBF=0                                                             
CX10639 CALL DOINX( 639 )                                               
CY10639 CALL DOINY( 639 )                                               
      DO 30 II=1,NPBF                                                   
      KK=KK+1                                                           
      NC=NCPBF(KK)                                                      
      NBF=NBF+NC                                                        
   30 CONTINUE                                                          
CX20640 CALL DOOUTX( 640 )                                              
CY20640 CALL DOOUTY( 640 )                                              
      KOLLE=.FALSE.                                                     
      KOLLL=.FALSE.                                                     
CX10641 CALL DOINX( 641 )                                               
CY10641 CALL DOINY( 641 )                                               
      DO 40 I=1,3                                                       
      IF(NEU.NE.EU(I)) GO TO 40                                         
      KOLLE=.TRUE.                                                      
      EFACT=FACTE(I)                                                    
   40 CONTINUE                                                          
CX20642 CALL DOOUTX( 642 )                                              
CY20642 CALL DOOUTY( 642 )                                              
CX10643 CALL DOINX( 643 )                                               
CY10643 CALL DOINY( 643 )                                               
      DO 50 I=1,2                                                       
      IF(NLU.NE.LU(I)) GO TO 50                                         
      KOLLL=.TRUE.                                                      
      RFACT=FACTL(I)                                                    
   50 CONTINUE                                                          
CX20644 CALL DOOUTX( 644 )                                              
CY20644 CALL DOOUTY( 644 )                                              
      IF(.NOT.KOLLE) GO TO 900                                          
      IF(.NOT.KOLLL) GO TO 910                                          
      WRITE(66,60)                                                      
   60 FORMAT(20X,' POTENTIAL PARAMETERS:'//)                            
      WRITE(66,65) NEU,NLU                                              
   65 FORMAT(20X,'      A (',A4,')      B (',A4,')')                    
CX10645 CALL DOINX( 645 )                                               
CY10645 CALL DOINY( 645 )                                               
      DO 70 I=1,NBF                                                     
      WRITE(66,80) A(I),B(I)                                            
      AC(I)=A(I)*EFACT                                                  
      BC(I)=B(I)*RFACT                                                  
      AD(I)=AC(I)/UNITE                                                 
      BD(I)=BC(I)*UNITL                                                 
      A(I)=AD(I)                                                        
      B(I)=BD(I)                                                        
   70 CONTINUE                                                          
CX20646 CALL DOOUTX( 646 )                                              
CY20646 CALL DOOUTY( 646 )                                              
      A(NBF+1)=0.D0                                                     
      B(NBF+1)=0.D0                                                     
   80 FORMAT(20X,2(1X,F12.4))                                           
CX40038 CALL SBOUTX (38)                                                
CY40038 CALL SBOUTY (38)                                                
      RETURN                                                            
  900 WRITE(66,901)                                                     
CX40038 CALL SBOUTX (38)                                                
CY40038 CALL SBOUTY (38)                                                
      STOP                                                              
  910 WRITE(66,911)                                                     
CX40038 CALL SBOUTX (38)                                                
CY40038 CALL SBOUTY (38)                                                
      STOP                                                              
  901 FORMAT(1X,']]] ERROR IN ENERGY CONVERSION FACTORS ')              
  911 FORMAT(1X,']]] ERROR IN LENGTH CONVERSION FACTORS ')              
CX40038 CALL SBOUTX (38)                                                
CY40038 CALL SBOUTY (38)                                                
      END                                                               
C================= PCINITK.N.S. 3/25 84 === REV.000  X.X.  X/X XX ======
C                                                                       
      SUBROUTINE PCINIT(X2,Y2,Z2,QA0,QB0,QC0,QD0,QA1,QB1,QC1,QD1        
     X                 ,AX0,AY0,AZ0,AX1,AY1,AZ1,FX,FY,FZ,TX,TY,TZ)      
C     X                 ,SX,SY,SZ,FSX,FSY,FSZ)                          
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION X2(*),Y2(*),Z2(*),QA0(*),QB0(*),QC0(*),QD0(*)           
      DIMENSION QA1(*),QB1(*),QC1(*),QD1(*),AX0(*),AY0(*),AZ0(*)        
      DIMENSION AX1(*),AY1(*),AZ1(*),FX(*),FY(*),FZ(*)                  
      DIMENSION TX(*),TY(*),TZ(*)                                       
C      DIMENSION TX(*),TY(*),TZ(*),SX(*),SY(*),SZ(*)                    
C      DIMENSION FSX(*),FSY(*),FSZ(*)                                   
      COMMON/CONTRL/ NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP    
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
CX30039 CALL SBINX (39)                                                 
CY30039 CALL SBINY (39)                                                 
      NSP=NSPECI(1)                                                     
      I=0                                                               
CX10647 CALL DOINX( 647 )                                               
CY10647 CALL DOINY( 647 )                                               
      DO 200 N=1,NTYPES                                                 
      NSP=NSPECI(N)                                                     
      NS=NSITES(N)                                                      
      DO 200 J=1,NSP                                                    
      I=I+1                                                             
      X2(I)=FX(I)*TSTEP**2/2.D0/DSUMM(N)                                
      Y2(I)=FY(I)*TSTEP**2/2.D0/DSUMM(N)                                
      Z2(I)=FZ(I)*TSTEP**2/2.D0/DSUMM(N)                                
      IF(NS.EQ.1) GO TO 200                                             
      QA1(I)=0.5D0*(-QC0(I)*AX0(I)-QD0(I)*AY0(I)+QB0(I)*AZ0(I))*TSTEP   
      QB1(I)=0.5D0*(+QD0(I)*AX0(I)-QC0(I)*AY0(I)-QA0(I)*AZ0(I))*TSTEP   
      QC1(I)=0.5D0*(+QA0(I)*AX0(I)+QB0(I)*AY0(I)+QD0(I)*AZ0(I))*TSTEP   
      QD1(I)=0.5D0*(-QB0(I)*AX0(I)+QA0(I)*AY0(I)-QC0(I)*AZ0(I))*TSTEP   
  200 CONTINUE                                                          
CX20648 CALL DOOUTX( 648 )                                              
CY20648 CALL DOOUTY( 648 )                                              
      CALL ROTATE(QA0,QB0,QC0,QD0,TX,TY,TZ,NSPECI(1),.TRUE.)            
      I=0                                                               
CX10649 CALL DOINX( 649 )                                               
CY10649 CALL DOINY( 649 )                                               
      DO 301 N=1,NTYPES                                                 
      NSP=NSPECI(N)                                                     
      NS =NSITES(N)                                                     
      IF(NS.EQ.1) GO TO 301                                             
      DXX=DRI(1,N)                                                      
      DYY=DRI(2,N)                                                      
      DZZ=DRI(3,N)                                                      
      DO 300 J=1,NSP                                                    
      I=I+1                                                             
      IF(DXX.NE.0.D0) AX1(I)=TX(I)*TSTEP/DXX                            
      IF(DYY.NE.0.D0) AY1(I)=TY(I)*TSTEP/DYY                            
      IF(DZZ.NE.0.D0) AZ1(I)=TZ(I)*TSTEP/DZZ                            
  300 CONTINUE                                                          
  301 CONTINUE                                                          
CX20650 CALL DOOUTX( 650 )                                              
CY20650 CALL DOOUTY( 650 )                                              
CX40040 CALL SBOUTX (40)                                                
CY40040 CALL SBOUTY (40)                                                
      RETURN                                                            
      END                                                               
C================= PREDICK.N.S. 3/25 84 === REV.000  X.X.  X/X XX ======
C                                                                       
      SUBROUTINE PREDIC(T)                                              
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION T(*)                                                    
      LOGICAL MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW               
     X      ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC          
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/WINDEX/ NOR,IQ, IW, IA,IM,IF,IT,IP,IO,IX(99)               
CX30041 CALL SBINX (41)                                                 
CY30041 CALL SBINY (41)                                                 
      NSP3=3*NSPECI(1)                                                  
      NSP4=4*NSPECI(1)                                                  
CX10651 CALL DOINX( 651 )                                               
CY10651 CALL DOINY( 651 )                                               
      DO 100 I=0,NOP-1                                                  
      T(IX(IM+0)+I)=T(IX(IM)+0)-T(IX(1)+I)                              
      T(IX(IM+1)+I)=T(IX(IM)+1)-T(IX(2)+I)                              
 100  T(IX(IM+2)+I)=T(IX(IM)+2)-T(IX(3)+I)                              
CX20652 CALL DOOUTX( 652 )                                              
CY20652 CALL DOOUTY( 652 )                                              
C                                                                       
C     CALL PRED56(T(IX(1)) ,NOP ,0,3*NOP,3*NATOMS)                      
C     CALL PRED56(T(IX(2)) ,NOP ,0,3*NOP,2*NATOMS)                      
C     CALL PRED56(T(IX(3)) ,NOP ,0,3*NOP,  NATOMS)                      
C     CALL PRED56(T(IX(IQ)),NSP4,1,NSP4,0)                              
      CALL PRED5(T(IX(1)),T(IX(4)),T(IX(7)),T(IX(10)),T(IX(13)),T(IX(16)
     1),NOP)                                                            
      CALL PRED5(T(IX(2)),T(IX(5)),T(IX(8)),T(IX(11)),T(IX(14)),T(IX(17)
     1),NOP)                                                            
      CALL PRED5(T(IX(3)),T(IX(6)),T(IX(9)),T(IX(12)),T(IX(15)),T(IX(18)
     1),NOP)                                                            
      CALL PRED6(T(IX(19)),T(IX(23)),T(IX(27)),T(IX(31)),T(IX(35)),     
     1           T(IX(39)),T(IX(43)),NSP4)                              
                                                                        
      IF(LQ2PC) THEN                                                    
CX40042 CALL SBOUTX (42)                                                
CY40042 CALL SBOUTY (42)                                                
      RETURN                                                            
      ENDIF                                                             
      CALL PRED6(T(IX(47)),T(IX(50)),T(IX(53)),T(IX(56)),T(IX(59)),     
     1           T(IX(62)),T(IX(65)),NSP3)                              
CX40042 CALL SBOUTX (42)                                                
CY40042 CALL SBOUTY (42)                                                
C     CALL PRED56(T(IX(IA)),NSP3,1,NSP3,0)                              
      RETURN                                                            
      END                                                               
C================= PRED56K.N.S. 3/25 84 === REV.000  X.X.  X/X XX ======
C                                                                       
      SUBROUTINE PRED56(T,N,NN,NNN,NNNN)                                
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION T(*)                                                    
CX30043 CALL SBINX (43)                                                 
CY30043 CALL SBINY (43)                                                 
CX10653 CALL DOINX( 653 )                                               
CY10653 CALL DOINY( 653 )                                               
      DO 100 N0=1,N                                                     
      N1=N0+NNN+NNNN                                                    
      N2=N1+NNN                                                         
      N3=N2+NNN                                                         
      N4=N3+NNN                                                         
      N5=N4+NNN                                                         
      T0=T(N0)                                                          
      T1=T(N1)                                                          
      T2=T(N2)                                                          
      T3=T(N3)                                                          
      T4=T(N4)                                                          
      T5=T(N5)                                                          
      T(N0)=T0+T1+     T2+     T3+     T4+      T5                      
      T(N1)=   T1+2.D0*T2+3.D0*T3+4.D0*T4+ 5.D0*T5                      
      T(N2)=           T2+3.D0*T3+6.D0*T4+10.D0*T5                      
      T(N3)=                   T3+4.D0*T4+10.D0*T5                      
      T(N4)=                           T4+ 5.D0*T5                      
  100 CONTINUE                                                          
CX20654 CALL DOOUTX( 654 )                                              
CY20654 CALL DOOUTY( 654 )                                              
      IF(NN.EQ.0) THEN                                                  
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
      RETURN                                                            
      ENDIF                                                             
CX10655 CALL DOINX( 655 )                                               
CY10655 CALL DOINY( 655 )                                               
      DO 200 N0=1,N                                                     
      N1=N0+NNN+NNNN                                                    
      N2=N1+NNN                                                         
      N3=N2+NNN                                                         
      N4=N3+NNN                                                         
      N5=N4+NNN                                                         
      N6=N5+NNN                                                         
      T6=T(N6)                                                          
      T(N0)=T(N0)+      T6                                              
      T(N1)=T(N1)+ 6.D0*T6                                              
      T(N2)=T(N2)+15.D0*T6                                              
      T(N3)=T(N3)+20.D0*T6                                              
      T(N4)=T(N4)+15.D0*T6                                              
      T(N5)=T(N5)+ 6.D0*T6                                              
  200 CONTINUE                                                          
CX20656 CALL DOOUTX( 656 )                                              
CY20656 CALL DOOUTY( 656 )                                              
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE PRED5(X0,X1,X2,X3,X4,X5,N)                             
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION X0(*),X1(*),X2(*),X3(*),X4(*),X5(*)                     
CX30045 CALL SBINX (45)                                                 
CY30045 CALL SBINY (45)                                                 
CX10657 CALL DOINX( 657 )                                               
CY10657 CALL DOINY( 657 )                                               
      DO 100 K=1,N                                                      
      X0(K)=X0(K)+X1(K)+X2(K)+X3(K)+X4(K)+X5(K)                         
      X1(K)=    X1(K)+2.D0*X2(K)+3.D0*X3(K)+4.D0*X4(K)+ 5.D0*X5(K)      
      X2(K)=           X2(K)+3.D0*X3(K)+6.D0*X4(K)+10.D0*X5(K)          
      X3(K)=                   X3(K)+4.D0*X4(K)+10.D0*X5(K)             
      X4(K)=                           X4(K)+ 5.D0*X5(K)                
  100 CONTINUE                                                          
CX20658 CALL DOOUTX( 658 )                                              
CY20658 CALL DOOUTY( 658 )                                              
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE PRED6(X0,X1,X2,X3,X4,X5,X6,N)                          
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION X0(*),X1(*),X2(*),X3(*),X4(*),X5(*),X6(*)               
CX30047 CALL SBINX (47)                                                 
CY30047 CALL SBINY (47)                                                 
CX10659 CALL DOINX( 659 )                                               
CY10659 CALL DOINY( 659 )                                               
      DO 100 K=1,N                                                      
      X0(K)=X0(K)+X1(K)+X2(K)+X3(K)+X4(K)+X5(K)+X6(K)                   
      X1(K)=X1(K)+2.D0*X2(K)+3.D0*X3(K)+4.D0*X4(K)+5.D0*X5(K)+6.D0*X6(K)
      X2(K)=        X2(K)+3.D0*X3(K)+6.D0*X4(K)+10.D0*X5(K)+15.D0*X6(K) 
      X3(K)=                 X3(K)+4.D0*X4(K)+10.D0*X5(K)+20.D0*X6(K)   
      X4(K)=                         X4(K)+ 5.D0*X5(K)+15.D0*X6(K)      
      X5(K)=                                     X5(K)+6.D0*X6(K)       
  100 CONTINUE                                                          
CX20660 CALL DOOUTX( 660 )                                              
CY20660 CALL DOOUTY( 660 )                                              
CX40048 CALL SBOUTX (48)                                                
CY40048 CALL SBOUTY (48)                                                
      RETURN                                                            
      END                                                               
C================= CORRECK.N.S. 3/25 84 === REV.001K.N.S. 4/26 84 ======
C                                                                       
      SUBROUTINE CORREC(T,C)                                            
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION T(*),C(*)                                               
      LOGICAL MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW               
     X      ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC          
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/WINDEX/ NOR,IQ, IW, IA,IM,IF,IT,IP,IO,IX(99)               
      COMMON/CORRCC/F02,F12,F22,F32,F42,F52,F62                         
     X             ,F01,F11,F21,F31,F41,F51,F61                         
     X             ,F03,F13,F23,F33,F43,F53,F63                         
CX30049 CALL SBINX (49)                                                 
CY30049 CALL SBINY (49)                                                 
      ALENGT=22.D-10/UNITL                                              
      DT2=TSTEP**2/2.D0                                                 
      NSP=NSPECI(1)                                                     
      NSP3=3*NSP                                                        
      NSP4=4*NSP                                                        
      IAC=IX(7)                                                         
      IFO=IX(IF)                                                        
      I=-1                                                              
CX10661 CALL DOINX( 661 )                                               
CY10661 CALL DOINY( 661 )                                               
      DO 110 IJK=1,3                                                    
      DO 105 N=1,NTYPES                                                 
      DO 100 J=1,NSPECI(N)                                              
      I=I+1                                                             
      NA=IAC+I                                                          
      NF=IFO+I                                                          
      C(I+1)=T(NA)-T(NF)*DT2/DSUMM(N)                                   
  100 CONTINUE                                                          
  105 CONTINUE                                                          
  110 CONTINUE                                                          
CX20662 CALL DOOUTX( 662 )                                              
CY20662 CALL DOOUTY( 662 )                                              
      I1=1                                                              
      I2=NOP+I1                                                         
      I3=NOP+I2                                                         
      CALL CORR56(T(IX(1)),C(I1),F02,F12,F22,F32,F42,F52,F62,NOP,0      
     X,3*NOP,3*NATOMS)                                                  
      CALL CORR56(T(IX(2)),C(I2),F02,F12,F22,F32,F42,F52,F62,NOP,0      
     X,3*NOP,2*NATOMS)                                                  
      CALL CORR56(T(IX(3)),C(I3),F02,F12,F22,F32,F42,F52,F62,NOP,0      
     X,3*NOP,  NATOMS)                                                  
C*****************************************************************      
C*****  COMPUTE THE MEAN SQUARE DISPLACEMENT (M.S.D)         *****      
C*****  THIS WILL NOT GIVE CORRECT M.S.D. ONE HAS TO         *****      
C*****  REANALYZE THE SIMULATION DATA TO GET M.S.D.          *****      
C*****************************************************************      
      CALL ZERO(NOP3,C(1))                                              
CX10663 CALL DOINX( 663 )                                               
CY10663 CALL DOINY( 663 )                                               
      DO 200 I=0,NOP-1                                                  
  200 C(I+1)=T(IX(IM+0)+I)-T(IX(1)+I)                                   
CX20664 CALL DOOUTX( 664 )                                              
CY20664 CALL DOOUTY( 664 )                                              
      I1=NOP+1                                                          
CX10665 CALL DOINX( 665 )                                               
CY10665 CALL DOINY( 665 )                                               
      DO 210 I=0,NOP-1                                                  
  210 C(I1+I)=T(IX(IM+1)+I)-T(IX(2)+I)                                  
CX20666 CALL DOOUTX( 666 )                                              
CY20666 CALL DOOUTY( 666 )                                              
      I2=I1+NOP                                                         
CX10667 CALL DOINX( 667 )                                               
CY10667 CALL DOINY( 667 )                                               
      DO 220 I=0,NOP-1                                                  
  220 C(I2+I)=T(IX(IM+2)+I)-T(IX(3)+I)                                  
CX20668 CALL DOOUTX( 668 )                                              
CY20668 CALL DOOUTY( 668 )                                              
      I=0                                                               
CX10669 CALL DOINX( 669 )                                               
CY10669 CALL DOINY( 669 )                                               
      DO 700 N=1,NTYPES                                                 
      R2(N)=0.D0                                                        
      NSP=NSPECI(N)                                                     
      RFACT=DBLE(NSP)/UNITL**2                                          
      DO 750 J=1,NSP                                                    
      I=I+1                                                             
      J1=I+NOP                                                          
      J2=J1+NOP                                                         
      RM2=C(I)*C(I)+C(J1)*C(J1)+C(J2)*C(J2)                             
      R2(N)=R2(N)+RM2/RFACT                                             
  750 CONTINUE                                                          
  700 CONTINUE                                                          
CX20670 CALL DOOUTX( 670 )                                              
CY20670 CALL DOOUTY( 670 )                                              
      CALL ZERO(NOP3,C(1))                                              
C*****************************************************************      
CX10671 CALL DOINX( 671 )                                               
CY10671 CALL DOINY( 671 )                                               
      DO 711 I=0,NOP-1                                                  
      T(IX(1)+I)=T(IX(1)+I)                                             
     *          -2.D0*ALENGT*DBLE(INT(T(IX(1)+I)/ALENGT))               
      T(IX(2)+I)=T(IX(2)+I)                                             
     *          -2.D0*ALENGT*DBLE(INT(T(IX(2)+I)/ALENGT))               
      T(IX(3)+I)=T(IX(3)+I)-2.D0*DBLE(INT(T(IX(3)+I)))                  
  711 CONTINUE                                                          
C*****************************************************************      
C***** IF NOT SECOND ORDER QUATERNION METHOD THEN            *****      
C*****************************************************************      
CX20672 CALL DOOUTX( 672 )                                              
CY20672 CALL DOOUTY( 672 )                                              
      IF(.NOT.LQ2PC) THEN                                               
C************************                                               
      CALL ROTATE(T(IX(IQ)),T(IX(IQ+1)),T(IX(IQ+2)),T(IX(IQ+3)),        
     X            T(IX(IT)),T(IX(IT+1)),T(IX(IT+2)),NSPECI(1),.TRUE.)   
      NSP=NSPECI(1)                                                     
      IQN=IX(IQ)-1                                                      
      IAN=IX(IA)-1                                                      
CX10673 CALL DOINX( 673 )                                               
CY10673 CALL DOINY( 673 )                                               
      DO 300 I=1,NSP                                                    
      IQ1=IQN+I                                                         
      Q1=T(IQ1)                                                         
      II=I+NSP                                                          
      IQ2=IQN+II                                                        
      Q2=T(IQ2)                                                         
      III=II+NSP                                                        
      IQ3=IQN+III                                                       
      Q3=T(IQ3)                                                         
      IIII=III+NSP                                                      
      IQ4=IQN+IIII                                                      
      Q4=T(IQ4)                                                         
      IA1=IAN+I                                                         
      A1=T(IA1)                                                         
      IA2=IAN+II                                                        
      A2=T(IA2)                                                         
      IA3=IAN+III                                                       
      A3=T(IA3)                                                         
      CQ1=0.5D0*(-Q3*A1-Q4*A2+Q2*A3)                                    
      CQ2=0.5D0*(+Q4*A1-Q3*A2-Q1*A3)                                    
      CQ3=0.5D0*(+Q1*A1+Q2*A2+Q4*A3)                                    
      CQ4=0.5D0*(-Q2*A1+Q1*A2-Q3*A3)                                    
      JQ1=IQ1+NSP4                                                      
      JQ2=JQ1+NSP4                                                      
      JQ3=JQ2+NSP4                                                      
      JQ4=JQ3+NSP4                                                      
      C(I)=   T(JQ1)-CQ1*TSTEP                                          
      C(II)=  T(JQ2)-CQ2*TSTEP                                          
      C(III)= T(JQ3)-CQ3*TSTEP                                          
      C(IIII)=T(JQ4)-CQ4*TSTEP                                          
  300 CONTINUE                                                          
CX20674 CALL DOOUTX( 674 )                                              
CY20674 CALL DOOUTY( 674 )                                              
      CALL CORR56(T(IX(IQ)),C,F01,F11,F21,F31,F41,F51,F61,NSP4,1,NSP4,0)
      ITO=IX(IT)-1                                                      
      IAN=IX(IA+3)-1                                                    
      I=0                                                               
CX10675 CALL DOINX( 675 )                                               
CY10675 CALL DOINY( 675 )                                               
      DO 401 N=1,NTYPES                                                 
      NSP=NSPECI(N)                                                     
      NS=NSITES(N)                                                      
      IF(NS.EQ.1) GO TO 401                                             
      DXX=DRI(1,N)                                                      
      DYY=DRI(2,N)                                                      
      DZZ=DRI(3,N)                                                      
      DYMZOX=(DYY-DZZ)/DXX                                              
      DZMXOY=(DZZ-DXX)/DYY                                              
      DXMYOZ=(DXX-DYY)/DZZ                                              
      DO 400 J=0,NSP-1                                                  
      I=I+1                                                             
      NT=ITO+I                                                          
      NA=IAN+I                                                          
      C(I)=0.D0                                                         
      IF(DXX.NE.0.D0) C(I)=T(NA)-T(NT)*TSTEP/DXX                        
     X-DYMZOX*T(IX(IA+1)+J)*T(IX(IA+2)+J)*TSTEP                         
      II=I+NSP                                                          
      NT=ITO+II                                                         
      NA=IAN+II                                                         
      C(II)=0.D0                                                        
      IF(DYY.NE.0.D0) C(II)=T(NA)-T(NT)*TSTEP/DYY                       
     X-DZMXOY*T(IX(IA+2)+J)*T(IX(IA+0)+J)*TSTEP                         
      III=II+NSP                                                        
      NT=ITO+III                                                        
      NA=IAN+III                                                        
      C(III)=0.D0                                                       
      IF(DZZ.NE.0.D0) C(III)=T(NA)-T(NT)*TSTEP/DZZ                      
     X-DXMYOZ*T(IX(IA+0)+J)*T(IX(IA+1)+J)*TSTEP                         
  400 CONTINUE                                                          
  401 CONTINUE                                                          
CX20676 CALL DOOUTX( 676 )                                              
CY20676 CALL DOOUTY( 676 )                                              
      CALL CORR56(T(IX(IA)),C,F01,F11,F21,F31,F41,F51,F61,NSP3,1,NSP3,0)
C******************************************************************     
C***** IF SECOND ORDER QUATERNION METHOD THEN                  ****     
C******************************************************************     
      ELSE                                                              
C*******************                                                    
      CALL ROTATE(T(IX(IQ)),T(IX(IQ+1)),T(IX(IQ+2)),T(IX(IQ+3)),        
     X            T(IX(IT)),T(IX(IT+1)),T(IX(IT+2)),NSPECI(1),.TRUE.)   
      DXX=DRI(1,1)                                                      
      DYY=DRI(2,1)                                                      
      DZZ=DRI(3,1)                                                      
      DYMZOX=(DYY-DZZ)/DXX                                              
      DZMXOY=(DZZ-DXX)/DYY                                              
      DXMYOZ=(DXX-DYY)/DZZ                                              
CX10677 CALL DOINX( 677 )                                               
CY10677 CALL DOINY( 677 )                                               
      DO 900 J=0,NSPECI(1)-1                                            
      Q1=T(IX(IQ)+J)                                                    
      Q2=T(IX(IQ+1)+J)                                                  
      Q3=T(IX(IQ+2)+J)                                                  
      Q4=T(IX(IQ+3)+J)                                                  
      QA=T(IX(IQ+4)+J)                                                  
      QB=T(IX(IQ+5)+J)                                                  
      QC=T(IX(IQ+6)+J)                                                  
      QD=T(IX(IQ+7)+J)                                                  
      A14=-Q3*QA+Q4*QB+Q1*QC-Q2*QD                                      
      A24=-Q4*QA-Q3*QB+Q2*QC+Q1*QD                                      
      A34= Q2*QA-Q1*QB+Q4*QC-Q3*QD                                      
      T(IX(IA)+J)=  2.D0*A14/TSTEP                                      
      T(IX(IA+1)+J)=2.D0*A24/TSTEP                                      
      T(IX(IA+2)+J)=2.D0*A34/TSTEP                                      
      A1=DYMZOX*A24*A34+0.5D0*DT2*T(IX(IT)+J)/DXX                       
      A2=DZMXOY*A14*A34+0.5D0*DT2*T(IX(IT+1)+J)/DYY                     
      A3=DXMYOZ*A14*A24+0.5D0*DT2*T(IX(IT+2)+J)/DZZ                     
      A4=-0.5D0*(QA*QA+QB*QB+QC*QC+QD*QD)                               
      C(J+1)=T(IX(IQ+8)+J)-(-A1*Q3-A2*Q4+A3*Q2+A4*Q1)                   
      K=J+NSPECI(1)                                                     
      C(K+1)=T(IX(IQ+9)+J)-(A1*Q4-A2*Q3-A3*Q1+A4*Q2)                    
      L=K+NSPECI(1)                                                     
      C(L+1)=T(IX(IQ+10)+J)-(A1*Q1+A2*Q2+A3*Q4+A4*Q3)                   
      M=L+NSPECI(1)                                                     
      C(M+1)=T(IX(IQ+11)+J)-(-A1*Q2+A2*Q1-A3*Q3+A4*Q4)                  
  900 CONTINUE                                                          
CX20678 CALL DOOUTX( 678 )                                              
CY20678 CALL DOOUTY( 678 )                                              
      CALL CORR56(T(IX(IQ)),C(1),F03,F13,F23,F33,F43,F53,F63,           
     X            NSP4,1,NSP4,0)                                        
      ENDIF                                                             
C************ THIS COMPLETES SECOND ORDER QUATERNION METHOD  *****      
C*****************************************************************      
      NSP=NSPECI(1)                                                     
      IQN=IX(IQ)-1                                                      
CX10679 CALL DOINX( 679 )                                               
CY10679 CALL DOINY( 679 )                                               
      DO 350 I=1,NSP                                                    
      IQ1=IQN+I                                                         
      Q1=T(IQ1)                                                         
      IQ2=IQ1+NSP                                                       
      Q2=T(IQ2)                                                         
      IQ3=IQ2+NSP                                                       
      Q3=T(IQ3)                                                         
      IQ4=IQ3+NSP                                                       
      Q4=T(IQ4)                                                         
      QNORM=SQRT(Q1*Q1+Q2*Q2+Q3*Q3+Q4*Q4)                               
      IF(QNORM.NE.0.D0) QNORM=1.D0/QNORM                                
      T(IQ1)=QNORM*Q1                                                   
      T(IQ2)=QNORM*Q2                                                   
      T(IQ3)=QNORM*Q3                                                   
      T(IQ4)=QNORM*Q4                                                   
  350 CONTINUE                                                          
C*****************************************************************      
CX20680 CALL DOOUTX( 680 )                                              
CY20680 CALL DOOUTY( 680 )                                              
CX10681 CALL DOINX( 681 )                                               
CY10681 CALL DOINY( 681 )                                               
      DO 500 I=0,NOP-1                                                  
  500 C(I+1)=T(IX(IW)+I)**2+T(IX(IW+1)+I)**2+T(IX(IW+2)+I)**2           
CX20682 CALL DOOUTX( 682 )                                              
CY20682 CALL DOOUTY( 682 )                                              
      NS=1                                                              
CX10683 CALL DOINX( 683 )                                               
CY10683 CALL DOINY( 683 )                                               
      DO 550 N=1,NTYPES                                                 
      NSP=NSPECI(N)                                                     
      TE(N)=(SSUM(NSP,C(NS),1)/DT2)*0.25D0*DSUMM(N)                     
      NS=NS+NSP                                                         
  550 CONTINUE                                                          
CX20684 CALL DOOUTX( 684 )                                              
CY20684 CALL DOOUTY( 684 )                                              
      NSP=NSPECI(1)                                                     
CX10685 CALL DOINX( 685 )                                               
CY10685 CALL DOINY( 685 )                                               
      DO 600 I=0,NSP-1                                                  
      C(I+1)=T(IX(IA)+I)**2                                             
      J=I+NSP                                                           
      C(J+1)=T(IX(IA+1)+I)**2                                           
      K=J+NSP                                                           
      C(K+1)=T(IX(IA+2)+I)**2                                           
  600 CONTINUE                                                          
CX20686 CALL DOOUTX( 686 )                                              
CY20686 CALL DOOUTY( 686 )                                              
      IS=1                                                              
      JS=IS+NSP                                                         
      KS=JS+NSP                                                         
CX10687 CALL DOINX( 687 )                                               
CY10687 CALL DOINY( 687 )                                               
      DO 650 N=1,NTYPES                                                 
      IF(NSITES(N).EQ.1) GO TO 650                                      
      NSP=NSPECI(N)                                                     
      RE(N)=0.5D0*(DRI(1,N)*SSUM(NSP,C(IS),1)                           
     X            +DRI(2,N)*SSUM(NSP,C(JS),1)                           
     X            +DRI(3,N)*SSUM(NSP,C(KS),1))                          
      IS=IS+NSP                                                         
      JS=JS+NSP                                                         
      KS=KS+NSP                                                         
  650 CONTINUE                                                          
CX20688 CALL DOOUTX( 688 )                                              
CY20688 CALL DOOUTY( 688 )                                              
CX40050 CALL SBOUTX (50)                                                
CY40050 CALL SBOUTY (50)                                                
      RETURN                                                            
      END                                                               
C================= CORR56K.N.S. 3/25 84 === REV.000  X.X.  X/X XX ======
C                                                                       
      SUBROUTINE CORR56(T,C,C0,C1,C2,C3,C4,C5,C6,N,NN,NNN,NNNN)         
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION T(*),C(*)                                               
CX30051 CALL SBINX (51)                                                 
CY30051 CALL SBINY (51)                                                 
CX10689 CALL DOINX( 689 )                                               
CY10689 CALL DOINY( 689 )                                               
      DO 100 N0=1,N                                                     
      N1=N0+NNN+NNNN                                                    
      N2=N1+NNN                                                         
      N3=N2+NNN                                                         
      N4=N3+NNN                                                         
      N5=N4+NNN                                                         
      T(N0)=T(N0)-C(N0)*C0                                              
      T(N1)=T(N1)-C(N0)*C1                                              
      T(N2)=T(N2)-C(N0)*C2                                              
      T(N3)=T(N3)-C(N0)*C3                                              
      T(N4)=T(N4)-C(N0)*C4                                              
      T(N5)=T(N5)-C(N0)*C5                                              
  100 CONTINUE                                                          
CX20690 CALL DOOUTX( 690 )                                              
CY20690 CALL DOOUTY( 690 )                                              
      IF(NN.EQ.0) THEN                                                  
CX40052 CALL SBOUTX (52)                                                
CY40052 CALL SBOUTY (52)                                                
      RETURN                                                            
      ENDIF                                                             
CX10691 CALL DOINX( 691 )                                               
CY10691 CALL DOINY( 691 )                                               
      DO 200 N0=1,N                                                     
      N1=N0+NNN+NNNN                                                    
      N2=N1+NNN                                                         
      N3=N2+NNN                                                         
      N4=N3+NNN                                                         
      N5=N4+NNN                                                         
      N6=N5+NNN                                                         
      T(N6)=T(N6)-C(N0)*C6                                              
  200 CONTINUE                                                          
CX20692 CALL DOOUTX( 692 )                                              
CY20692 CALL DOOUTY( 692 )                                              
CX40052 CALL SBOUTX (52)                                                
CY40052 CALL SBOUTY (52)                                                
      RETURN                                                            
      END                                                               
C================= SSUM    =============================================
C                                                                       
      FUNCTION SSUM(N,A,NJ)                                             
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION A(*)                                                    
CX30053 CALL SBINX (53)                                                 
CY30053 CALL SBINY (53)                                                 
      SUM=0.D0                                                          
CX10693 CALL DOINX( 693 )                                               
CY10693 CALL DOINY( 693 )                                               
      DO 10 I=1,N                                                       
      J=I-1                                                             
      SUM=SUM+A(J*NJ+1)                                                 
   10 CONTINUE                                                          
CX20694 CALL DOOUTX( 694 )                                              
CY20694 CALL DOOUTY( 694 )                                              
      SSUM=SUM                                                          
CX40054 CALL SBOUTX (54)                                                
CY40054 CALL SBOUTY (54)                                                
      RETURN                                                            
      END                                                               
C================= ZERO    =============================================
C                                                                       
      SUBROUTINE ZERO(M,A)                                              
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION A(M)                                                    
CX30055 CALL SBINX (55)                                                 
CY30055 CALL SBINY (55)                                                 
CX10695 CALL DOINX( 695 )                                               
CY10695 CALL DOINY( 695 )                                               
      DO 10 I=1,M                                                       
      A(I)=0.D0                                                         
   10 CONTINUE                                                          
CX20696 CALL DOOUTX( 696 )                                              
CY20696 CALL DOOUTY( 696 )                                              
CX40056 CALL SBOUTX (56)                                                
CY40056 CALL SBOUTY (56)                                                
      RETURN                                                            
      END                                                               
C================= SCALETK.N.S. 4/2  84 === REV.000  X.X.  X/X XX ======
C                                                                       
      SUBROUTINE SCALET(VX,VY,VZ,AX,AY,AZ)                              
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION VX(*),VY(*),VZ(*),AX(*),AY(*),AZ(*)                     
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      COMMON/AVERIZ/ V(13),VV(13)                                       
CX30057 CALL SBINX (57)                                                 
CY30057 CALL SBINY (57)                                                 
      IF(NSTEP.GE.NSSTS) WRITE(66,10) NSTEP                             
   10 FORMAT(1X,'*** TEMPERATURE SCALING SWITHED OFF AT TIME STEP ',I6) 
      IF(NSTEP.GE.NSSTS) LTS=.FALSE.                                    
      IF(MOD(NSTEP,NSBTS).NE.0) THEN                                    
CX40058 CALL SBOUTX (58)                                                
CY40058 CALL SBOUTY (58)                                                
      RETURN                                                            
      ENDIF                                                             
      I=0                                                               
CX10697 CALL DOINX( 697 )                                               
CY10697 CALL DOINY( 697 )                                               
      DO 30 N=1,NTYPES                                                  
      NS=NSITES(N)                                                      
      NSP=NSPECI(N)                                                     
      TSCAL=1.D0                                                        
      IF(TTS(N).NE.0.D0) TSCAL=SQRT(TTRAN(N)/TTS(N))                    
      DO 30 J=1,NSP                                                     
      I=I+1                                                             
      VX(I)=VX(I)*TSCAL                                                 
      VY(I)=VY(I)*TSCAL                                                 
      VZ(I)=VZ(I)*TSCAL                                                 
   30 CONTINUE                                                          
CX20698 CALL DOOUTX( 698 )                                              
CY20698 CALL DOOUTY( 698 )                                              
      I=0                                                               
CX10699 CALL DOINX( 699 )                                               
CY10699 CALL DOINY( 699 )                                               
      DO 36 N=1,NTYPES                                                  
      NS=NSITES(N)                                                      
      NSP=NSPECI(N)                                                     
      RSCAL=1.D0                                                        
      IF(NS.EQ.1) GO TO 36                                              
      IF(RTS(N).NE.0.D0) RSCAL=SQRT(TROT(N)/ RTS(N))                    
      DO 35 J=1,NSP                                                     
      I=I+1                                                             
      AX(I)=AX(I)*RSCAL                                                 
      AY(I)=AY(I)*RSCAL                                                 
      AZ(I)=AZ(I)*RSCAL                                                 
   35 CONTINUE                                                          
   36 CONTINUE                                                          
CX20700 CALL DOOUTX( 700 )                                              
CY20700 CALL DOOUTY( 700 )                                              
      NAVRS=0                                                           
CX10701 CALL DOINX( 701 )                                               
CY10701 CALL DOINY( 701 )                                               
      DO 40 I=1,13                                                      
      V(I)=0.D0                                                         
      VV(I)=0.D0                                                        
   40 CONTINUE                                                          
CX20702 CALL DOOUTX( 702 )                                              
CY20702 CALL DOOUTY( 702 )                                              
CX40058 CALL SBOUTX (58)                                                
CY40058 CALL SBOUTY (58)                                                
      RETURN                                                            
      END                                                               
C================= RAND    =============================================
C                                                                       
       FUNCTION RAND(N)                                                 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      COMMON /RANDNO/ R3(127),R1,I2                                     
      SAVE S, T, RMC, IW                                                
      DATA S,T,RMC/0.D0, 1.D0, 1.D0/                                    
      DATA IW/-1/                                                       
CX30059 CALL SBINX (59)                                                 
CY30059 CALL SBINY (59)                                                 
      IF((R1.LT.1.D0).AND.(N.EQ.0)) GO TO 60                            
       IF (IW.GT.0) GO TO 30                                            
   10  IW = IW + 1                                                      
       T = 0.5D0*T                                                      
       R1 = S                                                           
       S = S + T                                                        
       IF ((S.GT.R1).AND.(S.LT.1D0)) GO TO 10                           
       IKT = (IW - 1)/12                                                
       IC = IW - 12*IKT                                                 
       ID = 2**(13 - IC)                                                
CX10703 CALL DOINX( 703 )                                               
CY10703 CALL DOINY( 703 )                                               
       DO 20 I = 1, IC                                                  
   20  RMC = 0.5D0*RMC                                                  
CX20704 CALL DOOUTX( 704 )                                              
CY20704 CALL DOOUTY( 704 )                                              
       RM = 0.015625D0*0.015625D0                                       
   30  I2 = 127                                                         
       IR = MOD(ABS(N), 8190) + 1                                       
   40  R1 = 0.D0                                                        
CX10705 CALL DOINX( 705 )                                               
CY10705 CALL DOINY( 705 )                                               
       DO 50 I = 1, IKT                                                 
       IR = MOD(17*IR, 8191)                                            
   50  R1 = (R1 + DBLE(IR/2))*RM                                        
CX20706 CALL DOOUTX( 706 )                                              
CY20706 CALL DOOUTY( 706 )                                              
       IR = MOD (17*IR, 8191)                                           
       R1 = (R1 + DBLE(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.1D0) T = (R1 - 0.5D0) + (R3(I2) - 0.5D0)                
       R1 = T                                                           
       R3(I2) = R1                                                      
       I2 = I2 - 1                                                      
       RAND = R1                                                        
CX40060 CALL SBOUTX (60)                                                
CY40060 CALL SBOUTY (60)                                                
       RETURN                                                           
       END                                                              
C================== JACOBI =============================================
C                                                                       
      SUBROUTINE JACOBI(A,V,N)                                          
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION A(3,3),V(3,3)                                           
CX30061 CALL SBINX (61)                                                 
CY30061 CALL SBINY (61)                                                 
      N=3                                                               
      RHO=1.0D-8                                                        
      TES=0.D0                                                          
      SCL=0.D0                                                          
CX10707 CALL DOINX( 707 )                                               
CY10707 CALL DOINY( 707 )                                               
      DO 10 I=1,N                                                       
      SCL=SCL+A(I,I)*A(I,I)                                             
   10 CONTINUE                                                          
CX20708 CALL DOOUTX( 708 )                                              
CY20708 CALL DOOUTY( 708 )                                              
      SCL=SQRT(SCL)/DBLE(N)                                             
CX10709 CALL DOINX( 709 )                                               
CY10709 CALL DOINY( 709 )                                               
      DO 20 I=1,N                                                       
      DO 20 J=1,N                                                       
   20 A(I,J)=A(I,J)/SCL                                                 
CX20710 CALL DOOUTX( 710 )                                              
CY20710 CALL DOOUTY( 710 )                                              
CX10711 CALL DOINX( 711 )                                               
CY10711 CALL DOINY( 711 )                                               
      DO 30 I=1,N                                                       
      DO 30 J=1,N                                                       
      V(I,J)=0.D0                                                       
      IF(I.EQ.J)V(I,J)=1.D0                                             
   30 CONTINUE                                                          
CX20712 CALL DOOUTX( 712 )                                              
CY20712 CALL DOOUTY( 712 )                                              
CX10713 CALL DOINX( 713 )                                               
CY10713 CALL DOINY( 713 )                                               
      DO 100 I=2,N                                                      
      IU1=I-1                                                           
      DO 100 J=1,IU1                                                    
  100 TES=TES+2.D0*A(I,J)*A(I,J)                                        
CX20714 CALL DOOUTX( 714 )                                              
CY20714 CALL DOOUTY( 714 )                                              
      TES=SQRT(TES)                                                     
      M=0                                                               
  105 TES=TES/DBLE(N)                                                   
      IF(TES.LT.RHO)TES=RHO                                             
  110 CONTINUE                                                          
CX10715 CALL DOINX( 715 )                                               
CY10715 CALL DOINY( 715 )                                               
      DO      165      I=2,N                                            
      IUU1=I-1                                                          
      DO 165 J=1,IUU1                                                   
      IF(ABS(A(I,J))-TES)165,115,115                                    
  115 M=1                                                               
      V1=A(J,J)                                                         
      V2=A(I,J)                                                         
      V3=A(I,I)                                                         
      U=0.5D0*(V1-V3)                                                   
      IF(ABS(U)-RHO)120,125,125                                         
  120 OMG=-1.D0                                                         
      GO TO 130                                                         
  125 SQVVUU=SQRT(V2*V2+U*U)                                            
      OMG=-V2/SQVVUU                                                    
      IF(U.LT.0.D0)OMG=-OMG                                             
  130 SQOOO=SQRT(2.D0*(1.D0+SQRT(1.D0-OMG*OMG)))                        
      S=OMG/SQOOO                                                       
      C=SQRT(1.D0-S*S)                                                  
      DO 160 K=1,N                                                      
      IF(K-I)140,135,135                                                
  135 TEM=A(K,J)*C-A(K,I)*S                                             
      A(K,I)=A(K,J)*S+A(K,I)*C                                          
      A(K,J)=TEM                                                        
      GO TO 155                                                         
  140 IF(K-J)145,150,150                                                
  145 TEM=A(J,K)*C-A(I,K)*S                                             
      A(I,K)=A(J,K)*S+A(I,K)*C                                          
      A(J,K)=TEM                                                        
      GO TO 155                                                         
  150 TEM=A(K,J)*C-A(I,K)*S                                             
      A(I,K)=A(K,J)*S+A(I,K)*C                                          
      A(K,J)=TEM                                                        
  155 TEM=V(K,J)*C-V(K,I)*S                                             
      V(K,I)=V(K,J)*S+V(K,I)*C                                          
      V(K,J)=TEM                                                        
  160 CONTINUE                                                          
      A(J,J)=V1*C*C+V3*S*S-2.D0*V2*S*C                                  
      A(I,I)=V1*S*S+V3*C*C+2.D0*V2*S*C                                  
      A(I,J)=(V1-V3)*S*C+V2*(C*C-S*S)                                   
  165 CONTINUE                                                          
CX20716 CALL DOOUTX( 716 )                                              
CY20716 CALL DOOUTY( 716 )                                              
      IF(M-1)175,170,170                                                
  170 M=0                                                               
      GO TO 110                                                         
  175 IF(TES-RHO)180,180,105                                            
  180 CONTINUE                                                          
CX10717 CALL DOINX( 717 )                                               
CY10717 CALL DOINY( 717 )                                               
      DO      190      I=1,N                                            
      DO 190 J=1,N                                                      
  190 A(I,J)=SCL*A(I,J)                                                 
CX20718 CALL DOOUTX( 718 )                                              
CY20718 CALL DOOUTY( 718 )                                              
CX40062 CALL SBOUTX (62)                                                
CY40062 CALL SBOUTY (62)                                                
      RETURN                                                            
      END                                                               
C================= VALUESK.N.S. 3/30 84 === REV.001K.N.S. 4/27 84 ======
C                                                                       
      SUBROUTINE VALUES                                                 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
      COMMON/LRCCEV/ ELRCNC,VLRCNC,ELRCEW,VLRCEW,PELRC(11),VILRC(11)    
      COMMON/CONSTS/ AVSNO,BOLTZ,PLANCK,PI,ELECHG,EPS0                  
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      COMMON/AVERIZ/ SUWW ,SUWI ,SUWA ,SUII ,SUIA ,SUTOT ,STKEW ,SRKEW  
     X              ,STKEI      ,SETOT      ,STTW        ,SRTW  ,STTI   
     X             ,SSUWW,SSUWI,SSUWA,SSUII,SSUIA,SSUTOT,SSTKEW,SSRKEW  
     X             ,SSTKEI     ,SSETOT     ,SSTTW       ,SSRTW ,SSTTI   
CX30063 CALL SBINX (63)                                                 
CY30063 CALL SBINY (63)                                                 
      FACE=AVSNO*UNITE                                                  
      NSP=NSPECI(1)                                                     
      NOP=NSP+NION                                                      
C      NTOT=NOP+NATOMS                                                  
      EFACW= FACE/DBLE(NSP)                                             
      IF(NION.GT.0)                                                     
     XEFACI= FACE/DBLE(NION)                                            
C      EFACWI=FACE/DBLE(NSP+NION)                                       
C      EFACWA=FACE/DBLE(NSP+NATOMS)                                     
C      EFACIA=FACE/DBLE(NION+NATOMS)                                    
C      EFACT= FACE/DBLE(NTOT)                                           
      TKEW=TE(1)*EFACW                                                  
      TKEI=TE(2)*EFACI                                                  
      TTW=TKEW*TEMPF                                                    
      TTI=TKEI*TEMPF                                                    
      RKEW=RE(1)*EFACW                                                  
      RTW=RKEW*TEMPF                                                    
      UWW=EWW*EFACW                                                     
      UWI=EWI*EFACW                                                     
      UWA=EWA*EFACW                                                     
      UII=EII*EFACI                                                     
      UIA=EIA*EFACI                                                     
      UTOT=UWW+UWI+UWA+UII+UIA                                          
      PTOT=(EWW+EWI+EWA+EII+EIA)*EFACW                                  
      ETOT=UTOT+TKEW+TKEI+RKEW                                          
      TTS(1)=TTW                                                        
      TTS(2)=TTI                                                        
      RTS(1)=RTW                                                        
      RTS(2)=0.D0                                                       
      IF(NSTEP.EQ.0) GO TO 555                                          
      IF(MOD(NSTEP,NSBPO).NE.0) GO TO 170                               
  555 CONTINUE                                                          
      WRITE(66,50)                                                      
   50 FORMAT(1X,128('-'))                                               
      WRITE(66,100) NSTEP,UWW,UWI,UWA,UII,UIA                           
  100 FORMAT(/1X,'STEP',I6,' ENERGIES (KJ/MOL) - WW: ',-3PF12.4,' WI: ',
     X-3PF12.4,' WA: ',-3PF12.4,' II: ',-3PF12.4,' IA: ',-3PF12.4)      
      WRITE(66,120) TKEW,RKEW,TKEI                                      
  120 FORMAT(/1X,'KINETIC ENERGIES (KJ/MOL)  H2O/TR: ',-3PF12.4,' H2O/RO
     XT: ',-3PF12.4,' ION/TR: ',-3PF12.4)                               
      WRITE(66,140) TTW,RTW,TTI                                         
  140 FORMAT(/1X,'TEMPERATURES     (K)       H2O/TR: ',0PF12.2,'  H2O/RO
     XT: ',0PF12.2,' ION/TR: ',0PF12.2)                                 
      WRITE(66,160) UTOT,PTOT,ETOT,R2(1),R2(2)                          
  160 FORMAT(/1X,'TOT.ENERG.(KJ/MOL) U :',-3PF12.4,' (',-3PF12.4,') U+KI
     XN: ',-3PF12.4,' >>> MSD (M**2) W: ',1PD10.4,' I: ',1PD10.4)       
  170 CONTINUE                                                          
      NAVRS=NAVRS+1                                                     
      SUWW=SUWW+UWW                                                     
      SUWI=SUWI+UWI                                                     
      SUWA=SUWA+UWA                                                     
      SUII=SUII+UII                                                     
      SUIA=SUIA+UIA                                                     
      SUTOT=SUTOT+UTOT                                                  
      STKEW=STKEW+TKEW                                                  
      SRKEW=SRKEW+RKEW                                                  
      STKEI=STKEI+TKEI                                                  
      SETOT=SETOT+ETOT                                                  
      STTW=STTW+TTW                                                     
      SRTW=SRTW+RTW                                                     
      STTI=STTI+TTI                                                     
      SSUWW=SSUWW+UWW**2                                                
      SSUWI=SSUWI+UWI**2                                                
      SSUWA=SSUWA+UWA**2                                                
      SSUII=SSUII+UII**2                                                
      SSUIA=SSUIA+UIA**2                                                
      SSUTOT=SSUTOT+UTOT**2                                             
      SSTKEW=SSTKEW+TKEW**2                                             
      SSRKEW=SSRKEW+RKEW**2                                             
      SSTKEI=SSTKEI+TKEI**2                                             
      SSETOT=SSETOT+ETOT**2                                             
      SSTTW=SSTTW+TTW**2                                                
      SSRTW=SSRTW+RTW**2                                                
      SSTTI=SSTTI+TTI**2                                                
CX40064 CALL SBOUTX (64)                                                
CY40064 CALL SBOUTY (64)                                                
      RETURN                                                            
      END                                                               
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C================= INPUT K.N.S. 3/17 84 === REV.001K.N.S. 4/26 84 ===== 
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C                                                                       
      SUBROUTINE INPUT                                                  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      INTEGER SITES                                                     
      DOUBLE PRECISION MASS                                             
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONTRL/ NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP    
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/CONSTS/ AVSNO,BOLTZ,PLANCK,PI,ELECHG,EPS0                  
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      DIMENSION      CHARGE(10),DI(3,3),MASS(30)                        
      DIMENSION      R(3,10),SITE(3,30),AMASS(30),PMOI(3,3)             
     X              ,SUMM(3),Q(30),RRI(3,3),RMOI(3,3)                   
CX30065 CALL SBINX (65)                                                 
CY30065 CALL SBINY (65)                                                 
C     NAMELIST /SETUP /RHO,NOP,NTYPES,NENT,NION,NATOMS                  
C    X                ,MD,LSAVE,LRES,INMC,NPOCS                         
C     NAMELIST /SPECIE/NUMBER,SITES,CHARGE,MASS,ROTEMP,TRTEMP           
C     NAMELIST /MOLDYN/NSTEPS,NSBPO,NSBTS,NSSTS,DT,BREAK                
C    X                ,LTS,LLF,LPC,NOTR                                 
      CALL TYPE (1)                                                     
      READ(5,*)MD                                                       
      READ(5,*)LRES                                                     
      READ(5,*)LSAVE                                                    
      READ(5,*)INMC                                                     
      READ(5,1881)NTYPES,NOP,NION,NATOMS,NENT                           
      IF(MD) READ(5,*)LTS                                               
      IF(MD) READ(5,*)LLF                                               
      IF(MD) READ(5,*)LPC                                               
      IF(MD) READ(5,*)LQ2PC                                             
      IF(MD) READ(5,1685)DT                                             
 1685 FORMAT(D10.5)                                                     
      IF(MD) READ(5,1881)NSTEPS,NSBPO,NSSTS                             
 1881 FORMAT(6I10)                                                      
C     READ  (5,SETUP)                                                   
      IF(NION.GT.50) GO TO 900                                          
      FNOP=DBLE(NOP)                                                    
      K=0                                                               
      NSP=0                                                             
CX10719 CALL DOINX( 719 )                                               
CY10719 CALL DOINY( 719 )                                               
      DO 75 I=1,NTYPES                                                  
C     READ  (5,SPECIE)                                                  
      READ(5,1881) NUMBER,SITES                                         
      READ(5,1889) TRTEMP,ROTEMP                                        
 1889 FORMAT(2F10.5)                                                    
      NSITES(I)=SITES                                                   
      NSPECI(I)=NUMBER                                                  
      TTRAN(I)=TRTEMP                                                   
      TROT(I) =ROTEMP                                                   
      CALL COORDS(SITES,R)                                              
      DO 50 J=1,SITES                                                   
      K=K+1                                                             
      DO 10 LL=1,3                                                      
   10 SITE(LL,K)=R(LL,J)                                                
      READ(5,1898) MASS(J),CHARGE(J)                                    
 1898 FORMAT(2D10.5)                                                    
      Q(K)=CHARGE(J)                                                    
      AMASS(K)=MASS(J)                                                  
   50 CONTINUE                                                          
      L=K+1-SITES                                                       
      IF(SITES.NE.1) CALL FIXCOM(K,L,AMASS,SITE,PMOI,SUM,RMOI)          
      IF(SITES.EQ.1) SUM=AMASS(K)                                       
      SUMM(I)=SUM                                                       
      IF(SITES.EQ.1) GO TO 70                                           
      DO 60 J=1,3                                                       
   60 RRI(J,I)=PMOI(J,J)                                                
   70 NSP=NSP+NUMBER                                                    
   75 CONTINUE                                                          
CX20720 CALL DOOUTX( 720 )                                              
CY20720 CALL DOOUTY( 720 )                                              
      NSITET=K                                                          
      IF(NSP.NE.NOP) GO TO 910                                          
C     IF(MD) READ(5,MOLDYN)                                             
      CALL UNITS(SITE,Q,AMASS,SUMM,RRI)                                 
      WRITE(66,100) NTYPES                                              
      WRITE(66,150)                                                     
      K=0                                                               
CX10721 CALL DOINX( 721 )                                               
CY10721 CALL DOINY( 721 )                                               
      DO 95 I=1,NTYPES                                                  
      NS=NSITES(I)                                                      
      WRITE(66,110) I                                                   
      WRITE(66,120) NSPECI(I),NSITES(I)                                 
      DO 80 J=1,3                                                       
   80 DI(J,I)=DM(J,I)/3.33564D-30                                       
      WRITE(66,220) SUMM(I)                                             
      WRITE(66,230) TTRAN(I)                                            
      IF(NS.EQ.1) GO TO 85                                              
      WRITE(66,240) TROT(I)                                             
      WRITE(66,250) (RRI(J,I),J=1,3)                                    
      WRITE(66,260) (DM(J,I),J=1,3)                                     
      WRITE(66,270) (DI(J,I),J=1,3)                                     
C...  WRITE(66,280) (QM(J,I),J=1,3)                                     
C...  WRITE(66,285) (QM(J,I),J=4,9)                                     
   85 CONTINUE                                                          
      WRITE(66,290)                                                     
      DO 95 J=1,NS                                                      
      K=K+1                                                             
      WRITE(66,300) AMASS(K),DQ(K),SITE(1,K),SITE(2,K),SITE(3,K)        
      WRITE(6,*) AMASS(K),DQ(K),SITE(1,K),SITE(2,K),SITE(3,K),'******'  
   95 CONTINUE                                                          
CX20722 CALL DOOUTX( 722 )                                              
CY20722 CALL DOOUTY( 722 )                                              
      WRITE(66,310) NOP,NATOMS                                          
  100 FORMAT(/1X,'***            INPUT DATA PROCESSED FOR ',I2,' TYPE(S)
     X OF SPECII'//)                                                    
  110 FORMAT(/1X,'*** T Y P E : ',I2//)                                 
  120 FORMAT(/1X,'*** ... CONSISTS OF ',I5,'   PARTICLES       -    AND 
     XTHERE ARE ',I2,' POINT MASSES / CHARGES (SITES) IN EACH PARTICLE')
  150 FORMAT(//1X,'***          PHYSICAL / THERMODYNAMICAL SPECIFICATION
     X OF EACH SPECIES:'////)                                           
  220 FORMAT(1X,'*** TOTAL MASS: (KG)                  ',3(1PE15.5)/)   
  230 FORMAT(1X,'*** TRANSLATIONAL TEMPERATURE: (K)       ',F10.5/)     
  240 FORMAT(1X,'*** ROTATIONAL TEMPERATURE   : (K)       ',F10.5/)     
  250 FORMAT(1X,'*** MOMENT OF INERTIA: (KG*M**2)      ',3(1PE15.5)/)   
  260 FORMAT(1X,'*** DIPOLE MOMENT    : (C*M)          ',3(1PE15.5)/)   
  270 FORMAT(1X,'                       (DEBYE)        ',3(1PE15.5)/)   
C 280 FORMAT(1X,'*** QUADRUPOLE MOMENT: (C*M**2)       ',3(1PE15.5)/)   
C 285 FORMAT(1X,'                                      ',3(1PE15.5)/)   
  290 FORMAT(1X,'***      MASS (KG)           CHARGE (E)    X (NM)      
     X Y (NM)      Z (NM)'/)                                            
C 300 FORMAT(1X,'***',1PE15.5,0PE20.5,9P3(2X,F10.6))                    
  300 FORMAT(1X,'***',1PE15.5,0PE20.5,3(2X,9PF10.6))                    
  310 FORMAT(//1X,'*** MOVING PARTICLES: ',I5,' *** RIGID POINTS: ',I8/)
CX40066 CALL SBOUTX (66)                                                
CY40066 CALL SBOUTY (66)                                                
      RETURN                                                            
  900 WRITE(66,960)                                                     
CX40066 CALL SBOUTX (66)                                                
CY40066 CALL SBOUTY (66)                                                
      STOP                                                              
  910 WRITE(66,970) NSP,NOP                                             
CX40066 CALL SBOUTX (66)                                                
CY40066 CALL SBOUTY (66)                                                
      STOP                                                              
  960 FORMAT(/1X,']]] ONLY 50 IONS POSSIBLE IN THIS VERSION ')          
  970 FORMAT(/1X,']]] WRONG NR OF PARTICLES ( = ',I5,') - SHOULD BE',I5)
CX40066 CALL SBOUTX (66)                                                
CY40066 CALL SBOUTY (66)                                                
      END                                                               
C================= TYPE  K.N.S. 3/21 84 === REV.001K.N.S. 4/26 84 ======
C                                                                       
      SUBROUTINE TYPE(NNN)                                              
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/CONSTS/ AVSNO,BOLTZ,PLANCK,PI,ELECHG,EPS0                  
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      COMMON/MOUNTE/ NY1,NY2,NY3,KONFS                                  
      COMMON/CHARLO/ TRANP,ROTP,TEMP                                    
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
CX30067 CALL SBINX (67)                                                 
CY30067 CALL SBINY (67)                                                 
      IF(NNN.NE.1) GO TO 500                                            
      WRITE(66,1)                                                       
      WRITE(66,1)                                                       
      WRITE(66,3)                                                       
      WRITE(66,1)                                                       
      WRITE(66,1)                                                       
    1 FORMAT(1X,125('*'))                                               
    2 FORMAT(1X,62('*'),'  BY  SWAMY  KANDADAI    IBM  POUGHKEEPSIE/KING
     XSTON 1984 ','   ***')                                             
    3 FORMAT(1X,62('*'),'*************    MCMOLDYN BIO/1     ***********
     X**********','******')                                             
    5 FORMAT(///)                                                       
      WRITE(66,10)                                                      
   10 FORMAT(1X,15('*'),' M O L E C U L A R  D Y N A M I C S  FOR  BIOMO
     XLECULAR SYSTEMS               ',14X,20('*'))                      
      WRITE(66,1)                                                       
      WRITE(66,1)                                                       
      WRITE(66,2)                                                       
      WRITE(66,1)                                                       
      WRITE(66,1)                                                       
  500 IF(NNN.NE.2) GO TO 1000                                           
 1000 IF(NNN.NE.3) GO TO 1500                                           
      WRITE(66,1)                                                       
      WRITE(66,1)                                                       
      WRITE(66,5)                                                       
      WRITE(66,1009)                                                    
      WRITE(66,1010)                                                    
      WRITE(66,1011)                                                    
      WRITE(66,1012)                                                    
      WRITE(66,1013)                                                    
      WRITE(66,1014)                                                    
      WRITE(66,1015)                                                    
      WRITE(66,1016)                                                    
 1009 FORMAT(/24X,' P R O G R A M     I S    W O R K I N G      O N'//) 
 1010 FORMAT(1X,'   ',24X,'      MMM   MMM           DDDDDDD         ') 
 1011 FORMAT(1X,'   ',24X,'     MMM MMM MMM          DDDDDDDDD       ') 
 1012 FORMAT(1X,'   ',24X,'    MMM   M   MMM         DDD    DDD      ') 
 1013 FORMAT(1X,'   ',24X,'   MMM         MMM        DDD    DDD      ') 
 1014 FORMAT(1X,'   ',24X,'  MMM           MMM       DDDDDDDDD       ') 
 1015 FORMAT(1X,'   ',24X,' MMMMM         MMMMM      DDDDDDD         ') 
 1016 FORMAT(/1X,'   ',37X,'   M   O   D   E                       '//) 
      WRITE(66,1)                                                       
      WRITE(66,1)                                                       
      WRITE(66,1350) NSTEPS,DT                                          
      IF(LTS)WRITE(66,1360) NSSTS                                       
 1350 FORMAT(//1X,'*** NR OF STEPS PLANNED FOR THIS RUN IS ',I8,'       
     X AND THE TIME STEP TO BE USED IS     ',15PF10.4,' FS')            
 1360 FORMAT(/1X,'*** TEMPERATURE SCALING WILL BE USED UNTIL ',I5,'TH   
     XTIME STEP')                                                       
      IF(LLF) WRITE(66,1400)                                            
      IF(LPC) WRITE(66,1410)                                            
      IF(LQ2PC) WRITE(66,1420)                                          
      WRITE(66,1430)                                                    
      WRITE(66,1)                                                       
      WRITE(66,1)                                                       
      WRITE(66,5)                                                       
 1400 FORMAT(//1X,'***  BOTH TRANSLATIONAL AND ROTATIONAL EQUATIONS OF  
     XMOTION ARE SOLVED USING LEAPFROG ALGORITHMS')                     
 1410 FORMAT(//1X,'***  TRANSLATIONAL EQUATIONS OF MOTION ARE SOLVED USI
     XNG 6TH ORDER PREDICTOR-CORRECTOR ALGORITHM')                      
 1420 FORMAT(/1X,'***  ROTATIONAL    EQUATIONS OF MOTION ARE SOLVED USIN
     XG SECOND ORDER QUATERNION METHOD ***********')                    
 1430 FORMAT(/1X,'***  MOLECULAR REORIENTATIONS ARE DESCRIBED USING QUAT
     XERNION ALGEBRA'///)                                               
 1500 IF(NNN.NE.4) GO TO 2000                                           
      CALL CALAVS(NAVRS)                                                
 2000 CONTINUE                                                          
CX40068 CALL SBOUTX (68)                                                
CY40068 CALL SBOUTY (68)                                                
      RETURN                                                            
      END                                                               
C================= CALAVSK.N.S. 3/22 84 === REV.001K.N.S. 4/26 84 ======
C                                                                       
      SUBROUTINE CALAVS(NAVRS)                                          
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      CHARACTER*4 J,K,L,M,N                                             
      COMMON /AVERIZ/ V(13),VV(13)                                      
      DIMENSION AV(13),AVV(13)                                          
      DIMENSION J(13),K(13),L(13),M(13),N(13)                           
      SAVE J, K, L, M, N                                                
      DATA J(1),K(1),L(1),M(1),N(1)/'WW  ','POT.','EN. ','(KJ/','MOL)'/ 
     X J(2),K(2),L(2),M(2),N(2)/'WI  ','POT.','EN. ','(KJ/','MOL)'/     
     X J(3),K(3),L(3),M(3),N(3)/'WGA ','POT.','EN. ','(KJ/','MOL)'/     
     X J(4),K(4),L(4),M(4),N(4)/'II  ','POT.','EN. ','(KJ/','MOL)'/     
     X J(5),K(5),L(5),M(5),N(5)/'IGA ','POT.','EN. ','(KJ/','MOL)'/     
     X J(6),K(6),L(6),M(6),N(6)/'TOT.','POT.','EN. ','(KJ/','MOL)'/     
     X J(7),K(7),L(7),M(7),N(7)/'H2O ','TR. ',' KE ','(KJ/','MOL)'/     
     X J(8),K(8),L(8),M(8),N(8)/'H2O ','ROT.',' KE ','(KJ/','MOL)'/     
     X J(9),K(9),L(9),M(9),N(9)/'ION ','TR. ',' KE ','(KJ/','MOL)'/     
     XJ(10),K(10),L(10),M(10),N(10)/'TOT-','ENER','GY  ','(KJ/','MOL)'/ 
     XJ(11),K(11),L(11),M(11),N(11)/'H2O ','TR. ','TEMP','  ( ',' K )'/ 
     XJ(12),K(12),L(12),M(12),N(12)/'H2O ','ROT ','TEMP','  ( ',' K )'/ 
     XJ(13),K(13),L(13),M(13),N(13)/'ION ','TR. ','TEMP','  ( ',' K )'/ 
CX30069 CALL SBINX (69)                                                 
CY30069 CALL SBINY (69)                                                 
      FNAVS=DBLE(NAVRS)                                                 
CX10723 CALL DOINX( 723 )                                               
CY10723 CALL DOINY( 723 )                                               
      DO 100 I=1,13                                                     
      AV(I)=V(I)/FNAVS                                                  
      AVV(I)=SQRT(ABS(VV(I)/FNAVS-AV(I)*AV(I)))                         
  100 CONTINUE                                                          
CX20724 CALL DOOUTX( 724 )                                              
CY20724 CALL DOOUTY( 724 )                                              
      WRITE(66,150) NAVRS                                               
  150 FORMAT(//1X,'*** AVERAGE VALUES AFTER ',I8,' STEPS:'//)           
CX10725 CALL DOINX( 725 )                                               
CY10725 CALL DOINY( 725 )                                               
      DO 200 I=1,10                                                     
  200 WRITE(66,250) J(I),K(I),L(I),M(I),N(I),AV(I),AVV(I)               
CX20726 CALL DOOUTX( 726 )                                              
CY20726 CALL DOOUTY( 726 )                                              
  250 FORMAT(/1X,5A4,4X,-3PF12.4,' +/-  ',-3PF12.4)                     
      WRITE(1,2)                                                        
      WRITE(1,5) AV(10)                                                 
C AV(10) -  THE AVERAGE ENERGY - IS OUTPUT AS                           
C KILOJOULES PER MOLE AND IS STORED AS JOULES PER MOLE                  
C INTERNALLY. NOTE THE -3P FORMAT SPECIFER.                             
      IF( AV(10) .GT. -7495500.0   .AND.   AV(10) .LT. -7495300.0 ) THEN
        WRITE(1,3)                                                      
      ELSE                                                              
        WRITE(1,4)                                                      
      ENDIF                                                             
CX10727 CALL DOINX( 727 )                                               
CY10727 CALL DOINY( 727 )                                               
C                                                                       
      DO 300 I=11,13                                                    
  300 WRITE(66,350) J(I),K(I),L(I),M(I),N(I),AV(I),AVV(I)               
CX20728 CALL DOOUTX( 728 )                                              
CY20728 CALL DOOUTY( 728 )                                              
  350 FORMAT(/1X,5A4,4X, 0PF12.2,' +/-  ', 0PF12.2)                     
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,-3PF12.4)                                               
CX40070 CALL SBOUTX (70)                                                
CY40070 CALL SBOUTY (70)                                                
      RETURN                                                            
      END                                                               
C================= SCRTCHK.N.S. 3/22 84 === REV.001K.N.S. 4/26 84 ======
C                                                                       
      SUBROUTINE SCRTCH(QA0,QB0,QC0,QD0,WX0,WY0,WZ0,AX0,AY0,AZ0)        
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION QA0(*),QB0(*),QC0(*),QD0(*),WX0(*),WY0(*),WZ0(*)        
     X         ,AX0(*),AY0(*),AZ0(*)                                    
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
      COMMON/CONSTS/ AVSNO,BOLTZ,PLANCK,PI,ELECHG,EPS0                  
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      COMMON /RANDNO/ R3(127),R1,I2                                     
      RANDOM(LY)=2.D0*RAND(0)-1.D0+0.0001D0                             
CX30071 CALL SBINX (71)                                                 
CY30071 CALL SBINY (71)                                                 
C      XXXXX=RAND(1)                                                    
      LY=0                                                              
      IF(LSPA) GO TO 20                                                 
      I=0                                                               
CX10729 CALL DOINX( 729 )                                               
CY10729 CALL DOINY( 729 )                                               
      DO 11 N=1,NTYPES                                                  
      NS=NSITES(N)                                                      
C      FAC=1.D0                                                         
      IF(NS.EQ.1) GO TO 11                                              
      NSP=NSPECI(N)                                                     
      DO 10 J=1,NSP                                                     
      I=I+1                                                             
      QA0(I)=RANDOM(LY)                                                 
      QB0(I)=RANDOM(LY)                                                 
      QC0(I)=RANDOM(LY)                                                 
      QD0(I)=RANDOM(LY)                                                 
      QNORM=SQRT(QA0(I)**2+QB0(I)**2+QC0(I)**2+QD0(I)**2)               
      QA0(I)=QA0(I)/QNORM                                               
      QB0(I)=QB0(I)/QNORM                                               
      QC0(I)=QC0(I)/QNORM                                               
      QD0(I)=QD0(I)/QNORM                                               
   10 CONTINUE                                                          
   11 CONTINUE                                                          
CX20730 CALL DOOUTX( 730 )                                              
CY20730 CALL DOOUTY( 730 )                                              
   20 CONTINUE                                                          
      I=0                                                               
CX10731 CALL DOINX( 731 )                                               
CY10731 CALL DOINY( 731 )                                               
      DO 100 N=1,NTYPES                                                 
      NSP=NSPECI(N)                                                     
      FSP=DBLE(NSP)                                                     
      NS=NSITES(N)                                                      
      BE=BOLTZ/UNITE                                                    
      TT=BE*TTRAN(N)                                                    
      IF(LLF) TERMT=DSUMM(N)                                            
      IF(LPC) TERMT=TSTEP**2/DSUMM(N)                                   
      W =2.D0*SQRT(TERMT*TT)                                            
      DO 30 J=1,NSP                                                     
      I=I+1                                                             
      WX0(I)=W*RANDOM(LY)                                               
      WY0(I)=W*RANDOM(LY)                                               
      WZ0(I)=W*RANDOM(LY)                                               
   30 CONTINUE                                                          
      I=I-NSP                                                           
      W1=SSUM(NSP,WX0(I+1),1)/FSP                                       
      W2=SSUM(NSP,WY0(I+1),1)/FSP                                       
      W3=SSUM(NSP,WZ0(I+1),1)/FSP                                       
      WS=0.D0                                                           
      DO 40 J=1,NSP                                                     
      I=I+1                                                             
      WX0(I)=WX0(I)-W1                                                  
      WY0(I)=WY0(I)-W2                                                  
      WZ0(I)=WZ0(I)-W3                                                  
      WS=WS+WX0(I)**2+WY0(I)**2+WZ0(I)**2                               
   40 CONTINUE                                                          
      IF(.NOT.LLF) GO TO 50                                             
      WS=0.5D0*W/SQRT(WS/(3.D0*FSP))                                    
      GO TO 60                                                          
   50 IF(.NOT.LPC) GO TO 60                                             
      WS=SQRT(3.D0*TSTEP*TSTEP*FSP*TT/WS/DSUMM(N))                      
   60 CONTINUE                                                          
      I=I-NSP                                                           
      DO 70 J=1,NSP                                                     
      I=I+1                                                             
      WX0(I)=WS*WX0(I)                                                  
      WY0(I)=WS*WY0(I)                                                  
      WZ0(I)=WS*WZ0(I)                                                  
   70 CONTINUE                                                          
  100 CONTINUE                                                          
CX20732 CALL DOOUTX( 732 )                                              
CY20732 CALL DOOUTY( 732 )                                              
      WRITE(66,110)                                                     
  110 FORMAT(/1X,'*** START LINEAR VELOCITIES RANDOMIZED')              
      I=0                                                               
CX10733 CALL DOINX( 733 )                                               
CY10733 CALL DOINY( 733 )                                               
      DO 200 N=1,NTYPES                                                 
      IS=NSITES(N)                                                      
      NSP=NSPECI(N)                                                     
      IF(IS.EQ.1) GO TO 200                                             
      NSP=NSPECI(N)                                                     
      FSP=DBLE(NSP)                                                     
      BE=BOLTZ/UNITE                                                    
      TR=BE*TROT (N)                                                    
      AX=2.D0*SQRT(DRI(1,N)*TR)                                         
      AY=2.D0*SQRT(DRI(2,N)*TR)                                         
      AZ=2.D0*SQRT(DRI(3,N)*TR)                                         
      IF(.NOT.LPC) GO TO 125                                            
      AX=AX/DRI(1,N)                                                    
      AY=AY/DRI(2,N)                                                    
      AZ=AZ/DRI(3,N)                                                    
  125 CONTINUE                                                          
      DO 130 J=1,NSP                                                    
      I=I+1                                                             
      AX0(I)=AX*RANDOM(LY)                                              
      AY0(I)=AY*RANDOM(LY)                                              
      AZ0(I)=AZ*RANDOM(LY)                                              
  130 CONTINUE                                                          
      IF(NSP.EQ.1) GO TO 200                                            
      I=I-NSP                                                           
      A1=SSUM(NSP,AX0(I+1),1)/FSP                                       
      A2=SSUM(NSP,AY0(I+1),1)/FSP                                       
      A3=SSUM(NSP,AZ0(I+1),1)/FSP                                       
      ASX=0.D0                                                          
      ASY=0.D0                                                          
      ASZ=0.D0                                                          
      DO 140 J=1,NSP                                                    
      I=I+1                                                             
      AX0(I)=AX0(I)-A1                                                  
      AY0(I)=AY0(I)-A2                                                  
      AZ0(I)=AZ0(I)-A3                                                  
      ASX=ASX+AX0(I)**2                                                 
      ASY=ASY+AY0(I)**2                                                 
      ASZ=ASZ+AZ0(I)**2                                                 
  140 CONTINUE                                                          
      IF(LLF) THEN                                                      
      ASX=0.5D0*AX/SQRT(ASX/FSP )                                       
      ASY=0.5D0*AY/SQRT(ASY/FSP )                                       
      ASZ=0.5D0*AZ/SQRT(ASZ/FSP )                                       
      ELSE                                                              
      ASX=SQRT(FSP*TR/DRI(1,N)/ASX)                                     
      ASY=SQRT(FSP*TR/DRI(2,N)/ASY)                                     
      ASZ=SQRT(FSP*TR/DRI(3,N)/ASZ)                                     
      ENDIF                                                             
      I=I-NSP                                                           
      DO 170 J=1,NSP                                                    
      I=I+1                                                             
C     AX0(I)=ASX*AX0(I)*BREAK                                           
C     AY0(I)=ASY*AY0(I)*BREAK                                           
C     AZ0(I)=ASZ*AZ0(I)*BREAK                                           
      AX0(I)=ASX*AX0(I)                                                 
      AY0(I)=ASY*AY0(I)                                                 
      AZ0(I)=ASZ*AZ0(I)                                                 
  170 CONTINUE                                                          
  200 CONTINUE                                                          
CX20734 CALL DOOUTX( 734 )                                              
CY20734 CALL DOOUTY( 734 )                                              
      WRITE(66,210)                                                     
      WRITE(66,220)                                                     
  210 FORMAT(/1X,'*** START ANGULAR VELOCITIES RANDOMIZED')             
  220 FORMAT(/1X,'*** VELOCITIES SCALED TO TEMPERATURES')               
CX40072 CALL SBOUTX (72)                                                
CY40072 CALL SBOUTY (72)                                                
      RETURN                                                            
      END                                                               
C================= RESTARK.N.S. 3/24 84 === REV.001  S.K. 2/04 86 ======
C                                                                       
      SUBROUTINE RESTAR(T,INPUT,LU)                                     
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION T(*)                                                    
      CHARACTER*4 INTEGR,TYPE,INTE,GRATOR                               
      LOGICAL INPUT,ICHECK                                              
      DIMENSION INTEGR(3),TYPE(3)                                       
      SAVE INTEGR, TYPE                                                 
      DIMENSION OTTRAN(3),OTROT(3)                                      
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      COMMON /AVERIZ/ V(13),VV(13)                                      
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
      COMMON/WINDEX/ NOR,IQ, IW, IA,IM,IF,IT,IP,IO,IX(99)               
      DATA INTEGR/'LEAP','PRED','MONT'/                                 
      DATA TYPE  /'FROG','CORR','CARL'/                                 
CX30073 CALL SBINX (73)                                                 
CY30073 CALL SBINY (73)                                                 
      IF(LLF) ITYPE=1                                                   
      IF(LPC) ITYPE=2                                                   
      NOR=IM-1                                                          
      IF(.NOT.INPUT) GO TO 500                                          
C     REWIND LU                                                         
C     READ(LU) INTE,GRATOR                                              
      READ(LU,870) INTE,GRATOR                                          
      IF(INTE.NE.INTEGR(ITYPE)) WRITE(66,800) INTE,GRATOR               
C     READ(LU) OTTRAN,OTROT,OLDT,NNSTEP,NNOP,NNOR,NNTYPE                
      READ(LU,871) OTTRAN,OTROT                                         
      READ(LU,8711) OLDT,NNSTEP,NNOP,NNOR,NNTYPE                        
      IF(NNTYPE.NE.NTYPES) GO TO 910                                    
      IF  (NNOP.NE.NOP) GO TO 920                                       
      IF  (NNOR.NE.NOR) WRITE(66,810)                                   
CX10735 CALL DOINX( 735 )                                               
CY10735 CALL DOINY( 735 )                                               
      DO 10 I=1,3                                                       
      JBEG=IX(I)                                                        
      JEND=IX(I+1)-1                                                    
C  10 READ(LU) (T(J),J=JBEG,JEND)                                       
   10 READ(LU,872) (T(J),J=JBEG,JEND)                                   
CX20736 CALL DOOUTX( 736 )                                              
CY20736 CALL DOOUTY( 736 )                                              
CX10737 CALL DOINX( 737 )                                               
CY10737 CALL DOINY( 737 )                                               
      DO 15 I=4,NNOR                                                    
      JBEG=IX(I)                                                        
      JEND=IX(I+1)-1                                                    
C  15 READ(LU) (T(J),J=JBEG,JEND)                                       
   15 READ(LU,872) (T(J),J=JBEG,JEND)                                   
CX20738 CALL DOOUTX( 738 )                                              
CY20738 CALL DOOUTY( 738 )                                              
C     READ(LU) NAV,V,VV                                                 
      READ(LU,873) NAV                                                  
      READ(LU,872) V                                                    
      READ(LU,872) VV                                                   
      NAVRS=NAV                                                         
C                                                                       
  870 FORMAT(1X,A4,1X,A4)                                               
  871 FORMAT(1X,6(D12.5,1X))                                            
 8711 FORMAT(1X,D12.5,1X,4(I5,2X))                                      
  872 FORMAT(3(D21.14,3X))                                              
  873 FORMAT(I7)                                                        
C                                                                       
      IF(INTE.NE.INTEGR(ITYPE)) GO TO 300                               
      SDT=DT/OLDT                                                       
      IF(LPC) SDT=TSTEP/OLDT                                            
      ICHECK=.FALSE.                                                    
CX10739 CALL DOINX( 739 )                                               
CY10739 CALL DOINY( 739 )                                               
      DO 20 N=1,NTYPES                                                  
      IF(OTTRAN(N).NE.TTRAN(N)) ICHECK=.TRUE.                           
   20 IF (OTROT(N).NE.TROT(N))  ICHECK=.TRUE.                           
CX20740 CALL DOOUTX( 740 )                                              
CY20740 CALL DOOUTY( 740 )                                              
      IF(ICHECK) THEN                                                   
CX10741 CALL DOINX( 741 )                                               
CY10741 CALL DOINY( 741 )                                               
      DO 30 N=1,NTYPES                                                  
      TTS(N)=OTTRAN(N)                                                  
   30 RTS(N)=OTROT(N)                                                   
CX20742 CALL DOOUTX( 742 )                                              
CY20742 CALL DOOUTY( 742 )                                              
      CALL SCALET(T(IX(IW)),T(IX(IW+1)),T(IX(IW+2))                     
     X           ,T(IX(IA)),T(IX(IA+1)),T(IX(IA+2)))                    
      ENDIF                                                             
      ICHECK=.FALSE.                                                    
      IF(OLDT.NE.DT) ICHECK=.TRUE.                                      
      IF(ICHECK) THEN                                                   
CX10743 CALL DOINX( 743 )                                               
CY10743 CALL DOINY( 743 )                                               
      DO 40 I=0,NOP-1                                                   
      T(IX(IW+0)+I)=T(IX(IW+0)+I)*SDT                                   
      T(IX(IW+1)+I)=T(IX(IW+1)+I)*SDT                                   
      T(IX(IW+2)+I)=T(IX(IW+2)+I)*SDT                                   
   40 CONTINUE                                                          
CX20744 CALL DOOUTX( 744 )                                              
CY20744 CALL DOOUTY( 744 )                                              
      NSP = NSPECI(1)                                                   
CX10745 CALL DOINX( 745 )                                               
CY10745 CALL DOINY( 745 )                                               
      DO 41 I=0,NSP-1                                                   
      T(IX(IA+0)+I)=T(IX(IA+0)+I)*SDT                                   
      T(IX(IA+1)+I)=T(IX(IA+1)+I)*SDT                                   
      T(IX(IA+2)+I)=T(IX(IA+2)+I)*SDT                                   
   41 CONTINUE                                                          
CX20746 CALL DOOUTX( 746 )                                              
CY20746 CALL DOOUTY( 746 )                                              
      ENDIF                                                             
      IF(.NOT.LPC) GO TO 210                                            
C*** SET INITIAL POSITIONS TO COMPUTE M.S.D FOR PRED-CORR ***           
CX10747 CALL DOINX( 747 )                                               
CY10747 CALL DOINY( 747 )                                               
      DO 45 J=0,NOP-1                                                   
      T(IX(IM+0)+J)=T(IX(1)+J)                                          
      T(IX(IM+1)+J)=T(IX(2)+J)                                          
      T(IX(IM+2)+J)=T(IX(3)+J)                                          
   45 CONTINUE                                                          
CX20748 CALL DOOUTX( 748 )                                              
CY20748 CALL DOOUTY( 748 )                                              
      SDT2=SDT*SDT                                                      
      SDT3=SDT*SDT2                                                     
      SDT4=SDT*SDT3                                                     
      SDT5=SDT*SDT4                                                     
      SDT6=SDT*SDT5                                                     
      I=0                                                               
CX10749 CALL DOINX( 749 )                                               
CY10749 CALL DOINY( 749 )                                               
      DO 200 N=1,NTYPES                                                 
      NSP=NSPECI(N)                                                     
      DO 150 J=1,NSP                                                    
      DO 70 M=3,5                                                       
      T(IX(IW+M)+I)=T(IX(IW+M)+I)*SDT2                                  
   70 CONTINUE                                                          
      DO 75 M=6,8                                                       
      T(IX(IW+M)+I)=T(IX(IW+M)+I)*SDT3                                  
   75 CONTINUE                                                          
      DO 80 M=9,11                                                      
      T(IX(IW+M)+I)=T(IX(IW+M)+I)*SDT4                                  
   80 CONTINUE                                                          
      DO 85 M=12,14                                                     
      T(IX(IW+M)+I)=T(IX(IW+M)+I)*SDT5                                  
   85 CONTINUE                                                          
      I=I+1                                                             
  150 CONTINUE                                                          
  200 CONTINUE                                                          
CX20750 CALL DOOUTX( 750 )                                              
CY20750 CALL DOOUTY( 750 )                                              
      I=0                                                               
      NSP=NSPECI(1)                                                     
CX10751 CALL DOINX( 751 )                                               
CY10751 CALL DOINY( 751 )                                               
      DO 700 J=1,NSP                                                    
      DO 701 M=3,5                                                      
  701 T(IX(IA+M)+I)=T(IX(IA+M)+I)*SDT                                   
      DO 702 M=6,8                                                      
  702 T(IX(IA+M)+I)=T(IX(IA+M)+I)*SDT2                                  
      DO 703 M=9,11                                                     
  703 T(IX(IA+M)+I)=T(IX(IA+M)+I)*SDT3                                  
      DO 704 M=12,14                                                    
  704 T(IX(IA+M)+I)=T(IX(IA+M)+I)*SDT4                                  
      DO 705 M=15,17                                                    
  705 T(IX(IA+M)+I)=T(IX(IA+M)+I)*SDT5                                  
      DO 706 M=18,20                                                    
      T(IX(IA+M)+I)=T(IX(IA+M)+I)*SDT6                                  
  706 CONTINUE                                                          
      DO 100 M=4,7                                                      
      T(IX(IQ+M)+I)=T(IX(IQ+M)+I)*SDT                                   
  100 CONTINUE                                                          
      DO 105 M=8,11                                                     
      T(IX(IQ+M)+I)=T(IX(IQ+M)+I)*SDT2                                  
  105 CONTINUE                                                          
      DO 110 M=12,15                                                    
      T(IX(IQ+M)+I)=T(IX(IQ+M)+I)*SDT3                                  
  110 CONTINUE                                                          
      DO 115 M=16,19                                                    
      T(IX(IQ+M)+I)=T(IX(IQ+M)+I)*SDT4                                  
  115 CONTINUE                                                          
      DO 120 M=20,23                                                    
      T(IX(IQ+M)+I)=T(IX(IQ+M)+I)*SDT5                                  
  120 CONTINUE                                                          
      DO 125 M=24,27                                                    
      T(IX(IQ+M)+I)=T(IX(IQ+M)+I)*SDT6                                  
  125 CONTINUE                                                          
      I=I+1                                                             
  700 CONTINUE                                                          
CX20752 CALL DOOUTX( 752 )                                              
CY20752 CALL DOOUTY( 752 )                                              
  210 CONTINUE                                                          
      WRITE(66,250) INTE,GRATOR                                         
  250 FORMAT(//1X,'*** ',A4,'-',A4,' RESTART CONFIGURATION READ IN '///)
CX40074 CALL SBOUTX (74)                                                
CY40074 CALL SBOUTY (74)                                                
      RETURN                                                            
  300 CONTINUE                                                          
      LRES=.FALSE.                                                      
      LSPA=.TRUE.                                                       
      IF(INTE.EQ.INTEGR(2)) GO TO 330                                   
CX10753 CALL DOINX( 753 )                                               
CY10753 CALL DOINY( 753 )                                               
      DO 320 I=0,NOP-1                                                  
  320 T(IX(IQ)+I)=T(IX(4)+I)                                            
CX20754 CALL DOOUTX( 754 )                                              
CY20754 CALL DOOUTY( 754 )                                              
CX40074 CALL SBOUTX (74)                                                
CY40074 CALL SBOUTY (74)                                                
      RETURN                                                            
  330 CONTINUE                                                          
CX10755 CALL DOINX( 755 )                                               
CY10755 CALL DOINY( 755 )                                               
      DO 340 I=0,NOP-1                                                  
  340 T(IX(IQ)+I)=T(IX(19)+I)                                           
CX20756 CALL DOOUTX( 756 )                                              
CY20756 CALL DOOUTY( 756 )                                              
CX40074 CALL SBOUTX (74)                                                
CY40074 CALL SBOUTY (74)                                                
      RETURN                                                            
  500 CONTINUE                                                          
      WRITE(LU) INTEGR(ITYPE),TYPE(ITYPE)                               
      IF(LLF)WRITE(LU) TTRAN,TROT,DT   ,NSTEP,NOP,NOR,NTYPES            
      IF(LPC)WRITE(LU) TTRAN,TROT,TSTEP,NSTEP,NOP,NOR,NTYPES            
CX10757 CALL DOINX( 757 )                                               
CY10757 CALL DOINY( 757 )                                               
      DO 550 I=1,3                                                      
      JBEG=IX(I)                                                        
      JEND=IX(I+1)-1                                                    
  550 WRITE(LU) (T(J),J=JBEG,JEND)                                      
CX20758 CALL DOOUTX( 758 )                                              
CY20758 CALL DOOUTY( 758 )                                              
CX10759 CALL DOINX( 759 )                                               
CY10759 CALL DOINY( 759 )                                               
      DO 560 I=4,IM-1                                                   
      JBEG=IX(I)                                                        
      JEND=IX(I+1)-1                                                    
  560 WRITE(LU) (T(J),J=JBEG,JEND)                                      
CX20760 CALL DOOUTX( 760 )                                              
CY20760 CALL DOOUTY( 760 )                                              
      WRITE(LU) NAVRS,V,VV                                              
      REWIND LU                                                         
      IF(NSTEP.EQ.NSTEPS) WRITE(66,600) INTEGR(ITYPE),TYPE(ITYPE)       
  600 FORMAT(//1X,'*** ',A4,'-',A4,' RESTART CONFIGURATION DUMPED ')    
CX40074 CALL SBOUTX (74)                                                
CY40074 CALL SBOUTY (74)                                                
      RETURN                                                            
  800 FORMAT(//1X,'??? RESTART FILE IS OF ',A4,' - ',A4,' TYPE'/)       
  810 FORMAT(//1X,'??? NR OF RECORDS IN RESTART FILE DOES NOT FIT THE   
     XACTUAL CASE'/)                                                    
C 900 WRITE(66,950)                                                     
C     STOP                                                              
  910 WRITE(66,960) NNTYPE,NTYPES                                       
CX40074 CALL SBOUTX (74)                                                
CY40074 CALL SBOUTY (74)                                                
      STOP                                                              
  920 WRITE(66,970) NNOP,NOP                                            
CX40074 CALL SBOUTX (74)                                                
CY40074 CALL SBOUTY (74)                                                
      STOP                                                              
  950 FORMAT(//1X,']]] WRONG RESTART FILE ')                            
  960 FORMAT(//1X,']]] CANT READ THE RESTART FILE  -   NR OF MOLECULE TY
     XPES IN THE RESTART FILE IS ',I8,'    - - - SHOULD BE ',I8)        
  970 FORMAT(//1X,']]] CANT READ THE RESTART FILE  -  NR OF PARTICLES IN
     X RESTART FILE IS ',I5,' WHILE THE CURRENT NR IS ',I5/)            
CX40074 CALL SBOUTX (74)                                                
CY40074 CALL SBOUTY (74)                                                
      END                                                               
C****************************************************************       
      SUBROUTINE DUMP(X0,Y0,Z0,X1,Y1,Z1,SX,SY,SZ,WX,WY,WZ,              
     *Q1,Q2,Q3,Q4)                                                      
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      DIMENSION X0(*),Y0(*),Z0(*),SX(*),SY(*),SZ(*)                     
      DIMENSION X1(*),Y1(*),Z1(*),WX(*),WY(*),WZ(*)                     
      DIMENSION Q1(*),Q2(*),Q3(*),Q4(*)                                 
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
CX30075 CALL SBINX (75)                                                 
CY30075 CALL SBINY (75)                                                 
      NSP1=NSPECI(1)                                                    
      LLU=20                                                            
      WRITE(LLU) NSTEP                                                  
      WRITE(LLU) (X0(I),I=1,NOP)                                        
      WRITE(LLU) (Y0(I),I=1,NOP)                                        
      WRITE(LLU) (Z0(I),I=1,NOP)                                        
      WRITE(LLU) (SX(I),I=1,2000)                                       
      WRITE(LLU) (SX(I),I=2001,4000)                                    
      WRITE(LLU) (SX(I),I=4001,6000)                                    
      WRITE(LLU) (SY(I),I=1,2000)                                       
      WRITE(LLU) (SY(I),I=2001,4000)                                    
      WRITE(LLU) (SY(I),I=4001,6000)                                    
      WRITE(LLU) (SZ(I),I=1,2000)                                       
      WRITE(LLU) (SZ(I),I=2001,4000)                                    
      WRITE(LLU) (SZ(I),I=4001,6000)                                    
      WRITE(LLU) (X1(I),I=1,NOP)                                        
      WRITE(LLU) (Y1(I),I=1,NOP)                                        
      WRITE(LLU) (Z1(I),I=1,NOP)                                        
      WRITE(LLU) (WX(I),I=1,NSP1)                                       
      WRITE(LLU) (WY(I),I=1,NSP1)                                       
      WRITE(LLU) (WZ(I),I=1,NSP1)                                       
      WRITE(LLU) (Q1(I),I=1,NSP1)                                       
      WRITE(LLU) (Q2(I),I=1,NSP1)                                       
      WRITE(LLU) (Q3(I),I=1,NSP1)                                       
      WRITE(LLU) (Q4(I),I=1,NSP1)                                       
CX40076 CALL SBOUTX (76)                                                
CY40076 CALL SBOUTY (76)                                                
      RETURN                                                            
      END                                                               
C================= ACTION  S.K. 9/22 85 === REV.001  S.K.2/04 86 ====== 
C                                                                       
      SUBROUTINE ACTFOR(X0,Y0,Z0,FX,FY,FZ,TX,TY,TZ,SX,SY,SZ             
     X                 ,FSX,FSY,FSZ)                                    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION X0(*),Y0(*),Z0(*),FX(*),FY(*),FZ(*)                     
      DIMENSION TX(*),TY(*),TZ(*),SX(*),SY(*),SZ(*)                     
      DIMENSION FSX(*),FSY(*),FSZ(*)                                    
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON /SCRTZR/A(3),B(3),C(3),AA(3,650),BB(3,650),CC(3,650)       
     X              ,AAA(650),BBB(650),CCC(650),AAAA,BBBB,CCCC          
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
C     COMMON/TABELS/ QQ(11,11),AD(5),BD(5)                              
      COMMON/TABELS/ QQ(11,11),                                         
     1A1,A2,A3,A4,A5,B1,B2,B3,B4,B5                                     
      COMMON/INDEXQ/ NWI,NWJ,NWK,INDEX(11,11),INDS(3,10),NONE           
      DIMENSION AA1(4),AA2(4),BB1(4),BB2(4),SFX(4),SFY(4),SFZ(4)        
      DIMENSION MDEX(4),NDEX(4),JJJ(4),IXX(4)                           
      SAVE MDEX, NDEX                                                   
      DIMENSION FOX(2000),FOY(2000),FOZ(2000),F1X(2000),F1Y(2000),      
     1F1Z(2000),F2X(2000),F2Y(2000),F2Z(2000),                          
     2FPX(2000),FPY(2000),FPZ(2000),                                    
     3FAX(2000),FAY(2000),FAZ(2000),FIX(50),FIY(50),FIZ(50),            
     4FOXP(2000),FOYP(2000),FOZP(2000),F1XP(2000),F1YP(2000),           
     1F1ZP(2000),F2XP(2000),F2YP(2000),F2ZP(2000),                      
     2FPXP(2000),FPYP(2000),FPZP(2000),                                 
     4XDT(2000),YDT(2000),ZDT(2000),IND(2000)                           
      DATA MDEX/5,3,5,5/                                                
      DATA NDEX/1,2,4,5/                                                
CX30077 CALL SBINX (77)                                                 
CY30077 CALL SBINY (77)                                                 
      ALENGT=22.D-10/UNITL                                              
      ALENGM=1.D0/ALENGT                                                
      B1M=-B1                                                           
      B2M=-B2                                                           
      B3M=-B3                                                           
      B4M=-B4                                                           
      EFFE=0.D0                                                         
      EWW=0.D0                                                          
      EWI=0.D0                                                          
      EWA=0.D0                                                          
      EII=0.D0                                                          
      EIA=0.D0                                                          
      NSP=NSPECI(1)                                                     
      NSST=NSS+NION+NATOMS                                              
CX10761 CALL DOINX( 761 )                                               
CY10761 CALL DOINY( 761 )                                               
      DO 10 I=1,NSST                                                    
      FSX(I)=0.D0                                                       
      FSY(I)=0.D0                                                       
      FSZ(I)=0.D0                                                       
  10  CONTINUE                                                          
CX20762 CALL DOOUTX( 762 )                                              
CY20762 CALL DOOUTY( 762 )                                              
C     CALL ZERO(NSST,FSX(1))                                            
C     CALL ZERO(NSST,FSY(1))                                            
C     CALL ZERO(NSST,FSZ(1))                                            
C     CALL ZERO(NOP,FX(1))                                              
C     CALL ZERO(NOP,FY(1))                                              
C     CALL ZERO(NOP,FZ(1))                                              
C     CALL ZERO(NSP,TX(1))                                              
C     CALL ZERO(NSP,TY(1))                                              
C     CALL ZERO(NSP,TZ(1))                                              
      ISIT=NSITES(1)                                                    
C************************************************************           
C     DO 121 I=1,4                                                      
C     IXX(I)=0                                                          
C     JJJ(I)=0                                                          
C     SFX(I)=0.D0                                                       
C     SFY(I)=0.D0                                                       
C     SFZ(I)=0.D0                                                       
C     AA1(I)=0.D0                                                       
C     AA2(I)=0.D0                                                       
C     BB1(I)=0.D0                                                       
C     BB2(I)=0.D0                                                       
C 121 CONTINUE                                                          
C***  DEFINE THE CUTOFF AS OXYGEN-OXYGEN DISTANCE ********              
C***   WE HAVE USED A CUTOFF OF 8.0 ANGSTROMS HERE *******              
      RCUTS=64.D-20/UNITL/UNITL                                         
CX10763 CALL DOINX( 763 )                                               
CY10763 CALL DOINY( 763 )                                               
      DO 100 I=1,NSP                                                    
      FPX(I)=0.D0                                                       
      FPY(I)=0.D0                                                       
      FPZ(I)=0.D0                                                       
      FOX(I)=0.D0                                                       
      FOY(I)=0.D0                                                       
      FOZ(I)=0.D0                                                       
      F1X(I)=0.D0                                                       
      F1Y(I)=0.D0                                                       
      F1Z(I)=0.D0                                                       
      F2X(I)=0.D0                                                       
      F2Y(I)=0.D0                                                       
      F2Z(I)=0.D0                                                       
 100  CONTINUE                                                          
C***********************************************************            
C*****    WATER - WATER INTERACTIONS                   *****            
C***********************************************************            
CX20764 CALL DOOUTX( 764 )                                              
CY20764 CALL DOOUTY( 764 )                                              
CX10765 CALL DOINX( 765 )                                               
CY10765 CALL DOINY( 765 )                                               
      DO 240 I=2,NSP                                                    
      INS=(I-1)*ISIT                                                    
      FXO=0.D0                                                          
      FYO=0.D0                                                          
      FZO=0.D0                                                          
      FX1=0.D0                                                          
      FY1=0.D0                                                          
      FZ1=0.D0                                                          
      FX2=0.D0                                                          
      FY2=0.D0                                                          
      FZ2=0.D0                                                          
      FXP=0.D0                                                          
      FYP=0.D0                                                          
      FZP=0.D0                                                          
      DO 235 J=1,I-1                                                    
      IND(J)=0                                                          
      JNS=(J-1)*ISIT                                                    
      XD=X0(I)-X0(J)                                                    
      YD=Y0(I)-Y0(J)                                                    
      ZD=Z0(I)-Z0(J)                                                    
      XDT(J)=XD-2.D0*ALENGT*DBLE(INT(XD*ALENGM))                        
      YDT(J)=YD-2.D0*ALENGT*DBLE(INT(YD*ALENGM))                        
      ZDT(J)=ZD-2.D0*DBLE(INT(ZD))                                      
C O-O                                                                   
      DXS=XDT(J)+SX(INS+1)                                              
      DYS=YDT(J)+SY(INS+1)                                              
      DZS=ZDT(J)+SZ(INS+1)                                              
      RX=DXS-SX(JNS+1)                                                  
      RY=DYS-SY(JNS+1)                                                  
      RZ=DZS-SZ(JNS+1)                                                  
      RSQ=RX*RX+RY*RY+RZ*RZ                                             
      IF(RSQ.GE.RCUTS) GO TO  235                                       
      IND(J)=1                                                          
 235  CONTINUE                                                          
      L=0                                                               
      DO 236 J=1,I-1                                                    
      IF(IND(J).EQ.0) GO TO 236                                         
      L=L+1                                                             
      IND(L)=J                                                          
 236  CONTINUE                                                          
      IF(L.EQ.0) GO TO 240                                              
      DO 237 J=1,L                                                      
      K=IND(J)                                                          
      JNS=(K-1)*ISIT                                                    
      XD=XDT(K)                                                         
      YD=YDT(K)                                                         
      ZD=ZDT(K)                                                         
C O-O                                                                   
      DXS=XD+SX(INS+1)                                                  
      DYS=YD+SY(INS+1)                                                  
      DZS=ZD+SZ(INS+1)                                                  
      RX=DXS-SX(JNS+1)                                                  
      RY=DYS-SY(JNS+1)                                                  
      RZ=DZS-SZ(JNS+1)                                                  
      RSQ=RX*RX+RY*RY+RZ*RZ                                             
      R=SQRT(RSQ)                                                       
      EWWT  =A1*EXP(B1M*R)                                              
      FOR=(B1*EWWT  )/R                                                 
C     EWW=EWW+ETERM1                                                    
      FFX=RX*FOR                                                        
      FFY=RY*FOR                                                        
      FFZ=RZ*FOR                                                        
      FXOT  =       FFX                                                 
      FYOT  =       FFY                                                 
      FZOT  =       FFZ                                                 
      FOXJ=       FFX                                                   
      FOYJ=       FFY                                                   
      FOZJ=       FFZ                                                   
C O-H                                                                   
      RX=DXS-SX(JNS+2)                                                  
      RY=DYS-SY(JNS+2)                                                  
      RZ=DZS-SZ(JNS+2)                                                  
      RSQ=RX*RX+RY*RY+RZ*RZ                                             
      R=SQRT(RSQ)                                                       
      ETERM1=A2*EXP(B2M*R)                                              
      ETERM2=A3*EXP(B3M*R)                                              
      FOR=(B2*ETERM1+B3*ETERM2)/R                                       
      EWWT=EWWT+ETERM1+ETERM2                                           
      FFX=RX*FOR                                                        
      FFY=RY*FOR                                                        
      FFZ=RZ*FOR                                                        
      FXOT  =FXOT  +FFX                                                 
      FYOT  =FYOT  +FFY                                                 
      FZOT  =FZOT  +FFZ                                                 
      F1XJ=       FFX                                                   
      F1YJ=       FFY                                                   
      F1ZJ=       FFZ                                                   
C O-H                                                                   
      RX=DXS-SX(JNS+3)                                                  
      RY=DYS-SY(JNS+3)                                                  
      RZ=DZS-SZ(JNS+3)                                                  
      RSQ=RX*RX+RY*RY+RZ*RZ                                             
      R=SQRT(RSQ)                                                       
      ETERM1=A2*EXP(B2M*R)                                              
      ETERM2=A3*EXP(B3M*R)                                              
      FOR=(B2*ETERM1+B3*ETERM2)/R                                       
      EWWT=EWWT+ETERM1+ETERM2                                           
      FFX=RX*FOR                                                        
      FFY=RY*FOR                                                        
      FFZ=RZ*FOR                                                        
      FXOT  =FXOT  +FFX                                                 
      FYOT  =FYOT  +FFY                                                 
      FZOT  =FZOT  +FFZ                                                 
      F2XJ=       FFX                                                   
      F2YJ=       FFY                                                   
      F2ZJ=       FFZ                                                   
C H-H                                                                   
      DXS=XD+SX(INS+2)                                                  
      DYS=YD+SY(INS+2)                                                  
      DZS=ZD+SZ(INS+2)                                                  
      RX=DXS-SX(JNS+2)                                                  
      RY=DYS-SY(JNS+2)                                                  
      RZ=DZS-SZ(JNS+2)                                                  
      RSQ=RX*RX+RY*RY+RZ*RZ                                             
      R=SQRT(RSQ)                                                       
      RM1=1.D0/R                                                        
      TEX=QQ(2,3)*RM1                                                   
      EFFET=      TEX                                                   
      FOR=TEX/RSQ                                                       
      ETERM1=A4*EXP(B4M*R)                                              
      FOR=(B4*ETERM1)*RM1+FOR                                           
      EWWT=EWWT+ETERM1+TEX                                              
      FFX=RX*FOR                                                        
      FFY=RY*FOR                                                        
      FFZ=RZ*FOR                                                        
      FX1T  =       FFX                                                 
      FY1T  =       FFY                                                 
      FZ1T  =       FFZ                                                 
      F1XJ=F1XJ+FFX                                                     
      F1YJ=F1YJ+FFY                                                     
      F1ZJ=F1ZJ+FFZ                                                     
C H-H                                                                   
      RX=DXS-SX(JNS+3)                                                  
      RY=DYS-SY(JNS+3)                                                  
      RZ=DZS-SZ(JNS+3)                                                  
      RSQ=RX*RX+RY*RY+RZ*RZ                                             
      R=SQRT(RSQ)                                                       
      RM1=1.D0/R                                                        
      TEX=QQ(2,3)*RM1                                                   
      EFFET=EFFET+TEX                                                   
      FOR=TEX/RSQ                                                       
      ETERM1=A4*EXP(B4M*R)                                              
      FOR=(B4*ETERM1)*RM1+FOR                                           
      EWWT=EWWT+ETERM1+TEX                                              
      FFX=RX*FOR                                                        
      FFY=RY*FOR                                                        
      FFZ=RZ*FOR                                                        
      FX1T  =FX1T  +FFX                                                 
      FY1T  =FY1T  +FFY                                                 
      FZ1T  =FZ1T  +FFZ                                                 
      F2XJ=F2XJ+FFX                                                     
      F2YJ=F2YJ+FFY                                                     
      F2ZJ=F2ZJ+FFZ                                                     
C O-H                                                                   
      RX=DXS-SX(JNS+1)                                                  
      RY=DYS-SY(JNS+1)                                                  
      RZ=DZS-SZ(JNS+1)                                                  
      RSQ=RX*RX+RY*RY+RZ*RZ                                             
      R=SQRT(RSQ)                                                       
      ETERM1=A2*EXP(B2M*R)                                              
      ETERM2=A3*EXP(B3M*R)                                              
      FOR=(B2*ETERM1+B3*ETERM2)/R                                       
      EWWT=EWWT+ETERM1+ETERM2                                           
      FFX=RX*FOR                                                        
      FFY=RY*FOR                                                        
      FFZ=RZ*FOR                                                        
      FX1T  =FX1T  +FFX                                                 
      FY1T  =FY1T  +FFY                                                 
      FZ1T  =FZ1T  +FFZ                                                 
      FOXJ=FOXJ+FFX                                                     
      FOYJ=FOYJ+FFY                                                     
      FOZJ=FOZJ+FFZ                                                     
C H-P                                                                   
      RX=DXS-SX(JNS+4)                                                  
      RY=DYS-SY(JNS+4)                                                  
      RZ=DZS-SZ(JNS+4)                                                  
      RSQ=RX*RX+RY*RY+RZ*RZ                                             
      R=SQRT(RSQ)                                                       
      RM1=1.D0/R                                                        
      TEX=QQ(2,4)*RM1                                                   
      EFFET=EFFET+TEX                                                   
      FOR=TEX/RSQ                                                       
      EWWT=EWWT+TEX                                                     
      FFX=RX*FOR                                                        
      FFY=RY*FOR                                                        
      FFZ=RZ*FOR                                                        
      FX1T  =FX1T  +FFX                                                 
      FY1T  =FY1T  +FFY                                                 
      FZ1T  =FZ1T  +FFZ                                                 
      FPXJ=       FFX                                                   
      FPYJ=       FFY                                                   
      FPZJ=       FFZ                                                   
C H-H                                                                   
      DXS=XD+SX(INS+3)                                                  
      DYS=YD+SY(INS+3)                                                  
      DZS=ZD+SZ(INS+3)                                                  
      RX=DXS-SX(JNS+2)                                                  
      RY=DYS-SY(JNS+2)                                                  
      RZ=DZS-SZ(JNS+2)                                                  
      RSQ=RX*RX+RY*RY+RZ*RZ                                             
      R=SQRT(RSQ)                                                       
      RM1=1.D0/R                                                        
      TEX=QQ(2,3)*RM1                                                   
      EFFET=EFFET+TEX                                                   
      FOR=TEX/RSQ                                                       
      ETERM1=A4*EXP(B4M*R)                                              
      FOR=(B4*ETERM1)*RM1+FOR                                           
      EWWT=EWWT+ETERM1+TEX                                              
      FFX=RX*FOR                                                        
      FFY=RY*FOR                                                        
      FFZ=RZ*FOR                                                        
      FX2T  =       FFX                                                 
      FY2T  =       FFY                                                 
      FZ2T  =       FFZ                                                 
      F1XJ=F1XJ+FFX                                                     
      F1YJ=F1YJ+FFY                                                     
      F1ZJ=F1ZJ+FFZ                                                     
C H-H                                                                   
      RX=DXS-SX(JNS+3)                                                  
      RY=DYS-SY(JNS+3)                                                  
      RZ=DZS-SZ(JNS+3)                                                  
      RSQ=RX*RX+RY*RY+RZ*RZ                                             
      R=SQRT(RSQ)                                                       
      RM1=1.D0/R                                                        
      TEX=QQ(2,3)*RM1                                                   
      EFFET=EFFET+TEX                                                   
      FOR=TEX/RSQ                                                       
      ETERM1=A4*EXP(B4M*R)                                              
      FOR=(B4*ETERM1)*RM1+FOR                                           
      EWWT=EWWT+ETERM1+TEX                                              
      FFX=RX*FOR                                                        
      FFY=RY*FOR                                                        
      FFZ=RZ*FOR                                                        
      FX2T  =FX2T  +FFX                                                 
      FY2T  =FY2T  +FFY                                                 
      FZ2T  =FZ2T  +FFZ                                                 
      F2XJ=F2XJ+FFX                                                     
      F2YJ=F2YJ+FFY                                                     
      F2ZJ=F2ZJ+FFZ                                                     
C O-H                                                                   
      RX=DXS-SX(JNS+1)                                                  
      RY=DYS-SY(JNS+1)                                                  
      RZ=DZS-SZ(JNS+1)                                                  
      RSQ=RX*RX+RY*RY+RZ*RZ                                             
      R=SQRT(RSQ)                                                       
      ETERM1=A2*EXP(B2M*R)                                              
      ETERM2=A3*EXP(B3M*R)                                              
      FOR=(B2*ETERM1+B3*ETERM2)/R                                       
      EWWT=EWWT+ETERM1+ETERM2                                           
      FFX=RX*FOR                                                        
      FFY=RY*FOR                                                        
      FFZ=RZ*FOR                                                        
      FX2T  =FX2T  +FFX                                                 
      FY2T  =FY2T  +FFY                                                 
      FZ2T  =FZ2T  +FFZ                                                 
      FOXJ=FOXJ+FFX                                                     
      FOYJ=FOYJ+FFY                                                     
      FOZJ=FOZJ+FFZ                                                     
C H-P                                                                   
      RX=DXS-SX(JNS+4)                                                  
      RY=DYS-SY(JNS+4)                                                  
      RZ=DZS-SZ(JNS+4)                                                  
      RSQ=RX*RX+RY*RY+RZ*RZ                                             
      R=SQRT(RSQ)                                                       
      RM1=1.D0/R                                                        
      TEX=QQ(2,4)*RM1                                                   
      EFFET=EFFET+TEX                                                   
      FOR=TEX/RSQ                                                       
      EWWT=EWWT+TEX                                                     
      FFX=RX*FOR                                                        
      FFY=RY*FOR                                                        
      FFZ=RZ*FOR                                                        
      FX2T  =FX2T  +FFX                                                 
      FY2T  =FY2T  +FFY                                                 
      FZ2T  =FZ2T  +FFZ                                                 
      FPXJ=FPXJ+FFX                                                     
      FPYJ=FPYJ+FFY                                                     
      FPZJ=FPZJ+FFZ                                                     
C H-P                                                                   
      DXS=XD+SX(INS+4)                                                  
      DYS=YD+SY(INS+4)                                                  
      DZS=ZD+SZ(INS+4)                                                  
      RX=DXS-SX(JNS+2)                                                  
      RY=DYS-SY(JNS+2)                                                  
      RZ=DZS-SZ(JNS+2)                                                  
      RSQ=RX*RX+RY*RY+RZ*RZ                                             
      R=SQRT(RSQ)                                                       
      RM1=1.D0/R                                                        
      TEX=QQ(2,4)*RM1                                                   
      EFFET=EFFET+TEX                                                   
      FOR=TEX/RSQ                                                       
      EWWT=EWWT+TEX                                                     
      FFX=RX*FOR                                                        
      FFY=RY*FOR                                                        
      FFZ=RZ*FOR                                                        
      FXPT  =       FFX                                                 
      FYPT  =       FFY                                                 
      FZPT  =       FFZ                                                 
      F1XJ=F1XJ+FFX                                                     
      F1YJ=F1YJ+FFY                                                     
      F1ZJ=F1ZJ+FFZ                                                     
C H-P                                                                   
      RX=DXS-SX(JNS+3)                                                  
      RY=DYS-SY(JNS+3)                                                  
      RZ=DZS-SZ(JNS+3)                                                  
      RSQ=RX*RX+RY*RY+RZ*RZ                                             
      R=SQRT(RSQ)                                                       
      RM1=1.D0/R                                                        
      TEX=QQ(2,4)*RM1                                                   
      EFFET=EFFET+TEX                                                   
      FOR=TEX/RSQ                                                       
      EWWT=EWWT+TEX                                                     
      FFX=RX*FOR                                                        
      FFY=RY*FOR                                                        
      FFZ=RZ*FOR                                                        
      FXPT  =FXPT  +FFX                                                 
      FYPT  =FYPT  +FFY                                                 
      FZPT  =FZPT  +FFZ                                                 
      F2XJ=F2XJ+FFX                                                     
      F2YJ=F2YJ+FFY                                                     
      F2ZJ=F2ZJ+FFZ                                                     
C P-P                                                                   
      RX=DXS-SX(JNS+4)                                                  
      RY=DYS-SY(JNS+4)                                                  
      RZ=DZS-SZ(JNS+4)                                                  
      RSQ=RX*RX+RY*RY+RZ*RZ                                             
      R=SQRT(RSQ)                                                       
      RM1=1.D0/R                                                        
      TEX=QQ(4,4)*RM1                                                   
      EFFET=EFFET+TEX                                                   
      FOR=TEX/RSQ                                                       
      EWWT=EWWT+TEX                                                     
      EWW =EWWT+EWW                                                     
      EFFE=EFFET+EFFE                                                   
      FFX=RX*FOR                                                        
      FFY=RY*FOR                                                        
      FFZ=RZ*FOR                                                        
      FXPT  =FXPT  +FFX                                                 
      FYPT  =FYPT  +FFY                                                 
      FZPT  =FZPT  +FFZ                                                 
      FPXJ=FPXJ+FFX                                                     
      FPYJ=FPYJ+FFY                                                     
      FPZJ=FPZJ+FFZ                                                     
      FXO   =FXOT  +FXO                                                 
      FYO   =FYOT  +FYO                                                 
      FZO   =FZOT  +FZO                                                 
      FX1   =FX1T  +FX1                                                 
      FY1   =FY1T  +FY1                                                 
      FZ1   =FZ1T  +FZ1                                                 
      FX2   =FX2T  +FX2                                                 
      FY2   =FY2T  +FY2                                                 
      FZ2   =FZ2T  +FZ2                                                 
      FXP   =FXPT  +FXP                                                 
      FYP   =FYPT  +FYP                                                 
      FZP   =FZPT  +FZP                                                 
      FOXP(J)=       FOXJ                                               
      FOYP(J)=       FOYJ                                               
      FOZP(J)=       FOZJ                                               
      F1XP(J)=       F1XJ                                               
      F1YP(J)=       F1YJ                                               
      F1ZP(J)=       F1ZJ                                               
      F2XP(J)=       F2XJ                                               
      F2YP(J)=       F2YJ                                               
      F2ZP(J)=       F2ZJ                                               
      FPXP(J)=       FPXJ                                               
      FPYP(J)=       FPYJ                                               
      FPZP(J)=       FPZJ                                               
  237 CONTINUE                                                          
      DO 238 J=1,L                                                      
      K=IND(J)                                                          
      FOX(K)=FOX(K)-FOXP(J)                                             
      FOY(K)=FOY(K)-FOYP(J)                                             
      FOZ(K)=FOZ(K)-FOZP(J)                                             
      F1X(K)=F1X(K)-F1XP(J)                                             
      F1Y(K)=F1Y(K)-F1YP(J)                                             
      F1Z(K)=F1Z(K)-F1ZP(J)                                             
      F2X(K)=F2X(K)-F2XP(J)                                             
      F2Y(K)=F2Y(K)-F2YP(J)                                             
      F2Z(K)=F2Z(K)-F2ZP(J)                                             
      FPX(K)=FPX(K)-FPXP(J)                                             
      FPY(K)=FPY(K)-FPYP(J)                                             
      FPZ(K)=FPZ(K)-FPZP(J)                                             
 238  CONTINUE                                                          
      FOX(I)=FOX(I)+FXO                                                 
      FOY(I)=FOY(I)+FYO                                                 
      FOZ(I)=FOZ(I)+FZO                                                 
      F1X(I)=F1X(I)+FX1                                                 
      F1Y(I)=F1Y(I)+FY1                                                 
      F1Z(I)=F1Z(I)+FZ1                                                 
      F2X(I)=F2X(I)+FX2                                                 
      F2Y(I)=F2Y(I)+FY2                                                 
      F2Z(I)=F2Z(I)+FZ2                                                 
      FPX(I)=FPX(I)+FXP                                                 
      FPY(I)=FPY(I)+FYP                                                 
      FPZ(I)=FPZ(I)+FZP                                                 
  240 CONTINUE                                                          
C*************************************************************          
C*****        WATER - ION INTERACTIONS                   *****          
C*************************************************************          
CX20766 CALL DOOUTX( 766 )                                              
CY20766 CALL DOOUTY( 766 )                                              
      ISIT=NSITES(1)                                                    
      NSP=NSPECI(1)                                                     
      IF(NION.LE.0) GO TO 325                                           
      ION=0                                                             
  250  ION=ION+1                                                        
C     DO 250 ION=1,NION                                                 
      K=NSP+ION                                                         
      JSP=NSS+ION                                                       
      XX=0.D0                                                           
      YY=0.D0                                                           
      ZZ=0.D0                                                           
CX10767 CALL DOINX( 767 )                                               
CY10767 CALL DOINY( 767 )                                               
      DO 320 I=1,NSP                                                    
      INS=(I-1)*ISIT                                                    
      XD=X0(I)-X0(K)                                                    
      YD=Y0(I)-Y0(K)                                                    
      ZD=Z0(I)-Z0(K)                                                    
      XD=XD-2.D0*ALENGT*DBLE(INT(XD*ALENGM))                            
      YD=YD-2.D0*ALENGT*DBLE(INT(YD*ALENGM))                            
      ZD=ZD-2.D0*DBLE(INT(ZD))                                          
      DX=XD+SX(INS+1)                                                   
      DY=YD+SY(INS+1)                                                   
      DZ=ZD+SZ(INS+1)                                                   
      DD=DX*DX+DY*DY+DZ*DZ                                              
      DDI=1.D0/DD                                                       
      DSQI=SQRT(DDI)                                                    
      D6I=DDI*DDI*DDI                                                   
      D12I=D6I*D6I                                                      
      T1=A(1)*DSQI                                                      
      T2=B(1)*D6I                                                       
      T3=C(1)*D12I                                                      
      EWIT=  T1+T2+T3                                                   
      GGG=T1+6.D0*T2+12.D0*T3                                           
      GGG=GGG*DDI                                                       
      XO      =         DX*GGG                                          
      YO      =         DY*GGG                                          
      ZO      =         DZ*GGG                                          
      DX=XD+SX(INS+2)                                                   
      DY=YD+SY(INS+2)                                                   
      DZ=ZD+SZ(INS+2)                                                   
      DD=DX*DX+DY*DY+DZ*DZ                                              
      DDI=1.D0/DD                                                       
      DSQI=SQRT(DDI)                                                    
      D6I=DDI*DDI*DDI                                                   
      D12I=D6I*D6I                                                      
      T1=A(2)*DSQI                                                      
      T2=B(2)*D6I                                                       
      T3=C(2)*D12I                                                      
      EWIT=EWIT+T1+T2+T3                                                
      GGG=T1+6.D0*T2+12.D0*T3                                           
      GGG=GGG*DDI                                                       
      X1      =         DX*GGG                                          
      Y1      =         DY*GGG                                          
      Z1      =         DZ*GGG                                          
      DX=XD+SX(INS+3)                                                   
      DY=YD+SY(INS+3)                                                   
      DZ=ZD+SZ(INS+3)                                                   
      DD=DX*DX+DY*DY+DZ*DZ                                              
      DDI=1.D0/DD                                                       
      DSQI=SQRT(DDI)                                                    
      D6I=DDI*DDI*DDI                                                   
      D12I=D6I*D6I                                                      
      T1=A(3)*DSQI                                                      
      T2=B(3)*D6I                                                       
      T3=C(3)*D12I                                                      
      EWIT=EWIT+T1+T2+T3                                                
      EWI=EWI+EWIT                                                      
      GGG=T1+6.D0*T2+12.D0*T3                                           
      GGG=GGG*DDI                                                       
      X2      =         DX*GGG                                          
      Y2      =         DY*GGG                                          
      Z2      =         DZ*GGG                                          
      FOX(I)=FOX(I)+XO                                                  
      FOY(I)=FOY(I)+YO                                                  
      FOZ(I)=FOZ(I)+ZO                                                  
      F1X(I)=F1X(I)+X1                                                  
      F1Y(I)=F1Y(I)+Y1                                                  
      F1Z(I)=F1Z(I)+Z1                                                  
      F2X(I)=F2X(I)+X2                                                  
      F2Y(I)=F2Y(I)+Y2                                                  
      F2Z(I)=F2Z(I)+Z2                                                  
      XX=XX+XO+X1+X2                                                    
      YY=YY+YO+Y1+Y2                                                    
      ZZ=ZZ+ZO+Z1+Z2                                                    
  320 CONTINUE                                                          
CX20768 CALL DOOUTX( 768 )                                              
CY20768 CALL DOOUTY( 768 )                                              
      FSX(JSP)=FSX(JSP)-XX                                              
      FSY(JSP)=FSY(JSP)-YY                                              
      FSZ(JSP)=FSZ(JSP)-ZZ                                              
      IF(ION.LT.NION) GO TO 250                                         
C 250 CONTINUE                                                          
  325 CONTINUE                                                          
C*****************************************************************      
C*****        WATER - Z-DNA INTERACTIONS                     *****      
C*****************************************************************      
CX10769 CALL DOINX( 769 )                                               
CY10769 CALL DOINY( 769 )                                               
      DO 112 IA=1,NATOMS                                                
      FAX(IA)=0.D0                                                      
      FAY(IA)=0.D0                                                      
      FAZ(IA)=0.D0                                                      
 112  CONTINUE                                                          
CX20770 CALL DOOUTX( 770 )                                              
CY20770 CALL DOOUTY( 770 )                                              
      ISIT=NSITES(1)                                                    
      NSP=NSPECI(1)                                                     
CX10771 CALL DOINX( 771 )                                               
CY10771 CALL DOINY( 771 )                                               
      DO 500 I=1,NSP                                                    
      INS=(I-1)*ISIT                                                    
      DO 350 IA=1,NATOMS                                                
      K=IA+NOP                                                          
C     JSP=NSS+NION+IA                                                   
      XD=X0(I)-X0(K)                                                    
      YD=Y0(I)-Y0(K)                                                    
      ZD=Z0(I)-Z0(K)                                                    
      XD=XD-2.D0*ALENGT*DBLE(INT(XD*ALENGM))                            
      YD=YD-2.D0*ALENGT*DBLE(INT(YD*ALENGM))                            
      ZD=ZD-2.D0*DBLE(INT(ZD))                                          
      DX=XD+SX(INS+1)                                                   
      DY=YD+SY(INS+1)                                                   
      DZ=ZD+SZ(INS+1)                                                   
      DD=DX*DX+DY*DY+DZ*DZ                                              
      DDI=1.D0/DD                                                       
      DSQI=SQRT(DDI)                                                    
      D6I=DDI*DDI*DDI                                                   
      D12I=D6I*D6I                                                      
      T1=AA(1,IA)*DSQI                                                  
      T2=BB(1,IA)*D6I                                                   
      T3=CC(1,IA)*D12I                                                  
      EWAT=    T1+T2+T3                                                 
      GGG=T1+6.D0*T2+12.D0*T3                                           
      GGG=GGG*DDI                                                       
      FOX(I)=FOX(I)+(DX*GGG)                                            
      FOY(I)=FOY(I)+(DY*GGG)                                            
      FOZ(I)=FOZ(I)+(DZ*GGG)                                            
      FAX(IA)=FAX(IA)+(DX*GGG)                                          
      FAY(IA)=FAY(IA)+(DY*GGG)                                          
      FAZ(IA)=FAZ(IA)+(DZ*GGG)                                          
      DX=XD+SX(INS+2)                                                   
      DY=YD+SY(INS+2)                                                   
      DZ=ZD+SZ(INS+2)                                                   
      DD=DX*DX+DY*DY+DZ*DZ                                              
      DDI=1.D0/DD                                                       
      DSQI=SQRT(DDI)                                                    
      D6I=DDI*DDI*DDI                                                   
      D12I=D6I*D6I                                                      
      T1=AA(2,IA)*DSQI                                                  
      T2=BB(2,IA)*D6I                                                   
      T3=CC(2,IA)*D12I                                                  
      EWAT=EWAT+T1+T2+T3                                                
      GGG=T1+6.D0*T2+12.D0*T3                                           
      GGG=GGG*DDI                                                       
      F1X(I)=F1X(I)+(DX*GGG)                                            
      F1Y(I)=F1Y(I)+(DY*GGG)                                            
      F1Z(I)=F1Z(I)+(DZ*GGG)                                            
      FAX(IA)=FAX(IA)+(DX*GGG)                                          
      FAY(IA)=FAY(IA)+(DY*GGG)                                          
      FAZ(IA)=FAZ(IA)+(DZ*GGG)                                          
      DX=XD+SX(INS+3)                                                   
      DY=YD+SY(INS+3)                                                   
      DZ=ZD+SZ(INS+3)                                                   
      DD=DX*DX+DY*DY+DZ*DZ                                              
      DDI=1.D0/DD                                                       
      DSQI=SQRT(DDI)                                                    
      D6I=DDI*DDI*DDI                                                   
      D12I=D6I*D6I                                                      
      T1=AA(3,IA)*DSQI                                                  
      T2=BB(3,IA)*D6I                                                   
      T3=CC(3,IA)*D12I                                                  
      EWAT=EWAT+T1+T2+T3                                                
      GGG=T1+6.D0*T2+12.D0*T3                                           
      GGG=GGG*DDI                                                       
      F2X(I)=F2X(I)+(DX*GGG)                                            
      F2Y(I)=F2Y(I)+(DY*GGG)                                            
      F2Z(I)=F2Z(I)+(DZ*GGG)                                            
      FAX(IA)=FAX(IA)+(DX*GGG)                                          
      FAY(IA)=FAY(IA)+(DY*GGG)                                          
      FAZ(IA)=FAZ(IA)+(DZ*GGG)                                          
      EWA=EWA+EWAT                                                      
  350 CONTINUE                                                          
  500 CONTINUE                                                          
CX20772 CALL DOOUTX( 772 )                                              
CY20772 CALL DOOUTY( 772 )                                              
CX10773 CALL DOINX( 773 )                                               
CY10773 CALL DOINY( 773 )                                               
      DO 110 I=1,NSP                                                    
      INS=(I-1)*ISIT                                                    
      FSX(INS+1)=FSX(INS+1)+FOX(I)                                      
      FSY(INS+1)=FSY(INS+1)+FOY(I)                                      
      FSZ(INS+1)=FSZ(INS+1)+FOZ(I)                                      
      FSX(INS+2)=FSX(INS+2)+F1X(I)                                      
      FSY(INS+2)=FSY(INS+2)+F1Y(I)                                      
      FSZ(INS+2)=FSZ(INS+2)+F1Z(I)                                      
      FSX(INS+3)=FSX(INS+3)+F2X(I)                                      
      FSY(INS+3)=FSY(INS+3)+F2Y(I)                                      
      FSZ(INS+3)=FSZ(INS+3)+F2Z(I)                                      
      FSX(INS+4)=FSX(INS+4)+FPX(I)                                      
      FSY(INS+4)=FSY(INS+4)+FPY(I)                                      
      FSZ(INS+4)=FSZ(INS+4)+FPZ(I)                                      
 110  CONTINUE                                                          
C*****************************************************************      
C*****************************************************************      
C*****      ION - Z-DNA INTERACTIONS                         *****      
C*****************************************************************      
CX20774 CALL DOOUTX( 774 )                                              
CY20774 CALL DOOUTY( 774 )                                              
      IF(NION.LE.0) GO TO 999                                           
CX10775 CALL DOINX( 775 )                                               
CY10775 CALL DOINY( 775 )                                               
      DO 700 I=1,NION                                                   
      K=I+NSP                                                           
      ISP=NSS+I                                                         
      XX=0.D0                                                           
      YY=0.D0                                                           
      ZZ=0.D0                                                           
      XD=X0(K)                                                          
      YD=Y0(K)                                                          
      ZD=Z0(K)                                                          
      DO 650 J=1,NATOMS                                                 
      L=NOP+J                                                           
C     JSP=NSS+NION+J                                                    
      DX=XD-X0(L)                                                       
      DY=YD-Y0(L)                                                       
      DZ=ZD-Z0(L)                                                       
      DX=DX-2.D0*ALENGT*DBLE(INT(DX*ALENGM))                            
      DY=DY-2.D0*ALENGT*DBLE(INT(DY*ALENGM))                            
      DZ=DZ-2.D0*DBLE(INT(DZ))                                          
      DD=DX*DX+DY*DY+DZ*DZ                                              
      DDI=1.D0/DD                                                       
      DSQI=SQRT(DDI)                                                    
      D6I=DDI*DDI*DDI                                                   
      D12I=D6I*D6I                                                      
      T1=AAA(J)*DSQI                                                    
      T2=BBB(J)*D6I                                                     
      T3=CCC(J)*D12I                                                    
      EIA=EIA+T1+T2+T3                                                  
      GGG=T1+6.D0*T2+12.D0*T3                                           
      GGG=GGG*DDI                                                       
      XX      =XX      +(GGG*DX)                                        
      YY      =YY      +(GGG*DY)                                        
      ZZ      =ZZ      +(GGG*DZ)                                        
      FAX(J)=FAX(J)+(GGG*DX)                                            
      FAY(J)=FAY(J)+(GGG*DY)                                            
      FAZ(J)=FAZ(J)+(GGG*DZ)                                            
  650 CONTINUE                                                          
      FSX(ISP)=FSX(ISP)+XX                                              
      FSY(ISP)=FSY(ISP)+YY                                              
      FSZ(ISP)=FSZ(ISP)+ZZ                                              
  700 CONTINUE                                                          
CX20776 CALL DOOUTX( 776 )                                              
CY20776 CALL DOOUTY( 776 )                                              
CX10777 CALL DOINX( 777 )                                               
CY10777 CALL DOINY( 777 )                                               
      DO 111 IA=1,NATOMS                                                
      JSP=NSS+NION+IA                                                   
      FSX(JSP)=FSX(JSP)-FAX(IA)                                         
      FSY(JSP)=FSY(JSP)-FAY(IA)                                         
      FSZ(JSP)=FSZ(JSP)-FAZ(IA)                                         
 111  CONTINUE                                                          
C************************************************************           
CX20778 CALL DOOUTX( 778 )                                              
CY20778 CALL DOOUTY( 778 )                                              
C                                                           *           
C         ION  -  ION INTERACTIONS                          *           
C************************************************************           
      IF(NION.EQ.1) GO TO 999                                           
CX10779 CALL DOINX( 779 )                                               
CY10779 CALL DOINY( 779 )                                               
      DO 721 I=1,NION                                                   
      FIX(I)=0.D0                                                       
      FIY(I)=0.D0                                                       
      FIZ(I)=0.D0                                                       
 721  CONTINUE                                                          
CX20780 CALL DOOUTX( 780 )                                              
CY20780 CALL DOOUTY( 780 )                                              
CX10781 CALL DOINX( 781 )                                               
CY10781 CALL DOINY( 781 )                                               
      DO 720 I=2,NION                                                   
      XX=0.D0                                                           
      YY=0.D0                                                           
      ZZ=0.D0                                                           
      K=I+NSP                                                           
      XD=X0(K)                                                          
      YD=Y0(K)                                                          
      ZD=Z0(K)                                                          
      DO 730 J=1,I-1                                                    
      L=NSP+J                                                           
      JSP=NSS+J                                                         
      DX=XD-X0(L)                                                       
      DY=YD-Y0(L)                                                       
      DZ=ZD-Z0(L)                                                       
      DX=DX-2.D0*ALENGT*DBLE(INT(DX*ALENGM))                            
      DY=DY-2.D0*ALENGT*DBLE(INT(DY*ALENGM))                            
      DZ=DZ-2.D0*DBLE(INT(DZ))                                          
      DD=DX*DX+DY*DY+DZ*DZ                                              
      DDI=1.D0/DD                                                       
      DSQI=SQRT(DDI)                                                    
      D6I=DDI*DDI*DDI                                                   
      D12I=D6I*D6I                                                      
      T1=AAAA*DSQI                                                      
      T2=BBBB*D6I                                                       
      T3=CCCC*D12I                                                      
      EII=EII+T1+T2+T3                                                  
      GGG=T1+6.D0*T2+12.D0*T3                                           
      GGG=GGG*DDI                                                       
      XX      =XX      +(GGG*DX)                                        
      YY      =YY      +(GGG*DY)                                        
      ZZ      =ZZ      +(GGG*DZ)                                        
      FSX(JSP)=FSX(JSP)-(GGG*DX)                                        
      FSY(JSP)=FSY(JSP)-(GGG*DY)                                        
      FSZ(JSP)=FSZ(JSP)-(GGG*DZ)                                        
  730 CONTINUE                                                          
      FIX(I)=FIX(I)+XX                                                  
      FIY(I)=FIY(I)+YY                                                  
      FIZ(I)=FIZ(I)+ZZ                                                  
      GO TO 720                                                         
  720 CONTINUE                                                          
CX20782 CALL DOOUTX( 782 )                                              
CY20782 CALL DOOUTY( 782 )                                              
CX10783 CALL DOINX( 783 )                                               
CY10783 CALL DOINY( 783 )                                               
      DO 722 I=1,NION                                                   
      ISP=NSS+I                                                         
      FSX(ISP)=FSX(ISP)+FIX(I)                                          
      FSY(ISP)=FSY(ISP)+FIY(I)                                          
      FSZ(ISP)=FSZ(ISP)+FIZ(I)                                          
  722 CONTINUE                                                          
CX20784 CALL DOOUTX( 784 )                                              
CY20784 CALL DOOUTY( 784 )                                              
  999 CONTINUE                                                          
C*****************************************************************      
CX10785 CALL DOINX( 785 )                                               
CY10785 CALL DOINY( 785 )                                               
C     CALL ZERO(NOP,FX(1))                                              
C     CALL ZERO(NOP,FY(1))                                              
C     CALL ZERO(NOP,FZ(1))                                              
      DO 11 I=1,NOP                                                     
      FX(I)=0.D0                                                        
      FY(I)=0.D0                                                        
      FZ(I)=0.D0                                                        
  11  CONTINUE                                                          
CX20786 CALL DOOUTX( 786 )                                              
CY20786 CALL DOOUTY( 786 )                                              
CX10787 CALL DOINX( 787 )                                               
CY10787 CALL DOINY( 787 )                                               
      DO 800 I=1,NSPECI(1)                                              
      ISP=(I-1)*NSITES(1)                                               
      FX(I)=FX(I)+FSX(ISP+1)+FSX(ISP+2)+FSX(ISP+3)+FSX(ISP+4)           
      FY(I)=FY(I)+FSY(ISP+1)+FSY(ISP+2)+FSY(ISP+3)+FSY(ISP+4)           
      FZ(I)=FZ(I)+FSZ(ISP+1)+FSZ(ISP+2)+FSZ(ISP+3)+FSZ(ISP+4)           
      TX(I)=0.D0                                                        
      TY(I)=0.D0                                                        
      TZ(I)=0.D0                                                        
      TX(I)=      FSZ(ISP+1)*SY(ISP+1)-FSY(ISP+1)*SZ(ISP+1)             
     1     +      FSZ(ISP+2)*SY(ISP+2)-FSY(ISP+2)*SZ(ISP+2)             
     2     +      FSZ(ISP+3)*SY(ISP+3)-FSY(ISP+3)*SZ(ISP+3)             
     3     +      FSZ(ISP+4)*SY(ISP+4)-FSY(ISP+4)*SZ(ISP+4)             
      TY(I)=      FSX(ISP+1)*SZ(ISP+1)-FSZ(ISP+1)*SX(ISP+1)             
     1     +      FSX(ISP+2)*SZ(ISP+2)-FSZ(ISP+2)*SX(ISP+2)             
     2     +      FSX(ISP+3)*SZ(ISP+3)-FSZ(ISP+3)*SX(ISP+3)             
     3     +      FSX(ISP+4)*SZ(ISP+4)-FSZ(ISP+4)*SX(ISP+4)             
      TZ(I)=      FSY(ISP+1)*SX(ISP+1)-FSX(ISP+1)*SY(ISP+1)             
     1     +      FSY(ISP+2)*SX(ISP+2)-FSX(ISP+2)*SY(ISP+2)             
     2     +      FSY(ISP+3)*SX(ISP+3)-FSX(ISP+3)*SY(ISP+3)             
     3     +      FSY(ISP+4)*SX(ISP+4)-FSX(ISP+4)*SY(ISP+4)             
  800 CONTINUE                                                          
CX20788 CALL DOOUTX( 788 )                                              
CY20788 CALL DOOUTY( 788 )                                              
      N=NSS                                                             
CX10789 CALL DOINX( 789 )                                               
CY10789 CALL DOINY( 789 )                                               
      DO 810 I=NSPECI(1)+1,NOP                                          
      N=N+1                                                             
      FX(I)=FSX(N)                                                      
      FY(I)=FSY(N)                                                      
      FZ(I)=FSZ(N)                                                      
  810 CONTINUE                                                          
CX20790 CALL DOOUTX( 790 )                                              
CY20790 CALL DOOUTY( 790 )                                              
      CALL FCHECK(FSX,FSY,FSZ)                                          
CX40078 CALL SBOUTX (78)                                                
CY40078 CALL SBOUTY (78)                                                
      RETURN                                                            
      END                                                               
C*******************************************************************    
C****************************************************************       
C*****  THIS SUBROUTINE SCALES THE ANGULAR VELOCITIES IN THE ****       
C*****  SECOND ORDER QUARTERNION METHOD. THE QUARTERNION     ****       
C*****  VELOCITIES ARE COMPUTED FROM 2*DQ=Q*W                ****       
C****************************************************************       
      SUBROUTINE SCALER(QA0,QB0,QC0,QD0,QA,QB,QC,QD,AX,AY,AZ)           
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION QA0(*),QB0(*),QC0(*),QD0(*),AX(*),AY(*),AZ(*)           
      DIMENSION QA(*),QB(*),QC(*),QD(*)                                 
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
      COMMON/CONSTS/ AVSNO,BOLTZ,PLANCK,PI,ELECHG,EPS0                  
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      COMMON /AVERIZ/ V(13),VV(13)                                      
CX30079 CALL SBINX (79)                                                 
CY30079 CALL SBINY (79)                                                 
      IF(NSTEP.GE.NSSTS) WRITE(66,10) NSTEP                             
   10 FORMAT(1X,'*** TEMPERATURE SCALING SWITHED OFF AT TIME STEP ',I6) 
      IF(NSTEP.GE.NSSTS) LTS=.FALSE.                                    
      IF(MOD(NSTEP,NSBTS).NE.0) THEN                                    
CX40080 CALL SBOUTX (80)                                                
CY40080 CALL SBOUTY (80)                                                
      RETURN                                                            
      ENDIF                                                             
      K=0                                                               
CX10791 CALL DOINX( 791 )                                               
CY10791 CALL DOINY( 791 )                                               
      DO 30 N=1,NTYPES                                                  
      NSP=NSPECI(N)                                                     
      NSIT=NSITES(N)                                                    
      IF(NSIT.EQ.1) GO TO 30                                            
      DO 60 L=1,NSP                                                     
      K=K+1                                                             
      QA(K)=0.5D0*(-QC0(K)*AX(K)-QD0(K)*AY(K)+QB0(K)*AZ(K))*TSTEP       
      QB(K)=0.5D0*(+QD0(K)*AX(K)-QC0(K)*AY(K)-QA0(K)*AZ(K))*TSTEP       
      QC(K)=0.5D0*(+QA0(K)*AX(K)+QB0(K)*AY(K)+QD0(K)*AZ(K))*TSTEP       
      QD(K)=0.5D0*(-QB0(K)*AX(K)+QA0(K)*AY(K)-QC0(K)*AZ(K))*TSTEP       
   60 CONTINUE                                                          
   30 CONTINUE                                                          
CX20792 CALL DOOUTX( 792 )                                              
CY20792 CALL DOOUTY( 792 )                                              
      NAVRS=0                                                           
CX10793 CALL DOINX( 793 )                                               
CY10793 CALL DOINY( 793 )                                               
      DO 40 I=1,13                                                      
      V(I)=0.D0                                                         
      VV(I)=0.D0                                                        
   40 CONTINUE                                                          
CX20794 CALL DOOUTX( 794 )                                              
CY20794 CALL DOOUTY( 794 )                                              
CX40080 CALL SBOUTX (80)                                                
CY40080 CALL SBOUTY (80)                                                
      RETURN                                                            
      END                                                               
C***************************************************************        
C*****  THIS PROGRAM SUBTRACTS OUT BULK  VELOCITY           ****        
C***************************************************************        
      SUBROUTINE RNMTRN(WX0,WY0,WZ0)                                    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION WX0(*),WY0(*),WZ0(*)                                    
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
      COMMON/CONSTS/ AVSNO,BOLTZ,PLANCK,PI,ELECHG,EPS0                  
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      COMMON /AVERIZ/ V(13),VV(13)                                      
CX30081 CALL SBINX (81)                                                 
CY30081 CALL SBINY (81)                                                 
      IF(NSTEP.GE.NSSTS) WRITE(66,10) NSTEP                             
   10 FORMAT(1X,'*** TEMPERATURE SCALING SWITHED OFF AT TIME STEP ',I6) 
      IF(NSTEP.GE.NSSTS) LTS=.FALSE.                                    
      IF(MOD(NSTEP,NSBTS).NE.0) THEN                                    
CX40082 CALL SBOUTX (82)                                                
CY40082 CALL SBOUTY (82)                                                
      RETURN                                                            
      ENDIF                                                             
      I=0                                                               
CX10795 CALL DOINX( 795 )                                               
CY10795 CALL DOINY( 795 )                                               
      DO 100 N=1,NTYPES                                                 
      NSP=NSPECI(N)                                                     
      FSP=DBLE(NSP)                                                     
C      NS=NSITES(N)                                                     
      BE=BOLTZ/UNITE                                                    
      TT=BE*TTRAN(N)                                                    
      IF(LLF) TERMT=DSUMM(N)                                            
      IF(LPC) TERMT=TSTEP**2/DSUMM(N)                                   
      W =2.D0*SQRT(TERMT*TT)                                            
      W1=SSUM(NSP,WX0(I+1),1)/FSP                                       
      W2=SSUM(NSP,WY0(I+1),1)/FSP                                       
      W3=SSUM(NSP,WZ0(I+1),1)/FSP                                       
      WS=0.D0                                                           
      DO 40 J=1,NSP                                                     
      I=I+1                                                             
      WX0(I)=WX0(I)-W1                                                  
      WY0(I)=WY0(I)-W2                                                  
      WZ0(I)=WZ0(I)-W3                                                  
      WS=WS+WX0(I)**2+WY0(I)**2+WZ0(I)**2                               
   40 CONTINUE                                                          
      IF(.NOT.LLF) GO TO 50                                             
      WS=0.5D0*W/SQRT(WS/(3.D0*FSP))                                    
      GO TO 60                                                          
   50 IF(.NOT.LPC) GO TO 60                                             
      WS=SQRT(3.D0*TSTEP*TSTEP*FSP*TT/DSUMM(N)/WS)                      
   60 CONTINUE                                                          
      I=I-NSP                                                           
      DO 70 J=1,NSP                                                     
      I=I+1                                                             
      WX0(I)=WS*WX0(I)                                                  
      WY0(I)=WS*WY0(I)                                                  
      WZ0(I)=WS*WZ0(I)                                                  
   70 CONTINUE                                                          
  100 CONTINUE                                                          
CX20796 CALL DOOUTX( 796 )                                              
CY20796 CALL DOOUTY( 796 )                                              
      NAVRS=0                                                           
CX10797 CALL DOINX( 797 )                                               
CY10797 CALL DOINY( 797 )                                               
      DO 45 I=1,13                                                      
      V(I)=0.D0                                                         
      VV(I)=0.D0                                                        
   45 CONTINUE                                                          
CX20798 CALL DOOUTX( 798 )                                              
CY20798 CALL DOOUTY( 798 )                                              
CX40082 CALL SBOUTX (82)                                                
CY40082 CALL SBOUTY (82)                                                
      RETURN                                                            
      END                                                               
C****************************************************************       
C**** THIS PROGRAM SUBTRACTS BULK MOTION                      ***       
C****************************************************************       
      SUBROUTINE RNMROT(AX0,AY0,AZ0)                                    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
      DIMENSION AX0(*),AY0(*),AZ0(*)                                    
      LOGICAL        MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/LCONTR/ MC,MD,LRES,LSAVE,LSPA,INMC,LTS,NOTR,LGR,LWW        
     X              ,LWC1,LWA1,LWCN,LWAN,LWC1A1,LQ2PC,LCA,LWWW,LLF,LPC  
      COMMON/CONTRL/NSTEPS,NSTEP,NSBPO,NSBTS,NSSTS,NSTCK,NTYPES,NOP     
     X              ,NSITES(3),NSPECI(3),NSITET,NENT,NATOMS,NION        
     X              ,NDIM,NSS,NAVRS,NOP3,NOP4,NSS3,NPOCS                
      COMMON/CONVRT/ DSUMM(3),DRI(3,3),DQ(30),DSITE(3,30)               
     X              ,UNITM,UNITL,UNITE,COULF,ENERF,TEMPF,TSTEP          
      COMMON/CONSTS/ AVSNO,BOLTZ,PLANCK,PI,ELECHG,EPS0                  
      COMMON/PHYSIC/ TTRAN(3),TROT(3),RHO,VOLM,DT,FNOP,BREAK            
     X              ,DM(3,3),QM(9,3),TE(3),RE(3),TTS(3),RTS(3)          
     X              ,R2(3),EWW,EWI,EWA,EII,EIA                          
      COMMON /AVERIZ/ V(13),VV(13)                                      
CX30083 CALL SBINX (83)                                                 
CY30083 CALL SBINY (83)                                                 
      IF(NSTEP.GE.NSSTS) WRITE(66,10) NSTEP                             
   10 FORMAT(1X,'*** TEMPERATURE SCALING SWITHED OFF AT TIME STEP ',I6) 
      IF(NSTEP.GE.NSSTS) LTS=.FALSE.                                    
      IF(MOD(NSTEP,NSBTS).NE.0) THEN                                    
CX40084 CALL SBOUTX (84)                                                
CY40084 CALL SBOUTY (84)                                                
      RETURN                                                            
      ENDIF                                                             
      I=0                                                               
CX10799 CALL DOINX( 799 )                                               
CY10799 CALL DOINY( 799 )                                               
      DO 200 N=1,NTYPES                                                 
      IS=NSITES(N)                                                      
      NSP=NSPECI(N)                                                     
      IF(IS.EQ.1) GO TO 200                                             
      NSP=NSPECI(N)                                                     
      FSP=DBLE(NSP)                                                     
      BE=BOLTZ/UNITE                                                    
      TR=BE*TROT (N)                                                    
      AX=2.D0*SQRT(DRI(1,N)*TR)                                         
      AY=2.D0*SQRT(DRI(2,N)*TR)                                         
      AZ=2.D0*SQRT(DRI(3,N)*TR)                                         
      IF(.NOT.LPC) GO TO 125                                            
      AX=AX/DRI(1,N)                                                    
      AY=AY/DRI(2,N)                                                    
      AZ=AZ/DRI(3,N)                                                    
  125 CONTINUE                                                          
      IF(NSP.EQ.1) GO TO 200                                            
      A1=SSUM(NSP,AX0(I+1),1)/FSP                                       
      A2=SSUM(NSP,AY0(I+1),1)/FSP                                       
      A3=SSUM(NSP,AZ0(I+1),1)/FSP                                       
      ASX=0.D0                                                          
      ASY=0.D0                                                          
      ASZ=0.D0                                                          
      DO 140 J=1,NSP                                                    
      I=I+1                                                             
      AX0(I)=AX0(I)-A1                                                  
      AY0(I)=AY0(I)-A2                                                  
      AZ0(I)=AZ0(I)-A3                                                  
      ASX=ASX+AX0(I)**2                                                 
      ASY=ASY+AY0(I)**2                                                 
      ASZ=ASZ+AZ0(I)**2                                                 
  140 CONTINUE                                                          
      IF(LLF) THEN                                                      
      ASX=0.5D0*AX/SQRT(ASX/FSP )                                       
      ASY=0.5D0*AY/SQRT(ASY/FSP )                                       
      ASZ=0.5D0*AZ/SQRT(ASZ/FSP )                                       
      ELSE                                                              
      ASX=SQRT(FSP*TR/DRI(1,N)/ASX)                                     
      ASY=SQRT(FSP*TR/DRI(2,N)/ASY)                                     
      ASZ=SQRT(FSP*TR/DRI(3,N)/ASZ)                                     
      ENDIF                                                             
      I=I-NSP                                                           
      DO 170 J=1,NSP                                                    
      I=I+1                                                             
      AX0(I)=ASX*AX0(I)                                                 
      AY0(I)=ASY*AY0(I)                                                 
      AZ0(I)=ASZ*AZ0(I)                                                 
  170 CONTINUE                                                          
  200 CONTINUE                                                          
CX20800 CALL DOOUTX( 800 )                                              
CY20800 CALL DOOUTY( 800 )                                              
      NAVRS=0                                                           
CX10801 CALL DOINX( 801 )                                               
CY10801 CALL DOINY( 801 )                                               
      DO 45 I=1,13                                                      
      V(I)=0.D0                                                         
      VV(I)=0.D0                                                        
   45 CONTINUE                                                          
CX20802 CALL DOOUTX( 802 )                                              
CY20802 CALL DOOUTY( 802 )                                              
CX40084 CALL SBOUTX (84)                                                
CY40084 CALL SBOUTY (84)                                                
      RETURN                                                            
      END                                                               
