1 //
2 // comp2e.cc
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Curtis Janssen <cljanss@limitpt.com>
7 // Maintainer: LPS
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27 
28 #include <stdarg.h>
29 
30 #include <util/misc/formio.h>
31 #include <chemistry/qc/intv3/macros.h>
32 #include <chemistry/qc/intv3/flags.h>
33 #include <chemistry/qc/intv3/types.h>
34 #include <chemistry/qc/intv3/int2e.h>
35 #include <chemistry/qc/intv3/utils.h>
36 #include <chemistry/qc/intv3/tformv3.h>
37 
38 using namespace std;
39 using namespace sc;
40 
41 #undef DER_TIMING
42 #undef EREP_TIMING
43 
44 #if defined(DER_TIMING)||defined(EREP_TIMING)
45 #  include <util/misc/timer.h>
46 #endif
47 
48 static inline void
swtch(GaussianBasisSet * & i,GaussianBasisSet * & j)49 swtch(GaussianBasisSet* &i,GaussianBasisSet* &j)
50 {
51   GaussianBasisSet *tmp;
52   tmp = i;
53   i = j;
54   j = tmp;
55 }
56 
57 static inline void
sswtch(GaussianShell ** i,GaussianShell ** j)58 sswtch(GaussianShell**i,GaussianShell**j)
59 {
60   GaussianShell*tmp;
61   tmp = *i;
62   *i = *j;
63   *j = tmp;
64 }
65 
66 static inline void
iswtch(int * i,int * j)67 iswtch(int *i,int *j)
68 {
69   int tmp;
70   tmp = *i;
71   *i = *j;
72   *j = tmp;
73 }
74 
75 static void
fail()76 fail()
77 {
78   ExEnv::errn() << scprintf("failing module:\n%s",__FILE__) << endl;
79   abort();
80 }
81 
82 /* This computes the 2erep integrals for a shell quartet
83  * specified by psh1, psh2, psh3, psh4.
84  * The routine int_initialize_2erep must be called before
85  * any integrals can be computed.
86  * This routine may decide to change the shell ordering.
87  * The new ordering is placed in *psh1,4 on exit.
88  * for the derivatives.
89  */
90 void
erep(int & psh1,int & psh2,int & psh3,int & psh4)91 Int2eV3::erep(int &psh1, int &psh2, int &psh3, int &psh4)
92 {
93   compute_erep(0,&psh1,&psh2,&psh3,&psh4,0,0,0,0);
94   }
95 
96 /* This is an alternate interface to int_erep.  It takes
97  * as arguments the flags, an integer vector of shell numbers
98  * and an integer vector which will be filled in with size
99  * information, if it is non-NULL. */
100 void
erep(int * shells,int * sizes)101 Int2eV3::erep(int *shells, int  *sizes)
102 {
103   erep(shells[0],shells[1],shells[2],shells[3]);
104   if (sizes) {
105     sizes[0] = bs1_->shell(shells[0]).nfunction();
106     sizes[1] = bs2_->shell(shells[1]).nfunction();
107     sizes[2] = bs3_->shell(shells[2]).nfunction();
108     sizes[3] = bs4_->shell(shells[3]).nfunction();
109     }
110   }
111 
112 /* If we need a computation with adjusted angular momentum, then
113  * this lower level routine can be called instead of int_erep.
114  * The dam{1,2,3,4} arguments given the amount by which the
115  * angular momentum is to adjusted.  This differs from libint version
116  * 1 in that it used total angular momentum here.
117  */
118 void
compute_erep(int flags,int * psh1,int * psh2,int * psh3,int * psh4,int dam1,int dam2,int dam3,int dam4)119 Int2eV3::compute_erep(int flags, int *psh1, int *psh2, int *psh3, int *psh4,
120                       int dam1, int dam2, int dam3, int dam4)
121 {
122 #ifdef EREP_TIMING
123   char section[30];
124 #endif
125   GaussianBasisSet *pbs1=bs1_.pointer();
126   GaussianBasisSet *pbs2=bs2_.pointer();
127   GaussianBasisSet *pbs3=bs3_.pointer();
128   GaussianBasisSet *pbs4=bs4_.pointer();
129   int size;
130   int ii;
131   int size1, size2, size3, size4;
132   int tam1,tam2,tam3,tam4;
133   int i,j,k,l;
134   int ogc1,ogc2,ogc3,ogc4;
135   int sh1,sh2,sh3,sh4;
136   int am1,am2,am3,am4,am12, am34;
137   int minam1,minam2,minam3,minam4;
138   int redundant_index;
139   int e12,e13e24,e34;
140   int p12,p34,p13p24;
141   int eAB;
142 
143   /* Compute the offset shell numbers. */
144   osh1 = *psh1 + bs1_shell_offset_;
145   osh2 = *psh2 + bs2_shell_offset_;
146   osh3 = *psh3 + bs3_shell_offset_;
147   osh4 = *psh4 + bs4_shell_offset_;
148 
149   sh1 = *psh1;
150   sh2 = *psh2;
151   sh3 = *psh3;
152   sh4 = *psh4;
153 
154   /* Test the arguments to make sure that they are sensible. */
155   if (   sh1 < 0 || sh1 >= bs1_->nbasis()
156       ||( !int_unit2 && (sh2 < 0 || sh2 >= bs2_->nbasis()))
157       || sh3 < 0 || sh3 >= bs3_->nbasis()
158       ||( !int_unit4 && (sh4 < 0 || sh4 >= bs4_->nbasis()))) {
159     ExEnv::errn() << scprintf("compute_erep has been incorrectly used\n");
160     ExEnv::errn() << scprintf("shells (bounds): %d (%d), %d (%d), %d (%d), %d (%d)\n",
161             sh1,bs1_->nbasis()-1,
162             sh2,(bs2_.null()?0:bs2_->nbasis())-1,
163             sh3,bs3_->nbasis()-1,
164             sh4,(bs4_.null()?0:bs4_->nbasis())-1);
165     fail();
166     }
167 
168   /* Set up pointers to the current shells. */
169   int_shell1 = &bs1_->shell(sh1);
170   if (!int_unit2) int_shell2 = &bs2_->shell(sh2);
171   else int_shell2 = int_unit_shell;
172   int_shell3 = &bs3_->shell(sh3);
173   if (!int_unit4) int_shell4 = &bs4_->shell(sh4);
174   else int_shell4 = int_unit_shell;
175 
176 
177   /* Compute the maximum angular momentum on each centers to
178    * determine the most efficient way to invoke the building and shifting
179    * routines.  The minimum angular momentum will be computed at the
180    * same time. */
181   minam1 = int_shell1->min_am();
182   minam2 = int_shell2->min_am();
183   minam3 = int_shell3->min_am();
184   minam4 = int_shell4->min_am();
185   am1 = int_shell1->max_am();
186   am2 = int_shell2->max_am();
187   am3 = int_shell3->max_am();
188   am4 = int_shell4->max_am();
189 
190   am1 += dam1; minam1 += dam1;
191   am2 += dam2; minam2 += dam2;
192   am3 += dam3; minam3 += dam3;
193   am4 += dam4; minam4 += dam4;
194   am12 = am1 + am2;
195   am34 = am3 + am4;
196 
197   /* if no angular momentum remains on one of the centers return */
198   if (am1 < 0 || am2 < 0 || am3 < 0 || am4 < 0) {
199     return;
200     }
201 
202 #ifdef EREP_TIMING
203   sprintf(section,"erep am=%02d",am12+am34);
204   tim_enter(section);
205   tim_enter("setup");
206 #endif
207 
208   /* Convert the integral to the most efficient form. */
209   p12 = 0;
210   p34 = 0;
211   p13p24 = 0;
212 
213   if (am2 > am1) {
214     p12 = 1;
215     iswtch(&am1,&am2);iswtch(&sh1,&sh2);iswtch(psh1,psh2);iswtch(&osh1,&osh2);
216     iswtch(&dam1,&dam2);
217     iswtch(&minam1,&minam2);
218     sswtch(&int_shell1,&int_shell2);
219     swtch(pbs1,pbs2);
220     }
221   if (am4 > am3) {
222     p34 = 1;
223     iswtch(&am3,&am4);iswtch(&sh3,&sh4);iswtch(psh3,psh4);iswtch(&osh3,&osh4);
224     iswtch(&dam3,&dam4);
225     iswtch(&minam3,&minam4);
226     sswtch(&int_shell3,&int_shell4);
227     swtch(pbs3,pbs4);
228     }
229   if ((osh1 == osh4) && (osh2 == osh3) && (osh1 != osh2)) {
230     /* Don't make the permutation unless we won't override what was
231      * decided above about p34. */
232     if (am4 == am3) {
233       p34 = 1;
234       iswtch(&am3,&am4);iswtch(&sh3,&sh4);iswtch(psh3,psh4);iswtch(&osh3,&osh4);
235       iswtch(&dam3,&dam4);
236       iswtch(&minam3,&minam4);
237       sswtch(&int_shell3,&int_shell4);
238       swtch(pbs3,pbs4);
239       }
240     }
241   if ((am34 > am12)||((am34 == am12)&&(minam1 > minam3))) {
242     p13p24 = 1;
243     iswtch(&am1,&am3);iswtch(&sh1,&sh3);iswtch(psh1,psh3);iswtch(&osh1,&osh3);
244     iswtch(&am2,&am4);iswtch(&sh2,&sh4);iswtch(psh2,psh4);iswtch(&osh2,&osh4);
245     iswtch(&int_unit2,&int_unit4);
246     iswtch(&am12,&am34);
247     iswtch(&dam1,&dam3);
248     iswtch(&minam1,&minam3);
249     sswtch(&int_shell1,&int_shell3);
250     swtch(pbs1,pbs3);
251     iswtch(&dam2,&dam4);
252     iswtch(&minam2,&minam4);
253     sswtch(&int_shell2,&int_shell4);
254     swtch(pbs2,pbs4);
255     }
256   /* This tries to make centers A and B equivalent, if possible. */
257   else if (  (am3 == am1)
258            &&(am4 == am2)
259            && !int_unit2
260            && !int_unit4
261            &&(minam1 == minam3)
262            &&(!(  (bs1_ == bs2_)
263                 &&(bs1_->shell_to_center(sh1)==bs2_->shell_to_center(sh2))))
264            &&(   (bs3_ == bs4_)
265                &&(bs3_->shell_to_center(sh3)==bs4_->shell_to_center(sh4)))) {
266     p13p24 = 1;
267     iswtch(&am1,&am3);iswtch(&sh1,&sh3);iswtch(psh1,psh3);iswtch(&osh1,&osh3);
268     iswtch(&am2,&am4);iswtch(&sh2,&sh4);iswtch(psh2,psh4);iswtch(&osh2,&osh4);
269     iswtch(&am12,&am34);
270     iswtch(&dam1,&dam3);
271     iswtch(&minam1,&minam3);
272     sswtch(&int_shell1,&int_shell3);
273     swtch(pbs1,pbs3);
274     iswtch(&dam2,&dam4);
275     iswtch(&minam2,&minam4);
276     sswtch(&int_shell2,&int_shell4);
277     swtch(pbs2,pbs4);
278     }
279 
280   if ((pbs1 == pbs2)
281       &&(pbs1->shell_to_center(sh1)==pbs2->shell_to_center(sh2))) {
282     eAB = 1;
283     }
284   else {
285     eAB = 0;
286     }
287 
288   /* If the centers were permuted, then the int_expweighted variable may
289    * need to be changed. */
290   if (p12) {
291     iswtch(&int_expweight1,&int_expweight2);
292     }
293   if (p34) {
294     iswtch(&int_expweight3,&int_expweight4);
295     }
296   if (p13p24) {
297     iswtch(&int_expweight1,&int_expweight3);
298     iswtch(&int_expweight2,&int_expweight4);
299     }
300 
301   pbs1_ = pbs1;
302   pbs2_ = pbs2;
303   pbs3_ = pbs3;
304   pbs4_ = pbs4;
305 
306   int nc1 = int_shell1->ncontraction();
307   int nc2 = int_shell2->ncontraction();
308   int nc3 = int_shell3->ncontraction();
309   int nc4 = int_shell4->ncontraction();
310 
311   /* Compute the shell sizes. */
312   for (ii=size1=0; ii<nc1; ii++)
313     size1 += INT_NCART(int_shell1->am(ii)+dam1);
314   for (ii=size2=0; ii<nc2; ii++)
315     size2 += INT_NCART(int_shell2->am(ii)+dam2);
316   for (ii=size3=0; ii<nc3; ii++)
317     size3 += INT_NCART(int_shell3->am(ii)+dam3);
318   for (ii=size4=0; ii<nc4; ii++)
319     size4 += INT_NCART(int_shell4->am(ii)+dam4);
320   size = size1*size2*size3*size4;
321 
322   if (int_integral_storage) {
323 #ifdef EREP_TIMING
324       tim_change("check storage");
325 #endif
326     if (dam1 || dam2 || dam3 || dam4) {
327       ExEnv::errn() << scprintf("cannot use integral storage and dam\n");
328       fail();
329       }
330     if (    !int_unit2
331          && !int_unit4
332          && int_have_stored_integral(sh1,sh2,sh3,sh4,p12,p34,p13p24)) {
333       goto post_computation;
334       }
335     }
336 
337   /* Buildam up on center 1 and 3. */
338 #ifdef EREP_TIMING
339   tim_change("build");
340 #endif
341   int_buildgcam(minam1,minam2,minam3,minam4,
342                 am1,am2,am3,am4,
343                 dam1,dam2,dam3,dam4,
344                 sh1,sh2,sh3,sh4,
345                 eAB);
346 #ifdef EREP_TIMING
347   tim_change("cleanup");
348 #endif
349 
350   /* Begin loop over generalized contractions. */
351   ogc1 = 0;
352   for (i=0; i<nc1; i++) {
353     tam1 = int_shell1->am(i) + dam1;
354     if (tam1 < 0) continue;
355     int tsize1 = INT_NCART_NN(tam1);
356     ogc2 = 0;
357     for (j=0; j<nc2; j++) {
358       tam2 = int_shell2->am(j) + dam2;
359       if (tam2 < 0) continue;
360       int tsize2 = INT_NCART_NN(tam2);
361       ogc3 = 0;
362       for (k=0; k<nc3; k++) {
363         tam3 = int_shell3->am(k) + dam3;
364         if (tam3 < 0) continue;
365         int tsize3 = INT_NCART_NN(tam3);
366         ogc4 = 0;
367         for (l=0; l<nc4; l++) {
368           tam4 = int_shell4->am(l) + dam4;
369           if (tam4 < 0) continue;
370           int tsize4 = INT_NCART_NN(tam4);
371 
372 #ifdef EREP_TIMING
373   tim_change("shift");
374 #endif
375   /* Shift angular momentum from 1 to 2 and from 3 to 4. */
376   double *shiftbuffer = int_shiftgcam(i,j,k,l,tam1,tam2,tam3,tam4, eAB);
377 #ifdef EREP_TIMING
378   tim_change("cleanup");
379 #endif
380 
381   /* Place the integrals in the integral buffer. */
382   /* If permute_ is not set, then repack the integrals while copying. */
383   if ((!permute_)&&(p12||p34||p13p24)) {
384     int pam1,pam2,pam3,pam4;
385     int psize234,psize34;
386     int pogc1,pogc2,pogc3,pogc4;
387     int psize1,psize2,psize3,psize4;
388     pam1 = tam1;
389     pam2 = tam2;
390     pam3 = tam3;
391     pam4 = tam4;
392     pogc1 = ogc1;
393     pogc2 = ogc2;
394     pogc3 = ogc3;
395     pogc4 = ogc4;
396     psize1 = size1;
397     psize2 = size2;
398     psize3 = size3;
399     psize4 = size4;
400     if (p13p24) {
401       iswtch(&pam1,&pam3);
402       iswtch(&pam2,&pam4);
403       iswtch(&pogc1,&pogc3);
404       iswtch(&pogc2,&pogc4);
405       iswtch(&psize1,&psize3);
406       iswtch(&psize2,&psize4);
407       }
408     if (p34) {
409       iswtch(&pam3,&pam4);
410       iswtch(&pogc3,&pogc4);
411       iswtch(&psize3,&psize4);
412       }
413     if (p12) {
414       iswtch(&pam1,&pam2);
415       iswtch(&pogc1,&pogc2);
416       iswtch(&psize1,&psize2);
417       }
418     psize34 = psize4 * psize3;
419     psize234 = psize34 * psize2;
420     redundant_index = 0;
421     int newindexoffset = pogc1*psize234 + pogc2*psize34 + pogc3*psize4 + pogc4;
422     if (p13p24||p34) {
423       int stride1=psize234;
424       int stride2=psize34;
425       int stride3=psize4;
426       int stride4=1;
427       int tmp;
428       if (p12) {
429         tmp=stride1; stride1=stride2; stride2=tmp;
430         }
431       if (p34) {
432         tmp=stride3; stride3=stride4; stride4=tmp;
433         }
434       if (p13p24) {
435         tmp=stride1; stride1=stride3; stride3=tmp;
436         tmp=stride2; stride2=stride4; stride4=tmp;
437         }
438       int newindex1 = newindexoffset;
439       for (int ci1=0; ci1<tsize1; ci1++) {
440         int newindex12 = newindex1;
441         for (int ci2=0; ci2<tsize2; ci2++) {
442           int newindex123 = newindex12;
443           for (int ci3=0; ci3<tsize3; ci3++) {
444             double *tmp_shiftbuffer = &shiftbuffer[redundant_index];
445             int newindex1234 = newindex123;
446             for (int ci4=0; ci4<tsize4; ci4++) {
447               int_buffer[newindex1234] = tmp_shiftbuffer[ci4];
448               newindex1234 += stride4;
449               }
450             redundant_index+=tsize4;
451             newindex123 += stride3;
452             }
453           newindex12 += stride2;
454           }
455         newindex1 += stride1;
456         }
457       }
458     else if (nc3 == 1 && nc4 == 1) {
459       // this is the p12 only case w/o gen contractions on 3 & 4
460       // this special case collapses the 3rd and 4th indices together
461       for (int ci1=0; ci1<tsize1; ci1++) {
462         for (int ci2=0; ci2<tsize2; ci2++) {
463           int pci1=ci1;
464           int pci2=ci2;
465           if (p12) {
466             int tmp=pci1; pci1=pci2; pci2=tmp;
467             }
468           int newindex123 = newindexoffset + pci1*psize234 + pci2*psize34;
469           double *tmp_int_buffer = &int_buffer[newindex123];
470           double *tmp_shiftbuffer = &shiftbuffer[redundant_index];
471           for (int ci34=0; ci34<psize34; ci34++)
472             tmp_int_buffer[ci34] = tmp_shiftbuffer[ci34];
473           redundant_index += psize34;
474           }
475         }
476       }
477     else {
478       // this is the p12 only case w/gen. contr. on 3 & 4
479       for (int ci1=0; ci1<tsize1; ci1++) {
480         for (int ci2=0; ci2<tsize2; ci2++) {
481           int pci1=ci1;
482           int pci2=ci2;
483           if (p12) {
484             int tmp=pci1; pci1=pci2; pci2=tmp;
485             }
486           int newindex123 = newindexoffset + pci1*psize234 + pci2*psize34;
487           for (int ci3=0; ci3<tsize3; ci3++) {
488             double *tmp_int_buffer = &int_buffer[newindex123];
489             double *tmp_shiftbuffer = &shiftbuffer[redundant_index];
490             for (int ci4=0; ci4<tsize4; ci4++) {
491               tmp_int_buffer[ci4] = tmp_shiftbuffer[ci4];
492               }
493             redundant_index += tsize4;
494             newindex123 += psize4;
495             }
496           }
497         }
498       }
499     }
500   else if (nc3 == 1 && nc4 == 1) {
501     // this special case collapses the 3rd and 4th indices together
502     int size34 =  size3 * size4;
503     int size234 = size2 * size34;
504     double* redund_ints = shiftbuffer;
505     redundant_index = 0;
506     int newindex1 = ogc1*size234 + ogc2*size34 + ogc3*size4 + ogc4;
507     for (int ci1=0; ci1<tsize1; ci1++) {
508       int newindex12 = newindex1;
509       for (int ci2=0; ci2<tsize2; ci2++) {
510         double *tmp_int_buffer = &int_buffer[newindex12];
511         double *tmp_redund_ints = &redund_ints[redundant_index];
512         for (int ci34=0; ci34<size34; ci34++)
513           tmp_int_buffer[ci34] = tmp_redund_ints[ci34];
514         redundant_index += size34;
515         newindex12 += size34;
516         }
517       newindex1 += size234;
518       }
519     }
520   else {
521     int size34 =  size3 * size4;
522     int size234 = size2 * size34;
523     double* redund_ints = shiftbuffer;
524     redundant_index = 0;
525     int newindex1 = ogc1*size234 + ogc2*size34 + ogc3*size4 + ogc4;
526     for (int ci1=0; ci1<tsize1; ci1++) {
527       int newindex12 = newindex1;
528       for (int ci2=0; ci2<tsize2; ci2++) {
529         int newindex123 = newindex12;
530         for (int ci3=0; ci3<tsize3; ci3++) {
531           double *tmp_int_buffer = &int_buffer[newindex123];
532           double *tmp_redund_ints = &redund_ints[redundant_index];
533           for (int ci4=0; ci4<tsize4; ci4++) {
534             tmp_int_buffer[ci4] = tmp_redund_ints[ci4];
535             }
536           redundant_index += tsize4;
537           newindex123 += size4;
538           }
539         newindex12 += size34;
540         }
541       newindex1 += size234;
542       }
543     }
544 
545     /* End loop over generalized contractions. */
546           ogc4 += tsize4;
547           }
548         ogc3 += tsize3;
549         }
550       ogc2 += tsize2;
551       }
552     ogc1 += tsize1;
553     }
554 
555   if (   !int_unit2
556       && !int_unit4
557       && int_integral_storage) {
558 #ifdef EREP_TIMING
559       tim_change("maybe store");
560 #endif
561       int_store_integral(sh1,sh2,sh3,sh4,p12,p34,p13p24,size);
562     }
563 
564   /* We branch here if an integral was precomputed and the int_buffer
565    * has been already filled. */
566   post_computation:
567 
568 #ifdef EREP_TIMING
569   tim_change("post");
570 #endif
571 
572   /* Unpermute all of the permuted quantities. */
573   if ((!permute_)&&(p12||p34||p13p24)) {
574     if (p13p24) {
575       iswtch(&sh1,&sh3);iswtch(psh1,psh3);iswtch(&osh1,&osh3);
576       iswtch(&sh2,&sh4);iswtch(psh2,psh4);iswtch(&osh2,&osh4);
577       iswtch(&int_unit2,&int_unit4);
578       iswtch(&am1,&am3);
579       iswtch(&am2,&am4);
580       iswtch(&am12,&am34);
581       sswtch(&int_shell1,&int_shell3);
582       swtch(pbs1,pbs3);
583       sswtch(&int_shell2,&int_shell4);
584       swtch(pbs2,pbs4);
585       iswtch(&int_expweight1,&int_expweight3);
586       iswtch(&int_expweight2,&int_expweight4);
587       }
588     if (p34) {
589       iswtch(&sh3,&sh4);iswtch(psh3,psh4);iswtch(&osh3,&osh4);
590       iswtch(&am3,&am4);
591       sswtch(&int_shell3,&int_shell4);
592       swtch(pbs3,pbs4);
593       iswtch(&int_expweight3,&int_expweight4);
594       }
595     if (p12) {
596       iswtch(&sh1,&sh2);iswtch(psh1,psh2);iswtch(&osh1,&osh2);
597       iswtch(&am1,&am2);
598       sswtch(&int_shell1,&int_shell2);
599       swtch(pbs1,pbs2);
600       iswtch(&int_expweight1,&int_expweight2);
601       }
602     }
603 
604   pbs1_ = 0;
605   pbs2_ = 0;
606   pbs3_ = 0;
607   pbs4_ = 0;
608 
609   /* Transform to pure am (if requested in the centers structure). */
610   if (!(flags&INT_NOPURE)) {
611       transform_2e(integral_, int_buffer, int_buffer,
612                    &bs1_->shell(sh1),
613                    int_unit2?int_unit_shell:&bs2_->shell(sh2),
614                    &bs3_->shell(sh3),
615                    int_unit4?int_unit_shell:&bs4_->shell(sh4));
616     }
617 
618   /* Remove the redundant integrals, unless redundant_ is set. */
619   if (!redundant_) {
620     int redundant_offset = 0;
621     int nonredundant_offset = 0;
622     if ((osh1 == osh4)&&(osh2 == osh3)&&(osh1 != osh2)) {
623       ExEnv::errn() << scprintf("nonredundant integrals cannot be generated\n");
624       fail();
625       }
626     e12 = (int_unit2?0:(osh1 == osh2));
627     e13e24 = ((osh1 == osh3)
628               && ((int_unit2 && int_unit4)
629                   || ((int_unit2||int_unit4)?0:(osh2 == osh4))));
630     e34 = (int_unit4?0:(osh3 == osh4));
631     if (e12||e34||e13e24) {
632       nonredundant_erep(int_buffer,e12,e34,e13e24,
633                         int_shell1->nfunction(),
634                         int_shell2->nfunction(),
635                         int_shell3->nfunction(),
636                         int_shell4->nfunction(),
637                         &redundant_offset,
638                         &nonredundant_offset);
639       }
640     }
641 
642 #ifdef EREP_TIMING
643   tim_exit("post");
644   tim_exit(section);
645 #endif
646   }
647 
648 /* This computes the two electron derivatives for all unique
649  * centers in the passed shell quartet.  One center in
650  * the set of unique centers is not included.  This can
651  * be computed as minus the sum of the other derivatives.
652  * The list of centers for which integrals were computed can
653  * be determined from the contents of der_centers.
654  * The results are put into the global integral buffer in the
655  * format:
656  * +------------------+
657  * | dercenter1 +---+ |
658  * |            | x | |
659  * |            +---+ |
660  * |            | y | |
661  * |            +---+ |
662  * |            | z | |
663  * |            +---+ |
664  * +------------------+
665  * | dercenter2 +---+ |
666  * |            | x | |
667  * |            +---+ |
668  * |            | y | |
669  * |            +---+ |
670  * |            | z | |
671  * |            +---+ |
672  * +------------------+
673  * | dercenter3 +---+ |
674  * |            | x | |
675  * |            +---+ |
676  * |            | y | |
677  * |            +---+ |
678  * |            | z | |
679  * |            +---+ |
680  * +------------------+
681  */
682 
683 void
erep_all1der(int & psh1,int & psh2,int & psh3,int & psh4,der_centersv3_t * der_centers)684 Int2eV3::erep_all1der(int &psh1, int &psh2, int &psh3, int &psh4,
685                       der_centersv3_t *der_centers)
686 {
687   double *current_buffer;
688   int nints;
689   double *user_int_buffer;
690   int omit;
691   GaussianBasisSet *cs[4];
692   int sh[4];
693   int n_unique;
694   int i,j;
695   GaussianShell *shell1,*shell2,*shell3,*shell4;
696   GaussianBasisSet *ucs[4]; /* The centers struct for the unique centers. */
697   int unum[4];        /* The number of times that this unique center occurs. */
698   int uam[4];         /* The total angular momentum on each unique center. */
699   int am[4];
700   int osh[4];
701   int cen[4];
702   int ucen[4];
703   int ncart;
704   double *current_pure_buffer;
705 
706   cs[0] = bs1_.pointer();
707   cs[1] = bs2_.pointer();
708   cs[2] = bs3_.pointer();
709   cs[3] = bs4_.pointer();
710 
711   sh[0] = psh1;
712   sh[1] = psh2;
713   sh[2] = psh3;
714   sh[3] = psh4;
715 
716   /* Set up pointers to the current shells. */
717   shell1 = &bs1_->shell(psh1);
718   shell2 = &bs2_->shell(psh2);
719   shell3 = &bs3_->shell(psh3);
720   shell4 = &bs4_->shell(psh4);
721 
722   /* Number of cartesian and pure integrals. */
723   ncart = shell1->ncartesian()*shell2->ncartesian()
724          *shell3->ncartesian()*shell4->ncartesian();
725   nints = shell1->nfunction()*shell2->nfunction()
726          *shell3->nfunction()*shell4->nfunction();
727 
728   am[0] = shell1->max_am();
729   am[1] = shell2->max_am();
730   am[2] = shell3->max_am();
731   am[3] = shell4->max_am();
732 
733   /* Compute the offset shell numbers. */
734   osh[0] = psh1 + bs1_shell_offset_;
735   osh[1] = psh2 + bs2_shell_offset_;
736   osh[2] = psh3 + bs3_shell_offset_;
737   osh[3] = psh4 + bs4_shell_offset_;
738 
739   for (i=0; i<4; i++) cen[i] = cs[i]->shell_to_center(sh[i]);
740 
741   /* This macro returns true if two shell centers are the same. */
742 #define SC(cs1,shc1,cs2,shc2) (((cs1)==(cs2))&&((shc1)==(shc2)))
743 
744   /* Build the list of unique centers structures and shells. */
745   n_unique = 0;
746   for (i=0; i<4; i++) {
747     int unique = 1;
748     for (j=0; j<n_unique; j++) {
749       if (SC(ucs[j],ucen[j],cs[i],cen[i])) {
750         unique = 0;
751         uam[j] += am[i];
752         unum[j]++;
753         break;
754         }
755       }
756     if (unique) {
757       ucs[n_unique] = cs[i];
758       ucen[n_unique] = cen[i];
759       uam[n_unique] = am[i];
760       unum[n_unique] = 1;
761       n_unique++;
762       }
763     }
764 
765   /* Choose which unique center is to be omitted from the calculation. */
766   if (n_unique == 1) {
767     omit = 0;
768     }
769   else if (n_unique == 2) {
770     if (unum[0]==3) omit = 0;
771     else if (unum[1]==3) omit = 1;
772     else if (uam[1]>uam[0]) omit = 1;
773     else omit = 0;
774     }
775   else if (n_unique == 3) {
776     if (unum[0]==2) omit = 0;
777     else if (unum[1]==2) omit = 1;
778     else omit = 2;
779     }
780   else {
781     int max = 0;
782     omit = 0;
783     for (i=0; i<4; i++) {
784       if (uam[i]>max) {
785         max = uam[i];
786         omit = i;
787         }
788       }
789     }
790 
791   /* Save the location of the int_buffer. */
792   user_int_buffer = int_buffer;
793   int_buffer = int_derint_buffer;
794 
795   /* Zero out the result integrals. */
796   for (i=0; i<3*(n_unique-1)*ncart; i++) user_int_buffer[i] = 0.0;
797 
798   /* Loop thru the unique centers, computing the integrals and
799    * skip the derivative on the unique center specified by omit. */
800   der_centers->n = 0;
801   current_buffer = user_int_buffer;
802   for (i=0; i<n_unique; i++) {
803     if (i==omit) continue;
804 
805     der_centers->cs[der_centers->n] = ucs[i];
806     der_centers->num[der_centers->n] = ucen[i];
807     der_centers->n++;
808 
809     for (j=0; j<4; j++) {
810       if (SC(ucs[i],ucen[i],cs[j],cen[j])) {
811         int old_perm = permute();
812         set_permute(0);
813         compute_erep_1der(0,current_buffer,
814                           &psh1,&psh2,&psh3,&psh4,j);
815         set_permute(old_perm);
816         }
817       }
818 
819     current_buffer = &current_buffer[3*ncart];
820     }
821 
822   /* Put the information about the omitted center into der_centers. */
823   der_centers->ocs = ucs[omit];
824   der_centers->onum = ucen[omit];
825 
826   /* Transform to pure am. */
827   current_buffer = user_int_buffer;
828   current_pure_buffer = user_int_buffer;
829   for (i=0; i<3*der_centers->n; i++) {
830       transform_2e(integral_, current_buffer, current_pure_buffer,
831                    shell1, shell2, shell3, shell4);
832       current_buffer = &current_buffer[ncart];
833       current_pure_buffer = &current_pure_buffer[nints];
834     }
835 
836   /* Eliminate redundant integrals, unless flags specifies otherwise. */
837   current_buffer = user_int_buffer;
838   if (!redundant_) {
839     int redundant_offset = 0;
840     int nonredundant_offset = 0;
841     int e12,e13e24,e34;
842     int i;
843 
844     if ((osh[0] == osh[3])&&(osh[1] == osh[2])&&(osh[0] != osh[1])) {
845       ExEnv::errn() << scprintf("nonredundant integrals cannot be generated (1der)\n");
846       fail();
847       }
848 
849     /* Shell equivalence information. */
850     e12 = (osh[0] == osh[1]);
851     e13e24 = ((osh[0] == osh[2]) && (osh[1] == osh[3]));
852     e34 = (osh[2] == osh[3]);
853     if (e12||e13e24||e34) {
854       /* Repack x, y, and z integrals. */
855       for (i=0; i<3*der_centers->n; i++) {
856         nonredundant_erep(current_buffer,e12,e34,e13e24,
857                           shell1->nfunction(),
858                           shell2->nfunction(),
859                           shell3->nfunction(),
860                           shell4->nfunction(),
861                           &redundant_offset,
862                           &nonredundant_offset);
863         }
864       }
865     }
866 
867   /* Return the integral buffers to their normal state. */
868   int_derint_buffer = int_buffer;
869   int_buffer = user_int_buffer;
870   }
871 
872 /* This will call compute_erep
873  * to compute the derivatives in terms of order == 0 integrals.
874  * flags are the flags to the int_comp_erep routine
875  * psh1-4 are pointers to the shell numbers
876  * dercenter is 0, 1, 2, or 3 -- the center within the integral
877  *           which we are taking the derivative with respect to.
878  * The results are accumulated in buffer, which cannot be the same
879  * as the current int_buffer.
880  */
881 void
compute_erep_1der(int flags,double * buffer,int * psh1,int * psh2,int * psh3,int * psh4,int dercenter)882 Int2eV3::compute_erep_1der(int flags, double *buffer,
883                            int *psh1, int *psh2, int *psh3, int *psh4,
884                            int dercenter)
885 {
886   int ii;
887   int index;
888   int size1,size2,size3,size4,size1234;
889   int sizem234,sizem34,sizem2,sizem3,sizem4;
890   int sizep234,sizep34,sizep2,sizep3,sizep4;
891   GaussianShell *shell1,*shell2,*shell3,*shell4;
892 
893 #ifdef DER_TIMING
894   tim_enter("erep_1der");
895 #endif
896 
897   /* Set up pointers to the current shells. */
898   shell1 = &bs1_->shell(*psh1);
899   shell2 = &bs2_->shell(*psh2);
900   shell3 = &bs3_->shell(*psh3);
901   shell4 = &bs4_->shell(*psh4);
902 
903   if ((dercenter<0) || (dercenter > 3)) {
904     ExEnv::errn() << scprintf("illegal derivative center -- must be 0, 1, 2, or 3\n");
905     fail();
906     }
907 
908   /* Offsets for the intermediates with original angular momentum. */
909   for (ii=size1=0; ii<shell1->ncontraction(); ii++)
910     size1 += INT_NCART(shell1->am(ii));
911   for (ii=size2=0; ii<shell2->ncontraction(); ii++)
912     size2 += INT_NCART(shell2->am(ii));
913   for (ii=size3=0; ii<shell3->ncontraction(); ii++)
914     size3 += INT_NCART(shell3->am(ii));
915   for (ii=size4=0; ii<shell4->ncontraction(); ii++)
916     size4 += INT_NCART(shell4->am(ii));
917 
918   size1234 = size1*size2*size3*size4;
919 
920 #define DCTEST(n) ((dercenter==n)?1:0)
921   /* Offsets for the intermediates with angular momentum decremented. */
922   for (ii=sizem2=0; ii<shell2->ncontraction(); ii++)
923     sizem2 += INT_NCART(shell2->am(ii)-DCTEST(1));
924   for (ii=sizem3=0; ii<shell3->ncontraction(); ii++)
925     sizem3 += INT_NCART(shell3->am(ii)-DCTEST(2));
926   for (ii=sizem4=0; ii<shell4->ncontraction(); ii++)
927     sizem4 += INT_NCART(shell4->am(ii)-DCTEST(3));
928   sizem34 = sizem4 * sizem3;
929   sizem234 = sizem34 * sizem2;
930 
931   /* Offsets for the intermediates with angular momentum incremented. */
932   for (ii=sizep2=0; ii<shell2->ncontraction(); ii++)
933     sizep2 += INT_NCART(shell2->am(ii)+DCTEST(1));
934   for (ii=sizep3=0; ii<shell3->ncontraction(); ii++)
935     sizep3 += INT_NCART(shell3->am(ii)+DCTEST(2));
936   for (ii=sizep4=0; ii<shell4->ncontraction(); ii++)
937     sizep4 += INT_NCART(shell4->am(ii)+DCTEST(3));
938   sizep34 = sizep4 * sizep3;
939   sizep234 = sizep34 * sizep2;
940 
941 #ifdef DER_TIMING
942   tim_enter("- erep");
943 #endif
944 
945   int old_perm = permute();
946   set_permute(0);
947   int old_red = redundant();
948   set_redundant(1);
949   compute_erep(flags|INT_NOPURE,
950                psh1,psh2,psh3,psh4,
951                    -DCTEST(0),
952                    -DCTEST(1),
953                    -DCTEST(2),
954                    -DCTEST(3));
955   set_permute(old_perm);
956   set_redundant(old_red);
957 
958   /* Trouble if cpp is nonANSI. */
959 #define DERLOOP(index,indexp1) \
960    oc##indexp1 = 0;\
961    for ( c##indexp1 =0; c##indexp1 <shell##indexp1->ncontraction(); c##indexp1 ++) {\
962      am[index] = shell##indexp1->am(c##indexp1);\
963      FOR_CART(i[index],j[index],k[index],am[index])
964 
965 #define END_DERLOOP(index,indexp1,sign) \
966        END_FOR_CART\
967      oc##indexp1 += INT_NCART(am[index] sign DCTEST(index));\
968      }
969 
970 #define ALLDERLOOPS\
971     DERLOOP(0,1)\
972       DERLOOP(1,2)\
973         DERLOOP(2,3)\
974           DERLOOP(3,4)
975 
976 #define END_ALLDERLOOPS(sign)\
977             END_DERLOOP(3,4,sign)\
978           END_DERLOOP(2,3,sign)\
979         END_DERLOOP(1,2,sign)\
980       END_DERLOOP(0,1,sign)
981 
982   /* Place the contributions into the user integral buffer. */
983   index = 0;
984 
985   if (dercenter==0) {
986     int ogc1,ogc1m,gc1,i1,k1,f234,size234;
987     size234=size2*size3*size4;
988 
989 #ifdef DER_TIMING
990   tim_change("- 0");
991 #endif
992     /* The center 0 d/dx integrals */
993     ogc1 = 0;
994     ogc1m = 0;
995     for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
996       int am1 = shell1->am(gc1);
997       // only integrals with x^n n>0 on center 0 contribute
998       // so skip over n==0
999       index += (am1+1)*size234;
1000       int c1 = am1+1;
1001       for (i1=1; i1<=am1; i1++) {
1002         double factor = -i1;
1003         for (k1=0; k1<=am1-i1; k1++) {
1004           int c1xm1 = c1-am1-1;//=INT_CARTINDEX(am1-1,i1-1,j1)
1005           double *tmp_buffer=&buffer[index];
1006           double *tmp_int_buffer=&int_buffer[(ogc1m+c1xm1)*size234];
1007           for (f234=0; f234<size234; f234++) {
1008             tmp_buffer[f234] += factor * tmp_int_buffer[f234];
1009             }
1010           index+=size234;
1011           c1++;
1012           }
1013         }
1014       ogc1 += c1;
1015       ogc1m += INT_NCART(am1-1);
1016       }
1017 
1018     /* The center 0 d/dy integrals */
1019     ogc1 = 0;
1020     ogc1m = 0;
1021     for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
1022       int am1 = shell1->am(gc1);
1023       // only integrals with y^n n>0 on center 0 contribute
1024       // so skip over n==0 (by making i1+k1<am1)
1025       int c1 = 0;
1026       for (i1=0; i1<=am1; i1++) {
1027         for (k1=0; k1<=am1-i1-1; k1++) {
1028           double factor = -(am1-i1-k1);
1029           int c1ym1 = c1-i1;//=INT_CARTINDEX(am1-1,i1,j1-1)
1030           double *tmp_buffer=&buffer[index];
1031           double *tmp_int_buffer=&int_buffer[(ogc1m+c1ym1)*size234];
1032           for (f234=0; f234<size234; f234++) {
1033             tmp_buffer[f234] += factor * tmp_int_buffer[f234];
1034             }
1035           index+=size234;
1036           c1++;
1037           }
1038         // account for the y^n n==0 case by an extra increment
1039         c1++;
1040         index+=size234;
1041         }
1042       ogc1 += c1;
1043       ogc1m += INT_NCART(am1-1);
1044       }
1045 
1046     /* The center 0 d/dz integrals */
1047     ogc1 = 0;
1048     ogc1m = 0;
1049     for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
1050       int am1 = shell1->am(gc1);
1051       int c1 = 0;
1052       for (i1=0; i1<=am1; i1++) {
1053         // only integrals with z^n n>0 on center 0 contribute
1054         // so skip over n==0
1055         c1++;
1056         index+=size234;
1057         for (k1=1; k1<=am1-i1; k1++) {
1058           double factor = -k1;
1059           int c1zm1 = c1-i1-1;//=INT_CARTINDEX(am1-1,i1,j1)
1060           double *tmp_buffer=&buffer[index];
1061           double *tmp_int_buffer=&int_buffer[(ogc1m+c1zm1)*size234];
1062           for (f234=0; f234<size234; f234++) {
1063             tmp_buffer[f234] += factor * tmp_int_buffer[f234];
1064             }
1065           index+=size234;
1066           c1++;
1067           }
1068         }
1069       ogc1 += c1;
1070       ogc1m += INT_NCART(am1-1);
1071       }
1072     }
1073   else if (dercenter == 1) {
1074     int ogc2,ogc2m,gc2,i2,k2,f1,f34,size34,size234;
1075     size34 = size3*size4;
1076     size234 = size2*size3*size4;
1077 
1078 #ifdef DER_TIMING
1079   tim_change("- 1");
1080 #endif
1081     /* The center 1 d/dx integrals */
1082     ogc2 = 0;
1083     ogc2m = 0;
1084     for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
1085       int am2 = shell2->am(gc2);
1086       // only integrals with x^n n>0 on center 1 contribute
1087       // so skip over n==0
1088       int c2 = am2+1;
1089       for (i2=1; i2<=am2; i2++) {
1090         double factor = -i2;
1091         for (k2=0; k2<=am2-i2; k2++) {
1092           int c2xm1 = c2-am2-1;//=INT_CARTINDEX(am2-1,i2-1,j2)
1093           int buffer_index = (ogc2+c2)*size34;
1094           int int_buffer_index = (ogc2m+c2xm1)*size34;
1095           for (f1=0; f1<size1; f1++) {
1096             double *tmp_buffer=&buffer[buffer_index];
1097             double *tmp_int_buffer=&int_buffer[int_buffer_index];
1098             for (f34=0; f34<size34; f34++) {
1099                 tmp_buffer[f34] += factor * tmp_int_buffer[f34];
1100               }
1101             buffer_index += size234;
1102             int_buffer_index += sizem234;
1103             }
1104           c2++;
1105           }
1106         }
1107       ogc2 += c2;
1108       ogc2m += INT_NCART(am2-1);
1109       }
1110     index += size1234;
1111 
1112      /* The center 1 d/dy integrals */
1113     ogc2 = 0;
1114     ogc2m = 0;
1115     for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
1116       int am2 = shell2->am(gc2);
1117       // only integrals with y^n n>0 on center 1 contribute
1118       // so skip over n==0
1119       int c2 = 0;
1120       for (i2=0; i2<=am2; i2++) {
1121         for (k2=0; k2<=am2-i2-1; k2++) {
1122           double factor = -(am2-k2-i2);
1123           int c2ym1 = c2-i2;//=INT_CARTINDEX(am2-1,i2,j2-1)
1124           int buffer_index = size1234 + (ogc2+c2)*size34;
1125           int int_buffer_index = (ogc2m+c2ym1)*size34;
1126           for (f1=0; f1<size1; f1++) {
1127             double *tmp_buffer=&buffer[buffer_index];
1128             double *tmp_int_buffer=&int_buffer[int_buffer_index];
1129             for (f34=0; f34<size34; f34++) {
1130                 tmp_buffer[f34] += factor * tmp_int_buffer[f34];
1131               }
1132             buffer_index += size234;
1133             int_buffer_index += sizem234;
1134             }
1135           c2++;
1136           }
1137         // account for the y^n n==0 case by an extra increment
1138         c2++;
1139         }
1140       ogc2 += c2;
1141       ogc2m += INT_NCART(am2-1);
1142       }
1143     index += size1234;
1144 
1145      /* The center 1 d/dz integrals */
1146     ogc2 = 0;
1147     ogc2m = 0;
1148     for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
1149       int am2 = shell2->am(gc2);
1150       // only integrals with z^n n>0 on center 1 contribute
1151       // so skip over n==0
1152       int c2 = 0;
1153       for (i2=0; i2<=am2; i2++) {
1154         // account for the z^n n==0 case by an extra increment
1155         c2++;
1156         for (k2=1; k2<=am2-i2; k2++) {
1157           double factor = -k2;
1158           int c2zm1 = c2-i2-1;//=INT_CARTINDEX(am2-1,i2,j2-1)
1159           int buffer_index = size1234+size1234+(ogc2+c2)*size34;
1160           int int_buffer_index = (ogc2m+c2zm1)*size34;
1161           for (f1=0; f1<size1; f1++) {
1162             double *tmp_buffer=&buffer[buffer_index];
1163             double *tmp_int_buffer=&int_buffer[int_buffer_index];
1164             for (f34=0; f34<size34; f34++) {
1165                 tmp_buffer[f34] += factor * tmp_int_buffer[f34];
1166               }
1167             buffer_index += size234;
1168             int_buffer_index += sizem234;
1169             }
1170           c2++;
1171           }
1172         }
1173       ogc2 += c2;
1174       ogc2m += INT_NCART(am2-1);
1175       }
1176     index += size1234;
1177     }
1178   else if (dercenter == 2) {
1179     int ogc3,ogc3m,gc3,i3,k3,f12,f4,size12,size34;
1180     size12 = size1*size2;
1181     size34 = size3*size4;
1182 
1183 #ifdef DER_TIMING
1184   tim_change("- 2");
1185 #endif
1186     /* The center 2 d/dx integrals */
1187     ogc3 = 0;
1188     ogc3m = 0;
1189     for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
1190       int am3 = shell3->am(gc3);
1191       // only integrals with x^n n>0 on center 2 contribute
1192       // so skip over n==0
1193       int c3 = am3+1;
1194       for (i3=1; i3<=am3; i3++) {
1195         double factor = -i3;
1196         for (k3=0; k3<=am3-i3; k3++) {
1197           int c3xm1 = c3-am3-1;//=INT_CARTINDEX(am3-1,i3-1,j3)
1198           int buffer_index = (ogc3+c3)*size4;
1199           int int_buffer_index = (ogc3m+c3xm1)*size4;
1200           for (f12=0; f12<size12; f12++) {
1201             double *tmp_buffer=&buffer[buffer_index];
1202             double *tmp_int_buffer=&int_buffer[int_buffer_index];
1203             for (f4=0; f4<size4; f4++) {
1204                 tmp_buffer[f4] += factor * tmp_int_buffer[f4];
1205               }
1206             buffer_index += size34;
1207             int_buffer_index += sizem34;
1208             }
1209           c3++;
1210           }
1211         }
1212       ogc3 += c3;
1213       ogc3m += INT_NCART(am3-1);
1214       }
1215     index += size1234;
1216 
1217      /* The center 2 d/dy integrals */
1218     ogc3 = 0;
1219     ogc3m = 0;
1220     for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
1221       int am3 = shell3->am(gc3);
1222       // only integrals with y^n n>0 on center 2 contribute
1223       // so skip over n==0
1224       int c3 = 0;
1225       for (i3=0; i3<=am3; i3++) {
1226         for (k3=0; k3<=am3-i3-1; k3++) {
1227           double factor = -(am3-k3-i3);
1228           int c3ym1 = c3-i3;//=INT_CARTINDEX(am3-1,i3,j3-1)
1229           int buffer_index = size1234 + (ogc3+c3)*size4;
1230           int int_buffer_index = (ogc3m+c3ym1)*size4;
1231           for (f12=0; f12<size12; f12++) {
1232             double *tmp_buffer=&buffer[buffer_index];
1233             double *tmp_int_buffer=&int_buffer[int_buffer_index];
1234             for (f4=0; f4<size4; f4++) {
1235                 tmp_buffer[f4] += factor * tmp_int_buffer[f4];
1236               }
1237             buffer_index += size34;
1238             int_buffer_index += sizem34;
1239             }
1240           c3++;
1241           }
1242         // account for the y^n n==0 case by an extra increment
1243         c3++;
1244         }
1245       ogc3 += c3;
1246       ogc3m += INT_NCART(am3-1);
1247       }
1248     index += size1234;
1249 
1250      /* The center 2 d/dz integrals */
1251     ogc3 = 0;
1252     ogc3m = 0;
1253     for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
1254       int am3 = shell3->am(gc3);
1255       // only integrals with z^n n>0 on center 2 contribute
1256       // so skip over n==0
1257       int c3 = 0;
1258       for (i3=0; i3<=am3; i3++) {
1259         // account for the z^n n==0 case by an extra increment
1260         c3++;
1261         for (k3=1; k3<=am3-i3; k3++) {
1262           double factor = -k3;
1263           int c3zm1 = c3-i3-1;//=INT_CARTINDEX(am3-1,i3,j3)
1264           int buffer_index = size1234+size1234+(ogc3+c3)*size4;
1265           int int_buffer_index = (ogc3m+c3zm1)*size4;
1266           for (f12=0; f12<size12; f12++) {
1267             double *tmp_buffer=&buffer[buffer_index];
1268             double *tmp_int_buffer=&int_buffer[int_buffer_index];
1269             for (f4=0; f4<size4; f4++) {
1270                 tmp_buffer[f4] += factor * tmp_int_buffer[f4];
1271               }
1272             buffer_index += size34;
1273             int_buffer_index += sizem34;
1274             }
1275           c3++;
1276           }
1277         }
1278       ogc3 += c3;
1279       ogc3m += INT_NCART(am3-1);
1280       }
1281     index += size1234;
1282     }
1283   else if (dercenter == 3) {
1284     int ogc4,ogc4m,gc4,i4,k4,f123,size123;
1285     size123 = size1*size2*size3;
1286 
1287 #ifdef DER_TIMING
1288   tim_change("- 3");
1289 #endif
1290     /* The center 3 d/dx integrals */
1291     ogc4 = 0;
1292     ogc4m = 0;
1293     for (gc4=0; gc4<shell4->ncontraction(); gc4++) {
1294       int am4 = shell4->am(gc4);
1295       // only integrals with x^n n>0 on center 3 contribute
1296       // so skip over n==0
1297       int c4 = am4+1;
1298       for (i4=1; i4<=am4; i4++) {
1299         double factor = -i4;
1300         for (k4=0; k4<=am4-i4; k4++) {
1301           int c4xm1 = c4-am4-1;//=INT_CARTINDEX(am4-1,i4-1,j4)
1302           int buffer_index = ogc4+c4;
1303           int int_buffer_index = ogc4m+c4xm1;
1304           for (f123=0; f123<size123; f123++) {
1305             buffer[buffer_index] += factor * int_buffer[int_buffer_index];
1306             buffer_index += size4;
1307             int_buffer_index += sizem4;
1308             }
1309           c4++;
1310           }
1311         }
1312       ogc4 += c4;
1313       ogc4m += INT_NCART(am4-1);
1314       }
1315     index += size1234;
1316 
1317      /* The center 3 d/dy integrals */
1318     ogc4 = 0;
1319     ogc4m = 0;
1320     for (gc4=0; gc4<shell4->ncontraction(); gc4++) {
1321       int am4 = shell4->am(gc4);
1322       // only integrals with y^n n>0 on center 3 contribute
1323       // so skip over n==0
1324       int c4 = 0;
1325       for (i4=0; i4<=am4; i4++) {
1326         for (k4=0; k4<=am4-i4-1; k4++) {
1327           double factor = -(am4-k4-i4);
1328           int c4ym1 = c4-i4;//=INT_CARTINDEX(am4-1,i4,j4-1)
1329           int buffer_index = size1234 + ogc4+c4;
1330           int int_buffer_index = ogc4m+c4ym1;
1331           for (f123=0; f123<size123; f123++) {
1332             buffer[buffer_index] += factor * int_buffer[int_buffer_index];
1333             buffer_index += size4;
1334             int_buffer_index += sizem4;
1335             }
1336           c4++;
1337           }
1338         // account for the y^n n==0 case by an extra increment
1339         c4++;
1340         }
1341       ogc4 += c4;
1342       ogc4m += INT_NCART(am4-1);
1343       }
1344     index += size1234;
1345 
1346     /* The center 3 d/dz integrals */
1347     ogc4 = 0;
1348     ogc4m = 0;
1349     for (gc4=0; gc4<shell4->ncontraction(); gc4++) {
1350       int am4 = shell4->am(gc4);
1351       // only integrals with z^n n>0 on center 3 contribute
1352       // so skip over n==0
1353       int c4 = 0;
1354       for (i4=0; i4<=am4; i4++) {
1355         // account for the z^n n==0 case by an extra increment
1356         c4++;
1357         for (k4=1; k4<=am4-i4; k4++) {
1358           double factor = -k4;
1359           int c4zm1 = c4-i4-1;//=INT_CARTINDEX(am4-1,i4,j4-1)
1360           int buffer_index = size1234+size1234+ogc4+c4;
1361           int int_buffer_index = ogc4m+c4zm1;
1362           for (f123=0; f123<size123; f123++) {
1363             buffer[buffer_index] += factor * int_buffer[int_buffer_index];
1364             buffer_index += size4;
1365             int_buffer_index += sizem4;
1366             }
1367           c4++;
1368           }
1369         }
1370       ogc4 += c4;
1371       ogc4m += INT_NCART(am4-1);
1372       }
1373     index += size1234;
1374     }
1375 
1376 #ifdef DER_TIMING
1377   tim_change("+ erep");
1378 #endif
1379 
1380   /* Compute the next contribution to the integrals. */
1381   /* Tell the build routine that we need an exponent weighted contraction
1382    * with the exponents taken from the dercenter and adjust the
1383    * angular momentum of dercenter to the needed value. */
1384   if (dercenter==0) int_expweight1 = 1;
1385   else if (dercenter==1) int_expweight2 = 1;
1386   else if (dercenter==2) int_expweight3 = 1;
1387   else if (dercenter==3) int_expweight4 = 1;
1388   old_perm = permute();
1389   set_permute(0);
1390   old_red = redundant();
1391   set_redundant(1);
1392   compute_erep(flags|INT_NOPURE,
1393                psh1,psh2,psh3,psh4,
1394                      DCTEST(0),
1395                      DCTEST(1),
1396                      DCTEST(2),
1397                      DCTEST(3));
1398   set_permute(old_perm);
1399   set_redundant(old_red);
1400   if (dercenter==0) int_expweight1 = 0;
1401   else if (dercenter==1) int_expweight2 = 0;
1402   else if (dercenter==2) int_expweight3 = 0;
1403   else if (dercenter==3) int_expweight4 = 0;
1404 
1405   /* Place the contributions into the user integral buffer. */
1406   index = 0;
1407   if (dercenter==0) {
1408     int ogc1,ogc1p,gc1,i1,k1,f234,size234;
1409     size234=size2*size3*size4;
1410 
1411 #ifdef DER_TIMING
1412   tim_change("+ 0");
1413 #endif
1414     /* The center 0 d/dx integrals */
1415     ogc1 = 0;
1416     ogc1p = 0;
1417     for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
1418       int am1 = shell1->am(gc1);
1419       int c1 = 0;
1420       for (i1=0; i1<=am1; i1++) {
1421         for (k1=0; k1<=am1-i1; k1++) {
1422           int c1xp1 = c1+am1+2;//=INT_CARTINDEX(am1+1,i1+1,j1)
1423           double *tmp_buffer=&buffer[index];
1424           double *tmp_int_buffer=&int_buffer[(ogc1p+c1xp1)*size234];
1425           for (f234=0; f234<size234; f234++) {
1426             tmp_buffer[f234] += tmp_int_buffer[f234];
1427             }
1428           index+=size234;
1429           c1++;
1430           }
1431         }
1432       ogc1 += c1;
1433       ogc1p += INT_NCART(am1+1);
1434       }
1435 
1436     /* The center 0 d/dy integrals */
1437     ogc1 = 0;
1438     ogc1p = 0;
1439     for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
1440       int am1 = shell1->am(gc1);
1441       int c1 = 0;
1442       for (i1=0; i1<=am1; i1++) {
1443         for (k1=0; k1<=am1-i1; k1++) {
1444           int c1yp1 = c1+i1;//=INT_CARTINDEX(am1+1,i1,j1+1)
1445           double *tmp_buffer=&buffer[index];
1446           double *tmp_int_buffer=&int_buffer[(ogc1p+c1yp1)*size234];
1447           for (f234=0; f234<size234; f234++) {
1448             tmp_buffer[f234] += tmp_int_buffer[f234];
1449             }
1450           index+=size234;
1451           c1++;
1452           }
1453         }
1454       ogc1 += c1;
1455       ogc1p += INT_NCART(am1+1);
1456       }
1457 
1458     /* The center 0 d/dz integrals */
1459     ogc1 = 0;
1460     ogc1p = 0;
1461     for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
1462       int am1 = shell1->am(gc1);
1463       int c1 = 0;
1464       for (i1=0; i1<=am1; i1++) {
1465         for (k1=0; k1<=am1-i1; k1++) {
1466           int c1zp1 = c1+i1+1;//=INT_CARTINDEX(am1+1,i1,j1)
1467           double *tmp_buffer=&buffer[index];
1468           double *tmp_int_buffer=&int_buffer[(ogc1p+c1zp1)*size234];
1469           for (f234=0; f234<size234; f234++) {
1470             tmp_buffer[f234] += tmp_int_buffer[f234];
1471             }
1472           index+=size234;
1473           c1++;
1474           }
1475         }
1476       ogc1 += c1;
1477       ogc1p += INT_NCART(am1+1);
1478       }
1479     }
1480   else if (dercenter == 1) {
1481     int ogc2,ogc2p,gc2,i2,k2,f1,f34,size34,size234;
1482     size34 = size3*size4;
1483     size234 = size2*size3*size4;
1484 
1485 #ifdef DER_TIMING
1486   tim_change("+ 1");
1487 #endif
1488     /* The center 1 d/dx integrals */
1489     ogc2 = 0;
1490     ogc2p = 0;
1491     for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
1492       int am2 = shell2->am(gc2);
1493       int c2=0;
1494       for (i2=0; i2<=am2; i2++) {
1495         for (k2=0; k2<=am2-i2; k2++) {
1496           int c2xp1 = c2+am2+2;//=INT_CARTINDEX(am2+1,i2+1,j2)
1497           int buffer_index = (ogc2+c2)*size34;
1498           int int_buffer_index = (ogc2p+c2xp1)*size34;
1499           for (f1=0; f1<size1; f1++) {
1500             double *tmp_buffer=&buffer[buffer_index];
1501             double *tmp_int_buffer=&int_buffer[int_buffer_index];
1502             for (f34=0; f34<size34; f34++) {
1503                 tmp_buffer[f34] += tmp_int_buffer[f34];
1504               }
1505             buffer_index += size234;
1506             int_buffer_index += sizep234;
1507             }
1508           c2++;
1509           }
1510         }
1511       ogc2 += c2;
1512       ogc2p += INT_NCART(am2+1);
1513       }
1514     index += size1234;
1515 
1516      /* The center 1 d/dy integrals */
1517     ogc2 = 0;
1518     ogc2p = 0;
1519     for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
1520       int am2 = shell2->am(gc2);
1521       int c2 = 0;
1522       for (i2=0; i2<=am2; i2++) {
1523         for (k2=0; k2<=am2-i2; k2++) {
1524           int c2yp1 = c2+i2;//=INT_CARTINDEX(am2+1,i2,j2+1)
1525           int buffer_index = size1234 + (ogc2+c2)*size34;
1526           int int_buffer_index = (ogc2p+c2yp1)*size34;
1527           for (f1=0; f1<size1; f1++) {
1528             double *tmp_buffer=&buffer[buffer_index];
1529             double *tmp_int_buffer=&int_buffer[int_buffer_index];
1530             for (f34=0; f34<size34; f34++) {
1531                 tmp_buffer[f34] += tmp_int_buffer[f34];
1532               }
1533             buffer_index += size234;
1534             int_buffer_index += sizep234;
1535             }
1536           c2++;
1537           }
1538         }
1539       ogc2 += c2;
1540       ogc2p += INT_NCART(am2+1);
1541       }
1542     index += size1234;
1543 
1544     /* The center 1 d/dz integrals */
1545     ogc2 = 0;
1546     ogc2p = 0;
1547     for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
1548       int am2 = shell2->am(gc2);
1549       int c2 = 0;
1550       for (i2=0; i2<=am2; i2++) {
1551         for (k2=0; k2<=am2-i2; k2++) {
1552           int c2zp1 = c2+i2+1;//=INT_CARTINDEX(am2+1,i2,j2+1)
1553           int buffer_index = size1234+size1234+(ogc2+c2)*size34;
1554           int int_buffer_index = (ogc2p+c2zp1)*size34;
1555           for (f1=0; f1<size1; f1++) {
1556             double *tmp_buffer=&buffer[buffer_index];
1557             double *tmp_int_buffer=&int_buffer[int_buffer_index];
1558             for (f34=0; f34<size34; f34++) {
1559                 tmp_buffer[f34] += tmp_int_buffer[f34];
1560               }
1561             buffer_index += size234;
1562             int_buffer_index += sizep234;
1563             }
1564           c2++;
1565           }
1566         }
1567       ogc2 += c2;
1568       ogc2p += INT_NCART(am2+1);
1569       }
1570     index += size1234;
1571     }
1572   else if (dercenter == 2) {
1573     int ogc3,ogc3p,gc3,i3,k3,f12,f4,size12,size34;
1574     size12 = size1*size2;
1575     size34 = size3*size4;
1576 
1577 #ifdef DER_TIMING
1578   tim_change("+ 2");
1579 #endif
1580     /* The center 2 d/dx integrals */
1581     ogc3 = 0;
1582     ogc3p = 0;
1583     for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
1584       int am3 = shell3->am(gc3);
1585       int c3 = 0;
1586       for (i3=0; i3<=am3; i3++) {
1587         for (k3=0; k3<=am3-i3; k3++) {
1588           int c3xp1 = c3+am3+2;//=INT_CARTINDEX(am3+1,i3+1,j3)
1589           int buffer_index = (ogc3+c3)*size4;
1590           int int_buffer_index = (ogc3p+c3xp1)*size4;
1591           for (f12=0; f12<size12; f12++) {
1592             double *tmp_buffer=&buffer[buffer_index];
1593             double *tmp_int_buffer=&int_buffer[int_buffer_index];
1594             for (f4=0; f4<size4; f4++) {
1595                 tmp_buffer[f4] += tmp_int_buffer[f4];
1596               }
1597             buffer_index += size34;
1598             int_buffer_index += sizep34;
1599             }
1600           c3++;
1601           }
1602         }
1603       ogc3 += c3;
1604       ogc3p += INT_NCART(am3+1);
1605       }
1606     index += size1234;
1607 
1608      /* The center 2 d/dy integrals */
1609     ogc3 = 0;
1610     ogc3p = 0;
1611     for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
1612       int am3 = shell3->am(gc3);
1613       int c3 = 0;
1614       for (i3=0; i3<=am3; i3++) {
1615         for (k3=0; k3<=am3-i3; k3++) {
1616           int c3yp1 = c3+i3;//=INT_CARTINDEX(am3+1,i3,j3+1)
1617           int buffer_index = size1234 + (ogc3+c3)*size4;
1618           int int_buffer_index = (ogc3p+c3yp1)*size4;
1619           for (f12=0; f12<size12; f12++) {
1620             double *tmp_buffer=&buffer[buffer_index];
1621             double *tmp_int_buffer=&int_buffer[int_buffer_index];
1622             for (f4=0; f4<size4; f4++) {
1623                 tmp_buffer[f4] += tmp_int_buffer[f4];
1624               }
1625             buffer_index += size34;
1626             int_buffer_index += sizep34;
1627             }
1628           c3++;
1629           }
1630         }
1631       ogc3 += c3;
1632       ogc3p += INT_NCART(am3+1);
1633       }
1634     index += size1234;
1635 
1636      /* The center 2 d/dz integrals */
1637     ogc3 = 0;
1638     ogc3p = 0;
1639     for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
1640       int am3 = shell3->am(gc3);
1641       int c3 = 0;
1642       for (i3=0; i3<=am3; i3++) {
1643         for (k3=0; k3<=am3-i3; k3++) {
1644           int c3zp1 = c3+i3+1;//=INT_CARTINDEX(am3+1,i3,j3)
1645           int buffer_index = size1234+size1234+(ogc3+c3)*size4;
1646           int int_buffer_index = (ogc3p+c3zp1)*size4;
1647           for (f12=0; f12<size12; f12++) {
1648             double *tmp_buffer=&buffer[buffer_index];
1649             double *tmp_int_buffer=&int_buffer[int_buffer_index];
1650             for (f4=0; f4<size4; f4++) {
1651                 tmp_buffer[f4] += tmp_int_buffer[f4];
1652               }
1653             buffer_index += size34;
1654             int_buffer_index += sizep34;
1655             }
1656           c3++;
1657           }
1658         }
1659       ogc3 += c3;
1660       ogc3p += INT_NCART(am3+1);
1661       }
1662     index += size1234;
1663     }
1664   else if (dercenter == 3) {
1665     int ogc4,ogc4p,gc4,i4,k4,f123,size123;
1666     size123 = size1*size2*size3;
1667 
1668 #ifdef DER_TIMING
1669   tim_change("+ 3");
1670 #endif
1671     /* The center 3 d/dx integrals */
1672     ogc4 = 0;
1673     ogc4p = 0;
1674     for (gc4=0; gc4<shell4->ncontraction(); gc4++) {
1675       int am4 = shell4->am(gc4);
1676       int c4 = 0;
1677       for (i4=0; i4<=am4; i4++) {
1678         for (k4=0; k4<=am4-i4; k4++) {
1679           int c4xp1 = c4+am4+2;//=INT_CARTINDEX(am4+1,i4+1,j4)
1680           int buffer_index = ogc4+c4;
1681           int int_buffer_index = ogc4p+c4xp1;
1682           for (f123=0; f123<size123; f123++) {
1683             buffer[buffer_index] += int_buffer[int_buffer_index];
1684             buffer_index += size4;
1685             int_buffer_index += sizep4;
1686             }
1687           c4++;
1688           }
1689         }
1690       ogc4 += c4;
1691       ogc4p += INT_NCART(am4+1);
1692       }
1693     index += size1234;
1694 
1695      /* The center 3 d/dy integrals */
1696     ogc4 = 0;
1697     ogc4p = 0;
1698     for (gc4=0; gc4<shell4->ncontraction(); gc4++) {
1699       int am4 = shell4->am(gc4);
1700       int c4 = 0;
1701       for (i4=0; i4<=am4; i4++) {
1702         for (k4=0; k4<=am4-i4; k4++) {
1703           int c4yp1 = c4+i4;//=INT_CARTINDEX(am4+1,i4,j4+1)
1704           int buffer_index = size1234 + ogc4+c4;
1705           int int_buffer_index = ogc4p+c4yp1;
1706           for (f123=0; f123<size123; f123++) {
1707             buffer[buffer_index] += int_buffer[int_buffer_index];
1708             buffer_index += size4;
1709             int_buffer_index += sizep4;
1710             }
1711           c4++;
1712           }
1713         }
1714       ogc4 += c4;
1715       ogc4p += INT_NCART(am4+1);
1716       }
1717     index += size1234;
1718 
1719     /* The center 3 d/dz integrals */
1720     ogc4 = 0;
1721     ogc4p = 0;
1722     for (gc4=0; gc4<shell4->ncontraction(); gc4++) {
1723       int am4 = shell4->am(gc4);
1724       int c4 = 0;
1725       for (i4=0; i4<=am4; i4++) {
1726         for (k4=0; k4<=am4-i4; k4++) {
1727           int c4zp1 = c4+i4+1;//=INT_CARTINDEX(am4+1,i4,j4)
1728           int buffer_index = size1234+size1234+ogc4+c4;
1729           int int_buffer_index = ogc4p+c4zp1;
1730           for (f123=0; f123<size123; f123++) {
1731             buffer[buffer_index] += int_buffer[int_buffer_index];
1732             buffer_index += size4;
1733             int_buffer_index += sizep4;
1734             }
1735           c4++;
1736           }
1737         }
1738       ogc4 += c4;
1739       ogc4p += INT_NCART(am4+1);
1740       }
1741     index += size1234;
1742     }
1743 #ifdef DER_TIMING
1744   tim_exit(0);
1745   tim_exit(0);
1746 #endif
1747   }
1748 
1749 void
nonredundant_erep(double * buffer,int e12,int e34,int e13e24,int n1,int n2,int n3,int n4,int * red_off,int * nonred_off)1750 Int2eV3::nonredundant_erep(double *buffer, int e12, int e34, int e13e24,
1751                            int n1, int n2, int n3, int n4,
1752                            int *red_off, int *nonred_off)
1753 {
1754   int nonredundant_index;
1755   int i,j,k,l;
1756 
1757   double *redundant_ptr = &buffer[*red_off];
1758   double *nonredundant_ptr = &buffer[*nonred_off];
1759 
1760   nonredundant_index = 0;
1761   int n34 = n3*n4;
1762   for (i=0; i<n1; i++) {
1763     int jmax = INT_MAX2(e12,i,n2);
1764     for (j=0; j<=jmax; j++) {
1765       int kmax = INT_MAX3(e13e24,i,n3);
1766       for (k=0; k<=kmax; k++) {
1767         int lmax = INT_MAX4(e13e24,e34,i,j,k,n4);
1768         for (l=0; l<=lmax; l++) {
1769           nonredundant_ptr[l] = redundant_ptr[l];
1770           }
1771         redundant_ptr += n4;
1772         nonredundant_index += lmax+1;
1773         nonredundant_ptr += lmax+1;
1774         }
1775       redundant_ptr += (n3-(kmax+1))*n4;
1776       }
1777     redundant_ptr += (n2-(jmax+1))*n34;
1778     }
1779   *red_off += n1*n2*n34;
1780   *nonred_off += nonredundant_index;
1781   }
1782 
1783 /* This is an alternate interface to int_erep_all1der.  It takes
1784  * as arguments the flags, an integer vector of shell numbers
1785  * and an integer vector which will be filled in with size
1786  * information, if it is non-NULL, and the dercenters pointer. */
1787 void
erep_all1der(int * shells,int * sizes,der_centersv3_t * dercenters)1788 Int2eV3::erep_all1der(int *shells, int  *sizes,
1789                       der_centersv3_t *dercenters)
1790 {
1791   erep_all1der(shells[0],shells[1],shells[2],shells[3],
1792                dercenters);
1793   if (sizes) {
1794     sizes[0] = bs1_->shell(shells[0]).nfunction();
1795     sizes[1] = bs2_->shell(shells[1]).nfunction();
1796     sizes[2] = bs3_->shell(shells[2]).nfunction();
1797     sizes[3] = bs4_->shell(shells[3]).nfunction();
1798     }
1799   }
1800 
1801 void
int_erep_bound1der(int flags,int bsh1,int bsh2,int * size)1802 Int2eV3::int_erep_bound1der(int flags, int bsh1, int bsh2, int *size)
1803 {
1804   double *current_buffer;
1805   int nints;
1806   double *user_int_buffer;
1807   int i;
1808   GaussianShell *shell1,*shell2,*shell3,*shell4;
1809   int osh[4];
1810   int sh1 = bsh1;
1811   int sh2 = bsh2;
1812   int sh3 = bsh1;
1813   int sh4 = bsh2;
1814   int *psh1 = &sh1;
1815   int *psh2 = &sh2;
1816   int *psh3 = &sh3;
1817   int *psh4 = &sh4;
1818   int ncart;
1819   double *current_pure_buffer;
1820 
1821   /* Set up pointers to the current shells. */
1822   shell1 = &bs1_->shell(*psh1);
1823   shell2 = &bs2_->shell(*psh2);
1824   shell3 = &bs3_->shell(*psh3);
1825   shell4 = &bs4_->shell(*psh4);
1826 
1827   /* Number of cartesian and pure integrals. */
1828   ncart = shell1->ncartesian()*shell2->ncartesian()
1829          *shell3->ncartesian()*shell4->ncartesian();
1830   nints = shell1->nfunction() * shell2->nfunction()
1831         * shell3->nfunction() * shell4->nfunction();
1832 
1833   /* Compute the offset shell numbers. */
1834   osh[0] = *psh1 + bs1_shell_offset_;
1835   osh[1] = *psh2 + bs2_shell_offset_;
1836   osh[2] = *psh3 + bs3_shell_offset_;
1837   osh[3] = *psh4 + bs4_shell_offset_;
1838 
1839   /* Save the location of the int_buffer. */
1840   user_int_buffer = int_buffer;
1841   int_buffer = int_derint_buffer;
1842 
1843   /* Zero out the result integrals. */
1844   for (i=0; i<3*ncart; i++) user_int_buffer[i] = 0.0;
1845 
1846   /* Set the size so it is available to the caller. */
1847   *size = nints * 3;
1848 
1849   current_buffer = user_int_buffer;
1850   int old_perm = permute();
1851   compute_erep_bound1der(flags|INT_NOPURE,current_buffer,
1852                           psh1,psh2,psh3,psh4);
1853   set_permute(old_perm);
1854 
1855   /* Transform to pure am. */
1856   current_buffer = user_int_buffer;
1857   current_pure_buffer = user_int_buffer;
1858   for (i=0; i<3; i++) {
1859       transform_2e(integral_, current_buffer, current_pure_buffer,
1860                    &bs1_->shell(sh1),
1861                    &bs2_->shell(sh2),
1862                    &bs3_->shell(sh3),
1863                    &bs4_->shell(sh4));
1864       current_buffer = &current_buffer[ncart];
1865       current_pure_buffer = &current_pure_buffer[nints];
1866     }
1867 
1868   /* Eliminate redundant integrals, unless flags specifies otherwise. */
1869   current_buffer = user_int_buffer;
1870   if (!redundant_) {
1871     int redundant_offset = 0;
1872     int nonredundant_offset = 0;
1873     int e12,e13e24,e34;
1874     int i;
1875 
1876     if ((osh[0] == osh[3])&&(osh[1] == osh[2])&&(osh[0] != osh[1])) {
1877       ExEnv::errn() << scprintf("nonredundant integrals cannot be generated (1der)\n");
1878       fail();
1879       }
1880 
1881     /* Shell equivalence information. */
1882     e12 = (osh[0] == osh[1]);
1883     e13e24 = ((osh[0] == osh[2]) && (osh[1] == osh[3]));
1884     e34 = (osh[2] == osh[3]);
1885     /* Repack x, y, and z integrals. */
1886     for (i=0; i<3; i++) {
1887       nonredundant_erep(current_buffer,e12,e34,e13e24,
1888                              shell1->nfunction(),
1889                              shell2->nfunction(),
1890                              shell3->nfunction(),
1891                              shell4->nfunction(),
1892                              &redundant_offset,
1893                              &nonredundant_offset);
1894       }
1895     }
1896 
1897   /* Return the integral buffers to their normal state. */
1898   int_derint_buffer = int_buffer;
1899   int_buffer = user_int_buffer;
1900   }
1901 
1902 /* This routine computes a quantity needed to compute the
1903  * derivative integral bounds.
1904  * It fills int_buffer with (sh1+i sh2|sh1+i sh2).
1905  */
1906 void
compute_erep_bound1der(int flags,double * buffer,int * psh1,int * psh2,int * psh3,int * psh4)1907 Int2eV3::compute_erep_bound1der(int flags, double *buffer,
1908                                 int *psh1, int *psh2, int *psh3, int *psh4)
1909 {
1910   int oc1,oc2,oc3,oc4;
1911   int ii;
1912   int c1,c2,c3,c4;
1913   int i[4],j[4],k[4],am[4];
1914   int index;
1915   int sizem234,sizem34,sizem2,sizem3,sizem4;
1916   int sizep234,sizep34,sizep2,sizep3,sizep4;
1917   int sizepm234,sizepm34,sizepm2,sizepm3,sizepm4;
1918   GaussianShell *shell1,*shell2,*shell3,*shell4;
1919 
1920   /* Set up pointers to the current shells. */
1921   shell1 = &bs1_->shell(*psh1);
1922   shell2 = &bs2_->shell(*psh2);
1923   shell3 = &bs3_->shell(*psh3);
1924   shell4 = &bs4_->shell(*psh4);
1925 
1926 #define DCT1(n) ((((n)==0)||((n)==2))?1:0)
1927 #define DCT2(n) ((((n)==0)||((n)==2))?((((n)==0))?1:-1):0)
1928   /* Offsets for the intermediates with angular momentum decremented. */
1929   for (ii=sizem2=0; ii<shell2->ncontraction(); ii++)
1930     sizem2 += INT_NCART(shell2->am(ii)-DCT1(1));
1931   for (ii=sizem3=0; ii<shell3->ncontraction(); ii++)
1932     sizem3 += INT_NCART(shell3->am(ii)-DCT1(2));
1933   for (ii=sizem4=0; ii<shell4->ncontraction(); ii++)
1934     sizem4 += INT_NCART(shell4->am(ii)-DCT1(3));
1935   sizem34 = sizem4 * sizem3;
1936   sizem234 = sizem34 * sizem2;
1937 
1938   /* Offsets for the intermediates with angular momentum incremented. */
1939   for (ii=sizep2=0; ii<shell2->ncontraction(); ii++)
1940     sizep2 += INT_NCART(shell2->am(ii)+DCT1(1));
1941   for (ii=sizep3=0; ii<shell3->ncontraction(); ii++)
1942     sizep3 += INT_NCART(shell3->am(ii)+DCT1(2));
1943   for (ii=sizep4=0; ii<shell4->ncontraction(); ii++)
1944     sizep4 += INT_NCART(shell4->am(ii)+DCT1(3));
1945   sizep34 = sizep4 * sizep3;
1946   sizep234 = sizep34 * sizep2;
1947 
1948   /* Offsets for the intermediates with angular momentum incremented and
1949    * decremented. */
1950   for (ii=sizepm2=0; ii<shell2->ncontraction(); ii++)
1951     sizepm2 += INT_NCART(shell2->am(ii)+DCT2(1));
1952   for (ii=sizepm3=0; ii<shell3->ncontraction(); ii++)
1953     sizepm3 += INT_NCART(shell3->am(ii)+DCT2(2));
1954   for (ii=sizepm4=0; ii<shell4->ncontraction(); ii++)
1955     sizepm4 += INT_NCART(shell4->am(ii)+DCT2(3));
1956   sizepm34 = sizepm4 * sizepm3;
1957   sizepm234 = sizepm34 * sizepm2;
1958 
1959   /* END_DERLOOP must be redefined here because it previously depended
1960    * on the DCTEST macro */
1961 #undef END_DERLOOP
1962 #define END_DERLOOP(index,indexp1,sign) \
1963        END_FOR_CART\
1964      oc##indexp1 += INT_NCART(am[index] sign DCT1(index));\
1965      }
1966 
1967 #undef END_ALLDERLOOPS
1968 #define END_ALLDERLOOPS(sign)\
1969             END_DERLOOP(3,4,sign)\
1970           END_DERLOOP(2,3,sign)\
1971         END_DERLOOP(1,2,sign)\
1972       END_DERLOOP(0,1,sign)
1973 
1974   int old_perm = permute();
1975   set_permute(0);
1976   int old_red = redundant();
1977   set_redundant(1);
1978   compute_erep(flags,psh1,psh2,psh3,psh4,
1979                    -DCT1(0),
1980                    -DCT1(1),
1981                    -DCT1(2),
1982                    -DCT1(3));
1983   set_permute(old_perm);
1984   set_redundant(old_red);
1985 
1986    /* Place the contributions into the user integral buffer. */
1987    index = 0;
1988    /* The d/dx integrals */
1989   ALLDERLOOPS
1990     if (i[0]>0 && i[2]>0) {
1991       buffer[index] += i[0] * i[2] * int_buffer[
1992         (oc1 + INT_CARTINDEX(am[0]-DCT1(0),i[0]-DCT1(0),j[0])) * sizem234
1993        +(oc2 + INT_CARTINDEX(am[1]-DCT1(1),i[1]-DCT1(1),j[1])) * sizem34
1994        +(oc3 + INT_CARTINDEX(am[2]-DCT1(2),i[2]-DCT1(2),j[2])) * sizem4
1995        +(oc4 + INT_CARTINDEX(am[3]-DCT1(3),i[3]-DCT1(3),j[3]))
1996        ];
1997       }
1998     index++;
1999     END_ALLDERLOOPS(-)
2000 
2001    /* The d/dy integrals */
2002   ALLDERLOOPS
2003     if (j[0]>0 && j[2]>0) {
2004     buffer[index] += j[0] * j[2] * int_buffer[
2005          (oc1 + INT_CARTINDEX(am[0]-DCT1(0),i[0],j[0]-DCT1(0))) * sizem234
2006         +(oc2 + INT_CARTINDEX(am[1]-DCT1(1),i[1],j[1]-DCT1(1))) * sizem34
2007         +(oc3 + INT_CARTINDEX(am[2]-DCT1(2),i[2],j[2]-DCT1(2))) * sizem4
2008         +(oc4 + INT_CARTINDEX(am[3]-DCT1(3),i[3],j[3]-DCT1(3)))
2009         ];
2010       }
2011     index++;
2012   END_ALLDERLOOPS(-)
2013 
2014    /* The d/dz integrals */
2015   ALLDERLOOPS
2016     if (k[0]>0 && k[2]>0) {
2017     buffer[index] += k[0] * k[2] * int_buffer[
2018          (oc1 + INT_CARTINDEX(am[0]-DCT1(0),i[0],j[0])) * sizem234
2019         +(oc2 + INT_CARTINDEX(am[1]-DCT1(1),i[1],j[1])) * sizem34
2020         +(oc3 + INT_CARTINDEX(am[2]-DCT1(2),i[2],j[2])) * sizem4
2021         +(oc4 + INT_CARTINDEX(am[3]-DCT1(3),i[3],j[3]))
2022         ];
2023       }
2024     index++;
2025   END_ALLDERLOOPS(-)
2026 
2027   /* Compute the next contribution to the integrals. */
2028   /* Tell the build routine that we need an exponent weighted contraction
2029    * with the exponents taken from the dercenter and adjust the
2030    * angular momentum of dercenter to the needed value. */
2031   int_expweight1 = 1;
2032   int_expweight3 = 1;
2033   old_perm = permute();
2034   set_permute(0);
2035   old_red = redundant();
2036   set_redundant(1);
2037   compute_erep(flags,psh1,psh2,psh3,psh4,
2038                      DCT1(0),
2039                      DCT1(1),
2040                      DCT1(2),
2041                      DCT1(3));
2042   set_permute(old_perm);
2043   set_redundant(old_red);
2044   int_expweight1 = 0;
2045   int_expweight3 = 0;
2046 
2047   /* Place the contributions into the user integral buffer. */
2048   index = 0;
2049   /* The d/dx integrals */
2050   ALLDERLOOPS
2051           buffer[index] += int_buffer[
2052              (oc1+INT_CARTINDEX(am[0]+DCT1(0),i[0]+DCT1(0),j[0]))*sizep234
2053             +(oc2+INT_CARTINDEX(am[1]+DCT1(1),i[1]+DCT1(1),j[1]))*sizep34
2054             +(oc3+INT_CARTINDEX(am[2]+DCT1(2),i[2]+DCT1(2),j[2]))*sizep4
2055             +(oc4+INT_CARTINDEX(am[3]+DCT1(3),i[3]+DCT1(3),j[3]))
2056             ];
2057     index++;
2058     END_ALLDERLOOPS(+)
2059 
2060   /* The d/dy integrals */
2061   ALLDERLOOPS
2062           buffer[index] += int_buffer[
2063              (oc1+INT_CARTINDEX(am[0]+DCT1(0),i[0],j[0]+DCT1(0)))*sizep234
2064             +(oc2+INT_CARTINDEX(am[1]+DCT1(1),i[1],j[1]+DCT1(1)))*sizep34
2065             +(oc3+INT_CARTINDEX(am[2]+DCT1(2),i[2],j[2]+DCT1(2)))*sizep4
2066             +(oc4+INT_CARTINDEX(am[3]+DCT1(3),i[3],j[3]+DCT1(3)))
2067             ];
2068           index++;
2069     END_ALLDERLOOPS(+)
2070 
2071   /* The d/dz integrals */
2072   ALLDERLOOPS
2073           buffer[index] += int_buffer[
2074                (oc1 + INT_CARTINDEX(am[0]+DCT1(0),i[0],j[0])) * sizep234
2075               +(oc2 + INT_CARTINDEX(am[1]+DCT1(1),i[1],j[1])) * sizep34
2076               +(oc3 + INT_CARTINDEX(am[2]+DCT1(2),i[2],j[2])) * sizep4
2077               +(oc4 + INT_CARTINDEX(am[3]+DCT1(3),i[3],j[3]))
2078               ];
2079           index++;
2080     END_ALLDERLOOPS(+)
2081 
2082   /* END_DERLOOP must be redefined here because it previously depended
2083    * on the DCT1 macro */
2084 #undef END_DERLOOP
2085 #define END_DERLOOP(index,indexp1,sign) \
2086        END_FOR_CART\
2087      oc##indexp1 += INT_NCART(am[index] sign DCT2(index));\
2088      }
2089 
2090 #undef END_ALLDERLOOPS
2091 #define END_ALLDERLOOPS(sign)\
2092             END_DERLOOP(3,4,sign)\
2093           END_DERLOOP(2,3,sign)\
2094         END_DERLOOP(1,2,sign)\
2095       END_DERLOOP(0,1,sign)
2096 
2097   /* Compute the next contribution to the integrals. */
2098   /* Tell the build routine that we need an exponent weighted contraction
2099    * with the exponents taken from the dercenter and adjust the
2100    * angular momentum of dercenter to the needed value. */
2101   int_expweight1 = 1;
2102   old_perm = permute();
2103   set_permute(0);
2104   old_red = redundant();
2105   set_redundant(1);
2106   compute_erep(flags,psh1,psh2,psh3,psh4,
2107                      DCT2(0),
2108                      DCT2(1),
2109                      DCT2(2),
2110                      DCT2(3));
2111   set_permute(old_perm);
2112   set_redundant(old_red);
2113   int_expweight1 = 0;
2114 
2115   /* Place the contributions into the user integral buffer. */
2116   index = 0;
2117   /* The d/dx integrals */
2118   ALLDERLOOPS
2119     if (i[2] > 0)  {
2120           buffer[index] -= 2.0 * int_buffer[
2121              (oc1+INT_CARTINDEX(am[0]+DCT2(0),i[0]+DCT2(0),j[0]))*sizepm234
2122             +(oc2+INT_CARTINDEX(am[1]+DCT2(1),i[1]+DCT2(1),j[1]))*sizepm34
2123             +(oc3+INT_CARTINDEX(am[2]+DCT2(2),i[2]+DCT2(2),j[2]))*sizepm4
2124             +(oc4+INT_CARTINDEX(am[3]+DCT2(3),i[3]+DCT2(3),j[3]))
2125             ];
2126        }
2127     index++;
2128     END_ALLDERLOOPS(+)
2129 
2130   /* The d/dy integrals */
2131   ALLDERLOOPS
2132     if (j[2] > 0)  {
2133           buffer[index] -= 2.0 * int_buffer[
2134              (oc1+INT_CARTINDEX(am[0]+DCT2(0),i[0],j[0]+DCT2(0)))*sizepm234
2135             +(oc2+INT_CARTINDEX(am[1]+DCT2(1),i[1],j[1]+DCT2(1)))*sizepm34
2136             +(oc3+INT_CARTINDEX(am[2]+DCT2(2),i[2],j[2]+DCT2(2)))*sizepm4
2137             +(oc4+INT_CARTINDEX(am[3]+DCT2(3),i[3],j[3]+DCT2(3)))
2138             ];
2139        }
2140     index++;
2141     END_ALLDERLOOPS(+)
2142 
2143   /* The d/dz integrals */
2144   ALLDERLOOPS
2145     if (k[2] > 0)  {
2146           buffer[index] -= 2.0 * int_buffer[
2147                (oc1 + INT_CARTINDEX(am[0]+DCT2(0),i[0],j[0])) * sizepm234
2148               +(oc2 + INT_CARTINDEX(am[1]+DCT2(1),i[1],j[1])) * sizepm34
2149               +(oc3 + INT_CARTINDEX(am[2]+DCT2(2),i[2],j[2])) * sizepm4
2150               +(oc4 + INT_CARTINDEX(am[3]+DCT2(3),i[3],j[3]))
2151               ];
2152        }
2153     index++;
2154     END_ALLDERLOOPS(+)
2155   }
2156 
2157 /////////////////////////////////////////////////////////////////////////////
2158 
2159 // Local Variables:
2160 // mode: c++
2161 // c-file-style: "CLJ-CONDENSED"
2162 // End:
2163