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_E0CD 16 + 17 + ( SHELLA,SHELLP,SHELLC,SHELLD, 18 + NGQP,NEXQ,NGQEXQ, 19 + NXYZET,NXYZP,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_E0CD 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 [E0|CD] , E = A to P, adding up the 34C contributions from all the 2D PCD derivative 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, which can be achieved by having 40C the outer loops run over all possible x and y monomial 41C parts and the inner loops over all allowed E shell 42C combinations. The price to pay for such loop ordering 43C is the scattered addressing of locations within the 44C batch array, which has its rows and columns ordered 45C such that the E shells are increasing and within each 46C E shell the monomials are ordered such that x>y>z in 47C exponents. 48C 49C An example follows: 50C ------------------- 51C 52C Let E = 0,2 and C = 0 and D = 1. Then we have the left 53C and right hand of the batch array ordered as follows: 54C 55C left xyz right xyz 56C 57C 000 000 100 58C --- 010 59C 100 001 60C 010 61C 001 62C --- 63C 200 -> -5 64C 110 65C 101 -> -3 66C 020 67C 011 68C 002 -> 0 69C 70C The batch would thus have dimensions 10 x 1 x 3. 71C For the left side the reduced multiplication scheme 72C would have its outermost loop run over x=0,2, followed 73C by the next loop y=0,2-x. The innermost loop would then 74C run over the allowed shells E=E(max),max(E(min),x+y). 75C In this case all x,y-triples can be reused for all 76C appropriate shell combinations. 77C 78C To find the address of a specific x,y,z,E combination 79C inside the batch array, we first note that the z-part 80C is dependent on the x,y-parts and is hence not needed. 81C Lets look at the E-part first. The E-part is evaluated 82C from its dimension formula (E+1)*(E+2)/2. Organizing 83C the inner E-loop to run from E(max) always, the 84C dimension for E(max) is passed as the argument NXYZP 85C and all lower E dimensions are calculated by the 86C formula relating dimensions between E and E+1: 87C 88C dim(E+1) = dim(E) + E + 2 89C 90C In this way multiplications in the E-part can be 91C entirely avoided. The x,y-part is defined as the 92C part which has to be subtracted from dim(E) to 93C reach the xyz monomial position inside the E shell. 94C It can be divided into an x-part and a y-part. The 95C x-part is given by the fomula: 96C 97C x-part = - x*E + x(x-3)/2 98C 99C and for the example above has been given for E=2 100C and marked with arrows ->. The last term of the x-part 101C involves 1 multiplication and division, however it 102C can be changed to: 103C 104C x-1 105C x-part = - x*E - x + sum i 106C i=0 107C 108C and clever additions inside the x-loop avoid the use 109C of multiplications and divisions. The y-part is trivial 110C and is simply equal to -y. The overall conclusion is 111C thus that the location of a specific x,y,z,E quadruple 112C inside the batch comes at the cost of one x*E(max) 113C multiplication in the outermost x-loops, since the 114C other x*E ones can again be reached via stepwise 115C subtraction of x from x*E(max). 116C 117C Due to the very computational intensive steps inside 118C the x,y,z,E loops, special sections of identical 119C x,y,z,E loop structures have been given for each 120C # of roots =< 9, thus saving considerable computing 121C time over the general case. 122C 123C For comments on how the x,y,z,E loop structures are 124C coded please refer to the general root case. 125C 126C 127C Input: 128C 129C SHELLx = shell types for individual csh 130C x=A,C,D and csh sum P=A+B 131C NGQP = # of gaussian quadrature points 132C (roots) 133C NEXQ = current # of exponent quadruplets 134C NGQEXQ = product of # of gaussian quadrature 135C points times exponent quadruplets 136C NXYZET = sum of # of cartesian monomials 137C for all shells in the range 138C E = A,...,P=A+B 139C NXYZy = # of cartesian monomials for 140C y = P,C,D shells 141C INT2Dx = all current 2D PCD derivative 142C integrals for each cartesian 143C component (x = X,Y,Z) 144C DIFFx = is true, if differentiation was 145C performed along the x=Y,Z direction 146C TEMP1(2) = scratch arrays holding intermediate 147C 2D PCD derivative integral products 148C SCALE = the NGQEXQ scaling factors 149C 150C 151C Output: 152C 153C BATCH = batch of primitive cartesian 154C [E0|CD] derivative integrals 155C corresponding to all current 156C exponent quadruplets 157C 158C 159C AUTHOR : Norbert Flocke 160C------------------------------------------------------------------------ 161C 162C 163C ...include files and declare variables. 164C 165C 166 IMPLICIT NONE 167 168 LOGICAL DIFFY,DIFFZ 169 170 INTEGER I,J,K,M,N,R 171 INTEGER NGQP,NEXQ,NGQEXQ 172 INTEGER NXYZE,NXYZET,NXYZP,NXYZC,NXYZD 173 INTEGER SE 174 INTEGER SEEND 175 INTEGER SHELLA,SHELLP,SHELLC,SHELLD 176 INTEGER XC,YC,ZC,XD,YD,ZD 177 INTEGER XCP,XDP 178 INTEGER XE,YE,ZE 179 INTEGER XEMAX 180 INTEGER XEP,XYEP 181 INTEGER XYE 182 INTEGER YCMAX,YDMAX 183 INTEGER YEEND 184 185 DOUBLE PRECISION SUM 186 DOUBLE PRECISION ZERO 187 188 DOUBLE PRECISION SCALE (1:NGQEXQ) 189 DOUBLE PRECISION TEMP1 (1:NGQEXQ) 190 DOUBLE PRECISION TEMP2 (1:NGQEXQ) 191 192 DOUBLE PRECISION BATCH (1:NEXQ,1:NXYZET,1:NXYZC,1:NXYZD) 193 194 DOUBLE PRECISION INT2DX (1:NGQEXQ,0:SHELLP,0:SHELLC,0:SHELLD) 195 DOUBLE PRECISION INT2DY (1:NGQEXQ,0:SHELLP,0:SHELLC,0:SHELLD) 196 DOUBLE PRECISION INT2DZ (1:NGQEXQ,0:SHELLP,0:SHELLC,0:SHELLD) 197 198 PARAMETER (ZERO = 0.D0) 199C 200C 201C------------------------------------------------------------------------ 202C 203C 204C ...jump according to number of roots. 205C 206C 207 GOTO (1,2,3,4,5,6,7,8,9,10) MIN (NGQP,10) 208C 209C 210C ******************** 211C * # of roots = 1 * 212C ******************** 213C 214C 215 1 XDP = 0 216 DO 100 XD = SHELLD,0,-1 217 YDMAX = SHELLD - XD 218 XCP = 0 219 DO 102 XC = SHELLC,0,-1 220 YCMAX = SHELLC - XC 221 222 XEP = NXYZET + 3 223 DO 104 XE = 0,SHELLP 224 XEP = XEP + XE - 2 225 XEMAX = XE * SHELLP 226 YEEND = SHELLP - XE 227 228 DO M = 1,NEXQ 229 TEMP1 (M) = SCALE (M) * INT2DX (M,XE,XC,XD) 230 END DO 231 232 K = XDP 233 DO 110 YD = YDMAX,0,-1 234 K = K + 1 235 ZD = YDMAX - YD 236 J = XCP 237 DO 112 YC = YCMAX,0,-1 238 J = J + 1 239 ZC = YCMAX - YC 240 241 XYEP = XEP - XEMAX 242 DO 114 YE = 0,YEEND 243 XYE = XE + YE 244 XYEP = XYEP - 1 245 SEEND = MAX0 (SHELLA,XYE) 246 247 IF (.NOT.DIFFY .AND. YE+YC+YD.EQ.0) THEN 248 DO M = 1,NEXQ 249 TEMP2 (M) = TEMP1 (M) 250 END DO 251 ELSE 252 DO M = 1,NEXQ 253 TEMP2 (M) = TEMP1 (M) 254 + * INT2DY (M,YE,YC,YD) 255 END DO 256 END IF 257 258 I = XYEP 259 NXYZE = NXYZP 260 DO 120 SE = SHELLP,SEEND,-1 261 ZE = SE - XYE 262 263 IF (.NOT.DIFFZ .AND. ZE+ZC+ZD.EQ.0) THEN 264 DO M = 1,NEXQ 265 BATCH (M,I,J,K) = TEMP2 (M) 266 END DO 267 ELSE 268 DO M = 1,NEXQ 269 BATCH (M,I,J,K) = 270 + TEMP2 (M) 271 + * INT2DZ (M,ZE,ZC,ZD) 272 END DO 273 END IF 274 275 I = I - NXYZE + XE 276 NXYZE = NXYZE - SE - 1 277 120 CONTINUE 278 114 CONTINUE 279 112 CONTINUE 280 110 CONTINUE 281 104 CONTINUE 282 XCP = J 283 102 CONTINUE 284 XDP = K 285 100 CONTINUE 286 287 RETURN 288C 289C 290C ******************** 291C * # of roots = 2 * 292C ******************** 293C 294C 295 2 XDP = 0 296 DO 200 XD = SHELLD,0,-1 297 YDMAX = SHELLD - XD 298 XCP = 0 299 DO 202 XC = SHELLC,0,-1 300 YCMAX = SHELLC - XC 301 302 XEP = NXYZET + 3 303 DO 204 XE = 0,SHELLP 304 XEP = XEP + XE - 2 305 XEMAX = XE * SHELLP 306 YEEND = SHELLP - XE 307 308 DO M = 1,NGQEXQ 309 TEMP1 (M) = SCALE (M) * INT2DX (M,XE,XC,XD) 310 END DO 311 312 K = XDP 313 DO 210 YD = YDMAX,0,-1 314 K = K + 1 315 ZD = YDMAX - YD 316 J = XCP 317 DO 212 YC = YCMAX,0,-1 318 J = J + 1 319 ZC = YCMAX - YC 320 321 XYEP = XEP - XEMAX 322 DO 214 YE = 0,YEEND 323 XYE = XE + YE 324 XYEP = XYEP - 1 325 SEEND = MAX0 (SHELLA,XYE) 326 327 IF (.NOT.DIFFY .AND. YE+YC+YD.EQ.0) THEN 328 DO N = 1,NGQEXQ 329 TEMP2 (N) = TEMP1 (N) 330 END DO 331 ELSE 332 DO N = 1,NGQEXQ 333 TEMP2 (N) = TEMP1 (N) 334 + * INT2DY (N,YE,YC,YD) 335 END DO 336 END IF 337 338 I = XYEP 339 NXYZE = NXYZP 340 DO 220 SE = SHELLP,SEEND,-1 341 ZE = SE - XYE 342 343 IF (.NOT.DIFFZ .AND. ZE+ZC+ZD.EQ.0) THEN 344 R = 1 345 DO M = 1,NEXQ 346 BATCH (M,I,J,K) = TEMP2 (R) 347 + + TEMP2 (R+1) 348 R = R + 2 349 END DO 350 ELSE 351 R = 1 352 DO M = 1,NEXQ 353 BATCH (M,I,J,K) = 354 + TEMP2 (R) 355 + * INT2DZ (R,ZE,ZC,ZD) 356 + + TEMP2 (R+1) 357 + * INT2DZ (R+1,ZE,ZC,ZD) 358 R = R + 2 359 END DO 360 END IF 361 362 I = I - NXYZE + XE 363 NXYZE = NXYZE - SE - 1 364 220 CONTINUE 365 214 CONTINUE 366 212 CONTINUE 367 210 CONTINUE 368 204 CONTINUE 369 XCP = J 370 202 CONTINUE 371 XDP = K 372 200 CONTINUE 373 374 RETURN 375C 376C 377C ******************** 378C * # of roots = 3 * 379C ******************** 380C 381C 382 3 XDP = 0 383 DO 300 XD = SHELLD,0,-1 384 YDMAX = SHELLD - XD 385 XCP = 0 386 DO 302 XC = SHELLC,0,-1 387 YCMAX = SHELLC - XC 388 389 XEP = NXYZET + 3 390 DO 304 XE = 0,SHELLP 391 XEP = XEP + XE - 2 392 XEMAX = XE * SHELLP 393 YEEND = SHELLP - XE 394 395 DO M = 1,NGQEXQ 396 TEMP1 (M) = SCALE (M) * INT2DX (M,XE,XC,XD) 397 END DO 398 399 K = XDP 400 DO 310 YD = YDMAX,0,-1 401 K = K + 1 402 ZD = YDMAX - YD 403 J = XCP 404 DO 312 YC = YCMAX,0,-1 405 J = J + 1 406 ZC = YCMAX - YC 407 408 XYEP = XEP - XEMAX 409 DO 314 YE = 0,YEEND 410 XYE = XE + YE 411 XYEP = XYEP - 1 412 SEEND = MAX0 (SHELLA,XYE) 413 414 IF (.NOT.DIFFY .AND. YE+YC+YD.EQ.0) THEN 415 DO N = 1,NGQEXQ 416 TEMP2 (N) = TEMP1 (N) 417 END DO 418 ELSE 419 DO N = 1,NGQEXQ 420 TEMP2 (N) = TEMP1 (N) 421 + * INT2DY (N,YE,YC,YD) 422 END DO 423 END IF 424 425 I = XYEP 426 NXYZE = NXYZP 427 DO 320 SE = SHELLP,SEEND,-1 428 ZE = SE - XYE 429 430 IF (.NOT.DIFFZ .AND. ZE+ZC+ZD.EQ.0) THEN 431 R = 1 432 DO M = 1,NEXQ 433 BATCH (M,I,J,K) = TEMP2 (R) 434 + + TEMP2 (R+1) 435 + + TEMP2 (R+2) 436 R = R + 3 437 END DO 438 ELSE 439 R = 1 440 DO M = 1,NEXQ 441 BATCH (M,I,J,K) = 442 + TEMP2 (R) 443 + * INT2DZ (R,ZE,ZC,ZD) 444 + + TEMP2 (R+1) 445 + * INT2DZ (R+1,ZE,ZC,ZD) 446 + + TEMP2 (R+2) 447 + * INT2DZ (R+2,ZE,ZC,ZD) 448 R = R + 3 449 END DO 450 END IF 451 452 I = I - NXYZE + XE 453 NXYZE = NXYZE - SE - 1 454 320 CONTINUE 455 314 CONTINUE 456 312 CONTINUE 457 310 CONTINUE 458 304 CONTINUE 459 XCP = J 460 302 CONTINUE 461 XDP = K 462 300 CONTINUE 463 464 RETURN 465C 466C 467C ******************** 468C * # of roots = 4 * 469C ******************** 470C 471C 472 4 XDP = 0 473 DO 400 XD = SHELLD,0,-1 474 YDMAX = SHELLD - XD 475 XCP = 0 476 DO 402 XC = SHELLC,0,-1 477 YCMAX = SHELLC - XC 478 479 XEP = NXYZET + 3 480 DO 404 XE = 0,SHELLP 481 XEP = XEP + XE - 2 482 XEMAX = XE * SHELLP 483 YEEND = SHELLP - XE 484 485 DO M = 1,NGQEXQ 486 TEMP1 (M) = SCALE (M) * INT2DX (M,XE,XC,XD) 487 END DO 488 489 K = XDP 490 DO 410 YD = YDMAX,0,-1 491 K = K + 1 492 ZD = YDMAX - YD 493 J = XCP 494 DO 412 YC = YCMAX,0,-1 495 J = J + 1 496 ZC = YCMAX - YC 497 498 XYEP = XEP - XEMAX 499 DO 414 YE = 0,YEEND 500 XYE = XE + YE 501 XYEP = XYEP - 1 502 SEEND = MAX0 (SHELLA,XYE) 503 504 IF (.NOT.DIFFY .AND. YE+YC+YD.EQ.0) THEN 505 DO N = 1,NGQEXQ 506 TEMP2 (N) = TEMP1 (N) 507 END DO 508 ELSE 509 DO N = 1,NGQEXQ 510 TEMP2 (N) = TEMP1 (N) 511 + * INT2DY (N,YE,YC,YD) 512 END DO 513 END IF 514 515 I = XYEP 516 NXYZE = NXYZP 517 DO 420 SE = SHELLP,SEEND,-1 518 ZE = SE - XYE 519 520 IF (.NOT.DIFFZ .AND. ZE+ZC+ZD.EQ.0) THEN 521 R = 1 522 DO M = 1,NEXQ 523 BATCH (M,I,J,K) = TEMP2 (R) 524 + + TEMP2 (R+1) 525 + + TEMP2 (R+2) 526 + + TEMP2 (R+3) 527 R = R + 4 528 END DO 529 ELSE 530 R = 1 531 DO M = 1,NEXQ 532 BATCH (M,I,J,K) = 533 + TEMP2 (R) 534 + * INT2DZ (R,ZE,ZC,ZD) 535 + + TEMP2 (R+1) 536 + * INT2DZ (R+1,ZE,ZC,ZD) 537 + + TEMP2 (R+2) 538 + * INT2DZ (R+2,ZE,ZC,ZD) 539 + + TEMP2 (R+3) 540 + * INT2DZ (R+3,ZE,ZC,ZD) 541 R = R + 4 542 END DO 543 END IF 544 545 I = I - NXYZE + XE 546 NXYZE = NXYZE - SE - 1 547 420 CONTINUE 548 414 CONTINUE 549 412 CONTINUE 550 410 CONTINUE 551 404 CONTINUE 552 XCP = J 553 402 CONTINUE 554 XDP = K 555 400 CONTINUE 556 557 RETURN 558C 559C 560C ******************** 561C * # of roots = 5 * 562C ******************** 563C 564C 565 5 XDP = 0 566 DO 500 XD = SHELLD,0,-1 567 YDMAX = SHELLD - XD 568 XCP = 0 569 DO 502 XC = SHELLC,0,-1 570 YCMAX = SHELLC - XC 571 572 XEP = NXYZET + 3 573 DO 504 XE = 0,SHELLP 574 XEP = XEP + XE - 2 575 XEMAX = XE * SHELLP 576 YEEND = SHELLP - XE 577 578 DO M = 1,NGQEXQ 579 TEMP1 (M) = SCALE (M) * INT2DX (M,XE,XC,XD) 580 END DO 581 582 K = XDP 583 DO 510 YD = YDMAX,0,-1 584 K = K + 1 585 ZD = YDMAX - YD 586 J = XCP 587 DO 512 YC = YCMAX,0,-1 588 J = J + 1 589 ZC = YCMAX - YC 590 591 XYEP = XEP - XEMAX 592 DO 514 YE = 0,YEEND 593 XYE = XE + YE 594 XYEP = XYEP - 1 595 SEEND = MAX0 (SHELLA,XYE) 596 597 IF (.NOT.DIFFY .AND. YE+YC+YD.EQ.0) THEN 598 DO N = 1,NGQEXQ 599 TEMP2 (N) = TEMP1 (N) 600 END DO 601 ELSE 602 DO N = 1,NGQEXQ 603 TEMP2 (N) = TEMP1 (N) 604 + * INT2DY (N,YE,YC,YD) 605 END DO 606 END IF 607 608 I = XYEP 609 NXYZE = NXYZP 610 DO 520 SE = SHELLP,SEEND,-1 611 ZE = SE - XYE 612 613 IF (.NOT.DIFFZ .AND. ZE+ZC+ZD.EQ.0) THEN 614 R = 1 615 DO M = 1,NEXQ 616 BATCH (M,I,J,K) = TEMP2 (R) 617 + + TEMP2 (R+1) 618 + + TEMP2 (R+2) 619 + + TEMP2 (R+3) 620 + + TEMP2 (R+4) 621 R = R + 5 622 END DO 623 ELSE 624 R = 1 625 DO M = 1,NEXQ 626 BATCH (M,I,J,K) = 627 + TEMP2 (R) 628 + * INT2DZ (R,ZE,ZC,ZD) 629 + + TEMP2 (R+1) 630 + * INT2DZ (R+1,ZE,ZC,ZD) 631 + + TEMP2 (R+2) 632 + * INT2DZ (R+2,ZE,ZC,ZD) 633 + + TEMP2 (R+3) 634 + * INT2DZ (R+3,ZE,ZC,ZD) 635 + + TEMP2 (R+4) 636 + * INT2DZ (R+4,ZE,ZC,ZD) 637 R = R + 5 638 END DO 639 END IF 640 641 I = I - NXYZE + XE 642 NXYZE = NXYZE - SE - 1 643 520 CONTINUE 644 514 CONTINUE 645 512 CONTINUE 646 510 CONTINUE 647 504 CONTINUE 648 XCP = J 649 502 CONTINUE 650 XDP = K 651 500 CONTINUE 652 653 RETURN 654C 655C 656C ******************** 657C * # of roots = 6 * 658C ******************** 659C 660C 661 6 XDP = 0 662 DO 600 XD = SHELLD,0,-1 663 YDMAX = SHELLD - XD 664 XCP = 0 665 DO 602 XC = SHELLC,0,-1 666 YCMAX = SHELLC - XC 667 668 XEP = NXYZET + 3 669 DO 604 XE = 0,SHELLP 670 XEP = XEP + XE - 2 671 XEMAX = XE * SHELLP 672 YEEND = SHELLP - XE 673 674 DO M = 1,NGQEXQ 675 TEMP1 (M) = SCALE (M) * INT2DX (M,XE,XC,XD) 676 END DO 677 678 K = XDP 679 DO 610 YD = YDMAX,0,-1 680 K = K + 1 681 ZD = YDMAX - YD 682 J = XCP 683 DO 612 YC = YCMAX,0,-1 684 J = J + 1 685 ZC = YCMAX - YC 686 687 XYEP = XEP - XEMAX 688 DO 614 YE = 0,YEEND 689 XYE = XE + YE 690 XYEP = XYEP - 1 691 SEEND = MAX0 (SHELLA,XYE) 692 693 IF (.NOT.DIFFY .AND. YE+YC+YD.EQ.0) THEN 694 DO N = 1,NGQEXQ 695 TEMP2 (N) = TEMP1 (N) 696 END DO 697 ELSE 698 DO N = 1,NGQEXQ 699 TEMP2 (N) = TEMP1 (N) 700 + * INT2DY (N,YE,YC,YD) 701 END DO 702 END IF 703 704 I = XYEP 705 NXYZE = NXYZP 706 DO 620 SE = SHELLP,SEEND,-1 707 ZE = SE - XYE 708 709 IF (.NOT.DIFFZ .AND. ZE+ZC+ZD.EQ.0) THEN 710 R = 1 711 DO M = 1,NEXQ 712 BATCH (M,I,J,K) = TEMP2 (R) 713 + + TEMP2 (R+1) 714 + + TEMP2 (R+2) 715 + + TEMP2 (R+3) 716 + + TEMP2 (R+4) 717 + + TEMP2 (R+5) 718 R = R + 6 719 END DO 720 ELSE 721 R = 1 722 DO M = 1,NEXQ 723 BATCH (M,I,J,K) = 724 + TEMP2 (R) 725 + * INT2DZ (R,ZE,ZC,ZD) 726 + + TEMP2 (R+1) 727 + * INT2DZ (R+1,ZE,ZC,ZD) 728 + + TEMP2 (R+2) 729 + * INT2DZ (R+2,ZE,ZC,ZD) 730 + + TEMP2 (R+3) 731 + * INT2DZ (R+3,ZE,ZC,ZD) 732 + + TEMP2 (R+4) 733 + * INT2DZ (R+4,ZE,ZC,ZD) 734 + + TEMP2 (R+5) 735 + * INT2DZ (R+5,ZE,ZC,ZD) 736 R = R + 6 737 END DO 738 END IF 739 740 I = I - NXYZE + XE 741 NXYZE = NXYZE - SE - 1 742 620 CONTINUE 743 614 CONTINUE 744 612 CONTINUE 745 610 CONTINUE 746 604 CONTINUE 747 XCP = J 748 602 CONTINUE 749 XDP = K 750 600 CONTINUE 751 752 RETURN 753C 754C 755C ******************** 756C * # of roots = 7 * 757C ******************** 758C 759C 760 7 XDP = 0 761 DO 700 XD = SHELLD,0,-1 762 YDMAX = SHELLD - XD 763 XCP = 0 764 DO 702 XC = SHELLC,0,-1 765 YCMAX = SHELLC - XC 766 767 XEP = NXYZET + 3 768 DO 704 XE = 0,SHELLP 769 XEP = XEP + XE - 2 770 XEMAX = XE * SHELLP 771 YEEND = SHELLP - XE 772 773 DO M = 1,NGQEXQ 774 TEMP1 (M) = SCALE (M) * INT2DX (M,XE,XC,XD) 775 END DO 776 777 K = XDP 778 DO 710 YD = YDMAX,0,-1 779 K = K + 1 780 ZD = YDMAX - YD 781 J = XCP 782 DO 712 YC = YCMAX,0,-1 783 J = J + 1 784 ZC = YCMAX - YC 785 786 XYEP = XEP - XEMAX 787 DO 714 YE = 0,YEEND 788 XYE = XE + YE 789 XYEP = XYEP - 1 790 SEEND = MAX0 (SHELLA,XYE) 791 792 IF (.NOT.DIFFY .AND. YE+YC+YD.EQ.0) THEN 793 DO N = 1,NGQEXQ 794 TEMP2 (N) = TEMP1 (N) 795 END DO 796 ELSE 797 DO N = 1,NGQEXQ 798 TEMP2 (N) = TEMP1 (N) 799 + * INT2DY (N,YE,YC,YD) 800 END DO 801 END IF 802 803 I = XYEP 804 NXYZE = NXYZP 805 DO 720 SE = SHELLP,SEEND,-1 806 ZE = SE - XYE 807 808 IF (.NOT.DIFFZ .AND. ZE+ZC+ZD.EQ.0) THEN 809 R = 1 810 DO M = 1,NEXQ 811 BATCH (M,I,J,K) = TEMP2 (R) 812 + + TEMP2 (R+1) 813 + + TEMP2 (R+2) 814 + + TEMP2 (R+3) 815 + + TEMP2 (R+4) 816 + + TEMP2 (R+5) 817 + + TEMP2 (R+6) 818 R = R + 7 819 END DO 820 ELSE 821 R = 1 822 DO M = 1,NEXQ 823 BATCH (M,I,J,K) = 824 + TEMP2 (R) 825 + * INT2DZ (R,ZE,ZC,ZD) 826 + + TEMP2 (R+1) 827 + * INT2DZ (R+1,ZE,ZC,ZD) 828 + + TEMP2 (R+2) 829 + * INT2DZ (R+2,ZE,ZC,ZD) 830 + + TEMP2 (R+3) 831 + * INT2DZ (R+3,ZE,ZC,ZD) 832 + + TEMP2 (R+4) 833 + * INT2DZ (R+4,ZE,ZC,ZD) 834 + + TEMP2 (R+5) 835 + * INT2DZ (R+5,ZE,ZC,ZD) 836 + + TEMP2 (R+6) 837 + * INT2DZ (R+6,ZE,ZC,ZD) 838 R = R + 7 839 END DO 840 END IF 841 842 I = I - NXYZE + XE 843 NXYZE = NXYZE - SE - 1 844 720 CONTINUE 845 714 CONTINUE 846 712 CONTINUE 847 710 CONTINUE 848 704 CONTINUE 849 XCP = J 850 702 CONTINUE 851 XDP = K 852 700 CONTINUE 853 854 RETURN 855C 856C 857C ******************** 858C * # of roots = 8 * 859C ******************** 860C 861C 862 8 XDP = 0 863 DO 800 XD = SHELLD,0,-1 864 YDMAX = SHELLD - XD 865 XCP = 0 866 DO 802 XC = SHELLC,0,-1 867 YCMAX = SHELLC - XC 868 869 XEP = NXYZET + 3 870 DO 804 XE = 0,SHELLP 871 XEP = XEP + XE - 2 872 XEMAX = XE * SHELLP 873 YEEND = SHELLP - XE 874 875 DO M = 1,NGQEXQ 876 TEMP1 (M) = SCALE (M) * INT2DX (M,XE,XC,XD) 877 END DO 878 879 K = XDP 880 DO 810 YD = YDMAX,0,-1 881 K = K + 1 882 ZD = YDMAX - YD 883 J = XCP 884 DO 812 YC = YCMAX,0,-1 885 J = J + 1 886 ZC = YCMAX - YC 887 888 XYEP = XEP - XEMAX 889 DO 814 YE = 0,YEEND 890 XYE = XE + YE 891 XYEP = XYEP - 1 892 SEEND = MAX0 (SHELLA,XYE) 893 894 IF (.NOT.DIFFY .AND. YE+YC+YD.EQ.0) THEN 895 DO N = 1,NGQEXQ 896 TEMP2 (N) = TEMP1 (N) 897 END DO 898 ELSE 899 DO N = 1,NGQEXQ 900 TEMP2 (N) = TEMP1 (N) 901 + * INT2DY (N,YE,YC,YD) 902 END DO 903 END IF 904 905 I = XYEP 906 NXYZE = NXYZP 907 DO 820 SE = SHELLP,SEEND,-1 908 ZE = SE - XYE 909 910 IF (.NOT.DIFFZ .AND. ZE+ZC+ZD.EQ.0) THEN 911 R = 1 912 DO M = 1,NEXQ 913 BATCH (M,I,J,K) = TEMP2 (R) 914 + + TEMP2 (R+1) 915 + + TEMP2 (R+2) 916 + + TEMP2 (R+3) 917 + + TEMP2 (R+4) 918 + + TEMP2 (R+5) 919 + + TEMP2 (R+6) 920 + + TEMP2 (R+7) 921 R = R + 8 922 END DO 923 ELSE 924 R = 1 925 DO M = 1,NEXQ 926 BATCH (M,I,J,K) = 927 + TEMP2 (R) 928 + * INT2DZ (R,ZE,ZC,ZD) 929 + + TEMP2 (R+1) 930 + * INT2DZ (R+1,ZE,ZC,ZD) 931 + + TEMP2 (R+2) 932 + * INT2DZ (R+2,ZE,ZC,ZD) 933 + + TEMP2 (R+3) 934 + * INT2DZ (R+3,ZE,ZC,ZD) 935 + + TEMP2 (R+4) 936 + * INT2DZ (R+4,ZE,ZC,ZD) 937 + + TEMP2 (R+5) 938 + * INT2DZ (R+5,ZE,ZC,ZD) 939 + + TEMP2 (R+6) 940 + * INT2DZ (R+6,ZE,ZC,ZD) 941 + + TEMP2 (R+7) 942 + * INT2DZ (R+7,ZE,ZC,ZD) 943 R = R + 8 944 END DO 945 END IF 946 947 I = I - NXYZE + XE 948 NXYZE = NXYZE - SE - 1 949 820 CONTINUE 950 814 CONTINUE 951 812 CONTINUE 952 810 CONTINUE 953 804 CONTINUE 954 XCP = J 955 802 CONTINUE 956 XDP = K 957 800 CONTINUE 958 959 RETURN 960C 961C 962C ******************** 963C * # of roots = 9 * 964C ******************** 965C 966C 967 9 XDP = 0 968 DO 900 XD = SHELLD,0,-1 969 YDMAX = SHELLD - XD 970 XCP = 0 971 DO 902 XC = SHELLC,0,-1 972 YCMAX = SHELLC - XC 973 974 XEP = NXYZET + 3 975 DO 904 XE = 0,SHELLP 976 XEP = XEP + XE - 2 977 XEMAX = XE * SHELLP 978 YEEND = SHELLP - XE 979 980 DO M = 1,NGQEXQ 981 TEMP1 (M) = SCALE (M) * INT2DX (M,XE,XC,XD) 982 END DO 983 984 K = XDP 985 DO 910 YD = YDMAX,0,-1 986 K = K + 1 987 ZD = YDMAX - YD 988 J = XCP 989 DO 912 YC = YCMAX,0,-1 990 J = J + 1 991 ZC = YCMAX - YC 992 993 XYEP = XEP - XEMAX 994 DO 914 YE = 0,YEEND 995 XYE = XE + YE 996 XYEP = XYEP - 1 997 SEEND = MAX0 (SHELLA,XYE) 998 999 IF (.NOT.DIFFY .AND. YE+YC+YD.EQ.0) THEN 1000 DO N = 1,NGQEXQ 1001 TEMP2 (N) = TEMP1 (N) 1002 END DO 1003 ELSE 1004 DO N = 1,NGQEXQ 1005 TEMP2 (N) = TEMP1 (N) 1006 + * INT2DY (N,YE,YC,YD) 1007 END DO 1008 END IF 1009 1010 I = XYEP 1011 NXYZE = NXYZP 1012 DO 920 SE = SHELLP,SEEND,-1 1013 ZE = SE - XYE 1014 1015 IF (.NOT.DIFFZ .AND. ZE+ZC+ZD.EQ.0) THEN 1016 R = 1 1017 DO M = 1,NEXQ 1018 BATCH (M,I,J,K) = TEMP2 (R) 1019 + + TEMP2 (R+1) 1020 + + TEMP2 (R+2) 1021 + + TEMP2 (R+3) 1022 + + TEMP2 (R+4) 1023 + + TEMP2 (R+5) 1024 + + TEMP2 (R+6) 1025 + + TEMP2 (R+7) 1026 + + TEMP2 (R+8) 1027 R = R + 9 1028 END DO 1029 ELSE 1030 R = 1 1031 DO M = 1,NEXQ 1032 BATCH (M,I,J,K) = 1033 + TEMP2 (R) 1034 + * INT2DZ (R,ZE,ZC,ZD) 1035 + + TEMP2 (R+1) 1036 + * INT2DZ (R+1,ZE,ZC,ZD) 1037 + + TEMP2 (R+2) 1038 + * INT2DZ (R+2,ZE,ZC,ZD) 1039 + + TEMP2 (R+3) 1040 + * INT2DZ (R+3,ZE,ZC,ZD) 1041 + + TEMP2 (R+4) 1042 + * INT2DZ (R+4,ZE,ZC,ZD) 1043 + + TEMP2 (R+5) 1044 + * INT2DZ (R+5,ZE,ZC,ZD) 1045 + + TEMP2 (R+6) 1046 + * INT2DZ (R+6,ZE,ZC,ZD) 1047 + + TEMP2 (R+7) 1048 + * INT2DZ (R+7,ZE,ZC,ZD) 1049 + + TEMP2 (R+8) 1050 + * INT2DZ (R+8,ZE,ZC,ZD) 1051 R = R + 9 1052 END DO 1053 END IF 1054 1055 I = I - NXYZE + XE 1056 NXYZE = NXYZE - SE - 1 1057 920 CONTINUE 1058 914 CONTINUE 1059 912 CONTINUE 1060 910 CONTINUE 1061 904 CONTINUE 1062 XCP = J 1063 902 CONTINUE 1064 XDP = K 1065 900 CONTINUE 1066 1067 RETURN 1068C 1069C 1070C ******************** 1071C * # of roots > 9 * 1072C ******************** 1073C 1074C ...outer loops over x,x,x-triples. No skipping of the 1075C x,x,x-contribution of 0,0,0-type can be done here, 1076C since the 2DX integrals carry the Rys weight! 1077C 1078C 1079 10 XDP = 0 1080 DO 1000 XD = SHELLD,0,-1 1081 YDMAX = SHELLD - XD 1082 XCP = 0 1083 DO 1002 XC = SHELLC,0,-1 1084 YCMAX = SHELLC - XC 1085 1086 XEP = NXYZET + 3 1087 DO 1004 XE = 0,SHELLP 1088 XEP = XEP + XE - 2 1089 XEMAX = XE * SHELLP 1090 YEEND = SHELLP - XE 1091 1092 DO M = 1,NGQEXQ 1093 TEMP1 (M) = SCALE (M) * INT2DX (M,XE,XC,XD) 1094 END DO 1095C 1096C 1097C ...middle loops over y,y,y-triples. Skip multiplication 1098C of y,y,y-contributions, if no y-coordinate derivative 1099C was formed and we have a 0,0,0-triple, as then the 2DY 1100C PCD integrals are equal to 1. 1101C 1102C 1103 K = XDP 1104 DO 1010 YD = YDMAX,0,-1 1105 K = K + 1 1106 ZD = YDMAX - YD 1107 J = XCP 1108 DO 1012 YC = YCMAX,0,-1 1109 J = J + 1 1110 ZC = YCMAX - YC 1111 1112 XYEP = XEP - XEMAX 1113 DO 1014 YE = 0,YEEND 1114 XYE = XE + YE 1115 XYEP = XYEP - 1 1116 SEEND = MAX0 (SHELLA,XYE) 1117 1118 IF (.NOT.DIFFY .AND. YE+YC+YD.EQ.0) THEN 1119 DO N = 1,NGQEXQ 1120 TEMP2 (N) = TEMP1 (N) 1121 END DO 1122 ELSE 1123 DO N = 1,NGQEXQ 1124 TEMP2 (N) = TEMP1 (N) 1125 + * INT2DY (N,YE,YC,YD) 1126 END DO 1127 END IF 1128C 1129C 1130C ...inner loop over E shells. Skip multiplication of 1131C z,z,z-contributions, if we have a 0,0,0-triple and 1132C no derivations were performed on the z-coordinate, 1133C as then the 2DZ PCD integrals are equal to 1. 1134C All info concerning all three x,x,x-, y,y,y- and 1135C z,z,z-triples have been collected for all exponent 1136C quadruplets at once. Sum up the 2D X,Y,Z integral 1137C products to the appropriate place of the [E0|CD] 1138C batch. 1139C 1140C 1141 I = XYEP 1142 NXYZE = NXYZP 1143 DO 1020 SE = SHELLP,SEEND,-1 1144 ZE = SE - XYE 1145 1146 IF (.NOT.DIFFZ .AND. ZE+ZC+ZD.EQ.0) THEN 1147 R = 0 1148 DO M = 1,NEXQ 1149 SUM = ZERO 1150 DO N = 1,NGQP 1151 SUM = SUM + TEMP2 (R+N) 1152 END DO 1153 R = R + NGQP 1154 BATCH (M,I,J,K) = SUM 1155 END DO 1156 ELSE 1157 R = 0 1158 DO M = 1,NEXQ 1159 SUM = ZERO 1160 DO N = 1,NGQP 1161 SUM = SUM + TEMP2 (R+N) 1162 + * INT2DZ (R+N,ZE,ZC,ZD) 1163 END DO 1164 R = R + NGQP 1165 BATCH (M,I,J,K) = SUM 1166 END DO 1167 END IF 1168C 1169C 1170C ...next z,z,z-triple. 1171C 1172C 1173 I = I - NXYZE + XE 1174 NXYZE = NXYZE - SE - 1 1175 1020 CONTINUE 1176C 1177C 1178C ...next y,y,y-triple and next x,x,x-triple. 1179C 1180C 1181 1014 CONTINUE 1182 1012 CONTINUE 1183 1010 CONTINUE 1184 1185 1004 CONTINUE 1186 XCP = J 1187 1002 CONTINUE 1188 XDP = K 1189 1000 CONTINUE 1190C 1191C 1192C ...ready! 1193C 1194C 1195 RETURN 1196 END 1197