	program error88
C MAIN PROGRAM FOR PROGRAM BONDAN
C
      READ(5,1)NAME
 1    FORMAT(I2)
      IF(NAME-2)3,4,5
 3    CALL BADTEA
      GO TO 6
 4    CALL ERR
      GO TO 6
 5    WRITE(6,2)NAME
 2    FORMAT(' WRONG NUMBER ON CALLING CARD, ACTUAL NUMBER IS',I5)
 6    STOP
      END
C  PROGRAM BADTEA
      SUBROUTINE BADTEA
C THIS ROUTINE REVISED AND/OR RESEQUENCED ON OCT. 28, 1968
C
C        UNIVERSITY OF MINNESOTA   BOND ANGLES, BOND DISTANCES, THERMAL
C           ELLIPSOIDS AND ERRORS
C                       BY  LARRY W. FINGER         NOV,1965
C        CONVERTED TO FORTRAN IV     MAY 1967
C
C
      INCLUDE 'RFINE.PRM'
      DIMENSION S(3),IT(3),IPOS(2,3),FPB(8),IPB(10),IPA(10),
     1 IIH(8),SIG(QQPRM),JCELL(3,6),JPOS(3),JBETA(6),SIGG(QQPRM)
C
      COMMON/ABLE/XYZ(3,QQATM),BETA(6,QQATM),A(3),ANG(3),AS(3)
     1 ,CS(3),ATOM(2,QQATM),SIGA(6,6),SIGB(6,6,QQATM),RS(3,3,48)
     2 ,TS(3,48),TP(4,3),ITF(QQATM),NA,NSS,VOL,DLIMIT(3),PI,G(3,3)
     3 ,NLP,IERR,TITLE(18),CORR(1)
      COMMON/BAKER/XYZP(3,1000),IATP(20*QQATM)
C
      EQUIVALENCE (XYZP,SIG),(XYZP(181),SIGG)
C
      DATA FPB/0.0,0.5,2*0.0,0.25,0.75,2*0.0/
      DATA JPOS/1HX,1HY,1HZ/
      DATA IIH/3H   ,3H1/2,3H1/3,3H2/3,3H1/4,3H3/4,3H1/6,3H5/6/
      DATA IPB/3H   ,3H  X,3H +X,3H -X,3H  Y,3H +Y,3H -Y,3H  Z,3H +Z
     1,3H -Z/
      DATA IPA/0,1,1,-1,2,2,-2,3,3,-3/
      DATA IN,IOUT,INAUX,JCELL/5,6,4,1HA,2*1H ,1HB,2*1H ,1HC,2*1H ,4HCOS
     1 ,4HALPH,1HA,4HCOS ,4HBETA,1H ,4HCOS ,4HGAMM,1HA/
      DATA JBETA/2H11,2H22,2H33,2H12,2H13,2H23/
      DATA IP1,IA1,IB1,IC1,II1,IF1/1HP,1HA,1HB,1HC,1HI,1HF/
C
    1 FORMAT(18A4)
    2 FORMAT(1H110X,18A4)
    3 FORMAT(6F10.5)
    5 FORMAT(5I4,3XA1,6X,3F10.5)
    8 FORMAT(3(3A3,1X)I4)
    9 FORMAT(5X3(3A3,4X))
   34 FORMAT(10F8.3)
   50 FORMAT(8F10.8)
   56 FORMAT(2A3,7X3F7.5,I2/6F8.5)
   57 FORMAT(/24X,8A4,A1,9A4/10H REAL CELL,7X,6F12.5/
     1 17H RECIPROCAL CELL ,6F12.5//13H CELL VOLUME ,F16.5)
   58 FORMAT(//19H SYMMETRY OPERATORS,5X,A1,5H CELL,5X,A1,7HCENTRIC)
  106 FORMAT(/29H0STANDARD DEVIATIONS FOR CELL10X8A4,A1,9A4
     1    /31X,6F12.6//24H CELL CORRELATION MATRIX/)
  109 FORMAT(//50X,4HB OR/5X,4HATOM,13X,2(A4,6X),A1,6X,6(4HBETAA2,4X)/)
  110 FORMAT(4X,2A3,4X9F10.5)
  111 FORMAT(4X,2A3,4X3F10.5,F10.2)
  201 FORMAT(1X,3A4,2X,6F10.6)
  202 FORMAT(//44H TEMPERATURE FACTOR STANDARD DEVIATIONS FOR ,2A3/9X,
     16(4HBETA,A2,4X)/5X,6F10.6//19H CORRELATION MATRIX /)
  204 FORMAT(5X,4HBETA,A2,6F10.6)
  206 FORMAT(//34H ATOM POSITION STANDARD DEVIATIONS/10X,4HATOM,13X,3(A1
     1 ,9X)/)
  207 FORMAT(9X,2A3,4X,3F10.5)
  208 FORMAT(//5X,32HATOM POSITION CORRELATION MATRIX/)
  209 FORMAT(/4X,2A3,2X,1HX,7X,10F10.5)
  210 FORMAT(/12XA1,7X10F10.5)
  211 FORMAT(20X10F10.5)
C
      FPB(3)=1./3.
      FPB(4)=2./3.
      FPB(7)=1./6.
      FPB(8)=5./6.
      PI=ACOS(-1.0)
      READ(IN,1)TITLE
      WRITE(IOUT,2)TITLE
      READ(IN,3)A,ANG
C        CALCULATE SINES AND COSINES OF DIRECT SPACE ANGLES
      DO 4 I=1,3
      IF(ABS(ANG(I)-90.0).GT.1.0E-3)GO TO 1000
      S(I)=1.0
      ANG(I)=0.0
      GO TO 4
 1000 ANG(I)=ANG(I)*PI/180.0
      S(I)=SIN(ANG(I))
      ANG(I)=COS(ANG(I))
C        MATRIX G IS THE REAL SPACE METRIC TENSOR
    4 G(I,I)=A(I)**2
      DO 144 I=1,2
      IJ=I+1
      DO 144 J=IJ,3
      IK=6-I-J
      G(I,J)=A(I)*A(J)*ANG(IK)
  144 G(J,I)=G(I,J)
C        COMPUTE REAL SPACE CELL VOLUME
      VOL=A(1)*A(2)*A(3)*SQRT (1.0-ANG(1)**2-ANG(2)**2-ANG(3)**2
     1 +2.0*ANG(1)*ANG(2)*ANG(3))
C        COMPUTE RECIPROCAL SPACE CELL PARAMETERS
      AS(1)=A(2)*A(3)*S(1)/VOL
      AS(2)=A(1)*A(3)*S(2)/VOL
      AS(3)=A(1)*A(2)*S(3)/VOL
      CS(1)=(ANG(2)*ANG(3)-ANG(1))/(S(2)*S(3))
      CS(2)=(ANG(1)*ANG(3)-ANG(2))/(S(1)*S(3))
      CS(3)=(ANG(1)*ANG(2)-ANG(3))/(S(1)*S(2))
      WRITE(IOUT,57)JCELL,A,ANG,AS,CS,VOL
C        READ CONTROL CARD
      READ(IN,5)NA,ICNT,IBOD,IELV,IERR,LTYPE,DLIMIT
      IF(DLIMIT(3).LE.0.0) DLIMIT(3)=DLIMIT(2)
      I=2-ICNT
      WRITE(IOUT,58)LTYPE,JCELL(I,1)
C        READ SYMMETRY CARDS
      DO 7 I=1,48
      READ(IN,8)(IT(K),(IPOS(L,K),L=1,2),K=1,3),IEF
      WRITE(IOUT,9)(IT(K),(IPOS(L,K),L=1,2),K=1,3)
C        CLEAR ROTATIONAL MATRICES
      DO 12 J=1,3
      DO 10 K=1,3
   10 RS(K,J,I)=0
C        CONVERT TRANSLATIONAL OPERATOR
      DO 11 K=1,8
      IF(IT(J)-IIH(K))11,13,11
   11 CONTINUE
C        CARD CHARACTERS ARE NOT IN THE ACCEPTABLE LIST
      STOP 11
C        SET TRANSLATIONAL VALUE
   13 TS(J,I)=FPB(K)
      DO 14 K=1,3
      DO 14 L=1,2
C        COMPARE POSITIONAL INFORMATION WITH ACCEPTABLE VALUES
      DO 15 N=1,10
      IF(IPB(N)-IPOS(L,K))15,16,15
   15 CONTINUE
C        IMPROPER CODE ON CARD
      STOP 11
   16 IF(J-IABS (IPA(N)))14,17,14
C        STORE ROTATIONAL MATRIX ELEMENT
   17 RS(K,J,I)=ISIGN (1,IPA(N))
   14 CONTINUE
   12 CONTINUE
C        CHECK IF THIS IS THE LAST SYMMETRY CARD
      IF(IEF)18,7,18
    7 CONTINUE
C        TOO MANY SYMMETRY CARDS
      STOP 22
   18 NSS=I
      IF(ICNT)38,39,38
C        EXPAND LIST FOR CENTRIC SPACE GROUP
   39 DO 40 I=1,NSS
      N=I+NSS
      DO 40 J=1,3
      TS(J,N)=-TS(J,I)
      DO 40 K=1,3
   40 RS(K,J,N)=-RS(K,J,I)
      NSS=2*NSS
C        FORM MATRICES FOR LATTICE POINTS
   38 NLP=1
      DO 19 I=1,3
   19 TP(1,I)=0
      IF(LTYPE-IP1)21,20,21
C        LATTICE IS NOT PRIMITIVE
   21 DO 22 I=1,3
   22 TP(2,I)=0.5
      NLP=2
      IF(LTYPE-IA1)23,24,23
C        A CENTERED LATTICE
   24 TP(2,1)=0
      GO TO 20
   23 IF(LTYPE-IB1)25,26,25
C        B CENTERED LATTICE
   26 TP(2,2)=0
      GO TO 20
   25 IF(LTYPE-IC1)28,27,28
C        C CENTERED LATTICE
   27 TP(2,3)=0
      GO TO 20
   28 IF(LTYPE-II1)29,20,29
   29 IF(LTYPE-IF1)31,30,31
C        F CENTERED LATTICE
   30 NLP=4
      DO 32 I=2,4
      DO 33 J=1,3
   33 TP(I,J)=0.5
   32 TP(I,I-1)=0
      GO TO 20
C        R CENTERED LATTICE
   31 NLP=3
      TP(2,1)=1./3.
      TP(2,2)=2./3.
      TP(2,3)=TP(2,2)
      TP(3,1)=TP(2,2)
      TP(3,2)=TP(2,1)
      TP(3,3)=TP(2,1)
   20 IF(IERR.EQ.0)GO TO 101
C        READ CELL SIGMAS AND CORRELATION MATRIX, CHANGE TO VAR-COV
  100 READ(IN,50)(SIG(I),I=1,6)
C        CONVERT SIGMA OF ANGLES TO SIGMA OF COSINE
      DO 102 I=4,6
  102 SIG(I)=S(I-3)*SIG(I)*PI/180.0
      WRITE(IOUT,106)JCELL,(SIG(I),I=1,6)
      DO 200 I=1,6
      READ(IN,34)(SIGA(I,J),J=1,I)
      SIGA(I,I)=1.0
      WRITE(IOUT,201)(JCELL(J,I),J=1,3),(SIGA(I,J),J=1,I)
      DO 200 J=1,I
      SIGA(I,J)=SIGA(I,J)*SIG(I)*SIG(J)
  200 SIGA(J,I)=SIGA(I,J)
  101 WRITE(IOUT,109)JPOS,JBETA
      IO=IN
      IF(IERR.NE.2) GO TO 300
C        SET TO READ AUXILLARY INPUT UNIT
      IO=INAUX
      REWIND IO
C        READ STRUCTURE PARAMETERS
  300 DO 6 I=1,NA
      READ(IO,56)(ATOM(J,I),J=1,2),(XYZ(J,I),J=1,3),ITF(I),(BETA(J,I),J=
     11,6)
      IF(ITF(I))112,113,112
  112 WRITE(IOUT,110)(ATOM(J,I),J=1,2),(XYZ(J,I),J=1,3),(BETA(J,I),J=1,6
     1)
      GO TO 6
  113 WRITE(IOUT,111)(ATOM(J,I),J=1,2),(XYZ(J,I),J=1,3),BETA(1,I)
    6 CONTINUE
      IF(IERR.EQ.0.OR.(IELV.EQ.0.AND.IERR.NE.2)) GO TO 103
C        READ TEMPERATURE FACTOR SIGMAS AND CORR. MATRIX
      I=0
      DO 41 II=1,NA
      I=I+1
      IF(ITF(I))42,41,42
   42 READ(IO,50)(SIG(J),J=1,6)
      WRITE(IOUT,202)(ATOM(J,I),J=1,2),JBETA,(SIG(J),J=1,6)
      DO 203 J=1,6
      READ(IO,34)(SIGB(J,K,I),K=1,J)
      SIGB(J,J,I)=1.0
      WRITE(IOUT,204)JBETA(J),(SIGB(J,K,I),K=1,J)
C        CONVERT CORRELATION MATRIX TO VARIANCE-COVARIANCE MATRIX
      DO 203 K=1,J
      SIGB(J,K,I)=SIGB(J,K,I)*SIG(J)*SIG(K)
  203 SIGB(K,J,I)=SIGB(J,K,I)
   41 CONTINUE
  103 IF(IERR.EQ.0.OR.(IBOD.EQ.0.AND.IERR.NE.2)) GO TO 105
C        READ ATOM POSITION SIGMAS AND CORRELATION MATRIX
  104 WRITE(IOUT,206)JPOS
      NAA=3*NA
      DO 205 I=1,NA
      IJ=3*(I-1)+1
      JI=IJ+2
      READ(IO,50)(SIG(J),J=IJ,JI)
  205 WRITE(IOUT,207)(ATOM(J,I),J=1,2),(SIG(J),J=IJ,JI)
      WRITE(IOUT,208)
      L=0
      DO 223 I=1,NA
      DO 223 J=1,3
      IJ=3*(I-1)+J
      READ(IO,34)(SIGG(K),K=1,IJ)
      SIGG(IJ)=1.0
      IF(J.NE.1.OR.I.GT.4)GO TO 215
      WRITE(IOUT,209)(ATOM(K,I),K=1,2),(SIGG(K),K=1,IJ)
      GO TO 216
  215 IF(J.NE.1)GO TO 218
      WRITE(IOUT,209)(ATOM(K,I),K=1,2),(SIGG(K),K=1,10)
      GO TO 221
  218 IF(I.GE.4)GO TO 220
      WRITE(IOUT,210)JPOS(J),(SIGG(K),K=1,IJ)
      GO TO 216
  220 WRITE(IOUT,210)JPOS(J),(SIGG(K),K=1,10)
  221 WRITE(IOUT,211)(SIGG(K),K=11,IJ)
C        CONVERT CORRELATION MATRIX TO VARIANCE-COVARIANCE MATRIX
  216 DO 222 K=1,IJ
  222 SIGG(K)=SIGG(K)*SIG(K)*SIG(IJ)
      SIGG(IJ)=SIG(IJ)**2
      DO 223 K=1,IJ
      L=L+1
  223 CORR(L)=SIGG(K)
  105 IF(IBOD.EQ.0)GO TO 44
C        COMPUTE BOND DISTANCES AND ANGLES
   43 CALL BODAN(IBOD)
   44 IF(IELV)45,46,45
C        COMPUTE THERMAL ELLIPSOIDS
   45 CALL ELVIB
C        END OF CALCULATION, TERMINATE PROGRAM
   46 CALL EXIT
      END
      SUBROUTINE ERR
C        PROGRAM ERROR TO READ THE FILE OF RFINE88 MODIFIED FOR GENERAL
C        CONSTRAINTS.     MODIFIED JAN. 18, 1974
C
C   PROGRAM TO READ THE DATA FILE PREPARED BY THE FINGER LEAST-SQUARES
C        PROGRAM, REFORMAT THE PARAMETER VARIANCE-COVARIANCE MATRIX
C        AND THEN CALL THE BOND DISTANCE AND ANGLE SUBROUTINE OR THE
C        THERMAL ELLIPSOID SUBROUTINE TO COMPUTE THESE FUNCTIONS
C        AND THEIR ERRORS.
C
C
C***********************************************************************
C
C                             NOTE TO USER
C
C     DIMENSION OF ARRAY SIGX MUST BE AT LEAST
C     29*NA + (MRANK*(MRANK+1))/2 + 2*NDEGN + (3*NA*(3*NA+1))/2
C     MRANK = NUMBER OF REFINED PARAMETERS
C     NDEGN = NUMBER OF CONSTRAINT RELATIONS
C     NA = NUMBER OF ATOMS IN ASYMMETRIC UNIT
C
C***********************************************************************
C
      INCLUDE 'RFINE.PRM'
      PARAMETER (QQS=(3*qqatm*(3*qqatm+1))/2+(QQPRM*(QQPRM+1)/2)
     1+29*qqatm,QQA=20*QQATM)
      COMMON/ABLE/XYZ(3,QQATM),BETA(6,QQATM),A(3),ANG(3),AS(3)
     1 ,CS(3),ATOM(2,QQATM),SIGA(6,6),SIGB(6,6,QQATM),RS(3,3,48)
     2 ,TS(3,48),TP(4,3),ITF(QQATM),NA,NSS,VOL,DLIMIT(3),PI,G(3,3)
     3 ,NLP,IERR,TITLE(18),SIGX(QQS)
      DIMENSION IPARA(QQA),KPAR(QQATM),KRPAR(QQATM),KLAB(3)
     1 ,LABEL(3),S(3),TEMP(20),TEMP1(9,20),TEMP2(2,20),TENP(3)
     2 ,TEOP(5),ZT(5,QQATM),JSF(10),CR(10,QQATM),DR(15,QQATM)
     3 ,PAR(QQATM),PARNAM(2,QQATM),TEMP4(6),ITEMP4(7),NCNT(4)
     4 ,EXTIN(7),IPNT(6,QQATM),JPNT(QQPRM),KNSTR(2),CV(2),CNSDGN(2)
      EQUIVALENCE(SIGX,IPARA),(SIGX(QQA+1),KRPAR)
     1 ,(SIGX(QQA+QQATM+1),IPNT,JPNT),
     2 (SIGX(QQA+7*QQATM+1),PAR),(SIGX(QQA+8*QQATM+1),KPAR)
     3 ,(SIGX(QQA+9*QQATM+1),CR,DR,PARNAM,
     4 TEMP,TEMP1,TEMP2,TENP,TEOP,ZT,TEMP4,KNSTR,CNSDGN,CV)
C
      DATA IN,IOUT,ITAPE,NCNT/5,6,10,1,6,16,31/
C
C        DEFINE TWO INDEXING FUNCTIONS FOR TRIANGULAR MATRICES
      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
     1+MAX0(I,J)+1-MIN0(I,J)+MDEGN
c     I1(K)=IAND(K/256**3,255)
      I2(K)=IAND(K/256**2,255)
      I3(K)=IAND(K/256,255)
      I4(K)=IAND(K,255)
C
    2 FORMAT(18A4/I2,3A4)
    4 FORMAT(2I1,3F8.0)
    6 FORMAT(6E12.6)
    8 FORMAT(' ACTUAL LABEL ON THE FILE WAS '3A4,' THE SEARCH LABEL WAS
     $'3A4/'0JOB TERMINATED.')
   10 FORMAT('0',18A4/'0CELL VARIANCES'/)
   12 FORMAT(6(2X,E12.5))
   14 FORMAT(' ONLY ',I6,' MEMORY LOCATIONS FOR SIGX MATRIX, '
     1 ,I6,' NEEDED.')
C
      READ(IN,2)TITLE,KEY,LABEL
      PI= ACOS(-1.0)
C        IF KEY=0, A DATA SET SEARCH WILL BE MADE
  100 READ(ITAPE)KCYL,KLAB,AS,CS,G,VOL,NSS,NA,NSCAT,NLP,MRANK,NPARA,NSF,
     1TENP,ICENT,(((RS(I,J,K),I=1,3),J=1,3),K=1,NSS),((TS(I,J),I=1,3),J=
     21,NSS),((TP(I,J),J=1,3),I=1,NLP),(IPARA(I),I=1,NPARA),((TEMP1(I,J)
     3,I=1,9),J=1,NSCAT),(TEMP(I),I=1,NSF),TEOP,(JSF(I),I=1,NSF),NRSF
     4,NPAR,(PAR(I),(PARNAM(J,I),J=1,2),KPAR(I),KRPAR(I),I=1,NPAR),FISGN
     5,IFSGN,TEMP4,ITEMP4
      READ(ITAPE)((TEMP2(I,J),I=1,2),J=1,NSCAT),IEXT,EXTIN,((ATOM(J,I),J
     1=1,2),ITF(I),(ZT(J,I),J=1,5),(XYZ(J,I),J=1,3),(BETA(J,I),J=1,6),
     2(CR(J,I),J=1,10),(DR(J,I),J=1,15),I=1,NA),MDEGN
     3,(KNSTR(I),I=1,MDEGN)
      IF(NPAR.EQ.0)GO TO 110
C     FOR RIGID BODY REFINEMENTS APPLY LIBRATION CORRECTIONS TO XYZS.
      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 READ(ITAPE)ITEMP
      IF(ITEMP)120,110,120
C        READ VARIANCE-COVARIANCE MATRIX
  120 READ(ITAPE)MDEGN,(KNSTR(I),I=1,MDEGN)
      NMXL=(MRANK*(MRANK+1))/2+MDEGN
      J=MDEGN+1
      READ(ITAPE)(CV(I),I=J,NMXL)
C        CHECK IF DESIRED DATA SET READ
      IF(KEY)140,130,140
  130 DO 135 I=1,3
      IF(KLAB(I)-LABEL(I))136,135,136
  135 CONTINUE
      GO TO 140
C        INCORRECT OR UNEXPECTED LABEL
  136 WRITE(IOUT,8)KLAB,LABEL
      CALL EXIT
  140 READ(IN,4)MBODAN,MELVIB,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        SKIP NEXT SECTION IF NO ELLIPSOIDS CALCULATED
142   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
      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(KNSTR(I))
      IF(IPAR.LT.5.OR.IPAR.GT.10)GO TO 235
C        DEP. PARAMETER IS AN ANISO. COEF.
      IPAR=IPAR-4
      IAT=I3(KNSTR(I))
      J=I2(KNSTR(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)+CNSDGN(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(KNSTR(KK)).NE.IAT.OR.I4(KNSTR(KK)).LT.5.OR.I4(KNSTR(KK)).GT.
     110)GO TO 1500
      IJPAR=I4(KNSTR(KK))-4
      IF(IPAR.GT.IJPAR)GO TO 1500
      L=I2(KNSTR(KK))
      IK=INDEXU(J,L)
 1200 SIGB(IPAR,IJPAR,IAT)=SIGB(IPAR,IJPAR,IAT)+CNSDGN(I+1)*
     1CNSDGN(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
      i=1
  236 IF(IPARA(I).EQ.0)THEN
        I=I+1
        GO TO 236
      ELSE
        IPAR=IPARA(I)
      END IF
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
C        CHECK IF SUFFICIENT MEMORY FOR XYZ VAR-COV, THEN CLEAR
      if(i.gt.qqs-29*qqatm)then
        write(iout,14)i+29*qqatm,qqs
        stop
      endif
      DO 240 J=KB,I
  240 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
      KPNT=IPARA(IPAR)
      JPNT(NP)=KPNT
C        TEST IF POSITION REFINED
      IF(KPNT.EQ.0)GO TO 480
      KA=INDEXU(KPNT,KPNT)
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
      IF(MDEGN.EQ.0) GO TO 570
      NB=3*NA
      DO 565 I=1,MDEGN,2
      IPAR=I4(KNSTR(I))
      IF(IPAR.LT.2.OR.IPAR.GT.4)GO TO 565
      IAT=I3(KNSTR(I))
      J=I2(KNSTR(I))
      IPAR=3*(IAT-1)+IPAR-1
      DO 560 K=1,NB
      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     COMPUTE COVARIANCES OF UNREFINED POSITIONS WITH EACH OTHER.
       DO 563 KK=1,MDEGN,2
      IJPAR=I4(KNSTR(KK))
      IF(IJPAR.LT.2.OR.IJPAR.GT.4)GO TO 563
      JAT=I3(KNSTR(KK))
      L=I2(KNSTR(KK))
      IJPAR=3*(JAT-1)+IJPAR-1
      IF(IPAR.GT.IJPAR)GO TO 563
      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
  790 CONTINUE
C        READ AND CONVERT CELL ERROR MATRIX
      READ(IN,6)((SIGA(J,I),J=I,6),I=1,6)
      WRITE(IOUT,10)TITLE
      DO 795 I=1,6
 795  WRITE(IOUT,12)(SIGA(I,J),J=1,I)
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        EXPAND SYMMETRY OPERATORS FOR A CENTRIC CRYSTAL
  860 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)
  900 IF(MELVIB)910,920,910
C        COMPUTE THERMAL FUNCTIONS
  910 CALL ELVIB
C        END OF JOB
  920 CALL EXIT
      END
      SUBROUTINE BODAN(IBOD)
C
C THIS ROUTINE REVISED AND/OR RESEQUENCED ON OCT. 28, 1968
C
C        BOND DISTANCE AND ANGLE SUBPROGRAM WITH ERRORS
C
      INCLUDE 'RFINE.PRM'
      COMMON/ABLE/XYZ(3,QQATM),BETA(6,QQATM),A(3),ANG(3),AS(3)
     1 ,CS(3),ATOM(2,QQATM),SIGA(6,6),SIGB(6,6,QQATM),RS(3,3,48)
     2 ,TS(3,48),TP(4,3),ITF(QQATM),NA,NSS,VOL,DLIMIT(3),PI,G(3,3)
     3 ,NLP,IERR,TITLE(18),CORR(1)
      COMMON/BAKER/XYZP(3,1000),IATP(20*QQATM)
C
      DIMENSION XP(3),XYP(3),XXP(3),SIGX(6,6),DLDX(6),DSDX(9,2),YP(3),
     1 DLDA(6,2),SIGY(9,9),DNDA(6),DNDX(9),DTHDA(6),DTHDX(9),XZP(3),
     2 X(3,3),NXX(6),DXX(3),DXY(3),DXZ(3),DBD(3)
C
      EQUIVALENCE (DLDX,DSDX),(NXX(1),NXL),(NXX(2),NYL),(NXX(3),NZL)
C
      DATA IOUT/6/
C
 1002 FORMAT(46X,2A3,4F10.5,10X,I8)
 8999 FORMAT(31H THE ATOM ARRAY IS OVERFLOWING. ///15X,54HREDUCE THE LIM
     1ITS OF THE BOND DISTANCES AND TRY AGAIN. )
 9000 FORMAT(1H1,18A4/25H  BOND DISTANCES BETWEEN ,F5.2,5H AND ,F5.2,10H
     1 ANGSTROMS//10X,6HATOM 1,7X1HX,9X,1HY,9X,1HZ,2X,6HATOM 2,7X,1HX,9X
     2,1HY,9X,1HZ,4X,8HDISTANCE,3X5HERROR,4X,6HNUMBER//10X2A3,3F10.5,63X
     3,1H1)
 9001 FORMAT(10X2A3,3F10.5,56XI8)
 9002 FORMAT(1H1,18A4//69H BOND ANGLES AND ERRORS FOR ATOMS WITH INTERAT
     1OMIC DISTANCES BETWEEN F6.2,5H AND F6.2,11H ANGSTROMS./102H  THE A
     2NGLE IS ATOM 1-CENTRAL ATOM-ATOM 2.  THE NUMBER POINTS TO AN ENTRY
     3 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
     1ER,10X6HATOM 2,2X,6HNUMBER,10X,5HANGLE,5X,5HERROR/)
C
      WRITE(IOUT,9000)TITLE,DLIMIT(1),DLIMIT(3),(ATOM(J,1),J=1,2),(XYZ(J
     1,1),J=1,3)
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 II=1,NA
C        SET UP SEARCH ARRAYS FOR THIS ATOM
      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
      SIGX(N,1)=CORR(L)
      IF(N.EQ.1)GO TO 3
      DO 2 J=2,N
      L=L+1
    2 SIGX(N,J)=CORR(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
      SIGX(K+3,4)=CORR(L)
      IF(K.EQ.1)GO TO 77
      DO 78 N=2,K
      L=L+1
   78 SIGX(K+3,N+3)=CORR(L)
   77 JI=JI+IJ+K-1
C        GET CROSS TERMS FOR VARIANCE-COVARIANCE MATRIX
      CALL CORR1(SIGX,6,II,IATP(J)/64,CORR,1,4,NA)
C        COMPUTE VARIANCE OF DISTANCE
      CALL VARIAN(DS,DLDX,DLDA,SIGX,SIGA,1)
      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
C
C        LOOP THROUGH THE ASYMMETRIC UNIT TO GENERATE ATOMS
   67 DO 1 JJ= 1,NA
C        LOOP THROUGH THE SYMMETRY OPERATIONS
      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)=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)
      RETURN
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
      SIGX(K+3,4)=CORR(L)
      IF(K.EQ.1)GO TO 23
      DO 22 N=2,K
      L=L+1
   22 SIGX(K+3,N+3)=CORR(L)
   23 JI=JI+IJ+K-1
C        GENERATE CROSS MATRIX ELEMENTS
      CALL CORR1(SIGX,6,II,JJ,CORR,1,4,NA)
      S=SQRT (S)
      CALL VARIAN(DS,DLDX,DLDA,SIGX,SIGA,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 I=1,NA
      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)=CORR(L)
      IF(N.EQ.1)GO TO 124
      DO 25 K=2,N
      L=L+1
   25 SIGY(N+3,K+3)=CORR(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)=CORR(L)
      IF(K.EQ.1)GO TO 26
      DO 27 N=2,K
      L=L+1
   27 SIGY(K,N)=CORR(L)
   26 JI=JI+IJ+K-1
      CALL CORR1(SIGY,9,I,NAT,CORR,1,4,NA)
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)=CORR(L)
      IF(M.EQ.1)GO TO 50
      DO 49 N=2,M
      L=L+1
   49 SIGY(M+6,N+6)=CORR(L)
   50 JI=JI+IJ+M-1
      CALL CORR1(SIGY,9,I,NTA,CORR,4,7,NA)
      CALL CORR1(SIGY,9,NAT,NTA,CORR,1,7,NA)
      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))+
     1 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)*
     1 (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,DT
     1H
   40 CONTINUE
   29 CONTINUE
   28 CONTINUE
      RETURN
      END
      SUBROUTINE ELVIB
C
C THIS ROUTINE REVISED AND/OR RESEQUENCED ON OCT. 28, 1968
C
C        THIS ROUTINE COMPUTES ELLIPSOIDS OF THERMAL MOTION AND THE
C         ERRORS ASSOCIATED WITH THE PARAMETERS OF THESE ELLIPSOIDS
C
      INCLUDE 'RFINE.PRM'
      COMMON/ABLE/XYZ(3,QQATM),BETA(6,QQATM),A(3),ANG(3),AS(3)
     1 ,CS(3),ATOM(2,QQATM),SIGA(6,6),SIGB(6,6,QQATM),RS(3,3,48)
     2 ,TS(3,48),TP(4,3),ITF(QQATM),NA,NSS,VOL,DLIMIT(3),PI,G(3,3)
     3 ,NLP,IERR,TITLE(18),CORR(1)
C
      DIMENSION DDDC(6,3,3),DDDB(6,3,3),DTDA(6),DTDB(6),DSDA(6),DSDB(6),
     1 DDETDA(6),DDETDB(6),DLDA(6),DLDB(6),DQDA(6,3),DQDB(6,3),
     2 DBQDA(6),B(3,3),D(3,3),ROOT(3),DBQDB(6),DASDA(6,6)
     3 ,THETA(3),STH(3),DTHDB(6,3),DTHDA(6,3),DXDA(6,3,3),DXDB(6,3,3),
     4 X(3,3),IAXIS(3),QI(3)
C
      DATA IAXIS,IOUT/1HA,1HB,1HC,6/
C
 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)
 9002 FORMAT(///10X2A3,25H IS NON POSITIVE DEFINITE)
 9003 FORMAT(///5X,27HELLIPSOID OF VIBRATION FOR 2A3//20X,9HPARAMETER,15
     1X,5HVALUE,10X,5HERROR//13X22HEQUIVALENT ISOTROPIC B,F14.3,F15.3)
 9004 FORMAT(/10X,5HAXIS ,I1/18X,13HRMS AMPLITUDE,F18.6,F15.6,10H ANGSTR
     1OMS )
 9005 FORMAT(14X,15HANGLE WITH THE ,A1,5H AXIS,F14.2,F15.2,8H DEGREES)
C
      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
 9001 FORMAT(1H1,18A4//22H ELLIPSOIDS AND ERRORS)
C        MATRIX D IS THE PRODUCT GB
C        B IS THE MATRIX OF ANISOTROPIC COEFFICIENTS
C        G IS THE REAL SPACE METRIC TENSOR
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)
      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
      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 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)-
     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
      IF(ABS(DTDB(2)).LT.1.0E-10.AND.ABS(DTDB(3)).LT.1.0E-10)GO TO 9990
C        CHECK IF ATOM POSITIVE-DEFINITE
      IF(0.25*DTDB(2)**2-DTDB(3)-1.0E-5)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
      GO TO 9991
C THREE EQUAL ROOTS
9990  ROOT(1)=DET**(1.0/3.0)
      ROOT(2)=ROOT(1)
      ROOT(3)=ROOT(1)
C        PLACE ROOTS IN INCREASING ORDER
9991  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
      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
      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(
     1L,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(
     1L,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
     1*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))
     1 *CS(JX2)
      DBQDA(K)=DBQDA(K)+AS(J)*AS(L)*(QI(J)*QI(L)*DASDA(JX2+3,K)+CS(JX2)*
     1(QI(J)*DQDA(K,L)+QI(L)*DQDA(K,J)))+QI(J)*QI(L)*CS(JX2)*(AS(J)*
     2 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
      SUBROUTINE CORR1(SIGZ,NN,II,JJ,CORR,KK,LL,MM)
C
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)
C
      DIMENSION SIGZ(NN,NN),CORR(7260)
C
      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
      SUBROUTINE CELLD(DASDA,A,ANG,AS,CS,VOL)
C
      DIMENSION A(3),ANG(3),AS(3),CS(3),DASDA(6,6),SN(3)
C
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.
      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
      SUBROUTINE DERDET(A,B,C)
C
C        THIS ROUTINE COMPUTES THE DERIVATIVE OF A DETERMINANT
C
      DIMENSION A(6),B(3,3),C(6,3,3)
C
      DO 4 I=1,6
      A(I)=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))
     1  +C(I,J,N)*(B(M,J)*B(N,M)-B(M,M)*B(N,J))
     2  +C(I,J,M)*(B(N,J)*B(M,N)-B(N,N)*B(M,J))
      RETURN
      END
      SUBROUTINE VARIAN(A,B,C,D,E,II)
C
      DIMENSION B(6),C(6),D(6,6)   ,E(6,6,40)
C
      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
    3 A=SQRT (F)
      RETURN
      END
      SUBROUTINE LIBCOR(XYZ,PAR,KPAR,NPAR,G,NA,MDEGN,KNSTR)
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.
      DIMENSION XYZ(3,60),PAR(60),KPAR(60),G(3,3),KNSTR(MDEGN),ATOX(6)
     1,XTOA(6),RC(3),R(3),EL(3,3),OXYZ(3),SXYZ(3),L(6),OXP(3)
     2,TRNS(3,3),ELTEM(3,3)
      EQUIVALENCE(OMEGA,OXP(1)),(CHI,OXP(2)),(PHI,OXP(3))
      I1(K)=IAND(K/256**3,255)
      I2(K)=IAND(K/256**2,255)
      I3(K)=IAND(K/256,255)
      I4(K)=IAND(K,255)
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)
     1 .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.
      OXP(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.4.AND.IDEN.LE.6)GO TO 122
      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
  122 OXP(IDEN-3)=PAR(J)
      GO TO 130
  124 OXYZ(IDEN)=PAR(J)
      GO TO 130
  126 IDEN=IDEN-15
      IF(IDEN.GT.3)GO TO 128
      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
      IF(OMEGA.EQ.0..AND.CHI.EQ.0..AND.PHI.EQ.0.)GO TO 139
C     THE COORDINATE SYSTEM TO WHICH L IS REFERRED IS NOT THE STANDARD
C     ORTHONORMAL SYSTEM. L MUST BE TRANSFORMED TO THE STANDARD SYSTEM.
      CALL EULER(PHI,CHI,OMEGA,TRNS)
      DO 131 J=1,3
      DO 131 K=1,3
      ELTEM(J,K)=EL(J,K)
  131 EL(J,K)=0.
      DO 135 J=1,3
      DO 135 K=1,3
      DO 135 LL=1,3
      DO 135 M=1,3
  135 EL(J,K)=EL(J,K)+TRNS(J,LL)*TRNS(K,M)*ELTEM(LL,M)
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.
  139 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)GO TO 140
      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
      SUBROUTINE UTMUL(A,B,N,C)
C     MULTIPLIES VECTOR B BY UPPER TRIANGULAR MATRIX A, STORING THE
C     RESULT IN C. A IS STORED BY ROWS.
      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
      SUBROUTINE CORLIB(R3,ELT,RC)
      DIMENSION R3(3,3),ELT(3,3),RC(3),ACON(3,3),BCON(3)
      DO 30 I=1,3 
      BCON(I)=RC(I) 
      DO 30 J=1,3 
30    ACON(I,J)=0.
40    DO 50 I=1,3 
      DO 50 J=1,3 
50    ACON(I,J)=ACON(I,J)+EIJ(ELT,I,J)
      CALL SOLVE(ACON,BCON,R3)
      RETURN
      END 
      FUNCTION DETERM(C1,C2,C3) 
      DIMENSION C1(3),C2(3),C3(3) 
      DETERM=C1(1)*C2(2)*C3(3)+C1(2)*C2(3)*C3(1)+C1(3)*C2(1)*C3(2)
     1 -C1(3)*C2(2)*C3(1)-C1(2)*C2(1)*C3(3)-C1(1)*C2(3)*C3(2) 
      RETURN
      END 
      FUNCTION EIJ(EL,I,J)
      DIMENSION EL(3,3) 
      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,
     1 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, 
     1 L,L)+3.*(EL(K,K)**2+EL(L,L)**2)+2.*EL4(EL,K,K,L,L))
      RETURN
      END 
      FUNCTION EL4(EL,I,J,K,L)
      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 EULER(PHI,CHI,OMEGA,TRANS) 
C     SUBROUTINE TO GENERATE THE MATRIX B-INVERSE FROM THE EULERIAN 
C     ANGLES PHI,CHI,AND OMEGA. 
      DATA RAD/57.2957792/
      DIMENSION TRANS(3,3)
      ASSIGN 10 TO KGO
      GO TO 100 
10    TRANS(1,1)=CP*CX*CO-SP*SO 
      TRANS(2,1)=CP*CX*SO+SP*CO 
      TRANS(3,1)=-CP*SX 
      TRANS(1,2)=-SP*CX*CO-CP*SO
      TRANS(2,2)=-SP*CX*SO+CP*CO
      TRANS(3,2)=SP*SX
      TRANS(1,3)=SX*CO
      TRANS(2,3)=SX*SO
      TRANS(3,3)=CX 
      RETURN
      ENTRY DTDPHI(PHI,CHI,OMEGA,TRANS) 
      ASSIGN 20 TO KGO
      GO TO 100 
20    TRANS(1,1)=(-SP*CX*CO-CP*SO)/RAD
      TRANS(2,1)=(-SP*CX*SO+CP*CO)/RAD
      TRANS(3,1)=SP*SX/RAD
      TRANS(1,2)=(-CP*CX*CO+SP*SO)/RAD
      TRANS(2,2)=(-CP*CX*SO-SP*CO)/RAD
      TRANS(3,2)=CP*SX/RAD
      DO 25 I=1,3 
25    TRANS(I,3)=0. 
      RETURN
      ENTRY DTDCHI(PHI,CHI,OMEGA,TRANS) 
      ASSIGN 30 TO KGO
      GO TO 100 
30    TRANS(1,1)=-CP*SX*CO/RAD
      TRANS(2,1)=-CP*SX*SO/RAD
      TRANS(3,1)=-CP*CX/RAD
      TRANS(1,2)=SP*SX*CO/RAD
      TRANS(2,2)=SP*SX*SO/RAD
      TRANS(3,2)=SP*CX/RAD
      TRANS(1,3)=CX*CO/RAD
      TRANS(2,3)=CX*SO/RAD
      TRANS(3,3)=-SX/RAD
      RETURN
      ENTRY DTDOMG(PHI,CHI,OMEGA,TRANS) 
      ASSIGN 40 TO KGO
      GO TO 100 
40    TRANS(1,1)=(-CP*CX*SO-SP*CO)/RAD
      TRANS(2,1)=(CP*CX*CO-SP*SO)/RAD
      TRANS(1,2)=(SP*CX*SO-CP*CO)/RAD
      TRANS(2,2)=(-SP*CX*CO-CP*SO)/RAD
      TRANS(1,3)=-SX*SO/RAD
      TRANS(2,3)=SX*CO/RAD
      DO 45 I=1,3 
45    TRANS(3,I)=0. 
      RETURN
100   CP=COS(PHI/RAD) 
      CX=COS(CHI/RAD) 
      CO=COS(OMEGA/RAD) 
      SP=SIN(PHI/RAD) 
      SX=SIN(CHI/RAD) 
      SO=SIN(OMEGA/RAD) 
      GO TO KGO(10,20,30,40)
      END 
      SUBROUTINE SOLVE(A,B,C) 
      DIMENSION A(3,3),B(3),C(3)
      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 
