mopac606_nbo/src/repp.f

   1       SUBROUTINE REPP(NI,NJ,RIJ,RI,CORE)
   2 C***********************************************************************
   3 C
   4 C..VECTOR VERSION WRITTEN BY ERNEST R. DAVIDSON, INDIANA UNIVERSITY
   5 C
   6 C
   7 C  REPP CALCULATES THE TWO-ELECTRON REPULSION INTEGRALS AND THE
   8 C       NUCLEAR ATTRACTION INTEGRALS.
   9 C
  10 C     ON INPUT RIJ     = INTERATOMIC DISTANCE
  11 C              NI      = ATOM NUMBER OF FIRST ATOM
  12 C              NJ      = ATOM NUMBER OF SECOND ATOM
  13 C    (REF)     ADD     = ARRAY OF GAMMA, OR TWO-ELECTRON ONE-CENTER,
  14 C                        INTEGRALS.
  15 C    (REF)     TORE    = ARRAY OF NUCLEAR CHARGES OF THE ELEMENTS
  16 C    (REF)     DD      = ARRAY OF DIPOLE CHARGE SEPARATIONS
  17 C    (REF)     QQ      = ARRAY OF QUADRUPOLE CHARGE SEPARATIONS
  18 C
  19 C     THE COMMON BLOCKS ARE INITIALIZED IN BLOCK-DATA, AND NEVER CHANGED
  20 C
  21 C    ON OUTPUT RI      = ARRAY OF TWO-ELECTRON REPULSION INTEGRALS
  22 C              CORE    = 4 X 2 ARRAY OF ELECTRON-CORE ATTRACTION
  23 C                        INTEGRALS
  24 C
  25 C
  26 C *** THIS ROUTINE COMPUTES THE TWO-CENTRE REPULSION INTEGRALS AND THE
  27 C *** NUCLEAR ATTRACTION INTEGRALS.
  28 C *** THE TWO-CENTRE REPULSION INTEGRALS (OVER LOCAL COORDINATES) ARE
  29 C *** STORED AS FOLLOWS (WHERE P-SIGMA = O,  AND P-PI = P AND P* )
  30 C     (SS/SS)=1,   (SO/SS)=2,   (OO/SS)=3,   (PP/SS)=4,   (SS/OS)=5,
  31 C     (SO/SO)=6,   (SP/SP)=7,   (OO/SO)=8,   (PP/SO)=9,   (PO/SP)=10,
  32 C     (SS/OO)=11,  (SS/PP)=12,  (SO/OO)=13,  (SO/PP)=14,  (SP/OP)=15,
  33 C     (OO/OO)=16,  (PP/OO)=17,  (OO/PP)=18,  (PP/PP)=19,  (PO/PO)=20,
  34 C     (PP/P*P*)=21,   (P*P/P*P)=22.
  35 C *** THE STORAGE OF THE NUCLEAR ATTRACTION INTEGRALS  CORE(KL/IJ) IS
  36 C     (SS/)=1,   (SO/)=2,   (OO/)=3,   (PP/)=4
  37 C     WHERE IJ=1 IF THE ORBITALS CENTRED ON ATOM I,  =2 IF ON ATOM J.
  38 C *** NI AND NJ ARE THE ATOMIC NUMBERS OF THE TWO ELEMENTS.
  39 C
  40 C***********************************************************************
  41       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  42       LOGICAL SI,SJ
  43       COMMON /MULTIP/ DD(107),QQ(107),ADD(107,3)
  44       COMMON /CORE/ TORE(107)
  45       COMMON /NATORB/ NATORB(107)
  46       DIMENSION RI(22),CORE(4,2)
  47       DIMENSION ARG(72),SQR(72)
  48       DATA  TD/2.D00/
  49       DATA PP/0.5D00/
  50       DATA A0/0.529167D0/ ,EV/27.21D0/, EV1/13.605D0/, EV2/6.8025D0/,
  51      1 EV3/3.40125D0/, EV4/1.700625D0/
  52 C
  53 C     ATOMIC UNITS ARE USED IN THE CALCULATION,
  54 C     FINAL RESULTS ARE CONVERTED TO EV
  55 C
  56       R=RIJ/A0
  57 C
  58       SI = (NATORB(NI).GE.3)
  59       SJ = (NATORB(NJ).GE.3)
  60 C
  61       IF ((.NOT.SI) .AND. (.NOT.SJ)) THEN
  62 C
  63 C     HYDROGEN - HYDROGEN  (SS/SS)
  64 C
  65          AEE = PP/ADD(NI,1) + PP/ADD(NJ,1)
  66          AEE = AEE * AEE
  67          RI(1) = EV/SQRT(R*R+AEE)
  68          CORE(1,1 )= TORE(NJ)*RI(1)
  69          CORE(1,2) = TORE(NI)*RI(1)
  70 C
  71       ELSE IF (SI .AND. (.NOT.SJ)) THEN
  72 C
  73 C     HEAVY ATOM - HYDROGEN
  74 C
  75          AEE = PP/ADD(NI,1) + PP/ADD(NJ,1)
  76          AEE = AEE * AEE
  77          DA=DD(NI)
  78          QA=QQ(NI) * TD
  79          ADE = PP/ADD(NI,2) + PP/ADD(NJ,1)
  80          ADE = ADE * ADE
  81          AQE = PP/ADD(NI,3) + PP/ADD(NJ,1)
  82          AQE = AQE * AQE
  83          RSQ = R*R
  84          ARG(1) = RSQ + AEE
  85          XXX = R+DA
  86          ARG(2) = XXX*XXX + ADE
  87          XXX = R-DA
  88          ARG(3) = XXX*XXX + ADE
  89          XXX = R+QA
  90          ARG(4) = XXX*XXX + AQE
  91          XXX = R-QA
  92          ARG(5) = XXX*XXX + AQE
  93          ARG(6) = RSQ + AQE
  94          ARG(7) = ARG(6) + QA*QA
  95 C$DOIT ASIS
  96          DO 10 I = 1,7
  97             SQR(I) = SQRT(ARG(I))
  98    10    CONTINUE
  99          EE = EV/SQR(1)
 100          RI(1) = EE
 101          RI(2) = EV1/SQR(2) - EV1/SQR(3)
 102          RI(3) = EE + EV2/SQR(4) + EV2/SQR(5) - EV1/SQR(6)
 103          RI(4) = EE + EV1/SQR(7) - EV1/SQR(6)
 104          CORE(1,1) = TORE(NJ)*RI(1)
 105          CORE(1,2) = TORE(NI)*RI(1)
 106          CORE(2,1) = TORE(NJ)*RI(2)
 107          CORE(3,1) = TORE(NJ)*RI(3)
 108          CORE(4,1) = TORE(NJ)*RI(4)
 109 C
 110       ELSE IF ((.NOT.SI).AND.SJ) THEN
 111 C
 112 C     HYDROGEN - HEAVY ATOM
 113 C
 114          AEE = PP/ADD(NI,1) + PP/ADD(NJ,1)
 115          AEE = AEE * AEE
 116          DB=DD(NJ)
 117          QB=QQ(NJ) * TD
 118          AED = PP/ADD(NI,1) + PP/ADD(NJ,2)
 119          AED = AED * AED
 120          AEQ = PP/ADD(NI,1) + PP/ADD(NJ,3)
 121          AEQ = AEQ * AEQ
 122          RSQ = R*R
 123          ARG(1) = RSQ + AEE
 124          XXX = R-DB
 125          ARG(2) = XXX*XXX + AED
 126          XXX = R+DB
 127          ARG(3) = XXX*XXX + AED
 128          XXX = R-QB
 129          ARG(4) = XXX*XXX + AEQ
 130          XXX = R+QB
 131          ARG(5) = XXX*XXX + AEQ
 132          ARG(6) = RSQ + AEQ
 133          ARG(7) = ARG(6) + QB*QB
 134 C$DOIT ASIS
 135          DO 20 I = 1,7
 136             SQR(I) = SQRT(ARG(I))
 137    20    CONTINUE
 138          EE = EV/SQR(1)
 139          RI(1) = EE
 140          RI(5) = EV1/SQR(2)  - EV1/SQR(3)
 141          RI(11) = EE + EV2/SQR(4) + EV2/SQR(5) - EV1/SQR(6)
 142          RI(12) = EE + EV1/SQR(7) - EV1/SQR(6)
 143          CORE(1,1) = TORE(NJ)*RI(1)
 144          CORE(1,2) = TORE(NI)*RI(1)
 145          CORE(2,2) = TORE(NI)*RI(5)
 146          CORE(3,2) = TORE(NI)*RI(11)
 147          CORE(4,2) = TORE(NI)*RI(12)
 148 C
 149       ELSE
 150 C
 151 C     HEAVY ATOM - HEAVY ATOM
 152 C
 153 C     DEFINE CHARGE SEPARATIONS.
 154          DA=DD(NI)
 155          DB=DD(NJ)
 156          QA=QQ(NI) * TD
 157          QB=QQ(NJ) * TD
 158 C
 159          AEE = PP/ADD(NI,1) + PP/ADD(NJ,1)
 160          AEE = AEE * AEE
 161 C
 162          ADE = PP/ADD(NI,2) + PP/ADD(NJ,1)
 163          ADE = ADE * ADE
 164          AQE = PP/ADD(NI,3) + PP/ADD(NJ,1)
 165          AQE = AQE * AQE
 166          AED = PP/ADD(NI,1) + PP/ADD(NJ,2)
 167          AED = AED * AED
 168          AEQ = PP/ADD(NI,1) + PP/ADD(NJ,3)
 169          AEQ = AEQ * AEQ
 170          AXX = PP/ADD(NI,2) + PP/ADD(NJ,2)
 171          AXX = AXX * AXX
 172          ADQ = PP/ADD(NI,2) + PP/ADD(NJ,3)
 173          ADQ = ADQ * ADQ
 174          AQD = PP/ADD(NI,3) + PP/ADD(NJ,2)
 175          AQD = AQD * AQD
 176          AQQ = PP/ADD(NI,3) + PP/ADD(NJ,3)
 177          AQQ = AQQ * AQQ
 178          RSQ = R * R
 179          ARG(1) = RSQ + AEE
 180          XXX = R + DA
 181          ARG(2) = XXX * XXX + ADE
 182          XXX = R - DA
 183          ARG(3) = XXX*XXX + ADE
 184          XXX = R - QA
 185          ARG(4) = XXX*XXX + AQE
 186          XXX = R + QA
 187          ARG(5) = XXX*XXX + AQE
 188          ARG(6) = RSQ + AQE
 189          ARG(7) = ARG(6) + QA*QA
 190          XXX = R-DB
 191          ARG(8) = XXX*XXX + AED
 192          XXX = R+DB
 193          ARG(9) = XXX*XXX + AED
 194          XXX = R - QB
 195          ARG(10) = XXX*XXX + AEQ
 196          XXX = R + QB
 197          ARG(11) = XXX*XXX + AEQ
 198          ARG(12) = RSQ + AEQ
 199          ARG(13) = ARG(12) + QB*QB
 200          XXX = DA-DB
 201          ARG(14) = RSQ + AXX + XXX*XXX
 202          XXX = DA+DB
 203          ARG(15) = RSQ + AXX + XXX*XXX
 204          XXX = R + DA - DB
 205          ARG(16) = XXX*XXX + AXX
 206          XXX = R - DA + DB
 207          ARG(17) = XXX*XXX + AXX
 208          XXX = R - DA - DB
 209          ARG(18) = XXX*XXX + AXX
 210          XXX = R + DA + DB
 211          ARG(19) = XXX*XXX + AXX
 212          XXX = R + DA
 213          ARG(20) = XXX*XXX + ADQ
 214          ARG(21) = ARG(20) + QB*QB
 215          XXX = R - DA
 216          ARG(22) = XXX*XXX + ADQ
 217          ARG(23) = ARG(22) + QB*QB
 218          XXX = R - DB
 219          ARG(24) = XXX*XXX + AQD
 220          ARG(25) = ARG(24) + QA*QA
 221          XXX = R + DB
 222          ARG(26) = XXX*XXX + AQD
 223          ARG(27) = ARG(26) + QA*QA
 224          XXX = R + DA - QB
 225          ARG(28) = XXX*XXX + ADQ
 226          XXX = R - DA - QB
 227          ARG(29) = XXX*XXX + ADQ
 228          XXX = R + DA + QB
 229          ARG(30) = XXX*XXX + ADQ
 230          XXX = R - DA + QB
 231          ARG(31) = XXX*XXX + ADQ
 232          XXX = R + QA - DB
 233          ARG(32) = XXX*XXX + AQD
 234          XXX = R + QA + DB
 235          ARG(33) = XXX*XXX + AQD
 236          XXX = R - QA - DB
 237          ARG(34) = XXX*XXX + AQD
 238          XXX = R - QA + DB
 239          ARG(35) = XXX*XXX + AQD
 240          ARG(36) = RSQ + AQQ
 241          XXX = QA - QB
 242          ARG(37) = ARG(36) + XXX*XXX
 243          XXX = QA + QB
 244          ARG(38) = ARG(36) + XXX*XXX
 245          ARG(39) = ARG(36) + QA*QA
 246          ARG(40) = ARG(36) + QB*QB
 247          ARG(41) = ARG(39) + QB*QB
 248          XXX = R - QB
 249          ARG(42) = XXX*XXX + AQQ
 250          ARG(43) = ARG(42) + QA*QA
 251          XXX = R + QB
 252          ARG(44) = XXX*XXX + AQQ
 253          ARG(45) = ARG(44) + QA*QA
 254          XXX = R + QA
 255          ARG(46) = XXX*XXX + AQQ
 256          ARG(47) = ARG(46) + QB*QB
 257          XXX = R - QA
 258          ARG(48) = XXX*XXX + AQQ
 259          ARG(49) = ARG(48) + QB*QB
 260          XXX = R + QA - QB
 261          ARG(50) = XXX*XXX + AQQ
 262          XXX = R + QA + QB
 263          ARG(51) = XXX*XXX + AQQ
 264          XXX = R - QA - QB
 265          ARG(52) = XXX*XXX + AQQ
 266          XXX = R - QA + QB
 267          ARG(53) = XXX*XXX + AQQ
 268          QA=QQ(NI)
 269          QB=QQ(NJ)
 270          XXX = DA - QB
 271          XXX = XXX*XXX
 272          YYY = R - QB
 273          YYY = YYY*YYY
 274          ZZZ = DA + QB
 275          ZZZ = ZZZ*ZZZ
 276          WWW = R + QB
 277          WWW = WWW*WWW
 278          ARG(54) = XXX + YYY + ADQ
 279          ARG(55) = XXX + WWW + ADQ
 280          ARG(56) = ZZZ + YYY + ADQ
 281          ARG(57) = ZZZ + WWW + ADQ
 282          XXX = QA - DB
 283          XXX = XXX*XXX
 284          YYY = QA + DB
 285          YYY = YYY*YYY
 286          ZZZ = R + QA
 287          ZZZ = ZZZ*ZZZ
 288          WWW = R - QA
 289          WWW = WWW*WWW
 290          ARG(58) = ZZZ + XXX + AQD
 291          ARG(59) = WWW + XXX + AQD
 292          ARG(60) = ZZZ + YYY + AQD
 293          ARG(61) = WWW + YYY + AQD
 294          XXX = QA - QB
 295          XXX = XXX*XXX
 296          ARG(62) = ARG(36) + TD*XXX
 297          YYY = QA + QB
 298          YYY = YYY*YYY
 299          ARG(63) = ARG(36) + TD*YYY
 300          ARG(64) = ARG(36) + TD*(QA*QA+QB*QB)
 301          ZZZ = R + QA - QB
 302          ZZZ = ZZZ*ZZZ
 303          ARG(65) = ZZZ + XXX + AQQ
 304          ARG(66) = ZZZ + YYY + AQQ
 305          ZZZ = R + QA + QB
 306          ZZZ = ZZZ*ZZZ
 307          ARG(67) = ZZZ + XXX + AQQ
 308          ARG(68) = ZZZ + YYY + AQQ
 309          ZZZ = R - QA - QB
 310          ZZZ = ZZZ*ZZZ
 311          ARG(69) = ZZZ + XXX + AQQ
 312          ARG(70) = ZZZ + YYY + AQQ
 313          ZZZ = R - QA + QB
 314          ZZZ = ZZZ*ZZZ
 315          ARG(71) = ZZZ + XXX + AQQ
 316          ARG(72) = ZZZ + YYY + AQQ
 317          DO 30 I = 1,72
 318             SQR(I) = SQRT(ARG(I))
 319    30    CONTINUE
 320          EE = EV/SQR(1)
 321          DZE = -EV1/SQR(2) + EV1/SQR(3)
 322          QZZE = EV2/SQR(4) + EV2/SQR(5) - EV1/SQR(6)
 323          QXXE = EV1/SQR(7) - EV1/SQR(6)
 324          EDZ = - EV1/SQR(8) + EV1/SQR(9)
 325          EQZZ  = EV2/SQR(10) + EV2/SQR(11) - EV1/SQR(12)
 326          EQXX  = EV1/SQR(13) - EV1/SQR(12)
 327          DXDX  = EV1/SQR(14) - EV1/SQR(15)
 328          DZDZ  = EV2/SQR(16) + EV2/SQR(17) - EV2/SQR(18) - EV2/SQR(19)
 329          DZQXX =  EV2/SQR(20) - EV2/SQR(21) - EV2/SQR(22) + EV2/SQR(23)
 330          QXXDZ =  EV2/SQR(24) - EV2/SQR(25) - EV2/SQR(26) + EV2/SQR(27)
 331          DZQZZ = -EV3/SQR(28) + EV3/SQR(29) - EV3/SQR(30) + EV3/SQR(31)
 332      1       - EV2/SQR(22) + EV2/SQR(20)
 333          QZZDZ = -EV3/SQR(32) + EV3/SQR(33) - EV3/SQR(34) + EV3/SQR(35)
 334      1       + EV2/SQR(24) - EV2/SQR(26)
 335          QXXQXX = EV3/SQR(37) + EV3/SQR(38) - EV2/SQR(39) - EV2/SQR(40)
 336      1       + EV2/SQR(36)
 337          QXXQYY = EV2/SQR(41) - EV2/SQR(39) - EV2/SQR(40) + EV2/SQR(36)
 338          QXXQZZ = EV3/SQR(43) + EV3/SQR(45) - EV3/SQR(42) - EV3/SQR(44)
 339      1       - EV2/SQR(39) + EV2/SQR(36)
 340          QZZQXX = EV3/SQR(47) + EV3/SQR(49) - EV3/SQR(46) - EV3/SQR(48)
 341      1       - EV2/SQR(40) + EV2/SQR(36)
 342          QZZQZZ = EV4/SQR(50) + EV4/SQR(51) + EV4/SQR(52) + EV4/SQR(53)
 343      1       - EV3/SQR(48) - EV3/SQR(46) - EV3/SQR(42) - EV3/SQR(44)
 344      2       + EV2/SQR(36)
 345          DXQXZ = -EV2/SQR(54) + EV2/SQR(55) + EV2/SQR(56) - EV2/SQR(57)
 346          QXZDX = -EV2/SQR(58) + EV2/SQR(59) + EV2/SQR(60) - EV2/SQR(61)
 347          QXZQXZ = EV3/SQR(65) - EV3/SQR(67) - EV3/SQR(69) + EV3/SQR(71)
 348      1       - EV3/SQR(66) + EV3/SQR(68) + EV3/SQR(70) - EV3/SQR(72)
 349          RI(1) = EE
 350          RI(2) = -DZE
 351          RI(3) = EE + QZZE
 352          RI(4) = EE + QXXE
 353          RI(5) = -EDZ
 354          RI(6) = DZDZ
 355          RI(7) = DXDX
 356          RI(8) = -EDZ -QZZDZ
 357          RI(9) = -EDZ -QXXDZ
 358          RI(10) = -QXZDX
 359          RI(11) =  EE + EQZZ
 360          RI(12) =  EE + EQXX
 361          RI(13) = -DZE -DZQZZ
 362          RI(14) = -DZE -DZQXX
 363          RI(15) = -DXQXZ
 364          RI(16) = EE +EQZZ +QZZE +QZZQZZ
 365          RI(17) = EE +EQZZ +QXXE +QXXQZZ
 366          RI(18) = EE +EQXX +QZZE +QZZQXX
 367          RI(19) = EE +EQXX +QXXE +QXXQXX
 368          RI(20) = QXZQXZ
 369          RI(21) = EE +EQXX +QXXE +QXXQYY
 370          RI(22) = PP * (QXXQXX -QXXQYY)
 371 C
 372 C     CALCULATE CORE-ELECTRON ATTRACTIONS.
 373 C
 374          CORE(1,1) = TORE(NJ)*RI(1)
 375          CORE(2,1) = TORE(NJ)*RI(2)
 376          CORE(3,1) = TORE(NJ)*RI(3)
 377          CORE(4,1) = TORE(NJ)*RI(4)
 378          CORE(1,2) = TORE(NI)*RI(1)
 379          CORE(2,2) = TORE(NI)*RI(5)
 380          CORE(3,2) = TORE(NI)*RI(11)
 381          CORE(4,2) = TORE(NI)*RI(12)
 382 C
 383       END IF
 384 C
 385       RETURN
 386 C
 387       END