mopac606_nbo/src/rotate.f

   1       SUBROUTINE ROTATE (NI,NJ,XI,XJ,W,KR,E1B,E2A,ENUC,CUTOFF)
   2       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
   3 C***********************************************************************
   4 C
   5 C..IMPROVED SCALAR VERSION
   6 C..WRITTEN BY ERNEST R. DAVIDSON, INDIANA UNIVERSITY.
   7 C
   8 C
   9 C   ROTATE CALCULATES THE TWO-PARTICLE INTERACTIONS.
  10 C
  11 C   ON INPUT  NI     = ATOMIC NUMBER OF FIRST ATOM.
  12 C             NJ     = ATOMIC NUMBER OF SECOND ATOM.
  13 C             XI     = COORDINATE OF FIRST ATOM.
  14 C             XJ     = COORDINATE OF SECOND ATOM.
  15 C
  16 C ON OUTPUT W      = ARRAY OF TWO-ELECTRON REPULSION INTEGRALS.
  17 C           E1B,E2A= ARRAY OF ELECTRON-NUCLEAR ATTRACTION INTEGRALS,
  18 C                    E1B = ELECTRON ON ATOM NI ATTRACTING NUCLEUS OF NJ.
  19 C           ENUC   = NUCLEAR-NUCLEAR REPULSION TERM.
  20 C
  21 C
  22 C *** THIS ROUTINE COMPUTES THE REPULSION AND NUCLEAR ATTRACTION
  23 C     INTEGRALS OVER MOLECULAR-FRAME COORDINATES.  THE INTEGRALS OVER
  24 C     LOCAL FRAME COORDINATES ARE EVALUATED BY SUBROUTINE REPP AND
  25 C     STORED AS FOLLOWS (WHERE P-SIGMA = O,   AND P-PI = P AND P* )
  26 C     IN RI
  27 C     (SS/SS)=1,   (SO/SS)=2,   (OO/SS)=3,   (PP/SS)=4,   (SS/OS)=5,
  28 C     (SO/SO)=6,   (SP/SP)=7,   (OO/SO)=8,   (PP/SO)=9,   (PO/SP)=10,
  29 C     (SS/OO)=11,  (SS/PP)=12,  (SO/OO)=13,  (SO/PP)=14,  (SP/OP)=15,
  30 C     (OO/OO)=16,  (PP/OO)=17,  (OO/PP)=18,  (PP/PP)=19,  (PO/PO)=20,
  31 C     (PP/P*P*)=21,   (P*P/P*P)=22.
  32 C
  33 C***********************************************************************
  34       COMMON /NUMCAL/ NUMCAL
  35       SAVE ANALYT, ICALCN
  36       COMMON /MOLMEC/ HTYPE(4),NHCO(4,20),NNHCO,ITYPE
  37       CHARACTER*241 KEYWRD
  38       LOGICAL SI,SJ, ANALYT
  39       COMMON /NATORB/ NATORB(107)
  40       COMMON /TWOEL3/ F03(107)
  41       COMMON /ALPHA3/ ALP3(153)
  42       COMMON /ALPHA / ALP(107)
  43       COMMON /CORE  / TORE(107)
  44       COMMON /IDEAS / FN1(107,10),FN2(107,10),FN3(107,10)
  45       COMMON /ALPTM / ALPTM(30), EMUDTM(30)
  46       COMMON /ROTDUM/ CSS1,CSP1,CPPS1,CPPP1,CSS2,CSP2,CPPS2,CPPP2
  47       COMMON /ROTDU2/ X(3),Y(3),Z(3)
  48       COMMON /KEYWRD/ KEYWRD
  49       COMMON /W4G   / PREA(107), PREB(107), DA(107), DB(107)
  50       DIMENSION XI(3),XJ(3),W(100),E1B(10),E2A(10)
  51       DIMENSION RI(22),CCORE(4,2), BORON1(3,4), BORON2(3,4), BORON3(3,4)
  52       EQUIVALENCE (CCORE(1,1),CSS1)
  53       DATA ICALCN/0/
  54       DATA BORON1/  0.182613D0,  0.118587D0, -0.073280D0,
  55      1              0.412253D0, -0.149917D0,  0.000000D0,
  56      2              0.261751D0,  0.050275D0,  0.000000D0,
  57      3              0.359244D0,  0.074729D0,  0.000000D0/
  58       DATA BORON2/  6.D0,  6.D0,  5.D0,
  59      1             10.D0,  6.D0,  0.D0,
  60      2              8.D0,  5.D0,  0.D0,
  61      3              9.D0,  9.D0,  0.D0/
  62       DATA BORON3/  0.727592D0,  1.466639D0,  1.570975D0,
  63      1              0.832586D0,  1.186220D0,  0.000000D0,
  64      2              1.063995D0,  1.936492D0,  0.000000D0,
  65      3              0.819351D0,  1.574414D0,  0.000000D0/
  66 C
  67       IF(ICALCN.NE.NUMCAL)THEN
  68          ICALCN=NUMCAL
  69          ANALYT=(INDEX(KEYWRD,'ANALYT') .NE. 0)
  70          IF(ANALYT)THEN
  71             OPEN(UNIT=2,STATUS='SCRATCH',FORM='UNFORMATTED')
  72             REWIND 2
  73          ENDIF
  74       ENDIF
  75 C
  76       X(1)=XI(1)-XJ(1)
  77       X(2)=XI(2)-XJ(2)
  78       X(3)=XI(3)-XJ(3)
  79       RIJ=X(1)*X(1)+X(2)*X(2)+X(3)*X(3)
  80       IF (RIJ.LT.0.00002D0) THEN
  81 C
  82 C     SMALL RIJ CASE
  83 C
  84          DO 10 I=1,10
  85             E1B(I)=0.D0
  86             E2A(I)=0.D0
  87    10    CONTINUE
  88          W(KR)=0.D0
  89          ENUC=0.D0
  90 C
  91       ELSE IF (ITYPE.EQ.4) THEN
  92 C
  93 C     MINDO CASE
  94 C
  95          SUM=14.399D0/SQRT(RIJ+(7.1995D0/F03(NI)+7.1995D0/F03(NJ))**2)
  96          W(1)=SUM
  97          KR=KR+1
  98          DO 20 L=1,10
  99             E1B(L)=0.D0
 100             E2A(L)=0.D0
 101    20    CONTINUE
 102          E1B(1) = -SUM*TORE(NJ)
 103          E1B(3) = E1B(1)
 104          E1B(6) = E1B(1)
 105          E1B(10)= E1B(1)
 106          E2A(1) = -SUM*TORE(NI)
 107          E2A(3) = E2A(1)
 108          E2A(6) = E2A(1)
 109          E2A(10)= E2A(1)
 110          II = MAX(NI,NJ)
 111          NBOND = (II*(II-1))/2+NI+NJ-II
 112          RIJ = SQRT(RIJ)
 113          IF(NBOND.LT.154) THEN
 114             IF(NBOND.EQ.22 .OR. NBOND .EQ. 29) THEN
 115 C              NBOND = 22 IS C-H CASE
 116 C              NBOND = 29 IS N-H CASE
 117                SCALE=ALP3(NBOND)*EXP(-RIJ)
 118             ELSE
 119 C              NBOND < 154  IS NI < 18 AND NJ < 18 CASE
 120                SCALE=EXP(-ALP3(NBOND)*RIJ)
 121             ENDIF
 122          ELSE
 123 C              NBOND > 154 INVOLVES NI OR NJ > 18
 124             SCALE = 0
 125             IF(NATORB(NI).EQ.0) SCALE=      EXP(-ALP(NI)*RIJ)
 126             IF(NATORB(NJ).EQ.0) SCALE=SCALE+EXP(-ALP(NI)*RIJ)
 127          ENDIF
 128          IF (ABS(TORE(NI)).GT.20.D0 .AND. ABS(TORE(NJ)).GT.20.D0) THEN
 129             ENUC=0.D0
 130          ELSE IF (RIJ.LT.1.D0 .AND. NATORB(NI)*NATORB(NJ).EQ.0) THEN
 131             ENUC=0.D0
 132          ELSE
 133             ENUC = TORE(NI)*TORE(NJ)*SUM
 134      1       + ABS(TORE(NI)*TORE(NJ)*(14.399D0/RIJ-SUM)*SCALE)
 135          ENDIF
 136 C
 137 C     MNDO AND AM1 CASES
 138 C
 139 C *** THE REPULSION INTEGRALS OVER MOLECULAR FRAME (W) ARE STORED IN THE
 140 C     ORDER IN WHICH THEY WILL LATER BE USED.  IE.  (I,J/K,L) WHERE
 141 C     J.LE.I  AND  L.LE.K     AND L VARIES MOST RAPIDLY AND I LEAST
 142 C     RAPIDLY.  (ANTI-NORMAL COMPUTER STORAGE)
 143 C
 144       ELSE
 145 C
 146          RIJX = SQRT(RIJ)
 147          RIJ = MIN(RIJX,CUTOFF)
 148 C
 149 C *** COMPUTE INTEGRALS IN DIATOMIC FRAME
 150 C
 151          CALL REPP(NI,NJ,RIJ,RI,CCORE)
 152          IF(ANALYT)WRITE(2)(RI(I),I=1,22)
 153 C
 154          GAM = RI(1)
 155          A=1.D0/RIJX
 156          X(1) = X(1)*A
 157          X(2) = X(2)*A
 158          X(3) = X(3)*A
 159          IF (ABS(X(3)).GT.0.99999999D0) THEN
 160             X(3) = SIGN(1.D0,X(3))
 161             Y(1) = 0.D0
 162             Y(2) = 1.D0
 163             Y(3) = 0.D0
 164             Z(1) = 1.D0
 165             Z(2) = 0.D0
 166             Z(3) = 0.D0
 167          ELSE
 168             Z(3)=SQRT(1.D0-X(3)*X(3))
 169             A=1.D0/Z(3)
 170             Y(1)=-A*X(2)*SIGN(1.D0,X(1))
 171             Y(2)=ABS(A*X(1))
 172             Y(3)=0.D0
 173             Z(1)=-A*X(1)*X(3)
 174             Z(2)=-A*X(2)*X(3)
 175          ENDIF
 176          SI = (NATORB(NI).GT.1)
 177          SJ = (NATORB(NJ).GT.1)
 178          IF ( SI .OR. SJ) THEN
 179             XX11 = X(1)*X(1)
 180             XX21 = X(2)*X(1)
 181             XX22 = X(2)*X(2)
 182             XX31 = X(3)*X(1)
 183             XX32 = X(3)*X(2)
 184             XX33 = X(3)*X(3)
 185             YY11 = Y(1)*Y(1)
 186             YY21 = Y(2)*Y(1)
 187             YY22 = Y(2)*Y(2)
 188             ZZ11 = Z(1)*Z(1)
 189             ZZ21 = Z(2)*Z(1)
 190             ZZ22 = Z(2)*Z(2)
 191             ZZ31 = Z(3)*Z(1)
 192             ZZ32 = Z(3)*Z(2)
 193             ZZ33 = Z(3)*Z(3)
 194             YYZZ11 = YY11+ZZ11
 195             YYZZ21 = YY21+ZZ21
 196             YYZZ22 = YY22+ZZ22
 197             XY11 = 2.D0*X(1)*Y(1)
 198             XY21 =      X(1)*Y(2)+X(2)*Y(1)
 199             XY22 = 2.D0*X(2)*Y(2)
 200             XY31 =      X(3)*Y(1)
 201             XY32 =      X(3)*Y(2)
 202             XZ11 = 2.D0*X(1)*Z(1)
 203             XZ21 =      X(1)*Z(2)+X(2)*Z(1)
 204             XZ22 = 2.D0*X(2)*Z(2)
 205             XZ31 =      X(1)*Z(3)+X(3)*Z(1)
 206             XZ32 =      X(2)*Z(3)+X(3)*Z(2)
 207             XZ33 = 2.D0*X(3)*Z(3)
 208             YZ11 = 2.D0*Y(1)*Z(1)
 209             YZ21 =      Y(1)*Z(2)+Y(2)*Z(1)
 210             YZ22 = 2.D0*Y(2)*Z(2)
 211             YZ31 =      Y(1)*Z(3)
 212             YZ32 =      Y(2)*Z(3)
 213          ENDIF
 214 C
 215 C     (S S/S S)
 216          W(1)=RI(1)
 217          KI = 1
 218          IF (SJ) THEN
 219 C     (S S/PX S)
 220             W(2)=RI(5)*X(1)
 221 C     (S S/PX PX)
 222             W(3)=RI(11)*XX11+RI(12)*YYZZ11
 223 C     (S S/PY S)
 224             W(4)=RI(5)*X(2)
 225 C     (S S/PY PX)
 226             W(5)=RI(11)*XX21+RI(12)*YYZZ21
 227 C     (S S/PY PY)
 228             W(6)=RI(11)*XX22+RI(12)*YYZZ22
 229 C     (S S/PZ S)
 230             W(7)=RI(5)*X(3)
 231 C     (S S/PZ PX)
 232             W(8)=RI(11)*XX31+RI(12)*ZZ31
 233 C     (S S/PZ PY)
 234             W(9)=RI(11)*XX32+RI(12)*ZZ32
 235 C     (S S/PZ PZ)
 236             W(10)=RI(11)*XX33+RI(12)*ZZ33
 237             KI = 10
 238          ENDIF
 239 C
 240          IF (SI) THEN
 241 C     (PX S/S S)
 242             W(11)=RI(2)*X(1)
 243             IF (SJ) THEN
 244 C     (PX S/PX S)
 245                W(12)=RI(6)*XX11+RI(7)*YYZZ11
 246 C     (PX S/PX PX)
 247                W(13)=X(1)*(RI(13)*XX11+RI(14)*YYZZ11)
 248      1           +RI(15)*(Y(1)*XY11+Z(1)*XZ11)
 249 C     (PX S/PY S)
 250                W(14)=RI(6)*XX21+RI(7)*YYZZ21
 251 C     (PX S/PY PX)
 252                W(15)=X(1)*(RI(13)*XX21+RI(14)*YYZZ21)
 253      1           +RI(15)*(Y(1)*XY21+Z(1)*XZ21)
 254 C     (PX S/PY PY)
 255                W(16)=X(1)*(RI(13)*XX22+RI(14)*YYZZ22)
 256      1           +RI(15)*(Y(1)*XY22+Z(1)*XZ22)
 257 C     (PX S/PZ S)
 258                W(17)=RI(6)*XX31+RI(7)*ZZ31
 259 C     (PX S/PZ PX)
 260                W(18)=X(1)*(RI(13)*XX31+RI(14)*ZZ31)
 261      1           +RI(15)*(Y(1)*XY31+Z(1)*XZ31)
 262 C     (PX S/PZ PY)
 263                W(19)=X(1)*(RI(13)*XX32+RI(14)*ZZ32)
 264      1           +RI(15)*(Y(1)*XY32+Z(1)*XZ32)
 265 C     (PX S/PZ PZ)
 266                W(20)=X(1)*(RI(13)*XX33+RI(14)*ZZ33)
 267      1           +RI(15)*(          Z(1)*XZ33)
 268 C     (PX PX/S S)
 269                W(21)=RI(3)*XX11+RI(4)*YYZZ11
 270 C     (PX PX/PX S)
 271                W(22)=X(1)*(RI(8)*XX11+RI(9)*YYZZ11)
 272      1           +RI(10)*(Y(1)*XY11+Z(1)*XZ11)
 273 C     (PX PX/PX PX)
 274                W(23) =
 275      1     (RI(16)*XX11+RI(17)*YYZZ11)*XX11+RI(18)*XX11*YYZZ11
 276      2     +RI(19)*(YY11*YY11+ZZ11*ZZ11)
 277      3     +RI(20)*(XY11*XY11+XZ11*XZ11)
 278      4     +RI(21)*(YY11*ZZ11+ZZ11*YY11)
 279      5     +RI(22)*YZ11*YZ11
 280 C     (PX PX/PY S)
 281                W(24)=X(2)*(RI(8)*XX11+RI(9)*YYZZ11)
 282      1           +RI(10)*(Y(2)*XY11+Z(2)*XZ11)
 283 C     (PX PX/PY PX)
 284                W(25) =
 285      1     (RI(16)*XX11+RI(17)*YYZZ11)*XX21+RI(18)*XX11*YYZZ21
 286      2     +RI(19)*(YY11*YY21+ZZ11*ZZ21)
 287      3     +RI(20)*(XY11*XY21+XZ11*XZ21)
 288      4     +RI(21)*(YY11*ZZ21+ZZ11*YY21)
 289      5     +RI(22)*YZ11*YZ21
 290 C     (PX PX/PY PY)
 291                W(26) =
 292      1     (RI(16)*XX11+RI(17)*YYZZ11)*XX22+RI(18)*XX11*YYZZ22
 293      2     +RI(19)*(YY11*YY22+ZZ11*ZZ22)
 294      3     +RI(20)*(XY11*XY22+XZ11*XZ22)
 295      4     +RI(21)*(YY11*ZZ22+ZZ11*YY22)
 296      5     +RI(22)*YZ11*YZ22
 297 C     (PX PX/PZ S)
 298                W(27)=X(3)*(RI(8)*XX11+RI(9)*YYZZ11)
 299      1           +RI(10)*(         +Z(3)*XZ11)
 300 C     (PX PX/PZ PX)
 301                W(28) =
 302      1      (RI(16)*XX11+RI(17)*YYZZ11)*XX31
 303      2     +(RI(18)*XX11+RI(19)*ZZ11+RI(21)*YY11)*ZZ31
 304      3     +RI(20)*(XY11*XY31+XZ11*XZ31)
 305      4     +RI(22)*YZ11*YZ31
 306 C     (PX PX/PZ PY)
 307                W(29) =
 308      1      (RI(16)*XX11+RI(17)*YYZZ11)*XX32
 309      2     +(RI(18)*XX11+RI(19)*ZZ11+RI(21)*YY11)*ZZ32
 310      3     +RI(20)*(XY11*XY32+XZ11*XZ32)
 311      4     +RI(22)*YZ11*YZ32
 312 C     (PX PX/PZ PZ)
 313                W(30) =
 314      1      (RI(16)*XX11+RI(17)*YYZZ11)*XX33
 315      2     +(RI(18)*XX11+RI(19)*ZZ11+RI(21)*YY11)*ZZ33
 316      3     +RI(20)*XZ11*XZ33
 317 C     (PY S/S S)
 318                W(31)=RI(2)*X(2)
 319 C     (PY S/PX S)
 320                W(32)=RI(6)*XX21+RI(7)*YYZZ21
 321 C     (PY S/PX PX)
 322                W(33)=X(2)*(RI(13)*XX11+RI(14)*YYZZ11)
 323      1           +RI(15)*(Y(2)*XY11+Z(2)*XZ11)
 324 C     (PY S/PY S)
 325                W(34)=RI(6)*XX22+RI(7)*YYZZ22
 326 C     (PY S/PY PX)
 327                W(35)=X(2)*(RI(13)*XX21+RI(14)*YYZZ21)
 328      1           +RI(15)*(Y(2)*XY21+Z(2)*XZ21)
 329 C     (PY S/PY PY)
 330                W(36)=X(2)*(RI(13)*XX22+RI(14)*YYZZ22)
 331      1           +RI(15)*(Y(2)*XY22+Z(2)*XZ22)
 332 C     (PY S/PZ S)
 333                W(37)=RI(6)*XX32+RI(7)*ZZ32
 334 C     (PY S/PZ PX)
 335                W(38)=X(2)*(RI(13)*XX31+RI(14)*ZZ31)
 336      1           +RI(15)*(Y(2)*XY31+Z(2)*XZ31)
 337 C     (PY S/PZ PY)
 338                W(39)=X(2)*(RI(13)*XX32+RI(14)*ZZ32)
 339      1           +RI(15)*(Y(2)*XY32+Z(2)*XZ32)
 340 C     (PY S/PZ PZ)
 341                W(40)=X(2)*(RI(13)*XX33+RI(14)*ZZ33)
 342      1           +RI(15)*(         +Z(2)*XZ33)
 343 C     (PY PX/S S)
 344                W(41)=RI(3)*XX21+RI(4)*YYZZ21
 345 C     (PY PX/PX S)
 346                W(42)=X(1)*(RI(8)*XX21+RI(9)*YYZZ21)
 347      1           +RI(10)*(Y(1)*XY21+Z(1)*XZ21)
 348 C     (PY PX/PX PX)
 349                W(43) =
 350      1     (RI(16)*XX21+RI(17)*YYZZ21)*XX11+RI(18)*XX21*YYZZ11
 351      2     +RI(19)*(YY21*YY11+ZZ21*ZZ11)
 352      3     +RI(20)*(XY21*XY11+XZ21*XZ11)
 353      4     +RI(21)*(YY21*ZZ11+ZZ21*YY11)
 354      5     +RI(22)*YZ21*YZ11
 355 C     (PY PX/PY S)
 356                W(44)=X(2)*(RI(8)*XX21+RI(9)*YYZZ21)
 357      1           +RI(10)*(Y(2)*XY21+Z(2)*XZ21)
 358 C     (PY PX/PY PX)
 359                W(45) =
 360      1     (RI(16)*XX21+RI(17)*YYZZ21)*XX21+RI(18)*XX21*YYZZ21
 361      2     +RI(19)*(YY21*YY21+ZZ21*ZZ21)
 362      3     +RI(20)*(XY21*XY21+XZ21*XZ21)
 363      4     +RI(21)*(YY21*ZZ21+ZZ21*YY21)
 364      5     +RI(22)*YZ21*YZ21
 365 C     (PY PX/PY PY)
 366                W(46) =
 367      1     (RI(16)*XX21+RI(17)*YYZZ21)*XX22+RI(18)*XX21*YYZZ22
 368      2     +RI(19)*(YY21*YY22+ZZ21*ZZ22)
 369      3     +RI(20)*(XY21*XY22+XZ21*XZ22)
 370      4     +RI(21)*(YY21*ZZ22+ZZ21*YY22)
 371      5     +RI(22)*YZ21*YZ22
 372 C     (PY PX/PZ S)
 373                W(47)=X(3)*(RI(8)*XX21+RI(9)*YYZZ21)
 374      1           +RI(10)*(         +Z(3)*XZ21)
 375 C      (PY PX/PZ PX)
 376                W(48) =
 377      1     (RI(16)*XX21+RI(17)*YYZZ21)*XX31
 378      2     +(RI(18)*XX21+RI(19)*ZZ21+RI(21)*YY21)*ZZ31
 379      3     +RI(20)*(XY21*XY31+XZ21*XZ31)
 380      4     +RI(22)*YZ21*YZ31
 381 C      (PY PX/PZ PY)
 382                W(49) =
 383      1     (RI(16)*XX21+RI(17)*YYZZ21)*XX32
 384      2     +(RI(18)*XX21+RI(19)*ZZ21+RI(21)*YY21)*ZZ32
 385      3     +RI(20)*(XY21*XY32+XZ21*XZ32)
 386      4     +RI(22)*YZ21*YZ32
 387 C      (PY PX/PZ PZ)
 388                W(50) =
 389      1     (RI(16)*XX21+RI(17)*YYZZ21)*XX33
 390      2     +(RI(18)*XX21+RI(19)*ZZ21+RI(21)*YY21)*ZZ33
 391      3     +RI(20)*XZ21*XZ33
 392 C     (PY PY/S S)
 393                W(51)=RI(3)*XX22+RI(4)*YYZZ22
 394 C     (PY PY/PX S)
 395                W(52)=X(1)*(RI(8)*XX22+RI(9)*YYZZ22)
 396      1           +RI(10)*(Y(1)*XY22+Z(1)*XZ22)
 397 C      (PY PY/PX PX)
 398                W(53) =
 399      1     (RI(16)*XX22+RI(17)*YYZZ22)*XX11+RI(18)*XX22*YYZZ11
 400      2     +RI(19)*(YY22*YY11+ZZ22*ZZ11)
 401      3     +RI(20)*(XY22*XY11+XZ22*XZ11)
 402      4     +RI(21)*(YY22*ZZ11+ZZ22*YY11)
 403      5     +RI(22)*YZ22*YZ11
 404 C     (PY PY/PY S)
 405                W(54)=X(2)*(RI(8)*XX22+RI(9)*YYZZ22)
 406      1           +RI(10)*(Y(2)*XY22+Z(2)*XZ22)
 407 C      (PY PY/PY PX)
 408                W(55) =
 409      1     (RI(16)*XX22+RI(17)*YYZZ22)*XX21+RI(18)*XX22*YYZZ21
 410      2     +RI(19)*(YY22*YY21+ZZ22*ZZ21)
 411      3     +RI(20)*(XY22*XY21+XZ22*XZ21)
 412      4     +RI(21)*(YY22*ZZ21+ZZ22*YY21)
 413      5     +RI(22)*YZ22*YZ21
 414 C      (PY PY/PY PY)
 415                W(56) =
 416      1     (RI(16)*XX22+RI(17)*YYZZ22)*XX22+RI(18)*XX22*YYZZ22
 417      2     +RI(19)*(YY22*YY22+ZZ22*ZZ22)
 418      3     +RI(20)*(XY22*XY22+XZ22*XZ22)
 419      4     +RI(21)*(YY22*ZZ22+ZZ22*YY22)
 420      5     +RI(22)*YZ22*YZ22
 421 C     (PY PY/PZ S)
 422                W(57)=X(3)*(RI(8)*XX22+RI(9)*YYZZ22)
 423      1           +RI(10)*(         +Z(3)*XZ22)
 424 C      (PY PY/PZ PX)
 425                W(58) =
 426      1     (RI(16)*XX22+RI(17)*YYZZ22)*XX31
 427      2     +(RI(18)*XX22+RI(19)*ZZ22+RI(21)*YY22)*ZZ31
 428      3     +RI(20)*(XY22*XY31+XZ22*XZ31)
 429      4     +RI(22)*YZ22*YZ31
 430 C      (PY PY/PZ PY)
 431                W(59) =
 432      1     (RI(16)*XX22+RI(17)*YYZZ22)*XX32
 433      2     +(RI(18)*XX22+RI(19)*ZZ22+RI(21)*YY22)*ZZ32
 434      3     +RI(20)*(XY22*XY32+XZ22*XZ32)
 435      4     +RI(22)*YZ22*YZ32
 436 C      (PY PY/PZ PZ)
 437                W(60) =
 438      1     (RI(16)*XX22+RI(17)*YYZZ22)*XX33
 439      2     +(RI(18)*XX22+RI(19)*ZZ22+RI(21)*YY22)*ZZ33
 440      3     +RI(20)*XZ22*XZ33
 441 C     (PZ S/SS)
 442                W(61)=RI(2)*X(3)
 443 C     (PZ S/PX S)
 444                W(62)=RI(6)*XX31+RI(7)*ZZ31
 445 C     (PZ S/PX PX)
 446                W(63)=X(3)*(RI(13)*XX11+RI(14)*YYZZ11)
 447      1           +RI(15)*(         +Z(3)*XZ11)
 448 C     (PZ S/PY S)
 449                W(64)=RI(6)*XX32+RI(7)*ZZ32
 450 C     (PZ S/PY PX)
 451                W(65)=X(3)*(RI(13)*XX21+RI(14)*YYZZ21)
 452      1           +RI(15)*(         +Z(3)*XZ21)
 453 C     (PZ S/PY PY)
 454                W(66)=X(3)*(RI(13)*XX22+RI(14)*YYZZ22)
 455      1           +RI(15)*(         +Z(3)*XZ22)
 456 C     (PZ S/PZ S)
 457                W(67)=RI(6)*XX33+RI(7)*ZZ33
 458 C     (PZ S/PZ PX)
 459                W(68)=X(3)*(RI(13)*XX31+RI(14)*ZZ31)
 460      1           +RI(15)*(         +Z(3)*XZ31)
 461 C     (PZ S/PZ PY)
 462                W(69)=X(3)*(RI(13)*XX32+RI(14)*ZZ32)
 463      1           +RI(15)*(         +Z(3)*XZ32)
 464 C     (PZ S/PZ PZ)
 465                W(70)=X(3)*(RI(13)*XX33+RI(14)*ZZ33)
 466      1           +RI(15)*(         +Z(3)*XZ33)
 467 C     (PZ PX/S S)
 468                W(71)=RI(3)*XX31+RI(4)*ZZ31
 469 C     (PZ PX/PX S)
 470                W(72)=X(1)*(RI(8)*XX31+RI(9)*ZZ31)
 471      1           +RI(10)*(Y(1)*XY31+Z(1)*XZ31)
 472 C      (PZ PX/PX PX)
 473                W(73) =
 474      1     (RI(16)*XX31+RI(17)*ZZ31)*XX11+RI(18)*XX31*YYZZ11
 475      2     +RI(19)*ZZ31*ZZ11
 476      3     +RI(20)*(XY31*XY11+XZ31*XZ11)
 477      4     +RI(21)*ZZ31*YY11
 478      5     +RI(22)*YZ31*YZ11
 479 C     (PZ PX/PY S)
 480                W(74)=X(2)*(RI(8)*XX31+RI(9)*ZZ31)
 481      1           +RI(10)*(Y(2)*XY31+Z(2)*XZ31)
 482 C      (PZ PX/PY PX)
 483                W(75) =
 484      1     (RI(16)*XX31+RI(17)*ZZ31)*XX21+RI(18)*XX31*YYZZ21
 485      2     +RI(19)*ZZ31*ZZ21
 486      3     +RI(20)*(XY31*XY21+XZ31*XZ21)
 487      4     +RI(21)*ZZ31*YY21
 488      5     +RI(22)*YZ31*YZ21
 489 C      (PZ PX/PY PY)
 490                W(76) =
 491      1     (RI(16)*XX31+RI(17)*ZZ31)*XX22+RI(18)*XX31*YYZZ22
 492      2     +RI(19)*ZZ31*ZZ22
 493      3     +RI(20)*(XY31*XY22+XZ31*XZ22)
 494      4     +RI(21)*ZZ31*YY22
 495      5     +RI(22)*YZ31*YZ22
 496 C     (PZ PX/PZ S)
 497                W(77)=X(3)*(RI(8)*XX31+RI(9)*ZZ31)
 498      1           +RI(10)*(         +Z(3)*XZ31)
 499 C     (PZ PX/PZ PX)
 500                W(78) =
 501      1      (RI(16)*XX31+RI(17)*ZZ31)*XX31
 502      2     +(RI(18)*XX31+RI(19)*ZZ31)*ZZ31
 503      3     +RI(20)*(XY31*XY31+XZ31*XZ31)
 504      4     +RI(22)*YZ31*YZ31
 505 C      (PZ PX/PZ PY)
 506                W(79) =
 507      1      (RI(16)*XX31+RI(17)*ZZ31)*XX32
 508      2     +(RI(18)*XX31+RI(19)*ZZ31)*ZZ32
 509      3     +RI(20)*(XY31*XY32+XZ31*XZ32)
 510      4     +RI(22)*YZ31*YZ32
 511 C      (PZ PX/PZ PZ)
 512                W(80) =
 513      1      (RI(16)*XX31+RI(17)*ZZ31)*XX33
 514      2     +(RI(18)*XX31+RI(19)*ZZ31)*ZZ33
 515      3     +RI(20)*XZ31*XZ33
 516 C     (PZ PY/S S)
 517                W(81)=RI(3)*XX32+RI(4)*ZZ32
 518 C     (PZ PY/PX S)
 519                W(82)=X(1)*(RI(8)*XX32+RI(9)*ZZ32)
 520      1           +RI(10)*(Y(1)*XY32+Z(1)*XZ32)
 521 C      (PZ PY/PX PX)
 522                W(83) =
 523      1     (RI(16)*XX32+RI(17)*ZZ32)*XX11+RI(18)*XX32*YYZZ11
 524      2     +RI(19)*ZZ32*ZZ11
 525      3     +RI(20)*(XY32*XY11+XZ32*XZ11)
 526      4     +RI(21)*ZZ32*YY11
 527      5     +RI(22)*YZ32*YZ11
 528 C     (PZ PY/PY S)
 529                W(84)=X(2)*(RI(8)*XX32+RI(9)*ZZ32)
 530      1           +RI(10)*(Y(2)*XY32+Z(2)*XZ32)
 531 C      (PZ PY/PY PX)
 532                W(85) =
 533      1     (RI(16)*XX32+RI(17)*ZZ32)*XX21+RI(18)*XX32*YYZZ21
 534      2     +RI(19)*ZZ32*ZZ21
 535      3     +RI(20)*(XY32*XY21+XZ32*XZ21)
 536      4     +RI(21)*ZZ32*YY21
 537      5     +RI(22)*YZ32*YZ21
 538 C      (PZ PY/PY PY)
 539                W(86) =
 540      1     (RI(16)*XX32+RI(17)*ZZ32)*XX22+RI(18)*XX32*YYZZ22
 541      2     +RI(19)*ZZ32*ZZ22
 542      3     +RI(20)*(XY32*XY22+XZ32*XZ22)
 543      4     +RI(21)*ZZ32*YY22
 544      5     +RI(22)*YZ32*YZ22
 545 C     (PZ PY/PZ S)
 546                W(87)=X(3)*(RI(8)*XX32+RI(9)*ZZ32)
 547      1           +RI(10)*(         +Z(3)*XZ32)
 548 C      (PZ PY/PZ PX)
 549                W(88) =
 550      1      (RI(16)*XX32+RI(17)*ZZ32)*XX31
 551      2     +(RI(18)*XX32+RI(19)*ZZ32)*ZZ31
 552      3     +RI(20)*(XY32*XY31+XZ32*XZ31)
 553      4     +RI(22)*YZ32*YZ31
 554 C      (PZ PY/PZ PY)
 555                W(89) =
 556      1      (RI(16)*XX32+RI(17)*ZZ32)*XX32
 557      2     +(RI(18)*XX32+RI(19)*ZZ32)*ZZ32
 558      3     +RI(20)*(XY32*XY32+XZ32*XZ32)
 559      4     +RI(22)*YZ32*YZ32
 560 C       (PZ PY/PZ PZ)
 561                W(90) =
 562      1      (RI(16)*XX32+RI(17)*ZZ32)*XX33
 563      2     +(RI(18)*XX32+RI(19)*ZZ32)*ZZ33
 564      3     +RI(20)*XZ32*XZ33
 565 C     (PZ PZ/S S)
 566                W(91)=RI(3)*XX33+RI(4)*ZZ33
 567 C     (PZ PZ/PX S)
 568                W(92)=X(1)*(RI(8)*XX33+RI(9)*ZZ33)
 569      1           +RI(10)*(          Z(1)*XZ33)
 570 C       (PZ PZ/PX PX)
 571                W(93) =
 572      1     (RI(16)*XX33+RI(17)*ZZ33)*XX11+RI(18)*XX33*YYZZ11
 573      2     +RI(19)*ZZ33*ZZ11
 574      3     +RI(20)*XZ33*XZ11
 575      4     +RI(21)*ZZ33*YY11
 576 C     (PZ PZ/PY S)
 577                W(94)=X(2)*(RI(8)*XX33+RI(9)*ZZ33)
 578      1           +RI(10)*(         +Z(2)*XZ33)
 579 C       (PZ PZ/PY PX)
 580                W(95) =
 581      1     (RI(16)*XX33+RI(17)*ZZ33)*XX21+RI(18)*XX33*YYZZ21
 582      2     +RI(19)*ZZ33*ZZ21
 583      3     +RI(20)*XZ33*XZ21
 584      4     +RI(21)*ZZ33*YY21
 585 C       (PZ PZ/PY PY)
 586                W(96) =
 587      1     (RI(16)*XX33+RI(17)*ZZ33)*XX22+RI(18)*XX33*YYZZ22
 588      2     +RI(19)*ZZ33*ZZ22
 589      3     +RI(20)*XZ33*XZ22
 590      4     +RI(21)*ZZ33*YY22
 591 C     (PZ PZ/PZ S)
 592                W(97)=X(3)*(RI(8)*XX33+RI(9)*ZZ33)
 593      1           +RI(10)*(         +Z(3)*XZ33)
 594 C       (PZ PZ/PZ PX)
 595                W(98) =
 596      1      (RI(16)*XX33+RI(17)*ZZ33)*XX31
 597      2     +(RI(18)*XX33+RI(19)*ZZ33)*ZZ31
 598      3     +RI(20)*XZ33*XZ31
 599 C       (PZ PZ/PZ PY)
 600                W(99) =
 601      1      (RI(16)*XX33+RI(17)*ZZ33)*XX32
 602      2     +(RI(18)*XX33+RI(19)*ZZ33)*ZZ32
 603      3     +RI(20)*XZ33*XZ32
 604 C       (PZ PZ/PZ PZ)
 605                W(100) =
 606      1      (RI(16)*XX33+RI(17)*ZZ33)*XX33
 607      2     +(RI(18)*XX33+RI(19)*ZZ33)*ZZ33
 608      3     +RI(20)*XZ33*XZ33
 609                KI = 100
 610             ELSE
 611 C     (PX S/S S)
 612                W(2)=RI(2)*X(1)
 613 C     (PX PX/S S)
 614                W(3)=RI(3)*XX11+RI(4)*YYZZ11
 615 C     (PY S/S S)
 616                W(4)=RI(2)*X(2)
 617 C     (PY PX/S S)
 618                W(5)=RI(3)*XX21+RI(4)*YYZZ21
 619 C     (PY PY/S S)
 620                W(6)=RI(3)*XX22+RI(4)*YYZZ22
 621 C     (PZ S/SS)
 622                W(7)=RI(2)*X(3)
 623 C     (PZ PX/S S)
 624                W(8)=RI(3)*XX31+RI(4)*ZZ31
 625 C     (PZ PY/S S)
 626                W(9)=RI(3)*XX32+RI(4)*ZZ32
 627 C     (PZ PZ/S S)
 628                W(10)=RI(3)*XX33+RI(4)*ZZ33
 629                KI = 10
 630             END IF
 631          END IF
 632 C
 633 C *** NOW ROTATE THE NUCLEAR ATTRACTION INTEGRALS.
 634 C *** THE STORAGE OF THE NUCLEAR ATTRACTION INTEGRALS  CORE(KL/IJ) IS
 635 C     (SS/)=1,   (SO/)=2,   (OO/)=3,   (PP/)=4
 636 C
 637          E1B(1)=-CSS1
 638          IF(NATORB(NI).EQ.4) THEN
 639             E1B(2) = -CSP1 *X(1)
 640             E1B(3) = -CPPS1*XX11-CPPP1*YYZZ11
 641             E1B(4) = -CSP1 *X(2)
 642             E1B(5) = -CPPS1*XX21-CPPP1*YYZZ21
 643             E1B(6) = -CPPS1*XX22-CPPP1*YYZZ22
 644             E1B(7) = -CSP1 *X(3)
 645             E1B(8) = -CPPS1*XX31-CPPP1*ZZ31
 646             E1B(9) = -CPPS1*XX32-CPPP1*ZZ32
 647             E1B(10)= -CPPS1*XX33-CPPP1*ZZ33
 648          END IF
 649          E2A(1)=-CSS2
 650          IF(NATORB(NJ).EQ.4) THEN
 651             E2A(2) = -CSP2 *X(1)
 652             E2A(3) = -CPPS2*XX11-CPPP2*YYZZ11
 653             E2A(4) = -CSP2 *X(2)
 654             E2A(5) = -CPPS2*XX21-CPPP2*YYZZ21
 655             E2A(6) = -CPPS2*XX22-CPPP2*YYZZ22
 656             E2A(7) = -CSP2 *X(3)
 657             E2A(8) = -CPPS2*XX31-CPPP2*ZZ31
 658             E2A(9) = -CPPS2*XX32-CPPP2*ZZ32
 659             E2A(10)= -CPPS2*XX33-CPPP2*ZZ33
 660          END IF
 661          IF(ABS(TORE(NI)).GT.20.D0.AND.ABS(TORE(NJ)).GT.20.D0) THEN
 662 C SPARKLE-SPARKLE INTERACTION
 663             ENUC=0.D0
 664             RETURN
 665          ELSEIF (RIJ.LT.1.D0.AND.NATORB(NI)*NATORB(NJ).EQ.0) THEN
 666             ENUC=0.D0
 667             RETURN
 668          ENDIF
 669          SCALE = EXP(-ALP(NI)*RIJ)+EXP(-ALP(NJ)*RIJ)
 670 C
 671          IF (NI.EQ.24.AND.NJ.EQ.24) THEN
 672             SCALE = EXP(-ALPTM(NI)*RIJ)+EXP(-ALPTM(NJ)*RIJ)
 673          ENDIF
 674 C
 675          NT=NI+NJ
 676          IF(NT.EQ.8.OR.NT.EQ.9) THEN
 677             IF(NI.EQ.7.OR.NI.EQ.8) SCALE=SCALE+(RIJ-1.D0)*EXP(-ALP(NI)*R
 678      1IJ)
 679             IF(NJ.EQ.7.OR.NJ.EQ.8) SCALE=SCALE+(RIJ-1.D0)*EXP(-ALP(NJ)*R
 680      1IJ)
 681          ENDIF
 682          ENUC = TORE(NI)*TORE(NJ)*GAM
 683          SCALE=ABS(SCALE*ENUC)
 684          IF(ITYPE.EQ.2.AND.(NI.EQ.5.OR.NJ.EQ.5))THEN
 685 C
 686 C   LOAD IN AM1 BORON GAUSSIANS
 687 C
 688             NK=NI+NJ-5
 689 C   NK IS THE ATOMIC NUMBER OF THE NON-BORON ATOM
 690             NL=1
 691             IF(NK.EQ.1)NL=2
 692             IF(NK.EQ.6)NL=3
 693             IF(NK.EQ.9.OR.NK.EQ.17.OR.NK.EQ.35.OR.NK.EQ.53)NL=4
 694             DO 30 I=1,3
 695                FN1(5,I)=BORON1(I,NL)
 696                FN2(5,I)=BORON2(I,NL)
 697    30       FN3(5,I)=BORON3(I,NL)
 698          ENDIF
 699          IF(ITYPE.EQ.2.OR.ITYPE.EQ.3.OR.ITYPE.EQ.5) THEN
 700             DO 40 IG=1,10
 701                IF(ABS(FN1(NI,IG)).GT.0.D0) THEN
 702                   AX = FN2(NI,IG)*(RIJ-FN3(NI,IG))**2
 703                   IF(AX .LE. 25.D0) THEN
 704                      SCALE=SCALE +TORE(NI)*TORE(NJ)/RIJ*FN1(NI,IG)*EXP(-
 705      1AX)
 706                   ENDIF
 707                ENDIF
 708                IF(ABS(FN1(NJ,IG)).GT.0.D0) THEN
 709                   AX = FN2(NJ,IG)*(RIJ-FN3(NJ,IG))**2
 710                   IF(AX .LE. 25.D0) THEN
 711                      SCALE=SCALE +TORE(NI)*TORE(NJ)/RIJ*FN1(NJ,IG)*EXP(-
 712      1AX)
 713                   ENDIF
 714                ENDIF
 715    40       CONTINUE
 716          ENDIF
 717          IF (ITYPE.EQ.5.OR.ITYPE.EQ.6) THEN
 718 C            WRITE(6,*) ITYPE,NI,PREA(NI), PREB(NI), DA(NI), DB(NI)
 719 C            WRITE(6,*) ITYPE,NJ,PREA(NJ), PREB(NJ), DA(NJ), DB(NJ)
 720 C
 721 C THE PDDG FUNCTION HAS BEEN ADDED FOR THE PDDG METHODS
 722 C
 723             QCORR = 0.0D0
 724             ZAF=TORE(NI)/(TORE(NI)+TORE(NJ))
 725             ZBF=TORE(NJ)/(TORE(NI)+TORE(NJ))
 726             QCORR =
 727      1 (ZAF*PREA(NI)+ZBF*PREA(NJ))*EXP(-10.0D0*(RIJ-DA(NI)-DA(NJ))**2)+
 728      2 (ZAF*PREA(NI)+ZBF*PREB(NJ))*EXP(-10.0D0*(RIJ-DA(NI)-DB(NJ))**2)+
 729      3 (ZAF*PREB(NI)+ZBF*PREA(NJ))*EXP(-10.0D0*(RIJ-DB(NI)-DA(NJ))**2)+
 730      4 (ZAF*PREB(NI)+ZBF*PREB(NJ))*EXP(-10.0D0*(RIJ-DB(NI)-DB(NJ))**2)
 731          ELSE
 732             QCORR = 0.0D0
 733          ENDIF
 734 C         WRITE(6,*) 'QCORR:: ', QCORR
 735          ENUC=ENUC+SCALE+QCORR
 736 C
 737          IF(NATORB(NI)*NATORB(NJ).EQ.0)KI=0
 738          KR=KR+KI
 739 C
 740 C
 741       ENDIF
 742       RETURN
 743       END