SUBROUTINE POLAR
C.. 6/13/91
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
C**********************************************************************
C
C   POLAR SETS UP THE CALCULATION OF THE MOLECULAR ELECTRIC RESPONSE
C   PROPERTIES BY FFHPOL.
C
C**********************************************************************
      CHARACTER*2 ELEMNT
      DIMENSION X1(MAXORB,MAXORB), X2(MAXORB,MAXORB), X3(MAXORB,MAXORB)
      DIMENSION X4(MAXORB,MAXORB), X5(MAXORB,MAXORB), X6(MAXORB,MAXORB)
      DIMENSION X7(MAXORB,MAXORB), X8(MAXORB,MAXORB), X9(MAXORB,MAXORB)
      DIMENSION X10(MAXORB,MAXORB),X11(MAXORB,MAXORB),X12(MAXORB,MAXORB)
      DIMENSION X13(MAXORB,MAXORB)
      COMMON /WORK1/ X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,
     1               XDUMY(21*MPACK-10*MAXORB*MAXORB)
      COMMON /WORK3/ X11,X12,XDUMY1(4*MPACK-2*MAXORB*MAXORB)
      COMMON /SCRACH/ X13, XDUMY2(MAXPAR*MAXPAR-MAXORB*MAXORB)
C..
      COMMON /TITLES/ KOMENT,TITLE
      COMMON /POLVOL/ POLVOL(107)
      COMMON /KEYWRD/ KEYWRD
      COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
     1                NA(NUMATM),NB(NUMATM),NC(NUMATM)
      COMMON /GEOVAR/ NVAR,LOC(2,MAXPAR),IDUMY,XPARAM(MAXPAR)
      COMMON /SCFTYP/ EMIN,LIMSCF
      COMMON /TIMCOM/ TIME0
      COMMON /ELEMTS/ ELEMNT(107)
      COMMON /CORE  / CORE(107)
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /GEOSYM/ NDEP, LOCPAR(MAXPAR), IDEPFN(MAXPAR),
     1                    LOCDEP(MAXPAR)
      COMMON /GEOM  / GEO(3,NUMATM), COORD(3,NUMATM)
      COMMON /LAST  / LAST
      COMMON /EULER / TVEC(3,3),IDTVEC
      COMMON /OMVAL/ OMEGA
      COMMON /CHANEL/ IFILES(30)
      EQUIVALENCE(IW,IFILES(6)), (IR,IFILES(5))
      DIMENSION GRAD(MAXPAR),ROTVEC(3,3), TEMPV(3,3), DATAEV(10)
     1,VALUE(40)
      CHARACTER  KEYWRD*241, TYPE*7, KOMENT*81, TITLE*81
      CHARACTER  POLKEY*241, LINE*80
      LOGICAL LET, LIMSCF
      LIMSCF=.FALSE.
      TYPE=' MNDO  '
      LET=(INDEX(KEYWRD,'LET').NE.0)
      IF(INDEX(KEYWRD,'MINDO') .NE. 0) TYPE='MINDO/3'
      IF(INDEX(KEYWRD,'AM1') .NE. 0)    TYPE='  AM1  '
      WRITE(6,10)
   10 FORMAT('1',20('*'),' TDHF POLARIZABILITIES ',
     1        20('*'),//)
      CALL GMETRY(GEO,COORD)
C
C  ORIENT THE MOLECULE WITH THE MOMENTS OF INERTIA.
C  THIS IS DONE TO ENSURE A UNIQUE, REPRODUCEABLE SET OF DIRECTIONS.
C  IF LET IS SPECIFIED, THE INPUT ORIENTATION WILL BE USED.
C
      IF (.NOT.LET) THEN
         MASS = 1
         CALL AXIS(COORD,NUMAT,A,B,C,SUMW,MASS,ROTVEC)
         WRITE(6,20)
   20    FORMAT (/' ROTATION MATRIX FOR ORIENTATION OF MOLECULE:'/)
         DO 40 I = 1,3
            WRITE(6,30) (ROTVEC(I,J),J=1,3)
   30       FORMAT(5X,3F12.6)
   40    CONTINUE
C
C  ROTATE ATOMS
C
         DO 70 I = 1,NUMAT
            DO 60 J = 1,3
               SUM = 0.0D00
               DO 50 K = 1,3
                  SUM = SUM + COORD(K,I)*ROTVEC(K,J)
   50          CONTINUE
               GEO(J,I) = SUM
   60       CONTINUE
   70    CONTINUE
         DO 90 I = 1,NUMAT
            DO 80 J = 1,3
               COORD(J,I) = GEO(J,I)
   80       CONTINUE
   90    CONTINUE
         WRITE(6,'(//10X,''CARTESIAN COORDINATES '',/)')
         WRITE(6,'(4X,''NO.'',7X,''ATOM'',9X,''X'',
     1  9X,''Y'',9X,''Z'',/)')
         L=0
         DO 100 I=1,NUMAT
            IF(NAT(I).EQ.99.OR.NAT(I).EQ.107) GOTO 100
            L=L+1
            WRITE(6,'(I6,8X,A2,4X,3F10.4)')
     1           L,ELEMNT(NAT(I)),(COORD(J,L),J=1,3)
  100    CONTINUE
C
C  IF POLYMER, ROTATE TVEC
C  (BEWARE:  THE POLYMER SECTIONS MAY NOT WORK YET)
C
         IF (IDTVEC.GT.0) THEN
            DO 130 I = 1,IDTVEC
               DO 120 J = 1,3
                  SUM = 0.0D00
                  DO 110 K = 1,3
                     SUM = SUM + TVEC(K,I)*ROTVEC(K,J)
  110             CONTINUE
                  TEMPV(J,I) = SUM
  120          CONTINUE
  130       CONTINUE
            DO 150 I = 1,3
               DO 140 J = 1,IDTVEC
                  TVEC(I,J) = TEMPV(I,J)
  140          CONTINUE
  150       CONTINUE
            WRITE(6,160)((TVEC(J,I),J=1,3),I=1,IDTVEC)
  160       FORMAT(/' NEW TRANSLATION VECTOR:'/,
     1           ' ',3(3F15.5))
         ENDIF
      ENDIF
C
      LAST=1
      NA(1)=99
C
C  SET UP THE VARIABLES IN XPARAM AND LOC, THESE ARE IN CARTESIAN
C  COORDINATES.
C
      NDEP=0
      NUMAT=0
      SUMX=0.D0
      SUMY=0.D0
      SUMZ=0.D0
      DO 180 I=1,NATOMS
         IF((LABELS(I).NE.99).AND.(LABELS(I).NE.107)) THEN
            NUMAT=NUMAT+1
            LABELS(NUMAT)=LABELS(I)
            SUMX=SUMX+COORD(1,NUMAT)
            SUMY=SUMY+COORD(2,NUMAT)
            SUMZ=SUMZ+COORD(3,NUMAT)
            DO 170 J=1,3
  170       GEO(J,NUMAT)=COORD(J,NUMAT)
         ENDIF
  180 CONTINUE
      SUMX=SUMX/NUMAT
      SUMY=SUMY/NUMAT
      SUMZ=SUMZ/NUMAT
      SUMMAX=0.D0
      ATPOL=0.D0
      DO 190 I=1,NUMAT
         IF (LABELS(I).NE.107) THEN
            ATPOL=ATPOL+POLVOL(LABELS(I))
         ENDIF
         GEO(1,I)=GEO(1,I)-SUMX
         IF(SUMMAX.LT.ABS(GEO(1,I))) SUMMAX=ABS(GEO(1,I))
         GEO(2,I)=GEO(2,I)-SUMY
         IF(SUMMAX.LT.ABS(GEO(2,I))) SUMMAX=ABS(GEO(2,I))
         GEO(3,I)=GEO(3,I)-SUMZ
         IF(SUMMAX.LT.ABS(GEO(3,I))) SUMMAX=ABS(GEO(3,I))
  190 CONTINUE
C
      NVAR=0
      NATOMS = NUMAT
      CALL COMPFG(GEO, .TRUE., HEAT0, .TRUE., GRAD, .FALSE.)
      WRITE(6,200) HEAT0
  200 FORMAT(//' ENERGY OF "REORIENTED" SYSTEM WITHOUT FIELD:',
     1        F20.10)
C...............................................................
C
C  VARIABLES USED FOR TIME-DEPENDENT CALCULATIONS
C
C    OMEGA .........  FREQUENCY OF LIGHT (ACTUALLY INPUT AS ENERGY
C                     IN EV'S.
C    IWFLA .........  TYPE OF ALPHA CALCULATION FOR STORING MATRICES
C                     0 = STATIC
C                     1 = OMEGA
C                     2 = 2*OMEGA
C                     3 = 3*OMEGA
C    IWFLB .........  TYPE OF BETA CALCULATION FOR STORING MATRICES
C                     0 = (0,0)
C                     1 = (W,W) (SHG)
C                     2 = (0,W) (EOPE)
C                     3 = (W,-W) (OR)
C
C  INPUT NUMBER OF FREQENCIES TO RUN
C
C     IBET = 0  NO BETA CALC
C            1  ITERATIVE BETA
C           -1  NONITER BETA (SHG)
C           -2  NONITER EOPE
C           -3  NONITER OR
C
C     IGAM = 0  NO GAMMA CALC
C            1 THIRD HARMONIC GENERATION INPUT N,0,1,1
C            2 DC-EFISHG INPUT N,0,1,2
C            3 IDRI N,0,1,3
C            4 OKE N,0,1,4
C            5 DC EFIOR (NOT AVAILABLE)
C
      READ(5,'(A)') LINE
      CALL NUCHAR(LINE,VALUE,NVALUE)
      IWFLB=VALUE(1)
      IBET=VALUE(2)
      IGAM=VALUE(3)
      ATOL=VALUE(4)
      MAXITU=VALUE(5)
      MAXITA=VALUE(6)
      BTOL=VALUE(7)
         DO 220 I=1,3
  220    DATAEV(I)=(I-1)*0.25D0
         NFREQ=3
  230 CONTINUE
      IF (IGAM.NE.0) THEN
         IBET = 1
      ENDIF
      WRITE(6,240) NFREQ,IWFLB,IBET,IGAM
  240 FORMAT(//'  NFREQ=',I3,'  IWFLB=',I3,'  IBET=',I3,'  IGAM=',I3)
C
C ATOL IS THE MAXIMUM TOLERANCE IN MAKEUF AND BTOL IS THAT IN BMAKUF
C MAXITU IS THE MAXIMUM ITERATION IN BETAF AND MAXITA IS THE MAXIMUM
C ITERATION IN ALPHAF
C
C#      READ(IR,*,END=99,ERR=99) ATOL,MAXITU,MAXITA,BTOL
      WRITE(6,250) ATOL,BTOL,MAXITU,MAXITA
  250 FORMAT('  ATOL=',D12.5,'  BTOL=',D12.5,'    MAXITU=',I5,
     1        '    MAXITA=',I5)
C
C SET UP DIRECT ACCESS FILE FOR T-D MATRICES
      CALL OPENDA(0)
C
C CALCULATE ALPHA AT STATIC VALUES
C
      IF ((IWFLB .EQ. 2).OR.(IGAM .EQ. 2).OR.(IGAM.EQ.4).OR.
     1   (IBET.LE.-2)) THEN
         IWFLA = 0
         OMEGA = 0.0D00
         CALL ALPHAF(IWFLA,ATOL,MAXITA,X1,X2,X3,X4,X5,X6,X7)
      ENDIF
      IF (IGAM.EQ.4) THEN
         IWFLB=0
         CALL BETAF(IWFLB,MAXITU,BTOL,X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,
     1              X11,X12,X13)
      ENDIF
C
C CALCULATE FREQUENCY DEPENDENT VALUES
C
      DO 280 I = 1, NFREQ
C
C  READ IN FREQ:  ACTUALLY READ IN AS ENERGY IN EV.
C
         OMEGA=DATAEV(I)
         OMEGAU = OMEGA/27.2113961D+00
         IF (OMEGA.LT.1.0D-8) THEN
            WAVLEN = 999999.99D0
C#           WRITE(6,401) OMEGA
            WRITE(6,260)
  260       FORMAT(//,' ',65(1H*),/,
     1            ' CALCULATION OF STATIC FIELD QUANTITIES',/,
     2            ' ',65('*'))
         ELSE
            WRITE(6,270) OMEGA,OMEGAU,1239.8424D0/OMEGA,
     1 8065.541D0*OMEGA
  270       FORMAT(//,' ',70(1H*),
     1          /' CALCULATION FOR A FREQUENCY OF ',F10.5,' EV  =',
     2           F14.5,' A.U. '/18X,'WAVELENGTH OF ',F10.2,' NM  =',
     3           F14.5,' CM(-1)',/,
     4           ' ',70('*'))
         ENDIF
C
C  CALCULATE ALPHA(W)
C
         IWFLA = 1
         CALL ALPHAF(IWFLA,ATOL,MAXITA,X1,X2,X3,X4,X5,X6,X7)
C
C  PERFORM NONITERATIVE BETA CALCULATIONS
C
C   OPTICAL RECTIFICATION
         IF (IBET.EQ.-3) THEN
            CALL NONOR(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12)
         ENDIF
C   ELECTROPTIC POCKELS EFFECT
         IF (IBET.EQ.-2) THEN
            CALL NONOPE(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12)
         ENDIF
C   SECOND HARMONIC GENERATION
         IF (IBET.EQ.-1) THEN
            IWFLA = 2
            OMEGA = OMEGA*2.0D00
            CALL ALPHAF(IWFLA,ATOL,MAXITA,X1,X2,X3,X4,X5,X6,X7)
            OMEGA = OMEGA/2.0D00
            CALL NONBET(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12)
         ENDIF
C
C  PERFORM ITERATIVE BETA (SHG AND STATIC)CALCULATIONS
C
C         IF ((IBET.GT.0) .AND.(IGAM .EQ. 0)) THEN
         IF ((IBET.EQ.1) .AND. (IWFLB .LE. 1) .AND. (IGAM .EQ. 0)) THEN
            CALL BETAF(IWFLB,MAXITU,BTOL,X1,X2,X3,X4,X5,X6,X7,X8,
     1                 X9,X10,X11,X12,X13)
C
C PERFORM ITERATIVE BETA (EOPE AND OR) CALCULATIONS
C
         ELSEIF ((IBET.EQ.1).AND.(IWFLB.GT.1).AND.(IGAM.EQ.0)) THEN
            CALL BEOPOR(IWFLB,MAXITU,BTOL,X1,X2,X3,X4,X5,X6,X7,X8,
     1                  X9,X10,X11,X12,X13)
         ENDIF
C.......................................................................
C CALCULATE GAMMA VALUES
C.......................................................................
         IF ((IBET.GT.0).AND.(IGAM .LE. 3).AND.(IGAM.NE.0)) THEN
            IWFLB=1
            CALL BETAF(IWFLB,MAXITU,BTOL,X1,X2,X3,X4,X5,X6,X7,X8,
     1                 X9,X10,X11,X12,X13)
         ENDIF
C THIRD HARMONIC GENRATION
         IF (IGAM.EQ.1) THEN
            IWFLA = 3
            OMEGA = OMEGA*3.0D00
            CALL ALPHAF(IWFLA,ATOL,MAXITA,X1,X2,X3,X4,X5,X6,X7)
            OMEGA = OMEGA/3.0D00
            CALL NGAMTG(IGAM,X1,X2,X3,X4,X5,X6,X7,X8,X9)
         ENDIF
C DC-EFISHG
         IF (IGAM.EQ.2) THEN
            IWFLA = 2
            OMEGA = 2.0D00*OMEGA
            CALL ALPHAF(IWFLA,ATOL,MAXITA,X1,X2,X3,X4,X5,X6,X7)
            OMEGA = OMEGA/2.0D00
            IWFLB = 2
            CALL BEOPOR(IWFLB,MAXITU,BTOL,X1,X2,X3,X4,X5,X6,X7,X8,
     1                  X9,X10,X11,X12,X13)
            CALL NGEFIS(IGAM,X1,X2,X3,X4,X5,X6,X7,X8,X9)
         ENDIF
C IDRI
         IF (IGAM.EQ.3) THEN
            IWFLB=3
            CALL BEOPOR(IWFLB,MAXITU,BTOL,X1,X2,X3,X4,X5,X6,X7,X8,
     1                  X9,X10,X11,X12,X13)
            CALL NGIDRI(IGAM,X1,X2,X3,X4,X5,X6,X7,X8,X9)
         ENDIF
C OKE
         IF (IGAM.EQ.4) THEN
            IWFLB=2
            CALL BEOPOR(IWFLB,MAXITU,BTOL,X1,X2,X3,X4,X5,X6,X7,X8,
     1                  X9,X10,X11,X12,X13)
            CALL NGOKE(IGAM,X1,X2,X3,X4,X5,X6,X7,X8,X9)
         ENDIF
  280 CONTINUE
C
      RETURN
  290 WRITE(6,'('' DATA FOR POLAR CALCULATION EITHER'',
     1'' MISSING OR FAULTY'')')
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TF(UA,GA,UB,GB,T,NORBS,NCLOSE,IWFLB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C  THIS SUBROUTINE CREATES THE NEW T MATRIX
C
      DIMENSION UA(NORBS,NORBS),UB(NORBS,NORBS),
     1          T(NORBS,NORBS),GA(NORBS,NORBS),
     2          GB(NORBS,NORBS)
C
C  ZERO MATRIX INITIALLY
C
      CALL ZEROM(T,NORBS)
C
C CALCULATE T (IJ)(W,W)= SUM(GA(IK)(W)*UB(KJ)(W)+
C GB(IK)(W)*UA(KJ)(W)-UA(IK)(W)GB(KJ)(W)-UB(IK)(W)GA(KJ)(W)
C
      DO 30 I = 1,NORBS
         DO 20 J = 1,NORBS
            SUM1=0.0D0
            SUM2=0.0D0
            DO 10 K = 1,NORBS
C CALCULATE FOR (W,W), (0,W) VALUES
C
               SUM1 = SUM1+GA(I,K)*UB(K,J)+GB(I,K)*UA(K,J)
     1         -UA(I,K)*GB(K,J)-UB(I,K)*GA(K,J)
               SUM2 = SUM2+GA(J,K)*UB(K,I)+GB(J,K)*UA(K,I)
     1         -UA(J,K)*GB(K,I)-UB(J,K)*GA(K,I)
   10       CONTINUE
            T(I,J) = SUM1
            T(J,I) = SUM2
   20    CONTINUE
   30 CONTINUE
C
      RETURN
      END
      SUBROUTINE TRANSF(F,G,C,NORB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C  THIS SUBROUTINE FORMS THE G MATRIX BY TRANSFORMING F WITH C
C
      DIMENSION C(NORB,NORB),F(NORB,NORB),G(NORB,NORB)
C
      DO 40 I = 1,NORB
         DO 30 J = 1,NORB
C            IJ = I*(I-1)/2 + J
            TERM2 = 0.0D00
            DO 20 K = 1,NORB
               TERM = 0.0D00
               DO 10 L = 1,NORB
C                  KI = MAX0(K,L)
C                  LI = MIN0(K,L)
C                  KL = KI*(KI-1)/2 + LI
                  TERM = TERM + F(K,L)*C(L,J)
   10          CONTINUE
               TERM2 = TERM2 + TERM*C(K,I)
   20       CONTINUE
            G(I,J) = TERM2
   30    CONTINUE
   40 CONTINUE
      RETURN
      END
      FUNCTION TRSUB(UL,X,UR,L1,LM,NDIM)
C THIS PROGRAM CALCULATES TRACES OF MATRICES
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION UL(NDIM,NDIM),X(NDIM,NDIM),UR(NDIM,NDIM)
C
      SUM = 0.0D00
      DO 30 I = 1,L1
         DO 20 K = 1,LM
            SUML = 0.0D00
            DO 10 L = 1,LM
               SUML = SUML + X(K,L)*UR(L,I)
   10       CONTINUE
            SUM = SUM + SUML*UL(I,K)
   20    CONTINUE
   30 CONTINUE
      TRSUB = 2.0D00*SUM
      RETURN
      END
      FUNCTION TRUDGU(UL,X,UR,L1,LM,NDIM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION UL(NDIM,NDIM),X(NDIM,NDIM),UR(NDIM,NDIM)
C
      SUM = 0.0D00
      DO 30 I = 1,L1
         DO 20 K = 1,LM
            SUML = 0.0D00
            DO 10 L = 1,LM
               SUML = SUML + X(K,L)*UR(L,I)
   10       CONTINUE
            SUM = SUM + SUML*UL(K,I)
   20    CONTINUE
   30 CONTINUE
      TRUDGU = 2.0D00*SUM
      RETURN
      END
      FUNCTION TRUGDU(UL,X,UR,L1,LM,NDIM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION UL(NDIM,NDIM),X(NDIM,NDIM),UR(NDIM,NDIM)
C
      SUM = 0.0D00
      DO 30 I = 1,L1
         DO 20 K = 1,LM
            SUML = 0.0D00
            DO 10 L = 1,LM
               SUML = SUML + X(L,K)*UR(L,I)
   10       CONTINUE
            SUM = SUM + SUML*UL(I,K)
   20    CONTINUE
   30 CONTINUE
      TRUGDU = 2.0D00*SUM
      RETURN
      END
      FUNCTION TRUGUD(UL,X,UR,L1,LM,NDIM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION UL(NDIM,NDIM),X(NDIM,NDIM),UR(NDIM,NDIM)
C
      SUM = 0.0D00
      DO 30 I = 1,L1
         DO 20 K = 1,LM
            SUML = 0.0D00
            DO 10 L = 1,LM
               SUML = SUML + X(K,L)*UR(I,L)
   10       CONTINUE
            SUM = SUM + SUML*UL(I,K)
   20    CONTINUE
   30 CONTINUE
      TRUGUD = 2.0D00*SUM
      RETURN
      END
      SUBROUTINE ZEROM(X,M)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C  ZEROM ZEROS THE MATRIX X
C
      DIMENSION X(M,M)
      DO 20 I = 1,M
         DO 10 J = 1,M
            X(I,J) = 0.0D00
   10    CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE ALPHAF (IWFLA,ATOL,MAXITA,U,F,G,UOLD,H1,D,DA)
C
C  SUBROUTINE FOR THE CALCULATION OF THE FREQUENCY DEPENDENT FIRST-ORDER
C  RESPONCE MATRICIES UA AND DENSITIES DA.
C  USED TO COMPUTE THE FREQUENCY DEPENDENT POLARIZABILITY AND FOR
C  SOLVING THE SECOND-ORDER PROBLEM.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL LAST
      INCLUDE 'SIZES'
C
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM),NORBS,NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /VECTOR/ C(MORB2),EIGS(MAXORB),CA(MORB2),DUMY(MAXORB)
      COMMON /WMATRX/ W(N2ELEC*2)
      COMMON /GEOM  / GEO(3,NUMATM), COORD(3,NUMATM)
      COMMON /OMVAL/ OMEGA
      COMMON /KEYWRD/ KEYWRD
      CHARACTER*241 KEYWRD
      COMMON /CHANEL/ IFILES(30)
      EQUIVALENCE(IW,IFILES(6))
C
      DIMENSION U(MAXORB,MAXORB), F(MAXORB,MAXORB), G(MAXORB,MAXORB)
      DIMENSION UOLD(MAXORB,MAXORB), H1(MAXORB,MAXORB),
     1          D(MAXORB,MAXORB),DA(MAXORB,MAXORB), ALLALP(3,3)
      CHARACTER*1 ALAB
      DIMENSION ALAB(3)
      SAVE ALAB
      DATA ALAB/'X','Y','Z'/
C
      NSQR = NORBS*NORBS
      ALPAVG = 0.0D00
C COMPUTE OFFSETS FOR U AND G MATRICES
      IPOSU = 1 + 6*IWFLA
      IPOSG = 4 + 6*IWFLA
      WRITE(6,10) OMEGA
   10 FORMAT (/,' +++++ ALPHA AT ',1F13.5,' EV.')
C
C  CHOOSE A  COMPONENT
C  X: ID=1   Y: ID=2   Z: ID=3
C
      DO 70 ID = 1,3
         CMPTIM = SECOND()
         LAST = .FALSE.
C
C  CALCULATE THE DIPOLE MATRIX.
C
         CALL HMUF (H1,ID,COORD,NFIRST,NLAST,NAT,NORBS,NUMAT)
         CALL COPYM (H1,F,NORBS)
C
C  INITIALIZE UOLD TO ZERO
C
         CALL ZEROM (UOLD,NORBS)
C.................................................................
C  LOOP STARTS HERE
C.................................................................
         ICOUNT = 0
         ALPOLD = 0.0D00
   20    CONTINUE
         ICOUNT = ICOUNT + 1
         IF (ICOUNT.GT.MAXITA) LAST = .TRUE.
C
C  CREATE G MATRIX.
C
         CALL TRANSF (F,G,C,NORBS)
C
C  FORM U MATRIX
C
         CALL MAKEUF (U,UOLD,G,EIGS,LAST,NORBS,NNORB,NCLOSE,DIFF,ATOL)
C
C  FORM NEW DENSITY MATRIX
C
         CALL DENSF (U,C,CA,D,DA,NORBS,NCLOSE)
C
C COMPUTE TEST ALPHA TO BE USED FOR A CONVERGENCE TEST
C
         ALPHAW = AVAL(H1,D,NORBS)
         DELA = DABS(ALPOLD-ALPHAW)
         ALPOLD = ALPHAW
C.      WRITE(6,1500) ALPHAW
C. 1500 FORMAT ('  TEST ALPHA = ',D12.5)
C
C  CREATE NEW FOCK MATRIX
C
         CALL ZEROM (F,NORBS)
         CALL FFREQ2 (F,D,W,NUMAT,NFIRST,NLAST,NORBS)
         CALL FFREQ1 (F,D,DA,DA,NORBS)
         CALL HPLUSF (F,H1,NORBS)
C..............................................................
         IF (.NOT.LAST) GO TO 20
         CMPTIM = SECOND() - CMPTIM
         WRITE(6,30) ICOUNT,CMPTIM,DIFF,DELA
   30    FORMAT (/' CONVERGED IN',I4,' ITERATIONS IN',F10.2,' SECONDS',
     1         /'           DENSITY CONVERG. TO ',1PD12.5,
     2         /'             ALPHA CONVERG. TO ',1PD12.5,/)
C
C COMPUTE ALPHA
C
         ALPHAW = AVAL(H1,D,NORBS)
         ALLALP(ID,ID)=ALPHAW
         WRITE(6,40) ALAB(ID),ALAB(ID),ALPHAW
   40    FORMAT ('      ALPHA(',A1,',',A1,') = ',1PD14.7)
         ALPAVG = ALPAVG + ALPHAW
C
C  WRITE OUT U AND G FOR FUTURE USE
C
         CALL DAWRIT (U,NSQR,IPOSU+ID)
         CALL DAWRIT (G,NSQR,IPOSG+ID)
C
C  COMPUTE OTHER COMPONENTS
C
         DO 60 IC = 1,3
            IF (IC.NE.ID) THEN
               CALL HMUF (H1,IC,COORD,NFIRST,NLAST,NAT,NORBS,NUMAT)
               ALPHAW = AVAL(H1,D,NORBS)
               ALLALP(IC,ID)=ALPHAW
               WRITE(6,50) ALAB(IC),ALAB(ID),ALPHAW
   50          FORMAT ('      ALPHA(',A1,',',A1,') = ',1PD14.7)
            ENDIF
   60    CONTINUE
   70 CONTINUE
      ALPAVG = ALPAVG/3.0D00
      WRITE(6,80) ALPAVG
   80 FORMAT (/,'  ISOTROPIC AVERAGE ALPHA = ',1F13.5,' A.U.')
C
      RETURN
      END
      FUNCTION AVAL (H,D,NORBS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C.................................................................
C  COMPUTE POLARIZABILITY AS TRACE OF H*D
C.................................................................
      DIMENSION H(NORBS,NORBS),D(NORBS,NORBS)
      SUM = 0.0D00
      DO 20 I = 1,NORBS
         DO 10 J = 1,NORBS
            SUM = SUM + H(I,J)*D(J,I)
   10    CONTINUE
   20 CONTINUE
      AVAL = -SUM
      RETURN
      END
      SUBROUTINE BDENSF (UA,UB,UAB,C,D,DA,NORBS,NCLOSE,IWFLB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /CHANEL/ IFILES(30)
      EQUIVALENCE(IW,IFILES(6))
C
C  THIS SUBROUTINE IS USED TO COMPUTE THE FIRST-ORDER DENSITY
C
      DIMENSION C(NORBS,NORBS),D(NORBS,NORBS),UA(NORBS,NORBS),
     1       UB(NORBS,NORBS),DA(NORBS,NORBS),UAB(NORBS,NORBS)
C
C FORM DENSITY MATRIX
C
C
      CALL ZEROM(D,NORBS)
C
C CALCULATE
C
      DO 100 I = 1,NORBS
         DO 90 J = 1,NORBS
            S1 = 0.0D00
            S2 = 0.0D00
            S3 = 0.0D00
            S4 = 0.0D00
            DO 20 K = 1,NORBS
               DO 10 L = 1,NCLOSE
                  S1 = S1+C(I,K)*UAB(K,L)*C(J,L)
                  S2 = S2+C(I,L)*UAB(L,K)*C(J,K)
   10          CONTINUE
   20       CONTINUE
C
            DO 50 K = 1,NCLOSE
               DO 40 L= NCLOSE+1,NORBS
                  DO 30 M = 1,NCLOSE
                     S3 = S3+C(I,K)*(UA(K,L)*UB(L,M)+UB(K,L)*UA(L,M))*
     1C(J,M)
   30             CONTINUE
   40          CONTINUE
   50       CONTINUE
C
            DO 80 K = NCLOSE+1,NORBS
               DO 70 L = 1,NCLOSE
                  DO 60 M = NCLOSE+1, NORBS
                     S4 = S4+C(I,K)*(UA(K,L)*UB(L,M)+UB(K,L)*UA(L,M))*
     1C(J,M)
   60             CONTINUE
   70          CONTINUE
   80       CONTINUE
            D(I,J) =  2.0D00*(S1-S2+S3-S4)
   90    CONTINUE
  100 CONTINUE
C      WRITE(6,*) 'INITIAL DENSITY MATRIX FINAL FORM'
C      CALL MATOUT(D,EIGS,NORBS,NORBS,NORBS)
C
C CREATE DA
C
      DO 120 I = 1,NORBS
         DO 110 J = 1,NORBS
            DA(I,J) = D(I,J)/2.0D00
  110    CONTINUE
  120 CONTINUE
C
      RETURN
      END
      SUBROUTINE BEOPOR(IWFLB,MAXITU,BTOL,UA,UB,F,GA,GB,T,H1,
     1                   D,DA,UAB,UOLD1,G,X)
C
C THIS SUBROUTINE CALCULATES ITERATIVE BETA VALUES FOR
C THE ELECTROOPTIC POCKELS EFFECT AND OPTICAL RECTIFICATION
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION MAXU
      LOGICAL LAST
      INCLUDE 'SIZES'
C
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM),NORBS,NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /VECTOR/ C(MORB2),EIGS(MAXORB),CA(MORB2),DUMY(MAXORB)
      COMMON /WMATRX/ W(N2ELEC*2)
      COMMON /GEOM /  GEO(3,NUMATM),  COORD(3,NUMATM)
      COMMON /OMVAL/ OMEGA
      COMMON /CHANEL/ IFILES(30)
      EQUIVALENCE(IW,IFILES(6))
C
      DIMENSION UA(MAXORB,MAXORB),UB(MAXORB,MAXORB),F(MAXORB,MAXORB),
     1          GA(MAXORB,MAXORB),GB(MAXORB,MAXORB),
     2          T(MAXORB,MAXORB), H1(MAXORB,MAXORB),
     3          D(MAXORB,MAXORB),DA(MAXORB,MAXORB),
     4          UAB(MAXORB,MAXORB),UOLD1(MAXORB,MAXORB),
     5          G(MAXORB,MAXORB),X(MAXORB,MAXORB)
      CHARACTER*1 ALAB
      DIMENSION ALAB(3),IDA(9),IDB(9)
      SAVE ALAB, IDA, IDB
      DATA ALAB/'X','Y','Z'/
      DATA IDA /1,1,1,2,2,2,3,3,3/
      DATA IDB /1,2,3,1,2,3,1,2,3/
      ONE=1.0D00
      BETOLD = 0.0D00
      MAXSQ = NORBS*NORBS
      IF (IWFLB .EQ. 2) THEN
         IPOSU = 73
      ELSE
         IPOSU = 109
      ENDIF
      IPOSG = IPOSU + 9
      IPOSE = IPOSG + 9
      IPOSUM = IPOSE + 9
      IF (IWFLB .EQ. 0) THEN
         WRITE(6,10) OMEGA
   10    FORMAT(/,' +++++ BETA (STATIC) AT ',1F15.5 ,' EV.'/)
      ELSEIF (IWFLB .EQ. 2) THEN
         WRITE(6,20) OMEGA
   20    FORMAT(/,' +++++ BETA',
     1          ' (ELECTROOPTIC POCKELS EFFECT) AT ',1F15.5 ,' EV.'/)
      ELSE
         WRITE(6,30) OMEGA
   30    FORMAT(/,' +++++ BETA',
     1          ' (OPTICAL RECTIFICATION) AT ',1F15.5 ,' EV.'/)
      ENDIF
C
C  LOOP OVER COMPONENTS
C
      BAVX = 0.0D+00
      BAVY = 0.0D+00
      BAVZ = 0.0D+00
      DO 90 ID = 1,9
         CMPTIM = SECOND()
         IA=IDA(ID)
         IB=IDB(ID)
         LAST = .FALSE.
C
C  CALCULATE THE DIPOLE MATRIX.
C
         CALL HMUF(H1,IA,COORD,NFIRST,NLAST,NAT,NORBS,NUMAT)
C
C  INITIALIZE ZERO ARRAYS
C
         CALL ZEROM(UOLD1,NORBS)
         CALL ZEROM(UAB,NORBS)
         CALL ZEROM(F,NORBS)
C
C  INPUT U AND GA FROM ALPHA CALCULATIONS
C
         IF ((IWFLB .EQ. 2) .OR. (IWFLB .EQ. 0)) THEN
C  UA CONTAINS UA(0)
            JPU = 1 + IA
            CALL DAREAD(UA,MAXSQ,JPU)
C  GA CONTAINS GA(0)
            JPG = 4 + IA
            CALL DAREAD(GA,MAXSQ,JPG)
         ELSE
C  UA CONTAINS UA(W)
            JPU = 7 + IA
            CALL DAREAD(UA,MAXSQ,JPU)
C  GA CONTAINS GA(W)
            JPG = 10 + IA
            CALL DAREAD(GA,MAXSQ,JPG)
         ENDIF
C
C READ VALUES FOR (W,-W) CALCULATION  :  OR
C
         IF (IWFLB .EQ. 3) THEN
C  UB CONTAINS UB(-W) = -UB+(W)
            JPU = 7 + IB
            CALL DAREAD(X,MAXSQ,JPU)
            CALL FHPATN(UB,X,NORBS,2,-ONE)
C  GB CONTAINS GB(-W) = GB+(W)
            JPG = 10 + IB
            CALL DAREAD(X,MAXSQ,JPG)
            CALL FHPATN(GB,X,NORBS,2,ONE)
C
C READ VALUES FOR (0,W) CALCULATION  :  OKE
C
         ELSEIF (IWFLB .EQ. 0) THEN
C  UB CONTAINS UB(0)
            JPU = 1 + IB
            CALL DAREAD(UB,MAXSQ,JPU)
C  GB CONTAINS GB(0)
            JPG = 4 + IB
            CALL DAREAD(GB,MAXSQ,JPG)
         ELSE
C  UB CONTAINS UB(W)
            JPU = 7 + IB
            CALL DAREAD(UB,MAXSQ,JPU)
C  GB CONTAINS GB(W)
            JPG = 10 + IB
            CALL DAREAD(GB,MAXSQ,JPG)
         ENDIF
C
C  CONSTRUCT T-MATRIX ONE TIME
C
         CALL TF(UA,GA,UB,GB,T,NORBS,NCLOSE,IWFLB)
C
C  CALCULATE INITIAL DENSITY AND BETA VALUE
C
         CALL BDENSF(UA,UB,UAB,C,D,DA,NORBS,NCLOSE,IWFLB)
         BETAW = AVAL(H1,D,NORBS)
         DELA = DABS(BETOLD-BETAW)
         BETOLD = BETAW
C
C INITIALIZE FOCK MATRIX
C
         CALL FFREQ2(F,D,W,NUMAT,NFIRST,NLAST,NORBS)
         CALL FFREQ1(F,D,DA,DA,NORBS)
         CALL ZEROM(DA,NORBS)
         CALL HPLUSF(F,DA,NORBS)
C.................................................................
C  LOOP STARTS HERE
C.................................................................
         ICOUNT = 0
   40    CONTINUE
         ICOUNT = ICOUNT + 1
         IF (ICOUNT .GE. MAXITU) LAST = .TRUE.
C
C  CREATE G MATRIX.
C
         CALL TRANSF(F,G,C,NORBS)
C
C  FORM U MATRIX
C
         CALL BMAKUF(UA,UB,UAB,T,UOLD1,G,EIGS,LAST,NORBS,
     1              NCLOSE,DIFF,IWFLB,MAXU,BTOL)
C
C  FORM NEW DENSITY MATRIX
C
         CALL BDENSF(UA,UB,UAB,C,D,DA,NORBS,NCLOSE,IWFLB)
C...
C COMPUTE TEST BETA
C
         BETAW = AVAL(H1,D,NORBS)
         DELA = DABS(BETOLD-BETAW)
         BETOLD = BETAW
C         IF (LAST.OR.(ICOUNT.GT.(MAXITU-5))) THEN
C             WRITE(6,1500) ICOUNT,DELA,MAXU,DIFF
C 1500        FORMAT(' ',I4,'  DELTA BETA = ', D12.5,
C     X       ' MAXU = ', D12.5, '  UDIFF = ', D12.5)
C         ENDIF
C
C  CREATE NEW FOCK MATRIX
C
         CALL ZEROM(F,NORBS)
         CALL FFREQ2(F,D,W,NUMAT,NFIRST,NLAST,NORBS)
         CALL FFREQ1(F,D,DA,DA,NORBS)
         CALL ZEROM(DA,NORBS)
         CALL HPLUSF(F,DA,NORBS)
C..............................................................
         IF (.NOT.LAST) GO TO 40
         CMPTIM = SECOND() - CMPTIM
         WRITE(6,50) ICOUNT,CMPTIM
   50    FORMAT(/' CONVERGED IN',I4,' ITERATIONS IN',F10.2,
     1            ' SECONDS')
         WRITE(6,60) MAXU,DIFF
   60    FORMAT(' MAXIMUM UAB ELEMENT =',1F15.5 ,
     1           ',  MAXIMUM DIFFERENCE =',1F15.5 ,/)
C
C  COMPUTE OTHER COMPONENTS
C
         DO 80 IC = 1,3
            CALL HMUF(H1,IC,COORD,NFIRST,NLAST,NAT,NORBS,NUMAT)
            BETAW = AVAL(H1,D,NORBS)
            WRITE(6,70) ALAB(IC),ALAB(IA),ALAB(IB),BETAW
   70       FORMAT('      BETA(',A1,',',A1,',',A1,') = ',1F15.5 )
C CALCULATES THE AVERAGE VALUE OF BETA
C
            IF ((ID .EQ. 1) .AND. (IC .EQ. 1)) THEN
               BAVX = BAVX + 3.0D0*BETAW
            ELSEIF (((ID.EQ.5).OR.(ID.EQ.9)).AND. (IC .EQ. 1)) THEN
               BAVX = BAVX + BETAW
            ELSEIF (((ID.EQ.2).OR.(ID .EQ. 4)) .AND. (IC .EQ. 2)) THEN
               BAVX = BAVX + BETAW
            ELSEIF (((ID.EQ.3).OR.(ID .EQ. 7)) .AND. (IC .EQ. 3)) THEN
               BAVX = BAVX + BETAW
            ENDIF
C CALCULATES AVERAGE BETA IN Y-DIRECTION
C
            IF ((ID .EQ. 5) .AND. (IC .EQ. 2)) THEN
               BAVY = BAVY + 3.0D0*BETAW
            ELSEIF (((ID.EQ.2).OR.(ID .EQ. 4)) .AND. (IC .EQ. 1)) THEN
               BAVY = BAVY + BETAW
            ELSEIF (((ID.EQ.1).OR.(ID .EQ. 9)) .AND. (IC .EQ. 2)) THEN
               BAVY = BAVY + BETAW
            ELSEIF (((ID.EQ.6).OR.(ID .EQ. 8)) .AND. (IC .EQ. 3)) THEN
               BAVY = BAVY + BETAW
            ENDIF
C CALCULATES AVERAGE BETA IN THE Z-DIRECTION
C
            IF ((ID .EQ. 9) .AND. (IC .EQ. 3)) THEN
               BAVZ = BAVZ + 3.0D0 * BETAW
            ELSEIF (((ID.EQ.3).OR.(ID .EQ. 7)) .AND. (IC .EQ. 1)) THEN
               BAVZ = BAVZ + BETAW
            ELSEIF (((ID.EQ.6).OR.(ID .EQ. 8)) .AND. (IC .EQ. 2)) THEN
               BAVZ = BAVZ + BETAW
            ELSEIF (((ID.EQ.1).OR.(ID .EQ. 5)) .AND. (IC .EQ. 3)) THEN
               BAVZ = BAVZ + BETAW
            ENDIF
   80    CONTINUE
C
C CALL SUBROUTINE TO CALCULATE EPSILON AND UMINUS OMEGA,OMEGA
C  EPSILON IN H1 AND UMINUS IN DA
         CALL EPSAB(H1,EIGS,G,GA,GB,UA,UB,UAB,DA,NORBS,NCLOSE,IWFLB)
         CALL DAWRIT(UAB,MAXSQ,IPOSU+ID)
         CALL DAWRIT(G,MAXSQ,IPOSG+ID)
         CALL DAWRIT(H1,MAXSQ,IPOSE+ID)
         CALL DAWRIT(DA,MAXSQ,IPOSUM+ID)
   90 CONTINUE
      BAVX = BAVX/5.0D+00
      BAVY = BAVY/5.0D+00
      BAVZ = BAVZ/5.0D+00
      BVEC = (BAVX*BAVX+BAVY*BAVY+BAVZ*BAVZ)**0.5D+00
C
      WRITE(6,100) OMEGA,BAVX
  100 FORMAT(//,' AVERAGE BETAX VALUE AT ',F10.5,' EV = ',
     1   1F15.5 )
      WRITE(6,110) OMEGA,BAVY
  110 FORMAT(' AVERAGE BETAY VALUE AT ',F10.5,' EV = ',
     1   1F15.5 )
      WRITE(6,120) OMEGA,BAVZ
  120 FORMAT(' AVERAGE BETAZ VALUE AT ',F10.5,' EV = ',
     1   1F15.5 )
      WRITE(6,130) OMEGA,BVEC
  130 FORMAT(//,'  AVERAGE BETA VALUE AT ',F10.5,' EV = ',
     1   1F15.5 ,//)
C
      RETURN
      END
      SUBROUTINE BETAF(IWFLB,MAXITU,BTOL,UA,UB,F,GA,GB,T,H1,D,DA,
     1                  UAB,UOLD1,G,X)
C
C THIS SUBROUTINE CALCULATES ITERATIVE BETA VALUES FOR SECOND HARMONIC
C GENERATION.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION MAXU
      LOGICAL LAST
      INCLUDE 'SIZES'
C
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM),NORBS,NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /VECTOR/ C(MORB2),EIGS(MAXORB),CA(MORB2),DUMY(MAXORB)
      COMMON /WMATRX/ W(N2ELEC*2)
      COMMON /GEOM / GEO(3,NUMATM), COORD(3,NUMATM)
      COMMON /OMVAL/ OMEGA
      COMMON /KEYWRD/ KEYWRD
      COMMON /CHANEL/ IFILES(30)
      EQUIVALENCE(IW,IFILES(6))
      CHARACTER*241 KEYWRD
C
      DIMENSION UA(MAXORB,MAXORB),UB(MAXORB,MAXORB),F(MAXORB,MAXORB),
     1          GA(MAXORB,MAXORB),GB(MAXORB,MAXORB),
     2          T(MAXORB,MAXORB), H1(MAXORB,MAXORB),
     3          D(MAXORB,MAXORB),DA(MAXORB,MAXORB),
     4          UAB(MAXORB,MAXORB),UOLD1(MAXORB,MAXORB),
     5          G(MAXORB,MAXORB),X(MAXORB,MAXORB)
      DIMENSION ALLBET(3,3,3)
      CHARACTER*1 ALAB
      DIMENSION ALAB(3),IDA(6),IDB(6)
      SAVE ALAB, IDA, IDB
      DATA ALAB/'X','Y','Z'/
      DATA IDA /1,1,1,2,2,3/
      DATA IDB /1,2,3,2,3,3/
C
      ONE=1.0D00
      BETOLD = 0.0D00
      MAXSQ = NORBS*NORBS
      IPOSU = 25 + 24*IWFLB
      IPOSG = IPOSU + 6
      IPOSE = IPOSG + 6
      IPOSUM = IPOSE + 6
C
      IF (IWFLB .EQ. 0) THEN
         WRITE(6,10) OMEGA
   10    FORMAT(/,' +++++ BETA (STATIC) AT ',1F15.5 ,' EV.'/)
      ELSE
         WRITE(6,20) OMEGA
   20    FORMAT(/,' +++++ BETA',
     1           ' (SECOND HARMONIC GENERATION) AT ',1F13.5,' EV.'/)
      ENDIF
C
C  CHOOSE A  COMPONENT
C  X: ID=1   Y: ID=2   Z: ID=3
C
      BAVX = 0.0D+00
      BAVY = 0.0D+00
      BAVZ = 0.0D+00
      DO 80 ID = 1,6
         CMPTIM = SECOND()
         IA=IDA(ID)
         IB=IDB(ID)
         LAST = .FALSE.
C
C  CALCULATE THE DIPOLE MATRIX.
C
         CALL HMUF(H1,IA,COORD,NFIRST,NLAST,NAT,NORBS,NUMAT)
C
C  INITIALIZE ZERO ARRAYS
C
         CALL ZEROM(UOLD1,NORBS)
         CALL ZEROM(UAB,NORBS)
         CALL ZEROM(F,NORBS)
C
C  INPUT U AND GA FROM ALPHA CALCULATIONS
C
         IF ((IWFLB .EQ. 2) .OR. (IWFLB .EQ. 0)) THEN
            JPU = 1 + IA
            CALL DAREAD(UA,MAXSQ,JPU)
            JPG = 4 + IA
            CALL DAREAD(GA,MAXSQ,JPG)
         ELSE
            JPU = 7 + IA
            CALL DAREAD(UA,MAXSQ,JPU)
            JPG = 10 + IA
            CALL DAREAD(GA,MAXSQ,JPG)
         ENDIF
C READ VALUES FOR (W,-W)
         IF (IWFLB .EQ. 3) THEN
            JPU = 7 + IB
            CALL DAREAD(X,MAXSQ,JPU)
            CALL FHPATN(UB,X,NORBS,2,-ONE)
            JPG = 10 + IB
            CALL DAREAD(X,MAXSQ,JPG)
            CALL FHPATN(GB,X,NORBS,2,ONE)
C READ VALUES FOR OKE
C
         ELSEIF (IWFLB .EQ. 0) THEN
            JPU = 1 + IB
            CALL DAREAD(UB,MAXSQ,JPU)
            JPG = 4 + IB
            CALL DAREAD(GB,MAXSQ,JPG)
         ELSE
            JPU = 7 + IB
            CALL DAREAD(UB,MAXSQ,JPU)
            JPG = 10 + IB
            CALL DAREAD(GB,MAXSQ,JPG)
         ENDIF
C
C  CONSTRUCT T-MATRIX ONE TIME
C
         CALL TF(UA,GA,UB,GB,T,NORBS,NCLOSE,IWFLB)
C
C  CALCULATE INITIAL DENSITY AND BETA VALUE
C
         CALL BDENSF(UA,UB,UAB,C,D,DA,NORBS,NCLOSE,IWFLB)
         BETAW = AVAL(H1,D,NORBS)
         DELA = DABS(BETOLD-BETAW)
         BETOLD = BETAW
C
C INITIALIZE FOCK MATRIX
C
         CALL FFREQ2(F,D,W,NUMAT,NFIRST,NLAST,NORBS)
         CALL FFREQ1(F,D,DA,DA,NORBS)
         CALL ZEROM(DA,NORBS)
         CALL HPLUSF(F,DA,NORBS)
C.................................................................
C  LOOP STARTS HERE
C.................................................................
         ICOUNT = 0
   30    CONTINUE
         ICOUNT = ICOUNT + 1
         IF (ICOUNT .GE. MAXITU) LAST = .TRUE.
C
C  CREATE G MATRIX.
C
         CALL TRANSF(F,G,C,NORBS)
C
C  FORM U MATRIX
C
         CALL BMAKUF(UA,UB,UAB,T,UOLD1,G,EIGS,LAST,NORBS,
     1              NCLOSE,DIFF,IWFLB,MAXU,BTOL)
C
C  FORM NEW DENSITY MATRIX
C
         CALL BDENSF(UA,UB,UAB,C,D,DA,NORBS,NCLOSE,IWFLB)
C...
C COMPUTE TEST BETA
C
         BETAW = AVAL(H1,D,NORBS)
         DELA = DABS(BETOLD-BETAW)
         BETOLD = BETAW
C         IF (LAST.OR.(ICOUNT.GT.(MAXITU-5))) THEN
C             WRITE(6,1500) ICOUNT,DELA,MAXU,DIFF
C 1500        FORMAT(' ',I4,'  DELTA BETA = ', D12.5,
C     X       ' MAXU = ', D12.5, '  UDIFF = ', D12.5)
C         ENDIF
C
C  CREATE NEW FOCK MATRIX
C
         CALL ZEROM(F,NORBS)
         CALL FFREQ2(F,D,W,NUMAT,NFIRST,NLAST,NORBS)
         CALL FFREQ1(F,D,DA,DA,NORBS)
         CALL ZEROM(DA,NORBS)
         CALL HPLUSF(F,DA,NORBS)
C..............................................................
         IF (.NOT.LAST) GO TO 30
         CMPTIM = SECOND() - CMPTIM
         WRITE(6,40) ICOUNT,CMPTIM
   40    FORMAT(/' CONVERGED IN',I4,' ITERATIONS IN',F10.2,
     1            ' SECONDS')
         WRITE(6,50) MAXU,DIFF
   50    FORMAT(' MAXIMUM UAB ELEMENT =',1F15.5 ,
     1           ',  MAXIMUM DIFFERENCE =',1F15.5 ,/)
C
C COMPUTE BETA
C
C        CALL HMUF(H1,ID,COORD,NFIRST,NLAST,NAT,NORBS,NUMAT)
C        BETAW = AVAL(H1,D,NORBS)
C        WRITE(6,2000) ALAB(ID),ALAB(ID),ALAB(ID),BETAW
C2000    FORMAT('BETA(',A1,',',A1,','A1,') = ',D12.5)
C
C  COMPUTE OTHER COMPONENTS
C
         DO 70 IC = 1,3
            CALL HMUF(H1,IC,COORD,NFIRST,NLAST,NAT,NORBS,NUMAT)
            BETAW = AVAL(H1,D,NORBS)
            ALLBET(IC,IA,IB)=BETAW
            WRITE(6,60) ALAB(IC),ALAB(IA),ALAB(IB),BETAW
   60       FORMAT('      BETA(',A1,',',A1,',',A1,') = ',1F15.5 )
C
C CALCULATE AVERAGE BETA IN THE X-DIRECTION
C
            IF ((ID .EQ. 1) .AND. (IC .EQ. 1)) THEN
               BAVX = BAVX + 3.0D0 * BETAW
            ELSEIF ((ID .EQ. 2) .AND. (IC .EQ. 2)) THEN
               BAVX = BAVX + 2.0D0 * BETAW
            ELSEIF ((ID .EQ. 3) .AND. (IC .EQ. 3)) THEN
               BAVX = BAVX + 2.0D0 * BETAW
            ELSEIF (((ID.EQ.4).OR.(ID .EQ. 6)) .AND. (IC .EQ. 1)) THEN
               BAVX = BAVX + BETAW
            ENDIF
C CALCULATES AVERAGE BETA IN THE Y-DIRECTION
            IF ((ID .EQ. 4) .AND. (IC .EQ. 2)) THEN
               BAVY = BAVY + 3.0D0 * BETAW
            ELSEIF ((ID .EQ. 2) .AND. (IC .EQ. 1)) THEN
               BAVY = BAVY + 2.0D0 * BETAW
            ELSEIF ((ID .EQ. 5) .AND. (IC .EQ. 3)) THEN
               BAVY = BAVY + 2.0D0 * BETAW
            ELSEIF (((ID.EQ.1).OR.(ID .EQ. 6)) .AND. (IC .EQ. 2)) THEN
               BAVY = BAVY + BETAW
            ENDIF
C CALCULATES AVERAGE BETA IN THE Z-DIRECTION
            IF ((ID .EQ. 6) .AND. (IC .EQ. 3)) THEN
               BAVZ = BAVZ + 3.0D0 * BETAW
            ELSEIF ((ID .EQ. 3) .AND. (IC .EQ. 1)) THEN
               BAVZ = BAVZ + 2.0D0 * BETAW
            ELSEIF ((ID .EQ. 5) .AND. (IC .EQ. 2)) THEN
               BAVZ = BAVZ + 2.0D0 * BETAW
            ELSEIF (((ID.EQ.4) .OR. (ID.EQ.1)) .AND. (IC .EQ. 3)) THEN
               BAVZ = BAVZ + BETAW
            ENDIF
   70    CONTINUE
C
C
C CALL SUBROUTINE TO CALCULATE EPSILON AND UMINUS OMEGA,OMEGA
C  EPSILON IN H1 AND UMINUS IN DA
         CALL EPSAB(H1,EIGS,G,GA,GB,UA,UB,UAB,DA,NORBS,NCLOSE,IWFLB)
         CALL DAWRIT(UAB,MAXSQ,IPOSU+ID)
         CALL DAWRIT(G,MAXSQ,IPOSG+ID)
         CALL DAWRIT(H1,MAXSQ,IPOSE+ID)
         CALL DAWRIT(DA,MAXSQ,IPOSUM+ID)
   80 CONTINUE
C
      BAVX = BAVX/5.0D+00
      BAVY = BAVY/5.0D+00
      BAVZ = BAVZ/5.0D+00
C CALCULATES AVERAGE BETA
      BVEC = (BAVX*BAVX+BAVY*BAVY+BAVZ*BAVZ)**0.5D+00
C
      WRITE(6,90) OMEGA,BAVX
   90 FORMAT(//,' AVERAGE BETAX(SHG) VALUE AT',F10.5, ' EV = ',
     1   1F11.5 )
      WRITE(6,100) OMEGA,BAVY
  100 FORMAT(' AVERAGE BETAY(SHG) VALUE AT',F10.5, ' EV = ',
     1   1F11.5 )
      WRITE(6,110) OMEGA,BAVZ
  110 FORMAT(' AVERAGE BETAZ(SHG) VALUE AT',F10.5, ' EV = ',
     1   1F11.5)
C
      WRITE(6,120) OMEGA,BVEC
  120 FORMAT(//,' AVERAGE BETA (SHG) VALUE AT',F10.5, ' EV = ',
     1   1F11.5 ,//)
      RETURN
      END
      SUBROUTINE BETAL1(U0A,G0A,U1B,G1B,U1C,G1C,NCLOSE,NORBS,TERM)
C
C THIS SUBROUTINE CALCULATES THE TRACE OF UGU MATRICES
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U0A(NORBS,NORBS),U1B(NORBS,NORBS),U1C(NORBS,NORBS),
     1          G0A(NORBS,NORBS),G1B(NORBS,NORBS),G1C(NORBS,NORBS)
      T1A = TRUGUD(U0A,G1B,U1C,NCLOSE,NORBS,NORBS)
      T2A = TRUDGU(U1C,G1B,U0A,NCLOSE,NORBS,NORBS)
      T3A = TRUGDU(U1B,G1C,U0A,NCLOSE,NORBS,NORBS)
      T4A = TRUGDU(U0A,G1C,U1B,NCLOSE,NORBS,NORBS)
      T5A = TRUDGU(U1C,G0A,U1B,NCLOSE,NORBS,NORBS)
      T6A = TRUGUD(U1B,G0A,U1C,NCLOSE,NORBS,NORBS)
      T1B = TRUGUD(U0A,G1B,U1C,NORBS,NCLOSE,NORBS)
      T2B = TRUDGU(U1C,G1B,U0A,NORBS,NCLOSE,NORBS)
      T3B = TRUGDU(U1B,G1C,U0A,NORBS,NCLOSE,NORBS)
      T4B = TRUGDU(U0A,G1C,U1B,NORBS,NCLOSE,NORBS)
      T5B = TRUDGU(U1C,G0A,U1B,NORBS,NCLOSE,NORBS)
      T6B = TRUGUD(U1B,G0A,U1C,NORBS,NCLOSE,NORBS)
      TERM = T1B-T1A+T2B-T2A+T3A-T3B+T4A-T4B+T5B-T5A+T6B-T6A
      RETURN
      END
      SUBROUTINE BETALL(U2A,G2A,U1B,G1B,U1C,G1C,NCLOSE,NORBS,TERM)
C
C THIS SUBROUTINE CALCULATES TRACE OF UGU MATRICES
C WHEN A,B,C DIRECTIONS ARE DIFFERENT
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U2A(NORBS,NORBS),U1B(NORBS,NORBS),U1C(NORBS,NORBS),
     1          G2A(NORBS,NORBS),G1B(NORBS,NORBS),G1C(NORBS,NORBS)
      T1A = TRUDGU(U2A,G1B,U1C,NCLOSE,NORBS,NORBS)
      T2A = TRUGUD(U1C,G1B,U2A,NCLOSE,NORBS,NORBS)
      T3A = TRUGUD(U1B,G1C,U2A,NCLOSE,NORBS,NORBS)
      T4A = TRUDGU(U2A,G1C,U1B,NCLOSE,NORBS,NORBS)
      T5A = TRUGDU(U1C,G2A,U1B,NCLOSE,NORBS,NORBS)
      T6A = TRUGDU(U1B,G2A,U1C,NCLOSE,NORBS,NORBS)
      T1B = TRUDGU(U2A,G1B,U1C,NORBS,NCLOSE,NORBS)
      T2B = TRUGUD(U1C,G1B,U2A,NORBS,NCLOSE,NORBS)
      T3B = TRUGUD(U1B,G1C,U2A,NORBS,NCLOSE,NORBS)
      T4B = TRUDGU(U2A,G1C,U1B,NORBS,NCLOSE,NORBS)
      T5B = TRUGDU(U1C,G2A,U1B,NORBS,NCLOSE,NORBS)
      T6B = TRUGDU(U1B,G2A,U1C,NORBS,NCLOSE,NORBS)
      TERM = T1B-T1A+T2B-T2A+T3B-T3A+T4B-T4A+T5A-T5B+T6A-T6B
      RETURN
      END
      SUBROUTINE BETCOM(U1,G1,U2,G2,NCLOSE,NORBS,TERM)
C
C THIS SUBROUTINE CALCULATES TRACE OF UGU MATRICES
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U1(NORBS,NORBS),U2(NORBS,NORBS),
     1          G1(NORBS,NORBS),G2(NORBS,NORBS)
      T1A = TRUDGU(U2,G1,U1,NCLOSE,NORBS,NORBS)
      T2A = TRUGUD(U1,G1,U2,NCLOSE,NORBS,NORBS)
      T3A = TRUGDU(U1,G2,U1,NCLOSE,NORBS,NORBS)
      T1B = TRUDGU(U2,G1,U1,NORBS,NCLOSE,NORBS)
      T2B = TRUGUD(U1,G1,U2,NORBS,NCLOSE,NORBS)
      T3B = TRUGDU(U1,G2,U1,NORBS,NCLOSE,NORBS)
      TERM = 2.0D0*(T1B-T1A+T2B-T2A+T3A-T3B)
      RETURN
      END
      SUBROUTINE BMAKUF(UA,UB,UAB,T,UOLD1,GAB,EIGS,LAST,NORBS,
     1                    NCLOSE,DIFF,IWFLB,MAXU,BTOL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION MAXU
      LOGICAL LAST
C
C  THIS SUBROUTINE CREATES THE NEW TRANSFORMATION MATRIX U
C  AND THEN CHECKS FOR CONVERGENCE
C
      DIMENSION UA(NORBS,NORBS),UOLD1(NORBS,NORBS),T(NORBS,NORBS),
     1          EIGS(NORBS),GAB(NORBS,NORBS),UAB(NORBS,NORBS),
     2          UB(NORBS,NORBS)
      COMMON /OMVAL/ OMEGA
C
C  ZERO MATRIX INITIALLY
C      CALL ZEROM(UAB,NORBS)
C
C  CREATE DIAGONAL BLOCKS (OCC,OCC) AND (UNOCC,UNOCC)
C
      DO 30 I = 1,NORBS
         DO 20 J = 1,I
            SUM = 0.0D00
            IF (I.LE.NCLOSE) THEN
               KLL=NCLOSE+1
               KUL=NORBS
            ELSEIF ((I.GT.NCLOSE).AND.(J.GT.NCLOSE)) THEN
               KLL=1
               KUL=NCLOSE
            ENDIF
            DO 10 K =KLL,KUL
               SUM = SUM+UA(I,K)*UB(K,J)+UB(I,K)*UA(K,J)
   10       CONTINUE
            UAB(I,J) = SUM*0.5D00
            UAB(J,I) = SUM*0.5D00
   20    CONTINUE
   30 CONTINUE
C
C  CREATE OFF-DIAGONAL BLOCKS
C
      DO 80 K = NCLOSE+1,NORBS
         DO 70 L = 1,NCLOSE
            GO TO (40,50,60), IWFLB
C CALCULATE FOR (W,W) VALUES
C
   40       UAB(K,L) = 27.2113961D0*((GAB(K,L)+T(K,L))/
     1                 ((EIGS(L)-EIGS(K))-2.0D00*OMEGA))
            UAB(L,K) = 27.2113961D0*((GAB(L,K)+T(L,K))/
     1                 ((EIGS(K)-EIGS(L))-2.0D00*OMEGA))
            GO TO 70
C CALCULATE FOR (0,W) VALUES
C
   50       UAB(K,L) = 27.2113961D0*((GAB(K,L)+T(K,L))/
     1                 ((EIGS(L)-EIGS(K))-OMEGA))
            UAB(L,K) = 27.2113961D0*((GAB(L,K)+T(L,K))/
     1                 ((EIGS(K)-EIGS(L))-OMEGA))
            GO TO 70
C CALCULATE FOR (W,-W) VALUES
C
   60       UAB(K,L) = 27.2113961D0*((GAB(K,L)+T(K,L))/
     1                (EIGS(L)-EIGS(K)))
            UAB(L,K) = 27.2113961D0*((GAB(L,K)+T(L,K))/
     1                 (EIGS(K)-EIGS(L)))
   70    CONTINUE
   80 CONTINUE
C
C  CHECK FOR CONVERGENCE
C
      DIFF = 0.0D00
      MAXU = -1000.00D0
      DO 100 I = 1,NORBS
         DO 90 J = 1,NORBS
            UDIF = UAB(I,J)-UOLD1(I,J)
            IF (DIFF.LT.ABS(UDIF)) DIFF = ABS(UDIF)
            IF (MAXU.LT.UAB(I,J)) MAXU = UAB(I,J)
   90    CONTINUE
  100 CONTINUE
      IF (DIFF.LT.BTOL) THEN
         LAST = .TRUE.
      ENDIF
C
      DO 120 I = 1,NORBS
         DO 110 J = 1,NORBS
            UOLD1(I,J) = UAB(I,J)
  110    CONTINUE
  120 CONTINUE
C
      RETURN
      END
      SUBROUTINE COPYM(H,F,M)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C  COPYM COPIES MATRIX H INTO F
C
      DIMENSION F(M,M),H(M,M)
      DO 20 I = 1,M
         DO 10 J = 1,M
            F(I,J) = H(I,J)
   10    CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE DAREA1(V,LEN,IDAF,NS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION V(LEN)
C
C       READ A PHYSICAL RECORD FROM THE DAF
C
      READ(UNIT=IDAF, REC=NS) V
      RETURN
      END
      SUBROUTINE DAREAD(V,LEN,NREC)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION V(LEN)
C
      COMMON /IODAF/ IDAF,IRECLN,IRECST,IFILEN(145),IODA(145)
      COMMON /CHANEL/ IFILES(30)
      EQUIVALENCE(IW,IFILES(6))
C
C         READ A LOGICAL RECORD FROM THE DAF DICTIONARY FILE
C         A LOGICAL RECORD MAY SPAN SEVERAL PHYSICAL RECORDS.
C
      N = IODA(NREC)
      IF(N.EQ.-1) GO TO 20
      IS = -IRECLN + 1
      NS = N
      LENT = LEN
   10 CONTINUE
      IS = IS + IRECLN
      IF = IS + LENT - 1
      IF ((IF-IS+1) .GT. IRECLN) IF = IS + IRECLN - 1
      NSP = NS
      LENW = IF - IS + 1
      CALL DAREA1(V(IS),LENW,IDAF,NSP)
      LENT = LENT - IRECLN
      NS = NS + 1
      N = NS
      IF (LENT .GE. 1) GO TO 10
      RETURN
C
   20 CONTINUE
      WRITE(6,30) NREC,LEN
      STOP
C
   30 FORMAT(1X,'*** ERROR ***, ATTEMPT TO READ A DAF RECORD',
     1         ' THAT WAS NEVER WRITTEN. NREC,LEN=',I5,I10)
      END
      SUBROUTINE DAWRIT(V,LEN,NREC)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL NEWREC
C
      DIMENSION V(LEN)
C
      COMMON /IODAF/ IDAF,IRECLN,IRECST,IFILEN(145),IODA(145)
      COMMON /CHANEL/ IFILES(30)
      EQUIVALENCE(IW,IFILES(6))
C
C         WRITE A LOGICAL RECORD ON THE DAF DICTIONARY FILE
C         A LOGICAL RECORD MAY SPAN SEVERAL PHYSICAL RECORDS
C
      N = IODA(NREC)
      IF (N .GT. 0 .AND. LEN .NE. IFILEN(NREC)) GO TO 30
      NEWREC = .FALSE.
      IF (N .GT. 0) GO TO 10
      IODA(NREC) = IRECST
      IFILEN(NREC) = LEN
      NEWREC = .TRUE.
      IRECST = IRECST + (LEN-1)/IRECLN + 1
      N = IODA(NREC)
   10 CONTINUE
      IST = -IRECLN + 1
      NS = N
      LENT = LEN
   20 CONTINUE
      IST = IST + IRECLN
      IF = IST + LENT - 1
      IF ((IF-IST+1) .GT. IRECLN) IF = IST+IRECLN-1
      NSP = NS
      LENW = IF - IST + 1
      CALL DAWRT1(V(IST),LENW,IDAF,NSP)
      LENT = LENT - IRECLN
      NS = NS + 1
      N = NS
      IF (LENT .GE. 1) GO TO 20
      IF (NEWREC) WRITE(UNIT=IDAF, REC=1) IRECST,IODA,IFILEN,IS,IPK
C
C     THE CRAY/CTSS SYSTEM HAS A BUG IN IT!  THIS CALL TO DAWRIT
C     DOES NOT ALWAYS SUCCEED IN TRANFERRING DATA TO THE DISK,
C     LEAVING THE DATA ONLY IN THE BUFFER.  SUBSEQUENT CALLS
C     TO DAREAD FOR OTHER LOGICAL RECORDS WILL DESTROY THE
C     BUFFER RESIDENT DATA, AND THE DATA WILL BE LOST FOREVER.
C     THE FOLLOWING CALL QUARANTEES THE BUFFER IS FLUSHED
C     TO DISK.  IT SHOULD BE REMOVED IF THIS BUG IS EVER FIXED.
C
*CTS  CALL EMPTY(IDAF)
      RETURN
C
   30 CONTINUE
      WRITE(6,40) NREC,LEN,IFILEN(NREC)
      STOP
C
   40 FORMAT(1X,'DAWRIT HAS REQUESTED A RECORD WITH LENGTH',
     1       1X,'DIFFERENT THAN BEFORE - ABORT FORCED.'/
     2       1X,'DAF RECORD ',I5,' NEW LENGTH =',I5,' OLD LENGTH =',I5)
      END
C*MODULE IOLIB   *DECK DAWRT1
      SUBROUTINE DAWRT1(V,LEN,IDAF,NS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION V(LEN)
C
C     ----- WRITE A PHYSICAL RECORD ON THE DAF -----
C
      WRITE(UNIT=IDAF, REC=NS) V
      RETURN
      END
      SUBROUTINE DENSF(U,C,CA,D,DA,NORBS,NCLOSE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C  THIS SUBROUTINE IS USED TO COMPUTE THE FIRST-ORDER DENSITY
C  FROM CA = C*U
C
      DIMENSION C(NORBS,NORBS),CA(NORBS,NORBS),D(NORBS,NORBS)
      DIMENSION U(NORBS,NORBS),DA(NORBS,NORBS)
C
C  FORM DENSITY MATRIX      CA*N*C+ + C*N*CA+
C
      DO 40 I = 1,NORBS
         DO 30 J = 1,NORBS
            SUM = 0.0D00
            DO 20 K = 1,NORBS
               SK1 = 0.0D00
               SK2 = 0.0D00
               DO 10 L = 1,NCLOSE
                  SK1 = SK1 + U(K,L)*C(J,L)
                  SK2 = SK2 + C(I,L)*U(L,K)
   10          CONTINUE
               SUM = SUM + C(I,K)*SK1 - SK2*C(J,K)
   20       CONTINUE
            D(I,J) = 2.0D00*SUM
            DA(I,J) = SUM
   30    CONTINUE
   40 CONTINUE
C
      RETURN
      END
      SUBROUTINE EPSAB(EIGSAB,EIGS,GAB,GA,GB,UA,UB,
     1           UAB,UDMS,NORBS,NCLOSE,IWFLB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /OMVAL/ OMEGA
C
C  THIS SUBROUTINE CREATES THE NEW EPSILON MATRIX AND UDMS MATRIX
C
      DIMENSION EIGS(NORBS),EIGSAB(NORBS,NORBS),
     1          GA(NORBS,NORBS),GB(NORBS,NORBS),
     2          GAB(NORBS,NORBS),UAB(NORBS,NORBS),
     3          UA(NORBS,NORBS),UB(NORBS,NORBS),UDMS(NORBS,NORBS)
C
C  ZERO EPSILON OMEGA OMEGA MATRIX INITIALLY
C
      CALL ZEROM(EIGSAB,NORBS)
C
C  ZERO UAB MINUS OMEGA,OMEGA MATRIX INITIALLY
C
      CALL ZEROM(UDMS,NORBS)
C
      IF ((IWFLB .EQ. 0) .OR. (IWFLB .EQ. 1)) THEN
         OMVAL = 2.0D00 * OMEGA
      ELSEIF (IWFLB .EQ. 3) THEN
         OMVAL = 0.0D00
      ELSEIF (IWFLB .EQ. 2) THEN
         OMVAL = OMEGA
      ENDIF
      DO 30 I = 1,NCLOSE
         DO 20 J = 1,NCLOSE
            S1=0.0D00
            DO 10 K = NCLOSE+1, NORBS
C
C CALCULATION FOR EPSAB
C
               S1 = S1+ GA(I,K)*UB(K,J)+GB(I,K)*UA(K,J)
C
   10       CONTINUE
            EIGSAB(I,J)=GAB(I,J)+S1+UAB(I,J)*(EIGS(I)-
     1                    EIGS(J) + OMVAL)/27.2113961D00
   20    CONTINUE
   30 CONTINUE
C
C CALCULATION FOR UMS
C
      DO 60 I=1,NORBS
         DO 50 J=1,NORBS
            S2=0.0D00
            DO 40 K= 1,NORBS
               S2 = S2 + UA(I,K)*UB(K,J)+UB(I,K)*UA(K,J)
   40       CONTINUE
C
            UDMS(I,J) = S2-UAB(I,J)
   50    CONTINUE
   60 CONTINUE
C
      RETURN
      END
      SUBROUTINE FFREQ1(F, PTOT, PA, PB, NDIM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION F(NDIM,NDIM), PTOT(NDIM,NDIM),
     1          PA(NDIM,NDIM), PB(NDIM,NDIM)
C *********************************************************************
C
C *** COMPUTE THE REMAINING CONTRIBUTIONS TO THE ONE-CENTRE ELEMENTS.
C
C *********************************************************************
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /GAUSS / FN1(107),FN2(107)
     1       /MOLORB/ USPD(MAXORB),DUMY(MAXORB)
      COMMON /TWOELE/ GSS(107),GSP(107),GPP(107),GP2(107),HSP(107)
     1                ,GSD(107),GPD(107),GDD(107)
      COMMON /KEYWRD/ KEYWRD
      CHARACTER*241 KEYWRD
      LOGICAL FIRST
      SAVE FIRST
      DATA FIRST /.TRUE./
      IF(FIRST)THEN
         FIRST=.FALSE.
      ENDIF
      DO 60 II=1,NUMAT
         IA=NFIRST(II)
         IB=NMIDLE(II)
         IC=NLAST(II)
         NI=NAT(II)
         DTPOP=0.D0
         DAPOP=0.D0
         PTPOP=0.D0
         PAPOP=0.D0
         GOTO (60,30,20,20,20,10,10,10,10,10)IC-IA+2
   10    DTPOP=PTOT(IC,IC)+PTOT(IC-1,IC-1)
     1        +PTOT(IC-2,IC-2)+PTOT(IC-3,IC-3)
     2        +PTOT(IC-4,IC-4)
         DAPOP=PA(IC,IC)+PA(IC-1,IC-1)
     1        +PA(IC-2,IC-2)+PA(IC-3,IC-3)
     2        +PA(IC-4,IC-4)
   20    PTPOP=PTOT(IB,IB)+PTOT(IB-1,IB-1)
     1        +PTOT(IB-2,IB-2)
         PAPOP=PA(IB,IB)+PA(IB-1,IB-1)
     1        +PA(IB-2,IB-2)
   30    DBPOP=DTPOP-DAPOP
C
C     F(S,S)
C
         F(IA,IA) = F(IA,IA) + PB(IA,IA)*GSS(NI) + PTPOP*GSP(NI)
     1                       - PAPOP*HSP(NI) + DTPOP*GSD(NI)
         IF (NI.LT.3) GO TO 60
         IPLUS=IA+1
         DO 40 J=IPLUS,IB
C
C     F(P,P)
C
            F(J,J) = F(J,J) + PTOT(IA,IA)*GSP(NI) - PA(IA,IA)*HSP(NI)
     1                   + PB(J,J)*GPP(NI) + (PTPOP-PTOT(J,J))*GP2(NI)
     2                   - 0.5D0*(PAPOP-PA(J,J))*(GPP(NI)-GP2(NI))
     3                   + DTPOP*GPD(NI)
C
C     F(S,P)
C
            F(IA,J) = F(IA,J) + 2.D0*PTOT(IA,J)*HSP(NI)
     1                     - PA(IA,J)*(HSP(NI)+GSP(NI))
            F(J,IA) = F(J,IA) + 2.D0*PTOT(J,IA)*HSP(NI)
     1                     - PA(J,IA)*(HSP(NI)+GSP(NI))
   40    CONTINUE
C
C     F(P,P*)
C
         IMINUS=IB-1
         DO 50 J=IPLUS,IMINUS
            ICC=J+1
            DO 50 L=ICC,IB
               F(J,L) = F(J,L) + PTOT(J,L)*(GPP(NI)-GP2(NI))
     1                         - 0.5D0*PA(J,L)*(GPP(NI)+GP2(NI))
               F(L,J) = F(L,J) + PTOT(L,J)*(GPP(NI)-GP2(NI))
     1                         - 0.5D0*PA(L,J)*(GPP(NI)+GP2(NI))
   50    CONTINUE
   60 CONTINUE
      RETURN
      END
      SUBROUTINE FFREQ2(F,PTOT,W,NUMAT,NFIRST,NLAST,NORBS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION F(NORBS,NORBS), PTOT(NORBS,NORBS),
     1           NFIRST(*),
     2          NLAST(*), W(*)
C*******************************************************************
C
C  TDHF FORMS TWO ELECTRON TWO CENTER REPULSION PART OF THE FOCK
C  MATRIX
C ON INPUT PTOT = TOTAL DENSITY MATRIX
C          P    = ALPHA OR BETA DENSITY MATRIX
C          W    = TWO ELECTRON INTEGRAL MATRIX
C
C  ON OUTPUT F = PARTIAL FOCK MATRIX
C
C********************************************************************
      COMMON /KEYWRD/ KEYWRD
      CHARACTER*241 KEYWRD
   10 KK=0
C
      DO 70 II=2,NUMAT
         IIM1 = II - 1
         IA=NFIRST(II)
         IB=NLAST(II)
         DO 60 JJ=1,IIM1
            JA=NFIRST(JJ)
            JB=NLAST(JJ)
            DO 50 I=IA,IB
               DO 40 J=IA,I
                  FIJ = 1.0D00
                  IF (I.EQ.J) FIJ = 0.5D00
                  DO 30 K=JA,JB
                     DO 20 L=JA,K
                        FKL = 1.0D00
                        IF (K.EQ.L) FKL = 0.5D00
                        KK=KK+1
                        A=W(KK)
                        AINT=A*FKL*FIJ
                        F(I,J)=F(I,J)+ AINT*(PTOT(K,L)+PTOT(L,K))
                        F(J,I)=F(J,I)+ AINT*(PTOT(K,L)+PTOT(L,K))
                        F(K,L)=F(K,L)+ AINT*(PTOT(I,J)+PTOT(J,I))
                        F(L,K)=F(L,K)+ AINT*(PTOT(I,J)+PTOT(J,I))
                        AINT=AINT*0.5D00
                        F(I,L)=F(I,L)-AINT*PTOT(J,K)
                        F(L,I)=F(L,I)-AINT*PTOT(K,J)
                        F(K,J)=F(K,J)-AINT*PTOT(L,I)
                        F(J,K)=F(J,K)-AINT*PTOT(I,L)
                        F(I,K)=F(I,K)-AINT*PTOT(J,L)
                        F(K,I)=F(K,I)-AINT*PTOT(L,J)
                        F(J,L)=F(J,L)-AINT*PTOT(I,K)
                        F(L,J)=F(L,J)-AINT*PTOT(K,I)
   20                CONTINUE
   30             CONTINUE
   40          CONTINUE
   50       CONTINUE
   60    CONTINUE
   70 CONTINUE
      RETURN
      END
      SUBROUTINE FHPATN(A,B,NORBS,ITW,SIGN)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C THIS SUBROUTINE CONVERTS THE MATRICES INTO ITS ADJOINTS
C
      DIMENSION A(NORBS,NORBS),B(NORBS,NORBS)
      GO TO (10,40,40,10) ITW
   10 CONTINUE
      DO 30 J=1,NORBS
         DO 20 I=1,NORBS
            A(I,J)=B(I,J)
   20    CONTINUE
   30 CONTINUE
      GO TO 70
   40 CONTINUE
      DO 60 J=1,NORBS
         DO 50 I=1,NORBS
            A(I,J)=SIGN*B(J,I)
   50    CONTINUE
   60 CONTINUE
   70 CONTINUE
      RETURN
      END
      SUBROUTINE HMUF(H1,ID,COORD,NFIRST,NLAST,NAT,NORBS,NUMAT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C  FORM THE DIPOLE MOMENT MATRIX FOR COMPONENT ID
C
      COMMON /MULTIP/ DD(107),QQ(107),AM(107),AD(107),AQ(107)
      DIMENSION H1(NORBS,NORBS), COORD(3,*)
      DIMENSION NFIRST(NUMAT),NLAST(NUMAT),NAT(NUMAT)
C
C  ZERO H1 MATRIX
C
      CALL ZEROM(H1,NORBS)
C
C  FORM DIPOLE MATRIX
C
      DO 30 I = 1,NUMAT
         IA = NFIRST(I)
         IB = NLAST(I)
         NI = NAT(I)
         DO 20 I1 = IA,IB
            DO 10 J1 = IA,I1
               H1(I1,J1) = 0.0D00
               IO1 = I1 - IA
               JO1 = J1 - IA
               IF ((ID.EQ.1).AND.((JO1.EQ.0).AND.(IO1.EQ.1))) THEN
                  H1(I1,J1) = DD(NI)
                  H1(J1,I1) = DD(NI)
               ENDIF
               IF ((ID.EQ.2).AND.((JO1.EQ.0).AND.(IO1.EQ.2))) THEN
                  H1(I1,J1) = DD(NI)
                  H1(J1,I1) = DD(NI)
               ENDIF
               IF ((ID.EQ.3).AND.((JO1.EQ.0).AND.(IO1.EQ.3))) THEN
                  H1(I1,J1) = DD(NI)
                  H1(J1,I1) = DD(NI)
               ENDIF
   10       CONTINUE
            H1(I1,I1) = 0.0D00
C.. ADDED FOR TRANSLATION OF CENTER FROM ORIGIN
            H1(I1,I1) = H1(I1,I1) + 1.8897262D0*COORD(ID,I)
   20    CONTINUE
   30 CONTINUE
C
      RETURN
      END
      SUBROUTINE HPLUSF(F,H,NORBS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C HPLUSF ADDS THE 1 AND 2-ELECTRON PARTS OF THE FOCK MATRIX
C
      DIMENSION F(NORBS,NORBS), H(NORBS,NORBS)
      DO 20 I=1,NORBS
         DO 10 J=1,NORBS
            TERM = F(I,J)/27.2113961D00
            F(I,J)=H(I,J)+TERM
   10    CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE MAKEUF(U,UOLD,G,EIGS,LAST,NORBS,NNORBS,NCLOSE,
     1                   DIFF,ATOL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL LAST
C
C  THIS SUBROUTINE CREATES THE NEW TRANSFORMATION MATRIX U
C  AND THEN CHECKS FOR CONVERGENCE
C
      DIMENSION U(NORBS,NORBS),UOLD(NORBS,NORBS),
     1          EIGS(NORBS),G(NORBS,NORBS)
      COMMON /OMVAL/ OMEGA
C
C  ZERO MATRIX INITIALLY
C
      CALL ZEROM(U,NORBS)
C
C  CREATE OFF-DIAGONAL BLOCKS
C
      DO 20 K = NCLOSE+1,NORBS
         DO 10 L = 1,NCLOSE
            U(L,K) = 27.2113961D0*G(L,K)/(EIGS(K)-EIGS(L)-OMEGA)
            U(K,L) = 27.2113961D0*G(K,L)/(EIGS(L)-EIGS(K)-OMEGA)
   10    CONTINUE
   20 CONTINUE
C
C  CHECK FOR CONVERGENCE
C
      DIFF = 0.0D00
      DO 40 I = 1,NORBS
         DO 30 J = 1,NORBS
            UDIF = ABS(U(I,J)-UOLD(I,J))
            IF (DIFF.LT.UDIF) DIFF = UDIF
   30    CONTINUE
   40 CONTINUE
      IF (DIFF.LT.ATOL) THEN
         LAST = .TRUE.
      ENDIF
C
      DO 60 I = 1,NORBS
         DO 50 J = 1,NORBS
            UOLD(I,J) = U(I,J)
   50    CONTINUE
   60 CONTINUE
C
      RETURN
      END
      SUBROUTINE NGAMTG(IGAM,X,GD3,UD3,G1,U1,GS,USMD,EPS,US)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
C.....................................................................
C  CALCULATE GAMMA(THG) IN A NONITERATIVE FASHION
C.....................................................................
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM),NORBS,NELECS,NALPHA,NBETA,NCLOSE,
     2                NOPEN,NDUMY,FRACT
      COMMON /OMVAL/ OMEGA
      COMMON /CHANEL/ IFILES(30)
      EQUIVALENCE(IW,IFILES(6))
      DIMENSION X(MAXORB,MAXORB),GD3(MAXORB,MAXORB),UD3(MAXORB,MAXORB),
     1          G1(MAXORB,MAXORB),U1(MAXORB,MAXORB),GS(MAXORB,MAXORB),
     2          USMD(MAXORB,MAXORB),EPS(MAXORB,MAXORB),
     3          GAMMA(9),US(MAXORB,MAXORB)
      CHARACTER*1 ALAB
      DIMENSION ALAB(3),IDA(9),IDB(9),IDC(9),IDD(9),
     1          IPAIR(3,3)
      SAVE ALAB, IDA, IDB, IDC, IDD, IPAIR
      DATA ALAB /'X','Y','Z'/
      DATA IDA /1,2,3,1,1,2,2,3,3/
      DATA IDB /1,2,3,1,1,2,2,3,3/
      DATA IDC /1,2,3,2,3,1,3,1,2/
      DATA IDD /1,2,3,2,3,1,3,1,2/
      DATA IPAIR /1,2,3,2,4,5,3,5,6/
      ONE=1.D0
      MSQ = NORBS*NORBS
   10 WRITE(6,20) OMEGA
   20 FORMAT (//,
     1   ' GAMMA (THIRD HARMONIC GENERATION) AT ',F10.5,' EV.'//)
C
C IGAM=1 (THIRD HARMONIC GENERATION)
C
      JGARC=22
      JUARC=19
      JUREC=07
      JGREC=10
      JG2REC=55
      JU2REC=49
      JU2MRC=67
      JEPREC=61
C
C LOOP BEGINS FOR THE CALCULATION OF GAMMA(ABCD)
C
      GAV = 0.0D+00
      DO 100 IE=1,9
C
         IA=IDA(IE)
         IB=IDB(IE)
         IC=IDC(IE)
         ID=IDD(IE)
         ICD=IPAIR(IC,ID)
         IBD=IPAIR(IB,ID)
         IBC=IPAIR(IB,IC)
C
C  READ IN THE FIRST ORDER U3 OMEGA AND G3 OMEGA IN THE DIRECTION A
C
C MAKE GD3 OMEGA MATRIX FROM G3 MATRIX
C
         CALL DAREAD(X,MSQ,JGARC+IA)
         CALL FHPATN(GD3,X,NORBS,2,ONE)
C
C MAKE UD3 OMEGA MATRIX FROM U3 OMEGA MATRIX
C
         CALL DAREAD(X,MSQ,JUARC+IA)
         CALL FHPATN(UD3,X,NORBS,2,-ONE)
C
         YY=0.0D00
         IMOVE=1
   30    CONTINUE
C
   40    GO TO (50,60,70), IMOVE
C
   50    J2=IB
         J34=ICD
         GO TO 80
   60    J2=IC
         J34=IBD
         GO TO 80
   70    J2=ID
         J34=IBC
   80    CONTINUE
C
C  READ IN G1,U1,GS,US,UMS,EPS
C
C  GET  UB
         CALL DAREAD(U1,MSQ,JUREC+J2)
C  GET  GB
         CALL DAREAD(G1,MSQ,JGREC+J2)
C  GET  GCD
         CALL DAREAD(GS,MSQ,JG2REC+J34)
C  GET  UCD
         CALL DAREAD(US,MSQ,JU2REC+J34)
C  GET  USMD
         CALL DAREAD(USMD,MSQ,JU2MRC+J34)
C  GET  EPCD
         CALL DAREAD(EPS,MSQ,JEPREC+J34)
C
C
C FIRST KIND
C
         YY = YY+TRSUB(UD3,G1,US,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(USMD,G1,UD3,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(UD3,G1,US,NORBS,NCLOSE,NORBS)
         YY = YY+TRSUB(USMD,G1,UD3,NORBS,NCLOSE,NORBS)
C
C SECOND KIND
C
         YY = YY+TRSUB(UD3,GS,U1,NCLOSE,NORBS,NORBS)
         YY = YY+TRSUB(U1,GS,UD3,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(UD3,EPS,U1,NORBS,NCLOSE,NORBS)
         YY = YY-TRSUB(U1,EPS,UD3,NORBS,NCLOSE,NORBS)
C
C THIRD KIND
C
         YY = YY+TRSUB(U1,GD3,US,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(USMD,GD3,U1,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(U1,GD3,US,NORBS,NCLOSE,NORBS)
         YY = YY+TRSUB(USMD,GD3,U1,NORBS,NCLOSE,NORBS)
C
         IMOVE=IMOVE+1
         IF (IMOVE .LE. 3) GO TO 30
C
         GAMMA(IE) = YY
C
C CALCULATE THE AVERAGE GAMMA VALUE
C
         GAV = GAV + YY
C
C WRITE GAMMA(ABCD)
C
         WRITE(6,90)ALAB(IA),ALAB(IB),ALAB(IC),ALAB(ID),GAMMA(IE)
   90    FORMAT(' GAMMA(',A1,',',A1,',',A1,',',A1,') = ',1F13.5)
C
  100 CONTINUE
      GAVE = GAV/5.0D+00
      WRITE(6,110) OMEGA,GAVE
  110 FORMAT(//,' AVERAGE GAMMA VALUE AT ',F10.5,' = ',1F13.5,//)
      RETURN
      END
      SUBROUTINE NGEFIS(IGAM,X,GD3,UD3,G1,U1,GS,USMD,EPS,US)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
C.....................................................................
C  CALCULATE GAMMA(DC-EFISHG) IN A NONITERATIVE FASHION
C.....................................................................
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM),NORBS,NELECS,NALPHA,NBETA,NCLOSE,
     2                NOPEN,NDUMY,FRACT
      COMMON /OMVAL/ OMEGA
      COMMON /IODAF/IDAF,IRECLN,IRECST,IFILEN(145),IODA(145)
      COMMON /CHANEL/ IFILES(30)
      EQUIVALENCE(IW,IFILES(6))
      DIMENSION X(MAXORB,MAXORB),GD3(MAXORB,MAXORB),UD3(MAXORB,MAXORB),
     1          G1(MAXORB,MAXORB),U1(MAXORB,MAXORB),GS(MAXORB,MAXORB),
     2          USMD(MAXORB,MAXORB),EPS(MAXORB,MAXORB),
     3          GAMMA(15),US(MAXORB,MAXORB)
      CHARACTER*1 ALAB
      DIMENSION ALAB(3),IDA(15),IDB(15),IDC(15),IDD(15),
     1          IP(3,3),IPAIR(3,3)
      SAVE ALAB, IDA, IDB, IDC, IDD, IP, IPAIR
      DATA ALAB /'X','Y','Z'/
C
      DATA IDA /1,2,3,1,1,2,2,3,3,1,1,2,2,3,3/
      DATA IDB /1,2,3,2,3,1,3,1,2,1,1,2,2,3,3/
      DATA IDC /1,2,3,1,1,2,2,3,3,2,3,1,3,1,2/
      DATA IDD /1,2,3,2,3,1,3,1,2,2,3,1,3,1,2/
      DATA IP /1,2,3,2,4,5,3,5,6/
      DATA IPAIR /1,4,7,2,5,8,3,6,9/
      ONE=1.D0
      MSQ = NORBS*NORBS
      WRITE(6,10) OMEGA
   10 FORMAT (//,
     1   ' GAMMA (DC-EFISHG) AT ',F10.5,' EV.'//)
C
C GET DATA FROM ALPHA  AND ITERATIVE BETA CALCULATIONS
C
C   REQUIRED RECORDS FROM POLARIZABILITY CALCULATIONS
C   -------------------------------------------------------
C        0    W    2W    3W
C
C       -02- -08-  -14-  -20-  -U- MATRIX FOR -X- DIRECTION
C       -03- -09-  -15-  -21-  -U- MATRIX FOR -Y- DIRECTION
C       -04- -10-  -16-  -22-  -U- MATRIX FOR -Z- DIRECTION
C       -05- -11-  -17-  -23-  -G- MATRIX FOR -X- DIRECTION
C       -06- -12-  -18-  -24-  -G- MATRIX FOR -Y- DIRECTION
C       -07- -13-  -19-  -25-  -G- MATRIX FOR -Z- DIRECTION
C   -------------------------------------------------------
C      (0,0)    (W,W)    (0,W)    (W,-W)
C
C      -26-     -50-     -74-     -110-   -U- MATRIX FOR -XX- DIRECTION
C      -27-     -51-     -75-     -111-   -U- MATRIX FOR -XY- DIRECTION
C      -28-     -52-     -76-     -112-   -U- MATRIX FOR -XZ- DIRECTION
C                        -77-     -113-   -U- MATRIX DOE -YX- DIRECTION
C      -29-     -53-     -78-     -114-   -U- MATRIX FOR -YY- DIRECTION
C      -30-     -54-     -79-     -115-   -U- MATRIX FOR -YZ- DIRECTION
C                        -80-     -116-   -U- MATRIX FOR -ZX- DIRECTION
C                        -81-     -117-   -U- MATRIX FOR -ZY- DIRECTION
C      -31-     -55-     -82-     -118-   -U- MATRIX FOR -ZZ- DIRECTION
C   ------------------------------------------------------------------
C
C      -32-     -56-     -83-     -119-   -G- MATRIX FOR -XX- DIRECTION
C      -33-     -57-     -84-     -120-   -G- MATRIX FOR -XY- DIRECTION
C      -34-     -58-     -85-     -121-   -G- MATRIX FOR -XZ- DIRECTION
C                        -86-     -122-   -G- MATRIX FOR -YX- DIRECTION
C      -35-     -59-     -87-     -123-   -G- MATRIX FOR -YY- DIRECTION
C      -36-     -60-     -88-     -124-   -G- MATRIX FOR -YZ- DIRECTION
C                        -89-     -125-   -G- MATRIX FOR -ZX- DIRECTION
C                        -90-     -126-   -G- MATRIX FOR -ZY- DIRECTION
C      -37-     -61-     -91-     -127-   -G- MATRIX FOR -ZZ- DIRECTION
C   ------------------------------------------------------------------
C
C      -38-     -62-     -92-     -128-   -E- MATRIX FOR -XX- DIRECTION
C      -39-     -63-     -93-     -129-   -E- MATRIX FOR -XY- DIRECTION
C      -40-     -64-     -94-     -130-   -E- MATRIX FOR -XZ- DIRECTION
C                        -95-     -131-   -E- MATRIX FOR -YX- DIRECTION
C      -41-     -65-     -96-     -132-   -E- MATRIX FOR -YY- DIRECTION
C      -42-     -66-     -97-     -133-   -E- MATRIX FOR -YZ- DIRECTION
C                        -98-     -134-   -E- MATRIX FOR -ZX- DIRECTION
C                        -99-     -135-   -E- MATRIX FOR -ZY- DIRECTION
C      -43-     -67-     -100-    -136-   -E- MATRIX FOR -ZZ- DIRECTION
C   ------------------------------------------------------------------
C
C      -44-     -68-     -101-    -137-   -UM- MATRIX FOR -XX- DIRECTION
C      -45-     -69-     -102-    -138-   -UM- MATRIX FOR -XY- DIRECTION
C      -46-     -70-     -103-    -139-   -UM- MATRIX FOR -XZ- DIRECTION
C                        -104-    -140-   -UM- MATRIX FOR -YX- DIRECTION
C      -47-     -71-     -105-    -141-   -UM- MATRIX FOR -YY- DIRECTION
C      -48-     -72-     -106-    -142-   -UM- MATRIX FOR -YZ- DIRECTION
C                        -107-    -143-   -UM- MATRIX FOR -ZX- DIRECTION
C                        -108-    -144-   -UM- MATRIX FOR -ZY- DIRECTION
C      -49-     -73-     -109-    -145-   -UM- MATRIX FOR -ZZ- DIRECTION
C   ------------------------------------------------------------------
C
C CALCULATION OF DIFFERENT GAMMA VALUES IN A NONITERATIVE METHOD.
C
C IGAM=2 (DC-ELECTIC FIELD INDUCED SECOND HARMONIC GENERATION)
C
      JGARC=16
      JUARC=13
      JUREC=01
      JGREC=04
      JG2REC=55
      JU2REC=49
      JU2MRC=67
      JEPREC=61
C LOOP BEGINS FOR THE CALCULATION OF GAMMA(ABCD)
C
      GAV = 0.0D00
      DO 90 IE=1,15
         IA=IDA(IE)
         IB=IDB(IE)
         IC=IDC(IE)
         ID=IDD(IE)
         ICD=IP(IC,ID)
         IBD=IPAIR(IB,ID)
         IBC=IPAIR(IB,IC)
C
C  READ IN THE FIRST ORDER U3 OMEGA AND G3 OMEGA IN THE DIRECTION A
C  MAKE GD3 OMEGA MATRIX FROM G3 MATRIX
C
         CALL DAREAD(X,MSQ,JGARC+IA)
         CALL FHPATN(GD3,X,NORBS,2,ONE)
C
C MAKE UD3 OMEGA MATRIX FROM U3 OMEGA MATRIX
C
         CALL DAREAD(X,MSQ,JUARC+IA)
         CALL FHPATN(UD3,X,NORBS,2,-ONE)
         YY=0.0D00
         IMOVE=1
   20    CONTINUE
C
C DC EFISHG
C
   30    GO TO (40,50,60), IMOVE
   40    J2=IB
         J34=ICD
         GO TO 70
   50    J2=IC+6
C                 J34=IBD+24
         J3U=IBD+24
         J3G=IBD+27
         J3E=IBD+30
         J3UM=IBD+33
         GO TO 70
   60    J2=ID+6
         J3U=IBC+24
         J3G=IBC+27
         J3E=IBC+30
         J3UM=IBC+33
   70    CONTINUE
C
C  READ IN G1,U1,GS,US,UMS,EPS
C
C  CALL UB
C
         CALL DAREAD(U1,MSQ,JUREC+J2)
C  CALL GB
         CALL DAREAD(G1,MSQ,JGREC+J2)
         IF (IMOVE .EQ. 1) THEN
C  CALL GCD
            CALL DAREAD(GS,MSQ,JG2REC+J34)
C  CALL UCD
            CALL DAREAD(US,MSQ,JU2REC+J34)
C  CALL USMD
            CALL DAREAD(USMD,MSQ,JU2MRC+J34)
C  CALL EPCD
            CALL DAREAD(EPS,MSQ,JEPREC+J34)
C
         ELSE
C  CALL GCD
            CALL DAREAD(GS,MSQ,JG2REC+J3G)
C  CALL UCD
            CALL DAREAD(US,MSQ,JU2REC+J3U)
C  CALL USMD
            CALL DAREAD(USMD,MSQ,JU2MRC+J3UM)
C  CALL EPCD
            CALL DAREAD(EPS,MSQ,JEPREC+J3E)
         ENDIF
C
C FIRST KIND
C
         YY = YY+TRSUB(UD3,G1,US,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(USMD,G1,UD3,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(UD3,G1,US,NORBS,NCLOSE,NORBS)
         YY = YY+TRSUB(USMD,G1,UD3,NORBS,NCLOSE,NORBS)
C
C SECOND KIND
C
         YY = YY+TRSUB(UD3,GS,U1,NCLOSE,NORBS,NORBS)
         YY = YY+TRSUB(U1,GS,UD3,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(UD3,EPS,U1,NORBS,NCLOSE,NORBS)
         YY = YY-TRSUB(U1,EPS,UD3,NORBS,NCLOSE,NORBS)
C
C THIRD KIND
C
         YY = YY+TRSUB(U1,GD3,US,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(USMD,GD3,U1,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(U1,GD3,US,NORBS,NCLOSE,NORBS)
         YY = YY+TRSUB(USMD,GD3,U1,NORBS,NCLOSE,NORBS)
C
         IMOVE=IMOVE+1
         IF (IMOVE .LE. 3) GO TO 20
C
         GAMMA(IE) = YY
C CALCULATE THE AVERAGE GAMMA VALUE
         IF (IE .LE. 3) THEN
            GAV = GAV + 3*YY
         ELSEIF (IE .GT. 9) THEN
            GAV = GAV + YY
         ELSE
            GAV = GAV + 2*YY
         ENDIF
C
C WRITE GAMMA(ABCD)
C
         WRITE(6,80) ALAB(IA),ALAB(IB),ALAB(IC),ALAB(ID),GAMMA(IE)
   80    FORMAT(' GAMMA(',A1,',',A1,',',A1,',',A1,') = ',1PD14.7)
C
   90 CONTINUE
      GAVE = GAV/15.0D+00
      WRITE(6,100) OMEGA,GAVE
  100 FORMAT(//,' AVERAGE GAMMA VALUE AT ',F10.5, ' EV = ',1PD14.7,//)
      RETURN
      END
      SUBROUTINE NGIDRI(IGAM,X,GD3,UD3,G1,U1,GS,USMD,EPS,US)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
C.....................................................................
C  CALCULATE GAMMA(IDRI) IN A NONITERATIVE FASHION
C.....................................................................
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM),NORBS,NELECS,NALPHA,NBETA,NCLOSE,
     2                NOPEN,NDUMY,FRACT
      COMMON /OMVAL/ OMEGA
      COMMON /IODAF/IDAF,IRECLN,IRECST,IFILEN(145),IODA(145)
      COMMON /CHANEL/ IFILES(30)
      EQUIVALENCE(IW,IFILES(6))
      DIMENSION X(MAXORB,MAXORB),GD3(MAXORB,MAXORB),UD3(MAXORB,MAXORB),
     1          G1(MAXORB,MAXORB),U1(MAXORB,MAXORB),GS(MAXORB,MAXORB),
     2          USMD(MAXORB,MAXORB),EPS(MAXORB,MAXORB),
     3          GAMMA(15),US(MAXORB,MAXORB)
      DIMENSION ALAB(3),IDA(15),IDB(15),IDC(15),IDD(15),
     1          IP(3,3),IPAIR(3,3)
      CHARACTER*1 ALAB
      SAVE ALAB, IDA, IDB, IDC, IDD, IP, IPAIR
      DATA ALAB /'X','Y','Z'/
C
      DATA IDA /1,2,3,1,1,2,2,3,3,1,1,2,2,3,3/
      DATA IDB /1,2,3,1,1,2,2,3,3,2,3,1,3,1,2/
      DATA IDC /1,2,3,2,3,1,3,1,2,2,3,1,3,1,2/
      DATA IDD /1,2,3,2,3,1,3,1,2,1,1,2,2,3,3/
      DATA IP /1,2,3,2,4,5,3,5,6/
      DATA IPAIR /1,4,7,2,5,8,3,6,9/
      ONE=1.0D00
      MSQ = NORBS*NORBS
C
      WRITE(6,10) OMEGA
   10 FORMAT (//,
     1              ' GAMMA (IDRI) AT ',F10.5,' EV.'//)
C
C GET DATA FROM ALPHA  AND ITERATIVE BETA CALCULATIONS
C
C
C IGAM=3 (INTENSITY DEPENDENT REFRACTIVE INDEX OR DEGENERATED FOUR
C WAVE MIXING)
C
      JGARC=10
      JUARC=07
      JUREC=07
      JGREC=10
C LOOP BEGINS FOR THE CALCULATION OF GAMMA(ABCD)
C
      GAV = 0.0D+00
      DO 80 IE=1,15
         IA=IDA(IE)
         IB=IDB(IE)
         IC=IDC(IE)
         ID=IDD(IE)
         ICD=IPAIR(IC,ID)
         IBD=IPAIR(IB,ID)
         IBC=IP(IB,IC)
C
C  READ IN THE FIRST ORDER U3 OMEGA AND G3 OMEGA IN THE DIRECTION A
C
C MAKE GD3 OMEGA MATRIX FROM G3 MATRIX
C
         CALL DAREAD(X,MSQ,JGARC+IA)
         CALL FHPATN(GD3,X,NORBS,2,ONE)
C
C MAKE UD3 OMEGA MATRIX FROM U3 OMEGA MATRIX
C
         CALL DAREAD(X,MSQ,JUARC+IA)
         CALL FHPATN(UD3,X,NORBS,2,-ONE)
C
         YY=0.0D00
         IMOVE=1
   20    CONTINUE
C
C
C IDRI
C
         GO TO (30,40,50), IMOVE
   30    J2= IB
         J34=ICD
         JG2REC=118
         JU2REC=109
         JU2MRC=136
         JEPREC=127
         GO TO 60
   40    J2=IC
         J34=IBD
         JG2REC=118
         JU2REC=109
         JU2MRC=136
         JEPREC=127
         GO TO 60
   50    J2=ID
         J34=IBC
         JG2REC= 55
         JU2REC= 49
         JU2MRC= 67
         JEPREC= 61
   60    CONTINUE
C  READ IN G1,U1,GS,US,UMS,EPS
C
C  CALL UB
         IF (IMOVE.EQ.3) THEN
            CALL DAREAD(X,MSQ,JUREC+J2)
            CALL FHPATN(U1,X,NORBS,2,-ONE)
         ELSE
            CALL DAREAD(U1,MSQ,JUREC+J2)
         ENDIF
C  CALL GB
         IF (IMOVE.EQ.3) THEN
            CALL DAREAD(X,MSQ,JGREC+J2)
            CALL FHPATN(G1,X,NORBS,2,ONE)
         ELSE
            CALL DAREAD(G1,MSQ,JGREC+J2)
         ENDIF
C  CALL GCD
         CALL DAREAD(GS,MSQ,JG2REC+J34)
C  CALL UCD
         CALL DAREAD(US,MSQ,JU2REC+J34)
C  CALL USMD
         CALL DAREAD(USMD,MSQ,JU2MRC+J34)
C  CALL EPCD
         CALL DAREAD(EPS,MSQ,JEPREC+J34)
C
C FIRST KIND
C
         YY = YY+TRSUB(UD3,G1,US,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(USMD,G1,UD3,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(UD3,G1,US,NORBS,NCLOSE,NORBS)
         YY = YY+TRSUB(USMD,G1,UD3,NORBS,NCLOSE,NORBS)
C
C SECOND KIND
C
         YY = YY+TRSUB(UD3,GS,U1,NCLOSE,NORBS,NORBS)
         YY = YY+TRSUB(U1,GS,UD3,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(UD3,EPS,U1,NORBS,NCLOSE,NORBS)
         YY = YY-TRSUB(U1,EPS,UD3,NORBS,NCLOSE,NORBS)
C
C THIRD KIND
C
         YY = YY+TRSUB(U1,GD3,US,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(USMD,GD3,U1,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(U1,GD3,US,NORBS,NCLOSE,NORBS)
         YY = YY+TRSUB(USMD,GD3,U1,NORBS,NCLOSE,NORBS)
C
         IMOVE=IMOVE+1
         IF (IMOVE .LE. 3) GO TO 20
C
         GAMMA(IE) = YY
C
C CALCULATE THE AVERAGE GAMMA VALUE
C
         IF (IE .LE. 3) THEN
            GAV = GAV + 3.0D0*YY
         ELSEIF (IE .GT. 9) THEN
            GAV = GAV + YY
         ELSE
            GAV = GAV + 2.0D0*YY
         ENDIF
C
C WRITE GAMMA(ABCD)
C
         WRITE(6,70) ALAB(IA),ALAB(IB),ALAB(IC),ALAB(ID),GAMMA(IE)
   70    FORMAT(' GAMMA(',A1,',',A1,',',A1,',',A1,') = ',1PD14.7)
C
   80 CONTINUE
      GAVE = GAV/15.0D+00
      WRITE(6,90) OMEGA,GAVE
   90 FORMAT(//,'  AVERAGE GAMMA VALUE AT ',F10.5,' = ',1PD14.7,//)
      RETURN
      END
      SUBROUTINE NGOKE(IGAM,X,GD3,UD3,G1,U1,GS,USMD,EPS,US)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
C.....................................................................
C  CALCULATE GAMMA(OKE) IN A NONITERATIVE FASHION
C.....................................................................
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM),NORBS,NELECS,NALPHA,NBETA,NCLOSE,
     2                NOPEN,NDUMY,FRACT
      COMMON /OMVAL/ OMEGA
      COMMON /IODAF/IDAF,IRECLN,IRECST,IFILEN(145),IODA(145)
      COMMON /CHANEL/ IFILES(30)
      EQUIVALENCE(IW,IFILES(6))
      DIMENSION X(MAXORB,MAXORB),GD3(MAXORB,MAXORB),UD3(MAXORB,MAXORB),
     1          G1(MAXORB,MAXORB),U1(MAXORB,MAXORB),GS(MAXORB,MAXORB),
     2          USMD(MAXORB,MAXORB),EPS(MAXORB,MAXORB),
     3          GAMMA(15),US(MAXORB,MAXORB)
      CHARACTER*1 ALAB
      DIMENSION ALAB(3),IDA(15),IDB(15),IDC(15),IDD(15),
     1          IP(3,3),IPAIR(3,3)
      SAVE ALAB, IDA, IDB, IDC, IDD, IP, IPAIR
      DATA ALAB /'X','Y','Z'/
C
      DATA IDA /1,2,3,1,1,2,2,3,3,1,1,2,2,3,3/
      DATA IDB /1,2,3,1,1,2,2,3,3,2,3,1,3,1,2/
      DATA IDC /1,2,3,2,3,1,3,1,2,2,3,1,3,1,2/
      DATA IDD /1,2,3,2,3,1,3,1,2,1,1,2,2,3,3/
      DATA IP /1,2,3,2,4,5,3,5,6/
      DATA IPAIR /1,4,7,2,5,8,3,6,9/
      ONE=1.0D00
      MSQ = NORBS*NORBS
C
      IF (IGAM .EQ. 3 ) THEN
         WRITE(6,10) OMEGA
   10    FORMAT (//,
     1              ' GAMMA (IDRI) AT ',F10.5,' EV.'//)
      ELSE
         WRITE(6,20) OMEGA
   20    FORMAT (//,
     1              ' GAMMA (OKE) AT ',F10.5,' EV.'//)
      ENDIF
C
C DATA INCLUDING YX, ZY, ZX DIRECTIONS
C GET DATA FROM ALPHA  AND ITERATIVE BETA CALCULATIONS
C
C   REQUIRED RECORDS FROM POLARIZABILITY CALCULATIONS
C
C   -------------------------------------------------------
C        0    W    2W    3W
C
C       -02- -08-  -14-  -20-  -U- MATRIX FOR -X- DIRECTION
C       -03- -09-  -15-  -21-  -U- MATRIX FOR -Y- DIRECTION
C       -04- -10-  -16-  -22-  -U- MATRIX FOR -Z- DIRECTION
C       -05- -11-  -17-  -23-  -G- MATRIX FOR -X- DIRECTION
C       -06- -12-  -18-  -24-  -G- MATRIX FOR -Y- DIRECTION
C       -07- -13-  -19-  -25-  -G- MATRIX FOR -Z- DIRECTION
C   -------------------------------------------------------
C      (0,0)    (W,W)    (0,W)    (W,-W)
C
C      -26-     -50-     -74-     -110-   -U- MATRIX FOR -XX- DIRECTION
C      -27-     -51-     -75-     -111-   -U- MATRIX FOR -XY- DIRECTION
C      -28-     -52-     -76-     -112-   -U- MATRIX FOR -XZ- DIRECTION
C                        -77-     -113-   -U- MATRIX DOE -YX- DIRECTION
C      -29-     -53-     -78-     -114-   -U- MATRIX FOR -YY- DIRECTION
C      -30-     -54-     -79-     -115-   -U- MATRIX FOR -YZ- DIRECTION
C                        -80-     -116-   -U- MATRIX FOR -ZX- DIRECTION
C                        -81-     -117-   -U- MATRIX FOR -ZY- DIRECTION
C      -31-     -55-     -82-     -118-   -U- MATRIX FOR -ZZ- DIRECTION
C   ------------------------------------------------------------------
C
C      -32-     -56-     -83-     -119-   -G- MATRIX FOR -XX- DIRECTION
C      -33-     -57-     -84-     -120-   -G- MATRIX FOR -XY- DIRECTION
C      -34-     -58-     -85-     -121-   -G- MATRIX FOR -XZ- DIRECTION
C                        -86-     -122-   -G- MATRIX FOR -YX- DIRECTION
C      -35-     -59-     -87-     -123-   -G- MATRIX FOR -YY- DIRECTION
C      -36-     -60-     -88-     -124-   -G- MATRIX FOR -YZ- DIRECTION
C                        -89-     -125-   -G- MATRIX FOR -ZX- DIRECTION
C                        -90-     -126-   -G- MATRIX FOR -ZY- DIRECTION
C      -37-     -61-     -91-     -127-   -G- MATRIX FOR -ZZ- DIRECTION
C   ------------------------------------------------------------------
C
C      -38-     -62-     -92-     -128-   -E- MATRIX FOR -XX- DIRECTION
C      -39-     -63-     -93-     -129-   -E- MATRIX FOR -XY- DIRECTION
C      -40-     -64-     -94-     -130-   -E- MATRIX FOR -XZ- DIRECTION
C                        -95-     -131-   -E- MATRIX FOR -YX- DIRECTION
C      -41-     -65-     -96-     -132-   -E- MATRIX FOR -YY- DIRECTION
C      -42-     -66-     -97-     -133-   -E- MATRIX FOR -YZ- DIRECTION
C                        -98-     -134-   -E- MATRIX FOR -ZX- DIRECTION
C                        -99-     -135-   -E- MATRIX FOR -ZY- DIRECTION
C      -43-     -67-     -100-    -136-   -E- MATRIX FOR -ZZ- DIRECTION
C   ------------------------------------------------------------------
C
C      -44-     -68-     -101-    -137-   -UM- MATRIX FOR -XX- DIRECTION
C      -45-     -69-     -102-    -138-   -UM- MATRIX FOR -XY- DIRECTION
C      -46-     -70-     -103-    -139-   -UM- MATRIX FOR -XZ- DIRECTION
C                        -104-    -140-   -UM- MATRIX FOR -YX- DIRECTION
C      -47-     -71-     -105-    -141-   -UM- MATRIX FOR -YY- DIRECTION
C      -48-     -72-     -106-    -142-   -UM- MATRIX FOR -YZ- DIRECTION
C                        -107-    -143-   -UM- MATRIX FOR -ZX- DIRECTION
C                        -108-    -144-   -UM- MATRIX FOR -ZY- DIRECTION
C      -49-     -73-     -109-    -145-   -UM- MATRIX FOR -ZZ- DIRECTION
C   ------------------------------------------------------------------
C
C GET DATA FROM ALPHA  AND ITERATIVE BETA CALCULATIONS
C
C
C IGAM=4 (OPTICAL KERR EFFECT)
C
      JGARC=10
      JUARC=07
      JUREC=01
      JGREC=04
C LOOP BEGINS FOR THE CALCULATION OF GAMMA(ABCD)
C
      GAV = 0.0D+00
      DO 90 IE=1,15
         IA=IDA(IE)
         IB=IDB(IE)
         IC=IDC(IE)
         ID=IDD(IE)
         ICD=IPAIR(IC,ID)
         IBD=IPAIR(IB,ID)
         IBC=IP(IB,IC)
C
C  READ IN THE FIRST ORDER U3 OMEGA AND G3 OMEGA IN THE DIRECTION A
C
C MAKE GD3 OMEGA MATRIX FROM G3 MATRIX
C
         CALL DAREAD(X,MSQ,JGARC+IA)
         CALL FHPATN(GD3,X,NORBS,2,ONE)
C
C MAKE UD3 OMEGA MATRIX FROM U3 OMEGA MATRIX
C
         CALL DAREAD(X,MSQ,JUARC+IA)
         CALL FHPATN(UD3,X,NORBS,2,-ONE)
C
         YY=0.0D00
         IMOVE=1
   30    CONTINUE
C
C OKE
         GO TO (40,50,60), IMOVE
   40    J2= IB
         J34=ICD
         JG2REC=82
         JU2REC=73
         JU2MRC=100
         JEPREC=91
         GO TO 70
   50    J2=IC
         J34=IBD
         JG2REC=82
         JU2REC=73
         JU2MRC=100
         JEPREC=91
         GO TO 70
   60    J2=ID
         J34=IBC
         JG2REC=31
         JU2REC=25
         JU2MRC=43
         JEPREC=37
   70    CONTINUE
C  READ IN G1,U1,GS,US,UMS,EPS
C
C  CALL UB
         IF (IMOVE.EQ.3) THEN
            CALL DAREAD(U1,MSQ,JUARC+J2)
         ELSE
            CALL DAREAD(U1,MSQ,JUREC+J2)
         ENDIF
C  CALL GB
         IF (IMOVE.EQ.3) THEN
            CALL DAREAD(G1,MSQ,JGARC+J2)
         ELSE
            CALL DAREAD(G1,MSQ,JGREC+J2)
         ENDIF
C  CALL GCD
         CALL DAREAD(GS,MSQ,JG2REC+J34)
C  CALL UCD
         CALL DAREAD(US,MSQ,JU2REC+J34)
C  CALL USMD
         CALL DAREAD(USMD,MSQ,JU2MRC+J34)
C  CALL EPCD
         CALL DAREAD(EPS,MSQ,JEPREC+J34)
C
C FIRST KIND
C
         YY = YY+TRSUB(UD3,G1,US,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(USMD,G1,UD3,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(UD3,G1,US,NORBS,NCLOSE,NORBS)
         YY = YY+TRSUB(USMD,G1,UD3,NORBS,NCLOSE,NORBS)
C
C SECOND KIND
C
         YY = YY+TRSUB(UD3,GS,U1,NCLOSE,NORBS,NORBS)
         YY = YY+TRSUB(U1,GS,UD3,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(UD3,EPS,U1,NORBS,NCLOSE,NORBS)
         YY = YY-TRSUB(U1,EPS,UD3,NORBS,NCLOSE,NORBS)
C
C THIRD KIND
C
         YY = YY+TRSUB(U1,GD3,US,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(USMD,GD3,U1,NCLOSE,NORBS,NORBS)
         YY = YY-TRSUB(U1,GD3,US,NORBS,NCLOSE,NORBS)
         YY = YY+TRSUB(USMD,GD3,U1,NORBS,NCLOSE,NORBS)
C
         IMOVE=IMOVE+1
         IF (IMOVE .LE. 3) GO TO 30
C
         GAMMA(IE) = YY
C
C CALCULATE THE AVERAGE GAMMA VALUE
C
         IF (IE .LE. 3) THEN
            GAV = GAV + 3.0D0*YY
         ELSEIF (IE .GT. 9) THEN
            GAV = GAV + YY
         ELSE
            GAV = GAV + 2.0D0*YY
         ENDIF
C
C WRITE GAMMA(ABCD)
C
         WRITE(6,80) ALAB(IA),ALAB(IB),ALAB(IC),ALAB(ID),GAMMA(IE)
   80    FORMAT(' GAMMA(',A1,',',A1,',',A1,',',A1,') = ',1PD14.7)
C
   90 CONTINUE
      GAVE = GAV/15.0D+00
      WRITE(6,100) OMEGA,GAVE
  100 FORMAT(//,'  AVERAGE GAMMA VALUE AT ',F10.5,' = ',1PD14.7,//)
      RETURN
      END
      SUBROUTINE NONBET (U1X,U1Y,U1Z,U2X,U2Y,U2Z,
     1                   G1X,G1Y,G1Z,G2X,G2Y,G2Z)
C
C THIS SUBROUTINE CALCULATES SECOND HARMONIC GENERATION IN A
C NONITERATIVE WAY.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
C.....................................................................
C  CALCULATE BETA IN A NONITERATIVE FASHION
C.....................................................................
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM),NORBS,NELECS,NALPHA,NBETA,NCLOSE,
     2                NOPEN,NDUMY,FRACT
      COMMON /OMVAL/ OMEGA
      COMMON /CHANEL/ IFILES(30)
      EQUIVALENCE(IW,IFILES(6))
      DIMENSION U1X(MAXORB,MAXORB),U1Y(MAXORB,MAXORB),U1Z(MAXORB,MAXORB)
      DIMENSION U2X(MAXORB,MAXORB),U2Y(MAXORB,MAXORB),U2Z(MAXORB,MAXORB)
      DIMENSION G1X(MAXORB,MAXORB),G1Y(MAXORB,MAXORB),G1Z(MAXORB,MAXORB)
      DIMENSION G2X(MAXORB,MAXORB),G2Y(MAXORB,MAXORB),G2Z(MAXORB,MAXORB)
C
C GET DATA FROM ALPHA CALCULATIONS
C
      MAXSQ = NORBS*NORBS
      BAVX = 0.0D+00
      BAVY = 0.0D+00
      BAVZ = 0.0D+00
      CALL DAREAD (U1X,MAXSQ,8)
      CALL DAREAD (U1Y,MAXSQ,9)
      CALL DAREAD (U1Z,MAXSQ,10)
      CALL DAREAD (G1X,MAXSQ,11)
      CALL DAREAD (G1Y,MAXSQ,12)
      CALL DAREAD (G1Z,MAXSQ,13)
      CALL DAREAD (U2X,MAXSQ,14)
      CALL DAREAD (U2Y,MAXSQ,15)
      CALL DAREAD (U2Z,MAXSQ,16)
      CALL DAREAD (G2X,MAXSQ,17)
      CALL DAREAD (G2Y,MAXSQ,18)
      CALL DAREAD (G2Z,MAXSQ,19)
C XXX
      CALL BETCOM (U1X,G1X,U2X,G2X,NCLOSE,NORBS,BXXX)
      BAVX = BAVX +3.0D0 * BXXX
C YXX
      CALL BETCOM (U1X,G1X,U2Y,G2Y,NCLOSE,NORBS,BYXX)
      BAVY = BAVY + BYXX
C ZXX
      CALL BETCOM (U1X,G1X,U2Z,G2Z,NCLOSE,NORBS,BZXX)
      BAVZ = BAVZ + BZXX
C XXY
      CALL BETALL (U2X,G2X,U1X,G1X,U1Y,G1Y,NCLOSE,NORBS,BXXY)
      BAVY = BAVY + BXXY
C YXY
      CALL BETALL (U2Y,G2Y,U1X,G1X,U1Y,G1Y,NCLOSE,NORBS,BYXY)
      BAVX = BAVX + BYXY
C ZXY
      CALL BETALL (U2Z,G2Z,U1X,G1X,U1Y,G1Y,NCLOSE,NORBS,BZXY)
C XXZ
      CALL BETALL (U2X,G2X,U1X,G1X,U1Z,G1Z,NCLOSE,NORBS,BXXZ)
      BAVZ = BAVZ + BXXZ
C YXZ
      CALL BETALL (U2Y,G2Y,U1X,G1X,U1Z,G1Z,NCLOSE,NORBS,BYXZ)
C ZXZ
      CALL BETALL (U2Z,G2Z,U1X,G1X,U1Z,G1Z,NCLOSE,NORBS,BZXZ)
      BAVX = BAVX + BZXZ
C XYX
      CALL BETALL (U2X,G2X,U1Y,G1Y,U1X,G1X,NCLOSE,NORBS,BXYX)
      BAVY = BAVY + BXYX
C YYX
      CALL BETALL (U2Y,G2Y,U1Y,G1Y,U1X,G1X,NCLOSE,NORBS,BYYX)
      BAVX = BAVX + BYYX
C ZYX
      CALL BETALL (U2Z,G2Z,U1Y,G1Y,U1X,G1X,NCLOSE,NORBS,BZYX)
C XYY
      CALL BETCOM (U1Y,G1Y,U2X,G2X,NCLOSE,NORBS,BXYY)
      BAVX = BAVX + BXYY
C YYY
      CALL BETCOM (U1Y,G1Y,U2Y,G2Y,NCLOSE,NORBS,BYYY)
      BAVY = BAVY + 3.0D0*BYYY
C ZYY
      CALL BETCOM (U1Y,G1Y,U2Z,G2Z,NCLOSE,NORBS,BZYY)
      BAVZ = BAVZ + BZYY
C XYZ
      CALL BETALL (U2X,G2X,U1Y,G1Y,U1Z,G1Z,NCLOSE,NORBS,BXYZ)
C YYZ
      CALL BETALL (U2Y,G2Y,U1Y,G1Y,U1Z,G1Z,NCLOSE,NORBS,BYYZ)
      BAVZ = BAVZ + BYYZ
C ZYZ
      CALL BETALL (U2Z,G2Z,U1Y,G1Y,U1Z,G1Z,NCLOSE,NORBS,BZYZ)
      BAVY = BAVY + BZYZ
C XZX
      CALL BETALL (U2X,G2X,U1Z,G1Z,U1X,G1X,NCLOSE,NORBS,BXZX)
      BAVZ = BAVZ + BXZX
C YZX
      CALL BETALL (U2Y,G2Y,U1Z,G1Z,U1X,G1X,NCLOSE,NORBS,BYZX)
C ZZX
      CALL BETALL (U2Z,G2Z,U1Z,G1Z,U1X,G1X,NCLOSE,NORBS,BZZX)
      BAVX = BAVX + BZZX
C XZY
      CALL BETALL (U2X,G2X,U1Z,G1Z,U1Y,G1Y,NCLOSE,NORBS,BXZY)
C YZY
      CALL BETALL (U2Y,G2Y,U1Z,G1Z,U1Y,G1Y,NCLOSE,NORBS,BYZY)
      BAVZ = BAVZ + BYZY
C ZZY
      CALL BETALL (U2Z,G2Z,U1Z,G1Z,U1Y,G1Y,NCLOSE,NORBS,BZZY)
      BAVY = BAVY + BZZY
C XZZ
      CALL BETCOM (U1Z,G1Z,U2X,G2X,NCLOSE,NORBS,BXZZ)
      BAVX = BAVX + BXZZ
C YZZ
      CALL BETCOM (U1Z,G1Z,U2Y,G2Y,NCLOSE,NORBS,BYZZ)
      BAVY = BAVY + BYZZ
C ZZZ
      CALL BETCOM (U1Z,G1Z,U2Z,G2Z,NCLOSE,NORBS,BZZZ)
      BAVZ = BAVZ + 3.0D0 * BZZZ
C
      BAVX = BAVX/5.0D+00
      BAVY = BAVY/5.0D+00
      BAVZ = BAVZ/5.0D+00
C
      BVEC = (BAVX*BAVX+BAVY*BAVY+BAVZ*BAVZ)**0.5D+00
      WRITE(6,10)
   10 FORMAT (//,   ' BETA (SECOND HARMONIC GENERATION)'//)
      WRITE(6,20) BXXX,BYXX,BZXX,BXXY,BYXY,BZXY,
     1              BXXZ,BYXZ,BZXZ,BXYX,BYYX,BZYX,
     2              BXYY,BYYY,BZYY,BXYZ,BYYZ,BZYZ,
     3              BXZX,BYZX,BZZX,BXZY,BYZY,BZZY,
     4              BXZZ,BYZZ,BZZZ
   20 FORMAT (//,'  BXXX  ',D15.8,'  BYXX ',D15.8,'  BZXX ',D15.8,/,
     1           '  BXXY  ',D15.8,'  BYXY ',D15.8,'  BZXY ',D15.8,/,
     2           '  BXXZ  ',D15.8,'  BYXZ ',D15.8,'  BZXZ ',D15.8,/,
     3           '  BXYX  ',D15.8,'  BYYX ',D15.8,'  BZYX ',D15.8,/,
     4           '  BXYY  ',D15.8,'  BYYY ',D15.8,'  BZYY ',D15.8,/,
     5           '  BXYZ  ',D15.8,'  BYYZ ',D15.8,'  BZYZ ',D15.8,/,
     6           '  BXZX  ',D15.8,'  BYZX ',D15.8,'  BZZX ',D15.8,/,
     7           '  BXZY  ',D15.8,'  BYZY ',D15.8,'  BZZY ',D15.8,/,
     8           '  BXZZ  ',D15.8,'  BYZZ ',D15.8,'  BZZZ ',D15.8)
C
      WRITE(6,30) OMEGA,BAVX
   30 FORMAT (//, ' AVERAGE BETAX(SHG) VALUE AT',F10.5, 'EV = ',
     1         1F15.5 )
C
      WRITE(6,40) OMEGA,BAVY
   40 FORMAT (' AVERAGE BETAY(SHG) VALUE AT',F10.5, 'EV = ',
     1         1F15.5 )
C
      WRITE(6,50) OMEGA,BAVZ
   50 FORMAT (' AVERAGE BETAZ(SHG) VALUE AT',F10.5, 'EV = ',
     1         1F15.5 ,//)
C
      WRITE(6,60) OMEGA,BVEC
   60 FORMAT (//, ' AVERAGE BETA (SHG) VALUE AT',F10.5, 'EV = ',
     1         1F15.5 ,//)
      RETURN
      END
      SUBROUTINE NONOPE (U0X,U1Y,U1Z,U1X,U0Y,U0Z,
     1                    G0X,G1Y,G1Z,G1X,G0Y,G0Z)
C
C THIS SUBROUTINE CALCULATES ELECTROOPTIC POCKEL'S EFFECT
C IN A NONITERATIVE WAY.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
C.....................................................................
C  CALCULATE BETA IN A NONITERATIVE FASHION
C.....................................................................
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM),NORBS,NELECS,NALPHA,NBETA,NCLOSE,
     2                NOPEN,NDUMY,FRACT
      COMMON /OMVAL/ OMEGA
      COMMON /CHANEL/ IFILES(30)
      EQUIVALENCE(IW,IFILES(6))
      DIMENSION U0X(MAXORB,MAXORB),U1Y(MAXORB,MAXORB),U1Z(MAXORB,MAXORB)
      DIMENSION U1X(MAXORB,MAXORB),U0Y(MAXORB,MAXORB),U0Z(MAXORB,MAXORB)
      DIMENSION G0X(MAXORB,MAXORB),G1Y(MAXORB,MAXORB),G1Z(MAXORB,MAXORB)
      DIMENSION G1X(MAXORB,MAXORB),G0Y(MAXORB,MAXORB),G0Z(MAXORB,MAXORB)
C
      MAXSQ = NORBS*NORBS
C
C READ DATA FROM ALPHA CALCULATION
C
      BAVX = 0.0D+00
      BAVY = 0.0D+00
      BAVZ = 0.0D+00
C
      CALL DAREAD (U0X,MAXSQ,2)
      CALL DAREAD (U0Y,MAXSQ,3)
      CALL DAREAD (U0Z,MAXSQ,4)
      CALL DAREAD (G0X,MAXSQ,5)
      CALL DAREAD (G0Y,MAXSQ,6)
      CALL DAREAD (G0Z,MAXSQ,7)
      CALL DAREAD (U1X,MAXSQ,8)
      CALL DAREAD (U1Y,MAXSQ,9)
      CALL DAREAD (U1Z,MAXSQ,10)
      CALL DAREAD (G1X,MAXSQ,11)
      CALL DAREAD (G1Y,MAXSQ,12)
      CALL DAREAD (G1Z,MAXSQ,13)
C XXX
      CALL BETALL (U1X,G1X,U0X,G0X,U1X,G1X,NCLOSE,NORBS,BXXX)
      BAVX = BAVX + 3.0D0 * BXXX
C YXX
      CALL BETALL (U1Y,G1Y,U0X,G0X,U1X,G1X,NCLOSE,NORBS,BYXX)
      BAVY = BAVY + BYXX
C ZXX
      CALL BETALL (U1Z,G1Z,U0X,G0X,U1X,G1X,NCLOSE,NORBS,BZXX)
      BAVZ = BAVZ + BZXX
C XXY
      CALL BETALL (U1X,G1X,U0X,G0X,U1Y,G1Y,NCLOSE,NORBS,BXXY)
      BAVY = BAVY + BXXY
C YXY
      CALL BETALL (U1Y,G1Y,U0X,G0X,U1Y,G1Y,NCLOSE,NORBS,BYXY)
      BAVX = BAVX + BYXY
C ZXY
      CALL BETALL (U1Z,G1Z,U0X,G0X,U1Y,G1Y,NCLOSE,NORBS,BZXY)
C XXZ
      CALL BETALL (U1X,G1X,U0X,G0X,U1Z,G1Z,NCLOSE,NORBS,BXXZ)
      BAVZ = BAVZ + BXXZ
C YXZ
      CALL BETALL (U1Y,G1Y,U0X,G0X,U1Z,G1Z,NCLOSE,NORBS,BYXZ)
C ZXZ
      CALL BETALL (U1Z,G1Z,U0X,G0X,U1Z,G1Z,NCLOSE,NORBS,BZXZ)
      BAVX = BAVX + BZXZ
C XYX
      CALL BETALL (U1X,G1X,U0Y,G0Y,U1X,G1X,NCLOSE,NORBS,BXYX)
      BAVY = BAVY + BXYX
C YYX
      CALL BETALL (U1Y,G1Y,U0Y,G0Y,U1X,G1X,NCLOSE,NORBS,BYYX)
      BAVX = BAVX + BYYX
C ZYX
      CALL BETALL (U1Z,G1Z,U0Y,G0Y,U1X,G1X,NCLOSE,NORBS,BZYX)
C XYY
      CALL BETALL (U1X,G1X,U0Y,G0Y,U1Y,G1Y,NCLOSE,NORBS,BXYY)
      BAVX = BAVX + BXYY
C YYY
      CALL BETALL (U1Y,G1Y,U0Y,G0Y,U1Y,G1Y,NCLOSE,NORBS,BYYY)
      BAVY = BAVY + 3.0D0 * BYYY
C ZYY
      CALL BETALL (U1Z,G1Z,U0Y,G0Y,U1Y,G1Y,NCLOSE,NORBS,BZYY)
      BAVZ = BAVZ + BZYY
C XYZ
      CALL BETALL (U1X,G1X,U0Y,G0Y,U1Z,G1Z,NCLOSE,NORBS,BXYZ)
C YYZ
      CALL BETALL (U1Y,G1Y,U0Y,G0Y,U1Z,G1Z,NCLOSE,NORBS,BYYZ)
      BAVZ = BAVZ + BYYZ
C ZYZ
      CALL BETALL (U1Z,G1Z,U0Y,G0Y,U1Z,G1Z,NCLOSE,NORBS,BZYZ)
      BAVY = BAVY + BZYZ
C XZX
      CALL BETALL (U1X,G1X,U0Z,G0Z,U1X,G1X,NCLOSE,NORBS,BXZX)
      BAVZ = BAVZ + BXZX
C YZX
      CALL BETALL (U1Y,G1Y,U0Z,G0Z,U1X,G1X,NCLOSE,NORBS,BYZX)
C ZZX
      CALL BETALL (U1Z,G1Z,U0Z,G0Z,U1X,G1X,NCLOSE,NORBS,BZZX)
      BAVX = BAVX + BZZX
C XZY
      CALL BETALL (U1X,G1X,U0Z,G0Z,U1Y,G1Y,NCLOSE,NORBS,BXZY)
C YZY
      CALL BETALL (U1Y,G1Y,U0Z,G0Z,U1Y,G1Y,NCLOSE,NORBS,BYZY)
      BAVZ = BAVZ + BYZY
C ZZY
      CALL BETALL (U1Z,G1Z,U0Z,G0Z,U1Y,G1Y,NCLOSE,NORBS,BZZY)
      BAVY = BAVY + BZZY
C XZZ
      CALL BETALL (U1X,G1X,U0Z,G0Z,U1Z,G1Z,NCLOSE,NORBS,BXZZ)
      BAVX = BAVX + BXZZ
C YZZ
      CALL BETALL (U1Y,G1Y,U0Z,G0Z,U1Z,G1Z,NCLOSE,NORBS,BYZZ)
      BAVY = BAVY + BYZZ
C ZZZ
      CALL BETALL (U1Z,G1Z,U0Z,G0Z,U1Z,G1Z,NCLOSE,NORBS,BZZZ)
      BAVZ = BAVZ + 3.0D0 * BZZZ
C
      BAVX = BAVX/5.0D+00
      BAVY = BAVY/5.0D+00
      BAVZ = BAVZ/5.0D+00
C
      BVEC = (BAVX*BAVX+BAVY*BAVY+BAVZ*BAVZ)**0.5D+00
      WRITE(6,*)   '  BETA (ELECTOPTIC POCKELS EFFECT) '
      WRITE(6,10) BXXX,BYXX,BZXX,BXXY,BYXY,BZXY,
     1              BXXZ,BYXZ,BZXZ,BXYX,BYYX,BZYX,
     2              BXYY,BYYY,BZYY,BXYZ,BYYZ,BZYZ,
     3              BXZX,BYZX,BZZX,BXZY,BYZY,BZZY,
     4              BXZZ,BYZZ,BZZZ
   10 FORMAT (//,'  BXXX  ',D15.8,'  BYXX ',D15.8,'  BZXX ',D15.8,/,
     1           '  BXXY  ',D15.8,'  BYXY ',D15.8,'  BZXY ',D15.8,/,
     2           '  BXXZ  ',D15.8,'  BYXZ ',D15.8,'  BZXZ ',D15.8,/,
     3           '  BXYX  ',D15.8,'  BYYX ',D15.8,'  BZYX ',D15.8,/,
     4           '  BXYY  ',D15.8,'  BYYY ',D15.8,'  BZYY ',D15.8,/,
     5           '  BXYZ  ',D15.8,'  BYYZ ',D15.8,'  BZYZ ',D15.8,/,
     6           '  BXZX  ',D15.8,'  BYZX ',D15.8,'  BZZX ',D15.8,/,
     7           '  BXZY  ',D15.8,'  BYZY ',D15.8,'  BZZY ',D15.8,/,
     8           '  BXZZ  ',D15.8,'  BYZZ ',D15.8,'  BZZZ ',D15.8)
C
      WRITE(6,20) OMEGA,BAVX
   20 FORMAT(//,' AVERAGE BETAX VALUE AT', F10.5, 'EV = ',
     1         1F15.5 )
C
      WRITE(6,30) OMEGA,BAVY
   30 FORMAT(' AVERAGE BETAY VALUE AT', F10.5, 'EV = ',
     1         1F15.5 )
C
      WRITE(6,40) OMEGA,BAVZ
   40 FORMAT(' AVERAGE BETAZ VALUE AT', F10.5, 'EV = ',
     1         1F15.5 ,//)
C
      WRITE(6,50) OMEGA,BVEC
   50 FORMAT(//,' AVERAGE BETA(EOPE) VALUE AT', F10.5, 'EV = ',
     1         1F15.5 ,//)
      RETURN
      END
      SUBROUTINE NONOR (U0X,U1Y,U1Z,U1X,U0Y,U0Z,G0X,G1Y,G1Z,
     1                  G1X,G0Y,G0Z)
C
C THIS SUBROUTINE CALCULATES OPTICAL RECTIFICATION IN A
C NONITERATIVE WAY
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
C.....................................................................
C  CALCULATE BETA IN A NONITERATIVE FASHION
C.....................................................................
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM),NORBS,NELECS,NALPHA,NBETA,NCLOSE,
     2                NOPEN,NDUMY,FRACT
      COMMON /OMVAL/ OMEGA
      COMMON /CHANEL/ IFILES(30)
      EQUIVALENCE(IW,IFILES(6))
      DIMENSION U0X(MAXORB,MAXORB),U1Y(MAXORB,MAXORB),U1Z(MAXORB,MAXORB)
      DIMENSION U1X(MAXORB,MAXORB),U0Y(MAXORB,MAXORB),U0Z(MAXORB,MAXORB)
      DIMENSION G0X(MAXORB,MAXORB),G1Y(MAXORB,MAXORB),G1Z(MAXORB,MAXORB)
      DIMENSION G1X(MAXORB,MAXORB),G0Y(MAXORB,MAXORB),G0Z(MAXORB,MAXORB)
      MAXSQ = NORBS*NORBS
C
C READ DATA FROM ALPHA CALCULATION
C
      BAVX = 0.0D+00
      BAVY = 0.0D+00
      BAVZ = 0.0D+00
C
      CALL DAREAD (U0X,MAXSQ,2)
      CALL DAREAD (U0Y,MAXSQ,3)
      CALL DAREAD (U0Z,MAXSQ,4)
      CALL DAREAD (G0X,MAXSQ,5)
      CALL DAREAD (G0Y,MAXSQ,6)
      CALL DAREAD (G0Z,MAXSQ,7)
      CALL DAREAD (U1X,MAXSQ,8)
      CALL DAREAD (U1Y,MAXSQ,9)
      CALL DAREAD (U1Z,MAXSQ,10)
      CALL DAREAD (G1X,MAXSQ,11)
      CALL DAREAD (G1Y,MAXSQ,12)
      CALL DAREAD (G1Z,MAXSQ,13)
C
C NONITERATIVE BETA CALCULATION
C
C XXX
      CALL BETAL1 (U0X,G0X,U1X,G1X,U1X,G1X,NCLOSE,NORBS,BXXX)
      BAVX = BAVX + 3.0D0* BXXX
C YXX
      CALL BETAL1 (U0Y,G0Y,U1X,G1X,U1X,G1X,NCLOSE,NORBS,BYXX)
      BAVY = BAVY + BYXX
C ZXX
      CALL BETAL1 (U0Z,G0Z,U1X,G1X,U1X,G1X,NCLOSE,NORBS,BZXX)
      BAVZ = BAVZ + BZXX
C XXY
      CALL BETAL1 (U0X,G0X,U1X,G1X,U1Y,G1Y,NCLOSE,NORBS,BXXY)
      BAVY = BAVY + BXXY
C YXY
      CALL BETAL1 (U0Y,G0Y,U1X,G1X,U1Y,G1Y,NCLOSE,NORBS,BYXY)
      BAVX = BAVX + BYXY
C ZXY
      CALL BETAL1 (U0Z,G0Z,U1X,G1X,U1Y,G1Y,NCLOSE,NORBS,BZXY)
C XXZ
      CALL BETAL1 (U0X,G0X,U1X,G1X,U1Z,G1Z,NCLOSE,NORBS,BXXZ)
      BAVZ = BAVZ + BXXZ
C YXZ
      CALL BETAL1 (U0Y,G0Y,U1X,G1X,U1Z,G1Z,NCLOSE,NORBS,BYXZ)
C ZXZ
      CALL BETAL1 (U0Z,G0Z,U1X,G1X,U1Z,G1Z,NCLOSE,NORBS,BZXZ)
      BAVX = BAVX + BZXZ
C XYX
      CALL BETAL1 (U0X,G0X,U1Y,G1Y,U1X,G1X,NCLOSE,NORBS,BXYX)
      BAVY = BAVY + BXYX
C YYX
      CALL BETAL1 (U0Y,G0Y,U1Y,G1Y,U1X,G1X,NCLOSE,NORBS,BYYX)
      BAVX = BAVX + BYYX
C ZYX
      CALL BETAL1 (U0Z,G0Z,U1Y,G1Y,U1X,G1X,NCLOSE,NORBS,BZYX)
C XYY
      CALL BETAL1 (U0X,G0X,U1Y,G1Y,U1Y,G1Y,NCLOSE,NORBS,BXYY)
      BAVX = BAVX + BXYY
C YYY
      CALL BETAL1 (U0Y,G0Y,U1Y,G1Y,U1Y,G1Y,NCLOSE,NORBS,BYYY)
      BAVY = BAVY + 3.0D0 * BYYY
C ZYY
      CALL BETAL1 (U0Z,G0Z,U1Y,G1Y,U1Y,G1Y,NCLOSE,NORBS,BZYY)
      BAVZ = BAVZ + BZYY
C XYZ
      CALL BETAL1 (U0X,G0X,U1Y,G1Y,U1Z,G1Z,NCLOSE,NORBS,BXYZ)
C YYZ
      CALL BETAL1 (U0Y,G0Y,U1Y,G1Y,U1Z,G1Z,NCLOSE,NORBS,BYYZ)
      BAVZ = BAVZ + BYYZ
C ZYZ
      CALL BETAL1 (U0Z,G0Z,U1Y,G1Y,U1Z,G1Z,NCLOSE,NORBS,BZYZ)
      BAVY = BAVY + BZYZ
C XZX
      CALL BETAL1 (U0X,G0X,U1Z,G1Z,U1X,G1X,NCLOSE,NORBS,BXZX)
      BAVZ = BAVZ + BXZX
C YZX
      CALL BETAL1 (U0Y,G0Y,U1Z,G1Z,U1X,G1X,NCLOSE,NORBS,BYZX)
C ZZX
      CALL BETAL1 (U0Z,G0Z,U1Z,G1Z,U1X,G1X,NCLOSE,NORBS,BZZX)
      BAVX = BAVX + BZZX
C XZY
      CALL BETAL1 (U0X,G0X,U1Z,G1Z,U1Y,G1Y,NCLOSE,NORBS,BXZY)
C YZY
      CALL BETAL1 (U0Y,G0Y,U1Z,G1Z,U1Y,G1Y,NCLOSE,NORBS,BYZY)
      BAVZ = BAVZ + BYZY
C ZZY
      CALL BETAL1 (U0Z,G0Z,U1Z,G1Z,U1Y,G1Y,NCLOSE,NORBS,BZZY)
      BAVY = BAVY + BZZY
C XZZ
      CALL BETAL1 (U0X,G0X,U1Z,G1Z,U1Z,G1Z,NCLOSE,NORBS,BXZZ)
      BAVX = BAVX + BXZZ
C YZZ
      CALL BETAL1 (U0Y,G0Y,U1Z,G1Z,U1Z,G1Z,NCLOSE,NORBS,BYZZ)
      BAVY = BAVY + BYZZ
C ZZZ
      CALL BETAL1 (U0Z,G0Z,U1Z,G1Z,U1Z,G1Z,NCLOSE,NORBS,BZZZ)
      BAVZ = BAVZ + 3.0D0 * BZZZ
C
      BAVX = BAVX/5.0D+00
      BAVY = BAVY/5.0D+00
      BAVZ = BAVZ/5.0D+00
C
      BVEC = (BAVX*BAVX+BAVY*BAVY+BAVZ*BAVZ)**0.5D+00
      WRITE(6,10)
   10 FORMAT  (//, ' BETA (OPTICAL RECTIFICATION) ')
      WRITE(6,20) BXXX,BYXX,BZXX,BXXY,BYXY,BZXY,
     1              BXXZ,BYXZ,BZXZ,BXYX,BYYX,BZYX,
     2              BXYY,BYYY,BZYY,BXYZ,BYYZ,BZYZ,
     3              BXZX,BYZX,BZZX,BXZY,BYZY,BZZY,
     4              BXZZ,BYZZ,BZZZ
   20 FORMAT (//,'  BXXX  ',D15.8,'  BYXX ',D15.8,'  BZXX ',D15.8,/,
     1           '  BXXY  ',D15.8,'  BYXY ',D15.8,'  BZXY ',D15.8,/,
     2           '  BXXZ  ',D15.8,'  BYXZ ',D15.8,'  BZXZ ',D15.8,/,
     3           '  BXYX  ',D15.8,'  BYYX ',D15.8,'  BZYX ',D15.8,/,
     4           '  BXYY  ',D15.8,'  BYYY ',D15.8,'  BZYY ',D15.8,/,
     5           '  BXYZ  ',D15.8,'  BYYZ ',D15.8,'  BZYZ ',D15.8,/,
     6           '  BXZX  ',D15.8,'  BYZX ',D15.8,'  BZZX ',D15.8,/,
     7           '  BXZY  ',D15.8,'  BYZY ',D15.8,'  BZZY ',D15.8,/,
     8           '  BXZZ  ',D15.8,'  BYZZ ',D15.8,'  BZZZ ',D15.8)
C
      WRITE(6,30) OMEGA,BAVX
   30 FORMAT(//,' AVERAGE BETAX VALUE AT ',F10.5, 'EV = ',
     1        1F15.5 )
C
      WRITE(6,40) OMEGA,BAVY
   40 FORMAT(' AVERAGE BETAY VALUE AT ',F10.5, 'EV = ',
     1        1F15.5 )
C
      WRITE(6,50) OMEGA,BAVZ
   50 FORMAT(' AVERAGE BETAZ VALUE AT ',F10.5, 'EV = ',
     1        1F15.5 ,//)
C
      WRITE(6,60) OMEGA,BVEC
   60 FORMAT(//,' AVERAGE BETA(OR) VALUE AT ',F10.5, 'EV = ',
     1        1F15.5 ,//)
C
      RETURN
      END
      SUBROUTINE OPENDA(IREST)
C
C     - - - - OPEN MASTER DICTIONARY FILE 10 - - - -
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C
      COMMON /IODAF/ IDAF,IRECLN,IRECST,IFILEN(145),IODA(145)
      COMMON /JOBNAM/ JOBNAM
      COMMON /CHANEL/ IFILES(30)
      CHARACTER*80 JOBNAM
      EQUIVALENCE(IW,IFILES(6))
C
      IDAF = 17
*VAX  IRECLN = 512
      IRECLN = 1023
C
C             GET OPEN PARAMETERS FROM INPUT (OPTIONAL)
C
C           OBTAIN DICTNRY FILE NAME FROM RUN COMMAND LINE
C
*VAX      OPEN (UNIT=IDAF, FILE='DICTNRY', STATUS='UNKNOWN',
*VAX     *      ACCESS='DIRECT', FORM='UNFORMATTED', RECL=2*IRECLN)
C
C---------- modified by I. Cserny,  June 21, 1995 -------------
C-    I=INDEX(JOBNAM,' ')-1
C-    OPEN(UNIT=IDAF, FILE=JOBNAM(:I)//'.POL', STATUS='UNKNOWN',
C-   1      ACCESS='DIRECT', FORM='UNFORMATTED', RECL=8*IRECLN)
      OPEN(UNIT=IDAF, FILE='DICTNRY', STATUS='UNKNOWN',
     1      ACCESS='DIRECT', FORM='UNFORMATTED', RECL=8*IRECLN)
C--------------------------------------------------------------
C
C     ----- IS THIS A NEW OR OLD DAF FILE -----
C
      IF (IREST .NE. 0) GO TO 20
C
C        ----- MARK THE NEW DAF RECORDS AS EMPTY -----
C
      IRECST = 1
      DO 10 I = 1,145
         IODA(I) = -1
   10 CONTINUE
      IRECST = IRECST + 1
      WRITE(UNIT=IDAF, REC=1) IRECST,IODA,IFILEN,IS,IPK
      RETURN
C
C     ----- LOAD THE OLD DAF DIRECTORY -----
C
   20 CONTINUE
      READ(UNIT=IDAF, REC=1) IRECST,IODA,IFILEN,IS,IPK
      RETURN
      END