              SUBROUTINE SFAC(MODELL,IZAEHL,IRAUS)                              
C                                                                               
C                  STRUCTURE FACTOR GENERATION                                  
C                                                                               
C   MODELL=1 STRUCTURE FACTOR EQUATION BASED ON THE FOURIER TRANSFORM           
C   OF THE EDGEWORTH  EXPANSION. EQUATION HAD BEEN SUGGESTED BY C. K.           
C     JOHNSON (1969). THE EQUATION WAS POGRAMMED BY LARRY W. FINGER.            
C   MODELL=2 STRUCTURE FACTOR EQUATION BASED ON THE FOURIER TRANSFORM           
C   OF THE GRAM-CHARLIER SERIES EXPANSION. SERIES HAS BEEN  PREFERRED           
C     BY UDO ZUCKER, MAINLY BECAUSE OF THE EXACT FOURIER TRANSFORM.             
C   MODELL=3 LIKE 1, BUT INSTEAD OF A PRODUCT OF MILLER INDICES SPEC-           
C   IAL POLYNOMS ARE USED. THESE FUNCTIONS  CORRESPONDS  TO THE THREE           
C   DIMENSIONAL HERMIT POLYNOMS. BUT ONLY  CORRESPONDENCE  IS  GIVEN.           
C   THAT IS WHY THE LABEL HERMIT POLYNOMS  USED BY  C.K. JOHNSON  AND           
C                      OTHERS IS NOT SUITABLE.                                  
C   THE MODELS 2-4 ARE CALCULATED AND PROGRAMMED BY UDO ZUCKER  (MPI)           
C   REMIND THAT THE MODELS NO. 3 AND 4 ARE USING ONE PARAMETER  MORE.           
C                                                                               
C   EXPANSION OF THERMAL VIBRATIONS TILL SIXTH ORDER ADDED. (ONLY FOR           
C   THE GRAM CHARLIER SERIES EXPANSION. ONLY FOR THIS  EXPANSION  THE           
C   CORRECT AND COMPLETE FOURIER TRANSFORM IS GIVEN AND PROGRAMMABLE.           
C                                                                               
C        PROGRAMME WRITTEN  BY UDO ZUCKER  OCT. 78 - DEC. 80                    
C                                                                               
C               ORIGIN WAS A PROGRAMME OF L. W. FINGER                          
C                                                                               
C                                                                               
C     INPUT VARIABLES                                                           
C     ---------------                                                           
C                                                                               
C     MODELL   = 1,2,3,4 ...  STRUCTURE FACTOR MODEL, SEE COMMENT ABOVE.        
C     IZAEHL   = N ...... N ATOMS WITH A ALPHA PARAMETER NOT EQUAL ZERO.        
C     IRAUS    = 0/1 REFLECTION DO /DO NOT NEED COMPUT.  OF DERIVATIVES.        
C                                                                               
      DIMENSION SAVE(60),KDEGN(2640),AA(2640), QTRANS(6,24),V(3), VV(6),        
     & Q(3,3),D3HERM(10),D4HERM(15),ER(21,6),FR(28,6), DAE(21), DBE(21),        
     & DAF(28),DBF(28),HHHHH(21),HHHHJ(28),DADE(21,6), DBDE(21,6), DADF(        
     & 28,6),DBDF(28,6),T1(15),T2(15)                                           
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/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,ISAVE, 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/LSMAT/B(200),IEXT,IPARA(1000),MRANK,KCYL,LABEL(3),JSCRF(10         
     & ),NRSF, SIGMA, MDEGN, MAXIMA, IALFA(60),  RDAMP(200), IDAMP(200),        
     & A(2200)                                                                  
      COMMON/C/     CELLP(3,4),ICENT,JCELL,NSYM,RMAT(3,3,24),TRANS(3,24)        
      COMMON/E/DMAX,DMIN,RECELL(6),EXANG(6),EXTRAN(12)                          
      COMMON/F/DELTA,FCALC,FOBS,FOMAX,FOMIN,IREJ,MREF,MREJ,NREJ,NUMBER(3        
     & 6),RDEN(36), RHO, RHOMN,RHOMX, RNUM(36),WT,EXBETA,EXX,EXY,GAMMAP,        
     & FISGN,IFSGN,EXTIN(6),KEXTRF(7),DELMAX,IANIGC,GAMMAS,IZACH,RNEUEX(        
     & 7),NEUEXT,NDISTR,NIEXT                                                   
      COMMON/LS/ DAB(6),DAX(3),DBB(6),DBX(3),HCP(4),HHJ(6),HP(3,24),HPR(        
     & 3),TPR(24), HHHJ(10), DAC(10), DBC(10), HHHH(15), DAD(15),DBD(15)        
      COMMON/TAPE/IN,INN,IOUT,ISC1,ISC2,ISC3,ISC4,ISC5,IPUN                     
C                                                                              
      EQUIVALENCE (SAVE,AIND), (KDEGN,AA,DADC(1,37)), (V,DAX), (VV,HHJ),        
     & (D3HERM,DAE),(D4HERM,DAE(11)),(D4HERM(12),DBE),                          
     & (QTRANS,DBE(5)), (QTRANS(5,3),DAF), (QTRANS(3,8),DBF), (QTRANS(4,        
     & 8),HHHHH),(QTRANS(5,12),HHHHJ)                                           
      EQUIVALENCE         (DR(1,41),ER),(DR(7,49),FR),(DADD(1,41),DADE),        
     & (DADD(7,49),DADF),(DBDD(1,41),DBDE),(DBDD(7,49),DBDF)                    
C                                                                               
C                                                                               
C---- FAK1=2/3*PI**4           FAK2=4/3*PI**3             FAK3=1/6*PI**2        
C     FAK5=(2*PI)**5/120       FAK6=(2*PI)**6/720                               
C                                                                               
      DATA TWOPI,FAK1,FAK2,FAK3,FAK5,FAK6/6.2831853,64.939394,41.341702,        
     & 0.01688686,81.605249,85.456817/                                          
C                                                                               
C---- BECAUSE OF DIFFERENT MEANING OF VALUE IN REFINE AND SFAC.                 
C                                                                               
      MFSEL = ISAVE                                                             
C                                                                               
C---- DO WE NEED DERIVATIVES FOR THIS STRUCT. FACTOR. IF REJECTED=NO.           
C                                                                               
      IF(MFSEL.NE.0 .AND. IRAUS.EQ.1 .AND. MREJ.EQ.1)  MFSEL=0                  
C                                                                               
C---- GENERATE TRANSFORMED INDICES                                              
C                                                                               
      DO 10 I=1,NSYM                                                            
      TPR(I)=0.0                                                                
                  DO 10 J=1,3                                                   
                  HP(J,I)=0.0                                                   
                  TPR(I)=TPR(I)+H(J)*TRANS(J,I)*FISGN                           
                          DO 10 K=1,3                                           
   10                     HP(J,I)=HP(J,I)+H(K)*RMAT(K,J,I)*FISGN                
C                                                                               
C---- INITIALIZE                                                                
C                                                                               
      ACALC=0.0                                                                 
      BCALC=0.0                                                                 
      AMINUS=0.0                                                                
      BMINUS=0.0                                                                
C                                                                               
C---- START LOOP THROUGH ATOMS                                                  
C                                                                               
                             DO 450 L=1,NATOM                                   
C                                                                               
C                                                                               
C---- IF MODELL GRETER THREE THE INVERSE MATRIX OF BETA TRANSFORM IT-           
C     SELF LIKE A PRODUCT OF TWO MILLER INDICES.                                
C     IS IS WELL KNOWN FOR  MYSELF (U.Z.) THAT IN THE CURRENT VERSION           
C     FOR EACH REFLECTION EVER AND EVER THIS PART IS CALCULATED.                
C     THE SECOND POSSIBILITY OF SAVING THE RESULTS OF QTRANS FOR EACH           
C     ATOM NEEDS  INSTEAD OF  QTRANS(6,24)  AN ARRAY QTRANS(60,6,24).           
C     I DO NO HAVE 60*6*24 = 8640   ELEMENTS INSTEAD OF    6*24 = 144           
C     WITHOUT INCREASING CORE IN A MEANINGLESS WAY.                             
C     THIS VERSION NEEDS TIME INSTEAD OF CORE. USE GRAM-CHARLIER  IN-           
C     STEAD OF ALPHA-FORMALISM.                                                 
C                                                                               
      IF(MODELL.LT.3 ) GOTO 13                                                  
C---- LAUF=L TO PREVENT COMPILER WARNINGS                                       
      LAUF=L                                                                    
C---- Q IS INVERSE MATRIX OF BETA                                               
      CALL UMKEHR(BETA,LAUF,Q)                                                  
C---- HAS AN ATOM A STRUCTUREFACTOR WITH AN ALPHA .NE. ZERO?                    
      IF(IZAEHL.LE.0) GOTO 13                                                   
C                                                                               
C---- YES TRANSFORM                                                             
      DO 12 I=1,NSYM                                                            
         DO 12 II=1,6                                                           
         IBETA=II                                                               
         IGAMA=II                                                               
         IF(II.EQ.4 .OR. II.EQ.5) IBETA=1                                       
         IF(II.EQ.4) IGAMA=2                                                    
         IF(II.EQ.5) IGAMA=3                                                    
         IF(II.EQ.6) IBETA=2                                                    
         IF(II.EQ.6) IGAMA=3                                                    
         QTRANS(II,I)=0.0                                                       
               DO 12 J=1,3                                                      
               DO 12 K=1,3                                                      
  12           QTRANS(II,I)=QTRANS(II,I)+Q(J,K)                                 
     &           *RMAT(J,IBETA,I)*RMAT(K,IGAMA,I)*FISGN*FISGN                   
C                                                                               
C---- INITIALIZE PARAMETERS OF EACH ATOM                                        
C                                                                               
  13  TF=JCELL*(2-ICENT)                                                        
      EXPIJ=1.0                                                                 
      AAIND=0.0                                                                 
      BBIND=0.0                                                                 
      IF(MFSEL)20,60,20                                                         
   20                        DO 30 M=1,3                                        
C---- DERIVATIVES WTH TO COORDINATES.                                           
      DAX(M)=0.0                                                                
   30 DBX(M)=0.0                                                                
      IF(ISOT(L))40,70,40                                                       
   40                        DO 50 N=1,6                                        
C---- WTH TO ANISOTROPIC TEMPERATURE FACTOR.                                    
      DAB(N)=0.0                                                                
   50 DBB(N)=0.0                                                                
      IF(ISOT(L)-1)80,80,52                                                     
   52 DO 54 N=1,10                                                              
C---- WTH TO THIRD ORDER CUMULANTS.                                             
      DAC(N)=0.0                                                                
   54 DBC(N)=0.0                                                                
C---- WTH TO ALPHA                                                              
      DAALF=0.0                                                                 
      DBALF=0.0                                                                 
      IF(ISOT(L)-2)80,80,55                                                     
   55 DO 56 N=1,15                                                              
C---- WTH TO FOURTH ORDER CUM.                                                  
      DAD(N)=0.0                                                                
   56 DBD(N)=0.0                                                                
      IF(ISOT(L)-4)80,58,59                                                     
   59 DO 591 N=1,28                                                             
C---- WTH TO SIXTH ORDER CUM.                                                   
      DAF(N)=0.0                                                                
  591 DBF(N)=0.0                                                                
   58 DO 581 N=1,21                                                             
C---- WTH TO FIFTH ORDER CUM.                                                   
      DAE(N)=0.0                                                                
  581 DBE(N)=0.0                                                                
      GO TO 80                                                                  
C                                                                               
   60 IF(ISOT(L))80,70,80                                                       
C                                                                               
C---- CALCULATE ISOTROPIC TEMPERATURE FACTOR.                                   
C                                                                               
   70 TF=EXP(-BETA(1,L)*RHO)*TF                                                 
C                                                                               
C                                                                               
C---- START LOOP THROUGH SYMMETRY OPERATIONS.                                   
C                                                                               
C                                                                               
   80                        DO 280 IA=1,NSYM                                   
      ARG=TPR(IA)                                                               
                             DO 90 IB=1,3                                       
      HPR(IB)=HP(IB,IA)                                                         
C---- ARG = 2*PI (HX+KY+LZ)                                                     
   90 ARG=ARG+HPR(IB)*XYZ(IB,L)                                                 
C                                                                               
C---- INITIALIZE   SUMMATION ( H(K)*H(L)*H(M) * CR(KLM)) = CUM3                 
C                  AND CUM5 = SUM (HHHHH * ER)                                  
C                  AND D3ALF= SUM (GERMIT* CR)                                  
C     GERMIT = POLYNOMS OF ALPHA FORMALISM SIMILAR TO HERMIT POLY.              
C                                                                               
      CUM3=0.0                                                                  
      CUM5=0.0                                                                  
      D3ALF=0.0                                                                 
C                                                                               
      IF(ISOT(L))100,130,100                                                    
C                                                                               
C---- GENERATE TRANSFORMED INDEX PRODUCTS FOR ANISO. ATOMS                      
C                                                                               
  100 DO 110 IC=1,3                                                             
  110 HHJ(IC)=HPR(IC)**2                                                        
      HHJ(4)=2.0*HPR(1)*HPR(2)                                                  
      HHJ(5)=2.0*HPR(1)*HPR(3)                                                  
      HHJ(6)=2.0*HPR(2)*HPR(3)                                                  
C                                                                               
C---- CALCULATE ANISOTROPIC TEMPERATURE FACTOR                                  
C                                                                               
      TEMP=0.0                                                                  
      DO 120 ID=1,6                                                             
  120 TEMP=TEMP-HHJ(ID)*BETA(ID,L)                                              
C                                                                               
C---- MEANING OF CUM4,CUM6,... SIMILAR TO CUM3,CUM5,... (SEE ABOVE)             
C                                                                               
      CUM4=0.0                                                                  
      CUM6=0.0                                                                  
      D4ALF=0.0                                                                 
C                                                                               
C---- ALL MODELS ARE THE SAME IN HARMONIC CASE.  THEREFORE THE HARMONIC         
C     PART  ALWAYS CALCULATED UNDER LABEL 129. ( THIS  IS THE  QUICKEST         
C     METHOD.)                                                                  
C                                                                               
      IF(ISOT(L)-1)129,129,122                                                  
C                                                                               
C---- IF MODELL.GT.2 THE CUMULANT PRODUCT FOLLOWING BELOW ARE NOT USED.         
C     INSTEAD OF THESE PRODUCTS, WE NEED THE TERMS       G(JKL)*C(JKL).         
C     G(JKL) ARE CALCULATED IN THE ROUTINE GERMIT. G(JKL) ARE THE POLY-         
C     NOMS SIMILAR TO THE HERMIT PLOYNOMS. FROM NOW ON CALLED *GERMIT*-         
C     POLYNOMS.                                                                 
C                                                                               
  122 IF(MODELL.GT.2) GOTO 128                                                  
C                                                                               
C---- CUMULANT PRODUCTS IN  FORM OF H(J)*H(K)*H(L)*C(JKL)*MULTIPLICITY.         
C     MULTIPLIED BY 1.E-3                                                       
      DO 123 ID=1,3                                                             
  123 HHHJ(ID) =   HPR(ID)*HHJ(ID)*1.E-3                                        
      HHHJ( 4) = 3.0*HHJ(1)*HPR(2)*1.E-3                                        
      HHHJ(5 ) = 3.0*HHJ(2)*HPR(1)*1.E-3                                        
      HHHJ(6 ) = 3.0*HHJ(1)*HPR(3)*1.E-3                                        
      HHHJ(7 ) = 3.0*HHJ(3)*HPR(1)*1.E-3                                        
      HHHJ(8 ) = 3.0*HHJ(2)*HPR(3)*1.E-3                                        
      HHHJ(9 ) = 3.0*HHJ(3)*HPR(2)*1.E-3                                        
      HHHJ(10) = 3.0*HHJ(4)*HPR(3)*1.E-3                                        
C                                                                               
      DO 124 ID=1,10                                                            
  124 CUM3 = CUM3 + HHHJ(ID)*CR(ID,L)                                           
C                                                                               
      IF(ISOT(L)-2)128,128,125                                                  
C                                                                               
C---- CALCULATE FOURTH CUMULANT INDEX PRODUCTS                                  
C     MULTIPLIED BY 1.E-4                                                       
C                                                                               
  125 DO 126 ID=1,3                                                             
  126 HHHH(ID) = HHJ(ID)**2 *1.E-4                                              
      HHHH( 4) = 4.0 *HHHJ(1) * HPR(2) * 0.1                                    
      HHHH( 5) = 4.0 *HHHJ(1) * HPR(3) * 0.1                                    
      HHHH( 6) = 6.0 *HHJ(1)  * HHJ(2) * 1.E-4                                  
      HHHH( 7) = 12. *HHJ(1)  * HPR(2) * HPR(3) * 1.E-4                         
      HHHH( 8) = 6.0 *HHJ(1)  * HHJ(3) * 1.E-4                                  
      HHHH( 9) = 4.0 *HHHJ(2) * HPR(1) * 0.1                                    
      HHHH(10) = 12. *HHJ(2)  * HPR(1) * HPR(3) * 1.E-4                         
      HHHH(11) = 12. *HHJ(3)  * HPR(1) * HPR(2) * 1.E-4                         
      HHHH(12) = 4.0 *HHHJ(3) * HPR(1) * 0.1                                    
      HHHH(13) = 4.0 *HHHJ(2) * HPR(3) * 0.1                                    
      HHHH(14) = 6.0 *HHJ(2)  * HHJ(3) * 1.E-4                                  
      HHHH(15) = 4.0 *HHHJ(3) * HPR(2) * 0.1                                    
C                                                                               
      DO 127 ID=1,15                                                            
  127 CUM4 = CUM4 + HHHH(ID)*DR(ID,L)                                           
C                                                                               
      IF(ISOT(L)-3) 128,128,1251                                                
C                                                                               
C---- CALCULATE FIFTH CUMULANT INDEX PRODUCTS                                   
C     ELEMENTS OF TENSOR HAD BEEN MULTIPLIED BY 1.E5                            
C                                                                               
 1251 DO 1252 ID=1,3                                                            
 1252 HHHHH(ID) = HHHH(ID)* HPR(ID)* 0.1                                        
      HHHHH(4 ) = HHHH(1) * HPR(2) * 0.5                                        
      HHHHH(5 ) = HHHH(1) * HPR(3) * 0.5                                        
      HHHHH(6 ) = HHHJ(1) * HHJ(2) * 0.1                                        
      HHHHH(7 ) = HHHJ(1) * HHJ(6) * 0.1                                        
      HHHHH(8 ) = HHHJ(1) * HHJ(3) * 0.1                                        
      HHHHH(9 ) = HHHJ(2) * HHJ(1) * 0.1                                        
      HHHHH(10) = HHHJ(8) * HHJ(1) * 0.1                                        
      HHHHH(11) = HHHJ(9) * HHJ(1) * 0.1                                        
      HHHHH(12) = HHHJ(3) * HHJ(1) * 0.1                                        
      HHHHH(13) = HHHH(2) * HPR(1) * 0.5                                        
      HHHHH(14) = HHHH(9) * HPR(3) * 0.5                                        
      HHHHH(15) = HHHJ(10)* HHJ(6) * 2.5 E-2                                    
      HHHHH(16) = HHHH(15)* HPR(1) * 0.5                                        
      HHHHH(17) = HHHH(3) * HPR(1) * 0.5                                        
      HHHHH(18) = HHHH(2) * HPR(3) * 0.5                                        
      HHHHH(19) = HHHJ(2) * HHJ(3) * 0.1                                        
      HHHHH(20) = HHHJ(3) * HHJ(2) * 0.1                                        
      HHHHH(21) = HHHH(3) * HPR(2) * 0.5                                        
C                                                                               
      DO 1253 ID=1,21                                                           
 1253 CUM5 = CUM5 + HHHHH(ID)*ER(ID,L)                                          
C                                                                               
      IF(ISOT(L)-4) 128,128,1254                                                
C                                                                               
C---- CALCULATE SIXTH CUMULANT INDEX PRODUCTS                                   
C     ELEMENTS OF TENSOR ARE MULTIPLIED BY 1.E6                                 
C                                                                               
 1254 DO 1255 ID=1,3                                                            
 1255 HHHHJ(ID) = HHHHH(ID) * HPR(ID) *0.1                                      
      HHHHJ( 4) = HHHHH(1 ) * HPR(2 ) *0.6                                      
      HHHHJ( 5) = HHHHH(1 ) * HPR(3 ) *0.6                                      
      HHHHJ( 6) = HHHHH(4 ) * HPR(2 ) *0.3                                      
      HHHHJ( 7) = HHHHH(4 ) * HPR(3 ) *0.6                                      
      HHHHJ( 8) = HHHHH(5 ) * HPR(3 ) *0.3                                      
      HHHHJ( 9) = HHHHH(6 ) * HPR(2 ) *0.2                                      
      HHHHJ(10) = HHHHH(6 ) * HPR(3 ) *0.6                                      
      HHHHJ(11) = HHHHH(7 ) * HPR(3 ) *0.3                                      
      HHHHJ(12) = HHHHH(8 ) * HPR(3 ) *0.2                                      
      HHHHJ(13) = HHHHH(9 ) * HPR(2 ) *0.15                                     
      HHHHJ(14) = HHHHH(9 ) * HPR(3 ) *0.6                                      
      HHHHJ(15) = HHHHH(10) * HPR(3 ) *0.3                                      
      HHHHJ(16) = HHHHH(11) * HPR(3 ) *0.2                                      
      HHHHJ(17) = HHHHH(12) * HPR(3 ) *0.15                                     
      HHHHJ(18) = HHHHH(2 ) * HPR(1 ) *0.6                                      
      HHHHJ(19) = HHHHH(18) * HPR(1 ) *0.6                                      
      HHHHJ(20) = HHHHH(19) * HPR(1 ) *0.6                                      
      HHHHJ(21) = HHHHH(20) * HPR(1 ) *0.6                                      
      HHHHJ(22) = HHHHH(21) * HPR(1 ) *0.6                                      
      HHHHJ(23) = HHHHH(3 ) * HPR(1 ) *0.6                                      
      HHHHJ(24) = HHHHH(2 ) * HPR(3 ) *0.6                                      
      HHHHJ(25) = HHHHH(18) * HPR(3 ) *0.3                                      
      HHHHJ(26) = HHHHH(19) * HPR(3 ) *0.2                                      
      HHHHJ(27) = HHHHH(20) * HPR(3 ) *0.15                                     
      HHHHJ(28) = HHHHH(3 ) * HPR(2 ) *0.6                                      
C                                                                               
      DO 1256 ID=1,28                                                           
 1256 CUM6 = CUM6 + HHHHJ(ID)*FR(ID,L)                                          
C                                                                               
C                                                                               
C---- CALCULATION OF THE TEMPERATURE FACTOR CORRESPONDING TO MODEL.             
C                                                                               
C                                                                               
  128 GOTO (129,1290,3290,3290),MODELL                                          
C                                                                               
C---- FOURIER TRANSFORM OF THE EDGEWORTH SERIES EXPANSION                       
C                                                                               
  129 EXPIJ = EXP(TEMP+CUM4)                                                    
  130 T     = TWOPI*ARG-CUM3                                                    
      COSIJ = COS (T)*EXPIJ                                                     
      SINIJ = SIN (T)*EXPIJ                                                     
      AAIND = AAIND+COSIJ                                                       
      IF(ICENT)140,150,140                                                      
  140 BBIND = BBIND+SINIJ                                                       
  150 IF(MFSEL)190,280,190                                                      
C                                                                               
C---- CALCULATE POSITIONAL DERIVATIVES                                          
C                                                                               
  190 DO 200 IAA=1,3                                                            
  200 DAX(IAA)=DAX(IAA)-HPR(IAA)*SINIJ                                          
      IF(ICENT)210,230,210                                                      
  210 DO 220 IAB=1,3                                                            
  220 DBX(IAB)=DBX(IAB)+HPR(IAB)*COSIJ                                          
  230 IF(ISOT(L))240,280,240                                                    
C                                                                               
C---- CALCULATE THERMAL DERIVATIVES FOR ANISOTROPIC ATOMS                       
C                                                                               
  240 DO 250 IAC=1,6                                                            
  250 DAB(IAC)=DAB(IAC)-HHJ(IAC)*COSIJ                                          
      IF(ICENT)260,270,260                                                      
  260 DO 265 ID=1,6                                                             
  265 DBB(ID)=DBB(ID)-HHJ(ID)*SINIJ                                             
  270 IF(ISOT(L)-1)280,280,271                                                  
C                                                                               
C---- CALCULATE THIRD CUM. DERIVS.                                              
C                                                                               
  271 DO 272 ID=1,10                                                            
  272 DAC(ID)=DAC(ID)+HHHJ(ID)*SINIJ                                            
      IF(ICENT)273,275,273                                                      
  273 DO 274 ID=1,10                                                            
  274 DBC(ID)=DBC(ID)-HHHJ(ID)*COSIJ                                            
  275 IF(ISOT(L)-2)3001,3001,276                                                
C                                                                               
C---- CALCULATE FOURTH CUMULANT DERIVS.                                         
C                                                                               
  276 DO 277 ID=1,15                                                            
  277 DAD(ID)=DAD(ID)+HHHH(ID)*COSIJ                                            
      IF(ICENT)278,3001,278                                                     
  278 DO 279 ID=1,15                                                            
  279 DBD(ID)=DBD(ID)+HHHH(ID)*SINIJ                                            
C                                                                               
 3001 IF(MODELL.EQ.3) GOTO 3800                                                 
      GOTO 280                                                                  
C                                                                               
C                                                                               
C---- FOURIER TRANSFORM OF THE GRAM-CHARLIER SERIES EXPANSION.                  
C                                                                               
 1290 EXPIJ = EXP(TEMP)                                                         
 1300 T     = TWOPI*ARG                                                         
      COSIJ = COS (T)*EXPIJ                                                     
      SINIJ = SIN (T)*EXPIJ                                                     
      TERM1 = 1.+FAK1*CUM4-FAK6*CUM6                                            
      TERM2 = FAK2*CUM3-FAK5*CUM5                                               
      AC    = COSIJ*TERM1+SINIJ*TERM2                                           
      BC    = SINIJ*TERM1-COSIJ*TERM2                                           
      AAIND = AAIND+AC                                                          
      IF(ICENT)1400,1500,1400                                                   
 1400 BBIND = BBIND+BC                                                          
C                                                                               
 1500 IF(MFSEL)1900,2800,1900                                                   
C                                                                               
C---- CALCULATE POSITIONAL DERIVATIVES                                          
C                                                                               
 1900 DO 2000 IAA=1,3                                                           
 2000 DAX(IAA)=DAX(IAA)+HPR(IAA)*(COSIJ*TERM2-SINIJ*TERM1)                      
      IF(ICENT)2100,2300,2100                                                   
 2100 DO 2200 IAB=1,3                                                           
 2200 DBX(IAB)=DBX(IAB)+HPR(IAB)*(COSIJ*TERM1+SINIJ*TERM2)                      
 2300 IF(ISOT(L))2400,2800,2400                                                 
C                                                                               
C---- CALCULATE THERMAL DERIVATIVES FOR ANISOTROPIC ATOMS                       
C                                                                               
 2400 DO 2500 IAC=1,6                                                           
 2500 DAB(IAC)=DAB(IAC)-HHJ(IAC)*AC                                             
      IF(ICENT)2600,2700,2600                                                   
 2600 DO 2650 ID=1,6                                                            
 2650 DBB(ID)=DBB(ID)-HHJ(ID)*BC                                                
 2700 IF(ISOT(L)-1)2800,2800,2710                                               
C                                                                               
C---- CALCULATE THIRD CUMULANT DERIVATIVES                                      
C                                                                               
 2710 DO 2720 ID=1,10                                                           
 2720 DAC(ID)=DAC(ID)+HHHJ(ID)*SINIJ*FAK2                                       
      IF(ICENT)2730,2750,2730                                                   
 2730 DO 2740 ID=1,10                                                           
 2740 DBC(ID)=DBC(ID)-HHHJ(ID)*COSIJ*FAK2                                       
 2750 IF(ISOT(L)-2)2800,2800,2760                                               
C                                                                               
C---- CALCULATE FOURTH CUMULANT DERIVATIVES                                     
C                                                                               
 2760 DO 2770 ID=1,15                                                           
 2770 DAD(ID)=DAD(ID)+HHHH(ID)*COSIJ*FAK1                                       
      IF(ICENT)2780,2779,2780                                                   
 2780 DO 2790 ID=1,15                                                           
 2790 DBD(ID)=DBD(ID)+HHHH(ID)*SINIJ*FAK1                                       
 2779 IF(ISOT(L)-3)2800,2800,2771                                               
C                                                                               
C---- CALCULATE FIFTH CUMULANT DERIVATIVES                                      
C                                                                               
 2771 DO 2778 ID=1,21                                                           
 2778 DAE(ID)=DAE(ID)-HHHHH(ID)*SINIJ*FAK5                                      
      IF(ICENT) 2777,2776,2777                                                  
 2777 DO 2775 ID=1,21                                                           
 2775 DBE(ID)=DBE(ID)+HHHHH(ID)*COSIJ*FAK5                                      
 2776 IF(ISOT(L)-4) 2800,2800,2772                                              
C                                                                               
C---- CALCULATE SIXTH CUMULANT DERIVATIVES                                      
C                                                                               
 2772 DO 2774 ID=1,28                                                           
 2774 DAF(ID)=DAF(ID)-HHHHJ(ID)*COSIJ*FAK6                                      
      IF(ICENT) 2773,2800,2773                                                  
 2773 DO 2769 ID=1,28                                                           
 2769 DBF(ID)=DBF(ID)-HHHHJ(ID)*SINIJ*FAK6                                      
 2800 IF(MODELL.EQ.4) GOTO 3800                                                 
      GOTO 280                                                                  
C                                                                               
C---- CALCULATES THE PRODUCTS OF *GERMIT* POLYNOMS AND CUMULANTS.               
C     (ONLY IF MODELL GREATER TWO.)                                             
C                                                                               
C---- FIRST CALCULATE THE INVERSE MATRIX OF BETA AND THE POLYNOMS.              
C                                                                               
 3290 KSYM=IA                                                                   
      CALL GERMIT(HPR,QTRANS,ALF(L),HHHJ,HHHH,ISOT(L),KSYM,T1,T2,               
     &  D3HERM)                                                                 
C                                                                               
C---- CALCULATE  G(KLM)*C(KLM)   AND   G(KLMN)*D(KLMN)                          
C                                                                               
      DO 3240 ID=1,10                                                           
 3240 CUM3=CUM3+HHHJ(ID)*CR(ID,L)                                               
C                                                                               
      IF(ISOT(L)-2) 3280,3280,3250                                              
 3250 DO 3251 ID=1,15                                                           
 3251 CUM4=CUM4+HHHH(ID)*DR(ID,L)                                               
C                                                                               
 3280 IF(MODELL.EQ.3) GOTO 129                                                  
      IF(MODELL.EQ.4) GOTO 1290                                                 
C                                                                               
C---- IN THIS PART THE DERIVATIVES WRT TO ALPHA ARE CALCULATED, IF              
C     MODELL GREATER 2  AND IF REFINEMENT OF ALPHA REQUIRED.                    
C     DERIVATIVES OF G(JKL) ARE ALREADY CALCULATED.                             
C                                                                               
 3800 IF(IALFA(L).EQ.0) GOTO 280                                                
      CALL ABALFA(D4HERM,ALF(L),T1,T2,ISOT(L))                                  
C                                                                               
C---- SUMMATION OVER ALL CUMULANTS                                              
C                                                                               
      DO 3801 ID=1,10                                                           
 3801 D3ALF=D3ALF+D3HERM(ID)*CR(ID,L)                                           
      IF(ISOT(L)-2) 3810,3810,3811                                              
 3811 DO 3802 ID=1,15                                                           
 3802 D4ALF=D4ALF+D4HERM(ID)*DR(ID,L)                                           
C                                                                               
 3810 IF(MODELL.EQ.3) GOTO 3900                                                 
C                                                                               
C---- DERIVATIVES WRT TO ALPHA IN CASE OF EXPANSION RELATED TO GRAM-            
C     CHARLIER SERIES EXPANSION.                                                
C                                                                               
      DAALF=COSIJ*FAK1*D4ALF+FAK2*SINIJ*D3ALF+DAALF                             
      IF(ICENT.NE.0) DBALF=SINIJ*FAK1*D4ALF-FAK2*COSIJ*D3ALF+DBALF              
      GOTO 280                                                                  
C                                                                               
C---- DERIVATIVE WRT TO ALPHA IN CASE OF AN EXPANSION RELATED TO AN             
C     EDGEWORTH EXPANSION.                                                      
C                                                                               
 3900 DAALF=SINIJ*D3ALF+COSIJ*D4ALF + DAALF                                     
      IF(ICENT.NE.0) DBALF=-COSIJ*D3ALF+SINIJ*D4ALF+DBALF                       
      GOTO 280                                                                  
C                                                                               
C                                                                               
C---- FROM HERE FOR ALL MODELS THE SAME.                                        
C                                                                               
  280 CONTINUE                                                                  
C                                                                               
  360 ITY1  = ISCAT(1,L)                                                        
      ITY2  = ISCAT(2,L)                                                        
      TEMP4 = TF*SITE(L)                                                        
C                                                                               
      IF(ITY2)370,380,370                                                       
C                                                                               
C---- EVALUATE SCATTERING FACTOR FOR MULTIPLE OCCUPANCY                         
C                                                                               
  370 OCCB =OCCUP(L)-OCCA(L)                                                    
      REAL1=SCAT(ITY1)                                                          
      REAL2=SCAT(ITY2)                                                          
      TEMP =OCCA(L)*REAL1+OCCB*REAL2                                            
      TEMP1=OCCA(L)*ANOM(2,ITY1)+OCCB*ANOM(2,ITY2)                              
      TEMP2=REAL1-REAL2                                                         
      TEMP3=ANOM(2,ITY1)-ANOM(2,ITY2)                                           
      GO TO 390                                                                 
C                                                                               
C---- EVALUATE SCATTERING FACTOR FOR SINGLE OCCUPANCY                           
C                                                                               
  380 REAL1=SCAT(ITY1)                                                          
      TEMP =OCCUP(L)*REAL1                                                      
      TEMP3=ANOM(2,ITY1)                                                        
      TEMP1=OCCUP(L)*TEMP3                                                      
      TEMP2=REAL1                                                               
 390  IF(MFSEL)392,396,392                                                      
C                                                                               
C---- EVALUATE DERIVATIVES WITH RESPECT TO OCCUPANCY                            
C                                                                               
  392 DADO(L)=TEMP4*(TEMP2*AAIND-TEMP3*BBIND)                                   
      DBDO(L)=TEMP4*(TEMP2*BBIND+TEMP3*AAIND)                                   
C                                                                               
      IF(IFSGN.EQ.0)GO TO 396                                                   
C                                                                               
C---- EVALUATE STRUCTURE FACTOR FOR ENANTIOMORPHIC ORIENTATION                  
C                                                                               
      AMINUS =TEMP4*(TEMP*AAIND+TEMP1*BBIND)+AMINUS                             
      BMINUS=TEMP4*(-TEMP*BBIND+TEMP1*AAIND)+BMINUS                             
C                                                                               
C---- EVALUATE INDIVIDUAL ATOM CONTRIBUTIONS                                    
C                                                                               
  396 AIND(L)=TEMP4*(TEMP*AAIND-TEMP1*BBIND)                                    
      BIND(L)=TEMP4*(TEMP*BBIND+TEMP1*AAIND)                                    
      ACALC=ACALC+AIND(L)                                                       
      BCALC=BCALC+BIND(L)                                                       
      IF(MFSEL)400,450,400                                                      
C                                                                               
C                                                                               
C---- MODIFY DERIVATIVES                                                        
C                                                                               
  400 TEMP5=TWOPI*TEMP4                                                         
      DO 410 IAG=1,3                                                            
      DADX(IAG,L)=TEMP5*(TEMP*DAX(IAG)-TEMP1*DBX(IAG))                          
  410 DBDX(IAG,L)=TEMP5*(TEMP*DBX(IAG)+TEMP1*DAX(IAG))                          
      IF(ISOT(L))430,420,430                                                    
C                                                                               
C---- DERIVATIVE WITH RESPECT TO ISOTROPIC TEMPERATURE FACTOR                   
C                                                                               
  420 DADB(1,L)=-RHO*AIND(L)                                                    
      DBDB(1,L)=-RHO*BIND(L)                                                    
      GO TO 450                                                                 
C                                                                               
C---- DERIVATIVES WITH RESPECT TO ANISOTROPIC TEMPERATURE FACTOR                
C                                                                               
  430 DO 440 IAH=1,6                                                            
      DADB(IAH,L)=TEMP4*(TEMP*DAB(IAH)-TEMP1*DBB(IAH))                          
  440 DBDB(IAH,L)=TEMP4*(TEMP*DBB(IAH)+TEMP1*DAB(IAH))                          
      IF(MODELL.LT.3 .OR. IALFA(L).EQ.0) GOTO 4420                              
C                                                                               
C---- DERIVATIVE WITH RESPECT TO ALPHA.                                         
C                                                                               
      DAALFA(L)=TEMP4*(TEMP*DAALF-TEMP1*DBALF)                                  
      DBALFA(L)=TEMP4*(TEMP*DBALF+TEMP1*DAALF)                                  
 4420 IF(ISOT(L)-1)450,450,441                                                  
C                                                                               
C---- DERIVATIVES WRT THE THIRD CUMULANT                                        
C                                                                               
  441 DO 442 ID=1,10                                                            
      DADC(ID,L)=TEMP4*(TEMP*DAC(ID)-TEMP1*DBC(ID))                             
  442 DBDC(ID,L)=TEMP4*(TEMP*DBC(ID)+TEMP1*DAC(ID))                             
      IF(ISOT(L)-2)450,450,443                                                  
C                                                                               
C---- DERIVATIVES WRT THE FOURTH CUMULANT                                       
C                                                                               
  443 DO 444 ID=1,15                                                            
      DADD(ID,L)=TEMP4*(TEMP*DAD(ID)-TEMP1*DBD(ID))                             
  444 DBDD(ID,L)=TEMP4*(TEMP*DBD(ID)+TEMP1*DAD(ID))                             
      IF(ISOT(L)-3)450,450,4441                                                 
C                                                                               
C---- DERIVATIVES WRT THE FIFTH CUMULANT                                        
C                                                                               
 4441 DO 4442 ID=1,21                                                           
      DADE(ID,L)=TEMP4*(TEMP*DAE(ID)-TEMP1*DBE(ID))                             
 4442 DBDE(ID,L)=TEMP4*(TEMP*DBE(ID)+TEMP1*DAE(ID))                             
      IF(ISOT(L)-4) 450,450,4443                                                
C                                                                               
C---- DERIVATIVES WRT THE SIXTH CUMULANT                                        
C                                                                               
 4443 DO 4444 ID=1,28                                                           
      DADF(ID,L)=TEMP4*(TEMP*DAF(ID)-TEMP1*DBF(ID))                             
 4444 DBDF(ID,L)=TEMP4*(TEMP*DBF(ID)+TEMP1*DAF(ID))                             
C                                                                               
  450 CONTINUE                                                                  
C                                                                               
  470 RETURN                                                                    
C                                                                               
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++            
C                                                                               
      ENTRY EXCALC(NONPOS)                                                      
C                                                                               
C     NONPOS = 1/2  FOR THIS REFLECTION IS NOT/IS AN IMAGINARY EX -             
C                   TINCTION FACTOR COMPUTED.                                   
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++            
C                                                                               
C---- CALCULATE EXTINCTION FACTOR OF THE ZACHARIASSEN MODEL.                    
C                                                                               
      IF(EXBETA.LE.0.0)EXBETA=1.0                                               
C                                                                               
      IF(IEXT)1010,1000,1010                                                    
C                                                                               
C---- ISOTROPIC FORM                                                            
C                                                                               
 1000 EXX=EXBETA*EXTIN(1)*FCALC                                                 
      DRDE(1)=1.0                                                               
C                                                                               
      RETURN                                                                    
C                                                                               
C---- ANISOTROPIC EXTINCTION - CALCULATE TRIG FUNCTIONS OF ANGLES               
C     (CAUTION|   THIS SYSTEM IS NOT A BUSING AND LEVY SYSTEM.  )               
C     (           SEE DESCRIBTION OF PROGRAMME BY  UDO ZUCKER   )               
C                                                                               
 1010 SCHI=SIN(EXANG(1))                                                        
      CCHI=COS(EXANG(1))                                                        
      SPHI=SIN(EXANG(2))                                                        
      CPHI=COS(EXANG(2))                                                        
      IF(IEXT-1)1030,1020,1030                                                  
C                                                                               
C---- TYPE 1 CORRECTION                                                         
C                                                                               
 1020 V(1)=SPHI*SCHI                                                            
      V(2)=-CPHI*SCHI                                                           
      V(3)=-CCHI                                                                
      GO TO 1040                                                                
C                                                                               
C---- TYPE 2 ANISOTROPIC                                                        
C                                                                               
 1030 SNU=SIN(EXANG(3))                                                         
      CNU=COS(EXANG(3))                                                         
      V(1)=-CCHI*SPHI*CNU-SNU*CPHI                                              
      V(2)=CCHI*CPHI*CNU-SNU*SPHI                                               
      V(3)=-SCHI*CNU                                                            
C                                                                               
C---- FORM PRODUCTS OF COEFFICIENTS                                             
C                                                                               
 1040 DO 1050 I=1,3                                                             
 1050 VV(I)=V(I)**2                                                             
      VV(4)=2.0*V(1)*V(2)                                                       
      VV(5)=2.0*V(1)*V(3)                                                       
      VV(6)=2.0*V(2)*V(3)                                                       
      SUM=0.0                                                                   
      DO 1060 I=1,6                                                             
 1060 SUM=SUM+VV(I)*EXTIN(I)                                                    
      IF(SUM.GT.1.E-20) GOTO 1061                                               
      SUM=1.0                                                                   
      NONPOS=2                                                                  
 1061 TEMP=-1.0/SQRT(SUM)                                                       
      EXX=EXBETA*FCALC*ABS(TEMP)                                                
C                                                                               
C---- CALCULATE DERIVATIVES OF THE FORM DR-STAR/DE-IJ                           
C                                                                               
      DO 1070 I=1,6                                                             
 1070 DRDE(I)=0.5*VV(I)*TEMP/SUM                                                
C                                                                               
      RETURN                                                                    
      END                                                                       
C                                                                               
C                                                                               
                      SUBROUTINE MATRIX (MODELL)                                
C                                                                               
C     THIS SUBROUTINE CALCULATES THE  DERIVATIVES  OF THE  RESIDUAL  FRO        
C     THE DERIVATIVES OF THE CALCULATED STRUCTURE FACTOR AND  ACCUMULATE        
C           THE SUMS REQUIRED IN THE NORMAL EQUATIONS MATRICES.                 
C                                                                               
C                     WRITTEN  BY  L. W. FINGER                                 
C                                                                               
C         MODIFIED  FOR  BECKER & COPPENS  EXTINCTION  BY UDO ZUCKER            
C                                                                               
C                                                                               
      DIMENSION IBDF(2),SAVE(60),KDEGN(2640),AA(2640),DADTL(60),DBDTL(60        
     & ), LDEGN(200), DADE(21,6), DBDE(21,6), DADF(28,6),  DBDF(28,6)           
      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)        
      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/LSMAT/B(200),IEXT,IPARA(1000),MRANK,KCYL,LABEL(3),JSCRF(10         
     & ), NRSF, SIGMA, MDEGN, MAXIMA, IALFA(60), RDAMP(200),IDAMP(200),         
     & A(2200)                                                                  
      COMMON/F/DELTA,FCALC,FOBS,FOMAX,FOMIN,IREJ,MREF,MREJ,NREJ,NUMBER(3        
     & 6),RDEN(36),RHO,RHOMN,RHOMX,RNUM(36),WT,EXBETA,EXX,EXY,GAMMAP            
     &,FISGN,IFSGN,EXTIN(6),KEXTRF(7),DELMAX,IANIGC,GAMMAS,IZACH,RNEUEX(        
     &7),NEUEXT,NDISTR,NIEXT                                                    
      COMMON/LS/DERIV(200)                                                      
      COMMON/TAPE/IN,INN,IOUT,ISC1,ISC2,ISC3,ISC4,ISC5,IPUN                     
      COMMON/COPPEN/YEXT,BDF(16),JUNK6                                          
	COMMON/LEVER/ILEVER
C                                                                               
      EQUIVALENCE (SAVE,AIND),(KDEGN,AA,DADC(1,37)),(IBDF,BDF),(LDEGN,A)        
      EQUIVALENCE (DADD(1,41),DADE),(DADD(7,49),DADF),(DBDD(1,41),DBDE)         
     & ,(DBDD(7,49),DBDF),(DADTL,AIND),(DBDTL,BIND)                             
C                                                                               
      I2(K) = MOD(K/256**2,256)                                                 
      I3(K) = MOD(K/256,256)                                                    
      I4(K) = MOD(K,256)                                                        
C                                                                               
C                                                                               
      IF(IZACH.EQ.0)GOTO 102                                                    
C                                                                               
C---- EXTINCTION CORRECTION OF ZACHARIASSEN                                     
C                                                                               
C     CONS  MODIFICATIONS  OF DERIVATIVES  DUE  TO EXTINCTION.                  
C     CONS2 COEFFICIENT FOR DERIVATIVE OF F WTH TO EXTINCTION-                  
C           PARAMETERS.                                                         
C     IN CASE OF BECKER&COPPENS EXTINCTION CONS IS VORFAK AND                   
C     CONS WILL BE CALCULATED IN SUBROUTINE EXTIC.                              
C                                                                               
      CONS=-0.5*(EXY+EXY**3)/FCALC                                              
      CONS2=0.25*FCALC**3*EXBETA*EXY                                            
C                                                                               
      IF(MREF)80,90,80                                                          
C                                                                               
C---- REFINEMENT ON F*F                                                         
C                                                                               
   80 CONS=2.0*FCALC*CONS                                                       
      CONS2=2.0*FCALC*CONS2                                                     
C                                                                               
C---- INCLUDE DERIVATIVES, IF PARAMETER SHOULD BE REFINED.                      
C                                                                               
   90 J=1                                                                       
      IF(IEXT.NE.0)J=6                                                          
      DO 100 I=1,J                                                              
      IJ=KEXTRF(I)                                                              
      IF(IJ.NE.0)DERIV(IJ)=CONS2*DRDE(I)                                        
  100 CONTINUE                                                                  
      GOTO 101                                                                  
  102 CONTINUE                                                                  
C                                                                               
C---- EXTINCTION CORRECTION OF BECKER AND COPPENS                               
C                                                                               
      CONS=-VORFAK                                                              
      DO 2000 I=1,IPAREX(IBDF(1)+1)                                             
      IJ=KEXTRF(I)                                                              
 2000 IF(IJ.NE.0)DERIV(IJ)=-DRDE(I)                                             
C---- NEXT TO LINES ADDED 20.02.85 W.F.K.
C     IF(MREF) 2001,101,2001
C2001 CONS = 2.0*FCALC*CONS
  101 ACALC=ACALC*CONS                                                          
      BCALC=BCALC*CONS                                                          
C                                                                               
C---- DERIVATIVE WRT  ALPHA (MODELL > 2)                                        
C                                                                               
      IF(MODELL.LT.3) GOTO 104                                                  
      DO 1400 I=1,NATOM                                                         
      IJ=IALFA(I)                                                               
      IF(IJ.LE.0) GOTO 1400                                                     
      DERIV(IJ)=ACALC*DAALFA(I)+BCALC*DBALFA(I)                                 
 1400 CONTINUE                                                                  
C                                                                               
C---- SET ALL SCALING FACTORS TO ZERO, WHICH ARE NOT VALID FOR                  
C     THIS HKL -VALUE.                                                          
C                                                                               
  104 IF(NRSF.EQ.0)GO TO 140                                                    
      DO 120 I=1,NRSF                                                           
      IJ=JSCRF(I)                                                               
      IF(IJ.NE.0)DERIV(IJ)=0.                                                   
  120 CONTINUE                                                                  
C                                                                               
C---- DERIVATIVE WITH RESPECT TO SCALE FACTOR (FOR CURRENT REFLEC.)             
C                                                                               
      IF(JSCRF(ISCL).EQ.0) GO TO 140                                            
      J=JSCRF(ISCL)                                                             
      DERIV(J)=-FCALC                                                           
      IF(MREF)130,140,130                                                       
C                                                                               
C---- REFINEMENT ON F**2                                                        
C                                                                               
  130 DERIV(J)=-2.0*FCALC **2                                        
  140 CONTINUE                                                                  
      II=1                                                                      
C                                                                               
      IF(NPAR.EQ.0)GO TO 149                                                    
C                                                                               
C---- SPECIAL PARAMETERS                                                        
C                                                                               
      DO 145 I=1,NPAR                                                           
      IJ=IRPAR(I)                                                               
      IF(IJ.LE.0)GO TO 145                                                      
      DERIV(IJ)=0.                                                              
  141 IF(I2(LDEGN(II))-IJ)1415,142,145                                          
 1415 II=II+2                                                                   
      IF(II-MDEGN)141,141,149                                                   
  142 CALL MODIFY(DERIV(IJ),I3(LDEGN(II)),I4(LDEGN(II)),A(II+1))                
      GO TO 1415                                                                
  145 CONTINUE                                                                  
  149 IK=0                                                                      
C                                                                               
C---- CONVENTIONAL CRYSTALLOGRAPHIC PARAMETERS                                  
C                                                                               
                             DO 240 J=1,NATOM                                   
      IK=IK+1                                                                   
      IJ=IPARA(IK)                                                              
      IF(IJ)150,160,150                                                         
C                                                                               
C---- DERIVATIVE WRT OCCUPANCY                                                  
C                                                                               
  150 DERIV(IJ)=ACALC*DADO(J)+BCALC*DBDO(J)                                     
  152 IF(II.GT.MDEGN)GO TO 160                                                  
      IF(I2(LDEGN(II))-IJ)154,156,160                                           
  154 II=II+2                                                                   
      GO TO 152                                                                 
  156 CALL MODIFY(DERIV(IJ),I3(LDEGN(II)),I4(LDEGN(II)),A(II+1))                
      GO TO 154                                                                 
  160                        DO 180 K=1,3                                       
      IK=IK+1                                                                   
      IJ=IPARA(IK)                                                              
      IF(IJ)170,180,170                                                         
C                                                                               
C---- DERIVATIVE WRT POSITION                                                   
C                                                                               
  170 DERIV(IJ)=ACALC*DADX(K,J)+BCALC*DBDX(K,J)                                 
  172 IF(II.GT.MDEGN)GO TO 180                                                  
      IF(I2(LDEGN(II))-IJ)174,176,180                                           
  174 II=II+2                                                                   
      GO TO 172                                                                 
  176 CALL MODIFY(DERIV(IJ),I3(LDEGN(II)),I4(LDEGN(II)),A(II+1))                
      GO TO 174                                                                 
  180 CONTINUE                                                                  
      IF(ISOT(J))210,190,210                                                    
  190 IK=IK+1                                                                   
      IJ=IPARA(IK)                                                              
      IF(IJ)200,240,200                                                         
C                                                                               
C---- DERIVATIVE WRT ISOTROPIC TEMPERATURE FACTOR                               
C                                                                               
  200 DERIV(IJ)=ACALC*DADB(1,J)+BCALC*DBDB(1,J)                                 
  202 IF(II.GT.MDEGN)GO TO 240                                                  
      IF(I2(LDEGN(II))-IJ)204,206,240                                           
  204 II=II+2                                                                   
      GO TO 202                                                                 
  206 CALL MODIFY(DERIV(IJ),I3(LDEGN(II)),I4(LDEGN(II)),A(II+1))                
      GO TO 204                                                                 
  210                        DO 230 L=1,6                                       
      IK=IK+1                                                                   
      IJ=IPARA(IK)                                                              
      IF(IJ)220,230,220                                                         
C                                                                               
C---- DERIVATIVE WRT ANISOTROPIC TEMPERATURE FACTOR                             
C                                                                               
  220 DERIV(IJ)=ACALC*DADB(L,J)+BCALC*DBDB(L,J)                                 
  222 IF(II.GT.MDEGN)GO TO 230                                                  
      IF(I2(LDEGN(II))-IJ)224,226,230                                           
  224 II=II+2                                                                   
      GO TO 222                                                                 
  226 CALL MODIFY(DERIV(IJ),I3(LDEGN(II)),I4(LDEGN(II)),A(II+1))                
      GO TO 224                                                                 
  230 CONTINUE                                                                  
      IF(ISOT(J)-1)240,240,231                                                  
C                                                                               
C---- DERIVATIVES WRT TO THE THIRD CUMULANT                                     
C                                                                               
  231 DO 233 L=1,10                                                             
      IK=IK+1                                                                   
      IJ=IPARA(IK)                                                              
      IF(IJ)232,233,232                                                         
  232 DERIV(IJ)=ACALC*DADC(L,J)+BCALC*DBDC(L,J)                                 
 2322 IF(II.GT.MDEGN)GO TO 233                                                  
      IF(I2(LDEGN(II))-IJ)2324,2326,233                                         
 2324 II=II+2                                                                   
      GO TO 2322                                                                
 2326 CALL MODIFY(DERIV(IJ),I3(LDEGN(II)),I4(LDEGN(II)),A(II+1))                
      GO TO 2324                                                                
  233 CONTINUE                                                                  
      IF(ISOT(J)-2)240,240,234                                                  
C                                                                               
C---- DERIVATIVES WRT TO THE FOURTH CUMULANT                                    
C                                                                               
  234 DO 236 L=1,15                                                             
      IK=IK+1                                                                   
      IJ=IPARA(IK)                                                              
      IF(IJ)235,236,235                                                         
  235 DERIV(IJ)=ACALC*DADD(L,J)+BCALC*DBDD(L,J)                                 
 2352 IF(II.GT.MDEGN)GO TO 236                                                  
      IF(I2(LDEGN(II))-IJ)2354,2356,236                                         
 2354 II=II+2                                                                   
      GO TO 2352                                                                
 2356 CALL MODIFY(DERIV(IJ),I3(LDEGN(II)),I4(LDEGN(II)),A(II+1))                
      GO TO 2354                                                                
  236 CONTINUE                                                                  
      IF(ISOT(J)-3) 240,240,331                                                 
C                                                                               
C---- DERIVATIVES WRT TO THE FIFTH CUMULANT                                     
C                                                                               
  331 DO 333 L=1,21                                                             
      IK=IK+1                                                                   
      IJ=IPARA(IK)                                                              
      IF(IJ) 332,333,332                                                        
  332 DERIV(IJ)=ACALC*DADE(L,J)+BCALC*DBDE(L,J)                                 
 3322 IF(II.GT.MDEGN) GOTO 333                                                  
      IF(I2(LDEGN(II))-IJ)3324,3326,333                                         
 3324 II=II+2                                                                   
      GOTO 3322                                                                 
 3326 CALL MODIFY(DERIV(IJ),I3(LDEGN(II)),I4(LDEGN(II)),A(II+1))                
      GOTO 3324                                                                 
  333 CONTINUE                                                                  
      IF(ISOT(J)-4) 240,240,334                                                 
C                                                                               
C---- DERIVATIVES WRT TO THE SIXTH CUMULANT                                     
C                                                                               
  334 DO 336 L=1,28                                                             
      IK=IK+1                                                                   
      IJ=IPARA(IK)                                                              
      IF(IJ) 335,336,335                                                        
  335 DERIV(IJ)=ACALC*DADF(L,J)+BCALC*DBDF(L,J)                                 
 3352 IF(II.GT.MDEGN) GOTO 336                                                  
      IF(I2(LDEGN(II))-IJ) 3354,3356,336                                        
 3354 II=II+2                                                                   
      GOTO 3352                                                                 
 3356 CALL MODIFY(DERIV(IJ),I3(LDEGN(II)),I4(LDEGN(II)),A(II+1))                
      GOTO 3354                                                                 
  336 CONTINUE                                                                  
C                                                                               
  240 CONTINUE
C                                                                               
C *********************************************************************         
C ***********************  FORM MATRICES ******************************         
C *********************************************************************         
C                                                                               
      INDEX=MDEGN+1                                                             
      DO 280 M=1,MRANK                                                          
      TEMP=DERIV(M)*WT                                                          
      IF(TEMP)260,250,260                                                       
 250  INDEX=INDEX+MRANK-M+1                                                     
      GO TO 280                                                                 
 260  DO 270 N=M,MRANK                                                          
      A(INDEX)=A(INDEX)+TEMP*DERIV(N)                                           
 270  INDEX=INDEX+1                                                             
      B(M)=B(M)-DELTA*TEMP                                                      
  280 CONTINUE                                                                  
C
C	IF ILEVER=1 WRITE DERIVATES & WEIGHTS OUT TO UNIT 19
C
	IF(ILEVER.NE.0)THEN
	 WRITE(19) IH,WT,(DERIV(I),I=1,MRANK)
	ENDIF
C                                                                               
      RETURN                                                                    
      END                                                                       
C                                                                               
                 SUBROUTINE MODIFY(DERIV,NA,NP,TEMP)                            
C                                                                               
C     THIS ROUTINE APPLIES THE EQUATIONS                                        
C                                                                               
C                                                                               
C     D(A')           D(A)     D(P )                                            
C     ----- =  SUM ( ------ * ----M- )                                          
C     D(P )     M     D(P )    D(P )                                            
C        N               M        N                                             
C                                                                               
C                                                                               
C     D(B')           D(B)     D(P )                                            
C     ----- =  SUM ( ------ * ----M- )                                          
C     D(P )     M     D(P )    D(P )                                            
C        N               M        N                                             
C                                                                               
C     TO THE CALCULATED DERIVATIVES. A AND B ARE THE REAL AND  IMAGINARY        
C     PORTIONS OF THE CALCULATED STRUCTURE FACTOR, P(N) IS THE  INDEPEN-        
C     DENT PARAMETER, P(M) REPRESENTS A  DEPENDENT  PARAMETER, THE  SUM-        
C     MATION IS OVER ALL THE DEPENDENCIES, THE UNPRIMED DERIVATIVES  DE-        
C     NOTES THE ONE CALCULATED WITHOUT REGARD FOR DEPENDENCIES  AND  THE        
C     PRIMED DERIVATIVE IS  NEEDED TO  FORM THE  LEASTSQUARES  MATRICES.        
C                                                                               
C                      WRITTEN  BY  L. W. FINGER                                
C                                                                               
C     ERROR CORRECTION AND ADDING OF  FIFTH  AND  SIXTH  ORDER  CUMULANT        
C              DERIVATIVES  BY  UDO ZUCKER,  FEBRUARY 1980.                     
C                                                                               
      DIMENSION  SAVE(60), KDEGN(2640), AA(2640), DADE(21,6),DBDE(21,6),        
     & DADF(28,6),DBDF(28,6)                                                    
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)                     
C                                                                               
      EQUIVALENCE  (SAVE,AIND),(KDEGN,AA,DADC(1,37)),(DADD(1,41), DADE),        
     & (DADD(7,49),DADF),(DBDD(1,41),DBDE),(DBDD(7,49),DBDF)                    
C                                                                               
      IF(NP.GE.2)GO TO 100                                                      
C---- OCCUPANCIES.                                                              
      DERIV=DERIV+TEMP*(ACALC*DADO(NA)+BCALC*DBDO(NA))                          
      RETURN                                                                    
 100  IF(NP.GE.5)GO TO 200                                                      
C---- COORDINATES.                                                              
      DERIV=DERIV+TEMP*(ACALC*DADX(NP-1,NA)+BCALC*DBDX(NP-1,NA))                
      RETURN                                                                    
 200  IF(NP.GE.11)GO TO 300                                                     
C---- ANISOTROPIC THERMAL PARAMETERS.                                           
      DERIV=DERIV+TEMP*(ACALC*DADB(NP-4,NA)+BCALC*DBDB(NP-4,NA))                
      RETURN                                                                    
 300  IF(NP.GE.21)GO TO 400                                                     
C---- THERMAL CUMULANTS OF THIRD ORDER.                                         
      DERIV=DERIV+TEMP*(ACALC*DADC(NP-10,NA)+BCALC*DBDC(NP-10,NA))              
      RETURN                                                                    
 400  IF(NP.GE.36) GOTO 500                                                     
C---- T. C. O. FOURTH ORDER.                                                    
      DERIV=DERIV+TEMP*(ACALC*DADD(NP-20,NA)+BCALC*DBDD(NP-20,NA))              
      RETURN                                                                    
 500  IF(NP.GE.57) GOTO 600                                                     
C---- T. C. O. FIFTH ORDER.                                                     
      DERIV=DERIV+TEMP*(ACALC*DADE(NP-35,NA)+BCALC*DBDE(NP-35,NA))              
      RETURN                                                                    
C---- T. C. O. SIXTH ORDER.                                                     
 600  DERIV=DERIV+TEMP*(ACALC*DADF(NP-56,NA)+BCALC*DBDF(NP-56,NA))              
C                                                                               
      RETURN                                                                    
      END                                                                       
C                                                                               
                    SUBROUTINE RCALC(SCFF,NWT,GOF)                              
C                                                                               
C     THIS ROUTINE ACCUMULATES THE TERMS NEEDED TO CALCUALATE  THE  DIS-        
C     CREPANCY FACTORS WHICH ARE COMPUTED FOR ALL REFLECTIONS AND SEPA -        
C     RATELY FOR THOSE DATA USED IN CALCULATION.     A WEIGHTED RESIDUAL        
C                                                                               
C                                 2                                             
C              SUM (W*(FOBS-FCALC) )  1/2                                       
C     R(W) = ( --------------------- )       FOR REFIENEMENT ON FOBS            
C                            2                                                  
C                 SUM (W*FOBS )                                                 
C                                                                               
C                                                                               
C     OR                                                                        
C                                                                               
C                                   2                                           
C              SUM (W*(IOBS -ICALC ) )  1/2                                     
C     R(W) = ( ----------------------- )     FOR REFINEMENT ON IOBS             
C                              2                                                
C                   SUM (W*IOBS )                                               
C                                                                               
C     WILL BE COMPUTED.                                                         
C                                                                               
C     IN ADDITION THE CONVENTIONAL DISCREPANCY INDEX                            
C                                                                               
C              SUM(/FOBS-FCALC/)                                                
C     R    =   -----------------                                                
C                 SUM (FOBS)                                                    
C                                                                               
C     IS COMPUTED.                                                              
C                                                                               
C     IN ADDITION, THE RANGES OF VARIATION OF FOBS AND SIN(THETA)/LAMBDA        
C     ARE EACH DIVIDED INTO EIGHT SUBRANGES AND  THE  SUMS  NEEDED   FOR        
C     R(INT) ARE ACCUMULATED. IF THE WEIGHTING SCHEME  IS  CORRECT,  THE        
C                                        2                                      
C     EXPRESSION (1/S)*SUM(W*(FOBS-FCALC)  WHERE THE SUM IS TAKEN OVER A        
C     SUBSET OF THE FULL DATA SET CONTAINING  S  REFLECTIONS,  SHOULD BE        
C     INDEPENDENT OF THE COMPOSITION OF THE SUBSET.    IF WE DIVIDE THIS        
C                                                 2                             
C     EXPRESSION BY THE CONSTANT  (1/T)*SUM(W*FOBS ) ,    WHERE T IS THE        
C     TOTAL NUMBER OF REFLECTIONS IN THE DATA SET AND THE  SUM IS  TAKEN        
C     OVER THE ENTIRE SET, AND TAKE THE SQUARE ROOT OF THE QUOTIENT, THE        
C     RESULT SHOULD ALSO BE INDEPENDENT OF THE COMPOSITION  OF  THE SUB-        
C     SET. THIS NUMBER R(INT), WHICH CAN BE CONSIDERED THE AVERAGE CON -        
C     TRIBUTION TO R(W) FOR REFLECTIONS IN THE SUBSET IS PRINTED IN  THE        
C     OUTPUT FOR EACH SUBSET.                                                   
C                                                                               
C                        WRITTEN  BY L. W. FINGER                               
C                                                                               
C           THIS  ROUTINE REVISED MARCH 10, 1970   (L.W. FINGER)                
C                                                                               
C         MODIFICATIONS OCTOBER 78 - DECEMBER 80  BY    UDO ZUCKER              
C                                                                               
C     INPUT VARIABLES                                                           
C     ---------------                                                           
C                                                                               
C     SCFF = SCALEFACTOR   OF   PROCESSED   REFLECTION.                         
C     NWT  = FLAG FOR REFINEMENT WITH WEIGHTS OR UNIT-WEIGHTS.                  
C     GOF  = SEE BELOW.                                                         
C                                                                               
      COMMON/F/DELTA,FCALC,FOBS,FOMAX,FOMIN,IREJ,MREF,MREJ,NREJ,NUMBER(3        
     & 6),RDEN(36), RHO, RHOMN, RHOMX,RNUM(36),WT,EXBETA,EXX,EXY,GAMMAP,        
     & FISGN,IFSGN,EXTIN(6),KEXTRF(7),DELMAX,IANIGC,GAMMAS,IZACH,RNEUEX(        
     & 7),NEUEXT,NDISTR,NIEXT                                                   
C                                                                               
C                                                                               
C---- CALCULATE RANGE OF FOBS AND RHO                                           
C                                                                               
      SCF=SCFF                                                                  
      IRNG1=(SCF*FOBS-FOMIN)/(FOMAX-FOMIN+0.00001)*8.0+1.0                      
      IF(IRNG1)98,98,99                                                         
   98 IRNG1=1                                                                   
      GO TO 110                                                                 
   99 IF(IRNG1-8)110,110,100                                                    
  100 IRNG1=8                                                                   
  110 IRNG2=    (RHO-RHOMN)/(RHOMX-RHOMN+0.00001)*8.0+1.0                       
      IF(IRNG2)102,102,104                                                      
  102 IRNG2=1                                                                   
      GO TO 130                                                                 
  104 IF(IRNG2-8)130,130,120                                                    
  120 IRNG2=8                                                                   
C                                                                               
C---- CALCULATE NUMERATOR AND DENOMINATOR FOR WEIGHTED R                        
C                                                                               
C---- REFINEMENT WITH UNIT WEIGHTS ? (YES   NWT=2)                              
  130 WT=SCF**2/WT**2                                                           
      IF(NWT.EQ.2) WT=SCF**2                                                    
      IF(MREF)150,140,150                                                       
C                                                                               
C---- REFINE ON F                                                               
C                                                                               
  140 DELTA=FOBS-FCALC                                                          
      FDEN=WT*FOBS**2                                                           
      GO TO 160                                                                 
C                                                                               
C---- REFINEMENT OF F**2 (GOF = SUM(W*(FOBS-FCALC)**2  FOR GODNESS OF FI        
C                                                                               
  150 TEMP = WT*(FOBS-FCALC)**2                                                 
      IF(NREJ.EQ.0) GOF  = GOF + TEMP                                                         
      WT=0.25*WT/FOBS**2                                                        
      DELTA=FOBS**2-FCALC**2                                                    
      FDEN=FOBS**4*WT                                                           
  160 FNUM=WT*DELTA**2                                                          
      DELF=ABS(FOBS-FCALC)                                                      
                             DO 200 I=1,19,18                                   
      RDEN(I)=RDEN(I)+FDEN                                                      
      RDEN(I+1)=RDEN(I+1)+FOBS                                                  
      RNUM(I)=RNUM(I)+FNUM                                                      
      RNUM(I+1)=RNUM(I+1)+DELF                                                  
      NUMBER(I)=NUMBER(I)+1                                                     
      J=I+IRNG1+1                                                               
      II=0                                                                      
  170 RNUM(J)=RNUM(J)+FNUM                                                      
      RDEN(J)=RDEN(J)+FDEN                                                      
      NUMBER(J)=NUMBER(J)+1                                                     
      II=II+1                                                                   
      IF(II-1)190,180,190                                                       
  180 J=I+IRNG2+9                                                               
      GO TO 170                                                                 
  190 IF(NREJ)210,200,210                                                       
  200 CONTINUE                                                                  
C                                                                               
  210 RETURN                                                                    
      END                                                                       
                                                                                
                                                                                
                                                                                
                                                                                
                                                                                
                                                                                
                                                                                
