C    *******************************************************************        
C    ** FICHE F.7                                                     **        
C    ** TWO SEPARATE PARTS: NONRIGID AND RIGID TRIATOMIC MOLECULES.   **        
C    *******************************************************************        
                                                                                
                                                                                
                                                                                
C    *******************************************************************        
C    ** FICHE F.7 - PART A                                            **        
C    ** CONSTRAINT DYNAMICS FOR A NONRIGID TRIATOMIC MOLECULE.        **        
C    *******************************************************************        
                                                                                
                                                                                
                                                                                
        SUBROUTINE MOVE ( DT, TOL, MAXIT, BOX, K, WC )                          
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, OX, OY, OZ, FX, FY, FZ                    
        COMMON / BLOCK2 / DSQ, M                                                
                                                                                
C    *******************************************************************        
C    ** UPDATES ATOMIC POSITIONS WITH BOND CONSTRAINTS APPLIED.       **        
C    **                                                               **        
C    ** THIS ROUTINE USES THE VERLET ALGORITHM AND APPLIES BOND       **        
C    ** LENGTH CONSTRAINTS TO TWO BOND LENGTHS ONLY, 1-2 AND 2-3.     **        
C    ** WE SOLVE A SYSTEM OF QUADRATIC EQUATIONS ITERATIVELY, USING   **        
C    ** MATRIX INVERSION TO SOLVE ASSOCIATED LINEAR EQUATIONS,        **        
C    ** AND ITERATING TO CONVERGENCE.                                 **        
C    **                                                               **        
C    ** REFERENCE:                                                    **        
C    **                                                               **        
C    ** RYCKAERT ET AL., J. COMP. PHYS, 23, 327, 1977.                **        
C    **                                                               **        
C    ** PRINCIPAL VARIABLES:                                          **        
C    **                                                               **        
C    ** INTEGER N                          NUMBER OF MOLECULES        **        
C    ** INTEGER NA                         ATOMS PER MOL. (3 HERE)    **        
C    ** REAL    DT                         TIMESTEP                   **        
C    ** INTEGER MAXIT                      MAX NUMBER OF ITERATIONS   **        
C    ** REAL    TOL                        PRESCRIBED TOLERANCE       **        
C    ** REAL    BOX                        BOX LENGTH                 **        
C    ** REAL    K                          KINETIC ENERGY             **        
C    ** REAL    WC                         CONSTRAINT VIRIAL          **        
C    ** REAL    RX(N,NA),RY(N,NA),RZ(N,NA) ATOMIC POSITIONS AT TIME T **        
C    ** REAL    OX(N,NA),OY(N,NA),OZ(N,NA) OLD POSITIONS AT TIME T-DT **        
C    ** REAL    FX(N,NA),FY(N,NA),FZ(N,NA) ATOMIC FORCES AT TIME T    **        
C    ** REAL    M(NA)                      ATOMIC MASSES.             **        
C    ** REAL    DSQ(NA)                    SQUARED BOND LENGTHS       **        
C    ** DSQ(A) IS THE SQUARED BOND LENGTH BETWEEN ATOMS A AND A+1.    **        
C    **                                                               **        
C    ** USAGE:                                                        **        
C    **                                                               **        
C    ** THE ROUTINE SHOULD BE CALLED AFTER COMPUTING FORCES.          **        
C    ** IT RETURNS THE NEW POSITIONS (TIME T+DT) IN RX,RY,RZ, AND THE **        
C    ** CURRENT (TIME T) POSITIONS IN OX,OY,OZ.                       **        
C    ** IT ALSO COMPUTES THE KINETIC ENERGY K AND THE CONSTRAINT      **        
C    ** CONTRIBUTION TO THE ATOMIC VIRIAL WC.                         **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        INTEGER     NA                                                          
        PARAMETER ( NA = 3 )                                                    
                                                                                
        INTEGER     MAXIT                                                       
        REAL        TOL, DT, K, WC, BOX                                         
        REAL        RX(N,NA), RY(N,NA), RZ(N,NA)                                
        REAL        OX(N,NA), OY(N,NA), OZ(N,NA)                                
        REAL        FX(N,NA), FY(N,NA), FZ(N,NA)                                
        REAL        M(NA), DSQ(NA)                                              
                                                                                
        INTEGER     I, IT                                                       
        REAL        RX1, RY1, RZ1, RX2, RY2, RZ2, RX3, RY3, RZ3                 
        REAL        PX1, PY1, PZ1, PX2, PY2, PZ2, PX3, PY3, PZ3                 
        REAL        RX12, RY12, RZ12, RX23, RY23, RZ23                          
        REAL        PX12, PY12, PZ12, PX23, PY23, PZ23                          
        REAL        VXIA, VYIA, VZIA, OM1, OM2, OM3, DTSQ                       
        REAL        R12SQ, R23SQ, R12R23                                        
        REAL        P12SQ, P23SQ                                                
        REAL        P12R12, P12R23, P23R12, P23R23                              
        REAL        L12, L23, L12NEW, L23NEW                                    
        REAL        MAT11, MAT12, MAT21, MAT22                                  
        REAL        INV11, INV12, INV21, INV22                                  
        REAL        CONST1, CONST2, QUAD1, QUAD2, VEC1, VEC2                    
        REAL        Q111, Q122, Q112, Q211, Q222, Q212                          
        REAL        DETERM, LTOL, BOXINV                                        
        LOGICAL     DONE                                                        
                                                                                
C    *******************************************************************        
                                                                                
        BOXINV = 1.0 / BOX                                                      
        DTSQ   = DT ** 2                                                        
        LTOL   = TOL * DTSQ                                                     
        R12SQ  = DSQ(1)                                                         
        R23SQ  = DSQ(2)                                                         
        OM1    = 1.0 / M(1)                                                     
        OM2    = 1.0 / M(2)                                                     
        OM3    = 1.0 / M(3)                                                     
                                                                                
        K  = 0.0                                                                
        WC = 0.0                                                                
                                                                                
C    ** LOOP OVER MOLECULES **                                                  
                                                                                
        DO 2000 I = 1, N                                                        
                                                                                
C       ** VERLET ALGORITHM **                                                  
                                                                                
           RX1 = RX(I,1)                                                        
           RY1 = RY(I,1)                                                        
           RZ1 = RZ(I,1)                                                        
           PX1 = 2.0 * RX1 - OX(I,1) + DTSQ * FX(I,1) * OM1                     
           PY1 = 2.0 * RY1 - OY(I,1) + DTSQ * FY(I,1) * OM1                     
           PZ1 = 2.0 * RZ1 - OZ(I,1) + DTSQ * FZ(I,1) * OM1                     
           RX2 = RX(I,2)                                                        
           RY2 = RY(I,2)                                                        
           RZ2 = RZ(I,2)                                                        
           PX2 = 2.0 * RX2 - OX(I,2) + DTSQ * FX(I,2) * OM2                     
           PY2 = 2.0 * RY2 - OY(I,2) + DTSQ * FY(I,2) * OM2                     
           PZ2 = 2.0 * RZ2 - OZ(I,2) + DTSQ * FZ(I,2) * OM2                     
           RX3 = RX(I,3)                                                        
           RY3 = RY(I,3)                                                        
           RZ3 = RZ(I,3)                                                        
           PX3 = 2.0 * RX3 - OX(I,3) + DTSQ * FX(I,3) * OM3                     
           PY3 = 2.0 * RY3 - OY(I,3) + DTSQ * FY(I,3) * OM3                     
           PZ3 = 2.0 * RZ3 - OZ(I,3) + DTSQ * FZ(I,3) * OM3                     
                                                                                
C       ** CALCULATE RELATIVE VECTORS **                                        
                                                                                
           RX12 = RX1 - RX2                                                     
           RX12 = RX12 - ANINT ( RX12 * BOXINV ) * BOX                          
           RY12 = RY1 - RY2                                                     
           RY12 = RY12 - ANINT ( RY12 * BOXINV ) * BOX                          
           RZ12 = RZ1 - RZ2                                                     
           RZ12 = RZ12 - ANINT ( RZ12 * BOXINV ) * BOX                          
           RX23 = RX2 - RX3                                                     
           RX23 = RX23 - ANINT ( RX23 * BOXINV ) * BOX                          
           RY23 = RY2 - RY3                                                     
           RY23 = RY23 - ANINT ( RY23 * BOXINV ) * BOX                          
           RZ23 = RZ2 - RZ3                                                     
           RZ23 = RZ23 - ANINT ( RZ23 * BOXINV ) * BOX                          
           PX12 = PX1 - PX2                                                     
           PX12 = PX12 - ANINT ( PX12 * BOXINV ) * BOX                          
           PY12 = PY1 - PY2                                                     
           PY12 = PY12 - ANINT ( PY12 * BOXINV ) * BOX                          
           PZ12 = PZ1 - PZ2                                                     
           PZ12 = PZ12 - ANINT ( PZ12 * BOXINV ) * BOX                          
           PX23 = PX2 - PX3                                                     
           PX23 = PX23 - ANINT ( PX23 * BOXINV ) * BOX                          
           PY23 = PY2 - PY3                                                     
           PY23 = PY23 - ANINT ( PY23 * BOXINV ) * BOX                          
           PZ23 = PZ2 - PZ3                                                     
           PZ23 = PZ23 - ANINT ( PZ23 * BOXINV ) * BOX                          
                                                                                
C       ** CALCULATE SCALAR PRODUCTS **                                         
                                                                                
           R12R23 = RX12 * RX23 + RY12 * RY23 + RZ12 * RZ23                     
           P12SQ  = PX12 ** 2 + PY12 ** 2 + PZ12 ** 2                           
           P23SQ  = PX23 ** 2 + PY23 ** 2 + PZ23 ** 2                           
           P12R12 = PX12 * RX12 + PY12 * RY12 + PZ12 * RZ12                     
           P12R23 = PX12 * RX23 + PY12 * RY23 + PZ12 * RZ23                     
           P23R12 = PX23 * RX12 + PY23 * RY12 + PZ23 * RZ12                     
           P23R23 = PX23 * RX23 + PY23 * RY23 + PZ23 * RZ23                     
           CONST1 = R12SQ - P12SQ                                               
           CONST2 = R23SQ - P23SQ                                               
                                                                                
C       ** EVALUATE MATRIX ELEMENTS AND QUADRATIC COEFFICIENTS **               
                                                                                
           MAT11 =   2.0 * P12R12 * ( OM1 + OM2 )                               
           MAT12 = - 2.0 * P12R23 *   OM2                                       
           MAT21 = - 2.0 * P23R12 *   OM2                                       
           MAT22 =   2.0 * P23R23 * ( OM2 + OM3 )                               
                                                                                
           Q111 = - R12SQ  * ( OM1 + OM2 ) ** 2                                 
           Q122 = - R23SQ  *   OM2 ** 2                                         
           Q112 = + 2.0 * R12R23 * ( OM1 + OM2 ) * OM2                          
           Q211 = - R12SQ  *   OM2 ** 2                                         
           Q222 = - R23SQ  * ( OM2 + OM3 ) ** 2                                 
           Q212 = + 2.0 * R12R23 * ( OM2 + OM3 ) * OM2                          
                                                                                
C       ** INVERT MATRIX **                                                     
                                                                                
           DETERM = 1.0 / ( MAT11 * MAT22 - MAT21 * MAT12 )                     
           INV11 =  MAT22 * DETERM                                              
           INV12 = -MAT12 * DETERM                                              
           INV21 = -MAT21 * DETERM                                              
           INV22 =  MAT11 * DETERM                                              
                                                                                
C       ** PREPARE FOR ITERATIVE LOOP **                                        
                                                                                
           L12  = 0.0                                                           
           L23  = 0.0                                                           
           DONE = .FALSE.                                                       
           IT   = 0                                                             
                                                                                
C       ** BEGIN ITERATIVE LOOP **                                              
                                                                                
1000       IF ( ( .NOT. DONE ) .AND. ( IT .LE. MAXIT ) ) THEN                   
                                                                                
              QUAD1 = Q111 * L12 ** 2                                           
     :              + Q122 * L23 ** 2                                           
     :              + Q112 * L12 * L23                                          
              QUAD2 = Q211 * L12 ** 2                                           
     :              + Q222 * L23 ** 2                                           
     :              + Q212 * L12 * L23                                          
                                                                                
              VEC1 = CONST1 + QUAD1                                             
              VEC2 = CONST2 + QUAD2                                             
                                                                                
              L12NEW = INV11 * VEC1 + INV12 * VEC2                              
              L23NEW = INV21 * VEC1 + INV22 * VEC2                              
                                                                                
              DONE = ( ( ABS ( L12NEW - L12 ) .LE. LTOL ) .AND.                 
     :                 ( ABS ( L23NEW - L23 ) .LE. LTOL ) )                     
                                                                                
              L12 = L12NEW                                                      
              L23 = L23NEW                                                      
              IT  = IT + 1                                                      
                                                                                
              GOTO 1000                                                         
                                                                                
           ENDIF                                                                
                                                                                
C       ** END OF ITERATION **                                                  
                                                                                
           PX1 = PX1 + OM1 * ( L12 * RX12              )                        
           PY1 = PY1 + OM1 * ( L12 * RY12              )                        
           PZ1 = PZ1 + OM1 * ( L12 * RZ12              )                        
           PX2 = PX2 + OM2 * ( L23 * RX23 - L12 * RX12 )                        
           PY2 = PY2 + OM2 * ( L23 * RY23 - L12 * RY12 )                        
           PZ2 = PZ2 + OM2 * ( L23 * RZ23 - L12 * RZ12 )                        
           PX3 = PX3 + OM3 * (            - L23 * RX23 )                        
           PY3 = PY3 + OM3 * (            - L23 * RY23 )                        
           PZ3 = PZ3 + OM3 * (            - L23 * RZ23 )                        
                                                                                
           IF ( .NOT. DONE ) THEN                                               
                                                                                
              WRITE(*,'('' TOO MANY CONSTRAINT ITERATIONS '')')                 
              WRITE(*,'('' MOLECULE '',I5)') I                                  
              STOP                                                              
                                                                                
           ENDIF                                                                
                                                                                
C       ** CALCULATE KINETIC ENERGY CONTRIBUTION **                             
                                                                                
           VXIA = 0.5 * ( PX1 - OX(I,1) ) / DT                                  
           VYIA = 0.5 * ( PY1 - OY(I,1) ) / DT                                  
           VZIA = 0.5 * ( PZ1 - OZ(I,1) ) / DT                                  
           K = K + ( VXIA ** 2 + VYIA ** 2 + VZIA ** 2 ) * M(1)                 
           VXIA = 0.5 * ( PX2 - OX(I,2) ) / DT                                  
           VYIA = 0.5 * ( PY2 - OY(I,2) ) / DT                                  
           VZIA = 0.5 * ( PZ2 - OZ(I,2) ) / DT                                  
           K = K + ( VXIA ** 2 + VYIA ** 2 + VZIA ** 2 ) * M(2)                 
           VXIA = 0.5 * ( PX3 - OX(I,3) ) / DT                                  
           VYIA = 0.5 * ( PY3 - OY(I,3) ) / DT                                  
           VZIA = 0.5 * ( PZ3 - OZ(I,3) ) / DT                                  
           K = K + ( VXIA ** 2 + VYIA ** 2 + VZIA ** 2 ) * M(3)                 
                                                                                
C       ** CALCULATE CONSTRAINT VIRIAL CONTRIBUTION **                          
                                                                                
           WC = WC + L12 * R12SQ + L23 * R23SQ                                  
                                                                                
C       ** STORE AWAY POSITIONS **                                              
                                                                                
           OX(I,1) = RX1                                                        
           OY(I,1) = RY1                                                        
           OZ(I,1) = RZ1                                                        
           OX(I,2) = RX2                                                        
           OY(I,2) = RY2                                                        
           OZ(I,2) = RZ2                                                        
           OX(I,3) = RX3                                                        
           OY(I,3) = RY3                                                        
           OZ(I,3) = RZ3                                                        
           RX(I,1) = PX1                                                        
           RY(I,1) = PY1                                                        
           RZ(I,1) = PZ1                                                        
           RX(I,2) = PX2                                                        
           RY(I,2) = PY2                                                        
           RZ(I,2) = PZ2                                                        
           RX(I,3) = PX3                                                        
           RY(I,3) = PY3                                                        
           RZ(I,3) = PZ3                                                        
                                                                                
2000    CONTINUE                                                                
                                                                                
C    ** END OF LOOP OVER MOLECULES **                                           
                                                                                
        K  = 0.5 * K                                                            
        WC = WC / DTSQ / 3.0                                                    
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
C    *******************************************************************        
C    ** FICHE F.7 - PART B                                            **        
C    ** CONSTRAINT DYNAMICS FOR A RIGID NONLINEAR TRIATOMIC MOLECULE. **        
C    *******************************************************************        
                                                                                
                                                                                
                                                                                
        SUBROUTINE MOVE ( DT, TOL, MAXIT, BOX, K, WC )                          
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, OX, OY, OZ, FX, FY, FZ                    
        COMMON / BLOCK2 / DSQ, M                                                
                                                                                
C    *******************************************************************        
C    ** UPDATES ATOMIC POSITIONS WITH BOND CONSTRAINTS APPLIED.       **        
C    **                                                               **        
C    ** THIS ROUTINE USES THE VERLET ALGORITHM AND APPLIES BOND       **        
C    ** LENGTH CONSTRAINTS TO ALL THREE BONDS, NAMELY 1-2, 2-3, 3-1.  **        
C    ** WE SOLVE A SYSTEM OF QUADRATIC EQUATIONS ITERATIVELY, USING   **        
C    ** MATRIX INVERSION TO SOLVE ASSOCIATED LINEAR EQUATIONS,        **        
C    ** AND ITERATING TO CONVERGENCE.                                 **        
C    **                                                               **        
C    ** REFERENCE:                                                    **        
C    **                                                               **        
C    ** RYCKAERT ET AL., J. COMP. PHYS, 23, 327, 1977.                **        
C    **                                                               **        
C    ** PRINCIPAL VARIABLES:                                          **        
C    **                                                               **        
C    ** INTEGER N                          NUMBER OF MOLECULES        **        
C    ** INTEGER NA                         ATOMS PER MOL. (3 HERE)    **        
C    ** REAL    DT                         TIMESTEP                   **        
C    ** INTEGER MAXIT                      MAX NUMBER OF ITERATIONS   **        
C    ** REAL    BOX                        BOX LENGTH                 **        
C    ** REAL    TOL                        PRESCRIBED TOLERANCE       **        
C    ** REAL    K                          KINETIC ENERGY             **        
C    ** REAL    WC                         CONSTRAINT VIRIAL          **        
C    ** REAL    RX(N,NA),RY(N,NA),RZ(N,NA) ATOMIC POSITIONS AT TIME T **        
C    ** REAL    OX(N,NA),OY(N,NA),OZ(N,NA) OLD POSITIONS AT TIME T-DT **        
C    ** REAL    FX(N,NA),FY(N,NA),FZ(N,NA) ATOMIC FORCES AT TIME T    **        
C    ** REAL    M(NA)                      ATOMIC MASSES.             **        
C    ** REAL    DSQ(NA)                    SQUARED BOND LENGTHS       **        
C    ** DSQ(A) IS THE SQUARED BOND LENGTH BETWEEN ATOMS A AND A+1.    **        
C    **                                                               **        
C    ** USAGE:                                                        **        
C    **                                                               **        
C    ** THE ROUTINE SHOULD BE CALLED AFTER COMPUTING FORCES.          **        
C    ** IT RETURNS THE NEW POSITIONS (TIME T+DT) IN RX,RY,RZ, AND THE **        
C    ** CURRENT (TIME T) POSITIONS IN OX,OY,OZ.                       **        
C    ** IT ALSO COMPUTES THE KINETIC ENERGY K AND THE CONSTRAINT      **        
C    ** CONTRIBUTION TO THE VIRIAL WC.                                **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        INTEGER     NA                                                          
        PARAMETER ( NA = 3 )                                                    
                                                                                
        INTEGER     MAXIT                                                       
        REAL        TOL, DT, K, WC, BOX                                         
        REAL        RX(N,NA), RY(N,NA), RZ(N,NA)                                
        REAL        OX(N,NA), OY(N,NA), OZ(N,NA)                                
        REAL        FX(N,NA), FY(N,NA), FZ(N,NA)                                
        REAL        M(NA), DSQ(NA)                                              
                                                                                
        INTEGER     IT, I                                                       
        REAL        RX1, RY1, RZ1, RX2, RY2, RZ2, RX3, RY3, RZ3                 
        REAL        PX1, PY1, PZ1, PX2, PY2, PZ2, PX3, PY3, PZ3                 
        REAL        RX12, RY12, RZ12, RX23, RY23, RZ23, RX31, RY31, RZ31        
        REAL        PX12, PY12, PZ12, PX23, PY23, PZ23, PX31, PY31, PZ31        
        REAL        VXIA, VYIA, VZIA, OM1, OM2, OM3, DTSQ                       
        REAL        R12SQ, R23SQ, R31SQ, R12R23, R23R31, R31R12                 
        REAL        P12SQ, P23SQ, P31SQ                                         
        REAL        P12R12, P12R23, P12R31                                      
        REAL        P23R12, P23R23, P23R31                                      
        REAL        P31R12, P31R23, P31R31                                      
        REAL        L12, L23, L31, L12NEW, L23NEW, L31NEW, LTOL, BOXINV         
        REAL        CONST1, CONST2, CONST3                                      
        REAL        QUAD1, QUAD2, QUAD3, VEC1, VEC2, VEC3                       
        REAL        Q111, Q122, Q133, Q112, Q123, Q131                          
        REAL        Q211, Q222, Q233, Q212, Q223, Q231                          
        REAL        Q311, Q322, Q333, Q312, Q323, Q331                          
        REAL        MAT(3,3), INV(3,3)                                          
        LOGICAL     DONE                                                        
                                                                                
C    *******************************************************************        
                                                                                
        BOXINV = 1.0 / BOX                                                      
        DTSQ   = DT ** 2                                                        
        LTOL   = TOL * DTSQ                                                     
        R12SQ  = DSQ(1)                                                         
        R23SQ  = DSQ(2)                                                         
        R31SQ  = DSQ(3)                                                         
        OM1    = 1.0 / M(1)                                                     
        OM2    = 1.0 / M(2)                                                     
        OM3    = 1.0 / M(3)                                                     
                                                                                
        K  = 0.0                                                                
        WC = 0.0                                                                
                                                                                
C    ** LOOP OVER MOLECULES **                                                  
                                                                                
        DO 2000 I = 1, N                                                        
                                                                                
C       ** VERLET ALGORITHM **                                                  
                                                                                
           RX1 = RX(I,1)                                                        
           RY1 = RY(I,1)                                                        
           RZ1 = RZ(I,1)                                                        
           PX1 = 2.0 * RX1 - OX(I,1) + DTSQ * FX(I,1) * OM1                     
           PY1 = 2.0 * RY1 - OY(I,1) + DTSQ * FY(I,1) * OM1                     
           PZ1 = 2.0 * RZ1 - OZ(I,1) + DTSQ * FZ(I,1) * OM1                     
           RX2 = RX(I,2)                                                        
           RY2 = RY(I,2)                                                        
           RZ2 = RZ(I,2)                                                        
           PX2 = 2.0 * RX2 - OX(I,2) + DTSQ * FX(I,2) * OM2                     
           PY2 = 2.0 * RY2 - OY(I,2) + DTSQ * FY(I,2) * OM2                     
           PZ2 = 2.0 * RZ2 - OZ(I,2) + DTSQ * FZ(I,2) * OM2                     
           RX3 = RX(I,3)                                                        
           RY3 = RY(I,3)                                                        
           RZ3 = RZ(I,3)                                                        
           PX3 = 2.0 * RX3 - OX(I,3) + DTSQ * FX(I,3) * OM3                     
           PY3 = 2.0 * RY3 - OY(I,3) + DTSQ * FY(I,3) * OM3                     
           PZ3 = 2.0 * RZ3 - OZ(I,3) + DTSQ * FZ(I,3) * OM3                     
                                                                                
C       ** CALCULATE RELATIVE VECTORS **                                        
                                                                                
           RX12 = RX1 - RX2                                                     
           RX12 = RX12 - ANINT ( RX12 * BOXINV ) * BOX                          
           RY12 = RY1 - RY2                                                     
           RY12 = RY12 - ANINT ( RY12 * BOXINV ) * BOX                          
           RZ12 = RZ1 - RZ2                                                     
           RZ12 = RZ12 - ANINT ( RZ12 * BOXINV ) * BOX                          
           RX23 = RX2 - RX3                                                     
           RX23 = RX23 - ANINT ( RX23 * BOXINV ) * BOX                          
           RY23 = RY2 - RY3                                                     
           RY23 = RY23 - ANINT ( RY23 * BOXINV ) * BOX                          
           RZ23 = RZ2 - RZ3                                                     
           RZ23 = RZ23 - ANINT ( RZ23 * BOXINV ) * BOX                          
           RX31 = RX3 - RX1                                                     
           RX31 = RX31 - ANINT ( RX31 * BOXINV ) * BOX                          
           RY31 = RY3 - RY1                                                     
           RY31 = RY31 - ANINT ( RY31 * BOXINV ) * BOX                          
           RZ31 = RZ3 - RZ1                                                     
           RZ31 = RZ31 - ANINT ( RZ31 * BOXINV ) * BOX                          
                                                                                
           PX12 = PX1 - PX2                                                     
           PX12 = PX12 - ANINT ( PX12 * BOXINV ) * BOX                          
           PY12 = PY1 - PY2                                                     
           PY12 = PY12 - ANINT ( PY12 * BOXINV ) * BOX                          
           PZ12 = PZ1 - PZ2                                                     
           PZ12 = PZ12 - ANINT ( PZ12 * BOXINV ) * BOX                          
           PX23 = PX2 - PX3                                                     
           PX23 = PX23 - ANINT ( PX23 * BOXINV ) * BOX                          
           PY23 = PY2 - PY3                                                     
           PY23 = PY23 - ANINT ( PY23 * BOXINV ) * BOX                          
           PZ23 = PZ2 - PZ3                                                     
           PZ23 = PZ23 - ANINT ( PZ23 * BOXINV ) * BOX                          
           PX31 = PX3 - PX1                                                     
           PX31 = PX31 - ANINT ( PX31 * BOXINV ) * BOX                          
           PY31 = PY3 - PY1                                                     
           PY31 = PY31 - ANINT ( PY31 * BOXINV ) * BOX                          
           PZ31 = PZ3 - PZ1                                                     
           PZ31 = PZ31 - ANINT ( PZ31 * BOXINV ) * BOX                          
                                                                                
C       ** CALCULATE SCALAR PRODUCTS **                                         
                                                                                
           R12R23 = RX12 * RX23 + RY12 * RY23 + RZ12 * RZ23                     
           R23R31 = RX23 * RX31 + RY23 * RY31 + RZ23 * RZ31                     
           R31R12 = RX31 * RX12 + RY31 * RY12 + RZ31 * RZ12                     
           P12SQ  = PX12 ** 2 + PY12 ** 2 + PZ12 ** 2                           
           P23SQ  = PX23 ** 2 + PY23 ** 2 + PZ23 ** 2                           
           P31SQ  = PX31 ** 2 + PY31 ** 2 + PZ31 ** 2                           
           P12R12 = PX12 * RX12 + PY12 * RY12 + PZ12 * RZ12                     
           P12R23 = PX12 * RX23 + PY12 * RY23 + PZ12 * RZ23                     
           P12R31 = PX12 * RX31 + PY12 * RY31 + PZ12 * RZ31                     
           P23R12 = PX23 * RX12 + PY23 * RY12 + PZ23 * RZ12                     
           P23R23 = PX23 * RX23 + PY23 * RY23 + PZ23 * RZ23                     
           P23R31 = PX23 * RX31 + PY23 * RY31 + PZ23 * RZ31                     
           P31R12 = PX31 * RX12 + PY31 * RY12 + PZ31 * RZ12                     
           P31R23 = PX31 * RX23 + PY31 * RY23 + PZ31 * RZ23                     
           P31R31 = PX31 * RX31 + PY31 * RY31 + PZ31 * RZ31                     
           CONST1 = R12SQ - P12SQ                                               
           CONST2 = R23SQ - P23SQ                                               
           CONST3 = R31SQ - P31SQ                                               
                                                                                
C       ** CALCULATE MATRIX AND QUADRATIC COEFFICIENTS **                       
                                                                                
           MAT(1,1) =  2.0 * ( OM1 + OM2 ) * P12R12                             
           MAT(1,2) = -2.0 *   OM2         * P12R23                             
           MAT(1,3) = -2.0 *   OM1         * P12R31                             
           MAT(2,1) = -2.0 *   OM2         * P23R12                             
           MAT(2,2) =  2.0 * ( OM2 + OM3 ) * P23R23                             
           MAT(2,3) = -2.0 *   OM3         * P23R31                             
           MAT(3,1) = -2.0 *   OM1         * P31R12                             
           MAT(3,2) = -2.0 *   OM3         * P31R23                             
           MAT(3,3) =  2.0 * ( OM1 + OM3 ) * P31R31                             
                                                                                
           Q111 = - R12SQ * ( OM1 + OM2 ) ** 2                                  
           Q122 = - R23SQ * OM2 ** 2                                            
           Q133 = - R31SQ * OM1 ** 2                                            
           Q112 = + 2.0 * R12R23 * ( OM1 + OM2 ) * OM2                          
           Q123 = - 2.0 * R23R31 * OM1 * OM2                                    
           Q131 = + 2.0 * R31R12 * ( OM1 + OM2 ) * OM1                          
                                                                                
           Q211 = - R12SQ * OM2 ** 2                                            
           Q222 = - R23SQ * ( OM2 + OM3 ) ** 2                                  
           Q233 = - R31SQ * OM3 ** 2                                            
           Q212 = + 2.0 * R12R23 * ( OM2 + OM3 ) * OM2                          
           Q223 = + 2.0 * R23R31 * ( OM2 + OM3 ) * OM3                          
           Q231 = - 2.0 * R31R12 * OM2 * OM3                                    
                                                                                
           Q311 = - R12SQ * OM1 ** 2                                            
           Q322 = - R23SQ * OM3 ** 2                                            
           Q333 = - R31SQ * ( OM1 + OM3 ) ** 2                                  
           Q312 = - 2.0 * R12R23 * OM1 * OM3                                    
           Q323 = + 2.0 * R23R31 * ( OM1 + OM3 ) * OM3                          
           Q331 = + 2.0 * R31R12 * ( OM1 + OM3 ) * OM1                          
                                                                                
C       ** NOW CALL ROUTINE TO INVERT THE MATRIX MAT **                         
                                                                                
           CALL MATINV ( MAT, INV )                                             
                                                                                
C       ** PREPARE FOR ITERATIVE LOOP **                                        
                                                                                
           DONE = .FALSE.                                                       
           IT   = 0                                                             
           L12  = 0.0                                                           
           L23  = 0.0                                                           
           L31  = 0.0                                                           
                                                                                
C       ** BEGIN ITERATIVE LOOP **                                              
                                                                                
1000       IF ( ( .NOT. DONE ) .AND. ( IT .LE. MAXIT ) ) THEN                   
                                                                                
              QUAD1 = Q111 * L12 ** 2 + Q112 * L12 * L23                        
     :              + Q122 * L23 ** 2 + Q123 * L23 * L31                        
     :              + Q133 * L31 ** 2 + Q131 * L31 * L12                        
              QUAD2 = Q211 * L12 ** 2 + Q212 * L12 * L23                        
     :              + Q222 * L23 ** 2 + Q223 * L23 * L31                        
     :              + Q233 * L31 ** 2 + Q231 * L31 * L12                        
              QUAD3 = Q311 * L12 ** 2 + Q312 * L12 * L23                        
     :              + Q322 * L23 ** 2 + Q323 * L23 * L31                        
     :              + Q333 * L31 ** 2 + Q331 * L31 * L12                        
                                                                                
              VEC1 = CONST1 + QUAD1                                             
              VEC2 = CONST2 + QUAD2                                             
              VEC3 = CONST3 + QUAD3                                             
                                                                                
C          ** OBTAIN SOLUTIONS OF LINEARIZED EQUATION **                        
                                                                                
              L12NEW = INV(1,1) * VEC1 + INV(1,2) * VEC2                        
     :               + INV(1,3) * VEC3                                          
              L23NEW = INV(2,1) * VEC1 + INV(2,2) * VEC2                        
     :               + INV(2,3) * VEC3                                          
              L31NEW = INV(3,1) * VEC1 + INV(3,2) * VEC2                        
     :               + INV(3,3) * VEC3                                          
                                                                                
              DONE = ( ( ABS ( L12NEW - L12 ) .LE. LTOL ) .AND.                 
     :                 ( ABS ( L23NEW - L23 ) .LE. LTOL ) .AND.                 
     :                 ( ABS ( L31NEW - L31 ) .LE. LTOL ) )                     
                                                                                
              L12 = L12NEW                                                      
              L23 = L23NEW                                                      
              L31 = L31NEW                                                      
              IT  = IT + 1                                                      
                                                                                
              GOTO 1000                                                         
                                                                                
           ENDIF                                                                
                                                                                
C       ** END OF ITERATION **                                                  
                                                                                
           PX1 = PX1 + OM1 * ( L12 * RX12 - L31 * RX31 )                        
           PY1 = PY1 + OM1 * ( L12 * RY12 - L31 * RY31 )                        
           PZ1 = PZ1 + OM1 * ( L12 * RZ12 - L31 * RZ31 )                        
           PX2 = PX2 + OM2 * ( L23 * RX23 - L12 * RX12 )                        
           PY2 = PY2 + OM2 * ( L23 * RY23 - L12 * RY12 )                        
           PZ2 = PZ2 + OM2 * ( L23 * RZ23 - L12 * RZ12 )                        
           PX3 = PX3 + OM3 * ( L31 * RX31 - L23 * RX23 )                        
           PY3 = PY3 + OM3 * ( L31 * RY31 - L23 * RY23 )                        
           PZ3 = PZ3 + OM3 * ( L31 * RZ31 - L23 * RZ23 )                        
                                                                                
           IF ( .NOT. DONE ) THEN                                               
                                                                                
              WRITE(*,'('' TOO MANY CONSTRAINT ITERATIONS '')')                 
              WRITE(*,'('' MOLECULE '',I5)') I                                  
              STOP                                                              
                                                                                
           ENDIF                                                                
                                                                                
C       ** CALCULATE KINETIC ENERGY CONTRIBUTION **                             
                                                                                
           VXIA = 0.5 * ( PX1 - OX(I,1) ) / DT                                  
           VYIA = 0.5 * ( PY1 - OY(I,1) ) / DT                                  
           VZIA = 0.5 * ( PZ1 - OZ(I,1) ) / DT                                  
           K = K + ( VXIA ** 2 + VYIA ** 2 + VZIA ** 2 ) * M(1)                 
           VXIA = 0.5 * ( PX2 - OX(I,2) ) / DT                                  
           VYIA = 0.5 * ( PY2 - OY(I,2) ) / DT                                  
           VZIA = 0.5 * ( PZ2 - OZ(I,2) ) / DT                                  
           K = K + ( VXIA ** 2 + VYIA ** 2 + VZIA ** 2 ) * M(2)                 
           VXIA = 0.5 * ( PX3 - OX(I,3) ) / DT                                  
           VYIA = 0.5 * ( PY3 - OY(I,3) ) / DT                                  
           VZIA = 0.5 * ( PZ3 - OZ(I,3) ) / DT                                  
           K = K + ( VXIA ** 2 + VYIA ** 2 + VZIA ** 2 ) * M(3)                 
                                                                                
C       ** CALCULATE CONSTRAINT VIRIAL CONTRIBUTION **                          
                                                                                
           WC = WC + L12 * R12SQ + L23 * R23SQ + L31 * R31SQ                    
                                                                                
C       ** STORE RESULTS **                                                     
                                                                                
           OX(I,1) = RX1                                                        
           OY(I,1) = RY1                                                        
           OZ(I,1) = RZ1                                                        
           RX(I,1) = PX1                                                        
           RY(I,1) = PY1                                                        
           RZ(I,1) = PZ1                                                        
           OX(I,2) = RX2                                                        
           OY(I,2) = RY2                                                        
           OZ(I,2) = RZ2                                                        
           RX(I,2) = PX2                                                        
           RY(I,2) = PY2                                                        
           RZ(I,2) = PZ2                                                        
           OX(I,3) = RX3                                                        
           OY(I,3) = RY3                                                        
           OZ(I,3) = RZ3                                                        
           RX(I,3) = PX3                                                        
           RY(I,3) = PY3                                                        
           RZ(I,3) = PZ3                                                        
                                                                                
2000    CONTINUE                                                                
                                                                                
C    ** END OF LOOP OVER MOLECULES **                                           
                                                                                
        K  = 0.5 * K                                                            
        WC = WC / DTSQ / 3.0                                                    
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE MATINV ( MAT, INV )                                          
                                                                                
C    *******************************************************************        
C    ** A SIMPLE ROUTINE TO INVERT THE 3X3 MATRIX MAT AND RETURN INV  **        
C    *******************************************************************        
                                                                                
        REAL        MAT(3,3), INV(3,3), DETERM                                  
        REAL        TOL                                                         
        PARAMETER ( TOL = 1.E-7 )                                               
                                                                                
C    *******************************************************************        
                                                                                
C    ** GET ALL SIGNED COFACTORS, TRANSPOSE, AND PUT IN INV **                  
                                                                                
        INV(1,1) = MAT(2,2) * MAT(3,3) - MAT(2,3) * MAT(3,2)                    
        INV(2,1) = MAT(3,1) * MAT(2,3) - MAT(3,3) * MAT(2,1)                    
        INV(3,1) = MAT(2,1) * MAT(3,2) - MAT(2,2) * MAT(3,1)                    
        INV(1,2) = MAT(3,2) * MAT(1,3) - MAT(3,3) * MAT(1,2)                    
        INV(2,2) = MAT(1,1) * MAT(3,3) - MAT(1,3) * MAT(3,1)                    
        INV(3,2) = MAT(3,1) * MAT(1,2) - MAT(3,2) * MAT(1,1)                    
        INV(1,3) = MAT(1,2) * MAT(2,3) - MAT(1,3) * MAT(2,2)                    
        INV(2,3) = MAT(2,1) * MAT(1,3) - MAT(2,3) * MAT(1,1)                    
        INV(3,3) = MAT(1,1) * MAT(2,2) - MAT(1,2) * MAT(2,1)                    
                                                                                
C    ** GET DETERMINANT AND GUARD AGAINST ZERO **                               
                                                                                
        DETERM = MAT(1,1) * INV(1,1) + MAT(1,2) * INV(2,1)                      
     :         + MAT(1,3) * INV(3,1)                                            
                                                                                
        IF ( ABS ( DETERM ) .LT. TOL ) STOP 'ZERO DETERM IN MATINV'             
                                                                                
        INV(1,1) = INV(1,1) / DETERM                                            
        INV(1,2) = INV(1,2) / DETERM                                            
        INV(1,3) = INV(1,3) / DETERM                                            
        INV(2,1) = INV(2,1) / DETERM                                            
        INV(2,2) = INV(2,2) / DETERM                                            
        INV(2,3) = INV(2,3) / DETERM                                            
        INV(3,1) = INV(3,1) / DETERM                                            
        INV(3,2) = INV(3,2) / DETERM                                            
        INV(3,3) = INV(3,3) / DETERM                                            
                                                                                
        RETURN                                                                  
        END                                                                     
