mopac606_nbo/src/dipind.f

   1       SUBROUTINE DIPIND (DIPVEC)
   2 C...............................................................
   3 C  MODIFICATION OF DIPOLE SUBROUTINE FOR USE IN THE CALCULATION
   4 C  OF THE INDUCED DIPOLES FOR POLARIZABILITIES.
   5 C...............................................................
   6       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
   7       INCLUDE 'SIZES'
   8       COMMON /CORE  / CORE(107)
   9       COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
  10       COMMON /MOLMEC/ HTYPE(4),NHCO(4,20),NNHCO,ITYPE
  11       COMMON /GEOM  / GEO(3,NUMATM)
  12       COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
  13      1                NLAST(NUMATM),NORBS,NELECS,NALPHA,NBETA,
  14      2                NCLOSE,NOPEN,NDUMY,FRACT
  15       COMMON /NUMCAL/ NUMCAL
  16       COMMON /KEYWRD/ KEYWRD
  17       COMMON /ISTOPE/ AMS(107)
  18       COMMON /MULTIP/ DD(107), QQ(107), AM(107), AD(107), AQ(107)
  19       DIMENSION Q(MAXORB),Q2(MAXORB),DIPVEC(3),CENTER(3),
  20      1          COORD(3,NUMATM)
  21       CHARACTER*241 KEYWRD
  22 C
  23 C***********************************************************************
  24 C     DIPOLE CALCULATES DIPOLE MOMENTS
  25 C
  26 C  ON INPUT P     = DENSITY MATRIX
  27 C           Q     = TOTAL ATOMIC CHARGES, (NUCLEAR + ELECTRONIC)
  28 C           NUMAT = NUMBER OF ATOMS IN MOLECULE
  29 C           NAT   = ATOMIC NUMBERS OF ATOMS
  30 C           NFIRST= START OF ATOM ORBITAL COUNTERS
  31 C           COORD = COORDINATES OF ATOMS
  32 C
  33 C  OUTPUT  DIPOLE = DIPOLE MOMENT
  34 C***********************************************************************
  35 C
  36 C     IN THE ZDO APPROXIMATION, ONLY TWO TERMS ARE RETAINED IN THE
  37 C     CALCULATION OF DIPOLE MOMENTS.
  38 C     1. THE POINT CHARGE TERM (INDEPENDENT OF PARAMETERIZATION).
  39 C     2. THE ONE-CENTER HYBRIDIZATION TERM, WHICH ARISES FROM MATRIX
  40 C     ELEMENTS OF THE FORM <NS/R/NP>. THIS TERM IS A FUNCTION OF
  41 C     THE SLATER EXPONENTS (ZS,ZP) AND IS THUS DEPENDENT ON PARAMETER-
  42 C     IZATION. THE HYBRIDIZATION FACTORS (HYF(I)) USED IN THIS SUB-
  43 C     ROUTINE ARE CALCULATED FROM THE FOLLOWING FORMULAE.
  44 C     FOR SECOND ROW ELEMENTS <2S/R/2P>
  45 C     HYF(I)= 469.56193322*(SQRT(((ZS(I)**5)*(ZP(I)**5)))/
  46 C           ((ZS(I) + ZP(I))**6))
  47 C     FOR THIRD ROW ELEMENTS <3S/R/3P>
  48 C     HYF(I)=2629.107682607*(SQRT(((ZS(I)**7)*(ZP(I)**7)))/
  49 C           ((ZS(I) + ZP(I))**8))
  50 C     FOR FOURTH ROW ELEMENTS AND UP :
  51 C     HYF(I)=2*(2.10716)*DD(I)
  52 C     WHERE DD(I) IS THE CHARGE SEPARATION IN ATOMIC UNITS
  53 C
  54 C
  55 C     REFERENCES:
  56 C     J.A.POPLE & D.L.BEVERIDGE: APPROXIMATE M.O. THEORY
  57 C     S.P.MCGLYNN, ET AL: APPLIED QUANTUM CHEMISTRY
  58 C
  59       DIMENSION DIP(4,3)
  60       DIMENSION HYF(107,2)
  61       SAVE ICALCN, HYF, WTMOL, CHARGD
  62       LOGICAL  CHARGD
  63       DATA HYF(1,1)     / 0.0D00           /
  64       DATA   HYF(1,2) /0.0D0     /
  65       DATA   HYF(5,2) /6.520587D0/
  66       DATA   HYF(6,2) /4.253676D0/
  67       DATA   HYF(7,2) /2.947501D0/
  68       DATA   HYF(8,2) /2.139793D0/
  69       DATA   HYF(9,2) /2.2210719D0/
  70       DATA   HYF(14,2)/6.663059D0/
  71       DATA   HYF(15,2)/5.657623D0/
  72       DATA   HYF(16,2)/6.345552D0/
  73       DATA   HYF(17,2)/2.522964D0/
  74       DATA ICALCN/0/
  75 C
  76 C  SETUP FOR DIPOLE CALCULATION
  77 C
  78       CALL CHRGE (P,Q2)
  79       DO 10 I = 1,NUMAT
  80          Q(I) = CORE(NAT(I)) - Q2(I)
  81    10 CONTINUE
  82       CALL GMETRY (GEO,COORD)
  83 C
  84       IF (ICALCN.NE.NUMCAL) THEN
  85          DO 20 I=2,107
  86    20    HYF(I,1)= 5.0832*DD(I)
  87          WTMOL=0.D0
  88          SUM=0.D0
  89          DO 30 I=1,NUMAT
  90             WTMOL=WTMOL+AMS(NAT(I))
  91    30    SUM=SUM+Q(I)
  92          CHARGD=(ABS(SUM).GT.0.5D0)
  93          ICALCN=NUMCAL
  94          KTYPE=1
  95          IF(ITYPE.EQ.4)KTYPE=2
  96       ENDIF
  97       IF(CHARGD)THEN
  98 C
  99 C   NEED TO RESET ION'S POSITION SO THAT THE CENTER OF MASS IS AT THE
 100 C   ORIGIN.
 101 C
 102 C$DOIT ASIS
 103          DO 40 I=1,3
 104    40    CENTER(I)=0.D0
 105          DO 50 I=1,3
 106 C$DOIT VBEST
 107             DO 50 J=1,NUMAT
 108    50    CENTER(I)=CENTER(I)+AMS(NAT(J))*COORD(I,J)
 109 C$DOIT ASIS
 110          DO 60 I=1,3
 111    60    CENTER(I)=CENTER(I)/WTMOL
 112          DO 70 I=1,3
 113 C$DOIT VBEST
 114             DO 70 J=1,NUMAT
 115    70    COORD(I,J)=COORD(I,J)-CENTER(I)
 116       ENDIF
 117 C$DOIT ASIS
 118       DO 80 I=1,4
 119 C$DOIT ASIS
 120          DO 80 J=1,3
 121    80 DIP(I,J)=0.0D00
 122 C$DOIT ASIS
 123       DO 100 I=1,NUMAT
 124          NI=NAT(I)
 125          IA=NFIRST(I)
 126          L=NLAST(I)-IA
 127 C$DOIT ASIS
 128          DO 90 J=1,L
 129             K=((IA+J)*(IA+J-1))/2+IA
 130    90    DIP(J,2)=DIP(J,2)-HYF(NI,KTYPE)*P(K)
 131 C$DOIT ASIS
 132          DO 100 J=1,3
 133   100 DIP(J,1)=DIP(J,1)+4.803D00*Q(I)*COORD(J,I)
 134 C$DOIT ASIS
 135       DO 110 J=1,3
 136   110 DIP(J,3)=DIP(J,2)+DIP(J,1)
 137 C$DOIT ASIS
 138       DO 120 J=1,3
 139   120 DIP(4,J)=SQRT(DIP(1,J)**2+DIP(2,J)**2+DIP(3,J)**2)
 140       DIPVEC(1)= -DIP(1,3)
 141       DIPVEC(2)= -DIP(2,3)
 142       DIPVEC(3)= -DIP(3,3)
 143 C      WRITE (6,60) ((DIP(I,J),I=1,4),J=1,3)
 144 C   60 FORMAT (3(4F10.3))
 145       RETURN
 146 C
 147       END