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