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