      SUBROUTINE FFHPOL (HEAT0,ATPOL,DIPVEC)
      implicit double precision (A-H,O-Z)
      character*1 AXIS(3)
      logical LARGE,POLDIP,DEBUG
      include 'SIZES'
      include 'MOP6FILE'
C***********************************************************************
C  SUBROUTINE FOR THE FINITE FIELD CALCULATION OF ELECTRIC RESPONSE
C  PROPERTIES (DIPOLE MOMENT, POLARIZABILITY, AND 1ST AND 2ND
C  HYPERPOLARIZABILITY.
C
C  HENRY A. KURTZ, DEPARTMENT OF CHEMISTRY
C                  MEMPHIS STATE UNIVERSITY
C                  MEMPHIS, TN   38152
C
C***********************************************************************
      common /CORE  / CORE(107)
      common /GEOM  / GEO(3,NUMATM)
      common /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM),NORS,NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      common /COORD / COORD(3,NUMATM)
      common /KEYWRD/ KEYWRD
      common /FIELD / EFIELD(3)
      common /EULER / TVEC(3,3),IDTVEC
      character*241 KEYWRD
C
C
C     DIPE4 AND DIPDP HOLD THE CALCULATED DIPOLE MOMENTS
C
C     APOLE4 AND APOLDP HOLD THE POLARIZABILITY TENSOR AS
C                                A PACKED ARRAY XX,XY,YY,XZ,YZ,ZZ
C
C     BETAE4 AND BETAEP HOLD THE FIRST HYPERPOLARIZABILITY
C                                1. XXX
C                                2. YYY     6. YXX
C                                3. ZZZ     7. YZZ
C                                4. XYY     8. ZXX
C                                5. XZZ     9. ZYY
C
      dimension HEATE(3,2),DIPVEC(3),
     1          DIPE4(3),APOLE4(6),BETAE4(9),GAMME4(6),
     2          DIPDP(3),APOLDP(6),BETADP(9),GAMMDP(6),
     3          DIP1P(3),DIP1M(3),DIP2P(3),DIP2M(3)
      dimension IPTBD(6)
      data IPTBD /5,7,4,9,6,8/
C Energy: a.u. to kcal/mole
      AUTOKC = 23.061D+00*27.2107D+00
C Length: a.u. to Angstrom
      AUTOA  = 0.529177D+00
C Dipole: a.u. to debye
      AUTODB = 2.541563D+00
C Electric Field: a.u. to volt/meter
      AUTOVM = 51.4257D+00
      NBDIP = 1
      NGDIP = 4
      NBCNT = 4
      NGCNT = 4
C
      data AXIS/'X','Y','Z'/
      LARGE = (INDEX(KEYWRD,'LARGE').NE.0)
      DEBUG = (INDEX(KEYWRD,'DEBUG').NE.0)
C
C  FIELD STRENGTH IN A.U.
C
      EFVAL=0.001D0
      IDIP=1
C modification for variable field strength
      if (INDEX(KEYWRD,'POLAR=').NE.0)
     1           EFVAL = READA(KEYWRD,INDEX(KEYWRD,'POLAR='))
      write(MIOB,10) EFVAL
   10 format(//' APPLIED ELECTRIC FIELD MAGNITUDE: ',F15.5)
      POLDIP = .FALSE.
      if (IDIP.NE.0) POLDIP = .TRUE.
      SFE = 1.D00/EFVAL
      write(MIOB,20) 6.74834*ATPOL
   20 format(//' ATOMIC CONTRIBUTION TO THE POLARIZABILITY: ',F15.6,/,
     1          '  (IT IS ONLY APPLIED TO THE E4 RESULT)')
C.......................................................................
C  CALCULATE THE POLARIZABILITY AND HYPERPOLARIZABILITIES ALONG
C  THE THREE PRINCIPLE AXES.  (THESE AXES DEPEND ON YOUR ARBITRARY
C  ORIENTATION AND MAY NOT BE THE TRUE PRINCIPLE AXES.)
C.......................................................................
      do 150 ID = 1,3
         if (DEBUG) then
            write(MIOB,30) AXIS(ID)
   30       format(//,' ****** FIELD IN ',A1,' DIRECTION *****',/)
         end if
C
C ZERO THE FIELD
C
         do 40 I = 1,3
            EFIELD(I) = 0.0D00
   40    continue
         HNUC = 0.0D00
         do 50 I = 1,NUMAT
            HNUC = HNUC + EFVAL*GEO(ID,I)*CORE(NAT(I))*AUTOVM
   50    continue
         HNUC = HNUC*23.061D00
C +E(ID)
         EFIELD(ID) = EFVAL
         call COMPFG(GEO,.TRUE.,HEAT1P,.TRUE.,GRAD,.FALSE.)
         call DIPIND (DIP1P)
         DIIP = DIP1P(ID)
C -E(ID)
         EFIELD(ID) = -EFVAL
         call COMPFG(GEO,.TRUE.,HEAT1M,.TRUE.,GRAD,.FALSE.)
         call DIPIND (DIP1M)
         DIIM = DIP1M(ID)
C +2E(ID)
         EFIELD(ID) = 2.0D00*EFVAL
         call COMPFG(GEO,.TRUE.,HEAT2P,.TRUE.,GRAD,.FALSE.)
         call DIPIND (DIP2P)
C -2E(ID)
         EFIELD(ID) = -2.0D00*EFVAL
         call COMPFG(GEO,.TRUE.,HEAT2M,.TRUE.,GRAD,.FALSE.)
         call DIPIND (DIP2M)
C
C  CORRECT FOR ELECTRIC FIELD - NUCLEAR INTERACTIONS
C
         HEAT1P = HEAT1P + HNUC
         HEATE(ID,1) = HEAT1P
         HEAT1M = HEAT1M - HNUC
         HEATE(ID,2) = HEAT1M
         HEAT2P = HEAT2P + HNUC*2.D00
         HEAT2M = HEAT2M - HNUC*2.D00
C
         if (DEBUG) then
            write(MIOB,60)
   60       format(' FIELDS OF: ',5X,'F',21X,'2F')
            write(MIOB,70) HEAT1P,HEAT2P,HEAT1M,HEAT2M,
     %                 DIP1P(ID),DIP2P(ID),DIP1M(ID),DIP2M(ID)
   70       format('  ENERGY:'/,
     1           '   + ',2(F20.10,3X),/,'   - ',2(F20.10,3X),/,
     2           '  DIPOLE:'/,
     3           '   + ',2(F20.10,3X),/,'   - ',2(F20.10,3X))
         end if
C
C DIPOLE
C
         ETERM = (1.0D00/12.D00)*(HEAT2P - HEAT2M)
     1        - (2.0D00/3.0D00)*(HEAT1P - HEAT1M)
         DIPE4(ID) = ETERM*SFE/AUTOKC
C
C ALPHA
C
         IVL = (ID*(ID+1))/2
         ETERM = 2.5D00*HEAT0 - (4.D00/3.D00)*(HEAT1P + HEAT1M)
     1        + (1.D00/12.0D00)*(HEAT2P + HEAT2M)
         APOLE4(IVL) = ETERM*SFE*SFE/AUTOKC + ATPOL*6.74834
C
C BETA
C
         ETERM = (HEAT1P - HEAT1M) - 0.5D00*(HEAT2P - HEAT2M)
         BETAE4(ID) = ETERM*SFE*SFE*SFE/AUTOKC
C
C GAMMA
C
         ETERM = 4.0D00*(HEAT1P + HEAT1M) - (HEAT2P + HEAT2M)
     1        - 6.0D00*HEAT0
         GAMME4(ID) = ETERM*SFE*SFE*SFE*SFE/AUTOKC
C
C DIPOLE CALCULATIONS
C
         DMU = (2.0D00/3.0D00)*(DIP1P(ID) + DIP1M(ID))
     1      - (1.D00/6.0D00)*(DIP2P(ID) + DIP2M(ID))
         DIPDP(ID) = DMU/AUTODB
         AE = (2.0D00/3.0D00)*(DIP1P(ID) - DIP1M(ID))
     1     - (1.0D00/12.D00)*(DIP2P(ID) - DIP2M(ID))
         APOLDP(IVL) = AE*SFE/AUTODB
         BE = (1.D00/3.0D00)*(DIP2P(ID) + DIP2M(ID)
     1                     - DIP1P(ID) - DIP1M(ID))
         BETADP(ID) = BE*SFE*SFE/AUTODB
         GE = 0.5D00*(DIP2P(ID) - DIP2M(ID))
     1     - (DIP1P(ID) - DIP1M(ID))
         GAMMDP(ID) = GE*SFE*SFE*SFE/AUTODB
         do 80 KD = 1,3
            if (KD.LT.ID) then
               KVL = (ID*(ID-1))/2 + KD
               AKI = (2.0D00/3.0D00)*(DIP1P(KD) - DIP1M(KD))
     1         - (1.0D00/12.0D00)*(DIP2P(KD) - DIP2M(KD))
               APOLDP(KVL) = AKI*SFE/AUTODB
            end if
            if (KD.NE.ID) then
               BKII = (1.0D00/3.0D00)*(DIP2P(KD) + DIP2M(KD)
     1                           - DIP1P(KD) - DIP1M(KD))
               NBD = IPTBD(NBDIP)
               BETADP(NBD) = BKII*SFE*SFE/AUTODB
               NBDIP = NBDIP + 1
            end if
   80    continue
C.......................................................................
C
C  NOW CALCULATE THE OFF AXIS RESULTS.
C
C.......................................................................
         IDM1 = ID - 1
         do 140 JD = 1,IDM1
            HNUCJ = 0.0D00
            do 90 I = 1,NUMAT
               HNUCJ = HNUCJ + EFVAL*GEO(JD,I)*CORE(NAT(I))*51.4257
   90       continue
            HNUCJ = HNUCJ*23.061
            do 100 I = 1,3
               EFIELD(I) = 0.0D00
  100       continue
C
C DIAGONAL FIELDS WITH COMPONENTS EQUAL TO EFVAL
C
            EFIELD(ID) = EFVAL
            EFIELD(JD) = EFVAL
            call COMPFG(GEO,.TRUE.,HPP,.TRUE.,GRAD,.FALSE.)
            call DIPIND (DIP1P)
            DPP = DIP1P(ID)
            EFIELD(JD) = -EFVAL
            call COMPFG(GEO,.TRUE.,HPM,.TRUE.,GRAD,.FALSE.)
            call DIPIND (DIP1P)
            DPM = DIP1P(ID)
            EFIELD(ID) = -EFVAL
            call COMPFG(GEO,.TRUE.,HMM,.TRUE.,GRAD,.FALSE.)
            call DIPIND (DIP1P)
            DMM = DIP1P(ID)
            EFIELD(JD) = EFVAL
            call COMPFG(GEO,.TRUE.,HMP,.TRUE.,GRAD,.FALSE.)
            call DIPIND (DIP1P)
            DMP = DIP1P(ID)
            HPP = HPP + HNUC + HNUCJ
            HPM = HPM + HNUC - HNUCJ
            HMM = HMM - HNUC - HNUCJ
            HMP = HMP - HNUC + HNUCJ
            if (DEBUG) then
               write(MIOB,110)
  110          format(/,' ',12X,'+,+',15X,'+,-',15X,'-,+',15X,'-,-')
               write(MIOB,120) HPP,HPM,HMP,HMM
  120          format('  E ',4F15.6)
            end if
C
C  DIAGONAL FIELDS WITH COMPONENTS EQUAL TO 2*EFVAL
C
            EFIELD(ID) = EFVAL*2.D00
            EFIELD(JD) = EFVAL*2.D00
            call COMPFG(GEO,.TRUE.,H2PP,.TRUE.,GRAD,.FALSE.)
            EFIELD(JD) = -EFVAL*2.D00
            call COMPFG(GEO,.TRUE.,H2PM,.TRUE.,GRAD,.FALSE.)
            EFIELD(ID) = -EFVAL*2.D00
            call COMPFG(GEO,.TRUE.,H2MM,.TRUE.,GRAD,.FALSE.)
            EFIELD(JD) = EFVAL*2.D00
            call COMPFG(GEO,.TRUE.,H2MP,.TRUE.,GRAD,.FALSE.)
            H2PP = H2PP + 2.0D00*(HNUC + HNUCJ)
            H2PM = H2PM + 2.0D00*(HNUC - HNUCJ)
            H2MM = H2MM - 2.0D00*(HNUC + HNUCJ)
            H2MP = H2MP - 2.0D00*(HNUC - HNUCJ)
            if (DEBUG) then
               write(MIOB,130) H2PP,H2PM,H2MP,H2MM
  130          format(' 2E ',4F15.6)
            end if
C
            ATERM = (1.0D00/48.0D00)*(H2PP - H2PM - H2MP + H2MM)
     1          - (1.0D00/3.0D00)*(HPP - HPM - HMP + HMM)
            AIJ = ATERM*SFE*SFE/AUTOKC
            IVL = (ID*(ID-1))/2 + JD
            APOLE4(IVL) = AIJ
            BTERM = 0.5D00*(HMM - HPP + HPM - HMP)
     1          + HEATE(JD,1) - HEATE(JD,2)
            BJII = BTERM*SFE*SFE*SFE/AUTOKC
            BETAE4(NBCNT) = BJII
            NBCNT = NBCNT + 1
            BTERM = 0.5D00*(HMM - HPP + HMP - HPM)
     1          + HEATE(ID,1) - HEATE(ID,2)
            BIJJ = BTERM*SFE*SFE*SFE/AUTOKC
            BETAE4(NBCNT) = BIJJ
            NBCNT = NBCNT + 1
C
            GTERM = -(HPP + HMM + HPM + HMP) - 4.0D00*HEAT0
     1           + 2.0D00*(HEATE(ID,1) + HEATE(ID,2))
     2           + 2.0D00*(HEATE(JD,1) + HEATE(JD,2))
            GIIJJ = GTERM*SFE*SFE*SFE*SFE/AUTOKC
            GAMME4(NGCNT) = GIIJJ
            GDIP = 0.5D00*(DPP - DMP + DPM - DMM) - (DIIP - DIIM)
            GAMMDP(NGCNT) = GDIP*SFE*SFE*SFE/AUTODB
            NGCNT = NGCNT + 1
  140    continue
C
  150 continue
C-----------------------------------------------------------------------
C  SUMMARIZE THE RESULTS
C-----------------------------------------------------------------------
      write(MIOB,160)
  160 format(//,' ',30('*'),' DIPOLE ',30('*'),//)
      DIPE4T = SQRT(DIPE4(1)*DIPE4(1) + DIPE4(2)*DIPE4(2)
     1              + DIPE4(3)*DIPE4(3))
      DIPE4D = DIPE4T*AUTODB
      DIPDPT = SQRT(DIPDP(1)*DIPDP(1) + DIPDP(2)*DIPDP(2)
     1              + DIPDP(3)*DIPDP(3))
      DIPDPD = DIPDPT*AUTODB
      write(MIOB,170)
  170 format(21X,'E4',13X,'DIP',/)
      write(MIOB,180) 'X',DIPE4(1),DIPDP(1)
      write(MIOB,180) 'Y',DIPE4(2),DIPDP(2)
      write(MIOB,180) 'Z',DIPE4(3),DIPDP(3)
  180 format(5X,A1,7X,2F15.6)
      write(MIOB,190) DIPE4T,DIPDPT,
     1               DIPE4D,DIPDPD
  190 format(//' MAGNITUDE:  ',2F15.6,'  (A.U.)',/,
     1          ' ',12X,2F15.6,'  (DEBYE)')
C
C FIND EIGENVALUES AND EIGENVECTORS OF POLARIZATION MATRIX.
C
      write(MIOB,200)
  200 format(//,' ',22('*'),' POLARIZABILITY (ALPHA)',21('*'),//)
      AVGPE4 = (APOLE4(1)+APOLE4(3)+APOLE4(6))/3.0D00
      AVGA3 = AVGPE4*0.14818D00
      AVGESU = AVGPE4*0.296352D-24
      AVGPDP = (APOLDP(1)+APOLDP(3)+APOLDP(6))/3.0D00
      AVGA3D = AVGPDP*0.14818D00
      AVGESD = AVGPDP*0.296352D-24
      write(MIOB,210)
  210 format('  COMPONENT',12X,'E4',13X,'DIP',/)
      write(MIOB,220) 'XX',APOLE4(1),APOLDP(1),
     1              'YY',APOLE4(3),APOLDP(3),
     2              'ZZ',APOLE4(6),APOLDP(6),
     3              'XY',APOLE4(2),APOLDP(2),
     4              'XZ',APOLE4(4),APOLDP(4),
     5              'YZ',APOLE4(5),APOLDP(5)
  220 format(' ',5X,A4,5X,2F15.6)
      write(MIOB,230) AVGPE4,AVGPDP,AVGA3,AVGA3D,AVGESU,AVGESD
  230 format(//,' AVERAGE POLARIZABILITY:',8X,'E4',13X,'DIP',/,
     1           ' ',24X,2F15.6,'  A.U.',/,
     2           ' ',24X,2F15.6,'  ANG.**3',/,
     3           ' ',24X,2(1PD15.6),'  ESU')
C
C  CALCULATE "EXPERIMENTAL" HYPERPOLARIZABILITIES
C
C   8.65710D-33 is a.u. to e.s.u. conversion
      write(MIOB,240)
  240 format(//,' ',30('*'),' SECOND-ORDER (BETA)',25('*'),//)
      BX4 = 0.6D00*(BETAE4(1) + BETAE4(4) + BETAE4(6))
      BY4 = 0.6D00*(BETAE4(2) + BETAE4(5) + BETAE4(8))
      BZ4 = 0.6D00*(BETAE4(3) + BETAE4(7) + BETAE4(9))
      B4MU = (BX4*DIPE4(1) + BY4*DIPE4(2) + BZ4*DIPE4(3))/DIPE4T
      B4ESU = B4MU*8.65710D-03
      BXD = 0.6D00*(BETADP(1) + BETADP(4) + BETADP(6))
      BYD = 0.6D00*(BETADP(2) + BETADP(5) + BETADP(8))
      BZD = 0.6D00*(BETADP(3) + BETADP(7) + BETADP(9))
      BDMU = (BXD*DIPDP(1) + BYD*DIPDP(2) + BZD*DIPDP(3))/DIPDPT
      BDESU = BDMU*8.65710D-03
C
      write(MIOB,'(29X,A2,25X,A6)')'1X','(1/2)X'
      write(MIOB,250)
  250 format('  COMPONENT',2(12X,'E4',10X,'DIP',2X),/)
      write(MIOB,260) 'XXX',BETAE4(1),BETADP(1),BETAE4(1)/2,BETADP(1)/2
      write(MIOB,260) 'XYY',BETAE4(4),BETADP(4),BETAE4(4)/2,BETADP(4)/2
      write(MIOB,260) 'XZZ',BETAE4(6),BETADP(6),BETAE4(6)/2,BETADP(6)/2
      write(MIOB,260) 'YYY',BETAE4(2),BETADP(2),BETAE4(2)/2,BETADP(2)/2
      write(MIOB,260) 'YXX',BETAE4(5),BETADP(5),BETAE4(5)/2,BETADP(5)/2
      write(MIOB,260) 'YZZ',BETAE4(8),BETADP(8),BETAE4(8)/2,BETADP(8)/2
      write(MIOB,260) 'ZZZ',BETAE4(3),BETADP(3),BETAE4(3)/2,BETADP(3)/2
      write(MIOB,260) 'ZXX',BETAE4(7),BETADP(7),BETAE4(7)/2,BETADP(7)/2
      write(MIOB,260) 'ZYY',BETAE4(9),BETADP(9),BETAE4(9)/2,BETADP(9)/2
  260 format(' ',5X,A4,2(5X,2F12.3))
      write(MIOB,270)
  270 format(//,' VECTOR COMPONENTS GIVEN BY:',/,
     1          '      BI = (2/5)*(BI11+BI22+BI33)'/)
      write(MIOB,280) 'BX',BX4,BXD,BX4/2,BXD/2
      write(MIOB,280) 'BY',BY4,BYD,BY4/2,BYD/2
      write(MIOB,280) 'BZ',BZ4,BZD,BZ4/2,BZD/2
  280 format(' ',6X,A2,2(6X,2F12.3))
      write(MIOB,290)
  290 format(//'  VALUE OF BETA ALONG THE DIPOLE MOMENT:'/)
      write(MIOB,300) B4MU,BDMU,B4MU/2,BDMU/2,B4ESU,BDESU,
     1                  B4ESU/2,BDESU/2
  300 format(' ',4X,'B(AU)',2(5X,2F12.3,2X),/,
     1        ' ',4X,'B(ESU)',4X,2F12.3,7X,2F12.3,' (X10-30)')
C
      write(MIOB,310)
  310 format(//' ',24('*'),' THIRD-ORDER (GAMMA)',24('*'),//)
      GAMVAL = (GAMME4(1) + GAMME4(2) + GAMME4(3))
      GAMVAL = GAMVAL + 2.0D00*(GAMME4(4) + GAMME4(5) + GAMME4(6))
      GAMVAL = GAMVAL/5.0D00
C  5.05116D-40 is the a.u. to e.s.u. conversion
      GAMESU = GAMVAL*5.05116D-04
      GAMDIP = (GAMMDP(1) + GAMMDP(2) + GAMMDP(3))
      GAMDIP = GAMDIP + 2.0D00*(GAMMDP(4) + GAMMDP(5) + GAMMDP(6))
      GAMDIP = GAMDIP/5.0D00
      GAMDES = GAMDIP*5.05116D-04
      write(MIOB,'(23X,A2,25X,A6)')'1X','(1/6)X'
      write(MIOB,320)
  320 format(' ',17X,'E4',8X,'DIP',16X,'E4',8X,'DIP',/)
      write(MIOB,330) 'XXXX',GAMME4(1),GAMMDP(1),GAMME4(1)/6,GAMMDP(1)/6
      write(MIOB,330) 'YYYY',GAMME4(2),GAMMDP(2),GAMME4(2)/6,GAMMDP(2)/6
      write(MIOB,330) 'ZZZZ',GAMME4(3),GAMMDP(3),GAMME4(3)/6,GAMMDP(3)/6
      write(MIOB,330) 'XXYY',GAMME4(4),GAMMDP(4),GAMME4(4)/6,GAMMDP(4)/6
      write(MIOB,330) 'XXZZ',GAMME4(5),GAMMDP(5),GAMME4(5)/6,GAMMDP(5)/6
      write(MIOB,330) 'YYZZ',GAMME4(6),GAMMDP(6),GAMME4(6)/6,GAMMDP(6)/6
  330 format(5X,A4,2F12.3,5X,2F12.3)
      write(MIOB,340)
  340 format(//' AVERAGE GAMMA GIVEN BY:',/,
     1 '    (1/5)*[GXXX + GYYY + GZZZ + 2.0*(GXXYY + GXXZZ + GYYZZ)]')
      write(MIOB,'(/,20X,A2,22X,A6)')'1X','(1/6)X'
      write(MIOB,350) GAMVAL,GAMDIP,GAMVAL/6,GAMDIP/6,
     1              GAMESU,GAMDES,GAMESU/6,GAMDES/6
  350 format(/' <GAMMA> ',1PD12.5,1PD12.5,5X,1PD12.5,1PD12.5,'  A.U.'/,
     1       ' ',8X,1PD12.5,1PD12.5,5X,1PD12.5,1PD12.5,'  ESU (X10-36)')
C
      return
      end
