C     ***  B M F I T   ***  BEST  MOLECULAR  FIT   ***  S. C. NYBURG  **        
C     *************************   UPDATE     OCT   1981      ***********        
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
      CHARACTER*8 TYPE                                                          
      CHARACTER*8  PTYPE                                                        
      INTEGER*2 ORD                                                             
      COMMON TYPE(2,99),N,J,DS,NA,NB                                            
      COMMON D(99),  U(3,3),W(99),X(2,99),Y(2,99),Z(2,99),PTYPE(99)             
      DIMENSION CQ(3,3),A(9),R(9),ORD(99),XP(99),YP(99),ZP(99),EPX(99)          
      DIMENSION      CA(2),CB(2),CG(2),SG(2),CR(2),CP(2),EX(2,99)               
      DIMENSION EY(2,99),EZ(2,99),FRX(2,99),FRY(2,99),FRZ(2,99),RLC(2,3)        
      DIMENSION XS(2,99),YS(2,99),ZS(2,99),RAN(2,3),EPY(99),EPZ(99)             
      DIMENSION TITLE(18)
      OPEN(1, FILE=' ')
      OPEN(2, FILE=' ')
      OPEN(7, FILE=' ')
      OPEN(8, FILE=' ')
 1000 READ(7,79,END=999)TITLE                                                   
   79 FORMAT(18A4)                                                              
      DO 100 NY=1,2                                                             
      DO 100 NZ=1,50                                                            
      X(NY,NZ)=0.                                                               
      Y(NY,NZ)=0.                                                               
100   Z(NY,NZ)=0.                                                               
      DS=0.                                                                     
      J=0                                                                       
C     ITAP IS TH INPUT FILE HNIT FOR MOLECULE A IF=1,FOR B IF = 2               
      READ(7,*)N,NA,NB,NFIX,NFIT,NW,ITEMP,NDISCA,NDISCB,NCUTA,NCUTB,NNN         
C** 1 FORMAT(12I5)                                                              
      IF(NA-NB)101,101,103                                                      
  101 NTOT=NB                                                                   
      GO TO 10                                                                  
  103 NTOT=NA                                                                   
   10 CONTINUE                                                                  
      DPR=57.2957795                                                            
      DO 33 K=1,2                                                               
C     ******  RAN IN EITHER COSINES OR DEGREES                                  
      READ(7,* )RLC(K,1),RLC(K,2),RLC(K,3),RAN(K,1),RAN(K,2),RAN(K,3)           
C**12 FORMAT(6F10.6)                                                            
      IF (RAN(K,1).LT.1.0) GO TO 102                                            
      RAN(K,1)=COS(RAN(K,1)/DPR)                                                
      RAN(K,2)=COS(RAN(K,2)/DPR)                                                
      RAN(K,3)=COS(RAN(K,3)/DPR)                                                
  102 CA(K) = RAN(K,1)                                                          
      CB(K) = RAN(K,2)                                                          
      CG(K) = RAN(K,3)                                                          
      SG(K) = SQRT(1.0 - RAN(K,3)**2)                                           
      CR(K)=(CA(K)-CB(K)*CG(K))/SG(K)                                           
      CP(K)=SQRT(1.-CB(K)*CB(K)-CR(K)*CR(K))                                    
      IF(K-1)22,35,22                                                           
   35 NOW=NA                                                                    
      NOWW=NDISCA                                                               
      ITAP=1                                                                    
      IF(NCUTA.EQ.0)GO TO 2333                                                  
      II=0                                                                      
 2033 II=II+1                                                                   
      READ(ITAP,*)                                                              
C2133 FORMAT(A4)                                                                
      IF(II.EQ.NCUTA)GO TO 2333                                                 
      GO TO 2033                                                                
 2333 CONTINUE                                                                  
      GO TO 25                                                                  
   22 NOW=NB                                                                    
      NOWW=NDISCB                                                               
      ITAP=2                                                                    
      IF(NCUTB.EQ.0)GO TO 3333                                                  
      II=0                                                                      
 3033 II=II+1                                                                   
      READ(ITAP,*)                                                              
C3133 FORMAT(A4)                                                                
      IF(II.EQ.NCUTB)GO TO 3333                                                 
      GO TO 3033                                                                
 3333 CONTINUE                                                                  
   25 DO 33 KJ=1,NOW                                                            
      IF(KJ.GT.NOWW)GO TO 4033                                                  
      READ(ITAP,*)TYPE(K,KJ),FRX(K,KJ),FRY(K,KJ),FRZ(K,KJ)                      
    3 FORMAT (A6,21X,3F9.6 )                                                    
      IF(ITEMP.NE.0) READ(ITAP,*)                                               
      GO TO 4032                                                                
 4033 READ(7,*)TYPE(K,KJ),FRX(K,KJ),FRY(K,KJ),FRZ(K,KJ)                         
      IF(ITEMP.NE.0) READ(7,3)                                                  
 4032 CONTINUE                                                                  
      X(K,KJ)=FRX(K,KJ)*RLC(K,1)+FRY(K,KJ)*RLC(K,2)*CG(K)+FRZ(K,KJ)*            
     1RLC(K,3)*CB(K)                                                            
      Y(K,KJ)=FRY(K,KJ)*RLC(K,2)*SG(K)+FRZ(K,KJ)*RLC(K,3)*CR(K)                 
      Z(K,KJ)=FRZ(K,KJ)*RLC(K,3)*CP(K)                                          
      IF(NW)44,44,45                                                            
   44 W(KJ)=1.                                                                  
      GO TO 33                                                                  
   45 READ(7,* )EX(K,KJ),EY(K,KJ),EZ(K,KJ)                                     
      EX(K,KJ)=EX(K,KJ)*RLC(K,1)+EY(K,KJ)*RLC(K,2)*CG(K)+EZ(K,KJ)*              
     1RLC(K,3)*CB(K)                                                            
      EY(K,KJ)=EY(K,KJ)*RLC(K,2)*SG(K)+EZ(K,KJ)*RLC(K,3)*CR(K)                  
      EZ(K,KJ)=EZ(K,KJ)*RLC(K,3)*CP(K)                                          
C**46 FORMAT( 27X,3F9.6)                                                        
      EX(K,KJ)=SQRT((EX(K,KJ)**2+EY(K,KJ)**2+EZ(K,KJ)**2)/3.)                   
   33 CONTINUE                                                                  
      IF(NNN.EQ.0)GO TO 7                                                       
      READ(7,*  )ORD                                                            
C*105 FORMAT(24I3)                                                              
      DO 8 NORD=1,NB                                                            
      IF(ORD(NORD).EQ.0)ORD(NORD)=NORD                                          
      PTYPE(NORD)=TYPE(2,ORD(NORD))                                             
      XP(NORD)=X(2,ORD(NORD))                                                   
      YP(NORD)=Y(2,ORD(NORD))                                                   
      ZP(NORD)=Z(2,ORD(NORD))                                                   
      EPX(NORD)=EX(2,ORD(NORD))                                                 
      EPY(NORD)=EY(2,ORD(NORD))                                                 
      EPZ(NORD)=EZ(2,ORD(NORD))                                                 
    8 CONTINUE                                                                  
      DO 9 M=1,NB                                                               
      TYPE(2,M)=PTYPE(M)                                                        
      X(2,M)=XP(M)                                                              
      Y(2,M)=YP(M)                                                              
      Z(2,M)=ZP(M)                                                              
      EX(2,M)=EPX(M)                                                            
      EY(2,M)=EPY(M)                                                            
      EZ(2,M)=EPZ(M)                                                            
    9 CONTINUE                                                                  
C     FIND WEIGHT FOR EACH MATCHED ATOM PAIR                                    
    7 IF(NW.EQ.0) GO TO 49                                                      
      SES=0.                                                                    
      DO 47 KJ=1,N                                                              
      EX(1,KJ)=SQRT((EX(1,KJ))**2+(EX(2,KJ))**2)                                
   47 SES=SES+ (1./(EX(1,KJ)**2))                                               
      AN=N                                                                      
      AN=AN/SES                                                                 
      DO 48 KJ=1,N                                                              
   48 W(KJ)=AN/(EX(1,KJ)**2)                                                    
   49 WRITE(8,179)TITLE                                                         
  179 FORMAT(1H1,20A4,//)                                                       
      IF(NW.EQ.0)  GO TO 75                                                     
      WRITE(8,76)                                                               
   76 FORMAT(31X,'THIS IS A WEIGHTED FIT'//)                                    
      GO TO 80                                                                  
   75 WRITE(8,77)                                                               
   77 FORMAT(30X,'THIS IS AN UNWEIGHTED FIT'//)                                 
   80 WRITE(8,23)                                                               
   23 FORMAT(72H                         ORTHOGONAL COORDINATES OF ATOMS        
     1 BEFORE MATCHING//)                                                       
      WRITE(8,24)                                                               
   24 FORMAT(21X,'FIRST MOLECULE',31X,'SECOND MOLECULE'//)                      
      WRITE(8,20)RLC(1,1),RLC(1,2),RLC(1,3),RLC(2,1),RLC(2,2),RLC(2,3)          
   20 FORMAT(3X,'A=',F9.4,'   B=',F9.4,'   C=',F9.4,10X,'A=',F9.4,'   B=        
     1',F9.4,'   C=',F9.4)                                                      
      WRITE(8,43)RAN(1,1),RAN(1,2),RAN(1,3),RAN(2,1),RAN(2,2),RAN(2,3)          
   43 FORMAT(3X,'COSA=',F9.6,',COSB=',F9.6,',COSG=',F9.6,5X,'COSA=',F9.6        
     1,',COSB=',F9.6,',COSG=',F9.6//)                                           
      CALL NATOMS                                                               
C      TYPE OF FORMAT                                                           
   60 IF((NA-N).EQ.0)GO TO 50                                                   
      IF((NB-N).EQ.0)GO TO 54                                                   
      IF(NA-NB)51,52,53                                                         
C     OUTPUT ROUTINE F5                                                         
   51 NP=N+1                                                                    
      NAP=NA+1                                                                  
      IF(J.NE.0 )GO TO 84                                                       
      WRITE(8,26)(MI,TYPE(1,MI),X(1,MI),Y(1,MI),Z(1,MI),TYPE(2,MI),X(2,M        
     1I),Y(2,MI),Z(2,MI),MI=NP,NA)                                              
   26 FORMAT(I4,A6,3F10.6,3X,A6,3F10.6)                                         
      WRITE(8,74)(MP,TYPE(2,MP),X(2,MP),Y(2,MP),Z(2,MP),MP=NAP,NB)              
   74 FORMAT(I4,39X,A6,3F10.6)                                                  
      GO TO 34                                                                  
   84 WRITE(8,82)(MM,TYPE(1,MM),X(1,MM),Y(1,MM),Z(1,MM),TYPE(2,MM),X(2,M        
     1M),Y(2,MM),Z(2,MM),D(MM),MM=NP,NA)                                        
   82 FORMAT(I4,A6,3F10.6,3X,A6,4F10.6)                                         
   85 WRITE(8,42)(MR,TYPE(2,MR),X(2,MR),Y(2,MR),Z(2,MR),MR=NAP,NB)              
   42 FORMAT(I4,39X,A6,3F10.6)                                                  
   34 CONTINUE                                                                  
      GO TO 55                                                                  
C     OUTPUT ROUTINE F6                                                         
   52 NP=N+1                                                                    
      IF(J.NE.0 )GO TO 86                                                       
      WRITE(8,27)(MS,TYPE(1,MS),X(1,MS),Y(1,MS),Z(1,MS),TYPE(2,MS),X(2,M        
     1S),Y(2,MS),Z(2,MS),MS=NP,NA)                                              
   27 FORMAT(I4,A6,3F10.6,3X,A6,3F10.6)                                         
      GO TO 87                                                                  
   86 WRITE(8,36)(MT,TYPE(1,MT),X(1,MT),Y(1,MT),Z(1,MT),TYPE(2,MT),X(2,M        
     1T),Y(2,MT),Z(2,MT),D(MT),MT=NP,NA)                                        
   36 FORMAT(I4,A6,3F10.6,3X,A6,4F10.6)                                         
   87 CONTINUE                                                                  
      GO TO 55                                                                  
C      OUTPUT ROUTINE F4                                                        
   53 NP=N+1                                                                    
      NBP=NB+1                                                                  
      IF(J.NE.0 )GO TO 81                                                       
      WRITE(8,37)(MG,TYPE(1,MG),X(1,MG),Y(1,MG),Z(1,MG),TYPE(2,MG),X(2,M        
     1G),Y(2,MG),Z(2,MG),MG=NP,NB)                                              
   37 FORMAT(I4,A6,3F10.6,3X,A6,3F10.6)                                         
      WRITE(8,58)(MH,TYPE(1,MH),X(1,MH),Y(1,MH),Z(1,MH),MH=NBP,NA)              
      GO TO 135                                                                 
   81 WRITE(8,39)(MN,TYPE(1,MN),X(1,MN),Y(1,MN),Z(1,MN),TYPE(2,MN),X(2,M        
     1N),Y(2,MN),Z(2,MN),D(MN),MN=NP,NB)                                        
   39 FORMAT(I4,A6,3F10.6,3X,A6,4F10.6)                                         
   83 WRITE(8,58)(ML,TYPE(1,ML),X(1,ML),Y(1,ML),Z(1,ML),ML=NBP,NA)              
   58 FORMAT(I4,A6,3F10.6)                                                      
  135 CONTINUE                                                                  
      GO TO 55                                                                  
   50 IF(NB-N)66,65,66                                                          
   65 WRITE(8,29)                                                               
   29 FORMAT(33H                             NONE//)                            
      GO TO 55                                                                  
C     OUTPUT ROUTINE F3                                                         
   66 NBP=N+1                                                                   
      WRITE(8,40)(MF,TYPE(2,MF),X(2,MF),Y(2,MF),Z(2,MF),MF=NBP,NB)              
   40 FORMAT(I4,39X,A6,3F10.6)                                                  
      GO TO 55                                                                  
C      OUTPUT ROUTINE F2                                                        
   54 NAP=N+1                                                                   
      WRITE(8,41)(ME,TYPE(1,ME),X(1,ME),Y(1,ME),Z(1,ME),ME=NAP,NA)              
   41 FORMAT(I4,A6,3F10.6)                                                      
   55 IF(J.NE.0)GO TO 61                                                        
      WRITE(8,90)                                                               
   90 FORMAT(1H0)                                                               
      WRITE(8,93)                                                               
   93 FORMAT(69H                      CENTROID                                  
     1     CENTROID/)                                                           
C      FIND PAIR-WEIGHTED CENTROIDS OF MOLECULES                                
      IF(NFIX.NE.0)GO TO 62                                                     
      XA=0.                                                                     
      YA=0.                                                                     
      ZA=0.                                                                     
      XB=0.                                                                     
      YB=0.                                                                     
      ZB=0.                                                                     
      DO 4 NF=1,N                                                               
      XA=XA+X(1,NF)*W(NF)                                                       
      YA=YA+Y(1,NF)*W(NF)                                                       
      ZA=ZA+Z(1,NF)*W(NF)                                                       
      XB=XB+X(2,NF)*W(NF)                                                       
      YB=YB+Y(2,NF)*W(NF)                                                       
    4 ZB=ZB+Z(2,NF)*W(NF)                                                       
      AN=N                                                                      
      XA=XA/AN                                                                  
      YA=YA/AN                                                                  
      ZA=ZA/AN                                                                  
      XB=XB/AN                                                                  
      YB=YB/AN                                                                  
      ZB=ZB/AN                                                                  
      GO TO 63                                                                  
   62 XA=X(1,NFIX)                                                              
      YA=Y(1,NFIX)                                                              
      ZA=Z(1,NFIX)                                                              
      XB=X(2,NFIX)                                                              
      YB=Y(2,NFIX)                                                              
      ZB=Z(2,NFIX)                                                              
   63 WRITE(8,94)XA,YA,ZA,XB,YB,ZB                                              
   94 FORMAT(10X,3F10.6,9X,3F10.6,//)                                           
      WRITE(8,78)                                                               
   78 FORMAT(1H1)                                                               
      IF(NW.EQ.0)  GO TO 175                                                    
      WRITE(8,76)                                                               
      GO TO 180                                                                 
  175 WRITE(8,77)                                                               
C     MOVE CENTROIDS TO ORIGIN                                                  
  180 DO 5 NG=1,NTOT                                                            
      X(1,NG)=X(1,NG)-XA                                                        
      Y(1,NG)=Y(1,NG)-YA                                                        
      Z(1,NG)=Z(1,NG)-ZA                                                        
      X(2,NG)=X(2,NG)-XB                                                        
      Y(2,NG)=Y(2,NG)-YB                                                        
      Z(2,NG)=Z(2,NG)-ZB                                                        
      XS(2,NG)=X(2,NG)                                                          
      YS(2,NG)=Y(2,NG)                                                          
    5 ZS(2,NG)=Z(2,NG)                                                          
      DO 997 I=1,3                                                              
      DO 997 JJJ=1,3                                                            
      U(I,JJJ)=0.                                                               
  997 CQ(I,JJJ)=0.                                                              
      DO 17 I=1,N                                                               
      CQ(1,1)=CQ(1,1)+X(1,I )*X(2,I )*W(I)                                      
      CQ(1,2)=CQ(1,2)+Y(1,I )*X(2,I )*W(I)                                      
      CQ(1,3)=CQ(1,3)+Z(1,I )*X(2,I )*W(I)                                      
      CQ(2,1)=CQ(2,1)+X(1,I )*Y(2,I )*W(I)                                      
      CQ(2,2)=CQ(2,2)+Y(1,I )*Y(2,I )*W(I)                                      
      CQ(2,3)=CQ(2,3)+Z(1,I )*Y(2,I )*W(I)                                      
      CQ(3,1)=CQ(3,1)+X(1,I )*Z(2,I )*W(I)                                      
      CQ(3,2)=CQ(3,2)+Y(1,I )*Z(2,I )*W(I)                                      
      CQ(3,3)=CQ(3,3)+Z(1,I )*Z(2,I )*W(I)                                      
   17 CONTINUE                                                                  
      A(1)   =CQ(1,1)*CQ(1,1)+CQ(1,2)*CQ(1,2)+CQ(1,3)*CQ(1,3)                   
      A(2)   =CQ(1,1)*CQ(2,1)+CQ(1,2)*CQ(2,2)+CQ(1,3)*CQ(2,3)                   
      A(4)   =CQ(1,1)*CQ(3,1)+CQ(1,2)*CQ(3,2)+CQ(1,3)*CQ(3,3)                   
      A(3)   =CQ(2,1)*CQ(2,1)+CQ(2,2)*CQ(2,2)+CQ(2,3)*CQ(2,3)                   
      A(5)   =CQ(2,1)*CQ(3,1)+CQ(2,2)*CQ(3,2)+CQ(2,3)*CQ(3,3)                   
      A(6)   =CQ(3,1)*CQ(3,1)+CQ(3,2)*CQ(3,2)+CQ(3,3)*CQ(3,3)                   
      J=1                                                                       
      CALL EIGEN(A,R)                                                           
      NWR=0                                                                     
      NMIR=0                                                                    
      CALL VECTOR(A,R,CQ,NWR,NMIR,DETA)                                         
  213 IF(NWR.NE.0) GO TO 204                                                    
      IF(NMIR.EQ.1) GO TO 209                                                   
      IF(DETA.LT.0.) GO TO 203                                                  
      WRITE(8,211)                                                              
  211 FORMAT(20X,'THE BEST FIT IS GIVEN WITH THE SECOND MOLECULE UN-MIRR        
     1ORED',///)                                                                
      GO TO 207                                                                 
  203 WRITE(8,212)                                                              
  212 FORMAT(20X,'THE BEST FIT IS GIVEN WITH THE SECOND MOLECULE MIRRORE        
     1D',///)                                                                   
      GO TO 207                                                                 
  204 IF(DDET.LT.0.) GO TO 205                                                  
      WRITE(8,206)                                                              
  206 FORMAT(1H1,35X,'THIS IS THE MIRRORED FIT',///)                            
      GO TO 2                                                                   
  205 WRITE(8,208)                                                              
  208 FORMAT(1H1,20X,'THIS IS THE FIT WITH THE SECOND MOLECULE UN-MIRROR        
     1ED',///)                                                                  
      GO TO 2                                                                   
  209 WRITE(8,210)                                                              
  210 FORMAT( 1X,'EITHER ONE OF THE MOLECULES IS FLAT OR THE FIT IS THE         
     1SAME WHETHER OR NOT THE 2ND MOLECULE IS MIRRORED',///)                    
  207 DDET=DETA                                                                 
    2 DO 994 I=1,NB                                                             
      XN=A(1)*XS(2,I)+A(2)*YS(2,I)+A(3)*ZS(2,I)                                 
      YN=A(4)*XS(2,I)+A(5)*YS(2,I)+A(6)*ZS(2,I)                                 
      ZN=A(7)*XS(2,I)+A(8)*YS(2,I)+A(9)*ZS(2,I)                                 
      X(2,I)=XN                                                                 
      Y(2,I)=YN                                                                 
      Z(2,I)=ZN                                                                 
  994 CONTINUE                                                                  
      CALL DSUM                                                                 
      WRITE(8,11)                                                               
   11 FORMAT(87H                                                                
     1                          DELTA/)                                         
      CALL NATOMS                                                               
      GO TO 60                                                                  
   61 CONTINUE                                                                  
C     THIS FINAL SECTION PERFORMS THE TRANSFORMATION OF 2ND MOLECULE INT        
C       THE FRACTIONAL COORDINATES OF THE FIRST.                                
C      IF THIS IS REQUIRED, PUT NFIT=1.  OTHERWISE NFIT=0.                      
      IF(NFIT.EQ.0) GO TO 161                                                   
C     BRING 2ND MOLECULE BACK TO SAME CRYSTAL FRAME AS 1ST MOLECULE             
C     MOVE 2ND MOLECULE TO CENTROID OF 1ST MOLECULE IN ITS ORIGINAL POSI        
      DO 16 NG=1,NB                                                             
      X(2,NG)=X(2,NG)+XA                                                        
      Y(2,NG)=Y(2,NG)+YA                                                        
      Z(2,NG)=Z(2,NG)+ZA                                                        
   16 CONTINUE                                                                  
C     INVERSE MATRIX TO BRING TO FRACTIONAL COORDS                              
      DO 30 K=1,NB                                                              
      X(2,K)=(X(2,K)-(Y(2,K)*CG(1)/SG(1))                                       
     1          +(Z(2,K)/CP(1))*((CR(1)*CG(1)/SG(1)) -CB(1)))/RLC(1,1)          
      Y(2,K)= (Y(2,K)/SG(1)-((Z(2,K)*CR(1))/(CP(1)*SG(1))))/RLC(1,2)            
   30 Z(2,K)=Z(2,K)/(CP(1)*RLC(1,3))                                            
      WRITE(8,250)                                                              
  250 FORMAT     ( 1H0,///,'   NO', 4X,' ATOM ',6X,'FRACT. COORDINATES',        
     1 11X,' ATOM   FRACT. COORDS. OF 2ND MOLECULE',/23X,'FIRST MOLECULE        
     2',26X,'IN CELL OF FIRST')                                                 
      IF(NA-NB)67,67,68                                                         
   67 NUM=NA                                                                    
      GO TO 70                                                                  
   68 NUM=NB                                                                    
   70 WRITE(8,59)(M,TYPE(1,M),FRX(1,M),FRY(1,M),FRZ(1,M),TYPE(2,M),X(2,M        
     1),Y(2,M),Z(2,M),M=1,NUM)                                                  
   59 FORMAT(I4,5X,A6,3F10.6,5X,A6,3F10.6)                                      
      IF(NA-NB)71,161,72                                                        
   71 NAP=NA+1                                                                  
      WRITE(8,99)(M,TYPE(2,M),X(2,M),Y(2,M),Z(2,M),M=NAP,NB)                    
   99 FORMAT(I4,46X,A6,3F10.6)                                                  
      GO TO 161                                                                 
   72 NBP=NB+1                                                                  
      WRITE(8,73)(M,TYPE(1,M),FRX(1,M),FRY(1,M),FRZ(1,M),M=NBP,NA)              
   73 FORMAT(I4,5X,A6,3F10.6)                                                   
  161 IF(NMIR.EQ.1) GO TO 1000                                                  
      NWR=NWR+1                                                                 
      IF(NWR.EQ.2)GO TO 1000                                                    
      CALL VECTOR(A,R,CQ,NWR,NMIR,DETA)                                         
      GO TO 204                                                                 
  999 STOP                                                                      
      END                                                                       
      SUBROUTINE EIGEN(A,R)                                                     
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
      DIMENSION A(6),R(9)                                                       
      RANGE=1.0E-6                                                              
      IQ=-3                                                                     
      DO 20 J=1,3                                                               
      IQ=IQ+3                                                                   
      DO 20 I=1,3                                                               
      IJ=IQ+I                                                                   
      R(IJ)=0.0                                                                 
      IF(I-J) 20,15,20                                                          
   15 R(IJ)=1.0                                                                 
   20 CONTINUE                                                                  
   25 ANORM=0.0                                                                 
      DO 35 I=1,3                                                               
      DO 35 J=I,3                                                               
      IF(I-J) 30,35,30                                                          
   30 IA=I+(J*J-J)/2                                                            
      ANORM=ANORM+A(IA)*A(IA)                                                   
   35 CONTINUE                                                                  
      IF(ANORM)165,165,40                                                       
   40 ANORM=1.414*SQRT(ANORM)                                                   
      ANRMX=ANORM*RANGE/3.                                                      
      IND=0                                                                     
      THR=ANORM                                                                 
   45 THR=THR/3.                                                                
   50 L=1                                                                       
   55 M=L+1                                                                     
   60 MQ=(M*M-M)/2                                                              
      LQ=(L*L-L)/2                                                              
      LM=L+MQ                                                                   
   62 IF(ABS(A(LM))-THR) 130,65,65                                              
   65 IND=1                                                                     
      LL=L+LQ                                                                   
      MM=M+MQ                                                                   
      X =0.5*(A(LL)-A(MM))                                                      
   68 Y=-A(LM)/SQRT(A(LM)*A(LM)+X*X)                                            
      IF(X) 70,75,75                                                            
   70 Y=-Y                                                                      
   75 SINX=Y/SQRT(2.0*(1.0+(SQRT(1.0-Y*Y))))                                    
      SINX2=SINX*SINX                                                           
   78 COSX=SQRT(1.0-SINX2)                                                      
      COSX2=COSX*COSX                                                           
      SINCS=SINX*COSX                                                           
      ILQ=3*(L-1)                                                               
      IMQ=3*(M-1)                                                               
      DO 125 I=1,3                                                              
      IQ=(I*I-I)/2                                                              
      IF(I-L) 80,115,80                                                         
   80 IF(I-M) 85,115,90                                                         
   85 IM=I+MQ                                                                   
      GO TO 95                                                                  
   90 IM=M+IQ                                                                   
   95 IF(I-L) 100,105,105                                                       
  100 IL=I+LQ                                                                   
      GO TO 110                                                                 
  105 IL=L+IQ                                                                   
  110 X =A(IL)*COSX-A(IM)*SINX                                                  
      A(IM)=A(IL)*SINX+A(IM)*COSX                                               
      A(IL)=X                                                                   
  115 ILR=ILQ+I                                                                 
      IMR=IMQ+I                                                                 
      X =R(ILR)*COSX-R(IMR)*SINX                                                
      R(IMR)=R(ILR)*SINX+R(IMR)*COSX                                            
      R(ILR)=X                                                                  
  125 CONTINUE                                                                  
      X =2.0*A(LM)*SINCS                                                        
      Y =A(LL)*COSX2+A(MM)*SINX2-X                                              
      X =A(LL)*SINX2+A(MM)*COSX2+X                                              
      A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2)                             
      A(LL)=Y                                                                   
      A(MM)=X                                                                   
  130 IF(M-3) 135,140,135                                                       
  135 M=M+1                                                                     
      GO TO 60                                                                  
  140 IF(L-2) 145,150,145                                                       
  145 L=L+1                                                                     
      GO TO 55                                                                  
  150 IF(IND-1) 160,155,160                                                     
  155 IND=0                                                                     
      GO TO 50                                                                  
  160 IF(THR-ANRMX) 165,165,45                                                  
  165 IQ=-3                                                                     
      DO 185 I=1,3                                                              
      IQ=IQ+3                                                                   
      LL=I+(I*I-I)/2                                                            
      JQ=3*(I-2)                                                                
      DO 185 J=I,3                                                              
      JQ=JQ+3                                                                   
      MM=J+(J*J-J)/2                                                            
      IF(A(LL)-A(MM)) 170,185,185                                               
  170 X =A(LL)                                                                  
      A(LL)=A(MM)                                                               
      A(MM)=X                                                                   
      DO 180 K=1,3                                                              
      ILR=IQ+K                                                                  
      IMR=JQ+K                                                                  
      X =R(ILR)                                                                 
      R(ILR)=R(IMR)                                                             
  180 R(IMR)=X                                                                  
  185 CONTINUE                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE VECTOR(A,R,CQ,NWR,NMIR,DETA)                                   
      IMPLICIT REAL*8 (A-H, O-Z)                                                
      DIMENSION CQ(3,3)                                                         
      DIMENSION A(9),R(9)                                                       
      DIMENSION U(3,3)                                                          
      DIMENSION WV(3),ZZ(3,3),RC(3,3),B(3,3)                                    
      DO 28 I=1,3                                                               
      DO 28 J=1,3                                                               
   28 U(I,J)=0.                                                                 
      IF(NWR.NE.0) GO TO 52                                                     
      WV(1)=A(1)                                                                
      WV(2)=A(3)                                                                
      WV(3)=A(6)                                                                
      ZZ(1,1)=R(1)                                                              
      ZZ(1,2)=R(2)                                                              
      ZZ(1,3)=R(3)                                                              
      ZZ(2,1)=R(4)                                                              
      ZZ(2,2)=R(5)                                                              
      ZZ(2,3)=R(6)                                                              
      ZZ(3,1)=R(7)                                                              
      ZZ(3,2)=R(8)                                                              
      ZZ(3,3)=R(9)                                                              
      DO 21 I=1,3                                                               
      DO 21 J=1,3                                                               
      B(I,J)=0.                                                                 
   21 RC(I,J)=CQ(J,I)                                                           
C     MAKE SURE DET OF RC IS NOT EXACTLY ZERO                                   
  999 FORMAT(/'   VALUE OF DETERMINANT',E14.4)                                  
      DETR=  RC(1,1)*RC(2,2)*RC(3,3)                                            
     $      -RC(1,1)*RC(2,3)*RC(3,2)                                            
     $      +RC(1,2)*RC(2,3)*RC(3,1)                                            
     $      -RC(1,2)*RC(2,1)*RC(3,3)                                            
     $      +RC(1,3)*RC(2,1)*RC(3,2)                                            
     $      -RC(1,3)*RC(2,2)*RC(3,1)                                            
      WRITE(8,999)DETR                                                          
      IF(ABS(DETR).LT.0.0001 )GO TO 17                                          
      DO 19 K=1,3                                                               
      DO 19 I=1,3                                                               
      DO 19 J=1,3                                                               
      B(K,I)=B(K,I)+RC(I,J)*ZZ(K,J)                                             
   19 CONTINUE                                                                  
      GO TO 38                                                                  
   17 NMIR=1                                                                    
      IF(WV(2).LT.1.0E-6)   GO TO 39                                            
      DO 40 K=1,2                                                               
      DO 40 I=1,3                                                               
      DO 40 J=1,3                                                               
      B(K,I)=B(K,I)+RC(I,J)*ZZ(K,J)                                             
   40 CONTINUE                                                                  
      B (3,1)=B (1,2)*B (2,3)-B (1,3)*B (2,2)                                   
      B (3,2)=B (1,3)*B (2,1)-B (2,3)*B (1,1)                                   
      B (3,3)=B (1,1)*B (2,2)-B (2,1)*B (1,2)                                   
   38 BS=      B(1,1)**2+B(1,2)**2+B(1,3)**2                                    
      IF (BS.EQ.0.0) GO TO 100                                                  
      FACS=1./(B(1,1)**2+B(1,2)**2+B(1,3)**2)                                   
      DO 13 J=1,3                                                               
   13 B(1,J)=SQRT(FACS)*B(1,J)                                                  
  100 BS=      B(2,1)**2+B(2,2)**2+B(2,3)**2                                    
      IF(BS.EQ.0.0) GO TO 101                                                   
      FACS=1./(B(2,1)**2+B(2,2)**2+B(2,3)**2)                                   
      DO 14 J=1,3                                                               
   14 B(2,J)=SQRT(FACS)*B(2,J)                                                  
  101 BS=      B(3,1)**2+B(3,2)**2+B(3,3)**2                                    
      IF(BS.EQ.0.0) GO TO 102                                                   
      FACS=1./(B(3,1)**2+B(3,2)**2+B(3,3)**2)                                   
      DO 15 J=1,3                                                               
   15 B(3,J)=SQRT(FACS)*B(3,J)                                                  
      GO TO 102                                                                 
   39 DO 41 I=1,3                                                               
      DO 41 J=1,3                                                               
      B(1,I)=B(1,I)+RC(I,J)*ZZ(1,J)                                             
   41 CONTINUE                                                                  
      FACS=1./(B(1,1)**2+B(1,2)**2+B(1,3)**2)                                   
      DO 42 J=1,3                                                               
      B(1,J)=SQRT(FACS)*B(1,J)                                                  
42    CONTINUE                                                                  
      IF(B(1,1).EQ.0.) GO TO 43                                                 
      IF(B(1,2).EQ.0.) GO TO 37                                                 
      B(2,2)=B(1,1)/SQRT(B(1,1)**2+B(1,2)**2)                                   
      B(2,1)=-(B(1,2)*B(2,2))/B(1,1)                                            
      B(2,3)=0.                                                                 
      GO TO 44                                                                  
   43 B(2,1)=1.                                                                 
      B(2,2)=0.                                                                 
      B(2,3)=0.                                                                 
      GO TO 44                                                                  
   37 B(2,1)=0.                                                                 
      B(2,2)=1.                                                                 
      B(2,3)=0.                                                                 
   44 B(3,1)=B(1,2)*B(2,3)-B(1,3)*B(2,2)                                        
      B(3,2)=B(1,3)*B(2,1)-B(1,1)*B(2,3)                                        
      B(3,3)=B(1,1)*B(2,2)-B(1,2)*B(2,1)                                        
      GO TO 102                                                                 
52      B(3,1)=-B(3,1)                                                          
      B(3,2)=-B(3,2)                                                            
      B(3,3)=-B(3,3)                                                            
 102  DO 20 I=1,3                                                               
      DO 20 J=1,3                                                               
      DO 20 K=1,3                                                               
   20 U(I,J)=U(I,J)+B(K,I)*ZZ(K,J)                                              
      A(1)=U(1,1)                                                               
      A(2)=U(1,2)                                                               
      A(3)=U(1,3)                                                               
      A(4)=U(2,1)                                                               
      A(5)=U(2,2)                                                               
      A(6)=U(2,3)                                                               
      A(7)=U(3,1)                                                               
      A(8)=U(3,2)                                                               
      A(9)=U(3,3)                                                               
C     IS A MATRIX ORTHOGONAL?                                                   
      A11=A(1)*A(1)+A(2)*A(2)+A(3)*A(3)                                         
      A12=A(1)*A(4)+A(2)*A(5)+A(3)*A(6)                                         
      A13=A(1)*A(7)+A(2)*A(8)+A(3)*A(9)                                         
      A21=A(4)*A(1)+A(5)*A(2)+A(6)*A(3)                                         
      A22=A(4)*A(4)+A(5)*A(5)+A(6)*A(6)                                         
      A23=A(4)*A(7)+A(5)*A(8)+A(6)*A(9)                                         
      A31=A(7)*A(1)+A(8)*A(2)+A(9)*A(3)                                         
      A32=A(7)*A(4)+A(8)*A(5)+A(9)*A(6)                                         
      A33=A(7)*A(7)+A(8)*A(8)+A(9)*A(9)                                         
      WRITE(8,998)A11,A12,A13                                                   
      WRITE(8,998)A21,A22,A23                                                   
      WRITE(8,998)A31,A32,A33                                                   
 998  FORMAT(/20X,3E14.4,/)                                                     
      DETA=A(1)*A(5)*A(9)-A(1)*A(6)*A(8)+A(2)*A(6)*A(7)-A(2)*A(4)*A(9)+A        
     1(3)*A(4)*A(8)-A(3)*A(5)*A(7)                                              
      RETURN                                                                    
      END                                                                       
      SUBROUTINE DSUM                                                           
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
      CHARACTER*8 TYPE                                                          
      COMMON TYPE(2,99),N,J,DS,NA,NB                                            
      COMMON D(99),  U(3,3),W(99),X(2,99),Y(2,99),Z(2,99)                       
      DIMENSION A(9),R(9)                                                       
      DS=0.                                                                     
      IF(NA-NB)96,96,98                                                         
   96 NSUM=NA                                                                   
      GO TO 99                                                                  
   98 NSUM=NB                                                                   
   99 DO 60 MB=1,NSUM                                                           
   60 D(MB)=     (X(1,MB)-X(2,MB))**2+(Y(1,MB)-Y(2,MB))**2+(Z(1,MB              
     1)-Z(2,MB))**2                                                             
      DO 160 MU=1,N                                                             
  160 DS=DS+D(MU)                                                               
      DO 161 MV=1,NSUM                                                          
  161 D(MV)=SQRT(D(MV))                                                         
      RETURN                                                                    
      END                                                                       
      SUBROUTINE NATOMS                                                         
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
      CHARACTER*8 TYPE                                                          
      COMMON TYPE(2,99),N,J,DS,NA,NB                                            
      COMMON D(99),  U(3,3),W(99),X(2,99),Y(2,99),Z(2,99)                       
      IF(J.NE.0 )GO TO 56                                                       
      WRITE(8,26)(MA,TYPE(1,MA),X(1,MA),Y(1,MA),Z(1,MA),TYPE(2,MA),X(2,M        
     1A),Y(2,MA),Z(2,MA),MA=1,N)                                                
   26 FORMAT(I4,A6,3F10.6,12X,A6,3F10.6)                                        
      GO TO 57                                                                  
   56 WRITE(8,27)(MD,TYPE(1,MD),X(1,MD),Y(1,MD),Z(1,MD),TYPE(2,MD),X(2,M        
     1D),Y(2,MD),Z(2,MD),D(MD),MD=1,N)                                          
   27 FORMAT(I4,A6,3F10.6,3X,A6,4F10.6)                                         
      WRITE(8,164) DS                                                           
  164 FORMAT(/,'                                SUM OF DIFFERENCES SQUAR        
     1ED                               ',F10.6,//)                              
   57 WRITE (8,29)                                                              
   29 FORMAT(2H0 ,///)                                                          
      WRITE (8,28)                                                              
   28 FORMAT(56H                                         UNMATCHED ATOMS        
     1//)                                                                       
      RETURN                                                                    
      END                                                                       
