C
C $Id: dipole.F,v 1.11 1998/12/08 20:47:08 arjan Exp arjan $
C
C------------------------------------------------------------------------
C234567890123456789012345678901234567890123345678901234567890123456789012<stop

      Subroutine dipole
      implicit double precision(a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

C  Dipole.F calculates the dipole of a molecule.  Th variables required by
C  it are all contained in divcon.h, thus none are explicitly passed.  Basic-
C  ally, it calculates the center of mass, moves the atom so the origin resides
C  at that center, then calculates the point charge contribution of each atom to
C  the dipole.  It then calculates the contribution of hybridization to the
C  the dipole, and finally includes all these values in divcon.out.  For more
C  information, please see McGlynn, et al., "Introduction to Applied Quantum
C  Chemistry," 1972.

C  Variable declaration.

      integer counter
      double precision mulpt,weight
      dimension com(3), ptchg(3), ptchg2(3), ptchg3(3)
      dimension hybfac(maxatm), hybdip(3)

C  Initialization block

      DO 5 I=1,3
         com(i)=0.d0
         ptchg(i)=0.d0
         ptchg2(i)=0.d0
         ptchg3(i)=0.d0
         hybdip(i)=0.d0
 5    CONTINUE
      DO 7 j=1,natoms
         hybfac(j)=0.d0
 7    CONTINUE

C  Calculate center of mass.
      weight = 0.0
      DO 10 j=1,natoms
         weight=weight+emass(iatnum(j))
         DO 10 i=1,3
            com(i)=com(i)+xyz(i,j)*emass(iatnum(j))
 10   CONTINUE

      DO 20 i=1,3
         com(i)=com(i)/weight
 20   CONTINUE

C  Move molecule to center of mass.

      DO 30 j=1,natoms
         DO 30 i=1,3
            xyz(i,j)=xyz(i,j)-com(i)
 30   CONTINUE

C  Calculate point charge contribution to dipole for all three charge systems.
C  As the three charge types are contained in atchg, atchg2, and atchg3 for
C  Mulliken, CM1, and CM2 charges, respectively, the point charge array
C  (ptchg) is named using a similar system.  The reason a 2D array was not used
C  is the name of the atchg arrays presented diffculties as they are only 1D
C  arrays.

      DO 40 j=1,natoms
         DO 40 i=1,3
            ptchg(i)=ptchg(i)+4.803d0*atchg(j)*xyz(i,j)
            ptchg2(i)=ptchg2(i)+4.803d0*atchg2(j)*xyz(i,j)
            ptchg3(i)=ptchg3(i)+4.803d0*atchg3(j)*xyz(i,j)
 40   CONTINUE

C  Calculate hybridization factors and hybridization contribution to
C  dipole for the mUlliken charge system.

      DO 50 j=1,natoms
         x=dble(nquant(iatnum(j)))
         const=((5.0832d0)*(2.d0*x+1.d0)*(2.d0**(2.d0*x+1.d0)))
     .    /sqrt(3.d0)

         IF(MNDO) THEN
            frac1=(expnt1(0,iatnum(j))*expnt1(1,iatnum(j)))**
     .       (x+.5d0)
            frac2=(expnt1(0,iatnum(j))+expnt1(1,iatnum(j)))**
     .       (2.d0*x+2.d0)
         ELSEIF (AM1) THEN
            frac1=(expnt2(0,iatnum(j))*expnt2(1,iatnum(j)))**
     .       (x+.5d0)
            frac2=(expnt2(0,iatnum(j))+expnt2(1,iatnum(j)))**
     .       (2.d0*x+2.d0)
         ELSEIF (RM1) THEN
            frac1=(expnt5(0,iatnum(j))*expnt5(1,iatnum(j)))**
     .       (x+.5d0)
            frac2=(expnt5(0,iatnum(j))+expnt5(1,iatnum(j)))**
     .       (2.d0*x+2.d0)
         ELSEIF (AM1D) THEN
            frac1=(expnt6(0,iatnum(j))*expnt6(1,iatnum(j)))**
     .       (x+.5d0)
            frac2=(expnt6(0,iatnum(j))+expnt6(1,iatnum(j)))**
     .       (2.d0*x+2.d0)
         ELSEIF (PM3PDDG) THEN
            frac1=(expnt7(0,iatnum(j))*expnt7(1,iatnum(j)))**
     .       (x+.5d0)
            frac2=(expnt7(0,iatnum(j))+expnt7(1,iatnum(j)))**
     .       (2.d0*x+2.d0)
         ELSEIF (PM3) THEN
            frac1=(expnt3(0,iatnum(j))*expnt3(1,iatnum(j)))**
     .       (x+.5d0)
            frac2=(expnt3(0,iatnum(j))+expnt3(1,iatnum(j)))**
     .       (2.d0*x+2.d0)
         ENDIF
         hybfac(j)=const*frac1/frac2
 50   CONTINUE

      counter=0
      DO 60 j=1,natoms
         IF (iatnum(j).eq.1) then
            counter=counter+1
         ELSE
            hybdip(1)=hybdip(1)-hybfac(j)*PDIAG(counter+2)
            hybdip(2)=hybdip(2)-hybfac(j)*PDIAG(counter+4)
            hybdip(3)=hybdip(3)-hybfac(j)*PDIAG(counter+7)
            counter=counter+10
         ENDIF
 60   CONTINUE

C  Sum all and report the the three dipoles (Mul., CM1, CM2.)
C  There was an error that after much work found.  The hybridization
C  factor should not be included with CM1 and CM2 charges.  In hindsight,
C  this is obvious, as the very nature of the CM1 and CM2 charges accounts
C  for those factors.

      mulpt=sqrt(ptchg(1)**2+ptchg(2)**2+ptchg(3)**2)
      CM1pt=sqrt(ptchg2(1)**2+ptchg2(2)**2+ptchg2(3)**2)
      CM2pt=sqrt(ptchg3(1)**2+ptchg3(2)**2+ptchg3(3)**2)
      hybrid=sqrt(hybdip(1)**2+hybdip(2)**2+hybdip(3)**2)
! GM: 2010-10-31 begin
! bug: correction for the calculation of dipole moment
!     dipol1=mulpt+hybrid
      dipol1=sqrt( (ptchg(1)+hybdip(1))**2
     .            +(ptchg(2)+hybdip(2))**2
     .            +(ptchg(3)+hybdip(3))**2)
! GM: 2010-10-31 end
      dipol2=CM1pt
      dipol3=CM2pt

! GM: 2010-10-31 begin
! add: print dipole moment if screen option is activated
      if (screen) then
        write (iout,*) ' '
        write (iout,66) dipol1,(ptchg(i)+hybdip(i),i=1,3)
   66 format(' DIPOLE Total=',f12.6,' X=',f12.6,' Y=',f12.6,' Z=',f12.6)
        write (iout,*) ' '
!       write (iout,*) ' '
!       write (iout,*) 'DIPOLE:'
!       IF (MNDO) write(iout, *) 'NOTE: MNDO USES ONLY MULLIKEN CHARGES'
!       write (iout,65) hybrid
!       write (iout,67)
!       write (iout,70) mulpt,CM1pt,CM2pt
!       write (iout,80) dipol1,dipol2,dipol3
!       write (iout,*) ' '
!       write (iout,*) ' '
      endif
 65   format(' ','HYBRIDIZATION FACTOR =',3x,F6.4)
 67   format(' ', 17x,'MULLIKEN',5x,'CM1',7x,'CM2')
 70   format(1x,'POINT CHARGE =',4x,F6.4,4x,F6.4,4x,F6.4)
 80   format(1x,'TOTAL =',11x,F6.4,4x,F6.4,4x,F6.4)
! GM: 2010-10-31 end

      return
      end