FUNCTION DIPOLE (P,Q,COORD,DIPVEC, MODE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM), NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /DIPSTO/ UX,UY,UZ,CH(NUMATM)
      COMMON /KEYWRD/ KEYWRD
      COMMON /MOLMEC/ HTYPE(4),NHCO(4,20),NNHCO,ITYPE
      COMMON /ISTOPE/ AMS(107)
      COMMON /MULTIP/ DD(107), QQ(107), AM(107), AD(107), AQ(107)
      COMMON /NUMCAL/ NUMCAL
      DIMENSION P(*),Q(*),COORD(3,*),DIPVEC(3), CENTER(3)
      CHARACTER*241 KEYWRD
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)
      SAVE FIRST, HYF, DIP, WTMOL, CHARGD, FORCE
      LOGICAL FIRST, FORCE, CHARGD
      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/
      FIRST=(ICALCN.NE.NUMCAL)
      ICALCN=NUMCAL
      IF (FIRST) THEN
         DO 10 I=2,107
   10    HYF(I,1)= 5.0832*DD(I)
         WTMOL=0.D0
         SUM=0.D0
         DO 20 I=1,NUMAT
            WTMOL=WTMOL+AMS(NAT(I))
   20    SUM=SUM+Q(I)
         CHARGD=(ABS(SUM).GT.0.5D0)
         FORCE=(INDEX(KEYWRD,'FORCE') +INDEX(KEYWRD,'IRC').NE. 0)
         KTYPE=1
         IF(ITYPE.EQ.4)KTYPE=2
      ENDIF
      IF(.NOT.FORCE.AND.CHARGD)THEN
C
C   NEED TO RESET ION'S POSITION SO THAT THE CENTER OF MASS IS AT THE
C   ORIGIN.
C
         DO 30 I=1,3
   30    CENTER(I)=0.D0
         DO 40 I=1,3
            DO 40 J=1,NUMAT
   40    CENTER(I)=CENTER(I)+AMS(NAT(J))*COORD(I,J)
         DO 50 I=1,3
   50    CENTER(I)=CENTER(I)/WTMOL
         DO 60 I=1,3
            DO 60 J=1,NUMAT
   60    COORD(I,J)=COORD(I,J)-CENTER(I)
      ENDIF
      DO 70 I=1,4
         DO 70 J=1,3
   70 DIP(I,J)=0.0D00
      DO 90 I=1,NUMAT
         NI=NAT(I)
         IA=NFIRST(I)
         L=NLAST(I)-IA
         DO 80 J=1,L
            K=((IA+J)*(IA+J-1))/2+IA
   80    DIP(J,2)=DIP(J,2)-HYF(NI,KTYPE)*P(K)
         DO 90 J=1,3
   90 DIP(J,1)=DIP(J,1)+4.803D00*Q(I)*COORD(J,I)
      DO 100 J=1,3
  100 DIP(J,3)=DIP(J,2)+DIP(J,1)
      DO 110 J=1,3
  110 DIP(4,J)=SQRT(DIP(1,J)**2+DIP(2,J)**2+DIP(3,J)**2)
      IF( FORCE) THEN
         DIPVEC(1)=DIP(1,3)
         DIPVEC(2)=DIP(2,3)
         DIPVEC(3)=DIP(3,3)
      ENDIF
      IF(MODE.EQ.1)WRITE (6,130) ((DIP(I,J),I=1,4),J=1,3)
C     STORE DIPOLE MOMENT COMPONENTS IN UX,UY,UZ FOR USE IN
C     ASSIGNING CHARGES DETERMINED FROM THE ESP. BHB
      UX=DIP(1,3)
      UY=DIP(2,3)
      UZ=DIP(3,3)
C
C     STORE CHARGES Q IN ARRAY CH FOR USE IN ASSIGNING SYMMETRY TO
C     CHARGES. BHB
      DO 120 I=1,NUMAT
  120 CH(I)=Q(I)
      DIPOLE = DIP(4,3)
      RETURN
C
  130 FORMAT (' DIPOLE',11X,'X         Y         Z       TOTAL',/,
     1' POINT-CHG.',4F10.3/,' HYBRID',4X,4F10.3/,' SUM',7X,4F10.3)
C
      END