source: ThirdParty/mpqc_open/src/lib/chemistry/qc/intv3/build2e.cc@ 1ca493a

Candidate_v1.7.0 stable
Last change on this file since 1ca493a was 860145, checked in by Frederik Heber <heber@…>, 9 years ago

Merge commit '0b990dfaa8c6007a996d030163a25f7f5fc8a7e7' as 'ThirdParty/mpqc_open'
  • Property mode set to 100644
File size: 49.5 KB
RevLine 
[0b990d]1//
2// build2e.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 <stdlib.h>
29#include <assert.h>
30#include <math.h>
31
32#include <scconfig.h>
33#include <util/misc/formio.h>
34#include <chemistry/qc/intv3/macros.h>
35#include <chemistry/qc/intv3/fjt.h>
36#include <chemistry/qc/intv3/utils.h>
37#include <chemistry/qc/intv3/int2e.h>
38
39using namespace std;
40using namespace sc;
41
42#define CHECK_STACK_ALIGNMENT 0
43#if CHECK_STACK_ALIGNMENT
44static void
45stack_alignment_error(void *ptr, const char *where)
46{
47 ExEnv::outn() << "UNALIGNED STACK: " << where << ": " << ptr << endl;
48}
49static inline void
50stack_alignment_check(void *ptr, const char *where)
51{
52 if ((unsigned)ptr & 7) stack_alignment_error(ptr,where);
53}
54#else
55# define stack_alignment_check(ptr,where)
56#endif
57
58 /* MG is the maximum angular momentum for which we will use
59 * the generated build routines. It is defined in oint3/build.h */
60#define MINA(x) (((x)<MG)?(x):MG)
61
62static inline void
63iswtch(int *i, int *j)
64{
65 int tmp;
66
67 tmp = *i;
68 *i = *j;
69 *j = tmp;
70}
71
72static void
73fail()
74{
75 ExEnv::errn() << scprintf("failing module:\n%s",__FILE__) << endl;
76 abort();
77 }
78
79/* This initializes the build routines. It is called from
80 * int_initialize_erep. This allocates storage for the
81 * intermediate integrals. */
82void
83Int2eV3::int_init_buildgc(int order,
84 int am1, int am2, int am3, int am4,
85 int nc1, int nc2, int nc3, int nc4)
86{
87 int *jmax_for_con;
88 int am12;
89 int am34;
90 int am;
91 int i,j,k,l,m;
92 int ci,cj,ck,cl;
93 int e0f0_con_int_bufsize;
94 double *e0f0_con_int_buf;
95 int int_v_bufsize, int_v0_bufsize;
96 double *int_v_buf, *int_v0_buf;
97
98 used_storage_build_ = 0;
99
100 /* Convert the am1-4 to their canonical ordering. */
101 if (am2>am1) {
102 iswtch(&am1,&am2);
103 iswtch(&nc1,&nc2);
104 }
105 if (am4>am3) {
106 iswtch(&am3,&am4);
107 iswtch(&nc3,&nc4);
108 }
109 if ((am3 > am1)||((am3 == am1)&&(am4 > am2))) {
110 iswtch(&am1,&am3);
111 iswtch(&nc1,&nc3);
112 iswtch(&am2,&am4);
113 iswtch(&nc2,&nc4);
114 }
115
116 /* If the center permutation 1<->3 and 2<->4 is performed, then
117 * we may need the am for center 2 to be as big as for center 4. */
118 if (am4 > am2) am2 = am4;
119
120 /* As far as this routine knows the biggest nc can end up anywhere. */
121 if (nc2>nc1) nc1 = nc2;
122 if (nc3>nc1) nc1 = nc3;
123 if (nc4>nc1) nc1 = nc4;
124 nc2 = nc3 = nc4 = nc1;
125
126 jmax_for_con = (int *) malloc(sizeof(int)*nc1);
127 // storage for jmax_for_con is not counted since it is freed below
128 for (i=0; i<nc1; i++) {
129 int tmp;
130 jmax_for_con[i] = bs1_->max_am_for_contraction(i);
131 if ( (bs2_ != bs1_)
132 &&((tmp=(int_unit2?0:bs2_->max_am_for_contraction(i)))>jmax_for_con[i]))
133 jmax_for_con[i] = tmp;
134 if ( (bs3_ != bs1_) && (bs3_ != bs2_)
135 &&((tmp=bs3_->max_am_for_contraction(i))>jmax_for_con[i]))
136 jmax_for_con[i] = tmp;
137 if ( (bs4_ != bs1_) && (bs4_ != bs2_) && (bs4_ != bs3_)
138 &&((tmp=(int_unit4?0:bs4_->max_am_for_contraction(i)))>jmax_for_con[i]))
139 jmax_for_con[i] = tmp;
140 }
141
142 /* If derivatives are needed, then am1 can be bigger. */
143 if (order==1) am1++;
144 /* To compute derivative integral bounds, am3 can be bigger also. */
145 if (order==1 && int_derivative_bounds) am3++;
146
147 am12 = am1 + am2;
148 am34 = am3 + am4;
149 am = am12 + am34;
150
151 /* Allocate the intlist. */
152 contract_length.set_dim(am12+1,am34+1,am34+1);
153 build.int_v_list.set_dim(am12+1,am34+1,am+1);
154 used_storage_build_ += contract_length.nbyte();
155 used_storage_build_ += build.int_v_list.nbyte();
156#if CHECK_INTEGRAL_ALGORITHM
157 ExEnv::outn() << "contract_length: " << contract_length.nbyte() << endl;
158 ExEnv::outn() << "int_v_list: " << build.int_v_list.nbyte() << endl;
159#endif
160
161 /* Set all slots to 0 */
162 for (i=0; i<=am12; i++) {
163 for (j=0; j<=am34; j++) {
164 for (k=0; k<=am12+am34; k++) {
165 build.int_v_list(i,j,k) = 0;
166 }
167 }
168 }
169
170 for (i=0; i<=am12; i++) {
171 for (j=0; j<=am34; j++) {
172 for (k=0; k<=am34; k++) {
173 contract_length(i,j,k) = 0;
174 for (l=j; l<=k; l++) {
175 contract_length(i,j,k) += INT_NCART(i)*INT_NCART(l);
176 }
177 }
178 }
179 }
180
181 /* Compute the size of the buffer for the primitive integrals. */
182 int_v_bufsize = 0;
183 int_v0_bufsize = 0;
184 for (i=0; i<=am12; i++) {
185 for (j=0; j<=am34; j++) {
186 int_v0_bufsize += INT_NCART(i)*INT_NCART(j);
187 for (k=0; k<=am12+am34-i-j; k++) {
188 int_v_bufsize += INT_NCART(i)*INT_NCART(j);
189 }
190 }
191 }
192
193 int_v0_buf = (double*) malloc(sizeof(double)*int_v_bufsize);
194 used_storage_build_ += sizeof(double)*int_v_bufsize;
195 if (!int_v0_buf) {
196 ExEnv::errn() << scprintf("couldn't allocate all integral intermediates\n");
197 fail();
198 }
199 add_store(int_v0_buf);
200 int_v_buf = &int_v0_buf[int_v0_bufsize];
201
202 /* Allocate storage for the needed slots. */
203 for (i=0; i<=am12; i++) {
204 for (j=0; j<=am34; j++) {
205 build.int_v_list(i,j,0) = int_v0_buf;
206 int_v0_buf += INT_NCART(i)*INT_NCART(j);
207 for (k=1; k<=am12+am34-i-j; k++) {
208 build.int_v_list(i,j,k) = int_v_buf;
209 int_v_buf += INT_NCART(i)*INT_NCART(j);
210 }
211 }
212 }
213
214
215 /* Allocate storage for the contracted integrals (these are the output
216 * of the build routines). */
217 /* The ci, etc, indices refer to which set of contraction
218 * coefficients we are using. */
219 e0f0_con_int_bufsize = 0;
220 e0f0_con_ints_array = new IntV3Arraydoublep2***[nc1];
221 used_storage_build_ += sizeof(IntV3Arraydoublep2***)*nc1;
222 for (ci=0; ci<nc1; ci++) {
223 e0f0_con_ints_array[ci] = new IntV3Arraydoublep2**[nc2];
224 used_storage_build_ += sizeof(IntV3Arraydoublep2**)*nc2;
225 for (cj=0; cj<nc2; cj++) {
226 e0f0_con_ints_array[ci][cj] = new IntV3Arraydoublep2*[nc3];
227 used_storage_build_ += sizeof(IntV3Arraydoublep2*)*nc3;
228 for (ck=0; ck<nc3; ck++) {
229 e0f0_con_ints_array[ci][cj][ck] = new IntV3Arraydoublep2[nc4];
230 used_storage_build_ += sizeof(IntV3Arraydoublep2)*nc4;
231 for (cl=0; cl<nc4; cl++) {
232 int am12_for_con;
233 int am34_for_con;
234
235 am12_for_con = jmax_for_con[ci] + jmax_for_con[cj] + order;
236 if ((jmax_for_con[ck]!=am3)||(jmax_for_con[cl]!=am4)) {
237 am34_for_con = jmax_for_con[ck] + jmax_for_con[cl] + order;
238 }
239 else {
240 am34_for_con = jmax_for_con[ck] + jmax_for_con[cl];
241 }
242
243#if CHECK_INTEGRAL_ALGORITHM
244 ExEnv::outn() << "am12_for_con: " << am12_for_con << endl;
245 ExEnv::outn() << "am34_for_con: " << am34_for_con << endl;
246#endif
247
248 e0f0_con_ints_array[ci][cj][ck][cl].set_dim(am12_for_con+1,am34_for_con+1);
249 used_storage_build_ += e0f0_con_ints_array[ci][cj][ck][cl].nbyte();
250#if CHECK_INTEGRAL_ALGORITHM
251 ExEnv::outn() << "e0f0_con_ints_array: "
252 << e0f0_con_ints_array[ci][cj][ck][cl].nbyte()
253 << endl;
254#endif
255
256 /* Count how much storage for the integrals is needed. */
257 for (i=0; i<=am12_for_con; i++) {
258 for (k=0; k<=am34_for_con; k++) {
259 int s = INT_NCART(i)
260 * INT_NCART(k);
261 e0f0_con_int_bufsize += s;
262 }
263 }
264 }
265 }
266 }
267 }
268 e0f0_con_int_buf = (double*) malloc(sizeof(double)*e0f0_con_int_bufsize);
269 used_storage_build_ += e0f0_con_int_bufsize * sizeof(double);
270#if CHECK_INTEGRAL_ALGORITHM
271 ExEnv::outn() << "e0f0_int_buf: " << e0f0_con_int_bufsize * sizeof(double) << endl;
272#endif
273 if (!e0f0_con_int_buf) {
274 ExEnv::errn() << scprintf("couldn't allocate contracted integral storage\n");
275 fail();
276 }
277 add_store(e0f0_con_int_buf);
278 /* Allocate storage for the integrals which will be used by the shift
279 * routine. */
280 for (ci=0; ci<nc1; ci++) {
281 for (cj=0; cj<nc2; cj++) {
282 for (ck=0; ck<nc3; ck++) {
283 for (cl=0; cl<nc4; cl++) {
284 int am12_for_con;
285 int am34_for_con;
286
287 am12_for_con = jmax_for_con[ci] + jmax_for_con[cj] + order;
288 if ((jmax_for_con[ck]!=am3)||(jmax_for_con[cl]!=am4)) {
289 am34_for_con = jmax_for_con[ck] + jmax_for_con[cl] + order;
290 }
291 else {
292 am34_for_con = jmax_for_con[ck] + jmax_for_con[cl];
293 }
294
295 for (i=0; i<=am12_for_con; i++) {
296 for (k=0; k<=am34_for_con; k++) {
297 e0f0_con_ints_array[ci][cj][ck][cl](i,k) = 0;
298 }
299 }
300 for (i=0; i<=am12_for_con; i++) {
301 for (k=0; k<=am34_for_con; k++) {
302/* If there are Pople style s=p shells and the shells are ordered
303 * first s and then p and there are no p or d shells on the molecule,
304 * then this algorithm would will allocate a little more storage
305 * than needed. General contraction should be ordered high to
306 * low angular momentum for this reason. */
307 e0f0_con_ints_array[ci][cj][ck][cl](i,k)
308 = e0f0_con_int_buf;
309 e0f0_con_int_buf += INT_NCART(i)
310 * INT_NCART(k);
311 }
312 }
313 }
314 }
315 }
316 }
317
318 /* Initialize the build_routine array. */
319 for (i=0; i<4; i++) {
320 for (j=0; j<4; j++) {
321 for (k=0; k<4; k++) {
322 for (l=0; l<4; l++) {
323 for (m=0; m<2; m++) {
324 build_routine[i][j][k][l][m] = &BuildIntV3::impossible_integral;
325 }
326 }
327 }
328 }
329 }
330
331#define ASSIGN_BUILD(ii,j,k,l) \
332 build_routine[ii][j][k][l][0]= &BuildIntV3::i ## ii ## j ## k ## l ;\
333 build_routine[ii][j][k][l][1]= &BuildIntV3::i ## ii ## j ## k ## l ## eAB;
334#if (MG == 1) || (MG == 2) || (MG == 3) || (MG == 4)
335 ASSIGN_BUILD(0,1,0,0)
336 ASSIGN_BUILD(0,1,0,1)
337 ASSIGN_BUILD(0,1,1,1)
338 ASSIGN_BUILD(1,1,0,0)
339 ASSIGN_BUILD(1,1,1,1)
340#endif
341
342#if (MG == 2) || (MG == 3) || (MG == 4)
343 ASSIGN_BUILD(0,2,0,0)
344 ASSIGN_BUILD(0,2,0,1)
345 ASSIGN_BUILD(0,2,0,2)
346 ASSIGN_BUILD(0,2,1,1)
347 ASSIGN_BUILD(0,2,1,2)
348 ASSIGN_BUILD(0,2,2,2)
349 ASSIGN_BUILD(1,2,0,0)
350 ASSIGN_BUILD(1,2,0,1)
351 ASSIGN_BUILD(1,2,1,1)
352 ASSIGN_BUILD(1,2,1,2)
353 ASSIGN_BUILD(1,2,2,2)
354 ASSIGN_BUILD(2,2,0,0)
355 ASSIGN_BUILD(2,2,0,1)
356 ASSIGN_BUILD(2,2,1,1)
357 ASSIGN_BUILD(2,2,2,2)
358#endif
359
360#if (MG == 3) || (MG == 4)
361 ASSIGN_BUILD(0,3,0,0)
362 ASSIGN_BUILD(0,3,0,1)
363 ASSIGN_BUILD(0,3,0,2)
364 ASSIGN_BUILD(0,3,0,3)
365 ASSIGN_BUILD(0,3,1,1)
366 ASSIGN_BUILD(0,3,1,2)
367 ASSIGN_BUILD(0,3,1,3)
368 ASSIGN_BUILD(0,3,2,2)
369 ASSIGN_BUILD(0,3,2,3)
370 ASSIGN_BUILD(0,3,3,3)
371 ASSIGN_BUILD(1,3,0,0)
372 ASSIGN_BUILD(1,3,0,1)
373 ASSIGN_BUILD(1,3,0,2)
374 ASSIGN_BUILD(1,3,1,1)
375 ASSIGN_BUILD(1,3,1,2)
376 ASSIGN_BUILD(1,3,1,3)
377 ASSIGN_BUILD(1,3,2,2)
378 ASSIGN_BUILD(1,3,2,3)
379 ASSIGN_BUILD(1,3,3,3)
380 ASSIGN_BUILD(2,3,0,0)
381 ASSIGN_BUILD(2,3,0,1)
382 ASSIGN_BUILD(2,3,0,2)
383 ASSIGN_BUILD(2,3,1,1)
384 ASSIGN_BUILD(2,3,1,2)
385 ASSIGN_BUILD(2,3,2,2)
386 ASSIGN_BUILD(2,3,2,3)
387 ASSIGN_BUILD(2,3,3,3)
388 ASSIGN_BUILD(3,3,0,0)
389 ASSIGN_BUILD(3,3,0,1)
390 ASSIGN_BUILD(3,3,0,2)
391 ASSIGN_BUILD(3,3,1,1)
392 ASSIGN_BUILD(3,3,1,2)
393 ASSIGN_BUILD(3,3,2,2)
394 ASSIGN_BUILD(3,3,3,3)
395#endif
396
397 free(jmax_for_con);
398 saved_am12 = am12;
399 saved_am34 = am34;
400 saved_ncon = nc1;
401
402 used_storage_ += used_storage_build_;
403#if CHECK_INTEGRAL_ALGORITHM
404 ExEnv::outn() << "used_storage_build: " << used_storage_build_ << endl;
405#endif
406 }
407
408void
409Int2eV3::int_done_buildgc()
410{
411 int ci,cj,ck;
412
413 used_storage_ -= used_storage_build_;
414 used_storage_build_ = 0;
415
416 free_store();
417
418 for (ci=0; ci<saved_ncon; ci++) {
419 for (cj=0; cj<saved_ncon; cj++) {
420 for (ck=0; ck<saved_ncon; ck++) {
421 delete[] e0f0_con_ints_array[ci][cj][ck];
422 }
423 delete[] e0f0_con_ints_array[ci][cj];
424 }
425 delete[] e0f0_con_ints_array[ci];
426 }
427 delete[] e0f0_con_ints_array;
428
429 }
430
431/* add_store maintains a list of free storage allocated by int_init_buildgc */
432void
433Int2eV3::add_store(void *p)
434{
435 if (!store) {
436 store = (store_list_t*) malloc(sizeof(store_list_t));
437 assert(store);
438 store->p = 0;
439 n_store_last = 0;
440 }
441 if (n_store_last >= STORAGE_CHUNK) {
442 store_list_t* tmp = (store_list_t*) malloc(sizeof(store_list_t));
443 assert(tmp);
444 tmp->p = store;
445 store = tmp;
446 n_store_last = 0;
447 }
448 store->data[n_store_last++] = p;
449 }
450
451/* free_store frees the memory that add_store keeps track of */
452void
453Int2eV3::free_store()
454{
455 _free_store(store,n_store_last);
456 store = 0;
457 }
458
459void
460Int2eV3::_free_store(store_list_t* s, int n)
461{
462 int i;
463 if (!s) return;
464 for (i=0; i<n; i++) {
465 free(s->data[i]);
466 }
467 _free_store(s->p,STORAGE_CHUNK);
468 free(s);
469 }
470
471
472void
473Int2eV3::int_buildgcam(int minam1, int minam2, int minam3, int minam4,
474 int maxam1, int maxam2, int maxam3, int maxam4,
475 int dam1, int dam2, int dam3, int dam4,
476 int sh1, int sh2, int sh3, int sh4,
477 int eAB)
478{
479 int k,m,n;
480 int ci,cj,ck,cl;
481 int maxam12,maxam34;
482 int nc1,nc2,nc3,nc4;
483
484 if (maxam1<0 || maxam2<0 || maxam3<0 || maxam4<0) return;
485 if (minam1<0) minam1=0;
486 if (minam2<0) minam2=0;
487 if (minam3<0) minam3=0;
488 if (minam4<0) minam4=0;
489
490 maxam12 = maxam1 + maxam2;
491 maxam34 = maxam3 + maxam4;
492
493 nc1 = pbs1_->shell(sh1).ncontraction();
494 if (pbs2_.null()) nc2 = 1;
495 else nc2 = pbs2_->shell(sh2).ncontraction();
496 nc3 = pbs3_->shell(sh3).ncontraction();
497 if (pbs4_.null()) nc4 = 1;
498 else nc4 = pbs4_->shell(sh4).ncontraction();
499
500 /* Zero the target contracted integrals that the build routine
501 * will accumulate into. */
502 for (m=minam1; m<=maxam12; m++) {
503 for (n=minam3; n<=maxam34; n++) {
504 int nm_cart = INT_NCART(m)*INT_NCART(n);
505 for (ci=0; ci<nc1; ci++) {
506 if (m < int_shell1->am(ci)+dam1) continue;
507 for (cj=0; cj<nc2; cj++) {
508 if (int_shell1->am(ci)+dam1+int_shell2->am(cj)+dam2 < m)
509 continue;
510 for (ck=0; ck<nc3; ck++) {
511 if (n < int_shell3->am(ck)+dam3) continue;
512 for (cl=0; cl<nc4; cl++) {
513 if (int_shell3->am(ck)+dam3 +int_shell4->am(cl)+dam4 < n)
514 continue;
515 double *tmp = e0f0_con_ints_array[ci][cj][ck][cl](m,n);
516 for (int ii=0; ii<nm_cart; ii++) tmp[ii] = 0.0;
517 }
518 }
519 }
520 }
521 }
522 }
523
524 gen_shell_intermediates(sh1,sh2,sh3,sh4);
525
526 /* If enough of the functions come from generalized contractions
527 * to make it workwhile, then don't do redundant primitives
528 * at the additional cost of slower normalization computations.
529 */
530 if (nc1 + nc2 + nc3 + nc4 > 4)
531 build_using_gcs(nc1,nc2,nc3,nc4,
532 minam1,minam3,maxam12,maxam34,dam1,dam2,dam3,dam4,eAB);
533 else
534 build_not_using_gcs(nc1,nc2,nc3,nc4,
535 minam1,minam3,maxam12,maxam34,dam1,dam2,dam3,dam4,eAB);
536 }
537
538void
539Int2eV3::build_not_using_gcs(int nc1, int nc2, int nc3, int nc4,
540 int minam1, int minam3, int maxam12, int maxam34,
541 int dam1, int dam2, int dam3, int dam4, int eAB)
542{
543 int i,j,k,l,m;
544 int ci,cj,ck,cl;
545 double *bufferprim;
546
547#if 0
548 ExEnv::outn() << scprintf("not_gcs: %d%d%d%d\n",
549 int_expweight1,
550 int_expweight2,
551 int_expweight3,
552 int_expweight4
553 );
554#endif
555
556 /* Sum thru all possible contractions. */
557 for (ci=0; ci<nc1; ci++) {
558 int mlower = int_shell1->am(ci) + dam1;
559 if (mlower < 0) continue;
560 IntV3Arraydoublep2 ***e0f0_i = e0f0_con_ints_array[ci];
561 for (cj=0; cj<nc2; cj++) {
562 int mupper = mlower + int_shell2->am(cj) + dam2;
563 if (mupper < mlower) continue;
564 if (mlower < minam1) mlower = minam1;
565 if (mupper > maxam12) mupper = maxam12;
566 IntV3Arraydoublep2 **e0f0_ij = e0f0_i[cj];
567 for (ck=0; ck<nc3; ck++) {
568 int nlower = int_shell3->am(ck) + dam3;
569 if (nlower < 0) continue;
570 IntV3Arraydoublep2 *e0f0_ijk = e0f0_ij[ck];
571 for (cl=0; cl<nc4; cl++) {
572 int nupper = nlower + int_shell4->am(cl) + dam4;
573 if (nupper < nlower) continue;
574 if (nlower < minam3) nlower = minam3;
575 if (nupper > maxam34) nupper = maxam34;
576
577 /* Loop over the primitives. */
578 for (i=0; i<int_shell1->nprimitive(); i++) {
579 double coef0;
580 coef0 = int_shell1->coefficient_unnorm(ci,i);
581 if (int_expweight1) coef0 = coef0
582 * int_shell1->exponent(i);
583 /* This factor of two comes from the derivative integral formula. */
584 if (int_expweight1) coef0 *= 2.0;
585 if (int_expweight2) coef0 *= 2.0;
586 if (int_expweight3) coef0 *= 2.0;
587 if (int_expweight4) coef0 *= 2.0;
588 if (int_store1) opr1 = int_shell_to_prim[osh1] + i;
589 for (j=0; j<int_shell2->nprimitive(); j++) {
590 double coef1;
591 coef1 = int_shell2->coefficient_unnorm(cj,j);
592 if (int_expweight2) coef1 *= coef0
593 * int_shell2->exponent(j);
594 else coef1 *= coef0;
595 if (int_store1) opr2 = int_shell_to_prim[osh2] + j;
596 for (k=0; k<int_shell3->nprimitive(); k++) {
597 double coef2;
598 coef2 = int_shell3->coefficient_unnorm(ck,k);
599 if (int_expweight3) coef2 *= coef1
600 * int_shell3->exponent(k);
601 else coef2 *= coef1;
602 if (int_store1) opr3 = int_shell_to_prim[osh3] + k;
603 for (l=0; l<int_shell4->nprimitive(); l++) {
604 double coef3;
605 coef3 = int_shell4->coefficient_unnorm(cl,l);
606 if (int_expweight4) coef3 *= coef2
607 * int_shell4->exponent(l);
608 else coef3 *= coef2;
609 if (int_store1) opr4 = int_shell_to_prim[osh4] + l;
610
611 /* Produce the remaining intermediates. */
612 gen_prim_intermediates_with_norm(i,j,k,l, maxam12+maxam34,coef3);
613
614 /* Generate the target integrals. */
615 if ((maxam12 == 0) && (maxam34 == 0)) {
616 /* Do nothing: gen_prim_intermediates has set everything up. */
617 }
618 else if ((minam1<=MG)&&(minam3<=MG)&&(maxam12<=MG)&&(maxam34<=MG)) {
619 if (build_routine[minam1]
620 [maxam12]
621 [minam3]
622 [maxam34][eAB]==&BuildIntV3::impossible_integral){
623 ExEnv::errn() << scprintf("trying to build with int2v%d%d%d%d (exact)\n",
624 minam1,maxam12,minam3,maxam34);
625 }
626 if (!(build.*build_routine[minam1]
627 [maxam12]
628 [minam3]
629 [maxam34][eAB])()) {
630 ExEnv::outn() << "build2e.cc: did not succeed in building all integrals"
631 << endl;
632 abort();
633 }
634 }
635 else {
636 blockbuildprim(minam1,maxam12,minam3,maxam34);
637 }
638
639 /* Contract the primitive target integrals. */
640 /* Throw out all unneeded contractions. */
641 if (i||j||k||l) {
642 for (m=mlower; m<=mupper; m++) {
643 int o;
644 int sizec = contract_length(m,nlower,nupper);
645 double *restrictxx con_ints = e0f0_ijk[cl](m,nlower);
646 bufferprim = build.int_v_list(m,nlower,0);
647
648 for (o=sizec; o!=0; o--) {
649 *con_ints++ += *bufferprim++;
650 }
651
652 }
653 }
654 else {
655 // for the first primitive write to con_ints rather
656 // than accumulate into it
657 for (m=mlower; m<=mupper; m++) {
658 int o;
659 int sizec = contract_length(m,nlower,nupper);
660 double *restrictxx con_ints = e0f0_ijk[cl](m,nlower);
661 bufferprim = build.int_v_list(m,nlower,0);
662
663 for (o=sizec; o!=0; o--) {
664 *con_ints++ = *bufferprim++;
665 }
666
667 }
668 }
669
670 }
671 }
672 }
673 }
674
675 }
676 }
677 }
678 }
679
680 }
681
682void
683Int2eV3::build_using_gcs(int nc1, int nc2, int nc3, int nc4,
684 int minam1, int minam3, int maxam12, int maxam34,
685 int dam1, int dam2, int dam3, int dam4, int eAB)
686{
687 int i,j,k,l,m;
688 int ci,cj,ck,cl;
689 int maxam1234=maxam12+maxam34;
690 double coef0,coef1,coef2,coef3;
691 double ishl1expi=1.0, ishl2expj=1.0, ishl3expk=1.0;
692 double *bufferprim;
693 double c0scale;
694
695 /* Loop over the primitives. */
696 for (i=0; i<int_shell1->nprimitive(); i++) {
697 if (int_store1) opr1 = int_shell_to_prim[osh1] + i;
698 if (int_expweight1) ishl1expi=2.0*int_shell1->exponent(i);
699
700 for (j=0; j<int_shell2->nprimitive(); j++) {
701 if (int_store1) opr2 = int_shell_to_prim[osh2] + j;
702 ishl2expj = (int_expweight2) ?
703 2.0*int_shell2->exponent(j)*ishl1expi : ishl1expi;
704
705 for (k=0; k<int_shell3->nprimitive(); k++) {
706 if (int_store1) opr3 = int_shell_to_prim[osh3] + k;
707 ishl3expk = (int_expweight3) ?
708 2.0*int_shell3->exponent(k)*ishl2expj : ishl2expj;
709
710 for (l=0; l<int_shell4->nprimitive(); l++) {
711 if (int_store1) opr4 = int_shell_to_prim[osh4] + l;
712 c0scale = (int_expweight4) ?
713 2.0*int_shell4->exponent(l)*ishl3expk : ishl3expk;
714
715 /* Produce the remaining intermediates. */
716 gen_prim_intermediates(i,j,k,l, maxam1234);
717
718 /* Generate the target integrals. */
719 if (!maxam1234) {
720 /* Do nothing: gen_prim_intermediates has set everything up. */
721 }
722 else if ((minam1<=MG)&&(minam3<=MG)&&(maxam12<=MG)&&(maxam34<=MG)) {
723 intfunc brptr=build_routine[minam1][maxam12][minam3][maxam34][eAB];
724 if (brptr == &BuildIntV3::impossible_integral) {
725 ExEnv::errn() << scprintf("trying to build with int2v%d%d%d%d (exact)\n",
726 minam1,maxam12,minam3,maxam34);
727 }
728 if (!(build.*brptr)()) {
729 ExEnv::outn() << "build2e.cc: did not succeed in building all integrals"
730 << endl;
731 abort();
732 }
733 }
734 else {
735 blockbuildprim(minam1,maxam12,minam3,maxam34);
736 }
737
738 /* Sum thru all possible contractions.
739 * Throw out all unneeded contractions. */
740
741 for (ci=0; ci<nc1; ci++) {
742 int mlower = int_shell1->am(ci) + dam1;
743 if (mlower < 0) continue;
744 coef0 = int_shell1->coefficient_unnorm(ci,i)*c0scale;
745 IntV3Arraydoublep2 ***e0f0_i = e0f0_con_ints_array[ci];
746 for (cj=0; cj<nc2; cj++) {
747 int mupper = mlower + int_shell2->am(cj) + dam2;
748 if (mupper < mlower) continue;
749 if (mlower < minam1) mlower = minam1;
750 if (mupper > maxam12) mupper = maxam12;
751 coef1 = int_shell2->coefficient_unnorm(cj,j)*coef0;
752 IntV3Arraydoublep2 **e0f0_ij = e0f0_i[cj];
753 for (ck=0; ck<nc3; ck++) {
754 int nlower = int_shell3->am(ck) + dam3;
755 if (nlower < 0) continue;
756 coef2 = int_shell3->coefficient_unnorm(ck,k)*coef1;
757 IntV3Arraydoublep2 *e0f0_ijk = e0f0_ij[ck];
758 for (cl=0; cl<nc4; cl++) {
759 int nupper = nlower + int_shell4->am(cl) + dam4;
760 if (nupper < nlower) continue;
761 if (nlower < minam3) nlower = minam3;
762 if (nupper > maxam34) nupper = maxam34;
763 coef3 = int_shell4->coefficient_unnorm(cl,l)*coef2;
764
765 /* Contract the primitive target integrals. */
766 if (i||j||k||l) {
767 for (m=mlower; m<=mupper; m++) {
768 int o;
769 int sizec = contract_length(m,nlower,nupper);
770 double *restrictxx con_ints = e0f0_ijk[cl](m,nlower);
771 bufferprim = build.int_v_list(m,nlower,0);
772 /* Sum the integrals into the contracted integrals. */
773#ifdef SUNMOS
774 for (o=0; o < sizec; o++) {
775 con_ints[o] += coef3 * bufferprim[o];
776 }
777#else
778 for (o=sizec; o; o--) {
779 *con_ints++ += coef3 * *bufferprim++;
780 }
781#endif
782 }
783
784 }
785 else {
786 for (m=mlower; m<=mupper; m++) {
787 int o;
788 int sizec = contract_length(m,nlower,nupper);
789 double *restrictxx con_ints = e0f0_ijk[cl](m,nlower);
790 bufferprim = build.int_v_list(m,nlower,0);
791 /* Write the integrals to the contracted integrals. */
792#ifdef SUNMOS
793 for (o=0; o < sizec; o++) {
794 con_ints[o] = coef3 * bufferprim[o];
795 }
796#else
797 for (o=sizec; o; o--) {
798 *con_ints++ = coef3 * *bufferprim++;
799 }
800#endif
801 }
802 }
803 }
804 }
805 }
806 }
807
808
809 }
810 }
811 }
812 }
813 }
814
815/* This routine constructs intermediates needed for each quartet of
816 * primitives. It is given the total angular momentum as the argument
817 * and requires that the global primitive offsets and other global
818 * constants be initialized. */
819void
820Int2eV3::gen_prim_intermediates(int pr1, int pr2, int pr3, int pr4, int am)
821{
822 int i;
823 double T;
824 double pmq,pmq2;
825 double AmB,AmB2;
826 /* This is 2^(1/2) * pi^(5/4) */
827 const double sqrt2pi54 = 5.9149671727956129;
828 double conv_to_s;
829
830 if (int_store2) {
831 double *tmp;
832 build.int_v_zeta12 = int_prim_zeta(opr1,opr2);
833 build.int_v_zeta34 = int_prim_zeta(opr3,opr4);
834 build.int_v_oo2zeta12 = int_prim_oo2zeta(opr1,opr2);
835 build.int_v_oo2zeta34 = int_prim_oo2zeta(opr3,opr4);
836 tmp = int_prim_p(opr1,opr2);
837 build.int_v_p120 = *tmp++;
838 build.int_v_p121 = *tmp++;
839 build.int_v_p122 = *tmp;
840 tmp = int_prim_p(opr3,opr4);
841 build.int_v_p340 = *tmp++;
842 build.int_v_p341 = *tmp++;
843 build.int_v_p342 = *tmp;
844 build.int_v_k12 = int_prim_k(opr1,opr2);
845 build.int_v_k34 = int_prim_k(opr3,opr4);
846 }
847 else {
848 build.int_v_zeta12 = int_shell1->exponent(pr1) + int_shell2->exponent(pr2);
849 build.int_v_zeta34 = int_shell3->exponent(pr3) + int_shell4->exponent(pr4);
850 build.int_v_oo2zeta12 = 1.0/build.int_v_zeta12;
851 build.int_v_oo2zeta34 = 1.0/build.int_v_zeta34;
852 build.int_v_p120 = build.int_v_oo2zeta12
853 * ( int_shell1->exponent(pr1) * build.int_v_r10
854 + int_shell2->exponent(pr2) * build.int_v_r20 );
855 build.int_v_p121 = build.int_v_oo2zeta12
856 * ( int_shell1->exponent(pr1) * build.int_v_r11
857 + int_shell2->exponent(pr2) * build.int_v_r21 );
858 build.int_v_p122 = build.int_v_oo2zeta12
859 * ( int_shell1->exponent(pr1) * build.int_v_r12
860 + int_shell2->exponent(pr2) * build.int_v_r22 );
861 build.int_v_p340 = build.int_v_oo2zeta34
862 * ( int_shell3->exponent(pr3) * build.int_v_r30
863 + int_shell4->exponent(pr4) * build.int_v_r40 );
864 build.int_v_p341 = build.int_v_oo2zeta34
865 * ( int_shell3->exponent(pr3) * build.int_v_r31
866 + int_shell4->exponent(pr4) * build.int_v_r41 );
867 build.int_v_p342 = build.int_v_oo2zeta34
868 * ( int_shell3->exponent(pr3) * build.int_v_r32
869 + int_shell4->exponent(pr4) * build.int_v_r42 );
870
871 /* Compute AmB^2 for shell 1 and 2. */
872 AmB = build.int_v_r20 - build.int_v_r10;
873 AmB2 = AmB*AmB;
874 AmB = build.int_v_r21 - build.int_v_r11;
875 AmB2 += AmB*AmB;
876 AmB = build.int_v_r22 - build.int_v_r12;
877 AmB2 += AmB*AmB;
878
879 build.int_v_k12 = sqrt2pi54
880 * build.int_v_oo2zeta12
881 * exp( - int_shell1->exponent(pr1)*int_shell2->exponent(pr2)
882 * build.int_v_oo2zeta12
883 * AmB2 );
884
885 /* Compute AmB^2 for shells 3 and 4. */
886 AmB = build.int_v_r40 - build.int_v_r30;
887 AmB2 = AmB*AmB;
888 AmB = build.int_v_r41 - build.int_v_r31;
889 AmB2 += AmB*AmB;
890 AmB = build.int_v_r42 - build.int_v_r32;
891 AmB2 += AmB*AmB;
892
893 build.int_v_k34 = sqrt2pi54
894 * build.int_v_oo2zeta34
895 * exp( - int_shell3->exponent(pr3)*int_shell4->exponent(pr4)
896 * build.int_v_oo2zeta34
897 * AmB2 );
898
899 build.int_v_oo2zeta12 *= 0.5;
900 build.int_v_oo2zeta34 *= 0.5;
901 }
902
903 build.int_v_ooze = 1.0/(build.int_v_zeta12 + build.int_v_zeta34);
904
905 build.int_v_W0 = build.int_v_ooze*( build.int_v_zeta12 * build.int_v_p120
906 + build.int_v_zeta34 * build.int_v_p340 );
907 build.int_v_W1 = build.int_v_ooze*( build.int_v_zeta12 * build.int_v_p121
908 + build.int_v_zeta34 * build.int_v_p341 );
909 build.int_v_W2 = build.int_v_ooze*( build.int_v_zeta12 * build.int_v_p122
910 + build.int_v_zeta34 * build.int_v_p342 );
911
912 pmq = build.int_v_p120 - build.int_v_p340;
913 pmq2 = pmq*pmq;
914 pmq = build.int_v_p121 - build.int_v_p341;
915 pmq2 += pmq*pmq;
916 pmq = build.int_v_p122 - build.int_v_p342;
917 pmq2 += pmq*pmq;
918
919 T = build.int_v_zeta12
920 * build.int_v_zeta34
921 * build.int_v_ooze * pmq2;
922
923 double *fjttable = fjt_->values(am,T);
924
925 /* Convert the fjttable produced by int_fjt into the S integrals */
926 conv_to_s = sqrt(build.int_v_ooze) * build.int_v_k12 * build.int_v_k34;
927 for (i=0; i<=am; i++) {
928 build.int_v_list(0,0,i)[0] = fjttable[i] * conv_to_s;
929 }
930
931 }
932
933/* This is like gen_prim_intermediates, except the normalization is
934 * put into the ssss integrals. */
935void
936Int2eV3::gen_prim_intermediates_with_norm(int pr1, int pr2, int pr3, int pr4,
937 int am, double norm)
938{
939 int i;
940 double T;
941 double pmq,pmq2;
942 double AmB,AmB2;
943 /* This is 2^(1/2) * pi^(5/4) */
944 const double sqrt2pi54 = 5.9149671727956129;
945 double conv_to_s;
946
947 if (int_store2) {
948 build.int_v_zeta12 = int_prim_zeta(opr1,opr2);
949 build.int_v_zeta34 = int_prim_zeta(opr3,opr4);
950 build.int_v_oo2zeta12 = int_prim_oo2zeta(opr1,opr2);
951 build.int_v_oo2zeta34 = int_prim_oo2zeta(opr3,opr4);
952 build.int_v_p120 = int_prim_p(opr1,opr2,0);
953 build.int_v_p121 = int_prim_p(opr1,opr2,1);
954 build.int_v_p122 = int_prim_p(opr1,opr2,2);
955 build.int_v_p340 = int_prim_p(opr3,opr4,0);
956 build.int_v_p341 = int_prim_p(opr3,opr4,1);
957 build.int_v_p342 = int_prim_p(opr3,opr4,2);
958 build.int_v_k12 = int_prim_k(opr1,opr2);
959 build.int_v_k34 = int_prim_k(opr3,opr4);
960 }
961 else {
962 build.int_v_zeta12 = int_shell1->exponent(pr1) + int_shell2->exponent(pr2);
963 build.int_v_zeta34 = int_shell3->exponent(pr3) + int_shell4->exponent(pr4);
964 build.int_v_oo2zeta12 = 1.0/build.int_v_zeta12;
965 build.int_v_oo2zeta34 = 1.0/build.int_v_zeta34;
966 build.int_v_p120 = build.int_v_oo2zeta12
967 * ( int_shell1->exponent(pr1) * build.int_v_r10
968 + int_shell2->exponent(pr2) * build.int_v_r20 );
969 build.int_v_p121 = build.int_v_oo2zeta12
970 * ( int_shell1->exponent(pr1) * build.int_v_r11
971 + int_shell2->exponent(pr2) * build.int_v_r21 );
972 build.int_v_p122 = build.int_v_oo2zeta12
973 * ( int_shell1->exponent(pr1) * build.int_v_r12
974 + int_shell2->exponent(pr2) * build.int_v_r22 );
975 build.int_v_p340 = build.int_v_oo2zeta34
976 * ( int_shell3->exponent(pr3) * build.int_v_r30
977 + int_shell4->exponent(pr4) * build.int_v_r40 );
978 build.int_v_p341 = build.int_v_oo2zeta34
979 * ( int_shell3->exponent(pr3) * build.int_v_r31
980 + int_shell4->exponent(pr4) * build.int_v_r41 );
981 build.int_v_p342 = build.int_v_oo2zeta34
982 * ( int_shell3->exponent(pr3) * build.int_v_r32
983 + int_shell4->exponent(pr4) * build.int_v_r42 );
984
985 /* Compute AmB^2 for shell 1 and 2. */
986 AmB = build.int_v_r20 - build.int_v_r10;
987 AmB2 = AmB*AmB;
988 AmB = build.int_v_r21 - build.int_v_r11;
989 AmB2 += AmB*AmB;
990 AmB = build.int_v_r22 - build.int_v_r12;
991 AmB2 += AmB*AmB;
992
993 build.int_v_k12 = sqrt2pi54
994 * build.int_v_oo2zeta12
995 * exp( - int_shell1->exponent(pr1)*int_shell2->exponent(pr2)
996 * build.int_v_oo2zeta12
997 * AmB2 );
998
999 /* Compute AmB^2 for shells 3 and 4. */
1000 AmB = build.int_v_r40 - build.int_v_r30;
1001 AmB2 = AmB*AmB;
1002 AmB = build.int_v_r41 - build.int_v_r31;
1003 AmB2 += AmB*AmB;
1004 AmB = build.int_v_r42 - build.int_v_r32;
1005 AmB2 += AmB*AmB;
1006
1007 build.int_v_k34 = sqrt2pi54
1008 * build.int_v_oo2zeta34
1009 * exp( - int_shell3->exponent(pr3)*int_shell4->exponent(pr4)
1010 * build.int_v_oo2zeta34
1011 * AmB2 );
1012
1013 build.int_v_oo2zeta12 *= 0.5;
1014 build.int_v_oo2zeta34 *= 0.5;
1015 }
1016
1017 build.int_v_ooze = 1.0/(build.int_v_zeta12 + build.int_v_zeta34);
1018
1019 build.int_v_W0 = build.int_v_ooze*( build.int_v_zeta12 * build.int_v_p120
1020 + build.int_v_zeta34 * build.int_v_p340 );
1021 build.int_v_W1 = build.int_v_ooze*( build.int_v_zeta12 * build.int_v_p121
1022 + build.int_v_zeta34 * build.int_v_p341 );
1023 build.int_v_W2 = build.int_v_ooze*( build.int_v_zeta12 * build.int_v_p122
1024 + build.int_v_zeta34 * build.int_v_p342 );
1025
1026 pmq = build.int_v_p120 - build.int_v_p340;
1027 pmq2 = pmq*pmq;
1028 pmq = build.int_v_p121 - build.int_v_p341;
1029 pmq2 += pmq*pmq;
1030 pmq = build.int_v_p122 - build.int_v_p342;
1031 pmq2 += pmq*pmq;
1032
1033 T = build.int_v_zeta12
1034 * build.int_v_zeta34
1035 * build.int_v_ooze * pmq2;
1036
1037 double *fjttable = fjt_->values(am,T);
1038
1039 /* Convert the fjttable produced by int_fjt into the S integrals */
1040 conv_to_s = sqrt(build.int_v_ooze)
1041 * build.int_v_k12 * build.int_v_k34 * norm;
1042 for (i=0; i<=am; i++) {
1043 build.int_v_list(0,0,i)[0] = fjttable[i] * conv_to_s;
1044 }
1045
1046 }
1047
1048
1049/* This routine computes the shell intermediates. */
1050void
1051Int2eV3::gen_shell_intermediates(int sh1, int sh2, int sh3, int sh4)
1052{
1053 if (int_store1) {
1054 build.int_v_r10 = int_shell_r(osh1,0);
1055 build.int_v_r11 = int_shell_r(osh1,1);
1056 build.int_v_r12 = int_shell_r(osh1,2);
1057 build.int_v_r20 = int_shell_r(osh2,0);
1058 build.int_v_r21 = int_shell_r(osh2,1);
1059 build.int_v_r22 = int_shell_r(osh2,2);
1060 build.int_v_r30 = int_shell_r(osh3,0);
1061 build.int_v_r31 = int_shell_r(osh3,1);
1062 build.int_v_r32 = int_shell_r(osh3,2);
1063 build.int_v_r40 = int_shell_r(osh4,0);
1064 build.int_v_r41 = int_shell_r(osh4,1);
1065 build.int_v_r42 = int_shell_r(osh4,2);
1066 }
1067 else {
1068 build.int_v_r10 = pbs1_->r(pbs1_->shell_to_center(sh1),0);
1069 build.int_v_r11 = pbs1_->r(pbs1_->shell_to_center(sh1),1);
1070 build.int_v_r12 = pbs1_->r(pbs1_->shell_to_center(sh1),2);
1071 if (pbs2_.null()) {
1072 build.int_v_r20 = 0.0;
1073 build.int_v_r21 = 0.0;
1074 build.int_v_r22 = 0.0;
1075 }
1076 else {
1077 build.int_v_r20 = pbs2_->r(pbs2_->shell_to_center(sh2),0);
1078 build.int_v_r21 = pbs2_->r(pbs2_->shell_to_center(sh2),1);
1079 build.int_v_r22 = pbs2_->r(pbs2_->shell_to_center(sh2),2);
1080 }
1081 build.int_v_r30 = pbs3_->r(pbs3_->shell_to_center(sh3),0);
1082 build.int_v_r31 = pbs3_->r(pbs3_->shell_to_center(sh3),1);
1083 build.int_v_r32 = pbs3_->r(pbs3_->shell_to_center(sh3),2);
1084 if (pbs4_.null()) {
1085 build.int_v_r40 = 0.0;
1086 build.int_v_r41 = 0.0;
1087 build.int_v_r42 = 0.0;
1088 }
1089 else {
1090 build.int_v_r40 = pbs4_->r(pbs4_->shell_to_center(sh4),0);
1091 build.int_v_r41 = pbs4_->r(pbs4_->shell_to_center(sh4),1);
1092 build.int_v_r42 = pbs4_->r(pbs4_->shell_to_center(sh4),2);
1093 }
1094 }
1095 }
1096
1097void
1098Int2eV3::blockbuildprim(int minam1,int maxam12,int minam3,int maxam34)
1099{
1100 int m, b;
1101 int l=maxam12+maxam34;
1102
1103 // the (0,0,m) integrals have already been initialized
1104
1105 // compute (0,b,m) integrals
1106 for (m=l-1; m>=0; m--) {
1107 int bmax = l-m;
1108 if (bmax>maxam34) bmax=maxam34;
1109 blockbuildprim_3(1,bmax,m);
1110 }
1111
1112 // compute (a,b,m) integrals
1113 for (m=maxam12-1; m>=0; m--) {
1114 for (b=0; b<=maxam34; b++) {
1115// This is how the code was for a long while,
1116// but at some point it started giving the wrong
1117// answers and seems wrong from inspection. Valgrind
1118// flags that uninitialized I10i integrals are being
1119// used, which results from amin > 1. I have switched
1120// to the correctly behaving amin = 1.
1121// int amin = minam1-m;
1122// if (amin<1) amin=1;
1123// int amax = maxam12-m;
1124// blockbuildprim_1(amin,amax,b,m);
1125 int amax = maxam12-m;
1126 blockbuildprim_1(1,amax,b,m);
1127 }
1128 }
1129}
1130
1131void
1132Int2eV3::blockbuildprim_1(int amin,int amax,int am34,int m)
1133{
1134 double *I00;
1135 double *I10; /* = [a0|c0](m) */
1136 double *I11; /* = [a0|c0](m+1) */
1137 double *I20; /* = [a-1 0|c0](m) */
1138 double *I21; /* = [a-1 0|c0](m+1) */
1139 double *I31; /* = [a0|c-1 0](m+1) */
1140 int cartindex12;
1141 int cartindex34;
1142 int cartindex1234;
1143 int size34=0,size34m1;
1144 int i12, j12, k12;
1145 int i34, j34, k34;
1146
1147 double ***vlist1;
1148 double **vlist10;
1149 double **vlist11;
1150 double ***vlist2;
1151 double **vlist20;
1152
1153 vlist1 = build.int_v_list(amin-1);
1154 vlist10 = vlist1[am34];
1155
1156 if (am34) {
1157 vlist11 = vlist1[am34-1];
1158 }
1159
1160 if (amin>1) {
1161 vlist2 = build.int_v_list(amin-2);
1162 vlist20 = vlist2[am34];
1163 }
1164
1165 /* The size of the am34 group of primitives. */
1166 size34 = INT_NCART_NN(am34);
1167 /* The size of the group of primitives with ang. mom. = am34 - 1 */
1168 size34m1 = INT_NCART_DEC(am34,size34);
1169
1170 // Some local intermediates
1171 double half_ooze = 0.5 * build.int_v_ooze;
1172 double zeta34_ooze = build.int_v_zeta34 * build.int_v_ooze;
1173 double W0_m_p120 = build.int_v_W0 - build.int_v_p120;
1174 double p120_m_r10 = build.int_v_p120 - build.int_v_r10;
1175 double oo2zeta12 = build.int_v_oo2zeta12;
1176 double p121_m_r11 = build.int_v_p121 - build.int_v_r11;
1177 double W1_m_p121 = build.int_v_W1 - build.int_v_p121;
1178 double p122_m_r12 = build.int_v_p122 - build.int_v_r12;
1179 double W2_m_p122 = build.int_v_W2 - build.int_v_p122;
1180
1181 stack_alignment_check(&half_ooze, "buildprim_1: half_ooze");
1182
1183 for (int am12=amin; am12<=amax; am12++) {
1184 /* Construct the needed intermediate integrals. */
1185 double ***vlist0 = build.int_v_list(am12);
1186 double **vlist00 = vlist0[am34];
1187 I00 = vlist00[m];
1188 I10 = vlist10[m];
1189 I11 = vlist10[m+1];
1190 //I00 = build.int_v_list(am12,am34,m);
1191 //I10 = build.int_v_list(am12-1,am34,m);
1192 //I11 = build.int_v_list(am12-1,am34,m+1);
1193 if (am34) {
1194 I31 = vlist11[m+1];
1195 //I31 = build.int_v_list(am12 - 1, am34 - 1, m + 1);
1196 vlist11 = vlist0[am34-1];
1197 }
1198 if (am12>1) {
1199 I20 = vlist20[m];
1200 I21 = vlist20[m+1];
1201 //I20 = build.int_v_list(am12 - 2, am34, m);
1202 //I21 = build.int_v_list(am12 - 2, am34, m + 1);
1203 }
1204 vlist20 = vlist10;
1205 vlist10 = vlist00;
1206
1207 /* Construct the new integrals. */
1208 cartindex12 = 0;
1209 cartindex1234 = 0;
1210 // the i12==0, k12==0, j12=am12 case (build on y)
1211 i12 = 0;
1212 j12 = am12;
1213 k12 = 0;
1214
1215 int i12y1 = 0; //= INT_CARTINDEX(am12-1,i12,j12-1);
1216 int i12y1s34 = i12y1*size34;
1217 int i12y1s34m1 = i12y1*size34m1;
1218 double *I10i = &I10[i12y1s34];
1219 double *I11i = &I11[i12y1s34];
1220 double *restrictxx I00i = &I00[cartindex1234];
1221 if (j12==1) {
1222 for (cartindex34=0; cartindex34<size34; cartindex34++) {
1223 I00i[cartindex34]
1224 = I10i[cartindex34] * p121_m_r11
1225 + I11i[cartindex34] * W1_m_p121;
1226 }
1227 }
1228 else { // j12 > 1
1229 int i12y2s34 = 0; // = INT_CARTINDEX(am12-2,i12,j12-2)*size34;
1230 double *I20i = &I20[i12y2s34];
1231 double *I21i = &I21[i12y2s34];
1232 for (cartindex34=0; cartindex34<size34; cartindex34++) {
1233 I00i[cartindex34]
1234 = I10i[cartindex34] * p121_m_r11
1235 + I11i[cartindex34] * W1_m_p121
1236 + (j12 - 1) * oo2zeta12 * (I20i[cartindex34]
1237 - I21i[cartindex34] * zeta34_ooze);
1238 }
1239 }
1240 if (am34) {
1241 double *I31i = &I31[i12y1s34m1];
1242 cartindex34 = 0;
1243 for (i34=0; i34<=am34; i34++) {
1244 //note: k34 < am34-i34 instead of <= am34-i34, so j34 > 0
1245 int i34y1 = cartindex34-i34;//=INT_CARTINDEX(am34-1,i34,j34-1)
1246 j34 = am34 - i34;
1247 double j34_half_ooze = j34 * half_ooze;
1248 for (k34=0; k34<am34-i34; k34++) {
1249
1250 I00i[cartindex34] += j34_half_ooze * I31i[i34y1];
1251
1252 j34_half_ooze -= half_ooze;
1253 i34y1++;
1254 /* cartindex34 == INT_CARTINDEX(am34,i34,j34) */
1255 cartindex34++;
1256 }
1257 // increment cartindex34 here since the last k was skipped
1258 cartindex34++;
1259 }
1260 }
1261 cartindex12++;
1262 cartindex1234+=size34;
1263
1264 // the i12==0, j12==am12-1, k12==1 case (build on z)
1265 i12 = 0;
1266 j12 = am12 - 1;
1267 k12 = 1;
1268 int i12z1 = 0;//= INT_CARTINDEX(am12-1,i12,j12);
1269 int i12z1s34 = i12z1*size34;
1270 int i12z1s34m1 = i12z1*size34m1;
1271 I10i = &I10[i12z1s34];
1272 I11i = &I11[i12z1s34];
1273 I00i = &I00[cartindex1234];
1274 for (cartindex34=0; cartindex34<size34; cartindex34++) {
1275 I00i[cartindex34]
1276 = I10i[cartindex34] * p122_m_r12
1277 + I11i[cartindex34] * W2_m_p122;
1278 }
1279 if (am34) {
1280 double *I31i = &I31[i12z1s34m1];
1281 cartindex34 = 0;
1282 for (i34=0; i34<=am34; i34++) {
1283 // skip k34 == 0
1284 cartindex34++;
1285 int i34z1 = cartindex34-i34-1;//=INT_CARTINDEX(am34-1,i34,j34)
1286 double k34_half_ooze = half_ooze;
1287 for (k34=1; k34<=am34-i34; k34++) {
1288 I00i[cartindex34] += k34_half_ooze * I31i[i34z1];
1289 k34_half_ooze += half_ooze;
1290 i34z1++;
1291 cartindex34++;
1292 }
1293 }
1294 }
1295 cartindex12++;
1296 cartindex1234+=size34;
1297 // the i12==0, j12==am12-k12, k12>1 case (build on z)
1298 double k12m1_oo2zeta12 = oo2zeta12;
1299 for (k12=2; k12<=am12-i12; k12++) {
1300 j12 = am12 - k12;
1301 i12z1 = cartindex12-i12-1;//=INT_CARTINDEX(am12-1,i12,j12);
1302 i12z1s34 = i12z1*size34;
1303 i12z1s34m1 = i12z1*size34m1;
1304 int i12z2s34 = (cartindex12-i12-i12-2)*size34;
1305 //=INT_CARTINDEX(am12-2,i12,j12)*size34;
1306 I10i = &I10[i12z1s34];
1307 I11i = &I11[i12z1s34];
1308 double *I20i = &I20[i12z2s34];
1309 double *I21i = &I21[i12z2s34];
1310 I00i = &I00[cartindex1234];
1311 for (cartindex34=0; cartindex34<size34; cartindex34++) {
1312 I00i[cartindex34]
1313 = I10i[cartindex34] * p122_m_r12
1314 + I11i[cartindex34] * W2_m_p122
1315 + k12m1_oo2zeta12 * (I20i[cartindex34]
1316 - I21i[cartindex34] * zeta34_ooze);
1317 }
1318 if (am34) {
1319 double *I31i = &I31[i12z1s34m1];
1320 cartindex34 = 0;
1321 for (i34=0; i34<=am34; i34++) {
1322 // skip k34 == 0
1323 cartindex34++;
1324 int i34z1 = cartindex34-i34-1;//=INT_CARTINDEX(am34-1,i34,j34)
1325 double k34_half_ooze = half_ooze;
1326 for (k34=1; k34<=am34-i34; k34++) {
1327 I00i[cartindex34]
1328 += k34_half_ooze * I31i[i34z1];
1329 k34_half_ooze += half_ooze;
1330 i34z1++;
1331 cartindex34++;
1332 }
1333 }
1334 }
1335 cartindex12++;
1336 cartindex1234+=size34;
1337 k12m1_oo2zeta12 += oo2zeta12;
1338 }
1339
1340 // the i12==1 case (build on x)
1341 i12 = 1;
1342 int i12x1 = cartindex12-am12-1;//=INT_CARTINDEX(am12-1,i12-1,am12-i12)
1343 int i12x1s34 = i12x1*size34;
1344 int i12x1s34m1 = i12x1*size34m1;
1345 I00i = &I00[cartindex1234];
1346 I10i = &I10[i12x1s34];
1347 I11i = &I11[i12x1s34];
1348 //for (k12=0; k12<=am12-i12; k12++)
1349 int k12_cartindex34;
1350 int nk12_size34 = am12*size34;
1351 for (k12_cartindex34=0; k12_cartindex34<nk12_size34; k12_cartindex34++) {
1352 *I00i++ = *I10i++ * p120_m_r10 + *I11i++ * W0_m_p120;
1353 }
1354 I00i = &I00[cartindex1234];
1355 if (am34) {
1356 double *I31i = &I31[i12x1s34m1];
1357 for (k12=0; k12<am12; k12++) {
1358 // skip over i34==0
1359 double *I00is=&I00i[am34+1];
1360 double i34_half_ooze = half_ooze;
1361 for (i34=1; i34<=am34; i34++) {
1362 for (k34=i34; k34<=am34; k34++) { // index_k34 = true_k34 + i34
1363 *I00is++ += i34_half_ooze * *I31i++;
1364 }
1365 i34_half_ooze += half_ooze;
1366 }
1367 I00i += size34;
1368 }
1369 }
1370 cartindex12 += am12;
1371 cartindex1234 += am12*size34;
1372 // the i12>1 case (build on x)
1373 if (am12<2) continue;
1374 double i12m1_oo2zeta12 = oo2zeta12;
1375 i12x1 = cartindex12-am12-1;
1376 i12x1s34 = i12x1*size34;
1377 i12x1s34m1 = i12x1*size34m1;
1378 int i12x2s34 = (cartindex12-am12-am12-1)*size34;
1379 I10i = &I10[i12x1s34];
1380 I11i = &I11[i12x1s34];
1381 double *I20i = &I20[i12x2s34];
1382 double *I21i = &I21[i12x2s34];
1383 I00i = &I00[cartindex1234];
1384 for (i12=2; i12<=am12; i12++) {
1385 int sizek12_size34 = (am12-i12+1)*size34;
1386 int k12_c34;
1387 for (k12_c34=0; k12_c34<sizek12_size34; k12_c34++) {
1388 *I00i++
1389 = *I10i++ * p120_m_r10
1390 + *I11i++ * W0_m_p120
1391 + i12m1_oo2zeta12 * (*I20i++
1392 - *I21i++
1393 * zeta34_ooze);
1394 }
1395 i12m1_oo2zeta12 += oo2zeta12;
1396 }
1397 if (am34) {
1398 double *I31i = &I31[i12x1s34m1];
1399 I00i = &I00[cartindex1234];
1400 for (i12=2; i12<=am12; i12++) {
1401 for (k12=0; k12<=am12-i12; k12++) {
1402 // skip over i34==0
1403 double *I00is=&I00i[am34+1];
1404 double i34_half_ooze = half_ooze;
1405 for (i34=1; i34<=am34; i34++) {
1406 for (k34=0; k34<=am34-i34; k34++) {
1407 *I00is++ += i34_half_ooze * *I31i++;
1408 }
1409 i34_half_ooze += half_ooze;
1410 }
1411 I00i += size34;
1412 }
1413 }
1414 }
1415 }
1416}
1417
1418void
1419Int2eV3::blockbuildprim_3(int bmin,int bmax,int m)
1420{
1421 double *I00;
1422 double *I10; /* = [a0|c0](m) */
1423 double *I11; /* = [a0|c0](m+1) */
1424 double *I20; /* = [a0|c-1 0](m) */
1425 double *I21; /* = [a0|c-1 0](m+1) */
1426 int ci34m1,ci34m2;
1427 int size34,size34m1,size34m2;
1428 int i34, k34;
1429
1430 // These temporaries point to subblocks within the integrals arrays.
1431 double *I10o,*I11o,*I20o,*I21o;
1432
1433 double ***vlist0;
1434 double **vlist01;
1435 double **vlist02;
1436
1437 vlist0 = build.int_v_list(0);
1438 vlist01 = vlist0[bmin-1];
1439 if (bmin>1) {
1440 vlist02 = vlist0[bmin-2];
1441 }
1442
1443 for (int am34=bmin; am34<=bmax; am34++) {
1444
1445 /* Construct the needed intermediate integrals. */
1446 double **vlist00 = vlist0[am34];
1447 I00 = vlist00[m];
1448 I10 = vlist01[m];
1449 I11 = vlist01[m+1];
1450 //I00 = build.int_v_list(0, am34, m);
1451 //I10 = build.int_v_list(0, am34 - 1, m);
1452 //I11 = build.int_v_list(0, am34 - 1, m + 1);
1453 if (am34>1) {
1454 I20 = vlist02[m];
1455 I21 = vlist02[m+1];
1456 //I20 = build.int_v_list(0, am34 - 2, m);
1457 //I21 = build.int_v_list(0, am34 - 2, m + 1);
1458 }
1459 vlist02 = vlist01;
1460 vlist01 = vlist00;
1461
1462 /* The size of the group of primitives with ang. mom. = am34 - 1 */
1463 size34 = INT_NCART_NN(am34);
1464 size34m1 = INT_NCART_DEC(am34,size34);
1465 size34m2 = INT_NCART(am34-2);
1466
1467 // Useful constants
1468 double p340_m_r30 = build.int_v_p340 - build.int_v_r30;
1469 double W0_m_p340 = build.int_v_W0 - build.int_v_p340;
1470 double p341_m_r31 = build.int_v_p341 - build.int_v_r31;
1471 double W1_m_p341 = build.int_v_W1 - build.int_v_p341;
1472 double p342_m_r32 = build.int_v_p342 - build.int_v_r32;
1473 double W2_m_p342 = build.int_v_W2 - build.int_v_p342;
1474 double oo2zeta34 = build.int_v_oo2zeta34;
1475 double zeta12_ooze = build.int_v_zeta12 * build.int_v_ooze;
1476
1477 stack_alignment_check(&p340_m_r30, "buildprim_3: p340_m_r30");
1478
1479 /* Construct the new integrals. */
1480 double *restrictxx I00o = I00; // points the current target integral
1481 I10o = I10;
1482 I11o = I11;
1483 //int cartindex34 = 0;
1484 // i34 == 0, k34 == 0, j34 = am34
1485 /* ------------------ Build from the y position. */
1486 /* I10 I11 and I21 */
1487 *I00o = *I10o * p341_m_r31 + *I11o * W1_m_p341;
1488 if (am34>1) {
1489 I20o = I20;
1490 I21o = I21;
1491 *I00o += (am34 - 1) * oo2zeta34 * (*I20o
1492 - *I21o * zeta12_ooze);
1493 }
1494 //cartindex34++;
1495 // i34 == 0, k34 >= 1
1496 // loop over a portion of the l=am34-1 integrals
1497 I00o = &I00o[1];
1498 for (ci34m1=0; ci34m1<am34; ci34m1++) {
1499 /* ------------------ Build from the z position. */
1500 //note: ci34m1 = cartindex34 - i34 - 1;//=INT_CARTINDEX(am34-1,i34,j34)
1501 /* I10 and I11 */
1502 I00o[ci34m1] = I10o[ci34m1] * p342_m_r32 + I11o[ci34m1] * W2_m_p342;
1503 }
1504 // skip over i34 == 0, k34 == 1
1505 //cartindex34++;
1506 // i34 == 0, k34 > 1
1507 I00o = &I00o[1];
1508 // loop over a portion of the l=am34-2 integrals
1509 double k34m1_oo2zeta34 = oo2zeta34;
1510 for (ci34m2=0; ci34m2<am34-1; ci34m2++) {
1511 //note: k34 = 2+ci34m2
1512 /* ------------------ Build from the z position. */
1513 /* I20 and I21 */
1514 I00o[ci34m2]
1515 += k34m1_oo2zeta34 * (I20o[ci34m2] - I21o[ci34m2] * zeta12_ooze);
1516 k34m1_oo2zeta34 += oo2zeta34;
1517 }
1518 //cartindex34+=am34-1;
1519 // i34 >= 1
1520 I00o = &I00o[am34-1];
1521 //note: ci34m1 = INT_CARTINDEX(am34-1,i34-1,j34)
1522 for (ci34m1=0; ci34m1<size34m1; ci34m1++) {
1523 /* I10 and I11 contrib */
1524 /* ------------------ Build from the x position. */
1525 I00o[ci34m1] = I10o[ci34m1] * p340_m_r30 + I11o[ci34m1] * W0_m_p340;
1526 }
1527 // skip past i34 == 1
1528 //cartindex34 += am34;
1529 // i34 > 1
1530 I00o = &I00o[am34];
1531 //note: ci34m2=INT_CARTINDEX(am34-2,i34-2,j34)
1532 ci34m2=0;
1533 double i34m1_oo2zeta34 = oo2zeta34;
1534 for (i34=2; i34<=am34; i34++) {
1535 for (k34=0; k34<=am34-i34; k34++) {
1536 /* I20 and I21 contrib */
1537 /* ------------------ Build from the x position. */
1538 I00o[ci34m2]
1539 += i34m1_oo2zeta34 * (I20o[ci34m2] - I21o[ci34m2] * zeta12_ooze);
1540 ci34m2++;
1541 }
1542 i34m1_oo2zeta34 += oo2zeta34;
1543 }
1544 //cartindex34 += size34m2;
1545
1546 I00o = &I00o[size34m2];
1547 }
1548}
1549
1550/////////////////////////////////////////////////////////////////////////////
1551
1552// Local Variables:
1553// mode: c++
1554// c-file-style: "CLJ-CONDENSED"
1555// End:
Note: See TracBrowser for help on using the repository browser.