              SUBROUTINE BONVIB(IN,IOUT,LINE,MAXLIN,MARK)                   
C                                                                           
C      PROGRAM  TO READ  THE  DATA FILE  PREPARED BY THE PROGRAM            
C      REFINE, REFORMAT THE PARAMETER VARIANCE-COVARIANCE MATRIX            
C      AND THEN CALL THE  BOND DISTANCE AND  ANGLE SUBROUTINE OR            
C      THERMAL ELLIPSOID  SUBROUTINE  TO COMPUTE THESE FUNCTIONS            
C                       AND THEIR ERRORS.                                   
C                                                                           
C                 WRITTEN  BY  L. W. FINGER                                 
C                 MODIFIED BY  U. H. ZUCKER                                 
C                 MODIFIED BY  W. F. KUHS
C                                                                           
C---- DECLARATIONS                                                          
C                                                                           
      DIMENSION  INTER(20), X(35,4), Y(35,4),BDF(16),IBDF(2),EXTRAN(12),    
     & WINKEL(3),IPARA(1000), KPAR(60), KRPAR(60),KLAB(3),LABEL(3),S(3),    
     & TEMP(20),TEMP1(9,20),TEMP2(2,20),TENP(3),TEOP(5),ZT(5,60),JSF(10)    
     & ,CR(10,60),DR(15,60),PAR(60),PARNAM(2,60),TEMP4(6),ITEMP4(7),NCNT    
     & (4), ICNSTR(2), EXTIN(6), IPNT(6,60), JPNT(180),CNSTR(2),CV(400),    
     & CNSDGN(2),ALF(60),IALFA(60),LINE(11)                                 
C                                                                           
C---- ******************************************************************    
C                                                                           
C     MAXSIG (SEE BELOW)  HAS  TO  BE  THE  DIMENSION OF THE ARRAY SIGX     
C                        (COMMON ABLE)                                      
C                                                                           
C                  D I M E N S I O N    S I G X (128000)                       
C                                                                           
C---- ******************************************************************    
C                                                                           
      COMMON/ABLE/  XYZ(3,60),BETA(6,60),A(3),ANG(3),AS(3),CS(3),ATOM(2,    
     & 60), SIGA(6,6), SIGB(6,6,60),RS(3,3,48),TS(3,48),TP(4,3),ITF(60),    
     & NA, NSS, VOL,  DLIMIT(3),  PI, G(3,3), NLP, IERR, TITLE(18)
     &,SIGX(128000),CORR(5000)     
      COMMON/TENSOR/ CR,DR,ER(21,6),FR(28,6),CRF(3,3,3),
     & DRF(3,3,3,3),ERF(3,3,3,3,3)
      COMMON/METRIC/ AG,BG,CG,WINKEL,ARE(3),WINRE(3)
C                                                                           
      EQUIVALENCE    (CV,X), (CV(141),Y), (CV(281),INTER), (CV(301),BDF)    
     & ,(IBDF(1),BDF),(CV(317),EXTRAN), (SIGX,IPARA), (SIGX(1001),KRPAR)    
     & ,(SIGX(1061),IPNT,JPNT),(SIGX(1421),PAR),(SIGX(1481),KPAR), (SIGX    
     & (1541),PARNAM,TEMP,TEMP1,TEMP2,TENP, TEOP,ZT, TEMP4, CNSTR,    
     & ICNSTR,CNSDGN,CV,KNSTR)                                              
C                                                                           
      DATA ITAPE,NCNT/10,1,6,16,31/                                         
C                                                                           
C                                                                           
C---- STATEMENT FUNCTIONS                                                   
C                                                                           
      I2(K) = MOD(K/256**2,256)                                             
      I3(K) = MOD(K/256,256)                                                
      I4(K) = MOD(K,256)
C      IDIM = 128000	                                                    
C                                                                           
C---- DEFINE TWO INDEXING FUNCTIONS FOR TRIANGULAR MATRICES                 
C                                                                           
      INDEXL(I,J)=( ( MAX0(I,J)-1 )* MAX0(I,J) )/2 + MIN0(I,J) + NMXL       
      INDEXU(I,J)=(MRANK+1)*(MIN0(I,J)-1)-((MIN0(I,J)-1)*MIN0(I,J))/2       
     &            +MAX0(I,J)+1-MIN0(I,J)+MDEGN                              
C                                                                           
    2 FORMAT(18A4)                                                          
    4 FORMAT(3I1,3F8.0)                                                     
    6 FORMAT(6E12.6)                                                        
    8 FORMAT(7H1BONVIB/1X,6(1H-)//)                                         
   81 FORMAT(13H1TITLE OF RUN/1X,12(1H-)//1X,18A4)                          
   82 FORMAT(/18H0TIME OF EXECUTION/1X,17(1H-)//)                           
  121 FORMAT(/16H0CORE ALLOCATION /1X,15(1H-)//                             
     & 42H DIMENSION OF ARRAY SIGX (COMMON ABLE ) =  ,2X,I5,15X,            
     & 19HCURRENT MAXIMUM IS ,I6/)                                          
  124 FORMAT(13H0ERRORMESSAGE/1X,12(1H-)// 22H CORE IS INSUFFICIENT.,       
     &70H0YOU HAVE TO CHANGE THE DIMENSION OF SIGX IN COMMON  ABLE  TO A    
     &T LEAST ,I10,1H./17H0CURRENT VALUE IS,I9,1H.//16H0JOB TERMINATED.)    
  240 FORMAT(///////, 43X,29H LABEL  OF THE UNIT MAIN  =     ,3A4,          
     &10(/),35X,10HSUBSTANCE=,15X,12HTEMPERATURE=,15X,8HSPECIAL= )          
 1111 FORMAT(20X,42H P R O M E T H E U S  HAS BEEN CALLED FOR               
     & ,47HCALCULATION OF BONDS AND/OR THERMAL ELLIPSOIDS //                
     &40X,50H            WRITTEN  BY L.W. FINGER                  //        
     &40X,50H                                                     ///       
     &40X,51H THIS PROGRAM CALCULATES BONDS AND ANGLES AND THEIR  /         
     &40X,51H STANDARD DEVIATIONS.    IN ADDITION THE MEANSQUARE  /         
     &40X,51H AMPLITUDES,  THEIR  ERRORS  AND   THE    PRINCIPAL  /         
     &40X,51H AXES  OF  THERMAL   ELLIPSOIDS  CAN  BE  COMPUTED.  /         
     &40X,51H INFORMATION   ABOUT  THIS  PROGRAM  IS   AVAILABLE  /         
     &40X, 6H FROM ,11A4)                                                   
C                                                                           
C                                                                           
C----*******************************************************************    
C                                                                           
                          MAXSIG = 128000                                     
C                                                                           
C-----******************************************************************    
C                                                                           
      READ(IN,2)TITLE                                                       
C                                                                           
C---- READ MAINUNIT OF PROGRAM REFINE                                       
C                                                                           
  100 READ(ITAPE)KCYL,KLAB,AS,CS,G,VOL,NSS,NA,NSCAT,NLP,MRANK,NPARA,NSF,    
     &TENP,ICENT,(((RS(I,J,K),I=1,3),J=1,3),K=1,NSS),((TS(I,J),I=1,3),J=    
     &1,NSS),((TP(I,J),J=1,3),I=1,NLP),(IPARA(I),I=1,NPARA),((TEMP1(I,J)    
     &,I=1,9),J=1,NSCAT),(TEMP(I),I=1,NSF),TEOP,(JSF(I),I=1,NSF),NRSF       
     &,WINKEL,AG,BG,CG                                                      
     &,NPAR,(PAR(I),(PARNAM(J,I),J=1,2),KPAR(I),KRPAR(I),I=1,NPAR),FISGN    
     &,IFSGN,TEMP4,ITEMP4,BDF,GAMMAP,WVLNTH,INTPO,(JUNK6,I=1,4)             
     &,(INTER(I),I=1,NA),((X(J,I),I=1,INTPO),J=1,35),((Y(J,I),I=1,INTPO)    
     &,J=1,35),IUB,EXTRAN,IZACH,IREJEC,UNTEN,OBEN,MODELL,(ALF(I),I=1,NA) 
     &,(IALFA(I),I=1,NA)                                                 
      READ(ITAPE)((TEMP2(I,J),I=1,2),J=1,NSCAT),IEXT,EXTIN,((ATOM(J,I),J
     &=1,2),ITF(I),(ZT(J,I),J=1,5),(XYZ(J,I),J=1,3),(BETA(J,I),J=1,6),  
     &(CR(J,I),J=1,10),(DR(J,I),J=1,15),I=1,NA),MDEGN                   
     &,(CNSTR(I),I=1,MDEGN)
      DO 102 I=1,3
      ARE(I) = AS(I)
  102 WINRE(I) = CS(I)                                            
C                                                                       
C---- WRITE COVER                                                       
C                                                                       
      IF(MARK.NE.0) GOTO 1050                                           
      CALL SPACE(1,13,31,IOUT,MAXLIN)                                   
      WRITE(IOUT,1111)LINE                                              
      WRITE(IOUT,240)KLAB                                               
      CALL SPACE(2,13,31,IOUT,MAXLIN)                                   
      GOTO 101                                                          
C                                                                       
 1050 WRITE(IOUT,8)                                                     
C                                                                       
  101 WRITE(IOUT,81)TITLE                                               
      WRITE(IOUT,82)                                                    
      CALL ANFANG(0,iout)                                                    
C                                                                       
      IF(NPAR.EQ.0)GO TO 110                                            
C                                                                       
C---- FOR RIGID BODY REFINEMENTS APPLY LIBRATION CORRECTIONS TO XYZS.   
C                                                                       
      CALL LIBCOR(XYZ,PAR,KPAR,NPAR,G,NA,MDEGN,KNSTR)                   
      J=NPAR                                                            
      NPAR=0                                                            
      DO 105 I=1,J                                                      
      IF(KRPAR(I).GT.0)NPAR=NPAR+1                                      
  105 CONTINUE                                                          
C---- READ AND DISCARD STRUCTURE FACTOR RECORDS                         
  110 CONTINUE                                                          
  111 READ(ITAPE)ITEMP                                                  
      IF(ITEMP)120,111,120                                              
C---- READ VARIANCE-COVARIANCE MATRIX                                   
  120 READ(ITAPE)MDEGN,(CNSTR(I),I=1,MDEGN)                             
C                                                                       
C---- ENOUGH CORE TO RUN JOB?                                           
C                                                                       
      NMXL=(MRANK*(MRANK+1))/2+MDEGN                                    
      ISOUT=NMXL+3720                                                   
      WRITE(IOUT,121)ISOUT,MAXSIG                                       
      IF(MAXSIG-3720-NMXL)122,123,123                                   
  122 WRITE(IOUT,124)                                                   
      CALL EXIT                                                         
C                                                                       
C---- YES.                                                              
C                                                                       
  123 J=MDEGN+1                                                         
      READ(ITAPE)(CV(I),I=J,NMXL)                                       
      READ(IN,4)MBODAN,MELVIB,MTENSA,DLIMIT                           
      IF(FISGN.GT.0)GO TO 142                                           
C---- GET NEGATIVE OF POSITIONS                                         
      DO 141 I=1,NA                                                     
      DO 141 J=1,3                                                      
  141 XYZ(J,I)=-XYZ(J,I)                                                
C---- CALCULATE NUMBER OF MATRIX ROWS FOR EXTINCTION CORR.              
  142 IEXT=NPAR                                                         
      DO 143 I=1,7                                                      
      IF(ITEMP4(I).GT.0)IEXT=IEXT+1                                     
  143 CONTINUE                                                          
C                                                                       
C---- CALCULATE ROWS OF ALPHA REFINEMENT                                
C                                                                       
      IF(MODELL.LT.3) GOTO 144                                          
      DO 145 I=1,NA                                                     
      IF(IALFA(I).GT.0) IEXT=IEXT+1                                     
  145 CONTINUE                                                          
C---- CONVERT PARAMETER SELECTION INFORMATION TO MATRIX ROW POINTER     
  144 IJ=NRSF+IEXT                                                      
      DO 160 I=1,NPARA                                                  
      IF(IPARA(I))150,160,150                                           
  150 IJ=IJ+1                                                           
      IPARA(I)=IJ                                                       
  160 CONTINUE                                                          
C---- SKIP NEXT SECTION IF NO ELLIPSOIDS CALCULATED                     
      IF(MELVIB.EQ.0) GO TO 239                                         
C---- CRACK OUT ANISOTROPIC TEMPERATURE FACTOR VAR-COV MATRICES         
      IL=5                                                              
      DO 232 I=1,NA                                                     
      IF(ITF(I))170,230,170                                             
C---- CLEAR ARRAY                                                       
  170 DO 180 J=1,6                                                      
      IPNT(J,I)=0                                                       
      DO 180 K=1,6                                                      
  180 SIGB(K,J,I)=0.0                                                   
      IK=IL-1                                                           
      DO 220 J=1,6                                                      
C---- COMPUTE ROW FOR THIS PARAMETER                                    
      IK=IK+1                                                           
      JA=IPARA(IK)                                                      
      IPNT(J,I)=JA                                                      
      PI = ACOS(-1.0)                                                   
      IF(JA)190,220,190                                                 
  190 IJ=IK-1                                                           
      KA=INDEXU(JA,JA)                                                  
      DO 210 K=J,6                                                      
      IJ=IJ+1                                                           
      JB=IPARA(IJ)                                                      
      IF(JB)200,210,200                                                 
  200 KB=KA+JB-JA                                                       
      SIGB(K,J,I)=CV(KB)                                                
      SIGB(J,K,I)=CV(KB)                                                
  210 CONTINUE                                                          
  220 CONTINUE                                                          
  230 IJ=ITF(I)                                                         
  232 IL=IL+NCNT(IJ+1)+4                                                
C---- MODIFY SIGB IF NECESSARY                                          
      IF(MDEGN.LE.0) GO TO 239                                          
C---- LOOP THROUGH DEPENDENCIES                                         
      DO 235 I=1,MDEGN,2                                                
      IPAR=I4(ICNSTR(I))                                                
      IF(IPAR.LT.5.OR.IPAR.GT.10)GO TO 235                              
C---- DEP. PARAMETER IS AN ANISO. COEF.                                 
      IPAR=IPAR-4                                                       
      IAT=I3(ICNSTR(I))                                                 
      J=I2(ICNSTR(I))                                                   
C---- IS THE INDEPENDENT PARAMETER A MEMBER OF THE ANISOTROPIC          
C     TEMPERATURE FACTOR MATRIX FOR THE SAME ATOM? IF SO, THE           
C     COVARIANCE OF THE INDEPENDENT PARAMETER WITH THE DEPENDENT        
C     ONE MUST BE COMPUTED.                                             
      DO 1010 JJ=1,6                                                    
      JI=IPNT(JJ,IAT)                                                   
      IF(JI.EQ.0)GO TO 1010                                             
      IK=INDEXU(J,JI)                                                   
      SIGB(JJ,IPAR,IAT)=SIGB(JJ,IPAR,IAT)+CNSTR(I+1)*CV(IK)             
      SIGB(IPAR,JJ,IAT)=SIGB(JJ,IPAR,IAT)                               
 1010 CONTINUE                                                          
C---- COMPUTE COVARIANCES OF DEPENDENT PARAMETERS WITH EACH OTHER.      
 1011 DO 1500 KK=1,MDEGN,2                                              
      IF(I3(ICNSTR(KK)).NE.IAT.OR.I4(ICNSTR(KK)).LT.5.OR.I4(ICNSTR(KK)) 
     & .GT.10)GO TO 1500                                                
      IJPAR=I4(ICNSTR(KK))-4                                            
      L=I2(ICNSTR(KK))                                                  
      IK=INDEXU(J,L)                                                    
 1200 SIGB(IPAR,IJPAR,IAT)=SIGB(IPAR,IJPAR,IAT)+CNSTR(I+1)*CNSTR(KK+1)* 
     &                     CV(IK)                                       
      IF(IPAR.NE.IJPAR)SIGB(IJPAR,IPAR,IAT)=SIGB(IPAR,IJPAR,IAT)        
 1500 CONTINUE                                                          
  235 CONTINUE                                                          
  239 LXYZ=0                                                            
C---- SKIP NEXT SECTION IF NO BOND FUNCTIONS                            
      IF(MBODAN.EQ.0) GO TO 790                                         
C---- START UNPACKING POSITION MATRIX                                   
      IPAR=NRSF+IEXT+1                                                  
C---- KA WILL BE THE INDEX OF THE NEXT LOCATION IN THE UNPACKED ARRAY   
C     INCR IS THE STEP NEEDED TO ADVANCE KA AN ENTIRE ROW               
      KA=INDEXU(IPAR,IPAR)                                              
      INCR=MRANK+1-IPAR                                                 
C---- KB WILL INDEX THE NEXT LOCATION IN THE PACKED MATRIX              
      KB=NMXL+1                                                         
      LXYZ=(3*NA*(3*NA+1))/2                                            
      I=LXYZ+NMXL                                                       
      ISOUT=I+3720                                                      
C---- REQUEST SUFFICIENT CORE FOR XYZ VAR-COV, THEN CLEAR               
      IF(MAXSIG-3720-I)352,353,353                                      
  352 WRITE(IOUT,124)ISOUT,MAXSIG                                       
      CALL EXIT                                                         
C                                                                       
  353 DO 2401 J=KB,I                                                    
 2401 CV(J)=0.0                                                         
      NP=1                                                              
      KB=KB-1                                                           
      IPAR=0                                                            
C---- LOOP THROUGH ATOMS                                                
      DO 540 I=1,NA                                                     
C---- ADVANCE PARAMETER COUNTER FOR OCCUPANCY                           
      IPAR=IPAR+1                                                       
      J=ITF(I)                                                          
      NBI=NCNT(J+1)                                                     
      IF(IPARA(IPAR).EQ.0)GO TO 250                                     
C---- STEP OVER VAR-COV ROW FOR OCCUPANCY                               
      KA=KA+INCR                                                        
      INCR=INCR-1                                                       
C---- LOOP THROUGH POSITIONS                                            
  250 DO 480 J=1,3                                                      
      KB=KB+NP                                                          
      IF(FISGN.LT.0.0)XYZ(J,I)=-XYZ(J,I)                                
      IPAR=IPAR+1                                                       
      JPNT(NP)=IPARA(IPAR)                                              
C---- TEST IF POSITION REFINED                                          
      IF(IPARA(IPAR))260,480,260                                        
C---- COPY VARIANCE OF POSITION                                         
  260 CV(KB)=CV(KA)                                                     
      INCR=INCR-1                                                       
      KA=KA+1                                                           
      KKPAR=IPAR                                                        
      JP=NP                                                             
      KC=KB                                                             
      IF(J-3)280,340,280                                                
  280 ITEMP=J+1                                                         
C---- LOOP THROUGH OTHER POSITIONS FOR ATOM I                           
      DO 300 K=ITEMP,3                                                  
      KKPAR=KKPAR+1                                                     
      KC=KC+JP                                                          
      JP=JP+1                                                           
      IF(IPARA(KKPAR))290,300,290                                       
C---- COPY COVARIANCE ELEMENT                                           
  290 CV(KC)=CV(KA)                                                     
      KA=KA+1                                                           
  300 CONTINUE                                                          
C---- LOOP THROUGH TEMPERATURE FACTOR(S) FOR ATOM I                     
  340 DO 360 K=1,NBI                                                    
      KKPAR=KKPAR+1                                                     
      IF(IPARA(KKPAR))350,360,350                                       
  350 KA=KA+1                                                           
  360 CONTINUE                                                          
      IF(I-NA)370,480,370                                               
  370 ITEMP=I+1                                                         
C---- LOOP THROUGH THE OTHER ATOMS                                      
      DO 470 K=ITEMP,NA                                                 
      KKPAR=KKPAR+1                                                     
      IF(IPARA(KKPAR))380,390,380                                       
C---- OCCUP REFINED                                                     
  380 KA=KA+1                                                           
C---- LOOP THROUGH POSITIONS FOR ATOM K                                 
  390 DO 410 L=1,3                                                      
      KC=KC+JP                                                          
      JP=JP+1                                                           
      KKPAR=KKPAR+1                                                     
      IF(IPARA(KKPAR))400,410,400                                       
C---- COPY COVARIANCE ELEMENT                                           
  400 CV(KC)=CV(KA)                                                     
      KA=KA+1                                                           
  410 CONTINUE                                                          
      L=ITF(K)                                                          
      NB=NCNT(L+1)                                                      
C---- LOOP THROUGH TEMPERATURE FACTOR(S)                                
  440 DO 460 L=1,NB                                                     
      KKPAR=KKPAR+1                                                     
      IF(IPARA(KKPAR))450,460,450                                       
C---- ADVANCE COUNTER                                                   
  450 KA=KA+1                                                           
  460 CONTINUE                                                          
  470 CONTINUE                                                          
  480 NP=NP+1                                                           
C---- LOOP THROUGH VAR-COV ROWS FOR ATOM I TEMP FACTORS                 
      DO 530 J=1,NBI                                                    
      IPAR=IPAR+1                                                       
      IF(IPARA(IPAR))520,530,520                                        
C---- ADVANCE INDEX AND REDUCE INCREMENT                                
  520 KA=KA+INCR                                                        
      INCR=INCR-1                                                       
  530 CONTINUE                                                          
  540 CONTINUE                                                          
C---- ATOM VARIANCE-COVARIANCE MATRIX IS EXTRACTED, NOW MODIFY FOR      
C     CONSTRAINTS                                                       
 9999 FORMAT(9E14.6)                                                    
      IF(MDEGN.EQ.0) GO TO 570                                          
      NB=3*NA                                                           
      DO 565 I=1,MDEGN,2                                                
      IPAR=I4(ICNSTR(I))                                                
      IF(IPAR.LT.2.OR.IPAR.GT.4)GO TO 565                               
      IAT=I3(ICNSTR(I))                                                 
      J=I2(ICNSTR(I))                                                   
      IPAR=3*(IAT-1)+IPAR                                               
      DO 560 KK=1,NB                                                     
      K = KK
      IF(JPNT(K).EQ.0)GO TO 560                                         
C---- COMPUTE COVARIANCE OF POSITION WITH A REFINED POSITION.           
      KB=INDEXL(K,IPAR)                                                 
      KA=INDEXU(J,JPNT(K))                                              
      CV(KB)=CV(KB)+CV(KA)*CNSDGN(I+1)                                  
  560 CONTINUE                                                          
C     WRITE(IOUT,9999)(CV(K+NMXL),K=1,9)                                
C---- COMPUTE COVARIANCES OF UNREFINED POSITIONS WITH EACH OTHER.       
      DO 563 KK=1,MDEGN,2                                               
      IJPAR=I4(ICNSTR(KK))                                              
      IF(IJPAR.LT.2.OR.IJPAR.GT.4)GO TO 563                             
      JAT=I3(ICNSTR(KK))                                                
      L=I2(ICNSTR(KK))                                                  
      IJPAR=3*(JAT-1)+IJPAR                                             
      KA=INDEXU(J,L)                                                    
      KB=INDEXL(IPAR,IJPAR)                                             
      CV(KB)=CV(KB)+CNSDGN(I+1)*CNSDGN(KK+1)*CV(KA)                     
  563 CONTINUE                                                          
  565 CONTINUE                                                          
C---- MATRIX NOW TRANSFORMED, COPY TO LOW CORE                          
  570 KB=NMXL+1                                                         
      DO 580 I=1,LXYZ                                                   
      SIGX(I)=CV(KB)                                                    
  580 KB=KB+1                                                           
C---- CONTRACT MATRIX TO THE MINIMUM REQUIRED                           
C                                                                       
C---- READ AND CONVERT CELL ERROR MATRIX                                
 790  READ(IN,6)((SIGA(J,I),J=I,6),I=1,6)                               
C---- GET REAL SPACE CELL QUANTITIES                                    
      IERR=1                                                            
      DO 800 I=1,3                                                      
  800 A(I)=SQRT(G(I,I))                                                 
      ANG(1)=G(2,3)/(A(2)*A(3))                                         
      ANG(2)=G(1,3)/(A(1)*A(3))                                         
      ANG(3)=G(1,2)/(A(1)*A(2))                                         
      DO 810 I=1,3                                                      
  810 S(I)=SQRT(1.0-ANG(I)**2)                                          
      DO 840 I=1,6                                                      
      DO 840 J=I,6                                                      
      IF(I.LE.3.AND.J.LE.3) GO TO 840                                   
      IF(I-3)820,820,830                                                
  820 SIGA(J,I)=SIGA(J,I)*S(I-3)                                        
      GO TO 840                                                         
  830 SIGA(J,I)=SIGA(J,I)*S(I-3)*S(J-3)                                 
  840 SIGA(I,J)=SIGA(J,I)                                               
  850 IF(ICENT)880,860,880                                              
C---- PRINT VARIANCE-COVARIANCE MATRIX OF CELL
  860 WRITE(IOUT,9998) ((SIGA(J,I),J=1,6),I=1,6)
 9998 FORMAT(/,35H0VARIANCE-COVARIANCE MATRIX OF CELL,/,1X,34(1H-),//,
     &6(6E14.6,/))
C---- EXPAND SYMMETRY OPERATORS FOR A CENTRIC CRYSTAL                   
      ITEMP=NSS                                                         
      DO 870 I=1,NSS                                                    
      ITEMP=ITEMP+1                                                     
      DO 870 J=1,3                                                      
      TS(J,ITEMP)=-TS(J,I)                                              
      DO 870 K=1,3                                                      
  870 RS(K,J,ITEMP)=-RS(K,J,I)                                          
      NSS=ITEMP                                                         
  880 IF(MBODAN)890,900,890                                             
C---- COMPUTE BOND FUNCTIONS                                            
  890 CALL BODAN(MBODAN,IOUT)                                           
  900 IF(MELVIB)910,920,910                                             
C---- COMPUTE THERMAL FUNCTIONS                                         
  910 CALL ELVIB(IOUT)                                                  
C---- ANALYSE THERMAL MOTION TENSORS
  920 IF(MTENSA)930,940,930
  930 CALL TENSAL(IOUT,NA)
C---- END OF JOB                                                        
  940 RETURN                                                            
                                                                        
      END                                                               
C                                                                       
C                                                                       
                  SUBROUTINE BODAN(IBOD,IOUT)                           
                                                                        
                                                                        
C         BOND DISTANCE AND ANGLE SUBPROGRAM WITH ERRORS                
                                                                        
C                          L.W. FINGER                                  
C                                                                       
      COMMON/ABLE/XYZ(3,60),BETA(6,60),A(3),ANG(3),AS(3),CS(3),ATOM(2,  
     & 60),SIGA(6,6),SIGB(6,6,60),RS(3,3,48),TS(3,48),TP(4,3),ITF(60),  
     & NA,NSS,VOL,DLIMIT(3),PI,G(3,3),NLP,IERR,TITLE(18),
     & SIGX(128000),CORR(5000)     
      COMMON/BAKER/XYZP(3,1000),IATP(1000)                              
                                                                        
      DIMENSION XP(3),XYP(3),XXP(3),SIGZ(6,6),DLDX(6),DSDX(9,2),YP(3),  
     & DLDA(6,2),SIGY(9,9),DNDA(6),DNDX(9),DTHDA(6),DTHDX(9),XZP(3),    
     & X(3,3),NXX(6),DXX(3),DXY(3),DXZ(3),DBD(3),SIGAA(6,6,1),DLDAA(6,1)                        
                                                                        
      EQUIVALENCE (DLDX,DSDX),(NXX(1),NXL),(NXX(2),NYL),(NXX(3),NZL)    
     &,(SIGAA,SIGA),(DLDAA,DLDA)
                                                                        
 1002 FORMAT(46X,2A3,4F10.5,10X,I8)                                     
 8999 FORMAT(31H THE ATOM ARRAY IS OVERFLOWING. ///15X,54HREDUCE THE LIM
     &ITS OF THE BOND DISTANCES AND TRY AGAIN. )                        
 9000 FORMAT(1H1,18A4//25H BOND DISTANCES BETWEEN  ,F5.2,5H AND ,F5.2,10
     &H ANGSTROMS//10X,6HATOM 1,7X1HX,9X,1HY,9X,1HZ,2X,6HATOM 2,7X,1HX,9
     &X,1HY,9X,1HZ,4X,8HDISTANCE,3X5HERROR,4X,6HNUMBER//10X2A3,3F10.5,63
     &X,1H1)                                                            
 9001 FORMAT(10X2A3,3F10.5,56XI8)                                       
 9002 FORMAT(1H1,18A4//69H BOND ANGLES AND ERRORS FOR ATOMS WITH INTERAT
     &OMIC DISTANCES BETWEEN F6.2,5H AND F6.2,11H ANGSTROMS./102H0THE  A
     &NGLE IS ATOM 1-CENTRAL ATOM-ATOM 2.  THE NUMBER POINTS TO AN ENTRY
     & IN THE BOND DISTANCE OUTPUT.)                                    
 9003 FORMAT(15X,2A3,I6,12X,2A3,I6,F17.2,F10.2)                         
 9005 FORMAT(46X2A3,5F10.5,I8)                                          
 9010 FORMAT(///10X,12HCENTRAL ATOM,5X2A3,3F10.5//15X,6HATOM 1,2X,6HNUMB
     &ER,10X6HATOM 2,2X,6HNUMBER,10X,5HANGLE,5X,5HERROR/)               
                                                                        
      WRITE(IOUT,9000)TITLE,DLIMIT(1),DLIMIT(3),(ATOM(J,1),J=1,2),(XYZ(J
     &,1),J=1,3)                                                        
      NAHELP=NA                                                         
C---- SET UP PRELIMINARY INFORMATION FOR SEARCHING FOR ATOMS            
      PROD=A(1)*A(2)*A(3)                                               
      DO 61 J=1,3                                                       
      XYZP(J,1)=XYZ(J,1)                                                
      DBD(J)=PROD/A(J)                                                  
   61 DXX(J)=AS(J)*DLIMIT(3)                                            
C---- INITIALIZE AND FINISH STORING ATOM 1 IN THE ASYMMETRIC UNIT       
C     INTO THE GENERATED ATOM LIST                                      
C     THE ELEMENTS OF THE ARRAY IATP ARE 64*ATOM NUMBER + SYMMETRY NO   
      IATP(1)=65                                                        
      NAA=3*NA                                                          
      DO 37 I=1,9                                                       
      DO 37 J=1,2                                                       
   37 DSDX(I,J)=0                                                       
      NATOM=1                                                           
      SN=DLIMIT(1)**2                                                   
      SX=DLIMIT(3)**2                                                   
C---- LOOP THROUGH THE ATOMS IN THE ASYMMETRIC UNIT                     
      DO 1 IIHELP=1,NA                                                  
C---- SET UP SEARCH ARRAYS FOR THIS ATOM                                
      II=IIHELP                                                         
      DO 62 J=1,3                                                       
      DXY(J)=XYZ(J,II)+DXX(J)                                           
   62 DXZ(J)=XYZ(J,II)-DXX(J)                                           
C---- EXTRACT SELF VARIANCE-COVARIANCE MATRIX                           
      IJ=3*II-2                                                         
      JI=((IJ+1)*IJ)/2                                                  
      DO 3 N=1,3                                                        
      L=JI                                                              
      SIGZ(N,1)=SIGX(L)                                                 
      IF(N.EQ.1)GO TO 3                                                 
      DO 2 J=2,N                                                        
      L=L+1                                                             
    2 SIGZ(N,J)=SIGX(L)                                                 
    3 JI=JI+IJ+N-1                                                      
C---- IF FIRST ATOM, JUMP AROUND NEXT CODE                              
      IF(II.EQ.1)GO TO 67                                               
C---- SEARCH FOR PREVIOUS GENERATION OF THIS ATOM                       
   68 DO 115 I=1,NATOM                                                  
      IF(IATP(I)/64.EQ.II)GO TO 116                                     
      GO TO 115                                                         
  116 DO 117 J=1,3                                                      
      IF(ABS (XYZP(J,I)-XYZ(J,II))-0.00002)117,117,115                  
  117 CONTINUE                                                          
C---- ABOVE LOOP COMPLETED IF NOT NEW POSITION                          
      IJ=I                                                              
      GO TO 118                                                         
  115 CONTINUE                                                          
C---- ADD ATOM TO LIST                                                  
      IJ=NATOM+1                                                        
      NATOM=IJ                                                          
      IATP(NATOM)=II*64+1                                               
      DO 119 J=1,3                                                      
  119 XYZP(J,NATOM)=XYZ(J,II)                                           
C---- PRINT HEADING FOR ATOM IN ASYMMETRIC UNIT                         
  118 WRITE(IOUT,9001)(ATOM(J,II),J=1,2),(XYZ(J,II),J=1,3),IJ           
C---- LOOP THROUGH GENERATED LIST TO FIND DISTANCES                     
   70 DO 71 J=1,NATOM                                                   
      S=0                                                               
      DO 72 K=1,3                                                       
   72 XP(K)=XYZ(K,II)-XYZP(K,J)                                         
      DO 73 N=1,3                                                       
      XZP(N)=0                                                          
      DO 74 K=1,3                                                       
   74 XZP(N)=XZP(N)+XP(K)*G(K,N)                                        
   73 S=S+XZP(N)*XP(N)                                                  
      IF((S-SX)*(SN-S))71,71,75                                         
   75 S=SQRT (S)                                                        
      IF(IERR.EQ.0)GO TO 1001                                           
C---- COMPUTE DERIVATIVES                                               
      L=IATP(J)-64*(IATP(J)/64)                                         
      DLDA(4,1)=DBD(1)*XP(2)*XP(3)                                      
      DLDA(5,1)=DBD(2)*XP(1)*XP(3)                                      
      DLDA(6,1)=DBD(3)*XP(1)*XP(2)                                      
      DO 76 K=1,3                                                       
      DLDA(K,1)=XP(K)*XZP(K)/A(K)                                       
      DLDX(K)=XZP(K)                                                    
      DLDX(K+3)=0.                                                      
      DO 76 N=1,3                                                       
   76 DLDX(K+3)=DLDX(K+3)-XZP(N)*RS(N,K,L)                              
C---- FILL IN SELF VARIANCE-COVARIANCE MATRIX                           
      IJ=3*(IATP(J)/64)-2                                               
      JI=((IJ+1)*IJ)/2                                                  
      DO 77 K=1,3                                                       
      L=JI                                                              
      SIGZ(K+3,4)=SIGX(L)                                               
      IF(K.EQ.1)GO TO 77                                                
      DO 78 N=2,K                                                       
      L=L+1                                                             
   78 SIGZ(K+3,N+3)=SIGX(L)                                             
   77 JI=JI+IJ+K-1                                                      
C---- GET CROSS TERMS FOR VARIANCE-COVARIANCE MATRIX                    
      CALL CORR1(SIGZ,6,II,IATP(J)/64,SIGX,1,4,NAHELP)                  
C---- COMPUTE VARIANCE OF DISTANCE                                      
      CALL VARIAN(DS,DLDX,DLDAA,SIGZ,SIGAA,1)                             
C     WRITE(IOUT,9997) DS,DLDX,DLDAA,SIGZ,SIGAA
C9997 FORMAT(1X,E14.6,/,1X,6E14.6,/,2(6E14.6,/),/,6(6E14.6,/),
C    &/,6(6E14.6,/))
      DS=DS/S                                                           
      IJ=IATP(J)/64                                                     
      WRITE(IOUT,9005)(ATOM(K,IJ),K=1,2),(XYZP(K,J),K=1,3),S,DS,J       
      GO TO 71                                                          
C---- PRINT DISTANCE FOR CALCULATION WITHOUT ERRORS                     
 1001 IJ=IATP(J)/64                                                     
      WRITE(IOUT,1002)(ATOM(K,IJ),K=1,2),(XYZP(K,J),K=1,3),S,J          
   71 CONTINUE                                                          
                                                                        
                                                                        
C---- LOOP THROUGH THE ASYMMETRIC UNIT TO GENERATE ATOMS                
   67 DO 1 JJHELP= 1,NA                                                 
C---- LOOP THROUGH THE SYMMETRY OPERATIONS                              
      JJ=JJHELP                                                         
      DO 12 J=1,NSS                                                     
C---- GENERATE TRANSFORMED POSITION                                     
      DO 13 L=1,3                                                       
      XP(L)=TS(L,J)                                                     
      DO 13 M=1,3                                                       
   13 XP(L)=XP(L)+RS(L,M,J)*XYZ(M,JJ)                                   
C---- LOOP THROUGH LATTICE POINTS                                       
      DO 12 M=1,NLP                                                     
      DO 14 L=1,3                                                       
   14 XXP(L)=XP(L)+TP(M,L)                                              
C---- COMPUTE TRANSLATION LIMITS                                        
      DO 63 I=1,3                                                       
      NXX(I+3)=INT(3.0+DXY(I)-XXP(I))                                   
   63 NXX(I)=3-IFIX (XXP(I)-DXZ(I))                                     
      NXU=NXX(4)                                                        
      NYU=NXX(5)                                                        
      NZU=NXX(6)                                                        
      DO 12 IA=NXL,NXU                                                  
      XYP(1)=XXP(1)+FLOAT (IA-3)                                        
      DO 12 IB=NYL,NYU                                                  
      XYP(2)=XXP(2)+FLOAT (IB-3)                                        
      DO 12 IC=NZL,NZU                                                  
      XYP(3)=XXP(3)+FLOAT (IC-3)                                        
C---- CHECK FOR ATOM DUPLICATION                                        
      DO 113 I=1,NATOM                                                  
      IF(IATP(I)/64.EQ.JJ)GO TO 24                                      
      GO TO 113                                                         
   24 DO 114 L=1,3                                                      
      IF(ABS (XYP(L)-XYZP(L,I))-0.00002)114,114,113                     
  114 CONTINUE                                                          
C---- THROUGH HERE IF THIS ATOM HAS ALREADY BEEN GENERATED              
      GO TO 12                                                          
  113 CONTINUE                                                          
C---- COMPUTE DISTANCE                                                  
      S=0                                                               
      DO 15 K=1,3                                                       
   15 YP(K)=XYZ(K,II)-XYP(K)                                            
      DO 17 N=1,3                                                       
      XZP(N)=0                                                          
      DO 16 K=1,3                                                       
   16 XZP(N)=XZP(N)+YP(K)*G(K,N)                                        
   17 S=S+XZP(N)*YP(N)                                                  
C---- CHECK IF DISTANCE WITHIN LIMITS                                   
      IF((S-SX)*(SN-S))12,12,18                                         
C---- DISTANCE OK, CHECK IF ROOM IN GENERATED LIST                      
   18 IF(NATOM-1000)19,20,20                                            
C---- ARRAYS OVERFLOWING,  PRINT MESSAGE AND TERMINATE                  
   20 WRITE(IOUT,8999)                                                  
      CALL EXIT                                                         
C                                                                       
C---- ADD GENERATED ATOM TO LIST AND COMPUTE DERIVATIVES                
   19 NATOM=NATOM+1                                                     
      IATP(NATOM)=JJ*64+J                                               
      DLDA(4,1)=DBD(1)*YP(2)*YP(3)                                      
      DLDA(5,1)=DBD(2)*YP(1)*YP(3)                                      
      DLDA(6,1)=DBD(3)*YP(1)*YP(2)                                      
      DO 21 K=1,3                                                       
      XYZP(K,NATOM)=XYP(K)                                              
      DLDX(K)=XZP(K)                                                    
      DLDX(K+3)=0                                                       
      DLDA(K,1)=YP(K)*XZP(K)/A(K)                                       
      DO 21 N=1,3                                                       
   21 DLDX(K+3)=DLDX(K+3)-XZP(N)*RS(N,K,J)                              
      IJ=3*JJ-2                                                         
      IF(IERR.EQ.0)GO TO 1004                                           
 1003 JI=((IJ+1)*IJ)/2                                                  
      DO 23 K=1,3                                                       
      L=JI                                                              
      SIGZ(K+3,4)=SIGX(L)                                               
      IF(K.EQ.1)GO TO 23                                                
      DO 22 N=2,K                                                       
      L=L+1                                                             
   22 SIGZ(K+3,N+3)=SIGX(L)                                             
   23 JI=JI+IJ+K-1                                                      
C---- GENERATE CROSS MATRIX ELEMENTS                                    
      CALL CORR1(SIGZ,6,II,JJ,SIGX,1,4,NAHELP)                          
      S=SQRT (S)                                                        
      CALL VARIAN(DS,DLDX,DLDAA,SIGZ,SIGAA,1)                             
      DS=DS/S                                                           
      WRITE(IOUT,9005)(ATOM(N,JJ),N=1,2),XYP,S,DS,NATOM                 
      GO TO 12                                                          
 1004 S=SQRT (S)                                                        
C---- PRINT DISTANCE FOR NO ERRORS CALCULATED                           
      WRITE(IOUT,1002)(ATOM(N,JJ),N=1,2),XYP,S,NATOM                    
   12 CONTINUE                                                          
    1 CONTINUE                                                          
C---- END OF BOND DISTANCE CALCULATION, RETURN IF NO ANGLES             
      IF(IBOD.EQ.1) RETURN                                              
      WRITE(IOUT,9002)TITLE,(DLIMIT(I),I=1,2)                           
C---- SET UPPER LIMIT OF DISTANCE TO CONTACT DISTANCE                   
      SX=DLIMIT(2)**2                                                   
C---- LOOP THROUGH ASYMMETRIC UNIT                                      
      DO 28 I33K=1,NA                                                   
      I=I33K                                                            
      IJ=3*I-2                                                          
      JI=((IJ+1)*IJ)/2                                                  
      WRITE(IOUT,9010)(ATOM(K,I),K=1,2),(XYZ(K,I),K=1,3)                
      IF(IERR.EQ.0)GO TO 1006                                           
C---- GET SELF VAR-COV MATRIX                                           
      DO 124 N=1,3                                                      
      L=JI                                                              
      SIGY(N+3,4)=SIGX(L)                                               
      IF(N.EQ.1)GO TO 124                                               
      DO 25 K=2,N                                                       
      L=L+1                                                             
   25 SIGY(N+3,K+3)=SIGX(L)                                             
  124 JI=JI+IJ+N-1                                                      
C---- OUTSIDE LOOP THROUGH GENERATED ATOM LIST                          
 1006 ITEMP=NATOM-1                                                     
      DO 29 J=1,ITEMP                                                   
      DO 30 N=1,3                                                       
   30 XP(N)=XYZP(N,J)-XYZ(N,I)                                          
C---- COMPUTE AND CHECK DISTANCE                                        
      S=0                                                               
      DO 31 N=1,3                                                       
      XXP(N)=0                                                          
      DO 32 K=1,3                                                       
   32 XXP(N)=XXP(N)+XP(K)*G(K,N)                                        
   31 S=S+XXP(N)*XP(N)                                                  
      IF((S-SX)*(SN-S))29,29,33                                         
   33 S=SQRT (S)                                                        
      NAT=IATP(J)/64                                                    
      L   =IATP(J)-64*NAT                                               
      L1=L                                                              
      IF(IERR.EQ.0)GO TO 1008                                           
C---- COMPUTE DERIVATIVES                                               
      DLDA(4,1)=DBD(1)*XP(2)*XP(3)/S                                    
      DLDA(5,1)=DBD(2)*XP(1)*XP(3)/S                                    
      DLDA(6,1)=DBD(3)*XP(1)*XP(2)/S                                    
      DO 34 N=1,3                                                       
      DSDX(N+3,1)=-XXP(N)/S                                             
      DSDX(N,1)=0.0                                                     
      DLDA(N,1)=XP(N)*XXP(N)/(S*A(N))                                   
      DO 34 K=1,3                                                       
   34 DSDX(N,1)=DSDX(N,1)+XXP(K)*RS(K,N,L)/S                            
      IJ=3*NAT-2                                                        
      JI=((IJ+1)*IJ)/2                                                  
      DO 26 K=1,3                                                       
      L=JI                                                              
      SIGY(K,1)=SIGX(L)                                                 
      IF(K.EQ.1)GO TO 26                                                
      DO 27 N=2,K                                                       
      L=L+1                                                             
   27 SIGY(K,N)=SIGX(L)                                                 
   26 JI=JI+IJ+K-1                                                      
      CALL CORR1(SIGY,9,I,NAT,SIGX,1,4,NAHELP)                          
C---- INNER LOOP THROUGH GENERATED ATOM LIST                            
 1008 IK=J+1                                                            
      DO 40 K=IK,NATOM                                                  
C---- COMPUTE AND CHECK DISTANCE                                        
      S2=0                                                              
      DO 41 N=1,3                                                       
   41 YP(N)=XYZP(N,K)-XYZ(N,I)                                          
      DO 42 N=1,3                                                       
      XYP(N)=0.                                                         
      DO 43M=1,3                                                        
   43 XYP(N)=XYP(N)+YP(M)*G(M,N)                                        
   42 S2=S2+XYP(N)*YP(N)                                                
      IF((S2-SX)*(SN-S2))40,40,44                                       
   44 S2=SQRT (S2)                                                      
      NTA=IATP(K)/64                                                    
      L2=  IATP(K)-NTA*64                                               
      L=L2                                                              
      IF(IERR.EQ.0)GO TO 1010                                           
C---- COMPUTE DERIVATIVES                                               
      DLDA(4,2)=DBD(1)*YP(2)*YP(3)/S2                                   
      DLDA(5,2)=DBD(2)*YP(1)*YP(3)/S2                                   
      DLDA(6,2)=DBD(3)*YP(1)*YP(2)/S2                                   
      DO 45 M=1,3                                                       
      DSDX(M+3,2)=-XYP(M)/S2                                            
      DSDX(M+6,2)=0.                                                    
      DLDA(M,2)=YP(M)*XYP(M)/(S2*A(M))                                  
      DO 45 N=1,3                                                       
   45 DSDX(M+6,2)=DSDX(M+6,2)+XYP(N)*RS(N,M,L)/S2                       
C---- COMPUTE COSINE THETA                                              
 1010 CTH=0.                                                            
      DO 46 M=1,3                                                       
   46 CTH=CTH+XXP(M)*YP(M)                                              
      CTH=CTH/(S*S2)                                                    
C---- CHECK FOR COSINE OUT OF RANGE                                     
      IF(ABS (CTH)-1.0)47,47,48                                         
   48 CTH=SIGN (1.0,CTH)                                                
   47 STH=SQRT (1.0-CTH**2)                                             
      CHT=ACOS (CTH)*180.0/PI                                           
      IF(IERR.NE.0) GO TO 1012                                          
C---- PRINT ANGLE FOR NO ERRORS                                         
      WRITE(IOUT,9003)(ATOM(M,NAT),M=1,2),J,(ATOM(M,NTA),M=1,2),K,CHT   
      GO TO 40                                                          
C---- IF ANGLE CLOSE TO N*PI/2, SKIP ERROR CALCULATION                  
 1012 IF(ABS (STH)-1.0E-6)64,64,65                                      
   64 DTH=0                                                             
      GO TO 66                                                          
C---- EXTRACT VAR-COV MATRICES AND COMPUTE DERIVATIVES                  
   65 IJ=3*NTA-2                                                        
      JI=((IJ+1)*IJ)/2                                                  
      DO 50 M=1,3                                                       
      L=JI                                                              
      SIGY(M+6,7)=SIGX(L)                                               
      IF(M.EQ.1)GO TO 50                                                
      DO 49 N=2,M                                                       
      L=L+1                                                             
   49 SIGY(M+6,N+6)=SIGX(L)                                             
   50 JI=JI+IJ+M-1                                                      
      CALL CORR1(SIGY,9,I,NTA,SIGX,4,7,NAHELP)                          
      CALL CORR1(SIGY,9,NAT,NTA,SIGX,1,7,NAHELP)                        
      DO 51 N=1,3                                                       
      DO 51 M=1,3                                                       
   51 X(N,M)=XP(N)*YP(M)                                                
      DO 52 N=1,3                                                       
      NM=MOD(N,3)+1                                                     
      MN=MOD(N+1,3)+1                                                   
C---- DNDA IS THE DERIVATIVE OF THE NUMERATOR WRT THE CELL PARAMS.      
      DNDA(N)=2.0*A(N)*X(N,N)+A(NM)*ANG(MN)*(X(N,NM)+X(NM,N))+          
     & A(MN)*ANG(NM)*(X(N,MN)+X(MN,N))                                  
      DNDA(N+3)=A(NM)*A(MN)*(X(NM,MN)+X(MN,NM))                         
      DNDX(N+3)=-A(N)**2*(XP(N)+YP(N))-G(N,NM)*(YP(NM)+XP(NM))-G(N,MN)* 
     & (XP(MN)+YP(MN))                                                  
      DNDX(N)=0.                                                        
      DNDX(N+6)=DNDX(N)                                                 
      DO 52 M=1,3                                                       
      DO 52 L=1,3                                                       
      DNDX(N)=DNDX(N)+G(M,L)*RS(L,N,L1)*YP(M)                           
   52 DNDX(N+6)=DNDX(N+6)+G(M,L)*RS(L,N,L2)*XP(M)                       
      SS=1.0/(S*S2)                                                     
      DO 56 N=1,6                                                       
      DTHDA(N)=SS*DNDA(N)-CTH*(DLDA(N,2)/S2+DLDA(N,1)/S)                
   56 DTHDX(N)=SS*DNDX(N)-CTH*(DSDX(N,2)/S2+DSDX(N,1)/S)                
      DTH=DTHDA(1)**2*SIGA(1,1)+DTHDX(1)**2*SIGY(1,1)                   
      DO 79 N=7,9                                                       
   79 DTHDX(N)=SS*DNDX(N)-CTH*(DSDX(N,2)/S2+DSDX(N,1)/S)                
      DO 57 N=2,6                                                       
      IJ=N-1                                                            
      DTH=DTH+DTHDA(N)**2*SIGA(N,N)                                     
      DO 57 M=1,IJ                                                      
   57 DTH=DTH+2.0*DTHDA(N)*DTHDA(M)*SIGA(N,M)                           
      DO 58 N=2,9                                                       
      IJ=N-1                                                            
      DTH=DTH+DTHDX(N)**2*SIGY(N,N)                                     
      DO 58 M=1,IJ                                                      
   58 DTH=DTH+DTHDX(N)*DTHDX(M)*SIGY(N,M)*2.0                           
      IF(DTH)59,60,60                                                   
   59 DTH=0.                                                            
   60 DTH=SQRT (DTH)/STH*180.0/PI                                       
   66 WRITE(IOUT,9003)(ATOM(M,NAT),M=1,2),J,(ATOM(M,NTA),M=1,2),K,CHT,  
     & DTH                                                              
   40 CONTINUE                                                          
   29 CONTINUE                                                          
   28 CONTINUE                                                          
      RETURN                                                            
      END                                                               
C                                                                       
              SUBROUTINE CELLD(DASDA,A,ANG,AS,CS,VOL)                   
                                                                        
C                   WRITTEN BY L.W. FINGER                              
                                                                        
C     ROUTINE TO COMPUTE THE DERIVATIVES OF THE RECIPROCAL LATTICE      
C         PARAMETERS WITH RESPECT TO THE REAL CELL PARAMETERS.          
C     DASDA(I,J)= DERIVATIVE OF RECIPROCAL PARAMETER I WITH RESPECT     
C                     TO REAL CELL PARAMETER J.                         
C                                                                       
      DIMENSION A(3),ANG(3),AS(3),CS(3),DASDA(6,6),SN(3)                
C                                                                       
      TEMP=(A(1)*A(2)*A(3)/VOL)**2                                      
      SNPROD=1.0                                                        
C---- COMPUTE SINES OF REAL SPACE ANGLES AND PRODUCT OF SINES           
      DO 100 I=1,3                                                      
      SN(I)=SQRT(1.0-ANG(I)**2)                                         
  100 SNPROD=SNPROD*SN(I)                                               
C---- CLEAR DERIVATIVE ARRAY                                            
      DO 110 I=1,6                                                      
      DO 110 J=1,6                                                      
  110 DASDA(I,J)=0.0                                                    
C---- COMPUTE DERIVATIVES                                               
      DO 120 I=1,3                                                      
      TEMP1=TEMP*AS(I)                                                  
      J=MOD(I,3)+1                                                      
      K=MOD(I+1,3)+1                                                    
      L=I+3                                                             
      DASDA(I,I)=-AS(I)/A(I)                                            
      DASDA(I,L)=-TEMP1*SN(J)*SN(K)*CS(I)-AS(I)*ANG(I)/SN(I)**2         
      DASDA(I+3,I+3)=-SN(I)/SNPROD                                      
      DASDA(I,J+3)=-TEMP1*SN(I)*SN(K)*CS(J)                             
      DASDA(I,K+3)=-TEMP1*SN(I)*SN(J)*CS(K)                             
      DASDA(L,J+3)=SN(I)*ANG(K)/SNPROD+CS(I)*ANG(J)/SN(J)**2            
  120 DASDA(L,K+3)=SN(I)*ANG(J)/SNPROD+CS(I)*ANG(K)/SN(K)**2            
      RETURN                                                            
      END                                                               
C                                                                       
             SUBROUTINE CORR1(SIGZ,NN,II,JJ,CORR,KK,LL,MM)              
C                      WRITTEN BY L.W.FINGER                            
C     THIS ROUTINE EXTRACTS THE COVARIANCE ELEMENTS BETWEEN TWO ATOMS   
C     SIGZ = OUTPUT ARRAY  OF SIZE NN BY NN                             
C     II = ATOM  NUMBER OF FIRST ATOM                                   
C     JJ = ATOM NUMBER OF SECOND NUMBER                                 
C     CORR IS THE WHOLE POSITION VAR-COV MATRIX                         
C     KK AND LL DESCRIBE THE BLOCK OF SIGZ TO BE FILLED                 
C     MM = THE NUMBER OF ATOMS (NEEDED TO COMPUTE THE RANK OF CORR)     
                                                                        
      DIMENSION SIGZ(NN,NN),CORR(7260)                                  
                                                                        
      IK=MM                                                             
      IL=3*II-2                                                         
      IJ=3*JJ-2                                                         
      KI=KK-1                                                           
      LI=LL-1                                                           
      IF(II-JJ)1,2,3                                                    
    1 JI=((IJ-1)*IJ)/2+IL                                               
      IX=IJ                                                             
      GO TO 4                                                           
    3 JI=((IL-1)*IL)/2+IJ                                               
      IX=IL                                                             
    4 DO 5 I=1,3                                                        
      L=JI                                                              
      I1=I+LI                                                           
      DO 6 J=1,3                                                        
      J1=J+KI                                                           
      SIGZ(I1,J1)=CORR(L)                                               
    6 L=L+1                                                             
    5 JI=JI+IX+I-1                                                      
      RETURN                                                            
    2 SIGZ(LI+1,KI+1)=SIGZ(KI+1,KI+1)                                   
      DO 7 I=2,3                                                        
      I1=I+LI                                                           
      J1=I+KI                                                           
      SIGZ(I1,J1)=SIGZ(J1,J1)                                           
      IX=I-1                                                            
      DO 7 J=1,IX                                                       
      J2=J+KI                                                           
      J3=J+LI                                                           
      SIGZ(I1,J2)=SIGZ(J1,J2)                                           
    7 SIGZ(J3,J1)=SIGZ(I1,J2)                                           
      RETURN                                                            
      END                                                               
C                                                                       
                 SUBROUTINE DERDET(A,B,C)                               
                                                                        
C                    BY L.W.FINGER                                      
C     THIS ROUTINE COMPUTES THE DERIVATIVE OF A DETERMINANT             
                                                                        
      DIMENSION A(6),B(3,3),C(6,3,3)                                    
                                                                        
      DO 4 I=1,6                                                        
      A(I)=0.0                                                          
      DO 4 J=1,3                                                        
      M=MOD(J,3)+1                                                      
      N=MOD(J+1,3)+1                                                    
    4 A(I)=A(I)+C(I,J,J)*(B(N,N)*B(M,M)-B(N,M)*B(M,N))                  
     &  +C(I,J,N)*(B(M,J)*B(N,M)-B(M,M)*B(N,J))                         
     &  +C(I,J,M)*(B(N,J)*B(M,N)-B(N,N)*B(M,J))                         
      RETURN                                                            
      END                                                               
C                                                                       
                     SUBROUTINE ELVIB(IOUT)                             
                                                                        
C                        BY L.W. FINGER                                 
                                                                        
C     THIS  ROUTINE COMPUTES  ELLIPSOIDS OF  THERMAL  MOTION  AND  THE  
C     ERRORS  ASSOCIATED  WITH  THE  PARAMETERS  OF  THESE  ELLIPSOIDS  
                                                                        
      COMMON/ABLE/XYZ(3,60),BETA(6,60),A(3),ANG(3),AS(3),CS(3),ATOM(2,  
     & 60),SIGA(6,6),SIGB(6,6,60),RS(3,3,48),TS(3,48),TP(4,3),ITF(60),  
     & NA,NSS,VOL,DLIMIT(3),PI,G(3,3),NLP,IERR,TITLE(18),CORR(5000),
     & SIGX(128000)     
       
C            XYZ(3,60),BETA(6,60),A(3),ANG(3),AS(3),CS(3),ATOM(2,    
C     & 60), SIGA(6,6), SIGB(6,6,60),RS(3,3,48),TS(3,48),TP(4,3),ITF(60),    
C     & NA, NSS, VOL,  DLIMIT(3),  PI, G(3,3), NLP, IERR, TITLE(18)
C     &,SIGX(128000)
C                                                                 
      DIMENSION DDDC(6,3,3),DDDB(6,3,3),DTDA(6),DTDB(6),DSDA(6),DSDB(6),
     & DDETDA(6),DDETDB(6),DLDA(6),DLDB(6),DQDA(6,3),DQDB(6,3),         
     & DBQDA(6),B(3,3),D(3,3),ROOT(3),DBQDB(6),DASDA(6,6)               
     & ,THETA(3),STH(3),DTHDB(6,3),DTHDA(6,3),DXDA(6,3,3),DXDB(6,3,3),  
     & X(3,3),IAXIS(3),QI(3)                                            
                                                                        
      DATA IAXIS/1HA,1HB,1HC/                                           
                                                                        
 1013 FORMAT(1H1,18A4/11H ELLIPSOIDS)                                   
 1020 FORMAT(/10X,5HAXIS ,I1/18X,13HRMS AMPLITUDE,F18.6,16X,9HANGSTROMS)
 1025 FORMAT(14X,15HANGLE WITH THE ,A1,5H AXIS,F14.2,16X,7HDEGREES)     
 9001 FORMAT(1H1,18A4//22H ELLIPSOIDS AND ERRORS)                       
 9002 FORMAT(///10X2A3,25H IS NON POSITIVE DEFINITE)                    
 9003 FORMAT(///5X,27HELLIPSOID OF VIBRATION FOR ,2A3//20X,9HPARAMETER,
     &15X,5HVALUE,10X,5HERROR//13X,22HEQUIVALENT ISOTROPIC B,F14.3,
     &F15.3)  
 9004 FORMAT(/10X,5HAXIS ,I1/18X,13HRMS AMPLITUDE,F18.6,F15.6,10H ANGSTR
     &OMS )                                                             
 9005 FORMAT(14X,15HANGLE WITH THE ,A1,5H AXIS,F14.2,F15.2,8H DEGREES)  
                                                                        
      TWOPI=2.0*PI                                                      
      TFC=0.5/PI**2                                                     
      IF(IERR.NE.0)GO TO 1001                                           
      WRITE(IOUT,1013)TITLE                                             
      GO TO 1012                                                        
C---- GET DERIVATIVES OF RECIPROCAL CELL PARAMS WRT DIRECT PARAMS       
 1001 CALL CELLD(DASDA,A,ANG,AS,CS,VOL)                                 
      WRITE(IOUT,9001)TITLE                                             
C---- MATRIX D IS THE PRODUCT GB                                        
C     B IS THE MATRIX OF ANISOTROPIC COEFFICIENTS                       
C     G IS THE REAL SPACE METRIC                                        
C     DDDB(I,J,K) IS DERIVATIVE OF D(I,J) WRT BETA(K)                   
C     DDDC(I,J,K) IS DERIVATIVE OF D(I,J) WRT CELL(K)                   
C                                                                       
      DO 120 I=1,3                                                      
      DO 120 J=1,3                                                      
      DO 120 K=1,6                                                      
      DDDB(K,J,I)=0.                                                    
  120 DDDC(K,J,I)=0.0                                                   
      DO 121 I=1,3                                                      
      DO 122 J=1,3                                                      
  122 DDDB(I,J,I)=G(I,J)                                                
      DDDB(4,I,1)=G(I,2)                                                
      DDDB(6,I,3)=DDDB(4,I,1)                                           
      DDDB(5,I,1)=G(I,3)                                                
      DDDB(6,I,2)=DDDB(5,I,1)                                           
      DDDB(4,I,2)=G(I,1)                                                
  121 DDDB(5,I,3)=DDDB(4,I,2)                                           
C---- START LOOP THROUGH ATOM LIST                                      
 1012 DO 2 I=1,NA                                                       
      II=I                                                              
C---- SKIP OUT IF ATOM ISOTROPIC                                        
      IF(ITF(I))3,2,3                                                   
C---- EXPAND BETA LIST INTO SYMMETRIC MATRIX                            
C     B IS THE MATRIX OF THERMAL COEF.                                  
    3 DO 4 K=1,2                                                        
      B(K,K)=BETA(K,I)                                                  
      L=K+1                                                             
      DO 4 J=L,3                                                        
      JX2=K+J+1                                                         
      B(J,K)=BETA(JX2,I)                                                
    4 B(K,J)=B(J,K)                                                     
      B(3,3)=BETA(3,I)                                                  
C---- COMPUTE D MATRIX                                                  
C     D IS THE MATRIX PRODUCT GB                                         
      DO 5 K=1,3                                                         
      DO 5 L=1,3                                                         
      D(K,L)=0.                                                          
      DO 5 M=1,3                                                         
    5 D(K,L)=D(K,L)+G(K,M)*B(M,L)                                        
C---- SOLVE THIRD ORDER EQUATION TO GET RMS VALUES                       
C     T = TRACE OF MATRIX D                                             
C     S = SUM OF PRINCIPAL MINORS OF ORDER 2                            
C     DET = DETERMINANT OF D                                            
C                                                                       
      T=D(1,1)+D(2,2)+D(3,3)                                            
      S=D(1,1)*(D(2,2)+D(3,3))+D(2,2)*D(3,3)-D(1,2)*(D(2,1))-D(1,3)*D(3,
     & 1)-D(2,3)*D(3,2)                                                 
      DET=D(1,1)*(D(2,2)*D(3,3)-D(2,3)*D(3,2))-D(1,2)*(D(3,3)*D(2,1)-   
     & D(2,3)*D(3,1))+D(1,3)*(D(3,2)*D(2,1)-D(2,2)*D(3,1))              
C---- COMPUTE EQUIVALENT ISOTROPIC B                                    
      BISO=4.0/3.0*T                                                    
      DTDB(1)=S-T**2/3.0                                                
      DTDB(2)=(-2.0*T**3+9.0*T*S-27.0*DET)/27.0                         
      DTDB(3)=-DTDB(1)**3/27.0                                          
C---- CHECK IF ATOM POSITIVE-DEFINITE                                   
      IF(0.25*DTDB(2)**2-DTDB(3)-1.0E-6)1002,1002,7                     
 1002 DTDB(3)=-0.5*DTDB(2)/SQRT (DTDB(3))                               
      IF(ABS (DTDB(3))-1.0)1003,1003,6                                  
    6 DTDB(3)=SIGN (1.0,DTDB(3))                                        
 1003 PHI=ACOS (DTDB(3))                                                
      C=PHI/3.0                                                         
      DTDB(1)=2.0*SQRT (-DTDB(1)/3.0)                                   
      ROOT(1)=DTDB(1)*COS (C)+T/3.0                                     
      ROOT(2)=DTDB(1)*COS (C+TWOPI/3.0)+T/3.0                           
      ROOT(3)=DTDB(1)*COS (C-TWOPI/3.0)+T/3.0                           
C---- PLACE ROOTS IN INCREASING ORDER                                   
      DO 1004 J=1,2                                                     
      IJ=J+1                                                            
      DO 1004 K=IJ,3                                                    
      IF(ROOT(J).LE.ROOT(K))GO TO 1004                                  
 1005 HX=ROOT(J)                                                        
      ROOT(J)=ROOT(K)                                                   
      ROOT(K)=HX                                                        
 1004 CONTINUE                                                          
      IF(IERR)8,1006,8                                                  
C---- PRINT EQUIVALENT B FOR NO ERRORS                                  
 1006 WRITE(IOUT,9003)(ATOM(K,I),K=1,2),BISO                            
      GO TO 1007                                                        
C---- PRINT NON-POSITIVE DEFINITE MESSAGE                               
    7 WRITE(IOUT,9002)(ATOM(K,I),K=1,2)                                 
      GO TO 2                                                           
    8 DO 10 J=1,3                                                       
      DO 10 K=1,3                                                       
      M=MOD(K,3)+1                                                      
      N=MOD(K+1,3)+1                                                    
      DDDC(K,K,J)=D(K,J)/A(K)+A(K)*B(J,K)                               
      DDDC(K,M,J)=A(M)*B(K,J)*ANG(N)                                    
      DDDC(K,N,J)=A(N)*B(K,J)*ANG(M)                                    
      DDDC(K+3,M,J)=A(M)*A(N)*B(J,N)                                    
      DDDC(K+3,N,J)=A(N)*A(M)*B(J,M)                                    
   10 CONTINUE                                                          
C---- DTDA = DERIV OF TRACE WRT CELL PARAMS                             
C     DTDA = DERIV OF TRACE WRT TEMP FACTORS                            
C                                                                       
      DO 20 K=1,6                                                       
      DTDA(K)=0.                                                        
      DTDB(K)=0.0                                                       
      DO 18 J=1,3                                                       
      DTDB(K)=DTDB(K)+DDDB(K,J,J)                                       
   18 DTDA(K)=DTDA(K)+DDDC(K,J,J)                                       
C---- DSDA AND DSDB ARE DERIVS OF S WRT CELL PARAMS AND TEMP FACTORS    
      DSDA(K)=0.                                                        
      DSDB(K)=0.0                                                       
      DO 20 J=1,3                                                       
      DO 20 L=1,3                                                       
      IF(L.EQ.J)GO TO 20                                                
   22 DSDA(K)=DSDA(K)+D(J,J)*DDDC(K,L,L)-D(J,L)*DDDC(K,L,J)             
   21 DSDB(K)=DSDB(K)+D(J,J)*DDDB(K,L,L)-D(J,L)*DDDB(K,L,J)             
   20 CONTINUE                                                          
C---- COMPUTE DERIVATIVES OF DETERMINANT WRT CELL AND THERMAL PARAMS    
C     DERDET EVALUATES THE DERIVATIVE OF THE DETERMINANT                
C                                                                       
      CALL DERDET(DDETDA,D,DDDC)                                        
      CALL DERDET(DDETDB,D,DDDB)                                        
C---- COMPUTE SIGMA FOR EQ ISO B                                        
      CALL VARIAN(DELTA,DTDA,DTDB,SIGA,SIGB,II)                         
      DELTA=4.0/3.0*DELTA                                               
      WRITE(IOUT,9003)(ATOM(K,I),K=1,2),BISO,DELTA                      
C---- START FINDING EIGENVECTORS AND ERRORS                             
 1007 DO 98 K=1,3                                                       
      DO 98 J=1,3                                                       
      X(K,J)=D(K,J)                                                     
      DO 98 L=1,6                                                       
      DXDA(L,J,K)=DDDC(L,J,K)                                           
   98 DXDB(L,J,K)=DDDB(L,J,K)                                           
      DO 19 N=1,3                                                       
      DEN=    (2.0*T*ROOT(N)-S-3.0*ROOT(N)**2)                          
      IF(ABS(DEN).GT.0.0)DEN=1./DEN                                     
      IF(IERR.EQ.0)GO TO 1018                                           
      DO 23 K=1,6                                                       
      DLDA(K)=DEN*(-DDETDA(K)+ROOT(N)*DSDA(K)-ROOT(N)**2*DTDA(K))       
   23 DLDB(K)=DEN*(-DDETDB(K)+ROOT(N)*DSDB(K)-ROOT(N)**2*DTDB(K))       
      CALL VARIAN(DELTA,DLDA,DLDB,SIGA,SIGB,II)                         
C---- DELTA IS THE SIGMA OF THE ROOT                                    
 1018 IF(ROOT(N))7,25,25                                                
   25 ROOT1=SQRT (ROOT(N)*0.5)/PI                                       
      IF(IERR.EQ.0)GO TO 1016                                           
      DROOT=DELTA/(4.0*PI**2*ROOT1)                                     
      WRITE(IOUT,9004)N,ROOT1,DROOT                                     
      GO TO 1030                                                        
 1016 WRITE(IOUT,1020)N,ROOT1                                           
 1030 DO 100 K=1,3                                                      
      X(K,K)=D(K,K)-ROOT(N)                                             
      DO 100 L=1,6                                                      
      DXDA(L,K,K)=DDDC(L,K,K)-DLDA(L)                                   
  100 DXDB(L,K,K)=DDDB(L,K,K)-DLDB(L)                                   
      QMAG=0.                                                           
      DO 104 L=1,3                                                      
      M=MOD(L,3)+1                                                      
      STH(1)=(X(L,2)*X(M,3)-X(L,3)*X(M,2))*AS(1)                        
      STH(2)=(X(L,3)*X(M,1)-X(L,1)*X(M,3))*AS(2)                        
      STH(3)=(X(L,1)*X(M,2)-X(L,2)*X(M,1))*AS(3)                        
      TEMP=STH(1)**2+STH(2)**2+STH(3)**2                                
      IF(TEMP-QMAG)104,102,102                                          
  102 QMAG=TEMP                                                         
      K=L                                                               
      DO 103 M=1,3                                                      
  103 QI(M)=STH(M)                                                      
  104 CONTINUE                                                          
      QMAG=(QMAG+2.0*(QI(1)*QI(2)*CS(3)+QI(1)*QI(3)*CS(2)+QI(2)*QI(3)*  
     & CS(1)))                                                          
      DO 7531 M=1,3                                                     
 7531 QI(M)=QI(M)/AS(M)                                                 
      IF(QMAG.GT.0.0)GO TO 2001                                         
      QMAG=1.0E-3                                                       
 2001 QMAG=SQRT(QMAG)                                                   
      IF(IERR.EQ.0) GO TO 1021                                          
      M=MOD(K,3)+1                                                      
      DO 107 L=1,6                                                      
      DO 107 NN=1,3                                                     
      NM=MOD(NN,3)+1                                                    
      MN=MOD(NM,3)+1                                                    
  110 DQDA(L,NN)=X(K,NM)*DXDA(L,M,MN)+X(M,MN)*DXDA(L,K,NM)-X(K,MN)*DXDA(
     & L,M,NM)-X(M,NM)*DXDA(L,K,MN)                                     
  107 DQDB(L,NN)=X(K,NM)*DXDB(L,M,MN)+X(M,MN)*DXDB(L,K,NM)-X(K,MN)*DXDB(
     & L,M,NM)-X(M,NM)*DXDB(L,K,MN)                                     
      K=1                                                               
   27 DBQDB(K)=0.                                                       
      DBQDA(K)=0.                                                       
      DO 26 J=1,3                                                       
      DBQDB(K)=DBQDB(K)+2.0*QI(J)*AS(J)**2*DQDB(K,J)                    
      DBQDA(K)=DBQDA(K)+2.0*QI(J)*AS(J)**2*DQDA(K,J)+2.0*AS(J)*QI(J)**2 
     & *DASDA(J,K)                                                      
      DO 26 L=1,3                                                       
      IF(L.EQ.J) GO TO 26                                               
      JX2=6-L-J                                                         
      DBQDB(K)=DBQDB(K)+AS(J)*AS(L)*(QI(L)*DQDB(K,J)+QI(J)*DQDB(K,L))   
     & *CS(JX2)                                                         
      DBQDA(K)=DBQDA(K)+AS(J)*AS(L)*(QI(J)*QI(L)*DASDA(JX2+3,K)+CS(JX2)*
     & (QI(J)*DQDA(K,L)+QI(L)*DQDA(K,J)))+QI(J)*QI(L)*CS(JX2)*(AS(J)*   
     & DASDA(L,K)+AS(L)*DASDA(J,K))                                     
   26 CONTINUE                                                          
      DBQDA(K)=0.5*DBQDA(K)/QMAG                                        
      DBQDB(K)=0.5*DBQDB(K)/QMAG                                        
      K=K+1                                                             
      IF(K.LE.6) GO TO 27                                               
 1021 DO 30 K=1,3                                                       
      THETA(K)=QI(K)/(A(K)*QMAG)                                        
      IF(ABS (THETA(K))-1.0)32,32,31                                    
   31 THETA(K)=SIGN (1.0,THETA(K))                                      
   32 STH(K)=SQRT (1.0-THETA(K)**2)                                     
   30 THETA(K)=180.*ACOS (THETA(K))/PI                                  
      IF(IERR.EQ.0)GO TO 1023                                           
 1024 DO 33 K=1,3                                                       
      DO 34 J=1,6                                                       
C---- THETA IS THE ANGLE WITH THE I TH REAL AXIS                        
      DTHDB(J,K)=DQDB(J,K)/(A(K)*QMAG)-QI(K)*DBQDB(J)/(A(K)*QMAG**2)    
   34 DTHDA(J,K)=DQDA(J,K)/(A(K)*QMAG)-QI(K)*DBQDA(J)/(A(K)*QMAG**2)    
      DTHDA(K,K)=DTHDA(K,K)-QI(K)/(A(K)**2*QMAG)                        
      DO 33 J=1,6                                                       
      IF(ABS (STH(K))-1.0E-3)150,150,151                                
 151  DTHDA(J,K)=-DTHDA(J,K)/STH(K)                                     
      DTHDB(J,K)=-DTHDB(J,K)/STH(K)                                     
      GO TO 33                                                          
  150 DTHDA(J,K)=0.                                                     
      DTHDB(J,K)=0.                                                     
   33 CONTINUE                                                          
      DO 35 K=1,3                                                       
      CALL VARIAN(STH(K),DTHDA(1,K),DTHDB(1,K),SIGA,SIGB,II)            
   35 STH(K)=STH(K)*180.0/PI                                            
      WRITE(IOUT,9005)(IAXIS(K),THETA(K),STH(K),K=1,3)                  
      GO TO 19                                                          
 1023 WRITE(IOUT,1025)(IAXIS(K),THETA(K),K=1,3)                         
   19 CONTINUE                                                          
    2 CONTINUE                                                          
      RETURN                                                            
      END                                                               
C                                                                       
            SUBROUTINE LIBCOR(XYZ,PAR,KPAR,NPAR,G,NA,MDEGN,KNSTR)       
C                                                                       
C                           WRITTEN BY L.W.FINGER                       
C                                                                       
C     SUBROUTINE  TO APPLY  THE  LIBRATION  CORRECTION TO THE POSITION  
C     PARAMETERS, IN PREPARATION FOR CALCULATION OF BOND DISTANCES AND  
C     ANGLES.    IT IS ASSUMED THAT THE CONTRIBUTIONS TO THE ERRORS IN  
C     DISTANCES  AND ANGLES  DUE TO THE ERRORS IN DETERMINATION OF THE  
C     LIBRATION PARAMETERS IS  NEGLEGIBLE  COMPARED WITH THE ERRORS IN  
C     THE DETERMINATION OF THE POSITION PARAMETERS.                     
C                                                                       
C                                                                       
      DIMENSION XYZ(3,60),PAR(60),KPAR(60),G(3,3),KNSTR(MDEGN),ATOX(6)  
     & ,XTOA(6),RC(3),R(3),EL(3,3),OXYZ(3),SXYZ(3),L(6)                 
C
      I1(K)=MOD(K/256**3,256)                                           
      I2(K)=MOD(K/256**2,256)                                           
      I3(K)=MOD(K/256,256)                                              
      I4(K)=MOD(K,256)                                                  
C                                                                       
C---- COMPUTE TRANS FORMATION FROM CRYSTAL TO ORTHONORMAL AND BACK.     
      ATOX(1)=SQRT(G(1,1))                                              
      ATOX(2)=G(1,2)/ATOX(1)                                            
      ATOX(3)=G(1,3)/ATOX(1)                                            
      ATOX(4)=SQRT(G(2,2)-ATOX(2)**2)                                   
      ATOX(5)=(G(2,3)-ATOX(2)*ATOX(3))/ATOX(4)                          
      ATOX(6)=SQRT(G(3,3)-ATOX(3)**2-ATOX(5)**2)                        
      XTOA(1)=1./ATOX(1)                                                
      XTOA(4)=1./ATOX(4)                                                
      XTOA(6)=1./ATOX(6)                                                
      XTOA(2)=-ATOX(2)/(ATOX(1)*ATOX(4))                                
      XTOA(5)=-ATOX(5)/(ATOX(4)*ATOX(6))                                
      XTOA(3)=(ATOX(2)*ATOX(5)/ATOX(4)-ATOX(3))/(ATOX(1)*ATOX(6))       
C---- DETERMINE WHETHER ANY SPECIAL PARAMETERS ARE L TENSOR ELEMENTS.   
      ISET=-99                                                          
      IST=1                                                             
   50 DO 100 I=IST,NPAR                                                 
      IF(I1(KPAR(I)).EQ.1.AND.(I4(KPAR(I)).GE.16.AND.I4(KPAR(I)).LE.21) 
     & .AND.I2(KPAR(I)).NE.ISET) GO TO 110                              
  100 CONTINUE                                                          
      RETURN                                                            
C---- THERE IS AT LEAST 1 SET OF L TENSOR ELEMENTS. FIND THE MEMBERS OF 
C     THE SET AND FILL THEM IN.                                         
  110 ISET=I2(KPAR(I))                                                  
      DO 120 J=1,3                                                      
      OXYZ(J)=0.                                                        
      DO 120 K=1,3                                                      
  120 EL(J,K)=0.                                                        
      IL=0                                                              
      DO 130 J=1,NPAR                                                   
      IF(I2(KPAR(J)).NE.ISET.OR.I1(KPAR(J)).NE.1) GO TO 130             
      IDEN=I4(KPAR(J))                                                  
      IF(IDEN.GE.1.AND.IDEN.LE.3)GO TO 124                              
      IF(IDEN.GE.16.AND.IDEN.LE.21)GO TO 126                            
      GO TO 130                                                         
  124 OXYZ(IDEN)=PAR(J)                                                 
      GO TO 130                                                         
  126 IDEN=IDEN-15                                                      
      IF(IDEN-3)127,127,128                                             
  127 EL(IDEN,IDEN)=PAR(J)                                              
      GO TO 129                                                         
  128 JJ=1                                                              
      KK=3                                                              
      IF(IDEN.EQ.6)JJ=2                                                 
      IF(IDEN.EQ.4)KK=2                                                 
      EL(JJ,KK)=PAR(J)                                                  
      EL(KK,JJ)=PAR(J)                                                  
  129 IL=IL+1                                                           
      L(IL)=J                                                           
  130 CONTINUE                                                          
C---- WE MUST NOW FIND WHICH ATOMS DEPEND ON L IN THIS SET. ARRAY L     
C     CONTAINS A LIST OF THE PARAMETERS WHICH ARE MEMBERS OF THE L      
C     TENSOR IN THIS SET.                                               
      DO 1000 I=1,NA                                                    
      DO 140 J=1,IL                                                     
      JJ=L(J)                                                           
      DO 140 K=1,MDEGN,2                                                
      IF(I1(KNSTR(K)).NE.0.OR.I2(KNSTR(K)).GT.L(IL))GO TO 1000          
      IF(I2(KNSTR(K)).EQ.JJ.AND.I3(KNSTR(K)).EQ.I)GO TO 150             
  140 CONTINUE                                                          
C---- NO MATCH FOUND FOR THIS ATOM                                      
      GO TO 1000                                                        
C---- ATOM I DEPENDS ON THE CURRENT SET. APPLY LIBRATION CORRECTION.    
  150 DO 160 J=1,3                                                      
  160 SXYZ(J)=XYZ(J,I)-OXYZ(J)                                          
      CALL UTMUL(ATOX,SXYZ,3,RC)                                        
      CALL CORLIB(R,EL,RC)                                              
      CALL UTMUL(XTOA,R,3,SXYZ)                                         
      DO 170 J=1,3                                                      
  170 XYZ(J,I)=SXYZ(J)+OXYZ(J)                                          
 1000 CONTINUE                                                          
      IST=L(IL)+1                                                       
      IF(IST.GT.NPAR)RETURN                                             
      GO TO 50                                                          
      END                                                               
C                                                                       
                       SUBROUTINE  UTMUL(A,B,N,C)                       
C                                                                       
C                        WRITTEN BY L.W.FINGER                          
C                                                                       
C     MULTIPLIES VECTOR B BY UPPER TRIANGULAR MATRIX A, STORING THE     
C                 RESULT IN C. A IS STORED BY ROWS.                     
C                                                                       
      DIMENSION A(1),B(1),C(1)                                          
      K=0                                                               
      DO 100 I=1,N                                                      
      C(I)=0.                                                           
      DO 100 J=I,N                                                      
      K=K+1                                                             
  100 C(I)=C(I)+A(K)*B(J)                                               
      RETURN                                                            
      END                                                               
C                                                                       
                SUBROUTINE VARIAN(A,B,C,D,E,II)                         
C                    WRITTEN BY L.W.FINGER                              
                                                                        
      DIMENSION B(6),C(6),D(6,6)   ,E(6,6,40)                           
                                                                        
      N=II                                                              
      F=B(1)**2*D(1,1)+C(1)**2*E(1,1,N)                                 
      DO 1 I=2,6                                                        
      F=F+B(I)**2*D(I,I)+C(I)**2*E(I,I,N)                               
      IJ=I-1                                                            
      DO 1 J=1,IJ                                                       
    1 F=F+2.*B(I)*B(J)*D(I,J)+2.0*C(I)*C(J)*E(I,J,N)                    
      IF(F)2,3,3                                                        
    2 F=0.0                                                             
    3 A=SQRT (F)                                                        
      RETURN                                                            
      END                                                               
C                                                                       
        SUBROUTINE CORECT(R3,R1,R2,ELT,THETA,RC,BIGR,BIGRP,RP2)         
C                                                                       
C                      WRITTEN BY L.W.FINGER                            
C                                                                       
      DIMENSION R3(3),R1(3),R2(3),RC(3),BIGR(3),BIGRP(3),ELT(3,3),      
     & RP(3),RPP(3),ACON(3,3),BCON(3)                                   
C                                                                       
      RP2=0.                                                            
      DO 10 I=1,3                                                       
      RP(I)=R2(I)-R1(I)                                                 
   10 RP2=RP2+RP(I)**2                                                  
      ISIG=1                                                            
      DO 20 I33K=1,3                                                    
      I=I33K                                                            
      BCON(I)=RC(I)                                                     
      DO 20 J33K=1,3                                                    
      J=J33K                                                            
      ACON(I,J)=FIJ(THETA,ELT,RP,RP2,I,J)                               
   20 BCON(I)=BCON(I)+ACON(I,J)*R2(J)                                   
      GO TO 40                                                          
C                                                                       
C                                                                       
      E N T R Y  CORLIB(R3,ELT,RC)                                      
C                                                                       
C                                                                       
      DO 30 I=1,3                                                       
      BCON(I)=RC(I)                                                     
      DO 30 J=1,3                                                       
   30 ACON(I,J)=0.                                                      
      ISIG=0                                                            
   40 DO 50 I33K=1,3                                                    
      I=I33K                                                            
      DO 50 J33K=1,3                                                    
      J=J33K                                                            
   50 ACON(I,J)=ACON(I,J)+EIJ(ELT,I,J)                                  
      CALL SOLVE(ACON,BCON,R3)                                          
      IF(ISIG)59,59,60                                                  
   59 RETURN                                                            
   60 DO 70 I=1,3                                                       
   70 RPP(I)=R3(I)-R2(I)                                                
      CALL VXV(RP,RPP,BIGR)                                             
      CALL VXV(RP,BIGR,BIGRP)                                           
      RETURN                                                            
      END                                                               
C                                                                       
                      SUBROUTINE  SOLVE(A,B,C)                          
C                                                                       
C                      WRITTEN BY L.W.FINGER                            
C                                                                       
      DIMENSION A(3,3),B(3),C(3)                                        
C                                                                       
      DELTA=DETERM(A(1,1),A(1,2),A(1,3))                                
      C(1)=DETERM(B(1),A(1,2),A(1,3))/DELTA                             
      C(2)=DETERM(A(1,1),B(1),A(1,3))/DELTA                             
      C(3)=DETERM(A(1,1),A(1,2),B(1))/DELTA                             
      RETURN                                                            
      END                                                               
C                                                                       
                  FUNCTION FIJ(THETA,EL,RP,RP2,I,J)                     
C                                                                       
C                      WRITTEN BY L.W.FINGER                            
C                                                                       
      DIMENSION EL(3,3),RP(3),CON(3)                                    
C                                                                       
      IF(RP2.LT.0.01.OR.THETA.LT.0.000001)GO TO 100                     
      K=MOD(I,3)+1                                                      
      L=MOD(I+1,3)+1                                                    
      KK=MOD(J,3)+1                                                     
      LL=MOD(J+1,3)+1                                                   
      M2=MOD(J-I+2,3)+1                                                 
      M3=MOD(J-I+3,3)+1                                                 
      M4=MOD(J-I+4,3)+1                                                 
      CON1=1.-THETA/4.-.5*(EL(K,K)+EL(L,L))                             
      CON2=.5*EL(I,K)                                                   
      CON3=.5*EL(I,L)                                                   
      CON(M3)=-(RP(KK)**2+RP(LL)**2)                                    
      CON(M4)=RP(J)*RP(KK)                                              
      CON(M2)=RP(J)*RP(LL)                                              
      FIJ=(CON1*CON(1)+CON2*CON(2)+CON3*CON(3))*THETA/(2.*RP2)          
      RETURN                                                            
  100 FIJ=0.                                                            
      RETURN                                                            
      END                                                               
C                                                                       
                       SUBROUTINE VXV(V1,V2,V3)                         
C                                                                       
C                       WRITTEN BY L.W.FINGER                           
C                                                                       
      DIMENSION V1(3),V2(3),V3(3)                                       
C                                                                       
      V3(1)=V1(2)*V2(3)-V1(3)*V2(2)                                     
      V3(2)=V1(3)*V2(1)-V1(1)*V2(3)                                     
      V3(3)=V1(1)*V2(2)-V1(2)*V2(1)                                     
      RETURN                                                            
      END                                                               
C                                                                       
                        FUNCTION EIJ(EL,I,J)                            
C                                                                       
C                      WRITTEN BY L.W. FINGER                           
C                                                                       
      DIMENSION EL(3,3)                                                 
C                                                                       
      IF(I.EQ.J)GO TO 20                                                
      EIJ=.5*EL(I,J)-(1./24.)*(EL4(EL,I,1,1,J)+EL4(EL,I,2,2,J)+EL4(EL,  
     & I,3,3,J))                                                        
      RETURN                                                            
   20 K=MOD(I,3)+1                                                      
      L=MOD(I+1,3)+1                                                    
      EIJ=1.-.5*(EL(K,K)+EL(L,L))+(1./24.)*(EL4(EL,I,I,K,K)+EL4(EL,I,I, 
     & L,L)+3.*(EL(K,K)**2+EL(L,L)**2)+2.*EL4(EL,K,K,L,L))              
      RETURN                                                            
      END                                                               
C                                                                       
                    FUNCTION  DETERM(C1,C2,C3)                          
C                                                                       
C                     WRITTEN BY L.W.FINGER                             
C                                                                       
      DIMENSION C1(3),C2(3),C3(3)                                       
C                                                                       
      DETERM=C1(1)*C2(2)*C3(3)+C1(2)*C2(3)*C3(1)+C1(3)*C2(1)*C3(2)      
     & -C1(3)*C2(2)*C3(1)-C1(2)*C2(1)*C3(3)-C1(1)*C2(3)*C3(2)           
      RETURN                                                            
      END                                                               
C                                                                       
                     FUNCTION EL4(EL,I,J,K,L)                           
C                                                                       
C                     WRITTEN BY L.W.FINGER                             
C                                                                       
      DIMENSION EL(3,3)                                                 
      EL4=EL(I,J)*EL(K,L)+EL(I,K)*EL(J,L)+EL(I,L)*EL(J,K)               
      RETURN                                                            
      END                                                               

      SUBROUTINE TENSAL(IOUT,NATOM)
C---- ANALYSE HIGHER ORDER THERMAL MOTION TENSORS
C
      COMMON /ABLE/ XYZ(3,60),BETA(6,60),A(3),ANG(3),AS(3),CS(3),
     &ATOM(2,60), SIGA(6,6),SIGB(6,6,60),RS(3,3,48),TS(3,48),
     &TP(4,3),ITF(60),NA,NSS,VOL,DLIMIT(3),PI,G(3,3),NLP,IERR,
     &TITLE(18),SIGX(128000),CORR(5000)
      COMMON /TENSOR/ CR(10,60),DR(15,60),ER(21,6),FR(28,6),
     &CRF(3,3,3),DRF(3,3,3,3),ERF(3,3,3,3,3)
      COMMON /METRIC/ AG,BG,CG,WINKEL(3),ARE(3),WINRE(3)
      PI = ACOS(-1.)
C
C---- PRINT CELL INFORMATION
C
      WRITE(IOUT,10) AG,BG,CG,WINKEL,ARE,WINRE,G
   10 FORMAT(' DIRECT CELL ',6F10.6,/, ' RECIP. CELL ',6F10.6,
     &//,' METRIC TENSOR ',3(F10.6,1X),/,2(15X,3(F10.6,1X),/)//)
C
C---- LOOP OVER ATOMS
C
      DO 1000 IATOM = 1,NATOM   
      WRITE(IOUT,20) (ATOM(J,IATOM),J=1,2)
   20 FORMAT(1H0,60(1H*),' ATOM ',A4,A2)
C
C---- TRANSFORM TENSOR COMPONENTS INTO DISPLACEMENTS IN AE
C
      CALL UTRA(IOUT)
C
C---- EXPAND LIST INTO FULL TENSOR
C
      CALL TENFUL(IATOM,ICHK)
      IF(ICHK.EQ.0) GOTO 1000
C
C---- CALCULATE SKEW DIVERGENCE FROM RANK 3
C
      CALL SKEWDI(IOUT,IATOM)
C
C---- CALCULATE TRACE OF RANK 4
C
      CALL TRACE4(IOUT,IATOM)
C
C---- TRANSFORM TO CARTESIAN SYSTEM (AND CALCULATE INVARIANTS)
C
      CALL TRACAR(IOUT,IATOM)
 1000 CONTINUE
      RETURN 
      END
      SUBROUTINE UTRA(IOUT)
      RETURN
      END
      SUBROUTINE TENFUL(I,ICHK)
      COMMON /TENSOR/ CR(10,60),DR(15,60),ER(21,6),FR(28,6),
     &CRF(3,3,3),DRF(3,3,3,3),ERF(3,3,3,3,3)
C
C---- EXPAND LIST INTO FULL TENSOR FOR RANK 3
C
      ICHK = 1
      CHKSUM = 0.0
      DO 300 J=1,10
      CHKSUM = CHKSUM + ABS(CR(J,I))
  300 CR(J,I) = CR(J,I)*0.001
      DO 3000 K=1,3
 3000 CRF(K,K,K) = CR(K,I)
      DO 3100 K=4,10
      K1 = 1
      IF(K.GE.8) K1 = 2
      IF(K.LE.9) K2 = K - 6
      IF(K.LE.6) K2 = K - 5
      IF(K.EQ.7) K2 = 3
      IF(K.LE.5) K2 = K - 3
      IF(K.LE.10) K3 = 3
      IF(K.LE.5) K3 = 2
      IF(K.EQ.10) K1 = 1
      IF(K.EQ.10) K2 = 2
      CRF(K1,K2,K3) = CR(K,I)
      CRF(K1,K3,K2) = CR(K,I)
      CRF(K3,K2,K1) = CR(K,I)
      CRF(K3,K1,K2) = CR(K,I)
      CRF(K2,K1,K3) = CR(K,I)
      CRF(K2,K3,K1) = CR(K,I)
 3100 CONTINUE
C
C---- EXPAND LIST INTO FULL TENSOR FOR RANK 4
C
      DO 400 J = 1,15
      CHKSUM = CHKSUM + ABS(DR(J,I))
  400 DR(J,I) = DR(J,I)*0.0001
      DO 4000 K = 1,3
 4000 DRF(K,K,K,K) = DR(K,I)
      DO 4100 K = 4,15
      K1 = 1
      IF(K.GT.12) K1 = 2
      K2 = 1
      IF(K.GT.8) K2 = 2
      IF(K.GT.11) K2 = 3
      IF(K.EQ.13.OR.K.EQ.14) K2 = 2
      K3 = 1
      IF(K.GT.5) K3 = 2
      IF(K.GT.13) K3 = 3
      IF(K.EQ.8.OR.K.EQ.11.OR.K.EQ.12) K3 = 3
      K4 = 2
      IF(K.GT.4) K4 = 3
      IF(K.EQ.6.OR.K.EQ.9) K4 = 2
      DO 4150 JJ = 1,2
      IF(JJ.EQ.1) GOTO 4140
      KS = K1
      K1 = K3
      K3 = KS
      KS = K2
      K2 = K4
      K4 = KS
 4140 CONTINUE
      DRF(K1,K2,K3,K4) = DR(K,I)
      DRF(K1,K2,K4,K3) = DR(K,I)
      DRF(K1,K3,K2,K4) = DR(K,I)
      DRF(K1,K3,K4,K2) = DR(K,I)
      DRF(K1,K4,K2,K3) = DR(K,I)
      DRF(K1,K4,K3,K2) = DR(K,I)
      DRF(K2,K1,K3,K4) = DR(K,I)
      DRF(K2,K1,K4,K3) = DR(K,I)
      DRF(K2,K4,K1,K3) = DR(K,I)
      DRF(K2,K4,K3,K1) = DR(K,I)
      DRF(K2,K3,K4,K1) = DR(K,I)
      DRF(K2,K3,K1,K4) = DR(K,I)
 4150 CONTINUE
 4100 CONTINUE
      IF(CHKSUM.LT.1.E-8) ICHK = 0
      RETURN
      END
      SUBROUTINE SKEWDI(IOUT,IATOM)
C
C---- CALCULATE VECTOR INVARIANT G(REAL)JK * C IJK
C
      DOUBLE PRECISION UI(3)
      COMMON /ABLE/ XYZ(3,60),BETA(6,60),A(3),ANG(3),AS(3),CS(3),
     &ATOM(2,60), SIGA(6,6),SIGB(6,6,60),RS(3,3,48),TS(3,48),
     &TP(4,3),ITF(60),NA,NSS,VOL,DLIMIT(3),PI,G(3,3),NLP,IERR,
     &TITLE(18),SIGX(128000),CORR(5000)
      COMMON /TENSOR/ CR(10,60),DR(15,60),ER(21,6),FR(28,6),
     &CRF(3,3,3),DRF(3,3,3,3),ERF(3,3,3,3,3)
      COMMON /METRIC/ AG,BG,CG,WINKEL(3),ARE(3),WINRE(3)
      DIMENSION VI(3)
      DO 10 I=1,3
   10 VI(I) = 0.0
      DO 1000 I = 1,3
      DO 1000 J = 1,3
      DO 1000 K = 1,3
      VI(I) = VI(I) + G(J,K) * CRF(I,J,K)
 1000 CONTINUE
      WRITE (IOUT,990) (ATOM(J,IATOM),J=1,2),VI
  990 FORMAT (' VI FOR ATOM ',A4,A2,3F9.6)
C
C---- NOW TRANSFORM INTO AE**3 UNITS AND MULTIPLY WITH 3/5
C
      DO 2000 I = 1,3
      UI(I) = DBLE(0.6 * VI(I)/ARE(I)) 
 2000 CONTINUE
      WRITE(IOUT,991) (ATOM(J,IATOM),J=1,2),UI
  991 FORMAT (' VECTOR OF SKEW DIV. (IN AE**3) ',A4,A2,3F9.6)
C
C---- GET EQUIVALENT IN AE
C
      DO 3000 I = 1,3
      UI(I) = DABS(UI(I)**2)**(DBLE(1./6.))
 3000 CONTINUE
      WRITE(IOUT,992) (ATOM(J,IATOM),J=1,2),UI
  992 FORMAT (' MAGNITUDE OF SKEW DIV. (IN AE) ',A4,A2,3F9.6)
      RETURN
      END
      SUBROUTINE TRACE4(IOUT,IATOM)
C
C---- CALCULATE TRACE OF 4TH RANK TENSOR     
C
      DOUBLE PRECISION TD0
      COMMON /ABLE/ XYZ(3,60),BETA(6,60),A(3),ANG(3),AS(3),CS(3),
     &ATOM(2,60), SIGA(6,6),SIGB(6,6,60),RS(3,3,48),TS(3,48),
     &TP(4,3),ITF(60),NA,NSS,VOL,DLIMIT(3),PI,G(3,3),NLP,IERR,
     &TITLE(18),SIGX(128000),CORR(5000)
      COMMON /TENSOR/ CR(10,60),DR(15,60),ER(21,6),FR(28,6),
     &CRF(3,3,3),DRF(3,3,3,3),ERF(3,3,3,3,3)
      COMMON /METRIC/ AG,BG,CG,WINKEL(3),ARE(3),WINRE(3)
      TD0 = 0.0
      DO 4000 I=1,3
      DO 4000 J=1,3
      DO 4000 K=1,3
      DO 4000 L=1,3
      TD0 = TD0 + G(I,J)*G(K,L)*DRF(I,J,K,L)
 4000 CONTINUE
      WRITE (IOUT,900) (ATOM(J,IATOM),J=1,2),TD0
  900 FORMAT(' TRACE OF 4TH RANK TENSOR FOR ',A4,A2,F9.6)
C
C---- CALCULATE EQUIVALENT ISOTROPIC RMSD FROM TD0
C
      TD0 = (DABS(TD0)/(3.*5.))**0.25
      WRITE (IOUT,910) (ATOM(J,IATOM),J=1,2),TD0
  910 FORMAT(' EQUIVALENT RMSD FROM TRACE FOR ',A4,A2,F9.6,' AE')
      RETURN
      END
      SUBROUTINE TRACAR(IOUT,IATOM)
C
C---- TRANSFORM TO CARTESIAN SYSTEM (S. WILLIS & PRYOR P.100)
C
      DOUBLE PRECISION U4E11R,U4E22R,U4E33R,U4E11,U4E22,U4E33,
     &U4E12,U4E13,U4E23
      DIMENSION U3(3,3,3),U4(3,3,3,3),F(3,3),ANGL(3),ANGRE(3)
      DIMENSION BRF(3,3),U2(3,3),U2E(3,3),D(3),V(3,3),U4E(3,3,3,3)
      COMMON /ABLE/ XYZ(3,60),BETA(6,60),A(3),ANG(3),AS(3),CS(3),
     &ATOM(2,60), SIGA(6,6),SIGB(6,6,60),RS(3,3,48),TS(3,48),
     &TP(4,3),ITF(60),NA,NSS,VOL,DLIMIT(3),PI,G(3,3),NLP,IERR,
     &TITLE(18),SIGX(128000),CORR(5000)
      COMMON /TENSOR/ CR(10,60),DR(15,60),ER(21,6),FR(28,6),
     &CRF(3,3,3),DRF(3,3,3,3),ERF(3,3,3,3,3)
      COMMON /METRIC/ AG,BG,CG,WINKEL(3),ARE(3),WINRE(3)
C
C---- FILL IN TENSOR BETAIJ
C
      DO 4 I=1,3
    4 BRF(I,I) = BETA(I,IATOM)/(2.*PI*PI)
      BRF(1,2) = BETA(4,IATOM)/(2.*PI*PI)
      BRF(2,1) = BRF(1,2)
      BRF(1,3) = BETA(5,IATOM)/(2.*PI*PI)
      BRF(3,1) = BRF(1,3)
      BRF(2,3) = BETA(6,IATOM)/(2.*PI*PI)
      BRF(3,2) = BRF(2,3)
C
C---- FIRST FILL IN TRANSFORMATION MATRIX
C
      DO 5 I=1,3
      ANGL(I) = WINKEL(I)*PI/180.
    5 ANGRE(I) = ACOS(WINRE(I))
      DO 10 I=1,3
      DO 10 J=1,3
   10 F(I,J) = 0.0
       F(1,1) = 1./ARE(1)
       F(2,1) = -1./(ARE(1)*TAN(ANGRE(3)))
       F(3,1) = AG * COS(ANGL(2))
       F(2,2) = 1./(ARE(2)*SIN(ANGRE(3)))
       F(3,2) = BG * COS(ANGL(1))
       F(3,3) = CG
C
C---- INITIALIZE
C
      DO 50 IT=1,3
      DO 50 JT=1,3
      U2(IT,JT) = 0.0
      DO 50 KT=1,3
      U3(IT,JT,LT) = 0.0
      DO 50 LT=1,3
      U4(IT,JT,KT,LT) = 0.0
   50 CONTINUE
C
C---- TRANSFORM TENSOR COMPONENTS
C
      DO 90 IT=1,3
      DO 90 JT=1,3
      DO 90 I=1,3
      DO 90 J=1,3
      U2(IT,JT) = U2(IT,JT) + F(IT,I)*F(JT,J) * BRF(I,J)
   90 CONTINUE
      DO 100 IT=1,3
      DO 100 JT=1,3
      DO 100 KT=1,3
      DO 100 I=1,3
      DO 100 J=1,3
      DO 100 K=1,3
      U3(IT,JT,KT) = U3(IT,JT,KT) + F(IT,I)*F(JT,J)*F(KT,K) * CRF(I,J,K)
  100 CONTINUE
      DO 110 IT=1,3
      DO 110 JT=1,3
      DO 110 KT=1,3
      DO 110 LT=1,3
      DO 110 I=1,3
      DO 110 J=1,3
      DO 110 K=1,3
      DO 110 L=1,3
      U4(IT,JT,KT,LT) = U4(IT,JT,KT,LT) +
     & F(IT,I)*F(JT,J)*F(KT,K)*F(LT,L) * DRF(I,J,K,L)
  110 CONTINUE
C
C---- CALCULATE TRACE OF 4TH RANK TENSOR
C
C      DO 999 I=1,3
C      DO 999 J=1,3
C      DO 999 K=1,3
C      DO 999 L=1,3
C      WRITE(IOUT,99)I,J,K,L,U4(I,J,K,L)
C   99 FORMAT(4I4,F10.7)
C  999 CONTINUE
      U40 = 0.0
      DO 200 I=1,3
  200 U40 = U40 + U4(I,I,I,I)
      U40 = U40 + 2.*(U4(1,1,2,2)+U4(1,1,3,3)+U4(2,2,3,3))
      WRITE(IOUT,900) (ATOM(I,IATOM),I=1,2),U40
  900 FORMAT(' TRACE IN CART.FRAME FOR ATOM ',A4,A2,F9.6)      
C
C---- CALCULATE EIGENVALUES AND EIGENVECTORS OF U2
C
      WRITE(IOUT,910) U2
  910 FORMAT(' U2: ',3(F9.6,2X),/,2(5X,3(F9.6,2X),/))
	CALL JACOBI(U2,3,3,D,V,NROT)
	WRITE(IOUT,*) ' Unsorted Eigenvectors:'
	DO 11 I=1,3
	WRITE(IOUT,'(/1X,A,I3,A,F12.6)') 'Eigenvalue',I,' =',D(I)
		WRITE(IOUT,*) 'Eigenvector:'
		WRITE(IOUT,'(10X,5F12.6)') (V(J,I),J=1,3)
11	CONTINUE
	WRITE(IOUT,'(//,A,//)') '****** sorting ******'
	CALL EIGSRT(D,V,3,3)
	WRITE(IOUT,*) 'Sorted Eigenvectors:'
	DO 12 I=1,3
		WRITE(IOUT,'(/1X,A,I3,A,F12.6)') 'Eigenvalue',I,' =',D(I)
		WRITE(IOUT,*) 'Eigenvector:'
		WRITE(IOUT,'(10X,5F12.6)') (V(J,I),J=1,3)
12	CONTINUE
C
C---- REFILL MISSING TERMS
C
      U2(1,2) = U2(2,1)
      U2(1,3) = U2(3,1)
      U2(2,3) = U2(3,2)
C
C---- ROTATE CARTESIAN FRAME TO LINE UP WITH EIGENVECTORS
C
      CALL ROTAT(V,U2,U2E,U4,U4E)
C
C---- CALCULATE PARTIAL INVARIANT DELTA*U4E
C
      U4E11 = DBLE(U4E(1,1,1,1)+U4E(1,1,2,2)+U4E(1,1,3,3))
      U4E22 = DBLE(U4E(2,2,2,2)+U4E(2,2,3,3)+U4E(2,2,1,1))
      U4E33 = DBLE(U4E(3,3,3,3)+U4E(3,3,1,1)+U4E(3,3,2,2))
      U4E12 = DBLE(U4E(1,2,1,1)+U4E(1,2,2,2)+U4E(1,2,3,3))
      U4E13 = DBLE(U4E(1,3,1,1)+U4E(1,3,2,2)+U4E(1,3,3,3))
      U4E23 = DBLE(U4E(2,3,1,1)+U4E(2,3,2,2)+U4E(2,3,3,3))
      U4E11R = (DABS(U4E11/3.))**0.25
      U4E22R = (DABS(U4E22/3.))**0.25
      U4E33R = (DABS(U4E33/3.))**0.25
      WRITE(IOUT,920) U4E11,U4E22,U4E33,U4E12,U4E13,U4E23,
     &U4E11R,U4E22R,U4E33R
  920 FORMAT(/,' U4E-CONTRACTED : ',6F9.6, 'AE**4',/,
     &' U4E EQUI.RMSD  :',3F9.6, 'AE',/)
C
C---- CALCULATE EIGENVALUES AND EIGENVECTORS OF U2E
C
C	CALL JACOBI(U2E,3,3,D,V,NROT)
C	WRITE(IOUT,*) 'Unsorted Eigenvectors:'
C	DO 13 I=1,3
C	WRITE(IOUT,'(/1X,A,I3,A,F12.6)') 'Eigenvalue',I,' =',D(I)
C		WRITE(IOUT,*) 'Eigenvector:'
C		WRITE(IOUT,'(10X,5F12.6)') (V(J,I),J=1,3)
C13	CONTINUE
C	WRITE(IOUT,'(//,A,//)') '****** sorting ******'
C	CALL EIGSRT(D,V,3,3)
C	WRITE(IOUT,*) 'Sorted Eigenvectors:'
C	DO 14 I=1,3
C		WRITE(IOUT,'(/1X,A,I3,A,F12.6)') 'Eigenvalue',I,' =',D(I)
C		WRITE(IOUT,*) 'Eigenvector:'
C		WRITE(IOUT,'(10X,5F12.6)') (V(J,I),J=1,3)
C14	CONTINUE
      RETURN
      END
      SUBROUTINE ROTAT(A,U2,U2E,U4,U4E)
      DIMENSION U2(3,3),U2E(3,3),AA(3,3),A(3,3),
     & U4(3,3,3,3),U4E(3,3,3,3)
C
C---- TRANSFORM TENSOR COMPONENTS
C
      DO 80 I=1,3
      DO 80 J=1,3
      U2E(I,J) = 0.0
      DO 80 K=1,3
      DO 80 L=1,3
      U4E(I,J,K,L) = 0.0
   80 CONTINUE
      DO 90 IT=1,3
      DO 90 JT=1,3
      DO 90 I=1,3
      DO 90 J=1,3
      U2E(IT,JT) = U2E(IT,JT) + A(I,IT)*A(J,JT) * U2(I,J)
   90 CONTINUE
      DO 100 IT=1,3
      DO 100 JT=1,3
      DO 100 KT=1,3
      DO 100 LT=1,3
      DO 100 I=1,3
      DO 100 J=1,3
      DO 100 K=1,3
      DO 100 L=1,3
      U4E(IT,JT,KT,LT) = U4E(IT,JT,KT,LT) + A(I,IT)*A(J,JT)
     & *A(K,KT)*A(L,LT) * U4(I,J,K,L)
  100 CONTINUE
      RETURN
      END   
      SUBROUTINE JACOBI(A,N,NP,D,V,NROT)
      PARAMETER (NMAX=100)
      DIMENSION A(NP,NP),D(NP),V(NP,NP),B(NMAX),Z(NMAX)
      DO 12 IP=1,N
        DO 11 IQ=1,N
          V(IP,IQ)=0.
11      CONTINUE
        V(IP,IP)=1.
12    CONTINUE
      DO 13 IP=1,N
        B(IP)=A(IP,IP)
        D(IP)=B(IP)
        Z(IP)=0.
13    CONTINUE
      NROT=0
      DO 24 I=1,50
        SM=0.
        DO 15 IP=1,N-1
          DO 14 IQ=IP+1,N
            SM=SM+ABS(A(IP,IQ))
14        CONTINUE
15      CONTINUE
        IF(SM.EQ.0.)RETURN
        IF(I.LT.4)THEN
          TRESH=0.2*SM/N**2
        ELSE
          TRESH=0.
        ENDIF
        DO 22 IP=1,N-1
          DO 21 IQ=IP+1,N
            G=100.*ABS(A(IP,IQ))
            IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP)))
     *         .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN
              A(IP,IQ)=0.
            ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN
              H=D(IQ)-D(IP)
              IF(ABS(H)+G.EQ.ABS(H))THEN
                T=A(IP,IQ)/H
              ELSE
                THETA=0.5*H/A(IP,IQ)
                T=1./(ABS(THETA)+SQRT(1.+THETA**2))
                IF(THETA.LT.0.)T=-T
              ENDIF
              C=1./SQRT(1+T**2)
              S=T*C
              TAU=S/(1.+C)
              H=T*A(IP,IQ)
              Z(IP)=Z(IP)-H
              Z(IQ)=Z(IQ)+H
              D(IP)=D(IP)-H
              D(IQ)=D(IQ)+H
              A(IP,IQ)=0.
              DO 16 J=1,IP-1
                G=A(J,IP)
                H=A(J,IQ)
                A(J,IP)=G-S*(H+G*TAU)
                A(J,IQ)=H+S*(G-H*TAU)
16            CONTINUE
              DO 17 J=IP+1,IQ-1
                G=A(IP,J)
                H=A(J,IQ)
                A(IP,J)=G-S*(H+G*TAU)
                A(J,IQ)=H+S*(G-H*TAU)
17            CONTINUE
              DO 18 J=IQ+1,N
                G=A(IP,J)
                H=A(IQ,J)
                A(IP,J)=G-S*(H+G*TAU)
                A(IQ,J)=H+S*(G-H*TAU)
18            CONTINUE
              DO 19 J=1,N
                G=V(J,IP)
                H=V(J,IQ)
                V(J,IP)=G-S*(H+G*TAU)
                V(J,IQ)=H+S*(G-H*TAU)
19            CONTINUE
              NROT=NROT+1
            ENDIF
21        CONTINUE
22      CONTINUE
        DO 23 IP=1,N
          B(IP)=B(IP)+Z(IP)
          D(IP)=B(IP)
          Z(IP)=0.
23      CONTINUE
24    CONTINUE
      PAUSE '50 iterations should never happen'
      RETURN
      END
      SUBROUTINE EIGSRT(D,V,N,NP)
      DIMENSION D(NP),V(NP,NP)
      DO 13 I=1,N-1
        K=I
        P=D(I)
        DO 11 J=I+1,N
          IF(D(J).GE.P)THEN
            K=J
            P=D(J)
          ENDIF
11      CONTINUE
        IF(K.NE.I)THEN
          D(K)=D(I)
          D(I)=P
          DO 12 J=1,N
            P=V(J,I)
            V(J,I)=V(J,K)
            V(J,K)=P
12        CONTINUE
        ENDIF
13    CONTINUE
      RETURN
      END
      
 
