1 //
2 // tformv3.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 #if defined(__GNUC__)
29 #pragma implementation
30 #endif
31 
32 #include <stdlib.h>
33 #include <string.h>
34 #include <math.h>
35 
36 #include <util/misc/formio.h>
37 #include <chemistry/qc/basis/integral.h>
38 #include <chemistry/qc/intv3/macros.h>
39 #include <chemistry/qc/intv3/tformv3.h>
40 #include <chemistry/qc/intv3/utils.h>
41 
42 using namespace sc;
43 
44 ////////////////////////////////////////////////////////////////////////////
45 
46 #define PRINT 0
47 
48 void
transform_init()49 Int2eV3::transform_init()
50 {
51   source = 0;
52   nsourcemax = 0;
53 }
54 
55 void
transform_done()56 Int2eV3::transform_done()
57 {
58   delete[] source;
59 }
60 
61 void
transform_init()62 Int1eV3::transform_init()
63 {
64   source = 0;
65   nsourcemax = 0;
66 }
67 
68 void
transform_done()69 Int1eV3::transform_done()
70 {
71   delete[] source;
72 }
73 
74 static void
do_copy1(double * source,double * target,int chunk,int n1,int s1,int offset1,int n2,int s2,int offset2)75 do_copy1(double *source, double *target, int chunk,
76          int n1, int s1, int offset1,
77          int n2, int s2, int offset2)
78 {
79   int i1, i2;
80 
81   for (i1=0; i1<n1; i1++) {
82       int off = ((offset1 + i1)*s2 + offset2)*chunk;
83       for (i2=0; i2<n2*chunk; i2++, off++) {
84           target[off] = source[off];
85         }
86     }
87 }
88 
89 static void
do_copy2(double * source,double * target,int n1,int s1,int offset1,int n2,int s2,int offset2,int n3,int s3,int offset3,int n4,int s4,int offset4)90 do_copy2(double *source, double *target,
91          int n1, int s1, int offset1,
92          int n2, int s2, int offset2,
93          int n3, int s3, int offset3,
94          int n4, int s4, int offset4)
95 {
96   int i1, i2, i3, i4;
97 
98   for (i1=0; i1<n1; i1++) {
99       for (i2=0; i2<n2; i2++) {
100           for (i3=0; i3<n3; i3++) {
101               int off = (((offset1 + i1)*s2 + offset2 + i2)
102                          *s3 + offset3 + i3)*s4 + offset4;
103               for (i4=0; i4<n4; i4++, off++) {
104                   target[off] = source[off];
105                 }
106             }
107         }
108     }
109 }
110 
111 static void
do_sparse_transform11(double * source,double * target,int chunk,SphericalTransformIter & trans,int offsetcart1,int offsetpure1,int n2,int s2,int offset2)112 do_sparse_transform11(double *source, double *target, int chunk,
113                       SphericalTransformIter& trans,
114                       int offsetcart1,
115                       int offsetpure1,
116                       int n2, int s2, int offset2)
117 {
118   int i2;
119 
120   for (trans.begin(); trans.ready(); trans.next()) {
121       double coef = trans.coef();
122       int pure = trans.pureindex();
123       int cart = trans.cartindex();
124       int offtarget = ((offsetpure1 + pure)*s2 + offset2)*chunk;
125       int offsource = ((offsetcart1 + cart)*s2 + offset2)*chunk;
126       for (i2=0; i2<n2*chunk; i2++) {
127           target[offtarget++] += coef * source[offsource++];
128         }
129     }
130 }
131 
132 static void
do_sparse_transform12(double * source,double * target,int chunk,SphericalTransformIter & trans,int n1,int offset1,int s2cart,int offsetcart2,int s2pure,int offsetpure2)133 do_sparse_transform12(double *source, double *target, int chunk,
134                       SphericalTransformIter& trans,
135                       int n1, int offset1,
136                       int s2cart, int offsetcart2,
137                       int s2pure, int offsetpure2)
138 {
139   int i1, ichunk;
140 
141   for (trans.begin(); trans.ready(); trans.next()) {
142       double coef = trans.coef();
143       int pure = trans.pureindex();
144       int cart = trans.cartindex();
145       for (i1=0; i1<n1; i1++) {
146           int offtarget = ((offset1 + i1)*s2pure + offsetpure2 + pure)*chunk;
147           int offsource = ((offset1 + i1)*s2cart + offsetcart2 + cart)*chunk;
148           for (ichunk=0; ichunk<chunk; ichunk++) {
149               target[offtarget++] += coef * source[offsource++];
150             }
151         }
152     }
153 }
154 
155 static void
do_sparse_transform2_1(double * source,double * target,SphericalTransformIter & trans,int stcart,int stpure,int ogctcart,int ogctpure,int n2,int s2,int ogc2,int n3,int s3,int ogc3,int n4,int s4,int ogc4)156 do_sparse_transform2_1(double *source, double *target,
157                        SphericalTransformIter& trans,
158                        int stcart, int stpure,
159                        int ogctcart, int ogctpure,
160                        int n2, int s2, int ogc2,
161                        int n3, int s3, int ogc3,
162                        int n4, int s4, int ogc4)
163 {
164   int i2, i3, i4;
165 
166   for (trans.begin(); trans.ready(); trans.next()) {
167       double coef = trans.coef();
168       int pure = trans.pureindex();
169       int cart = trans.cartindex();
170       int offtarget1 = ogctpure + pure;
171       int offsource1 = ogctcart + cart;
172       int offtarget2 = offtarget1*s2 + ogc2;
173       int offsource2 = offsource1*s2 + ogc2;
174       for (i2=0; i2<n2; i2++,offtarget2++,offsource2++) {
175           int offtarget3 = offtarget2*s3 + ogc3;
176           int offsource3 = offsource2*s3 + ogc3;
177           for (i3=0; i3<n3; i3++,offtarget3++,offsource3++) {
178               int offtarget4 = offtarget3*s4 + ogc4;
179               int offsource4 = offsource3*s4 + ogc4;
180               for (i4=0; i4<n4; i4++,offtarget4++,offsource4++) {
181                   target[offtarget4] += coef * source[offsource4];
182                 }
183             }
184         }
185     }
186 }
187 
188 static void
do_sparse_transform2_2(double * source,double * target,SphericalTransformIter & trans,int stcart,int stpure,int ogctcart,int ogctpure,int n1,int s1,int ogc1,int n3,int s3,int ogc3,int n4,int s4,int ogc4)189 do_sparse_transform2_2(double *source, double *target,
190                        SphericalTransformIter& trans,
191                        int stcart, int stpure,
192                        int ogctcart, int ogctpure,
193                        int n1, int s1, int ogc1,
194                        int n3, int s3, int ogc3,
195                        int n4, int s4, int ogc4)
196 {
197   int i1, i3, i4;
198 
199   for (trans.begin(); trans.ready(); trans.next()) {
200       double coef = trans.coef();
201       int pure = trans.pureindex();
202       int cart = trans.cartindex();
203       int offtarget1 = ogc1;
204       int offsource1 = ogc1;
205       for (i1=0; i1<n1; i1++,offtarget1++,offsource1++) {
206           int offtarget2 = offtarget1*stpure + ogctpure + pure;
207           int offsource2 = offsource1*stcart + ogctcart + cart;
208           int offtarget3 = offtarget2*s3 + ogc3;
209           int offsource3 = offsource2*s3 + ogc3;
210           for (i3=0; i3<n3; i3++,offtarget3++,offsource3++) {
211               int offtarget4 = offtarget3*s4 + ogc4;
212               int offsource4 = offsource3*s4 + ogc4;
213               for (i4=0; i4<n4; i4++,offtarget4++,offsource4++) {
214                   target[offtarget4] += coef * source[offsource4];
215                 }
216             }
217         }
218     }
219 }
220 
221 static void
do_sparse_transform2_3(double * source,double * target,SphericalTransformIter & trans,int stcart,int stpure,int ogctcart,int ogctpure,int n1,int s1,int ogc1,int n2,int s2,int ogc2,int n4,int s4,int ogc4)222 do_sparse_transform2_3(double *source, double *target,
223                        SphericalTransformIter& trans,
224                        int stcart, int stpure,
225                        int ogctcart, int ogctpure,
226                        int n1, int s1, int ogc1,
227                        int n2, int s2, int ogc2,
228                        int n4, int s4, int ogc4)
229 {
230   int i1, i2, i4;
231 
232   for (trans.begin(); trans.ready(); trans.next()) {
233       double coef = trans.coef();
234       int pure = trans.pureindex();
235       int cart = trans.cartindex();
236       int offtarget1 = ogc1;
237       int offsource1 = ogc1;
238       for (i1=0; i1<n1; i1++,offtarget1++,offsource1++) {
239           int offtarget2 = offtarget1*s2 + ogc2;
240           int offsource2 = offsource1*s2 + ogc2;
241           for (i2=0; i2<n2; i2++,offtarget2++,offsource2++) {
242               int offtarget3 = offtarget2*stpure + ogctpure + pure;
243               int offsource3 = offsource2*stcart + ogctcart + cart;
244               int offtarget4 = offtarget3*s4 + ogc4;
245               int offsource4 = offsource3*s4 + ogc4;
246               for (i4=0; i4<n4; i4++,offtarget4++,offsource4++) {
247                   target[offtarget4] += coef * source[offsource4];
248                 }
249             }
250         }
251     }
252 }
253 
254 static void
do_sparse_transform2_4(double * source,double * target,SphericalTransformIter & trans,int stcart,int stpure,int ogctcart,int ogctpure,int n1,int s1,int ogc1,int n2,int s2,int ogc2,int n3,int s3,int ogc3)255 do_sparse_transform2_4(double *source, double *target,
256                        SphericalTransformIter& trans,
257                        int stcart, int stpure,
258                        int ogctcart, int ogctpure,
259                        int n1, int s1, int ogc1,
260                        int n2, int s2, int ogc2,
261                        int n3, int s3, int ogc3)
262 {
263   int i1, i2, i3;
264 
265   for (trans.begin(); trans.ready(); trans.next()) {
266       double coef = trans.coef();
267       int pure = trans.pureindex();
268       int cart = trans.cartindex();
269       int offtarget1 = ogc1;
270       int offsource1 = ogc1;
271       for (i1=0; i1<n1; i1++,offtarget1++,offsource1++) {
272           int offtarget2 = offtarget1*s2 + ogc2;
273           int offsource2 = offsource1*s2 + ogc2;
274           for (i2=0; i2<n2; i2++,offtarget2++,offsource2++) {
275               int offtarget3 = offtarget2*s3 + ogc3;
276               int offsource3 = offsource2*s3 + ogc3;
277               int offtarget4 = offtarget3*stpure + ogctpure + pure;
278               int offsource4 = offsource3*stcart + ogctcart + cart;
279               for (i3=0; i3<n3; i3++,offtarget4+=stpure,offsource4+=stcart) {
280                   //for (i3=0; i3<n3; i3++,offtarget3++,offsource3++) {
281                   //int offtarget4 = offtarget3*stpure + ogctpure + pure;
282                   //int offsource4 = offsource3*stcart + ogctcart + cart;
283                   target[offtarget4] += coef * source[offsource4];
284                 }
285             }
286         }
287     }
288 }
289 
290 static void
do_sparse_transform2(double * source,double * target,int index,SphericalTransformIter & trans,int stcart,int stpure,int ogctcart,int ogctpure,int n1,int s1,int ogc1,int n2,int s2,int ogc2,int n3,int s3,int ogc3,int n4,int s4,int ogc4)291 do_sparse_transform2(double *source, double *target,
292                      int index, SphericalTransformIter& trans,
293                      int stcart, int stpure,
294                      int ogctcart, int ogctpure,
295                      int n1, int s1, int ogc1,
296                      int n2, int s2, int ogc2,
297                      int n3, int s3, int ogc3,
298                      int n4, int s4, int ogc4)
299 {
300   switch (index) {
301   case 0:
302       do_sparse_transform2_1(source, target, trans,
303                              stcart,stpure,
304                              ogctcart, ogctpure,
305                              n2, s2, ogc2,
306                              n3, s3, ogc3,
307                              n4, s4, ogc4);
308       break;
309   case 1:
310       do_sparse_transform2_2(source, target, trans,
311                              stcart,stpure,
312                              ogctcart, ogctpure,
313                              n1, s1, ogc1,
314                              n3, s3, ogc3,
315                              n4, s4, ogc4);
316       break;
317   case 2:
318       do_sparse_transform2_3(source, target, trans,
319                              stcart,stpure,
320                              ogctcart, ogctpure,
321                              n1, s1, ogc1,
322                              n2, s2, ogc2,
323                              n4, s4, ogc4);
324       break;
325   case 3:
326       do_sparse_transform2_4(source, target, trans,
327                              stcart,stpure,
328                              ogctcart, ogctpure,
329                              n1, s1, ogc1,
330                              n2, s2, ogc2,
331                              n3, s3, ogc3);
332       break;
333     }
334 }
335 
336 /* make sure enough space exists for the source integrals */
337 void
source_space(int nsource)338 Int2eV3::source_space(int nsource)
339 {
340   if (nsourcemax < nsource) {
341       delete[] source;
342       source = new double[nsource*3];
343       target = &source[nsource];
344       scratch = &source[nsource*2];
345       nsourcemax = nsource;
346     }
347 }
348 
349 void
copy_to_source(double * integrals,int nsource)350 Int2eV3::copy_to_source(double *integrals, int nsource)
351 {
352   int i;
353   double *tmp, *tmp2;
354 
355   /* Allocate more temporary space if needed. */
356   source_space(nsource);
357 
358   tmp = source;
359   tmp2 = integrals;
360   for (i=0; i<nsource; i++) *tmp++ = *tmp2++;
361 }
362 
363 /* make sure enough space exists for the source integrals */
364 void
source_space(int nsource)365 Int1eV3::source_space(int nsource)
366 {
367   if (nsourcemax < nsource) {
368       delete[] source;
369       source = new double[nsource];
370       nsourcemax = nsource;
371     }
372 }
373 
374 void
copy_to_source(double * integrals,int nsource)375 Int1eV3::copy_to_source(double *integrals, int nsource)
376 {
377   int i;
378   double *tmp, *tmp2;
379 
380   /* Allocate more temporary space if needed. */
381   source_space(nsource);
382 
383   tmp = source;
384   tmp2 = integrals;
385   for (i=0; i<nsource; i++) *tmp++ = *tmp2++;
386 }
387 
388 void
do_transform_1e(Integral * integ,double * integrals,GaussianShell * sh1,GaussianShell * sh2,int chunk)389 Int1eV3::do_transform_1e(Integral *integ,
390                          double *integrals,
391                          GaussianShell *sh1, GaussianShell *sh2,
392                          int chunk)
393 {
394   int i, j;
395   int ogc1, ogc2;
396   int ogc1pure, ogc2pure;
397   int am1, am2;
398   int pure1 = sh1->has_pure();
399   int pure2 = sh2->has_pure();
400   int ncart1 = sh1->ncartesian();
401   int ncart2 = sh2->ncartesian();
402   int nfunc1 = sh1->nfunction();
403   int nfunc2 = sh2->nfunction();
404   int nfunci, nfuncj;
405 
406   if (!pure1 && !pure2) return;
407 
408   /* Loop through the generalized general contractions,
409    * transforming the first index. */
410   if (pure1) {
411       copy_to_source(integrals, ncart1*ncart2*chunk);
412       memset(integrals, 0, sizeof(double)*sh1->nfunction()*ncart2*chunk);
413 
414       ogc1 = 0;
415       ogc1pure = 0;
416       for (i=0; i<sh1->ncontraction(); i++) {
417           am1 = sh1->am(i);
418           nfunci = sh1->nfunction(i);
419           ogc2 = 0;
420           for (j=0; j<sh2->ncontraction(); j++) {
421               am2 = sh2->am(j);
422               nfuncj = sh2->nfunction(j);
423 
424               if (sh1->is_pure(i)) {
425                   SphericalTransformIter
426                       trans(integ->spherical_transform(sh1->am(i)));
427                   do_sparse_transform11(source, integrals, chunk,
428                                         trans,
429                                         ogc1,
430                                         ogc1pure,
431                                         INT_NCART(am2), ncart2, ogc2);
432                 }
433               else {
434                   do_copy1(source, integrals, chunk,
435                            nfunci, nfunc1, ogc1pure,
436                            INT_NCART(am2), ncart2, ogc2);
437                 }
438               ogc2 += INT_NCART(am2);
439             }
440           ogc1 += INT_NCART(am1);
441           ogc1pure += INT_NPURE(am1);
442         }
443     }
444 
445   if (pure2) {
446       copy_to_source(integrals, nfunc1*ncart2*chunk);
447       memset(integrals, 0,
448              sizeof(double)*sh1->nfunction()*sh2->nfunction()*chunk);
449 
450       ogc1 = 0;
451       for (i=0; i<sh1->ncontraction(); i++) {
452           am1 = sh1->am(i);
453           nfunci = sh1->nfunction(i);
454           ogc2 = 0;
455           ogc2pure = 0;
456           for (j=0; j<sh2->ncontraction(); j++) {
457               am2 = sh2->am(j);
458               nfuncj = sh2->nfunction(j);
459 
460               if (sh2->is_pure(j)) {
461                   SphericalTransformIter
462                       trans(integ->spherical_transform(sh2->am(j)));
463                   do_sparse_transform12(source, integrals, chunk,
464                                         trans,
465                                         INT_NPURE(am1), ogc1,
466                                         ncart2, ogc2,
467                                         sh2->nfunction(), ogc2pure);
468                 }
469               else {
470                   do_copy1(source, integrals, chunk,
471                            nfunci, nfunc1, ogc1,
472                            nfuncj, nfunc2, ogc2pure);
473                 }
474               ogc2 += INT_NCART(am2);
475               ogc2pure += INT_NPURE(am2);
476             }
477           ogc1 += INT_NPURE(am1);
478         }
479     }
480 }
481 
482 /* it is ok for integrals and target to overlap */
483 void
transform_1e(Integral * integ,double * integrals,double * target,GaussianShell * sh1,GaussianShell * sh2,int chunk)484 Int1eV3::transform_1e(Integral *integ,
485                       double *integrals, double *target,
486                       GaussianShell *sh1, GaussianShell *sh2, int chunk)
487 {
488   int ntarget;
489 
490   do_transform_1e(integ, integrals, sh1, sh2, chunk);
491 
492   /* copy the integrals to the target, if necessary */
493   ntarget = sh1->nfunction() * sh2->nfunction();
494   if (integrals != target) {
495       memmove(target, integrals, ntarget*sizeof(double)*chunk);
496     }
497 }
498 
499 /* it is not ok for integrals and target to overlap */
500 void
accum_transform_1e(Integral * integ,double * integrals,double * target,GaussianShell * sh1,GaussianShell * sh2,int chunk)501 Int1eV3::accum_transform_1e(Integral *integ,
502                             double *integrals, double *target,
503                             GaussianShell *sh1, GaussianShell *sh2, int chunk)
504 {
505   int i, ntarget;
506 
507   do_transform_1e(integ, integrals, sh1, sh2, chunk);
508 
509   /* accum the integrals to the target */
510   ntarget = sh1->nfunction() * sh2->nfunction() * chunk;
511   for (i=0; i<ntarget; i++) target[i] += integrals[i];
512 }
513 
514 void
transform_1e(Integral * integ,double * integrals,double * target,GaussianShell * sh1,GaussianShell * sh2)515 Int1eV3::transform_1e(Integral*integ,
516                             double *integrals, double *target,
517                             GaussianShell *sh1, GaussianShell *sh2)
518 {
519   transform_1e(integ, integrals, target, sh1, sh2, 1);
520 }
521 
522 void
accum_transform_1e(Integral * integ,double * integrals,double * target,GaussianShell * sh1,GaussianShell * sh2)523 Int1eV3::accum_transform_1e(Integral*integ,
524                                   double *integrals, double *target,
525                                   GaussianShell *sh1, GaussianShell *sh2)
526 {
527   accum_transform_1e(integ, integrals, target, sh1, sh2, 1);
528 }
529 
530 void
transform_1e_xyz(Integral * integ,double * integrals,double * target,GaussianShell * sh1,GaussianShell * sh2)531 Int1eV3::transform_1e_xyz(Integral*integ,
532                           double *integrals, double *target,
533                           GaussianShell *sh1, GaussianShell *sh2)
534 {
535   transform_1e(integ, integrals, target, sh1, sh2, 3);
536 }
537 
538 void
accum_transform_1e_xyz(Integral * integ,double * integrals,double * target,GaussianShell * sh1,GaussianShell * sh2)539 Int1eV3::accum_transform_1e_xyz(Integral*integ,
540                                 double *integrals, double *target,
541                                 GaussianShell *sh1, GaussianShell *sh2)
542 {
543   accum_transform_1e(integ, integrals, target, sh1, sh2, 3);
544 }
545 
546 void
do_gencon_sparse_transform_2e(Integral * integ,double * integrals,double * target,int index,GaussianShell * sh1,GaussianShell * sh2,GaussianShell * sh3,GaussianShell * sh4)547 Int2eV3::do_gencon_sparse_transform_2e(Integral*integ,
548                                        double *integrals, double *target,
549                                        int index,
550                                        GaussianShell *sh1, GaussianShell *sh2,
551                                        GaussianShell *sh3, GaussianShell *sh4)
552 {
553   int ncart[4];
554   int nfunc[4];
555   int nfunci, nfuncj, nfunck, nfuncl;
556   int ncarti, ncartj, ncartk, ncartl;
557   int i, j, k, l;
558   int ogccart[4];
559   int ogcfunc[4];
560   int am1, am2, am3, am4;
561 
562   int ntarget1;
563   int ntarget2;
564   int ntarget3;
565   int ntarget4;
566 
567   int nsource1;
568   int nsource2;
569   int nsource3;
570   int nsource4;
571 
572   int *ni = &ncarti;
573   int *nj = &ncartj;
574   int *nk = &ncartk;
575   int *nl = &ncartl;
576 
577   int *ogc1 = &ogccart[0];
578   int *ogc2 = &ogccart[1];
579   int *ogc3 = &ogccart[2];
580   int *ogc4 = &ogccart[3];
581 
582   GaussianShell *shell;
583 
584   int *tgencon;
585 
586   ncart[0] = sh1->ncartesian();
587   ncart[1] = sh2->ncartesian();
588   ncart[2] = sh3->ncartesian();
589   ncart[3] = sh4->ncartesian();
590   nfunc[0] = sh1->nfunction();
591   nfunc[1] = sh2->nfunction();
592   nfunc[2] = sh3->nfunction();
593   nfunc[3] = sh4->nfunction();
594 
595   ntarget1 = ncart[0];
596   ntarget2 = ncart[1];
597   ntarget3 = ncart[2];
598   ntarget4 = ncart[3];
599 
600   nsource1 = ncart[0];
601   nsource2 = ncart[1];
602   nsource3 = ncart[2];
603   nsource4 = ncart[3];
604 
605   if (index >= 0) {
606       ntarget1 = nfunc[0];
607       if (index >= 1) {
608           ntarget2 = nfunc[1];
609           nsource1 = nfunc[0];
610           ni = &nfunci;
611           ogc1 = &ogcfunc[0];
612           if (index >= 2) {
613               ntarget3 = nfunc[2];
614               nsource2 = nfunc[1];
615               nj = &nfuncj;
616               ogc2 = &ogcfunc[1];
617               if (index >= 3) {
618                   ntarget4 = nfunc[3];
619                   nsource3 = nfunc[2];
620                   nk = &nfunck;
621                   ogc3 = &ogcfunc[2];
622                 }
623             }
624         }
625     }
626 
627   switch (index) {
628   case 0:
629       shell = sh1;
630       tgencon = &i;
631       break;
632   case 1:
633       shell = sh2;
634       tgencon = &j;
635       break;
636   case 2:
637       shell = sh3;
638       tgencon = &k;
639       break;
640   case 3:
641       shell = sh4;
642       tgencon = &l;
643       break;
644   default:
645       shell = 0;
646       tgencon = 0;
647       break;
648     }
649 
650 #if PRINT
651     {
652       double *tmp = integrals;
653       ExEnv::outn() << scprintf("Before transform of index %d (%dx%dx%dx%d)\n",
654              index, nsource1, nsource2, nsource3, nsource4);
655       for (i=0; i<nsource1; i++) {
656           for (j=0; j<nsource2; j++) {
657               for (k=0; k<nsource3; k++) {
658                   for (l=0; l<nsource4; l++) {
659                       if (fabs(*tmp)>1.e-15) {
660                           ExEnv::outn() << scprintf("(%d %d|%d %d) = %15.11lf\n",i,j,k,l,*tmp);
661                         }
662                       tmp++;
663                     }
664                 }
665             }
666         }
667     }
668 #endif
669 
670   copy_to_source(integrals, nsource1*nsource2*nsource3*nsource4);
671   memset(target, 0, sizeof(double)*ntarget1*ntarget2*ntarget3*ntarget4);
672 
673   ogccart[0] = 0;
674   ogcfunc[0] = 0;
675   for (i=0; i<sh1->ncontraction(); i++) {
676       am1 = sh1->am(i);
677       nfunci = sh1->nfunction(i);
678       ncarti = INT_NCART(am1);
679       ogccart[1] = 0;
680       ogcfunc[1] = 0;
681       for (j=0; j<sh2->ncontraction(); j++) {
682           am2 = sh2->am(j);
683           nfuncj = sh2->nfunction(j);
684           ncartj = INT_NCART(am2);
685           ogccart[2] = 0;
686           ogcfunc[2] = 0;
687           for (k=0; k<sh3->ncontraction(); k++) {
688               am3 = sh3->am(k);
689               nfunck = sh3->nfunction(k);
690               ncartk = INT_NCART(am3);
691               ogccart[3] = 0;
692               ogcfunc[3] = 0;
693               for (l=0; l<sh4->ncontraction(); l++) {
694                   am4 = sh4->am(l);
695                   nfuncl = sh4->nfunction(l);
696                   ncartl = INT_NCART(am4);
697 
698                   if (shell->is_pure(*tgencon)) {
699                       SphericalTransformIter
700                         trans(integ->spherical_transform(shell->am(*tgencon)));
701                       do_sparse_transform2(source, target,
702                                            index, trans,
703                                            ncart[index], nfunc[index],
704                                            ogccart[index], ogcfunc[index],
705                                            *ni, nsource1, *ogc1,
706                                            *nj, nsource2, *ogc2,
707                                            *nk, nsource3, *ogc3,
708                                            *nl, nsource4, *ogc4);
709                     }
710                   else {
711                       do_copy2(source, integrals,
712                                *ni, nsource1, *ogc1,
713                                *nj, nsource2, *ogc2,
714                                *nk, nsource3, *ogc3,
715                                *nl, nsource4, *ogc4);
716                     }
717                   ogccart[3] += ncartl;
718                   ogcfunc[3] += nfuncl;
719                 }
720               ogccart[2] += ncartk;
721               ogcfunc[2] += nfunck;
722             }
723           ogccart[1] += ncartj;
724           ogcfunc[1] += nfuncj;
725         }
726       ogccart[0] += ncarti;
727       ogcfunc[0] += nfunci;
728     }
729 
730 
731 #if PRINT
732     {
733       double *tmp = integrals;
734       ExEnv::outn() << scprintf("After transform of index %d (%dx%dx%dx%d)\n",
735              index, ntarget1, ntarget2, ntarget3, ntarget4);
736       for (i=0; i<ntarget1; i++) {
737           for (j=0; j<ntarget2; j++) {
738               for (k=0; k<ntarget3; k++) {
739                   for (l=0; l<ntarget4; l++) {
740                       if (fabs(*tmp)>1.e-15) {
741                           ExEnv::outn()
742                             << scprintf("(%d %d|%d %d) = %15.11lf\n",
743                                         i,j,k,l,*tmp);
744                         }
745                       tmp++;
746                     }
747                 }
748             }
749         }
750     }
751 #endif
752 }
753 
754 void
transform_2e_slow(Integral * integ,double * integrals,double * target,GaussianShell * sh1,GaussianShell * sh2,GaussianShell * sh3,GaussianShell * sh4)755 Int2eV3::transform_2e_slow(Integral *integ, double *integrals, double *target,
756                            GaussianShell *sh1, GaussianShell *sh2,
757                            GaussianShell *sh3, GaussianShell *sh4)
758 {
759   int pure1 = sh1->has_pure();
760   int pure2 = sh2->has_pure();
761   int pure3 = sh3->has_pure();
762   int pure4 = sh4->has_pure();
763 
764   if (pure1) {
765       do_gencon_sparse_transform_2e(integ,
766                                     integrals, target, 0, sh1, sh2, sh3, sh4);
767       integrals = target;
768     }
769   if (pure2) {
770       do_gencon_sparse_transform_2e(integ,
771                                     integrals, target, 1, sh1, sh2, sh3, sh4);
772       integrals = target;
773     }
774   if (pure3) {
775       do_gencon_sparse_transform_2e(integ,
776                                     integrals, target, 2, sh1, sh2, sh3, sh4);
777       integrals = target;
778     }
779   if (pure4) {
780       do_gencon_sparse_transform_2e(integ,
781                                     integrals, target, 3, sh1, sh2, sh3, sh4);
782       integrals = target;
783     }
784 
785   if (integrals != target) {
786       int nint = sh1->nfunction() * sh2->nfunction()
787                * sh3->nfunction() * sh4->nfunction();
788       memmove(target, integrals, sizeof(double)*nint);
789     }
790 }
791 
792 /////////////////////////////////////////////////////////////////////////////
793 
794 static void
do_sparse_transform2_1new(double * source,double * target,SphericalTransformIter & trans,int stcart,int stpure,int n2,int n3,int n4)795 do_sparse_transform2_1new(double *source, double *target,
796                           SphericalTransformIter& trans,
797                           int stcart, int stpure,
798                           int n2,
799                           int n3,
800                           int n4)
801 {
802   int i234, n234=n2*n3*n4;
803 
804   for (trans.begin(); trans.ready(); trans.next()) {
805     double coef = trans.coef();
806     int pure = trans.pureindex();
807     int cart = trans.cartindex();
808     int offtarget4 = pure*n234;
809     int offsource4 = cart*n234;
810     for (i234=0; i234<n234; i234++,offtarget4++,offsource4++) {
811       target[offtarget4] += coef * source[offsource4];
812       }
813     }
814 }
815 
816 static void
do_sparse_transform2_2new(double * source,double * target,SphericalTransformIter & trans,int stcart,int stpure,int n1,int n3,int n4)817 do_sparse_transform2_2new(double *source, double *target,
818                           SphericalTransformIter& trans,
819                           int stcart, int stpure,
820                           int n1,
821                           int n3,
822                           int n4)
823 {
824   int i1, i34, n34=n3*n4;
825   int n34stpure=n34*(stpure-1); // -1 because of the increment int the loop
826   int n34stcart=n34*(stcart-1); // ditto
827 
828   for (trans.begin(); trans.ready(); trans.next()) {
829     double coef = trans.coef();
830     int pure = trans.pureindex();
831     int cart = trans.cartindex();
832     int offtarget4 = pure*n34;
833     int offsource4 = cart*n34;
834     for (i1=0; i1<n1; i1++) {
835       for (i34=0; i34<n34; i34++,offtarget4++,offsource4++) {
836         target[offtarget4] += coef * source[offsource4];
837         }
838       offtarget4 += n34stpure;
839       offsource4 += n34stcart;
840       }
841     }
842 }
843 
844 static void
do_sparse_transform2_3new(double * source,double * target,SphericalTransformIter & trans,int stcart,int stpure,int n1,int n2,int n4)845 do_sparse_transform2_3new(double *source, double *target,
846                           SphericalTransformIter& trans,
847                           int stcart, int stpure,
848                           int n1,
849                           int n2,
850                           int n4)
851 {
852   int i12, i4, n12=n1*n2, n4stpure=n4*stpure, n4stcart=n4*stcart;
853 
854   for (trans.begin(); trans.ready(); trans.next()) {
855     double coef = trans.coef();
856     int pure = trans.pureindex();
857     int cart = trans.cartindex();
858     int offtarget3 = pure*n4;
859     int offsource3 = cart*n4;
860     for (i12=0; i12<n12; i12++) {
861       int offtarget4 = offtarget3;
862       int offsource4 = offsource3;
863       for (i4=0; i4<n4; i4++,offtarget4++,offsource4++) {
864         target[offtarget4] += coef * source[offsource4];
865         }
866       offtarget3 += n4stpure;
867       offsource3 += n4stcart;
868       }
869     }
870 }
871 
872 static void
do_sparse_transform2_4new(double * source,double * target,SphericalTransformIter & trans,int stcart,int stpure,int n1,int n2,int n3)873 do_sparse_transform2_4new(double *source, double *target,
874                           SphericalTransformIter& trans,
875                           int stcart, int stpure,
876                           int n1,
877                           int n2,
878                           int n3)
879 {
880   int n123=n1*n2*n3;
881 
882   for (trans.begin(); trans.ready(); trans.next()) {
883     double coef = trans.coef();
884     int pure = trans.pureindex();
885     int cart = trans.cartindex();
886     int offtarget4 = pure;
887     int offsource4 = cart;
888     for (int i123=0; i123<n123; i123++) {
889       target[offtarget4] += coef * source[offsource4];
890       offtarget4 += stpure;
891       offsource4 += stcart;
892       }
893     }
894 }
895 
896 // Cartint and pureint may overlap.  The must be enough space
897 // in pureint to hold all the cartints.  The cartint buffer
898 // will be overwritten in any case.
899 void
transform_2e(Integral * integ,double * cartint,double * pureint,GaussianShell * sh1,GaussianShell * sh2,GaussianShell * sh3,GaussianShell * sh4)900 Int2eV3::transform_2e(Integral *integ,
901                       double *cartint, double *pureint,
902                       GaussianShell *sh1, GaussianShell *sh2,
903                       GaussianShell *sh3, GaussianShell *sh4)
904 {
905   int pure1 = sh1->has_pure();
906   int pure2 = sh2->has_pure();
907   int pure3 = sh3->has_pure();
908   int pure4 = sh4->has_pure();
909 
910   int nfunc1=sh1->nfunction();
911   int nfunc2=sh2->nfunction();
912   int nfunc3=sh3->nfunction();
913   int nfunc4=sh4->nfunction();
914   int nfunc34 =  nfunc3 * nfunc4;
915   int nfunc234 = nfunc2 * nfunc34;
916   int nfunc1234 = nfunc1 * nfunc234;
917 
918   if (!pure1 && !pure2 && !pure3 && !pure4) {
919     if (pureint!=cartint) memmove(pureint, cartint, sizeof(double)*nfunc1234);
920     return;
921     }
922 
923   int ncart1=sh1->ncartesian();
924   int ncart2=sh2->ncartesian();
925   int ncart3=sh3->ncartesian();
926   int ncart4=sh4->ncartesian();
927   int ncart34 =  ncart3 * ncart4;
928   int ncart234 = ncart2 * ncart34;
929 
930   // allocate the scratch arrays, if needed
931   source_space(ncart1*ncart234);
932 
933   int ncon1 = sh1->ncontraction();
934   int ncon2 = sh2->ncontraction();
935   int ncon3 = sh3->ncontraction();
936   int ncon4 = sh4->ncontraction();
937 
938   if (ncon1==1 && ncon2==1 && ncon3==1 && ncon4==1) {
939     double *sourcebuf = cartint;
940     double *targetbuf = target;
941     // transform indices
942     if (pure1) {
943       SphericalTransformIter transi(integ->spherical_transform(sh1->am(0)));
944       memset(targetbuf,0,sizeof(double)*nfunc1*ncart2*ncart3*ncart4);
945       do_sparse_transform2_1new(sourcebuf, targetbuf, transi,
946                                 ncart1, nfunc1,
947                                 ncart2,
948                                 ncart3,
949                                 ncart4);
950       double*tmp=sourcebuf; sourcebuf=targetbuf; targetbuf=tmp;
951       }
952     if (pure2) {
953       SphericalTransformIter transj(integ->spherical_transform(sh2->am(0)));
954       memset(targetbuf,0,sizeof(double)*nfunc1*nfunc2*ncart3*ncart4);
955       do_sparse_transform2_2new(sourcebuf, targetbuf, transj,
956                                 ncart2, nfunc2,
957                                 nfunc1,
958                                 ncart3,
959                                 ncart4);
960       double*tmp=sourcebuf; sourcebuf=targetbuf; targetbuf=tmp;
961       }
962     if (pure3) {
963       SphericalTransformIter transk(integ->spherical_transform(sh3->am(0)));
964       memset(targetbuf,0,sizeof(double)*nfunc1*nfunc2*nfunc3*ncart4);
965       do_sparse_transform2_3new(sourcebuf, targetbuf, transk,
966                                 ncart3, nfunc3,
967                                 nfunc1,
968                                 nfunc2,
969                                 ncart4);
970       double*tmp=sourcebuf; sourcebuf=targetbuf; targetbuf=tmp;
971       }
972     if (pure4) {
973       SphericalTransformIter transl(integ->spherical_transform(sh4->am(0)));
974       memset(targetbuf,0,sizeof(double)*nfunc1234);
975       do_sparse_transform2_4new(sourcebuf, targetbuf, transl,
976                                 ncart4, nfunc4,
977                                 nfunc1,
978                                 nfunc2,
979                                 nfunc3);
980       double*tmp=sourcebuf; sourcebuf=targetbuf; targetbuf=tmp;
981       }
982     if (sourcebuf!=pureint)
983       memmove(pureint, sourcebuf, sizeof(double)*nfunc1234);
984     }
985   else {
986     // begin gc loop
987     int ogccart1 = 0;
988     int ogcfunc1 = 0;
989     for (int i=0; i<ncon1; i++) {
990       int am1 = sh1->am(i);
991       int nfunci = sh1->nfunction(i);
992       int ispurei = sh1->is_pure(i);
993       int ncarti = INT_NCART_NN(am1);
994       int ogccart2 = 0;
995       int ogcfunc2 = 0;
996       SphericalTransformIter transi(integ->spherical_transform(am1));
997       for (int j=0; j<ncon2; j++) {
998         int am2 = sh2->am(j);
999         int nfuncj = sh2->nfunction(j);
1000         int ispurej = sh2->is_pure(j);
1001         int ncartj = INT_NCART_NN(am2);
1002         int ogccart3 = 0;
1003         int ogcfunc3 = 0;
1004         SphericalTransformIter transj(integ->spherical_transform(am2));
1005         for (int k=0; k<ncon3; k++) {
1006           int am3 = sh3->am(k);
1007           int nfunck = sh3->nfunction(k);
1008           int ispurek = sh3->is_pure(k);
1009           int ncartk = INT_NCART_NN(am3);
1010           int ogccart4 = 0;
1011           int ogcfunc4 = 0;
1012           SphericalTransformIter transk(integ->spherical_transform(am3));
1013           for (int l=0; l<ncon4; l++) {
1014             int am4 = sh4->am(l);
1015             int nfuncl = sh4->nfunction(l);
1016             int ispurel = sh4->is_pure(l);
1017             int ncartl = INT_NCART_NN(am4);
1018 
1019     ;
1020     // copy to source buffer
1021     int cartindex1 = ogccart1*ncart234
1022                    + ogccart2*ncart34 + ogccart3*ncart4 + ogccart4;
1023     double *tmp_source = source;
1024     int is;
1025     for (is=0; is<ncarti; is++) {
1026       int cartindex12 = cartindex1;
1027       for (int js=0; js<ncartj; js++) {
1028         int cartindex123 = cartindex12;
1029         for (int ks=0; ks<ncartk; ks++) {
1030           double *tmp_cartint = &cartint[cartindex123];
1031           for (int ls=0; ls<ncartl; ls++) {
1032             *tmp_source++ = *tmp_cartint++;
1033             }
1034           cartindex123 += ncart4;
1035           }
1036         cartindex12 += ncart34;
1037         }
1038       cartindex1 += ncart234;
1039       }
1040 
1041 
1042     double *sourcebuf = source;
1043     double *targetbuf = target;
1044 
1045     // transform indices
1046     if (ispurei) {
1047       memset(targetbuf,0,sizeof(double)*nfunci*ncartj*ncartk*ncartl);
1048       do_sparse_transform2_1new(sourcebuf, targetbuf, transi,
1049                                 ncarti, nfunci,
1050                                 ncartj,
1051                                 ncartk,
1052                                 ncartl);
1053       double*tmp=sourcebuf; sourcebuf=targetbuf; targetbuf=tmp;
1054       }
1055     if (ispurej) {
1056       memset(targetbuf,0,sizeof(double)*nfunci*nfuncj*ncartk*ncartl);
1057       do_sparse_transform2_2new(sourcebuf, targetbuf, transj,
1058                                 ncartj, nfuncj,
1059                                 nfunci,
1060                                 ncartk,
1061                                 ncartl);
1062       double*tmp=sourcebuf; sourcebuf=targetbuf; targetbuf=tmp;
1063       }
1064     if (ispurek) {
1065       memset(targetbuf,0,sizeof(double)*nfunci*nfuncj*nfunck*ncartl);
1066       do_sparse_transform2_3new(sourcebuf, targetbuf, transk,
1067                                 ncartk, nfunck,
1068                                 nfunci,
1069                                 nfuncj,
1070                                 ncartl);
1071       double*tmp=sourcebuf; sourcebuf=targetbuf; targetbuf=tmp;
1072       }
1073     if (ispurel) {
1074       memset(targetbuf,0,sizeof(double)*nfunci*nfuncj*nfunck*nfuncl);
1075       SphericalTransformIter transl(integ->spherical_transform(am4));
1076       do_sparse_transform2_4new(sourcebuf, targetbuf, transl,
1077                                 ncartl, nfuncl,
1078                                 nfunci,
1079                                 nfuncj,
1080                                 nfunck);
1081       double*tmp=sourcebuf; sourcebuf=targetbuf; targetbuf=tmp;
1082       }
1083 
1084     // copy to scratch buffer
1085     int funcindex1 = ogcfunc1*nfunc234
1086                    + ogcfunc2*nfunc34 + ogcfunc3*nfunc4 + ogcfunc4;
1087     tmp_source = sourcebuf;
1088     for (is=0; is<nfunci; is++) {
1089       int funcindex12 = funcindex1;
1090       for (int js=0; js<nfuncj; js++) {
1091         int funcindex123 = funcindex12;
1092         for (int ks=0; ks<nfunck; ks++) {
1093           double *tmp_scratch = &scratch[funcindex123];
1094           for (int ls=0; ls<nfuncl; ls++) {
1095             *tmp_scratch++ = *tmp_source++;
1096             }
1097           funcindex123 += nfunc4;
1098           }
1099         funcindex12 += nfunc34;
1100         }
1101       funcindex1 += nfunc234;
1102       }
1103 
1104     // end gc loop
1105             ogccart4 += ncartl;
1106             ogcfunc4 += nfuncl;
1107             }
1108           ogccart3 += ncartk;
1109           ogcfunc3 += nfunck;
1110           }
1111         ogccart2 += ncartj;
1112         ogcfunc2 += nfuncj;
1113         }
1114       ogccart1 += ncarti;
1115       ogcfunc1 += nfunci;
1116       }
1117     memcpy(pureint, scratch, sizeof(double)*nfunc1234);
1118     }
1119 }
1120 
1121 /////////////////////////////////////////////////////////////////////////////
1122 
1123 // Local Variables:
1124 // mode: c++
1125 // c-file-style: "CLJ-CONDENSED"
1126 // End:
1127