      SUBROUTINE MATRIX(ILVG)
C
C             ROUTINE TO FORM LEAST SQUARES MATRIX
C
      INCLUDE 'RFINE.PRM'
      INCLUDE 'A.INC'
      INCLUDE 'XPAR.INC'
      INCLUDE 'D.INC'
      INCLUDE 'F.INC'
      INCLUDE 'LS.INC'
      INCLUDE 'LSMAT.INC'
      INCLUDE 'CNSTRS.INC'
      INCLUDE 'TWIN.INC'
      INCLUDE 'EXTN.INC'
      INCLUDE 'ABSORB.INC'
      DIMENSION LDEGN(200),DERIV(QQPRM),ADERIV(QQPRM+2)
      EQUIVALENCE (ADERIV,DAB)
     1,(ADERIV(1),APWT),(ADERIV(2),APDELT),(ADERIV(3),DERIV(1))
C
      IF(IZACH.NE.0)THEN
C
C---- EXTINCTION CORRECTION OF ZACHARIASSEN
C
C     CONS  MODIFICATIONS  OF DERIVATIVES  DUE  TO EXTINCTION.
C     CONS2 COEFFICIENT FOR DERIVATIVE OF F WTH EXTINCTION-
C           PARAMETERS.
C     IN CASE OF BECKER & COPPENS EXTINCTION CONS IS VORFAK AND
C     WILL BE CALCULATED IN SUBROUTINE EXTIC.
C
        FTEMP=FCALC
        IF(FTEMP.EQ.0.0)FTEMP=1.0E6
        CONS=-0.5*(EXY+EXY**3)*ABSFCT/FTEMP
        CONS2=0.25*FCALC**3*EXBETA*EXY/ABSFCT
C
        IF(MREF.NE.0) THEN
C---- REFINEMENT ON F*F
          CONS=2.0*FCALC*CONS
          CONS2=2.0*FCALC*CONS2
        END IF
C---- INCLUDE DERIVATIVES, IF PARAMETER SHOULD BE REFINED.
        IJ=KEXTRF(1)
        IF(IJ.NE.0)DERIV(IJ)=CONS2*DRDE(1)
      ELSE
C
C---- EXTINCTION CORRECTION OF BECKER AND COPPENS
C
        CONS=-VORFAK*ABSFCT
	  DO 330 I=1,7
            IJ=KEXTRF(I)
	    IF(IJ.NE.0)DERIV(IJ)=-DRDE(I)
330	  CONTINUE
	END IF
C --- ANALYTICAL ABSORPTION CORRECTION
      IF(NABSRN.GT.0)THEN
        N=NABSRN*(NABSRN+2)
        TEMP=FCALC*0.5/ABSFCT
        DO 348 I=1,N
          IJ=KABSRF(I)
          IF(IJ.NE.0)DERIV(IJ)=TEMP*DABSCF(I)
  348   CONTINUE
      END IF
      IF(NRSF.NE.0)THEN
        DO 350 I=1,NRSF
          IJ=JSCRF(I)
          IF(IJ.NE.0)DERIV(IJ)=0.
  350   CONTINUE
C DERIVATIVE WITH RESPECT TO SCALE FACTOR
        IF(JSCRF(ISCL).NE.0) THEN
          J=JSCRF(ISCL)
          IF(MREF.NE.0)THEN
C REFINEMENT ON F**2
  	    DERIV(J)=-2.0*FCALC**2/SCALE(ISCL)
	  ELSE
C REFINEMENT ON F
            DERIV(J)=-FCALC/SCALE(ISCL)
	  END IF
        END IF
      END IF
      ACALC=ACALC*CONS
      BCALC=BCALC*CONS
      II=1
      IF(NPAR.EQ.0)GO TO 390
      DO 380 I=1,NPAR
      IJ=IRPAR(I)
      IF(IJ.LE.0)GO TO 380
      DERIV(IJ)=0.
  380 CONTINUE
C REFINEMENT OF TWIN VOLUMES
  390	DO 400 I=1,NT
	  IJ=JTWNMT(I)
	  IF(IJ.GT.0)DERIV(IJ)=CONS*((ATWN(I)**2+BTWN(I)**2)/2.0
     1       -TFPAR*FPAREN)
  400	CONTINUE
	IK=0
                             DO 560 J=1,NATOM
      IK=IK+1
      IJ=IPARA(IK)
      IF(IJ)410,420,410
C             DERIVATIVE WRT OCCUPANCY
  410 DERIV(IJ)=ACALC*DADO(J)+BCALC*DBDO(J)
  420                        DO 440 K=1,3
      IK=IK+1
      IJ=IPARA(IK)
      IF(IJ)430,440,430
C             DERIVATIVE WRT POSITION
  430 DERIV(IJ)=ACALC*DADX(K,J)+BCALC*DBDX(K,J)
  440 CONTINUE
      IF(ISOT(J))470,450,470
  450 IK=IK+1
      IJ=IPARA(IK)
      IF(IJ)460,560,460
C             DERIVATIVE WRT ISOTROPIC TEMPERATURE FACTOR
  460 DERIV(IJ)=ACALC*DADB(1,J)+BCALC*DBDB(1,J)
      GO TO 560
  470                        DO 490 L=1,6
      IK=IK+1
      IJ=IPARA(IK)
      IF(IJ)480,490,480
C             DERIVATIVE WRT ANISOTROPIC TEMPERATURE FACTOR
  480 DERIV(IJ)=ACALC*DADB(L,J)+BCALC*DBDB(L,J)
  490 CONTINUE
      IF(ISOT(J)-1)560,560,500
C        DERIVATIVES WRT TO THE THIRD CUMULANT
  500 DO 520 L=1,10
      IK=IK+1
      IJ=IPARA(IK)
      IF(IJ)510,520,510
  510 DERIV(IJ)=ACALC*DADC(L,J)+BCALC*DBDC(L,J)
  520 CONTINUE
      IF(ISOT(J)-2)560,560,530
C        DERIVATIVES WRT TO THE FOURTH CUMULANT
  530 DO 550 L=1,15
      IK=IK+1
      IJ=IPARA(IK)
      IF(IJ)540,550,540
  540 DERIV(IJ)=ACALC*DADD(L,J)+BCALC*DBDD(L,J)
  550 CONTINUE
  560 CONTINUE
C
C MODIFY DERIVATIVES
C
	IF(MDEGN.GT.0)CALL MODIFY
C WRITE DATA FOR LEVERAGE CALCULATION
	IF(ILVG.NE.0)WRITE(50)IH,WT,(DERIV(I),I=1,MRANK)
C             FORM MATRICES
      INDEX=MDEGN+1
      DELTA=DELTA*BIW
                             DO 610 M=1,MRANK
      TEMP=DERIV(M)*WT
      IF(TEMP)590,580,590
  580 INDEX=INDEX+MRANK-M+1
      GO TO 610
  590                        DO 600 N=M,MRANK
      A(INDEX)=A(INDEX)+TEMP*DERIV(N)
  600 INDEX=INDEX+1
      B(M)=B(M)-DELTA*TEMP
  610 CONTINUE
      RETURN
      END
