C  2 FEB 06 - MWS - ONEFLD OPTION FIXED
C 14 NOV 05 - DGF - PAD COMMON BLOCK ENRGYS
C 19 SEP 05 - SG  - DO ALL 3 FIELD DIRECTIONS TO PRODUCE FULL TENSORS
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 19 MAY 04 - JAB - PRINT A MESSAGE ONLY ONCE IN PARALLEL JOBS
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C 26 MAR 03 - MWS - TAKE ENERGY FROM FUNCT COMMON BLOCK, NOT ENRGYS
C  6 SEP 01 - MWS - ADD DUMMY ARGUMENTS TO NAMEIO CALL
C 25 JUN 01 - MWS - ALTER COMMON BLOCK WFNOPT
C 13 FEB 99 - MAF - FFLDX: CORRECT SYMMETRY OFF/ON FOR SPHER HARMONICS
C 13 JUN 96 - VAG - ADD VARIABLE FOR CI TYPE TO SCFOPT COMMON
C  5 MAR 96 - MWS - CHANGE SHELL SYMMETRY PACKING COMMONS
C  3 JAN 96 - MWS - FFLDX: CALL NEW ONE ELECTRON INTEGRAL DRIVER
C 28 SEP 95 - MWS - FFLDX: USE DIPINT ROUTINE TO GET THESE INTEGRALS
C 10 NOV 94 - MWS - FFLDX: REMOVE FTNCHEK WARNINGS
C 27 OCT 94 - MWS - FFLDX: DELETE OLD MP2PAR COMMON
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C  2 MAY 94 - MWS - FFLDX: FIX CHECK RUNS
C 16 JUL 93 - MWS - INCREASE MAXIMUM CI ROOTS TO 100
C 31 OCT 92 - MWS - PREVENT SEMIEMPIRICAL FIELDS (NO DIPOLE INTEGRALS)
C 22 JUN 92 - MWS - ADD 2ND ARG TO GUESMO CALL
C  1 MAY 92 - MWS - INCLUDE HENRY KURTZ' FFIELD CALCULATION IN GAMESS
C
C*MODULE FFIELD  *DECK FFLDX
      SUBROUTINE FFLDX
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL ONEFLD,OFFDIA,GOPARR,DSKWRK,MASWRK,
     *        MCCI,SYM,PK,PKSAVE,PANDK,BLOCK,EFLDL
C
      PARAMETER (MXATM=2000, MXSH=5000)
C
      DIMENSION EFIELD(5),EFVAL(3),FFENER(13),KFLI(13),KFLJ(13),
     *          DIP(13,3),EONE(3),COM(3),ALPHAE(3,3),BETAE(3,3),
     *          GAMMAE(3,3),ALPHAD(3,3),BETAD(3,3),GAMMAD(3,3),
     *          FDIPOL(3),DDIPOL(3),ENSAV(3,2:3),DPSAV(3,2:3,3)
C
      COMMON /EFLDC / EVEC(3),EFLDL
      COMMON /ENRGMP/ EMP2,EMP3,EMP4
      COMMON /FMCOM / X(1)
      COMMON /FUNCT / E,EG(3,MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MP2PAR/ OSPT,TOL,METHOD,NWDMP2,MPPROP,
     *                NACORE,NBCORE,NOA,NOB,NORB,NBF,NOMIT,MOCPHF,MAXITC
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROCS,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PKFIL / PK,PANDK,BLOCK
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /TMVALS/ TI,TX,TIM
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      COMMON /XYZPRP/ XP,YP,ZP,
     *                DIPMX,DIPMY,DIPMZ,
     *                QXX,QYY,QZZ,QXY,QXZ,QYZ,
     *                QMXX,QMYY,QMZZ,QMXY,QMXZ,QMYZ,
     *                OXXX,OXXY,OXXZ,OXYY,OYYY,OYYZ,
     *                OXZZ,OYZZ,OZZZ,OXYZ,
     *                OMXXX,OMXXY,OMXXZ,OMXYY,OMYYY,
     *                OMYYZ,OMXZZ,OMYZZ,OMZZZ,OMXYZ
C
      PARAMETER (ZERO=0.0D+00, TWO=2.0D+00, DFAC=2.541766D+00,
     *           ZAPTOL=1.0D-04)
C
      PARAMETER (NNAM=5)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA FFCALC/8HFFCALC  /
      DATA QNAM/8HEFIELD  ,8HONEFLD  ,8HSYM     ,
     *          8HESTEP   ,8HOFFDIA  /
      DATA KQNAM/33,0,0,  3,0/
C
      DATA KFLI/1,3,3,1,2,2,2,1,3,5,5,4,4/
      DATA KFLJ/1,1,2,2,2,1,3,3,3,5,4,4,5/
C
      DATA RMC/8HMCSCF    /
      DATA CHECK/8HCHECK   /
      DATA RNONE/8HNONE    /
C
C        CALCULATION OF POLARIZABILITY QUANTITIES BY APPLIED FIELDS,
C        STATIC ALPHA, BETA, AND GAMMA COMPUTED BY FINITE DIFFERENCING
C
C        THE METHOD IS DESCRIBED IN THE PAPER
C       "CALCULATION OF NONLINEAR OPTICAL PROPERTIES OF MOLECULES"
C        H.A.KURTZ, J.J.P.STEWART, K.M.DIETER
C        J.COMPUT.CHEM.  11, 82-87 (1990)
C
C        THIS CODE CONTRIBUTED BY HENRY KURTZ TO GAMESS IN APRIL 1992
C        MODIFIED BY SANKHA GHOSH IN SEPT. 2005 TO COMPUTE FULL TENSORS
C
      IF(MASWRK) WRITE(IW,9000)
      MCCI = SCFTYP.EQ.RMC  .OR.  CITYP.NE.RNONE
C
      ESTEP = 1.0D-03
      ONEFLD = .FALSE.
      SYM    = .FALSE.
      OFFDIA = .TRUE.
      EONE(1) = ZERO
      EONE(2) = ZERO
      EONE(3) = ZERO
C
      JRET = 0
      CALL NAMEIO (IR,JRET,FFCALC,NNAM,QNAM,KQNAM,
     *             EONE,ONEFLD,SYM,ESTEP,OFFDIA,
     *             0,0,0,0,
     *   0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     *   0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     *   0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0)
      IF(JRET.GT.1) THEN
         IF(MASWRK) WRITE(IW,*) 'ERROR IN $FFCALC GROUP'
         CALL ABRT
      END IF
C
      IF(MASWRK) THEN
         WRITE(IW,9010)
         IF(ONEFLD) THEN
            WRITE(IW,9015) ONEFLD,SYM,(EONE(IE),IE=1,3)
         ELSE
            WRITE(IW,9016) OFFDIA,ESTEP
         END IF
      END IF
C
C     *******************************************
C     MP2 SHOULD GO FASTEST WITH JK INTEGRAL FILE
C     *******************************************
C
      PKSAVE = PK
      NOPKSV = NOPK
      IF(MPLEVL.EQ.2) PK=.FALSE.
      IF(MPLEVL.EQ.2) NOPK=1
C
C     --- ENERGY COMPUTATION USING ONLY ONE APPLIED FIELD ---
C     NOTE THAT DIPOLE INTEGRALS SHOULD PRE-EXIST FROM ABOVE.
C     RESET SYMMETRY TO C1 FOR ASYMMETRIC APPLIED FIELDS.
C
      IF(ONEFLD) THEN
         EFVAL(1) = EONE(1)
         EFVAL(2) = EONE(2)
         EFVAL(3) = EONE(3)
         MFLD=1
         IF(MASWRK) WRITE(IW,9090) MFLD,(EFVAL(J),J=1,3)
C
C                 SET UP FOR A TYPICAL $EFIELD COMPUTATION
C           EFLDL CAUSES DIPOLE INTS TO BE ADDED TO H INSIDE -ONEEI-,
C           IT WILL ALSO LEAD TO CORRECTING THE NUCLEAR REP. ENERGY.
C
         EFLDL = .TRUE.
         EVEC(1) = EONE(1)
         EVEC(2) = EONE(2)
         EVEC(3) = EONE(3)
C
         IF(.NOT.SYM) CALL SYMOFF
         CALL ONEEI
         CALL GUESMO(GUESS)
         CALL JANDK
         CALL WFN
         GO TO 888
      END IF
C
C     ***************************************
C     CALCULATE FIELD-FREE ENERGY OF MOLECULE
C     ***************************************
C
      IF(MASWRK) WRITE(IW,9020)
      CALL ENERGX
C
      IF(E.EQ.ZERO  .AND.  EXETYP.NE.CHECK) THEN
         IF(MASWRK) WRITE(IW,9030)
         CALL ABRT
         STOP
      END IF
      IF(TIM .GE. TIMLIM) THEN
         IF(MASWRK) WRITE(IW,9040)
         CALL ABRT
         STOP
      END IF
C
C     ***************************************************
C     EVALUATE AND SAVE DIPOLE INTEGRALS IN C.O.M. COORDS
C     ***************************************************
C
      CALL CALCOM(XP,YP,ZP)
      COM(1) = XP
      COM(2) = YP
      COM(3) = ZP
      CALL DIPINT(COM(1),COM(2),COM(3),.FALSE.)
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
C
      CALL VALFM(LOADFM)
      LHAM = LOADFM + 1
      LVEC = LHAM   + L2
      LAST = LVEC   + L3
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C     ***********************
C     SAVE FIELD FREE RESULTS
C     ***********************
      IF(MPLEVL .GT. 0) THEN
        FFENER(1) = EMP2
      ELSE
        FFENER(1) = E
      END IF
C
      IF(MPLEVL .GT. 0 .AND. MPPROP .EQ. 1) THEN
         DIP(1,1) = DIPMX/DFAC
         DIP(1,2) = DIPMY/DFAC
         DIP(1,3) = DIPMZ/DFAC
      ELSE IF(CCTYP .NE. RNONE .OR.
     *        (MPLEVL .GT. 0 .AND. MPPROP .EQ. 0))THEN
         DIP(1,1) = ZERO
         DIP(1,2) = ZERO
         DIP(1,3) = ZERO
      END IF
C
      IF(MASWRK) THEN
         WRITE(IW,9050) 1,ZERO,ZERO,ZERO
         WRITE(IW,9060) FFENER(1)
         WRITE(IW,9070) DIPMX,DIPMY,DIPMZ
      END IF
C     ************************************************
C     SAVE FIELD FREE ORBITALS, AND POSSIBLY CI VECTOR
C     ************************************************
      CALL FFSVMO(0,X(LVEC))
      IF(MCCI) CALL FFSVCI(0)
C     **********************************************
C     SAVE ORIGINAL ONE ELECTRON INTEGRALS FOR LATER
C     **********************************************
      CALL DAREAD(IDAF,IODA,X(LHAM),L2,11,0)
      CALL DAWRIT(IDAF,IODA,X(LHAM),L2,94,0)
      CALL RETFM(NEED)
C
C     *************************************************************
C     *           -----  SET UP THE BASE FIELDS -----             *
C     *  IF AXIAL COMPONENTS ONLY (JAXIS.EQ.0) THEN FIELDS ARE    *
C     *  1:    0                                                  *
C     *  2:    I      4:  2I                                      *
C     *  3:   -I      5: -2I                                      *
C     *                                                           *
C     *  IF OFF AXIS (JAXIS.NE.0) THEN GO AROUND THE PLANE IN     *
C     *  AS SYSTEMATIC WAY AS POSSIBLE TO HELP CONVERGENCE.       *
C     *  1:    0    0                                             *
C     *  2:   -I    0        6:    I    0          10:  -2I  -2J  *
C     *  3:   -I    J        7:    I   -J          11:  -2I   2J  *
C     *  4:    0    J        8:    O   -J          12:   2I   2J  *
C     *  5:    I    J        9:   -I   -J          13:   2I  -2J  *
C     *************************************************************
C
      EFIELD(1) =  ZERO
      EFIELD(2) =  ESTEP
      EFIELD(3) = -ESTEP
      EFIELD(4) =  TWO*ESTEP
      EFIELD(5) = -TWO*ESTEP
      NFLDS     =  5
      SYM       = .FALSE.
C
C     ----- RESET SYMMETRY TO C1 FOR ASYMMETRIC APPLIED FIELDS -----
C     IF NECESSARY, CALCULATE A GRANDE LIST OF TWO ELECTRON INTEGRALS
C     AND ALSO TURN OFF USE OF SYMMETRY IN FOCK DIAGONALIZATIONS.
C     USING SYMOFF FOR C1 LOOKS STUPID, BUT SETS UP SYMON CALL LATER.
C
      IF(NT.GT.1) THEN
         IF(MASWRK) WRITE(IW,9080)
         CALL SYMOFF
         CALL JANDK
      ELSE
         CALL SYMOFF
      END IF
C
      IF(TIM .GE. TIMLIM) THEN
         IF(MASWRK) WRITE(IW,9040)
         CALL ABRT
         STOP
      END IF
C
      IF(EXETYP.EQ.CHECK) THEN
         CALL VCLR(FFENER,1,13  )
         CALL VCLR(DIP   ,1,13*3)
      END IF
C
C     ***********************
C     INITIALIZE THE MATRICES
C     ***********************
C
      CALL VCLR(FDIPOL,1,3)
      CALL VCLR(DDIPOL,1,3)
      CALL VCLR(ALPHAE,1,9)
      CALL VCLR(ALPHAD,1,9)
      CALL VCLR( BETAE,1,9)
      CALL VCLR( BETAD,1,9)
      CALL VCLR(GAMMAE,1,9)
      CALL VCLR(GAMMAD,1,9)
C
C     ******************
C     LOOP OVER THE AXIS
C     ******************
      DO IAXIS = 1,3
C
C        ****************************
C        LOOP OVER ALL APPLIED FIELDS
C        ****************************
            NPRINT = -5
            DO 400 IFLD = 2,NFLDS
               FFENER(IFLD) = ZERO
               DIP(IFLD,1)  = ZERO
               DIP(IFLD,2)  = ZERO
               DIP(IFLD,3)  = ZERO
               EFVAL(1)     = ZERO
               EFVAL(2)     = ZERO
               EFVAL(3)     = ZERO
C
C            SETUP THIS ELECTRIC FIELD
C
               EFVAL(IAXIS) = EFIELD(IFLD)
C
               MFLD = IFLD-1
               IF(MASWRK) WRITE(IW,9090) MFLD,(EFVAL(J),J=1,3)
               IF(EXETYP.EQ.CHECK) GO TO 400
C           **********************************************
C           ADD ELECTRIC FIELD TO ONE ELECTRON HAMILTONIAN
C           **********************************************
               CALL VALFM (LOADFM)
               LMU  = LOADFM + 1
               LH   = LMU    + L2
               LAST = LH     + L2
               NEED = LAST - LOADFM - 1
               CALL GETFM(NEED)
C
               CALL DAREAD(IDAF,IODA,X(LH),L2,94,0)
C
               IF(EFVAL(1).NE.ZERO) THEN
                 CALL DAREAD(IDAF,IODA,X(LMU),L2,95,0)
                 CALL DAXPY(L2,EFVAL(1),X(LMU),1,X(LH),1)
               END IF
               IF(EFVAL(2).NE.ZERO) THEN
                 CALL DAREAD(IDAF,IODA,X(LMU),L2,96,0)
                 CALL DAXPY(L2,EFVAL(2),X(LMU),1,X(LH),1)
               END IF
               IF(EFVAL(3).NE.ZERO) THEN
                 CALL DAREAD(IDAF,IODA,X(LMU),L2,97,0)
                 CALL DAXPY(L2,EFVAL(3),X(LMU),1,X(LH),1)
               END IF
C
               CALL DAWRIT(IDAF,IODA,X(LH),L2,11,0)
C           *******************************************
C           RECOVER THE FIELD FREE ORBITALS AND DENSITY
C           *CAUTION, STORAGE AT -LMU- IS USED AS L3!*
C           *******************************************
               CALL FFSVMO(1,X(LMU))
               IF(MCCI) CALL FFSVCI(1)
               CALL ORTHDN
               CALL RETFM(NEED)
C           ***************************************
C           CALCULATE WAVEFUNCTION UNDER THIS FIELD
C           ***************************************
               CALL WFN
C
               IF(E.EQ.ZERO) THEN
               IF(MASWRK) WRITE(IW,9030)
                 CALL ABRT
                 STOP
               END IF
C
               IF(TIM .GE. TIMLIM) THEN
               IF(MASWRK) WRITE(IW,9040)
                 CALL ABRT
                 STOP
               END IF
C           ***********************************
C           ADD NUCLEAR CONTRIBUTION FROM FIELD
C           ***********************************
               ANUCF = ZERO
               DO JNUF = 1,NAT
                  DO INUF = 1,3
                     COMCOR = C(INUF,JNUF) - COM(INUF)
                     ANUCF = ANUCF - EFVAL(INUF)*COMCOR*ZAN(JNUF)
                  END DO
               END DO
C
               EWITHF = E + ANUCF
               FFENER(IFLD) = EWITHF
C           ****************
C           CALCULATE DIPOLE
C           ****************
               IF(CCTYP .NE. RNONE .OR.
     *           (MPLEVL .GT. 0 .AND. MPPROP .EQ. 0))THEN
                 DIP(IFLD,1) = ZERO
                 DIP(IFLD,2) = ZERO
                 DIP(IFLD,3) = ZERO
               ELSE IF(MPLEVL .GT. 0 .AND. MPPROP .EQ. 1)THEN
                 CALL ELMOMC
                 DIP(IFLD,1) = DIPMX/DFAC
                 DIP(IFLD,2) = DIPMY/DFAC
                 DIP(IFLD,3) = DIPMZ/DFAC
               ELSE
                 CALL ELMOMC
                 DIP(IFLD,1) = DIPMX/DFAC
                 DIP(IFLD,2) = DIPMY/DFAC
                 DIP(IFLD,3) = DIPMZ/DFAC
               END IF
C           ************************************
C           SAVE THE RESULTS FOR FILDS : I, & -I
C           ************************************
               IF(IFLD .EQ. 2 .OR. IFLD .EQ. 3) THEN
                 ENSAV(IAXIS,IFLD)   = FFENER(IFLD)
                 DPSAV(IAXIS,IFLD,1) = DIP(IFLD,1)
                 DPSAV(IAXIS,IFLD,2) = DIP(IFLD,2)
                 DPSAV(IAXIS,IFLD,3) = DIP(IFLD,3)
               END IF
C
               IF(MASWRK) THEN
                 MFLD = IFLD-1
                 WRITE(IW,9050) MFLD,(EFVAL(J),J=1,3)
                 WRITE(IW,9060) EWITHF
                 WRITE(IW,9070) DIPMX,DIPMY,DIPMZ
               END IF
  400       CONTINUE
C
C              NOTES ADDED BY MWS, JAN 2006.
C    THE DIFFERENTIATIONS ARE MORE THAN A BIT DODGY, FOR
C    A SYMMETRIC MOLECULE LIKE WATER WE END UP WITH SOME
C    OF THE BB BELOW BEING 1D-13.  THERE PROBABLY SHOULD
C    BE A TEST ZAPPING ALL THE NUMERATORS TO ZERO IF THEY
C    ARE VERY SMALL COMPARED TO THE ESTEP POWER BEING
C    USED IN THEIR DENOMINATOR.  SHOULD ESTEP BE BIGGER?
C
C    I'D GUESS THAT THERE ARE ONLY 3 SIGNIFICANT FIGURES
C    IN THE BETA AND GAMMA TENSOR COMPONENTS.  PRINT LESS?
C
C    ALSO, SHOULDN'T THERE BE GAMMA COMPONENTS LIKE XXXY?
C
               BD = -(2.0D+00/3.0D+00) *(FFENER(2)-FFENER(3))
     *              + (1.0D+00/12.0D+00)*(FFENER(4)-FFENER(5))
               FDIPOL(IAXIS) = BD/ESTEP
               BD = (2.0D+00/3.0D+00)*(DIP(2,IAXIS) + DIP(3,IAXIS))
     *             - (1.0D+00/6.0D+00)*(DIP(4,IAXIS) + DIP(5,IAXIS))
               DDIPOL(IAXIS) = BD
               BA =  (5.0D+00/2.0D+00) * FFENER(1)
     *              - (4.0D+00/3.0D+00) *(FFENER(2) + FFENER(3))
     *              + (1.0D+00/12.0D+00)*(FFENER(4) + FFENER(5))
               ALPHAE(IAXIS,IAXIS) = BA/(ESTEP*ESTEP)
               DO IB = 1,3
                  BA = (2.0D+00/3.0D+00) *(DIP(2,IB) - DIP(3,IB))
     *                 - (1.0D+00/12.0D+00)*(DIP(4,IB) - DIP(5,IB))
                  ALPHAD(IAXIS,IB) = BA/ESTEP
               END DO
               BB = (FFENER(2)-FFENER(3)) -
     *              (1.0D+00/2.0D+00)*(FFENER(4)-FFENER(5))
               BETAE(IAXIS,IAXIS) = BB/(ESTEP**3)
               DO IB = 1,3
                  BB = (1.0D+00/3.0D+00)*(DIP(4,IB) + DIP(5,IB) -
     *                                    DIP(2,IB) - DIP(3,IB))
                  BETAD(IAXIS,IB) = BB/(ESTEP*ESTEP)
               END DO
               BG = -6.0D+00*FFENER(1) + 4.0D+00*(FFENER(2)+FFENER(3))
     *              - (FFENER(4) + FFENER(5))
               GAMMAE(IAXIS,IAXIS) = BG/(ESTEP**4)
               BG = (1.0D+00/2.0D+00)*(DIP(4,IAXIS) - DIP(5,IAXIS)) -
     *              (DIP(2,IAXIS) - DIP(3,IAXIS))
               GAMMAD(IAXIS,IAXIS) = BG/(ESTEP**3)
      END DO
C
C     **************************************
C     OFF-AXIAL COMPONENTS ONLY IF REQUESTED
C     **************************************
      IF (OFFDIA) THEN
        NFLDS = 13
C       *******************
C       LOOP OVER THE PLANE
C       *******************
        DO IAXIS = 1,3
           DO JAXIS = IAXIS,3
              IF (IAXIS .NE. JAXIS) THEN
C          **********************************
C          BEGIN LOOP OVER ALL APPLIED FIELDS
C          **********************************
               NPRINT = -5
            DO 500 IFLD = 2,NFLDS
C
               FFENER(IFLD) = ZERO
               DIP(IFLD,1)  = ZERO
               DIP(IFLD,2)  = ZERO
               DIP(IFLD,3)  = ZERO
               EFVAL(1)     = ZERO
               EFVAL(2)     = ZERO
               EFVAL(3)     = ZERO
C
C            SETUP THIS ELECTRIC FIELD
C
               EFVAL(IAXIS) = EFIELD(KFLI(IFLD))
               EFVAL(JAXIS) = EFIELD(KFLJ(IFLD))
               MFLD = IFLD-1
               IF(MASWRK) WRITE(IW,9090) MFLD,(EFVAL(J),J=1,3)
               IF(EXETYP.EQ.CHECK) GO TO 500
C
               IF (IFLD .EQ. 2) THEN
                  FFENER(IFLD) = ENSAV(IAXIS,3)
                  DIP(IFLD,1)  = DPSAV(IAXIS,3,1)
                  DIP(IFLD,2)  = DPSAV(IAXIS,3,2)
                  DIP(IFLD,3)  = DPSAV(IAXIS,3,3)
               ELSE IF (IFLD .EQ. 4) THEN
                  FFENER(IFLD) = ENSAV(JAXIS,2)
                  DIP(IFLD,1)  = DPSAV(JAXIS,2,1)
                  DIP(IFLD,2)  = DPSAV(JAXIS,2,2)
                  DIP(IFLD,3)  = DPSAV(JAXIS,2,3)
               ELSE IF (IFLD .EQ. 6) THEN
                  FFENER(IFLD) = ENSAV(IAXIS,2)
                  DIP(IFLD,1)  = DPSAV(IAXIS,2,1)
                  DIP(IFLD,2)  = DPSAV(IAXIS,2,2)
                  DIP(IFLD,3)  = DPSAV(IAXIS,2,3)
               ELSE IF (IFLD .EQ. 8) THEN
                  FFENER(IFLD) = ENSAV(JAXIS,3)
                  DIP(IFLD,1)  = DPSAV(JAXIS,3,1)
                  DIP(IFLD,2)  = DPSAV(JAXIS,3,2)
                  DIP(IFLD,3)  = DPSAV(JAXIS,3,3)
               ELSE
C
C           **********************************************
C           ADD ELECTRIC FIELD TO ONE ELECTRON HAMILTONIAN
C           **********************************************
               CALL VALFM (LOADFM)
               LMU  = LOADFM + 1
               LH   = LMU    + L2
               LAST = LH     + L2
               NEED = LAST - LOADFM - 1
               CALL GETFM(NEED)
C
               CALL DAREAD(IDAF,IODA,X(LH),L2,94,0)
C
               IF(EFVAL(1).NE.ZERO) THEN
                 CALL DAREAD(IDAF,IODA,X(LMU),L2,95,0)
                 CALL DAXPY(L2,EFVAL(1),X(LMU),1,X(LH),1)
               END IF
               IF(EFVAL(2).NE.ZERO) THEN
                 CALL DAREAD(IDAF,IODA,X(LMU),L2,96,0)
                 CALL DAXPY(L2,EFVAL(2),X(LMU),1,X(LH),1)
               END IF
               IF(EFVAL(3).NE.ZERO) THEN
                 CALL DAREAD(IDAF,IODA,X(LMU),L2,97,0)
                 CALL DAXPY(L2,EFVAL(3),X(LMU),1,X(LH),1)
               END IF
C
               CALL DAWRIT(IDAF,IODA,X(LH),L2,11,0)
C           *******************************************
C           RECOVER THE FIELD FREE ORBITALS AND DENSITY
C           *CAUTION, STORAGE AT -LMU- IS USED AS L3!*
C           *******************************************
               CALL FFSVMO(1,X(LMU))
               IF(MCCI) CALL FFSVCI(1)
               CALL ORTHDN
               CALL RETFM(NEED)
C           ***************************************
C           CALCULATE WAVEFUNCTION UNDER THIS FIELD
C           ***************************************
               CALL WFN
C
               IF(E.EQ.ZERO) THEN
               IF(MASWRK) WRITE(IW,9030)
                 CALL ABRT
                 STOP
               END IF
C
               IF(TIM .GE. TIMLIM) THEN
               IF(MASWRK) WRITE(IW,9040)
                 CALL ABRT
                 STOP
               END IF
C           ***********************************
C           ADD NUCLEAR CONTRIBUTION FROM FIELD
C           ***********************************
               ANUCF = ZERO
               DO JNUF = 1,NAT
                  DO INUF = 1,3
                     COMCOR = C(INUF,JNUF) - COM(INUF)
                     ANUCF = ANUCF - EFVAL(INUF)*COMCOR*ZAN(JNUF)
                  END DO
               END DO
C
               EWITHF = E + ANUCF
               FFENER(IFLD) = EWITHF
C           ****************
C           CALCULATE DIPOLE
C           ****************
               IF(CCTYP .NE. RNONE .OR.
     *           (MPLEVL .GT. 0 .AND. MPPROP .EQ. 0))THEN
                 DIP(IFLD,1) = ZERO
                 DIP(IFLD,2) = ZERO
                 DIP(IFLD,3) = ZERO
               ELSE IF(MPLEVL .GT. 0 .AND. MPPROP .EQ. 1)THEN
                 CALL ELMOMC
                 DIP(IFLD,1) = DIPMX/DFAC
                 DIP(IFLD,2) = DIPMY/DFAC
                 DIP(IFLD,3) = DIPMZ/DFAC
               ELSE
                 CALL ELMOMC
                 DIP(IFLD,1) = DIPMX/DFAC
                 DIP(IFLD,2) = DIPMY/DFAC
                 DIP(IFLD,3) = DIPMZ/DFAC
               END IF
               END IF
C
               IF(MASWRK) THEN
                 MFLD = IFLD-1
                 WRITE(IW,9050) MFLD,(EFVAL(J),J=1,3)
                 WRITE(IW,9060) EWITHF
                 WRITE(IW,9070) DIPMX,DIPMY,DIPMZ
               END IF
  500 CONTINUE
C
               AIJ = (FFENER(10)-FFENER(11)-FFENER(13)+FFENER(12)) /
     *                              48.0D+00
     *              -(FFENER(5)-FFENER(7)-FFENER(3)+FFENER(9))/3.0D+00
               ALPHAE(IAXIS,JAXIS) = AIJ/(ESTEP*ESTEP)
               ALPHAE(JAXIS,IAXIS) = ALPHAE(IAXIS,JAXIS)
C
               BIJJ = 0.5D+00*(FFENER(9)-FFENER(5)+FFENER(3)-FFENER(7))
     *                       + FFENER(6) - FFENER(2)
               BETAE(IAXIS,JAXIS) = BIJJ/(ESTEP**3)
               BJII = 0.5D+00*(FFENER(9)-FFENER(5)+FFENER(7)-FFENER(3))
     *                       + FFENER(4) - FFENER(8)
               BETAE(JAXIS,IAXIS) = BJII/(ESTEP**3)
C
               GIIJJ = -4.0D+00*FFENER(1) - (FFENER(5)+FFENER(9) +
     *                                       FFENER(7)+FFENER(3))
     *               + 2.0D+00*(FFENER(6) + FFENER(2))
     *               + 2.0D+00*(FFENER(4) + FFENER(8))
               GAMMAE(IAXIS,JAXIS) = GIIJJ/(ESTEP**4)
               GAMMAE(JAXIS,IAXIS) = GAMMAE(IAXIS,JAXIS)
C
               GDIP = 0.5D+00*(DIP(5,IAXIS) - DIP(3,IAXIS) +
     *                         DIP(7,IAXIS) - DIP(9,IAXIS))
     *                      - (DIP(6,IAXIS) - DIP(2,IAXIS))
               GAMMAD(IAXIS,JAXIS) = GDIP/(ESTEP**3)
               GAMMAD(JAXIS,IAXIS) = GAMMAD(IAXIS,JAXIS)
C
C
           END IF
           END DO
        END DO
      END IF
C           ZAP OUT SMALL COMPONENTS
C
      DO J=1,3
         IF(ABS(FDIPOL(J)).LT.ZAPTOL) FDIPOL(J) = ZERO
         IF(ABS(DDIPOL(J)).LT.ZAPTOL) DDIPOL(J) = ZERO
         DO I=1,3
            IF(ABS(ALPHAE(I,J)).LT.ZAPTOL) ALPHAE(I,J) = ZERO
            IF(ABS( BETAE(I,J)).LT.ZAPTOL)  BETAE(I,J) = ZERO
            IF(ABS(GAMMAE(I,J)).LT.ZAPTOL) GAMMAE(I,J) = ZERO
            IF(ABS(ALPHAD(I,J)).LT.ZAPTOL) ALPHAD(I,J) = ZERO
            IF(ABS( BETAD(I,J)).LT.ZAPTOL)  BETAD(I,J) = ZERO
            IF(ABS(GAMMAD(I,J)).LT.ZAPTOL) GAMMAD(I,J) = ZERO
         ENDDO
      ENDDO
C
      IF(MASWRK) THEN
        WRITE(IW,8050)
        WRITE(IW,8010) (FDIPOL(K),K=1,3)
        WRITE(IW,8020) ((ALPHAE(I,J),J=1,3),I=1,3)
        WRITE(IW,8030) ((BETAE(I,J),I=1,3),J=1,3)
        WRITE(IW,8040) ((GAMMAE(I,J),J=1,3),I=1,3)
        WRITE(IW,8060)
        WRITE(IW,8010) (DDIPOL(K),K=1,3)
        WRITE(IW,8020) ((ALPHAD(I,J),J=1,3),I=1,3)
        WRITE(IW,8030) ((BETAD(I,J),J=1,3),I=1,3)
        WRITE(IW,8040) ((GAMMAD(I,J),J=1,3),I=1,3)
      END IF
C
  888 CONTINUE
      IF(.NOT.SYM) CALL SYMON
        PK = PKSAVE
        NOPK=NOPKSV
      IF(MASWRK) THEN
        WRITE(IW,*) '..... DONE WITH FINITE FIELD CALCULATION .....'
        CALL TIMIT(1)
      END IF
      RETURN
C
 9000 FORMAT(//'      ******************************'/
     *         '      *                            *'/
     *         '      *  FINITE FIELD CALCULATION  *'/
     *         '      *                            *'/
     *         '      ******************************')
 9010 FORMAT(/5X,13(1H-)/5X,'$FFCALC INPUT'/5X,13(1H-))
 9015 FORMAT(1X,'ONEFLD=',L1,'   SYM=',L1/
     *       1X,'USER SPECIFIED FIELD IS E(X,Y,Z)=',3F9.5)
 9016 FORMAT(1X,'FULL TENSOR OPTIONS: OFFDIA=',L1,'   ESTEP=',F9.5)
 9020 FORMAT(//1X,'CALCULATING ENERGY WITH NO FIELD PRESENT')
 9030 FORMAT(1X,'*** FFIELD RUN STOPPED BY CONVERGENCE FAILURE ***')
 9040 FORMAT(1X,'*** INSUFFICIENT TIME TO CONTINUE FFIELD JOB ***')
 9050 FORMAT(//1X,'APPLIED FIELD',I3,' IS E(X,Y,Z) = ',3F9.5)
 9060 FORMAT(1X,'ENERGY WITH THE FIELD ON ',F30.16)
 9070 FORMAT(1X,'X DIPOLE WITH FIELD ON   ',F30.16/
     *       1X,'Y DIPOLE WITH FIELD ON   ',F30.16/
     *       1X,'Z DIPOLE WITH FIELD ON   ',F30.16//)
 9080 FORMAT(1X,'RECALCULATING INTEGRALS W/O POINT GROUP SYMMETRY')
 9090 FORMAT(//1X,'--- EVALUATING WAVEFUNCTION FOR FIELD',I3,' ---'/
     *          1X,'APPLIED AT STRENGTH E(X,Y,Z) = ',3F9.5)
 8010 FORMAT(1X,'DIPOLE',1X,'#',8X,'X',18X,'Y',18X,'Z',3X,'(A.U.)'/,
     *    1X,'##########################################################
     ##########',/,
     *       8X,'#',1X,1P,E15.7,4X,E15.7,4X,E15.7,0P/)
 8020 FORMAT(2X,'ALPHA',1X,'#',8X,'X',18X,'Y',18X,'Z',3X,'(A.U.)'/,
     *    1X,'##########################################################
     ##########',/,
     *       3X,'X',4X,'#',1X,1P,E15.7,4X,E15.7,4X,E15.7,0P/,
     *       3X,'Y',4X,'#',1X,1P,E15.7,4X,E15.7,4X,E15.7,0P/,
     *       3X,'Z',4X,'#',1X,1P,E15.7,4X,E15.7,4X,E15.7,0P/)
 8030 FORMAT(2X,'BETA',2X,'#',8X,'X',18X,'Y',18X,'Z',3X,'(A.U.)'/,
     *    1X,'##########################################################
     ##########',/,
     *       2X,'XX',4X,'#',1X,1P,E15.7,4X,E15.7,4X,E15.7,0P/,
     *       2X,'YY',4X,'#',1X,1P,E15.7,4X,E15.7,4X,E15.7,0P/,
     *       2X,'ZZ',4X,'#',1X,1P,E15.7,4X,E15.7,4X,E15.7,0P/)
 8040 FORMAT(2X,'GAMMA',1X,'#',8X,'XX',17X,'YY',17X,'ZZ',3X,'(A.U.)'/,
     *    1X,'##########################################################
     ##########',/,
     *       2X,'XX',4X,'#',1X,1P,E15.7,4X,E15.7,4X,E15.7,0P/,
     *       2X,'YY',4X,'#',1X,1P,E15.7,4X,E15.7,4X,E15.7,0P/,
     *       2X,'ZZ',4X,'#',1X,1P,E15.7,4X,E15.7,4X,E15.7,0P//)
 8050 FORMAT(//,19X,'----ENERGY BASED RESULTS----',12X/)
 8060 FORMAT(//,19X,'----DIPOLE BASED RESULTS----',12X/)
      END
