      SUBROUTINE DIPIND (DIPVEC)
C...............................................................
C  MODIFICATION OF DIPOLE SUBROUTINE FOR USE IN THE CALCULATION
C  OF THE INDUCED DIPOLES FOR POLARIZABILITIES.
C...............................................................
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      common /CORE  / CORE(107)
      common /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
      common /MOLMEC/ HTYPE(4),NHCO(4,20),NNHCO,ITYPE
      common /GEOM  / GEO(3,NUMATM)
      common /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM),NORBS,NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      common /NUMCAL/ NUMCAL
      common /KEYWRD/ KEYWRD
      common /ISTOPE/ AMS(107)
      common /MULTIP/ DD(107), QQ(107), AM(107), AD(107), AQ(107)
      dimension Q(MAXORB),Q2(MAXORB),DIPVEC(3),CENTER(3),
     1          COORD(3,NUMATM)
      character*241 KEYWRD
C
C***********************************************************************
C     DIPOLE CALCULATES DIPOLE MOMENTS
C
C  ON INPUT P     = DENSITY MATRIX
C           Q     = TOTAL ATOMIC CHARGES, (NUCLEAR + ELECTRONIC)
C           NUMAT = NUMBER OF ATOMS IN MOLECULE
C           NAT   = ATOMIC NUMBERS OF ATOMS
C           NFIRST= START OF ATOM ORBITAL COUNTERS
C           COORD = COORDINATES OF ATOMS
C
C  OUTPUT  DIPOLE = DIPOLE MOMENT
C***********************************************************************
C
C     IN THE Zdo APPROXIMATION, ONLY TWO TERMS ARE RETAINED IN THE
C     CALCULATION OF DIPOLE MOMENTS.
C     1. THE POINT CHARGE TERM (INDEPENDENT OF PARAMETERIZATION).
C     2. THE ONE-CENTER HYBRIDIZATION TERM, WHICH ARISES FROM MATRIX
C     ELEMENTS OF THE FORM <NS/R/NP>. THIS TERM IS A FUNCTION OF
C     THE SLATER EXPONENTS (ZS,ZP) AND IS THUS DEPENDENT ON PARAMETER-
C     IZATION. THE HYBRIDIZATION FACTORS (HYF(I)) USED IN THIS SUB-
C     ROUTINE ARE CALCULATED FROM THE FOLLOWING FORMULAE.
C     FOR SECOND ROW ELEMENTS <2S/R/2P>
C     HYF(I)= 469.56193322*(SQRT(((ZS(I)**5)*(ZP(I)**5)))/
C           ((ZS(I) + ZP(I))**6))
C     FOR THIRD ROW ELEMENTS <3S/R/3P>
C     HYF(I)=2629.107682607*(SQRT(((ZS(I)**7)*(ZP(I)**7)))/
C           ((ZS(I) + ZP(I))**8))
C     FOR FOURTH ROW ELEMENTS AND UP :
C     HYF(I)=2*(2.10716)*DD(I)
C     WHERE DD(I) IS THE CHARGE SEPARATION IN ATOMIC UNITS
C
C
C     REFERENCES:
C     J.A.POPLE & D.L.BEVERIDGE: APPROXIMATE M.O. THEORY
C     S.P.MCGLYNN, ET AL: APPLIED QUANTUM CHEMISTRY
C
      dimension DIP(4,3)
      dimension HYF(107,2)
      logical FORCE, CHARGD
      save ICALCN, HYF, WTMOL, CHARGD, FORCE
      data HYF(1,1)     / 0.0D00           /
      data   HYF(1,2) /0.0D0     /
      data   HYF(5,2) /6.520587D0/
      data   HYF(6,2) /4.253676D0/
      data   HYF(7,2) /2.947501D0/
      data   HYF(8,2) /2.139793D0/
      data   HYF(9,2) /2.2210719D0/
      data   HYF(14,2)/6.663059D0/
      data   HYF(15,2)/5.657623D0/
      data   HYF(16,2)/6.345552D0/
      data   HYF(17,2)/2.522964D0/
      data ICALCN/0/
C
C  SETUP FOR DIPOLE CALCULATION
C
      call CHRGE (P,Q2)
      do 10 I = 1,NUMAT
         Q(I) = CORE(NAT(I)) - Q2(I)
   10 continue
      call GMETRY (GEO,COORD)
C
      if (ICALCN.NE.NUMCAL) then
         do 20 I=2,107
   20    HYF(I,1)= 5.0832*DD(I)
         WTMOL=0.D0
         SUM=0.D0
         do 30 I=1,NUMAT
            WTMOL=WTMOL+AMS(NAT(I))
   30    SUM=SUM+Q(I)
         CHARGD=(ABS(SUM).GT.0.5D0)
         ICALCN=NUMCAL
         FORCE=(INDEX(KEYWRD,'FORCE') +INDEX(KEYWRD,'IRC').NE. 0)
         KTYPE=1
         if (ITYPE.EQ.4)  KTYPE=2
      end if
      if (CHARGD) then
C
C   NEED TO RESET ION'S POSITION SO THAT THE CENTER OF MASS IS AT THE
C   ORIGIN.
C
C$DOIT ASIS
         do 40 I=1,3
   40    CENTER(I)=0.D0
         do 50 I=1,3
C$DOIT VBEST
            do 50 J=1,NUMAT
   50    CENTER(I)=CENTER(I)+AMS(NAT(J))*COORD(I,J)
C$DOIT ASIS
         do 60 I=1,3
   60    CENTER(I)=CENTER(I)/WTMOL
         do 70 I=1,3
C$DOIT VBEST
            do 70 J=1,NUMAT
   70    COORD(I,J)=COORD(I,J)-CENTER(I)
      end if
C$DOIT ASIS
      do 80 I=1,4
C$DOIT ASIS
         do 80 J=1,3
   80 DIP(I,J)=0.0D00
C$DOIT ASIS
      do 100 I=1,NUMAT
         NI=NAT(I)
         IA=NFIRST(I)
         L=NLAST(I)-IA
C$DOIT ASIS
         do 90 J=1,L
            K=((IA+J)*(IA+J-1))/2+IA
   90    DIP(J,2)=DIP(J,2)-HYF(NI,KTYPE)*P(K)
C$DOIT ASIS
         do 100 J=1,3
  100 DIP(J,1)=DIP(J,1)+4.803D00*Q(I)*COORD(J,I)
C$DOIT ASIS
      do 110 J=1,3
  110 DIP(J,3)=DIP(J,2)+DIP(J,1)
C$DOIT ASIS
      do 120 J=1,3
  120 DIP(4,J)=SQRT(DIP(1,J)**2+DIP(2,J)**2+DIP(3,J)**2)
      DIPVEC(1)= -DIP(1,3)
      DIPVEC(2)= -DIP(2,3)
      DIPVEC(3)= -DIP(3,3)
      DIPCOM = DIP(4,3)
C      write(MIOB,60) ((DIP(I,J),I=1,4),J=1,3)
C   60 format(3(4F10.3))
      return
C
      end
