                     FUNCTION IPAK(I1,I2,I3,I4)                         
C                                                                       
C     FUNCTION REQUIRES AN INTEGER WORD LENGTH OF AT LEAST 32 BITS.     
C     IT STORES FOUR INTEGER NUMBERS LESS THAN 255 IN  ONE COMPUTER     
C                              WORD.                                    
C                         BY L. W. FINGER                               
C                                                                       
      IPAK=((I1*256+I2)*256+I3)*256+I4                                  
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
                     SUBROUTINE RDEGN(KCYL,NEW)                         
C                                                                       
C     SUBROUTINE TO READ IN THE DEPENDENCY RELATIONSHIPS AND ARRANGE    
C              THEM IN ORDER OF INDEPENDENT PARAMETERS.                 
C                                                                       
C               ORIGIN IS A ROUTINE BY L. W. FINGER                     
C                                                                       
C         ERROR CORRECTION AND ADDITION OF THE CONSTRAINT ORDERS        
C                                                                       
C     (N)EW   =  READ NEW CONSTRAINTS.                                  
C     (O)LD   =  TAKE OLD CONSTRAINTS AS NEW ONES.                      
C     (P)RINT =  PRINT OLD CONSTRAINTS AND TAKE THESE AS NEW ONES.      
C     (A)DD   =  ADD THE CONSTRAIN TO THE OLD ONES.                     
C     (D)ELETE=  DELETE THIS CONSTRAIN AND TAKE REST AS NEW ONES.       
C     (E)ND   =  END OF INPUT LIST.                                     
C                                                                       
C                    BY UDO ZUCKER    DECEMBER 1980                     
C                                                                       
C     INPUT VARIABLE                                                    
C     --------------                                                    
C                                                                       
C     KCYL   NUMBER OF CURRENT CYCLE. (A PRIORI RUN=0)                  
C     NEW    IS NEW(7) INPUT FLAG FOR CONSTRAINTS                       
C                                                                       
C                                                                       
      DOUBLE PRECISION ICARD,NAME                                       
C                                                                       
      DIMENSION SAVE(60),KDEGN(2640),AA(2640),KSTORE(1000),ASTORE(1000) 
      DIMENSION IDPAR(2),INPAR(2,5),CNPAR(5),ITYP(7),ICARD(3),NAME(3)   
C                                                                       
      COMMON/D/ACALC,AIND(60),ANOM(2,20),BCALC,BIND(60),DADB(6,60),DADO(
     & 60),DADX(3,60),DBDB(6,60),DBDO(60),DBDX(3,60),H(3),IH(3),ISCAT(2,
     & 60),ISCL,MFSEL,OCCA(60),OCCUP(60), SCALE(10), SCAT(20), SITE(60),
     & DADC(10,60), DBDC(10,60), DADD(15,60), DBDD(15,60),AMINUS,BMINUS,
     & DRDE(7),NDEGN,VORFAK,IPAREX(7),DAALFA(60),DBALFA(60)             
      COMMON/TAPE/IN,INN,IOUT,ISC1,ISC2,ISC3,ISC4,ISC5,IPUN             
      COMMON/G/BUFF(1023),FCURVE(9,20),IFMT(18),IWT,NSCAT,BETFAC,WVLNTH,
     & WTFAC                                                            
C                                                                       
      EQUIVALENCE (SAVE,AIND),(KDEGN,AA,DADC(1,37)),(BUFF,ASTORE,KSTORE)
C                                                                       
      DATA ITYP/1H ,1HN,1HO,1HA,1HD,1HP,1HE/                            
      DATA ICARD/8HN(EW)   ,8HO(LD)   ,8HP(RINT) /                      
      DATA NAME /8HO(LD)   ,8HA(DD)   ,8HD(ELETE)/                      
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)
 2000 FORMAT(/13H0ERRORMESSAGE/1X,12(1H-)//1X,A8,15H MUST BE FIRST ,    
     & 23HCONSTRAINT-CARD SYMBOL.//16H0JOB TERMINATED.)                 
 2001 FORMAT(/8H0WARNING/1X,7(1H-)//30H A(DD) AFTER N(EW) IN DEPENDEN,  
     &19HCY LIST IS ILLEGAL./35H0CARD ASSUMED AS CONTINUATION-CARD.//   
     &1H )                                                              
 2002 FORMAT(/8H0WARNING/1X,7(1H-)//25H D(ELETE) AFTER N(EW) IN ,       
     &27HDEPENDENCY LIST IS ILLEGAL./27H0CARD ASSUMED AS NOT GIVEN.//   
     & 1H )                                                             
 2003 FORMAT(/13H0ERRORMESSAGE/1X,12(1H-)//24H TOO MANY CONSTRAINTS TO  
     &,7HDELETE./44H0USING N(EW) CONSTRAINTS SEEMS TO BE BETTER./       
     & 16H0JOB TERMINATED.)                                             
 2004 FORMAT(/8H0WARNING/1X,7(1H-)//25H P(RINT) SHOULD BE FIRST ,       
     &24HCARD OF DEPENDENCY LIST./27H0CARD ASSUMED AS NOT GIVEN.//1H )  
 2005 FORMAT(/13H0ERRORMESSAGE/1X,12(1H-)//23H IN AN * A PRIORI * RUN   
     & ,2X,A8,16H IS MEANINGLESS./16H0JOB TERMINATED.)                  
 9000 FORMAT(A1,I2,I3,5(2I3,F8.0))                                      
 9011 FORMAT(/13H0ERRORMESAGE/1X,12(1H-)//21H0TOO MANY DEPENDENCY ,     
     & 10HRELATIONS./16H0JOB TERMINATED.)                               
 9016 FORMAT(/13H0ERRORMESSAGE/1X,12(1H-)//21H0IMPROPER DEPENDENCY ,    
     & 12HINFORMATION./16H0JOB TERMINATED.)                             
C                                                                       
C---- INITIALIZE                                                        
C                                                                       
      NDEL  = 0                                                         
      KPDEL = 1                                                         
      KP    = 1                                                         
      MOLD  = 4                                                         
      NOT1  = 0                                                         
      NOT3  = 0                                                         
      KPMAX = 500                                                       
      NOLD  = NDEGN                                                     
      NDEGN = 0                                                         
      KMAX  = 2640                                                      
      M     = 0                                                         
      KP    = 1                                                         
      IERR  = 0                                                         
C                                                                       
C---- READ CONSTRAINT                                                   
C                                                                       
  10  READ(IN,9000)LABEL,IDPAR,((INPAR(I,J),I=1,2),CNPAR(J),J=1,5)      
C                                                                       
C---- LAST CARD?  (CAN BE MARKED WITH E(ND) TOO.)                       
C                                                                       
      IF(LABEL.NE.ITYP(1)) GOTO 900                                     
      IF(IDPAR(1).EQ.0.AND.IDPAR(2).EQ.0) GOTO 66                       
C                                                                       
C---- M IS NUMBER OF CONSTRAINT-CARDS READ IN.                          
  900 M = M + 1                                                         
C---- INTERPRET LABEL                                                   
C                                                                       
      DO 8 I = 1,7                                                      
      J = I                                                             
    8 IF(LABEL.EQ.ITYP(I)) GOTO 9                                       
C                                                                       
C---- BRANCH CORRESPONDING TO SYMBOL READ IN. FIRST COMPUTED  GOTO      
C     WHEN SYMBOL READ FIRST, SECOND WHEN A BLANK HAS TO BE INTER-      
C     PRET AS CONTINUATION CARD (WITH SAME SYMBOL)                      
C                                                                       
    9 GOTO (12,11,22,33,44,55,66),J                                     
   12 GOTO (111,112,44,11),MOLD                                         
C                                                                       
C---- A N(EW) WAS READ IN. N(EW) MUST BE FIRST CARD.                    
C                                                                       
   11 MOLD = 1                                                          
      IF(M.NE.1) GOTO 1000                                              
C                                                                       
C---- NOT1=1 MEANS N(EW) WAS READ BEFORE. MOLD IS DESIGNATOR FOR SYMBOL.
C                                                                       
      NOT1 = 1                                                          
  111 MOLD = 1                                                          
C                                                                       
  112 DO 20 I=1,5                                                       
      IF(INPAR(1,I).EQ.0.AND.INPAR(2,I).EQ.0.AND.ABS(CNPAR(I)).LT.1E-8)      
     & GO TO 10                                                         
      IF(INPAR(2,I).LT.0)GO TO 1222                                     
      KDEGN(KP)=IPAK(INPAR(1,I),INPAR(2,I),IDPAR(1),IDPAR(2))           
      AA(KP+1)=CNPAR(I)                                                 
      KP=KP+2                                                           
      NDEGN=NDEGN+1                                                     
      IF(KP.LT.KMAX)GO TO 20                                            
C---- TOO MANY RELATIONS                                                
 1111 WRITE(IOUT,9011)                                                  
      CALL EXIT                                                         
C                                                                       
C---- 'RUN' INDICATED ON DEPENDENCY CARD.                               
 1222 IF(MOD(I,2).NE.0.OR.INPAR(1,I).NE.0)GO TO 16                      
      JMAX=-INPAR(2,I)                                                  
      JMIN=INPAR(2,I-1)+1                                               
      DO 15 JX1=JMIN,JMAX                                               
      J=JX1                                                             
      KDEGN(KP)=IPAK(INPAR(1,I),J,IDPAR(1),IDPAR(2))                    
      AA(KP+1)=0.                                                       
      KP=KP+2                                                           
      NDEGN=NDEGN+1                                                     
      IF(KP-KMAX)15,1111,1111                                           
   15 CONTINUE                                                          
      GO TO 20                                                          
C---- ERROR ON DEPENDENCY CARD                                          
   16 WRITE(IOUT,9016)                                                  
      CALL EXIT                                                         
C                                                                       
   20 CONTINUE                                                          
      GO TO 10                                                          
C                                                                       
C---- O(LD) WAS READ IN.  IT MUST BE FIRST CARD.                        
C                                                                       
   22 MOLD  = 2                                                         
      IERR  = 1                                                         
      NEW   = 4                                                         
      IF(KCYL.EQ.0) GOTO 1005                                           
      IF(M.NE.1) GOTO 1000                                              
      NDEGN = NOLD                                                      
      GOTO 9999                                                         
C                                                                       
C---- A(DD) A CONSTRAIN                                                 
C                                                                       
   33 IF(NOT1.NE.0) GOTO 1001                                           
      MOLD  = 2                                                         
      IF(KCYL.EQ.0) GOTO 1005                                           
      IF(NOT3.NE.0) GOTO 112                                            
      NOT3  = 1                                                         
      NDEGN = NOLD                                                      
      KP    = NDEGN*2 + 1                                               
      GOTO 112                                                          
C                                                                       
C---- D(ELETE) CONSTRAINT                                               
C                                                                       
   44 IF(NOT1.NE.0) GOTO 1002                                           
      MOLD  = 3                                                         
      IF(KCYL.EQ.0) GOTO 1005                                           
      IF(NOT3.NE.0) GOTO 441                                            
      NDEGN = NOLD                                                      
      KP    = NDEGN*2 + 1                                               
  441 DO 450  I = 1,5                                                   
      IF(INPAR(1,I).EQ.0.AND.INPAR(2,I).EQ.0.AND.ABS(CNPAR(I)).LT.1E-8)      
     & GO TO 10                                                         
      IF(INPAR(2,I).LT.0)GO TO 451                                      
      KSTORE(KPDEL)=IPAK(INPAR(1,I),INPAR(2,I),IDPAR(1),IDPAR(2))       
      ASTORE(KPDEL+1)=CNPAR(I)                                          
      KPDEL=KPDEL+2                                                     
      NDEL =NDEL +1                                                     
      IF(KPDEL.LT.KPMAX)GO TO 450                                       
C---- TOO MANY RELATIONS                                                
  452 WRITE(IOUT,2003)                                                  
      CALL EXIT                                                         
C                                                                       
C---- 'RUN' INDICATED ON DEPENDENCY CARD.                               
  451 IF(MOD(I,2).NE.0.OR.INPAR(1,I).NE.0) GO TO 456                    
      JMAX=-INPAR(2,I)                                                  
      JMIN=INPAR(2,I-1)+1                                               
      DO 457 JX1=JMIN,JMAX                                              
      J=JX1                                                             
      KSTORE(KPDEL)=IPAK(INPAR(1,I),J,IDPAR(1),IDPAR(2))                
      ASTORE(KPDEL+1)=0.                                                
      KPDEL=KPDEL+2                                                     
      NDEL =NDEL +1                                                     
      IF(KPDEL-KPMAX)457,452,452                                        
  457 CONTINUE                                                          
      GO TO 450                                                         
C---- ERROR ON DEPENDENCY CARD                                          
  456 WRITE(IOUT,9016)                                                  
      CALL EXIT                                                         
C                                                                       
  450 CONTINUE                                                          
      GO TO 10                                                          
C                                                                       
C---- P(RINT) WAS READ IN.                                              
C                                                                       
   55 IF(M.NE.1) GOTO 1004                                              
      NDEGN=NOLD                                                        
      GOTO 9999                                                         
C                                                                       
C---- E(ND) WAS READ.                                                   
C                                                                       
   66 IERR=4                                                            
C                                                                       
C---- DELETE CONSTRAINTS IF DESIRED.                                    
      IF(NDEL.EQ.0) GOTO 100                                            
      KPDEL = KPDEL - 2                                                 
      KP    = KP    - 2                                                 
C                                                                       
      DO 77 I = 1,KPDEL,2                                               
         DO 78 J = 1,KP,2                                               
         IF(KSTORE(I)-KDEGN(J))78,79,78                                 
   79    IF(ABS(ASTORE(I+1)-AA(J+1)).GT.1.E-8) GOTO 78                  
         KDEGN(J) = 0                                                   
   78    CONTINUE                                                       
   77 CONTINUE                                                          
C                                                                       
      IP=0                                                              
      DO 80 I = 1,KP,2                                                  
      IF(KDEGN(I)) 82,81,82                                             
   81 IP = IP + 2                                                       
      GOTO 80                                                           
   82 KDEGN(I-IP) = KDEGN(I)                                            
      AA(I-IP+1)  = AA(I+1)                                             
   80 CONTINUE                                                          
C                                                                       
      NDEGN = NDEGN - IP/2                                              
C                                                                       
C---- RESHUFFLE CONSTRAINT INDICATORS TO GROUP THEM ACCORDING TO THE    
C     INDEPENDENT PARAMETER.                                            
C                                                                       
 100  II=NDEGN-1                                                        
      IF(II.EQ.0) GOTO 200                                              
      DO 110 I=1,II                                                     
      JJ=I+1                                                            
      DO 110 J=JJ,NDEGN                                                 
      IF(KDEGN(2*I-1)-KDEGN(2*J-1))110,110,105                          
  105 ITEMP=KDEGN(2*I-1)                                                
      TEMP=AA(2*I)                                                      
      KDEGN(2*I-1)=KDEGN(2*J-1)                                         
      AA(2*I)=AA(2*J)                                                   
      KDEGN(2*J-1)=ITEMP                                                
      AA(2*J)=TEMP                                                      
  110 CONTINUE                                                          
  200 CONTINUE                                                          
 9999 RETURN                                                            
C                                                                       
C-----------------------------------------------------------------------
C     ERROR PART                                                        
C                                                                       
C---- O(LD) OR N(EW) OR P(RINT) MUST BE FIRST CARD OF CONSTRAINTS.      
C                                                                       
 1000 WRITE(IOUT,2000)ICARD(MOLD)                                       
      CALL EXIT                                                         
C                                                                       
C---- A(DD) AFTER N(EW) IS ILLEGAL, USES AS N(EW).                      
C                                                                       
 1001 WRITE(IOUT,2001)                                                  
      GOTO 111                                                          
C                                                                       
C---- D(ELETE) AFTER N(EW) IS ILLEGAL, USES AS NOT GIVEN.               
C                                                                       
 1002 WRITE(IOUT,2002)                                                  
      GOTO 10                                                           
C                                                                       
C---- P(RINT) MUST BE FIRST CARD.                                       
C                                                                       
 1004 WRITE(IOUT,2004)                                                  
      GOTO 10                                                           
C                                                                       
C---- IN AN *A PRIORI RUN*  D(ELETE), A(DD) AND O(LD) IS MEANING-       
C     LESS.                                                             
 1005 IF(IERR.EQ.0) IERR=MOLD                                           
      WRITE(IOUT,2005)NAME(IERR)                                        
      CALL EXIT                                                         
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
         SUBROUTINE SYMINV (A,NRANK,SCR,NCB,B,EPS,DET,JRANK,IASIZ)            
C                                                                       
C                                                                       
C     THE METHOD USED BY THIS MATRIX INVERSION ROUTINE HAS BEEN         
C                  DESCRIBED BY BUSING AND LEVY                         
C     'A PRCEDURE FOR INVERTING LARGE MATRICES' COMM. ACM, 5,45         
C              AND WAS CODED BY MR.JAMES MUNDSTOCK                      
C     OF  THE  NUMERICAL  ANALYSIS  CENTER  OF  THE  UNIVERSITY         
C                         OF MINNESOTA.                                 
C                                                                       
C     REVISED AND/OR RESEQUENCED OCT. 28, 1968  BY L. W. FINGER         
C                                                                       
C     ARGUMENT LIST                                                     
C     -------------                                                     
C                                                                       
C     A     = INPUT AND OUTPUT ARRAY - STORED BY ROWS WITH ONLY         
C             THE  UNIQUE PORTION IN STORAGE                            
C     NRANK = INPUT SIZE OF SYSTEM OF EQUATIONS                         
C     SCR   = A SCRATCH ARRAY OF DIMENSION NRANK                        
C     B     = MATRIX  TO  BE PREMULTIPLIED  BY THE INVERSE OF A         
C     NCB   = NUMBER OF COLUMNS OF B                                    
C     EPS   = SIZE  OF SMALLEST  ELEMENT  TO BE USED AS A PIVOT         
C     DET   = OUTPUT VALUE OF THE DETERMINANT OF A                      
C     JRANK = ACTUAL  RANK  OF A    ( THIS IS AN OUTPUT VALUE.)         
C                                                                       
      DIMENSION A(IASIZ),B(NRANK,NCB),SCR(NRANK)                                  
C                                                                       
C---- INITIALIZE                                                        
C                                                                       
      DET=1.0                                                           
      DO 100 I=1,NRANK                                                  
  100 SCR(I)=1.0                                                        
C                                                                       
C---- START MASTER LOOP                                                 
C                                                                       
      DO 280 I=1,NRANK                                                  
C---- SEARCH FOR PIVOT TESTING ONLY THE DIAGONAL ELEMENTS               
      BIG=0.0                                                           
C---- PRESET INDEX AND STEP FOR DIAGONAL ELEMENTS                       
      J =-NRANK                                                         
      KJ=NRANK+1                                                        
      DO 130 K=1,NRANK                                                  
C---- AUGMENT INDEX AND REDUCE STEP                                     
      J=J+KJ                                                            
      KJ=KJ-1                                                           
C---- CHECK IF THIS ELEMENT ALREADY USED AS PIVOT                       
      IF(SCR(K))130,110,110                                             
C---- J INDEXES A(K,K)                                                  
  110 TEMP=ABS(A(J))                                                    
      IF(TEMP-BIG)130,130,120                                           
  120 BIG=TEMP                                                          
      KK=K                                                              
      KL=J                                                              
  130 CONTINUE                                                          
C---- COMPARE LARGEST DIAGONAL ELEMENT WITH EPS                         
      IF(BIG-EPS)140,150,150                                            
C---- SINGULARITY RETURN                                                
  140 JRANK=I-1                                                         
      RETURN                                                            
  150 SCR(KK)=-1.0                                                      
C---- KL INDEXES A(KK,KK)                                               
      DET=DET*A(KL)                                                     
      Q=1.0/A(KL)                                                       
C---- SET UP INDEXING FOR ELIMINATION STEPS                             
      J=KK-NRANK                                                        
      NJ=NRANK+1                                                        
      K=0                                                               
      IK=KK-1                                                           
      JK=KK+1                                                           
      IF(IK)210,210,160                                                 
  160 DO 200 N=1,IK                                                     
      NJ=NJ-1                                                           
      J=J+NJ                                                            
C---- J INDEXES A(KK,N)                                                 
      TEMP=Q*A(J)                                                       
      NI=J                                                              
      MJ=NJ                                                             
      DO 170 M=N,IK                                                     
      K=K+1                                                             
C---- J INDEXES A(KK,M)                                                 
C---- K INDEXES A(M,N)                                                  
      A(K)=A(K)-TEMP*SCR(M)*A(J)                                        
      MJ=MJ-1                                                           
  170 J=J+MJ                                                            
      J=NI                                                              
C---- J INDEXES A(KK,N)                                                 
      A(J)=TEMP                                                         
      K=K+1                                                             
      L=KL                                                              
      IF(JK-NRANK)180,180,200                                           
  180 DO 190 M=JK,NRANK                                                 
      L=L+1                                                             
      K=K+1                                                             
C---- K INDEXES A(M,N)                                                  
C---- L INDEXES A(M,KK)                                                 
  190 A(K)=A(K)-TEMP*A(L)                                               
  200 CONTINUE                                                          
  210 K=KL+NRANK-KK                                                     
      J=KL                                                              
      IF(JK-NRANK)220,220,250                                           
  220 DO 240 N=JK,NRANK                                                 
      J=J+1                                                             
C---- J INDEXES A(N,KK)                                                 
      TEMP=SCR(N)*A(J)*Q                                                
      NI=J                                                              
      DO 230 M=N,NRANK                                                  
      K=K+1                                                             
C---- K INDEXES A(M,N)                                                  
C---- J INDEXES A(M,KK)                                                 
      A(K)=A(K)-TEMP*A(J)                                               
  230 J=J+1                                                             
  240 J=NI                                                              
  250 K=KL                                                              
C---- K INDEXES A(KK,KK)                                                
      A(K)=Q                                                            
      IF(JK-NRANK)260,260,280                                           
  260 DO 270 M=JK,NRANK                                                 
      K=K+1                                                             
C---- K INDEXES A(M,KK)                                                 
  270 A(K)=-Q*A(K)                                                      
  280 CONTINUE                                                          
C---- INVERSION COMPLETE, EXIT IF INVERSION ONLY                        
      IF(NCB.EQ.1.AND.NRANK.EQ.3) GOTO 360
      GOTO 290
C---- PREMULTIPLY THE NCB COLUMNS OF B BY A INVERSE                     
  290 DO 350 K=1,NCB                                                    
      DO 330 I=1,NRANK                                                  
      Q=0.0                                                             
      L=I-NRANK                                                         
      KK=NRANK+1                                                        
      DO 300 J=1,I                                                      
      KK=KK-1                                                           
      L=L+KK                                                            
C---- L INDEXES A(I,J)                                                  
  300 Q=Q+A(L)*B(J,K)                                                   
      IF(I-NRANK)310,330,330                                            
  310 KJ=I+1                                                            
      DO 320 J=KJ,NRANK                                                 
      L=L+1                                                             
C---- L INDEXES A(J,I)                                                  
  320 Q=Q+A(L)*B(J,K)                                                   
  330 SCR(I)=Q                                                          
      DO 340 J=1,NRANK                                                  
  340 B(J,K)=SCR(J)                                                     
  350 CONTINUE                                                          
C---- SET NON-SINGULAR FLAG AND EXIT                                    
  360 JRANK=NRANK                                                       
      RETURN                                                            
      END                                                               
C                                                                       
                         SUBROUTINE TRNFRM                              
C                                                                       
C              ROUTINE FOR SPECIAL PARAMETER TREATMENT                  
C                                                                       
C                                                                       
C     SUBROUTINE  CALLED  ONCE PER  JOB,  AFTER CELL CONSTANTS HAVE BEEN
C     READ IN,  TO SET UP  TRANSFORMS  TO TAKE VECTORS FROM CRYSTAL AXES
C     TO AN  ORTHONORMAL  COORDINATE  SYSTEM, AND TO TAKE THE SECOND AND
C              THIRD CUMULANTS BACK TO THE CRYSTAL SYSTEM.              
C                                                                       
C                         BY  L. W. FINGER                              
C                                                                       
      COMMON/A/BETA(6,60),G(3,3),ISOT(60),NATOM,TAG(2,60),TITLE(18),XYZ(
     &3,60),PI,CR(10,60),DR(15,60),ALF(60),ITTYP(60)                             
      COMMON/XPAR/NPAR,ATOX(3,3),XTOA(3,3),PAR(60),PARNAM(2,60),IPAR(60)
     &,IRPAR(60),OXYZ(3),SXYZ(3),OXP(3),T(3,3),EL(3,3),S(3,3),TRNS(3,3),
     & TTRNS(3,3), TRANSI(3,3), TTRNSI(3,3), RTEM(3), BTEM(6),  GTEM(10)
C                                                                       
      ATOX(1,1) = SQRT(G(1,1))                                          
      ATOX(1,2) = G(1,2)/ATOX(1,1)                                      
      ATOX(1,3) = G(1,3)/ATOX(1,1)                                      
      ATOX(2,1) = 0.                                                    
      ATOX(2,2) = SQRT(G(2,2)-ATOX(1,2)**2)                             
      ATOX(2,3) = (G(2,3)-ATOX(1,2)*ATOX(1,3))/ATOX(2,2)                
      ATOX(3,1) = 0.                                                    
      ATOX(3,2) = 0.                                                    
      ATOX(3,3) = SQRT(G(3,3)-ATOX(1,3)**2-ATOX(2,3)**2)                
C                                                                       
      DO 10 I  =1,3                                                     
   10 XTOA(I,I)=1./ATOX(I,I)                                            
      XTOA(1,2)=-ATOX(1,2)/(ATOX(1,1)*ATOX(2,2))                        
      XTOA(2,3)=-ATOX(2,3)/(ATOX(2,2)*ATOX(3,3))                        
      XTOA(1,3)=(ATOX(1,2)*ATOX(2,3)/ATOX(2,2)-ATOX(1,3))/(ATOX(1,1)*ATO
     &X(3,3))                                                           
      XTOA(2,1)=0.                                                      
      XTOA(3,1)=0.                                                      
      XTOA(3,2)=0.                                                      
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
                           SUBROUTINE PAREXT                            
C                                                                       
C          SUBROUTINE TO READ IN SPECIAL PARAMETER INFORMATION.         
C                                                                       
C                           BY L. W. FINGER                             
C                                                                       
C     READ 1 CARD FOR EACH OF N SPECIAL PARAMETERS,      CONTAINING THE 
C     FOLLOWING INFORMATION.                                            
C                                                                       
C     COL 1-8     NAME OF PARAMETER  - UP TO 8 ALPHANUMERIC CHARACTERS. 
C     COL 9-10    TYPE OF PARAMETER. 1=SPECIAL COORDINATE SYSTEM (SHAPE 
C                 CONSTRAINT) OR RIGID BODY THERMAL MOTION,   2=DEFINED 
C                 BY  USER  IN   SUBROUTINE SPDERI (SUPPLIED AS DUMMY). 
C                 0 INDICATES A LINEAR FUNCTION, WITH COEFFICIENTS   AS 
C                 ON CONSTRAINT  CARDS.                                 
C     COL 11-12   NUMBER INDICATING WHICH OF SEVERAL SPECIAL COORDINATE 
C                 OR RIGID BODY SYSTEMS THIS PARAMETER BELONGS TO.      
C     COL 13-14   PARAMETER DESIGNATOR.  IF TYPE 1,  THE  CODE  IS   AS 
C                 FOLLOWS 1-X COORDINATE OF ORIGIN (IN CRYSTAL SYSTEM), 
C                 2-Y COORDINATE OF  ORIGIN, 3-Z  COORDINATE OF ORIGIN, 
C                 4-OMEGA ANGLE FOR TRANSFORM TO STANDARD, ORTHONORMAL  
C                 SYSTEM, 5-CHI  ANGLE  FOR TRANSFORM, 6-PHI ANGLE FOR  
C                 TRANSFORM,  7-X COORDINATE  IS SPECIAL SYSTEM,   8-Y  
C                 COORDINATE IS SPECIAL SYSTEM, 9-Z COORDINATE IN SPE-  
C                 CIAL SYSTEM,  10-T11, 11-T22, 12-T33,13-T12, 14-T13,  
C                 15-T23, 16-L11, 17-L22, 18-L33, 19-L12, 20-L13,  21-  
C                 L23, 22-S11, 23-S22, 24-S33, 25-S12, 26-S13, 27-S23,  
C                 28-S21,29-S31, 30-S32.                                
C     COL 15-24   INITIAL VALUE OF THIS PARAMETER.                      
C                                                                       
C                                                                       
      COMMON/XPAR/NPAR,ATOX(3,3),XTOA(3,3),PAR(60),PARNAM(2,60),IPAR(60)
     & ,IRPAR(60),OXYZ(3),SXYZ(3),OXP(3),T(3,3),EL(3,3),S(3,3),TRNS(3,3)
     & ,TTRNS(3,3), TRANSI(3,3), TTRNSI(3,3), RTEM(3), BTEM(6), GTEM(10)
      COMMON/TAPE/IN,INN,IOUT,ISC1,ISC2,ISC3,ISC4,ISC5,IPUN             
C                                                                       
      DIMENSION ITEM(3)                                                 
      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                                                                       
      IF(NPAR.EQ.0)RETURN                                               
      WRITE(IOUT,9001)                                                  
      DO 100 I=1,NPAR                                                   
      READ(IN,9000)(PARNAM(J,I),J=1,2),ITEM,PAR(I)                      
      IPAR(I)=IPAK(ITEM(1),ITEM(2),0,ITEM(3))                           
      WRITE(IOUT,9002)(PARNAM(J,I),J=1,2),ITEM,PAR(I)                   
 100  CONTINUE                                                          
      RETURN                                                            
C                                                                       
 9000 FORMAT(2A4,3I2,F10.4)                                             
 9001 FORMAT(19H0SPECIAL PARAMETERS)                                    
 9002 FORMAT(1X,2A4,3(2X,I2),2X,F10.4)                                  
C                                                                       
      ENTRY PUNPAR                                                      
C                                                                       
      DO 200 I=1,NPAR                                                   
      ITEM(1)=I1(IPAR(I))                                               
      ITEM(2)=I2(IPAR(I))                                               
      ITEM(3)=I4(IPAR(I))                                               
      WRITE(IPUN,9000)(PARNAM(J,I),J=1,2),ITEM,PAR(I)                   
 200  CONTINUE                                                          
C                                                                       
      RETURN                                                            
      END                                                               
