1C Copyright (c) 2003-2010 University of Florida 2C 3C This program is free software; you can redistribute it and/or modify 4C it under the terms of the GNU General Public License as published by 5C the Free Software Foundation; either version 2 of the License, or 6C (at your option) any later version. 7 8C This program is distributed in the hope that it will be useful, 9C but WITHOUT ANY WARRANTY; without even the implied warranty of 10C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11C GNU General Public License for more details. 12 13C The GNU General Public License is included in this distribution 14C in the file COPYRIGHT. 15 SUBROUTINE ERD__DERV_INT2D_TO_ABCD 16 + 17 + ( SHELLA,SHELLB,SHELLC,SHELLD, 18 + NGQP,NEXQ,NGQEXQ, 19 + NXYZA,NXYZB,NXYZC,NXYZD, 20 + INT2DX,INT2DY,INT2DZ, 21 + DIFFY,DIFFZ, 22 + TEMP1,TEMP2, 23 + SCALE, 24 + 25 + BATCH ) 26 + 27C------------------------------------------------------------------------ 28C OPERATION : ERD__DERV_INT2D_TO_ABCD 29C MODULE : ELECTRON REPULSION INTEGRALS DIRECT 30C MODULE-ID : ERD 31C SUBROUTINES : none 32C DESCRIPTION : This routine assembles the set of batches of cartesian 33C derivative eris [AB|CD], adding up the contributions 34C from all the derivative 2D ABCD integrals. 35C 36C The routine uses the reduced Rys multiplication scheme 37C as suggested in R.Lindh, U.Ryu and B.Liu, JCP 95, 5889. 38C This scheme reuses intermediate products between 39C 2DX and 2DY integrals. 40C 41C Due to the very computational intensive steps inside 42C the x,y,z loops, special sections of identical x,y,z 43C loop structures have been given for each # of roots 44C =< 9, thus saving considerable computing time over the 45C general case. 46C 47C For comments on how the x,y,z loop structures are 48C coded please refer to the general root case. 49C 50C 51C Input: 52C 53C SHELLx = shell types for individual csh 54C x = A,B,C,D 55C NGQP = # of gaussian quadrature points 56C (roots) 57C NEXQ = current # of exponent quadruplets 58C NGQEXQ = product of # of gaussian quadrature 59C points times exponent quadruplets 60C NXYZx = # of cartesian monomials for 61C x = A,B,C,D shells 62C INT2Dx = all current 2D ABCD derivative 63C integrals for each cartesian 64C component (x = X,Y,Z) 65C DIFFx = is true, if differentiation was 66C performed along the x=Y,Z direction 67C TEMP1(2) = scratch arrays holding intermediate 68C 2D ABCD derivative integral products 69C SCALE = the NGQEXQ scaling factors 70C 71C 72C Output: 73C 74C BATCH = batch of primitive cartesian 75C [AB|CD] derivative integrals 76C corresponding to all current 77C exponent quadruplets 78C 79C 80C AUTHOR : Norbert Flocke 81C------------------------------------------------------------------------ 82C 83C ...include files and declare variables. 84C 85C 86 IMPLICIT NONE 87 88 LOGICAL DIFFY,DIFFZ 89 90 INTEGER I,J,K,L,M,N,R 91 INTEGER NGQP,NEXQ,NGQEXQ 92 INTEGER NXYZA,NXYZB,NXYZC,NXYZD 93 INTEGER SHELLA,SHELLB,SHELLC,SHELLD 94 INTEGER XA,YA,ZA,XB,YB,ZB,XC,YC,ZC,XD,YD,ZD 95 INTEGER XAP,XBP,XCP,XDP 96 INTEGER YAMAX,YBMAX,YCMAX,YDMAX 97 98 DOUBLE PRECISION SUM 99 DOUBLE PRECISION ZERO 100 101 DOUBLE PRECISION SCALE (1:NGQEXQ) 102 DOUBLE PRECISION TEMP1 (1:NGQEXQ) 103 DOUBLE PRECISION TEMP2 (1:NGQEXQ) 104 105 DOUBLE PRECISION BATCH (1:NEXQ, 106 + 1:NXYZA,1:NXYZB,1:NXYZC,1:NXYZD) 107 108 DOUBLE PRECISION INT2DX (1:NGQEXQ, 109 + 0:SHELLA,0:SHELLB,0:SHELLC,0:SHELLD) 110 DOUBLE PRECISION INT2DY (1:NGQEXQ, 111 + 0:SHELLA,0:SHELLB,0:SHELLC,0:SHELLD) 112 DOUBLE PRECISION INT2DZ (1:NGQEXQ, 113 + 0:SHELLA,0:SHELLB,0:SHELLC,0:SHELLD) 114 115 PARAMETER (ZERO = 0.D0) 116C 117C 118C------------------------------------------------------------------------ 119C 120C 121C ...jump according to number of roots. 122C 123C 124 GOTO (1,2,3,4,5,6,7,8,9,10) MIN (NGQP,10) 125C 126C 127C ******************** 128C * # of roots = 1 * 129C ******************** 130C 131C 132 1 XDP = 0 133 DO 100 XD = SHELLD,0,-1 134 YDMAX = SHELLD - XD 135 XCP = 0 136 DO 102 XC = SHELLC,0,-1 137 YCMAX = SHELLC - XC 138 XBP = 0 139 DO 104 XB = SHELLB,0,-1 140 YBMAX = SHELLB - XB 141 XAP = 0 142 DO 106 XA = SHELLA,0,-1 143 YAMAX = SHELLA - XA 144 145 DO M = 1,NEXQ 146 TEMP1 (M) = SCALE (M) * INT2DX (M,XA,XB,XC,XD) 147 END DO 148 149 L = XDP 150 DO 110 YD = YDMAX,0,-1 151 L = L + 1 152 ZD = YDMAX - YD 153 K = XCP 154 DO 112 YC = YCMAX,0,-1 155 K = K + 1 156 ZC = YCMAX - YC 157 J = XBP 158 DO 114 YB = YBMAX,0,-1 159 J = J + 1 160 ZB = YBMAX - YB 161 I = XAP 162 DO 116 YA = YAMAX,0,-1 163 I = I + 1 164 ZA = YAMAX - YA 165 166 IF (.NOT.DIFFY .AND. YA+YB+YC+YD.EQ.0) THEN 167 DO M = 1,NEXQ 168 TEMP2 (M) = TEMP1 (M) 169 END DO 170 ELSE 171 DO M = 1,NEXQ 172 TEMP2 (M) = TEMP1 (M) 173 + * INT2DY (M,YA,YB,YC,YD) 174 END DO 175 END IF 176 177 IF (.NOT.DIFFZ .AND. ZA+ZB+ZC+ZD.EQ.0) THEN 178 DO M = 1,NEXQ 179 BATCH (M,I,J,K,L) = TEMP2 (M) 180 END DO 181 ELSE 182 DO M = 1,NEXQ 183 BATCH (M,I,J,K,L) = 184 + TEMP2 (M) 185 + * INT2DZ (M,ZA,ZB,ZC,ZD) 186 END DO 187 END IF 188 189 116 CONTINUE 190 114 CONTINUE 191 112 CONTINUE 192 110 CONTINUE 193 XAP = I 194 106 CONTINUE 195 XBP = J 196 104 CONTINUE 197 XCP = K 198 102 CONTINUE 199 XDP = L 200 100 CONTINUE 201 202 RETURN 203C 204C 205C ******************** 206C * # of roots = 2 * 207C ******************** 208C 209C 210 2 XDP = 0 211 DO 200 XD = SHELLD,0,-1 212 YDMAX = SHELLD - XD 213 XCP = 0 214 DO 202 XC = SHELLC,0,-1 215 YCMAX = SHELLC - XC 216 XBP = 0 217 DO 204 XB = SHELLB,0,-1 218 YBMAX = SHELLB - XB 219 XAP = 0 220 DO 206 XA = SHELLA,0,-1 221 YAMAX = SHELLA - XA 222 223 DO M = 1,NGQEXQ 224 TEMP1 (M) = SCALE (M) * INT2DX (M,XA,XB,XC,XD) 225 END DO 226 227 L = XDP 228 DO 210 YD = YDMAX,0,-1 229 L = L + 1 230 ZD = YDMAX - YD 231 K = XCP 232 DO 212 YC = YCMAX,0,-1 233 K = K + 1 234 ZC = YCMAX - YC 235 J = XBP 236 DO 214 YB = YBMAX,0,-1 237 J = J + 1 238 ZB = YBMAX - YB 239 I = XAP 240 DO 216 YA = YAMAX,0,-1 241 I = I + 1 242 ZA = YAMAX - YA 243 244 IF (.NOT.DIFFY .AND. YA+YB+YC+YD.EQ.0) THEN 245 DO N = 1,NGQEXQ 246 TEMP2 (N) = TEMP1 (N) 247 END DO 248 ELSE 249 DO N = 1,NGQEXQ 250 TEMP2 (N) = TEMP1 (N) 251 + * INT2DY (N,YA,YB,YC,YD) 252 END DO 253 END IF 254 255 IF (.NOT.DIFFZ .AND. ZA+ZB+ZC+ZD.EQ.0) THEN 256 R = 1 257 DO M = 1,NEXQ 258 BATCH (M,I,J,K,L) = TEMP2 (R) 259 + + TEMP2 (R+1) 260 R = R + 2 261 END DO 262 ELSE 263 R = 1 264 DO M = 1,NEXQ 265 BATCH (M,I,J,K,L) = 266 + TEMP2 (R) 267 + * INT2DZ (R,ZA,ZB,ZC,ZD) 268 + + TEMP2 (R+1) 269 + * INT2DZ (R+1,ZA,ZB,ZC,ZD) 270 R = R + 2 271 END DO 272 END IF 273 274 216 CONTINUE 275 214 CONTINUE 276 212 CONTINUE 277 210 CONTINUE 278 XAP = I 279 206 CONTINUE 280 XBP = J 281 204 CONTINUE 282 XCP = K 283 202 CONTINUE 284 XDP = L 285 200 CONTINUE 286 287 RETURN 288C 289C 290C ******************** 291C * # of roots = 3 * 292C ******************** 293C 294C 295 3 XDP = 0 296 DO 300 XD = SHELLD,0,-1 297 YDMAX = SHELLD - XD 298 XCP = 0 299 DO 302 XC = SHELLC,0,-1 300 YCMAX = SHELLC - XC 301 XBP = 0 302 DO 304 XB = SHELLB,0,-1 303 YBMAX = SHELLB - XB 304 XAP = 0 305 DO 306 XA = SHELLA,0,-1 306 YAMAX = SHELLA - XA 307 308 DO M = 1,NGQEXQ 309 TEMP1 (M) = SCALE (M) * INT2DX (M,XA,XB,XC,XD) 310 END DO 311 312 L = XDP 313 DO 310 YD = YDMAX,0,-1 314 L = L + 1 315 ZD = YDMAX - YD 316 K = XCP 317 DO 312 YC = YCMAX,0,-1 318 K = K + 1 319 ZC = YCMAX - YC 320 J = XBP 321 DO 314 YB = YBMAX,0,-1 322 J = J + 1 323 ZB = YBMAX - YB 324 I = XAP 325 DO 316 YA = YAMAX,0,-1 326 I = I + 1 327 ZA = YAMAX - YA 328 329 IF (.NOT.DIFFY .AND. YA+YB+YC+YD.EQ.0) THEN 330 DO N = 1,NGQEXQ 331 TEMP2 (N) = TEMP1 (N) 332 END DO 333 ELSE 334 DO N = 1,NGQEXQ 335 TEMP2 (N) = TEMP1 (N) 336 + * INT2DY (N,YA,YB,YC,YD) 337 END DO 338 END IF 339 340 IF (.NOT.DIFFZ .AND. ZA+ZB+ZC+ZD.EQ.0) THEN 341 R = 1 342 DO M = 1,NEXQ 343 BATCH (M,I,J,K,L) = TEMP2 (R) 344 + + TEMP2 (R+1) 345 + + TEMP2 (R+2) 346 R = R + 3 347 END DO 348 ELSE 349 R = 1 350 DO M = 1,NEXQ 351 BATCH (M,I,J,K,L) = 352 + TEMP2 (R) 353 + * INT2DZ (R,ZA,ZB,ZC,ZD) 354 + + TEMP2 (R+1) 355 + * INT2DZ (R+1,ZA,ZB,ZC,ZD) 356 + + TEMP2 (R+2) 357 + * INT2DZ (R+2,ZA,ZB,ZC,ZD) 358 R = R + 3 359 END DO 360 END IF 361 362 316 CONTINUE 363 314 CONTINUE 364 312 CONTINUE 365 310 CONTINUE 366 XAP = I 367 306 CONTINUE 368 XBP = J 369 304 CONTINUE 370 XCP = K 371 302 CONTINUE 372 XDP = L 373 300 CONTINUE 374 375 RETURN 376C 377C 378C ******************** 379C * # of roots = 4 * 380C ******************** 381C 382C 383 4 XDP = 0 384 DO 400 XD = SHELLD,0,-1 385 YDMAX = SHELLD - XD 386 XCP = 0 387 DO 402 XC = SHELLC,0,-1 388 YCMAX = SHELLC - XC 389 XBP = 0 390 DO 404 XB = SHELLB,0,-1 391 YBMAX = SHELLB - XB 392 XAP = 0 393 DO 406 XA = SHELLA,0,-1 394 YAMAX = SHELLA - XA 395 396 DO M = 1,NGQEXQ 397 TEMP1 (M) = SCALE (M) * INT2DX (M,XA,XB,XC,XD) 398 END DO 399 400 L = XDP 401 DO 410 YD = YDMAX,0,-1 402 L = L + 1 403 ZD = YDMAX - YD 404 K = XCP 405 DO 412 YC = YCMAX,0,-1 406 K = K + 1 407 ZC = YCMAX - YC 408 J = XBP 409 DO 414 YB = YBMAX,0,-1 410 J = J + 1 411 ZB = YBMAX - YB 412 I = XAP 413 DO 416 YA = YAMAX,0,-1 414 I = I + 1 415 ZA = YAMAX - YA 416 417 IF (.NOT.DIFFY .AND. YA+YB+YC+YD.EQ.0) THEN 418 DO N = 1,NGQEXQ 419 TEMP2 (N) = TEMP1 (N) 420 END DO 421 ELSE 422 DO N = 1,NGQEXQ 423 TEMP2 (N) = TEMP1 (N) 424 + * INT2DY (N,YA,YB,YC,YD) 425 END DO 426 END IF 427 428 IF (.NOT.DIFFZ .AND. ZA+ZB+ZC+ZD.EQ.0) THEN 429 R = 1 430 DO M = 1,NEXQ 431 BATCH (M,I,J,K,L) = TEMP2 (R) 432 + + TEMP2 (R+1) 433 + + TEMP2 (R+2) 434 + + TEMP2 (R+3) 435 R = R + 4 436 END DO 437 ELSE 438 R = 1 439 DO M = 1,NEXQ 440 BATCH (M,I,J,K,L) = 441 + TEMP2 (R) 442 + * INT2DZ (R,ZA,ZB,ZC,ZD) 443 + + TEMP2 (R+1) 444 + * INT2DZ (R+1,ZA,ZB,ZC,ZD) 445 + + TEMP2 (R+2) 446 + * INT2DZ (R+2,ZA,ZB,ZC,ZD) 447 + + TEMP2 (R+3) 448 + * INT2DZ (R+3,ZA,ZB,ZC,ZD) 449 R = R + 4 450 END DO 451 END IF 452 453 416 CONTINUE 454 414 CONTINUE 455 412 CONTINUE 456 410 CONTINUE 457 XAP = I 458 406 CONTINUE 459 XBP = J 460 404 CONTINUE 461 XCP = K 462 402 CONTINUE 463 XDP = L 464 400 CONTINUE 465 466 RETURN 467C 468C 469C ******************** 470C * # of roots = 5 * 471C ******************** 472C 473C 474 5 XDP = 0 475 DO 500 XD = SHELLD,0,-1 476 YDMAX = SHELLD - XD 477 XCP = 0 478 DO 502 XC = SHELLC,0,-1 479 YCMAX = SHELLC - XC 480 XBP = 0 481 DO 504 XB = SHELLB,0,-1 482 YBMAX = SHELLB - XB 483 XAP = 0 484 DO 506 XA = SHELLA,0,-1 485 YAMAX = SHELLA - XA 486 487 DO M = 1,NGQEXQ 488 TEMP1 (M) = SCALE (M) * INT2DX (M,XA,XB,XC,XD) 489 END DO 490 491 L = XDP 492 DO 510 YD = YDMAX,0,-1 493 L = L + 1 494 ZD = YDMAX - YD 495 K = XCP 496 DO 512 YC = YCMAX,0,-1 497 K = K + 1 498 ZC = YCMAX - YC 499 J = XBP 500 DO 514 YB = YBMAX,0,-1 501 J = J + 1 502 ZB = YBMAX - YB 503 I = XAP 504 DO 516 YA = YAMAX,0,-1 505 I = I + 1 506 ZA = YAMAX - YA 507 508 IF (.NOT.DIFFY .AND. YA+YB+YC+YD.EQ.0) THEN 509 DO N = 1,NGQEXQ 510 TEMP2 (N) = TEMP1 (N) 511 END DO 512 ELSE 513 DO N = 1,NGQEXQ 514 TEMP2 (N) = TEMP1 (N) 515 + * INT2DY (N,YA,YB,YC,YD) 516 END DO 517 END IF 518 519 IF (.NOT.DIFFZ .AND. ZA+ZB+ZC+ZD.EQ.0) THEN 520 R = 1 521 DO M = 1,NEXQ 522 BATCH (M,I,J,K,L) = TEMP2 (R) 523 + + TEMP2 (R+1) 524 + + TEMP2 (R+2) 525 + + TEMP2 (R+3) 526 + + TEMP2 (R+4) 527 R = R + 5 528 END DO 529 ELSE 530 R = 1 531 DO M = 1,NEXQ 532 BATCH (M,I,J,K,L) = 533 + TEMP2 (R) 534 + * INT2DZ (R,ZA,ZB,ZC,ZD) 535 + + TEMP2 (R+1) 536 + * INT2DZ (R+1,ZA,ZB,ZC,ZD) 537 + + TEMP2 (R+2) 538 + * INT2DZ (R+2,ZA,ZB,ZC,ZD) 539 + + TEMP2 (R+3) 540 + * INT2DZ (R+3,ZA,ZB,ZC,ZD) 541 + + TEMP2 (R+4) 542 + * INT2DZ (R+4,ZA,ZB,ZC,ZD) 543 R = R + 5 544 END DO 545 END IF 546 547 516 CONTINUE 548 514 CONTINUE 549 512 CONTINUE 550 510 CONTINUE 551 XAP = I 552 506 CONTINUE 553 XBP = J 554 504 CONTINUE 555 XCP = K 556 502 CONTINUE 557 XDP = L 558 500 CONTINUE 559 560 RETURN 561C 562C 563C ******************** 564C * # of roots = 6 * 565C ******************** 566C 567C 568 6 XDP = 0 569 DO 600 XD = SHELLD,0,-1 570 YDMAX = SHELLD - XD 571 XCP = 0 572 DO 602 XC = SHELLC,0,-1 573 YCMAX = SHELLC - XC 574 XBP = 0 575 DO 604 XB = SHELLB,0,-1 576 YBMAX = SHELLB - XB 577 XAP = 0 578 DO 606 XA = SHELLA,0,-1 579 YAMAX = SHELLA - XA 580 581 DO M = 1,NGQEXQ 582 TEMP1 (M) = SCALE (M) * INT2DX (M,XA,XB,XC,XD) 583 END DO 584 585 L = XDP 586 DO 610 YD = YDMAX,0,-1 587 L = L + 1 588 ZD = YDMAX - YD 589 K = XCP 590 DO 612 YC = YCMAX,0,-1 591 K = K + 1 592 ZC = YCMAX - YC 593 J = XBP 594 DO 614 YB = YBMAX,0,-1 595 J = J + 1 596 ZB = YBMAX - YB 597 I = XAP 598 DO 616 YA = YAMAX,0,-1 599 I = I + 1 600 ZA = YAMAX - YA 601 602 IF (.NOT.DIFFY .AND. YA+YB+YC+YD.EQ.0) THEN 603 DO N = 1,NGQEXQ 604 TEMP2 (N) = TEMP1 (N) 605 END DO 606 ELSE 607 DO N = 1,NGQEXQ 608 TEMP2 (N) = TEMP1 (N) 609 + * INT2DY (N,YA,YB,YC,YD) 610 END DO 611 END IF 612 613 IF (.NOT.DIFFZ .AND. ZA+ZB+ZC+ZD.EQ.0) THEN 614 R = 1 615 DO M = 1,NEXQ 616 BATCH (M,I,J,K,L) = TEMP2 (R) 617 + + TEMP2 (R+1) 618 + + TEMP2 (R+2) 619 + + TEMP2 (R+3) 620 + + TEMP2 (R+4) 621 + + TEMP2 (R+5) 622 R = R + 6 623 END DO 624 ELSE 625 R = 1 626 DO M = 1,NEXQ 627 BATCH (M,I,J,K,L) = 628 + TEMP2 (R) 629 + * INT2DZ (R,ZA,ZB,ZC,ZD) 630 + + TEMP2 (R+1) 631 + * INT2DZ (R+1,ZA,ZB,ZC,ZD) 632 + + TEMP2 (R+2) 633 + * INT2DZ (R+2,ZA,ZB,ZC,ZD) 634 + + TEMP2 (R+3) 635 + * INT2DZ (R+3,ZA,ZB,ZC,ZD) 636 + + TEMP2 (R+4) 637 + * INT2DZ (R+4,ZA,ZB,ZC,ZD) 638 + + TEMP2 (R+5) 639 + * INT2DZ (R+5,ZA,ZB,ZC,ZD) 640 R = R + 6 641 END DO 642 END IF 643 644 616 CONTINUE 645 614 CONTINUE 646 612 CONTINUE 647 610 CONTINUE 648 XAP = I 649 606 CONTINUE 650 XBP = J 651 604 CONTINUE 652 XCP = K 653 602 CONTINUE 654 XDP = L 655 600 CONTINUE 656 657 RETURN 658C 659C 660C ******************** 661C * # of roots = 7 * 662C ******************** 663C 664C 665 7 XDP = 0 666 DO 700 XD = SHELLD,0,-1 667 YDMAX = SHELLD - XD 668 XCP = 0 669 DO 702 XC = SHELLC,0,-1 670 YCMAX = SHELLC - XC 671 XBP = 0 672 DO 704 XB = SHELLB,0,-1 673 YBMAX = SHELLB - XB 674 XAP = 0 675 DO 706 XA = SHELLA,0,-1 676 YAMAX = SHELLA - XA 677 678 DO M = 1,NGQEXQ 679 TEMP1 (M) = SCALE (M) * INT2DX (M,XA,XB,XC,XD) 680 END DO 681 682 L = XDP 683 DO 710 YD = YDMAX,0,-1 684 L = L + 1 685 ZD = YDMAX - YD 686 K = XCP 687 DO 712 YC = YCMAX,0,-1 688 K = K + 1 689 ZC = YCMAX - YC 690 J = XBP 691 DO 714 YB = YBMAX,0,-1 692 J = J + 1 693 ZB = YBMAX - YB 694 I = XAP 695 DO 716 YA = YAMAX,0,-1 696 I = I + 1 697 ZA = YAMAX - YA 698 699 IF (.NOT.DIFFY .AND. YA+YB+YC+YD.EQ.0) THEN 700 DO N = 1,NGQEXQ 701 TEMP2 (N) = TEMP1 (N) 702 END DO 703 ELSE 704 DO N = 1,NGQEXQ 705 TEMP2 (N) = TEMP1 (N) 706 + * INT2DY (N,YA,YB,YC,YD) 707 END DO 708 END IF 709 710 IF (.NOT.DIFFZ .AND. ZA+ZB+ZC+ZD.EQ.0) THEN 711 R = 1 712 DO M = 1,NEXQ 713 BATCH (M,I,J,K,L) = TEMP2 (R) 714 + + TEMP2 (R+1) 715 + + TEMP2 (R+2) 716 + + TEMP2 (R+3) 717 + + TEMP2 (R+4) 718 + + TEMP2 (R+5) 719 + + TEMP2 (R+6) 720 R = R + 7 721 END DO 722 ELSE 723 R = 1 724 DO M = 1,NEXQ 725 BATCH (M,I,J,K,L) = 726 + TEMP2 (R) 727 + * INT2DZ (R,ZA,ZB,ZC,ZD) 728 + + TEMP2 (R+1) 729 + * INT2DZ (R+1,ZA,ZB,ZC,ZD) 730 + + TEMP2 (R+2) 731 + * INT2DZ (R+2,ZA,ZB,ZC,ZD) 732 + + TEMP2 (R+3) 733 + * INT2DZ (R+3,ZA,ZB,ZC,ZD) 734 + + TEMP2 (R+4) 735 + * INT2DZ (R+4,ZA,ZB,ZC,ZD) 736 + + TEMP2 (R+5) 737 + * INT2DZ (R+5,ZA,ZB,ZC,ZD) 738 + + TEMP2 (R+6) 739 + * INT2DZ (R+6,ZA,ZB,ZC,ZD) 740 R = R + 7 741 END DO 742 END IF 743 744 716 CONTINUE 745 714 CONTINUE 746 712 CONTINUE 747 710 CONTINUE 748 XAP = I 749 706 CONTINUE 750 XBP = J 751 704 CONTINUE 752 XCP = K 753 702 CONTINUE 754 XDP = L 755 700 CONTINUE 756 757 RETURN 758C 759C 760C ******************** 761C * # of roots = 8 * 762C ******************** 763C 764C 765 8 XDP = 0 766 DO 800 XD = SHELLD,0,-1 767 YDMAX = SHELLD - XD 768 XCP = 0 769 DO 802 XC = SHELLC,0,-1 770 YCMAX = SHELLC - XC 771 XBP = 0 772 DO 804 XB = SHELLB,0,-1 773 YBMAX = SHELLB - XB 774 XAP = 0 775 DO 806 XA = SHELLA,0,-1 776 YAMAX = SHELLA - XA 777 778 DO M = 1,NGQEXQ 779 TEMP1 (M) = SCALE (M) * INT2DX (M,XA,XB,XC,XD) 780 END DO 781 782 L = XDP 783 DO 810 YD = YDMAX,0,-1 784 L = L + 1 785 ZD = YDMAX - YD 786 K = XCP 787 DO 812 YC = YCMAX,0,-1 788 K = K + 1 789 ZC = YCMAX - YC 790 J = XBP 791 DO 814 YB = YBMAX,0,-1 792 J = J + 1 793 ZB = YBMAX - YB 794 I = XAP 795 DO 816 YA = YAMAX,0,-1 796 I = I + 1 797 ZA = YAMAX - YA 798 799 IF (.NOT.DIFFY .AND. YA+YB+YC+YD.EQ.0) THEN 800 DO N = 1,NGQEXQ 801 TEMP2 (N) = TEMP1 (N) 802 END DO 803 ELSE 804 DO N = 1,NGQEXQ 805 TEMP2 (N) = TEMP1 (N) 806 + * INT2DY (N,YA,YB,YC,YD) 807 END DO 808 END IF 809 810 IF (.NOT.DIFFZ .AND. ZA+ZB+ZC+ZD.EQ.0) THEN 811 R = 1 812 DO M = 1,NEXQ 813 BATCH (M,I,J,K,L) = TEMP2 (R) 814 + + TEMP2 (R+1) 815 + + TEMP2 (R+2) 816 + + TEMP2 (R+3) 817 + + TEMP2 (R+4) 818 + + TEMP2 (R+5) 819 + + TEMP2 (R+6) 820 + + TEMP2 (R+7) 821 R = R + 8 822 END DO 823 ELSE 824 R = 1 825 DO M = 1,NEXQ 826 BATCH (M,I,J,K,L) = 827 + TEMP2 (R) 828 + * INT2DZ (R,ZA,ZB,ZC,ZD) 829 + + TEMP2 (R+1) 830 + * INT2DZ (R+1,ZA,ZB,ZC,ZD) 831 + + TEMP2 (R+2) 832 + * INT2DZ (R+2,ZA,ZB,ZC,ZD) 833 + + TEMP2 (R+3) 834 + * INT2DZ (R+3,ZA,ZB,ZC,ZD) 835 + + TEMP2 (R+4) 836 + * INT2DZ (R+4,ZA,ZB,ZC,ZD) 837 + + TEMP2 (R+5) 838 + * INT2DZ (R+5,ZA,ZB,ZC,ZD) 839 + + TEMP2 (R+6) 840 + * INT2DZ (R+6,ZA,ZB,ZC,ZD) 841 + + TEMP2 (R+7) 842 + * INT2DZ (R+7,ZA,ZB,ZC,ZD) 843 R = R + 8 844 END DO 845 END IF 846 847 816 CONTINUE 848 814 CONTINUE 849 812 CONTINUE 850 810 CONTINUE 851 XAP = I 852 806 CONTINUE 853 XBP = J 854 804 CONTINUE 855 XCP = K 856 802 CONTINUE 857 XDP = L 858 800 CONTINUE 859 860 RETURN 861C 862C 863C ******************** 864C * # of roots = 9 * 865C ******************** 866C 867C 868 9 XDP = 0 869 DO 900 XD = SHELLD,0,-1 870 YDMAX = SHELLD - XD 871 XCP = 0 872 DO 902 XC = SHELLC,0,-1 873 YCMAX = SHELLC - XC 874 XBP = 0 875 DO 904 XB = SHELLB,0,-1 876 YBMAX = SHELLB - XB 877 XAP = 0 878 DO 906 XA = SHELLA,0,-1 879 YAMAX = SHELLA - XA 880 881 DO M = 1,NGQEXQ 882 TEMP1 (M) = SCALE (M) * INT2DX (M,XA,XB,XC,XD) 883 END DO 884 885 L = XDP 886 DO 910 YD = YDMAX,0,-1 887 L = L + 1 888 ZD = YDMAX - YD 889 K = XCP 890 DO 912 YC = YCMAX,0,-1 891 K = K + 1 892 ZC = YCMAX - YC 893 J = XBP 894 DO 914 YB = YBMAX,0,-1 895 J = J + 1 896 ZB = YBMAX - YB 897 I = XAP 898 DO 916 YA = YAMAX,0,-1 899 I = I + 1 900 ZA = YAMAX - YA 901 902 IF (.NOT.DIFFY .AND. YA+YB+YC+YD.EQ.0) THEN 903 DO N = 1,NGQEXQ 904 TEMP2 (N) = TEMP1 (N) 905 END DO 906 ELSE 907 DO N = 1,NGQEXQ 908 TEMP2 (N) = TEMP1 (N) 909 + * INT2DY (N,YA,YB,YC,YD) 910 END DO 911 END IF 912 913 IF (.NOT.DIFFZ .AND. ZA+ZB+ZC+ZD.EQ.0) THEN 914 R = 1 915 DO M = 1,NEXQ 916 BATCH (M,I,J,K,L) = TEMP2 (R) 917 + + TEMP2 (R+1) 918 + + TEMP2 (R+2) 919 + + TEMP2 (R+3) 920 + + TEMP2 (R+4) 921 + + TEMP2 (R+5) 922 + + TEMP2 (R+6) 923 + + TEMP2 (R+7) 924 + + TEMP2 (R+8) 925 R = R + 9 926 END DO 927 ELSE 928 R = 1 929 DO M = 1,NEXQ 930 BATCH (M,I,J,K,L) = 931 + TEMP2 (R) 932 + * INT2DZ (R,ZA,ZB,ZC,ZD) 933 + + TEMP2 (R+1) 934 + * INT2DZ (R+1,ZA,ZB,ZC,ZD) 935 + + TEMP2 (R+2) 936 + * INT2DZ (R+2,ZA,ZB,ZC,ZD) 937 + + TEMP2 (R+3) 938 + * INT2DZ (R+3,ZA,ZB,ZC,ZD) 939 + + TEMP2 (R+4) 940 + * INT2DZ (R+4,ZA,ZB,ZC,ZD) 941 + + TEMP2 (R+5) 942 + * INT2DZ (R+5,ZA,ZB,ZC,ZD) 943 + + TEMP2 (R+6) 944 + * INT2DZ (R+6,ZA,ZB,ZC,ZD) 945 + + TEMP2 (R+7) 946 + * INT2DZ (R+7,ZA,ZB,ZC,ZD) 947 + + TEMP2 (R+8) 948 + * INT2DZ (R+8,ZA,ZB,ZC,ZD) 949 R = R + 9 950 END DO 951 END IF 952 953 916 CONTINUE 954 914 CONTINUE 955 912 CONTINUE 956 910 CONTINUE 957 XAP = I 958 906 CONTINUE 959 XBP = J 960 904 CONTINUE 961 XCP = K 962 902 CONTINUE 963 XDP = L 964 900 CONTINUE 965 966 RETURN 967C 968C 969C ******************** 970C * # of roots > 9 * 971C ******************** 972C 973C ...outer loops over x,x,x,x-quadruples. No skipping of 974C x,x,x,x-contribution of 0,0,0,0-type can be done here, 975C since the 2DX integrals carry the Rys weight! 976C 977C 978 10 XDP = 0 979 DO 1000 XD = SHELLD,0,-1 980 YDMAX = SHELLD - XD 981 XCP = 0 982 DO 1002 XC = SHELLC,0,-1 983 YCMAX = SHELLC - XC 984 XBP = 0 985 DO 1004 XB = SHELLB,0,-1 986 YBMAX = SHELLB - XB 987 XAP = 0 988 DO 1006 XA = SHELLA,0,-1 989 YAMAX = SHELLA - XA 990 991 DO M = 1,NGQEXQ 992 TEMP1 (M) = SCALE (M) * INT2DX (M,XA,XB,XC,XD) 993 END DO 994C 995C 996C ...inner loops over y,y,y,y-quadruples. Skip the 997C multiplication of y,y,y,y-contributions, if no 998C y-coordinate derivative was formed and we have a 999C 0,0,0,0-quadruple, as then the 2DY ABCD integrals 1000C are equal to 1. 1001C 1002C 1003 L = XDP 1004 DO 1010 YD = YDMAX,0,-1 1005 L = L + 1 1006 ZD = YDMAX - YD 1007 K = XCP 1008 DO 1012 YC = YCMAX,0,-1 1009 K = K + 1 1010 ZC = YCMAX - YC 1011 J = XBP 1012 DO 1014 YB = YBMAX,0,-1 1013 J = J + 1 1014 ZB = YBMAX - YB 1015 I = XAP 1016 DO 1016 YA = YAMAX,0,-1 1017 I = I + 1 1018 ZA = YAMAX - YA 1019 1020 IF (.NOT.DIFFY .AND. YA+YB+YC+YD.EQ.0) THEN 1021 DO N = 1,NGQEXQ 1022 TEMP2 (N) = TEMP1 (N) 1023 END DO 1024 ELSE 1025 DO N = 1,NGQEXQ 1026 TEMP2 (N) = TEMP1 (N) 1027 + * INT2DY (N,YA,YB,YC,YD) 1028 END DO 1029 END IF 1030C 1031C 1032C ...skip multiplication of z,z,z,z-contributions, if we 1033C have a 0,0,0,0-quadruple and no derivations were 1034C performed on the z-coordinate, as then the 2DZ ABCD 1035C integrals are equal to 1. All info concerning all 1036C three x,x,x,x-, y,y,y,y- and z,z,z,z-quadruples have 1037C been collected for all exponent quadruplets at once. 1038C Sum up the 2D X,Y,Z integral products to the 1039C appropriate places of the [AB|CD] batch. 1040C 1041C 1042 IF (.NOT.DIFFZ .AND. ZA+ZB+ZC+ZD.EQ.0) THEN 1043 R = 0 1044 DO M = 1,NEXQ 1045 SUM = ZERO 1046 DO N = 1,NGQP 1047 SUM = SUM + TEMP2 (R+N) 1048 END DO 1049 R = R + NGQP 1050 BATCH (M,I,J,K,L) = SUM 1051 END DO 1052 ELSE 1053 R = 0 1054 DO M = 1,NEXQ 1055 SUM = ZERO 1056 DO N = 1,NGQP 1057 SUM = SUM + TEMP2 (R+N) 1058 + * INT2DZ (R+N,ZA,ZB,ZC,ZD) 1059 END DO 1060 R = R + NGQP 1061 BATCH (M,I,J,K,L) = SUM 1062 END DO 1063 END IF 1064C 1065C 1066C ...next z,z,z,z-quadruple and y,y,y,y-quadruple. 1067C 1068C 1069 1016 CONTINUE 1070 1014 CONTINUE 1071 1012 CONTINUE 1072 1010 CONTINUE 1073C 1074C 1075C ...next x,x,x,x-quadruple. 1076C 1077C 1078 XAP = I 1079 1006 CONTINUE 1080 XBP = J 1081 1004 CONTINUE 1082 XCP = K 1083 1002 CONTINUE 1084 XDP = L 1085 1000 CONTINUE 1086C 1087C 1088C ...ready! 1089C 1090C 1091 RETURN 1092 END 1093