source: ThirdParty/mpqc_open/src/lib/chemistry/qc/intv3/shift2e.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: 15.5 KB
Line 
1//
2// shift2e.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 <util/misc/formio.h>
29#include <chemistry/qc/intv3/macros.h>
30#include <chemistry/qc/intv3/int2e.h>
31
32using namespace std;
33using namespace sc;
34
35//#undef CHECK_INTEGRAL_ALGORITHM
36//#define CHECK_INTEGRAL_ALGORITHM 1
37
38static inline void
39iswtch(int *i, int *j)
40{
41 int tmp;
42
43 tmp = *i;
44 *i = *j;
45 *j = tmp;
46}
47
48/* This initializes the shift routines. It is called by int_initialize_erep.
49 * It is passed the maximum am to be found on each center.
50 */
51void
52Int2eV3::int_init_shiftgc(int order, int am1, int am2, int am3, int am4)
53{
54 /* The intermediate integral arrays are allocated by the
55 * build initialization routine. */
56
57 used_storage_shift_ = 0;
58
59 /* Convert the am1-4 to their canonical ordering. */
60 if (am2>am1) {
61 iswtch(&am1,&am2);
62 }
63 if (am4>am3) {
64 iswtch(&am3,&am4);
65 }
66 if ((am3 > am1)||((am3 == am1)&&(am4 > am2))) {
67 iswtch(&am1,&am3);
68 iswtch(&am2,&am4);
69 }
70
71 /* If the center permutation 1<->3 and 2<->4 is performed, then
72 * we may need the am for center 2 to be as big as for center 4. */
73 if (am4 > am2) am2 = am4;
74
75 /* If derivatives are needed am1 will need to be larger. */
76 if (order==1) am1++;
77 /* For derivative integral bounds am3 will need to be larger. */
78 if (order==1 && int_derivative_bounds) am3++;
79
80 // Set up the new intermediate arrays.
81 int e, c, d;
82 int ndata34_e = 0;
83 for (e=am1; e<=am1+am2; e++) {
84 int size_e = INT_NCART(e);
85 ndata34_e += size_e;
86 }
87 int ndata34_f = 0;
88 for (d=1; d<=am4; d++) {
89 int size_d = INT_NCART(d);
90 int size_dm1 = INT_NCART(d-1);
91 int off_cp1_dm1 = INT_NCART(am3) * size_dm1;
92 int off_c_d = 0;
93 for (c=am3; c<=am3+am4-d; c++) {
94 int size_c = INT_NCART(c);
95 int size_cp1 = INT_NCART(c+1);
96 off_c_d += size_c * size_d;
97 off_cp1_dm1 += size_cp1 * size_dm1;
98 }
99 if (off_c_d > ndata34_f) ndata34_f = off_c_d;
100 if (off_cp1_dm1 > ndata34_f) ndata34_f = off_cp1_dm1;
101 }
102 int ndata34 = ndata34_e * ndata34_f;
103
104 int ndata12 = 0;
105 int a, b;
106 int size_c_d = INT_NCART(am3)*INT_NCART(am4);
107 for (b=1; b<=am2; b++) {
108 int size_b = INT_NCART(b);
109 int size_bm1 = INT_NCART(b-1);
110 int off_a_b = 0;
111 int off_ap1_bm1 = INT_NCART(am1) * size_bm1 * size_c_d;
112 for (a=am1; a<=am1+am2-b; a++) {
113 int size_a = INT_NCART(a);
114 int size_ap1 = INT_NCART(a+1);
115 off_a_b += size_a * size_b * size_c_d;
116 off_ap1_bm1 += size_ap1 * size_bm1 * size_c_d;
117 }
118 if (off_a_b > ndata12) ndata12 = off_a_b;
119 if (off_ap1_bm1 > ndata12) ndata12 = off_ap1_bm1;
120 }
121 int ndatamax = (ndata12>ndata34?ndata12:ndata34);
122 buf34 = new double[ndata34];
123 buf12 = new double[ndata12];
124 bufshared = new double[ndatamax];
125
126 used_storage_shift_ += sizeof(double)*(ndata34+ndata12+ndatamax);
127
128 used_storage_ += used_storage_shift_;
129 }
130
131void
132Int2eV3::int_done_shiftgc()
133{
134 used_storage_ -= used_storage_shift_;
135 delete[] buf12;
136 delete[] buf34;
137 delete[] bufshared;
138 }
139
140/* This is the principle entry point for the am shifting routines.
141 * tam1-4 is the target angular momentum on centers 1-4
142 * sh1-4 are the shell numbers on centers 1-4
143 */
144double *
145Int2eV3::int_shiftgcam(int gc1, int gc2, int gc3, int gc4,
146 int tam1, int tam2, int tam3, int tam4, int peAB)
147{
148 int am1,am2,am3,am4;
149
150 /* Copy the gc{1,2,3,4} into g{1,2,3,4} (static globals). */
151 g1 = gc1;
152 g2 = gc2;
153 g3 = gc3;
154 g4 = gc4;
155
156 /* Compute the angular momentum quartet. */
157 am1 = tam1;
158 am2 = tam2;
159 am3 = tam3;
160 am4 = tam4;
161
162 // (a0|b0) does need shifting
163 if (am2==0 && am4==0) {
164 return e0f0_con_ints_array[g1][g2][g3][g4](am1,am3);
165 }
166
167 /* Copy the A B equivalency info into a static global variable. */
168 eAB = peAB;
169
170 /* Compute the intermediates. */
171 AmB[0] = build.int_v_r10 - build.int_v_r20;
172 AmB[1] = build.int_v_r11 - build.int_v_r21;
173 AmB[2] = build.int_v_r12 - build.int_v_r22;
174 CmD[0] = build.int_v_r30 - build.int_v_r40;
175 CmD[1] = build.int_v_r31 - build.int_v_r41;
176 CmD[2] = build.int_v_r32 - build.int_v_r42;
177
178#if CHECK_INTEGRAL_ALGORITHM > 1
179 ExEnv::outn() << "generating ("
180 << am1 << "," << am2 << "," << am3 << "," << am4 << ")"
181 << ":" << endl;
182#endif
183
184 // the (e0|f0) integrals have been initialized
185 IntV3Arraydoublep2 &e0f0 = e0f0_con_ints_array[g1][g2][g3][g4];
186
187 // generate (e0|cd) for each needed e
188 int e, c, d;
189 int off_e = 0;
190 int size34 = INT_NCART(am3)*INT_NCART(am4);
191 double *buf34_1 = buf34;
192 double *buf34_2 = bufshared;
193 for (e=am1; e<=am1+am2; e++) {
194 int size_e = INT_NCART(e);
195 for (d=1; d<=am4; d++) {
196 int size_d = INT_NCART(d);
197 int size_dm1 = INT_NCART(d-1);
198 int off_c_dm1 = 0;
199 int off_cp1_dm1 = size_e * INT_NCART(am3) * size_dm1;
200 int off_c_d = 0;
201 for (c=am3; c<=am3+am4-d; c++) {
202 int size_c = INT_NCART(c);
203 int size_cp1 = INT_NCART(c+1);
204 double *I0001, *I0010, *I0000;
205 if (d==am4) {
206 I0001 = &buf12[off_e];
207 }
208 else I0001 = &buf34_1[off_c_d];
209 if (d==1) {
210 I0010 = e0f0(e,c+1);
211 I0000 = e0f0(e,c);
212 }
213 else {
214 I0010 = &buf34_2[off_cp1_dm1];
215 I0000 = &buf34_2[off_c_dm1];
216 }
217 shiftam_34(I0001,I0010,I0000,e,0,c,d);
218 off_c_d += size_e * size_c * size_d;
219 off_c_dm1 = off_cp1_dm1;
220 off_cp1_dm1 += size_e * size_cp1 * size_dm1;
221 }
222 // swap the buffers.
223 double *tmp = buf34_1;
224 buf34_1 = buf34_2;
225 buf34_2 = tmp;
226 }
227 off_e += size_e * size34;
228 }
229
230 // generate (ab|cd)
231 int a, b;
232 int size_c_d = size34;
233 double *buf12_1 = bufshared;
234 double *buf12_2 = buf12;
235 for (b=1; b<=am2; b++) {
236 int size_b = INT_NCART(b);
237 int size_bm1 = INT_NCART(b-1);
238 int off_a_b = 0;
239 int off_ap1_bm1 = INT_NCART(am1) * size_bm1 * size_c_d;
240 int off_a_bm1 = 0;
241 for (a=am1; a<=am1+am2-b; a++) {
242 int size_a = INT_NCART(a);
243 int size_ap1 = INT_NCART(a+1);
244 double *I0100 = &buf12_1[off_a_b];
245 double *I1000;
246 double *I0000;
247 if (b==1 && am4 == 0) {
248 I1000 = e0f0(a+1,am3);
249 if (eAB) I0000 = 0;
250 else I0000 = e0f0(a,am3);
251 }
252 else {
253 I1000 = &buf12_2[off_ap1_bm1];
254 if (eAB) I0000 = 0;
255 else I0000 = &buf12_2[off_a_bm1];
256 }
257 if (eAB) shiftam_12eAB(I0100,I1000,I0000,a,b,am3,am4);
258 else shiftam_12(I0100,I1000,I0000,a,b,am3,am4);
259 off_a_b += size_a * size_b * size_c_d;
260 off_a_bm1 = off_ap1_bm1;
261 off_ap1_bm1 += size_ap1 * size_bm1 * size_c_d;
262 }
263 // swap the buffers.
264 double *tmp = buf12_1;
265 buf12_1 = buf12_2;
266 buf12_2 = tmp;
267 }
268
269 /* Construct the target integrals. */
270 return buf12_2;
271 }
272
273/* Shift angular momentum from center 1 to center 2.
274 * I0100 are the target integrals.
275 * am1-4 is the angular momentum on each of the centers in the target set.
276 */
277void
278Int2eV3::shiftam_12(double *I0100, double *I1000, double *I0000,
279 int am1, int am2, int am3, int am4)
280{
281 int i;
282 int i1,k1;
283 int size2, size2m134, size34;
284
285#if CHECK_INTEGRAL_ALGORITHM > 1
286 ExEnv::outn() << "(" << am1 << "," << am2 << "," << am3 << "," << am4 << ")"
287 << " <- "
288 << "(" << am1+1 << "," << am2-1 << "," << am3 << "," << am4 << ")"
289 << "(" << am1 << "," << am2-1 << "," << am3 << "," << am4 << ")"
290 << endl;
291#endif
292
293 size2m134 = INT_NCART(am2-1)*INT_NCART(am3)*INT_NCART(am4);
294 size34 = INT_NCART(am3)*INT_NCART(am4);
295 size2 = INT_NCART(am2);
296
297 int size_zcontrib = am2*size34;
298 int size_xcontrib = (size2-(am2+1))*size34;
299
300 double AmB0 = AmB[0];
301 double AmB1 = AmB[1];
302 double AmB2 = AmB[2];
303
304 /* Loop over the target integrals. */
305 double *restrictxx I0100i=I0100;
306 int cartindex1 = 0;
307 for (i1=0; i1<=am1; i1++) {
308 for (k1=0; k1<=am1-i1; k1++) {
309 //int j1 = am1 - i1 - k1;
310 int ci1x1 = (cartindex1 + am1 + 2) * size2m134;
311 int ci1y1 = (cartindex1 + i1) * size2m134;
312 int ci1z1 = (cartindex1 + i1 + 1) * size2m134;
313 //note:
314 //ci1x1 = INT_CARTINDEX(am1+1,i1+1,j1) * size2m134;
315 //ci1y1 = INT_CARTINDEX(am1+1,i1,j1+1) * size2m134;
316 //ci1z1 = INT_CARTINDEX(am1+1,i1,j1) * size2m134;
317 int ci1 = cartindex1 * size2m134;
318 // i2 == 0, k2 == 0, j2 == am2 (>0)
319 double *I1000i=&I1000[ci1y1];
320 double *I0000i=&I0000[ci1];
321 for (i=0; i<size34; i++) {
322 I0100i[i] = I1000i[i] + I0000i[i] * AmB1;
323 }
324 I0100i=&I0100i[size34];
325 // i2 == 0, k2 > 0
326 I1000i=&I1000[ci1z1];
327 I0000i=&I0000[ci1];
328 for (i=0; i<size_zcontrib; i++) {
329 I0100i[i] = I1000i[i] + I0000i[i] * AmB2;
330 }
331 I0100i=&I0100i[size_zcontrib];
332 // i2 >= 1
333 I1000i=&I1000[ci1x1];
334 I0000i=&I0000[ci1];
335 for (i=0; i<size_xcontrib; i++) {
336 I0100i[i] = I1000i[i] + I0000i[i] * AmB0;
337 }
338 I0100i=&I0100i[size_xcontrib];
339
340 cartindex1++;
341 }
342 }
343 }
344
345
346/* Shift angular momentum from center 1 to center 2 when centers
347 * one and two are the same.
348 * I0100 are the target integrals.
349 * am1-4 is the angular momentum on each of the centers in the target set.
350 */
351void
352Int2eV3::shiftam_12eAB(double *I0100, double *I1000, double *I0000,
353 int am1, int am2, int am3, int am4)
354{
355 int i;
356 int i1,k1;
357 int size2, size2m134, size34;
358
359#if CHECK_INTEGRAL_ALGORITHM > 1
360 ExEnv::outn() << "(" << am1 << "," << am2 << "," << am3 << "," << am4 << ")"
361 << " <- "
362 << "(" << am1+1 << "," << am2-1 << "," << am3 << "," << am4 << ")"
363 << "(" << am1 << "," << am2-1 << "," << am3 << "," << am4 << ")"
364 << endl;
365#endif
366
367 size2m134 = INT_NCART(am2-1)*INT_NCART(am3)*INT_NCART(am4);
368 size34 = INT_NCART(am3)*INT_NCART(am4);
369 size2 = INT_NCART(am2);
370
371 int size_zcontrib = am2*size34;
372 int size_xcontrib = (size2-(am2+1))*size34;
373
374 /* Loop over the target integrals. */
375 double *restrictxx I0100i=I0100;
376 int cartindex1 = 0;
377 for (i1=0; i1<=am1; i1++) {
378 for (k1=0; k1<=am1-i1; k1++) {
379 //int j1 = am1 - i1 - k1;
380 int ci1x1 = (cartindex1 + am1 + 2) * size2m134;
381 int ci1y1 = (cartindex1 + i1) * size2m134;
382 int ci1z1 = (cartindex1 + i1 + 1) * size2m134;
383 //note:
384 //ci1x1 = INT_CARTINDEX(am1+1,i1+1,j1) * size2m134;
385 //ci1y1 = INT_CARTINDEX(am1+1,i1,j1+1) * size2m134;
386 //ci1z1 = INT_CARTINDEX(am1+1,i1,j1) * size2m134;
387 // i2 == 0, k2 == 0, j2 == am2 (>0)
388 double *I1000i=&I1000[ci1y1];
389 for (i=0; i<size34; i++) {
390 I0100i[i] = I1000i[i];
391 }
392 I0100i=&I0100i[size34];
393 // i2 == 0, k2 > 0
394 I1000i=&I1000[ci1z1];
395 for (i=0; i<size_zcontrib; i++) {
396 I0100i[i] = I1000i[i];
397 }
398 I0100i=&I0100i[size_zcontrib];
399 // i2 >= 1
400 I1000i=&I1000[ci1x1];
401 for (i=0; i<size_xcontrib; i++) {
402 I0100i[i] = I1000i[i];
403 }
404 I0100i=&I0100i[size_xcontrib];
405
406 cartindex1++;
407 }
408 }
409 }
410
411void
412Int2eV3::shiftam_34(double *restrictxx I0001, double *I0010, double *I0000,
413 int am1, int am2, int am3, int am4)
414{
415 int i1,k1,cartindex1;
416 int i2,k2,cartindex2;
417 int i3,k3,cartindex3;
418 int i4,k4,cartindex4;
419 int cartindex1234;
420 int size23p14m1,size3p14m1,size4m1,size234m1,size34m1;
421
422#if CHECK_INTEGRAL_ALGORITHM > 1
423 ExEnv::outn() << "(" << am1 << "," << am2 << "," << am3 << "," << am4 << ")"
424 << " <- "
425 << "(" << am1 << "," << am2 << "," << am3+1 << "," << am4-1 << ")"
426 << "(" << am1 << "," << am2 << "," << am3 << "," << am4-1 << ")"
427 << endl;
428#endif
429
430 size23p14m1 = INT_NCART(am2)*INT_NCART(am3+1)*INT_NCART(am4-1);
431 size3p14m1 = INT_NCART(am3+1)*INT_NCART(am4-1);
432 size4m1 = INT_NCART(am4-1);
433
434 size234m1 = INT_NCART(am2)*INT_NCART(am3)*INT_NCART(am4-1);
435 size34m1 = INT_NCART(am3)*INT_NCART(am4-1);
436
437 double CmD0 = CmD[0];
438 double CmD1 = CmD[1];
439 double CmD2 = CmD[2];
440
441 /* Loop over the target integrals. */
442 cartindex1234 = 0;
443 cartindex1 = 0;
444 for (i1=0; i1<=am1; i1++) {
445 for (k1=0; k1<=am1-i1; k1++) {
446 //int j1 = am1 - i1 - k1;
447 int ci1_I0010 = cartindex1 * size23p14m1;
448 int ci1_I0000 = cartindex1 * size234m1;
449 cartindex2 = 0;
450 for (i2=0; i2<=am2; i2++) {
451 for (k2=0; k2<=am2-i2; k2++) {
452 //int j2 = am2 - i2 - k2;
453 int ci2_I0010 = ci1_I0010 + cartindex2 * size3p14m1;
454 int ci2_I0000 = ci1_I0000 + cartindex2 * size34m1;
455 cartindex3 = 0;
456 for (i3=0; i3<=am3; i3++) {
457 for (k3=0; k3<=am3-i3; k3++) {
458 //int j3 = am3 - i3 - k3;
459 //note: cartindex3 + am3 + 2 = INT_CARTINDEX(am3+1,i3+1,j3)
460 int ci3_I0010 = ci2_I0010
461 + (cartindex3 + am3 + 2)*size4m1;
462 int ci3_I0000 = ci2_I0000 + cartindex3*size4m1;
463 //cartindex4 = 0;
464 // this routine called only when am4 > 0
465 ///// CASE 1: i4 = 0 k4 = 0 j4 = am4; shift on y
466 //note: j4 = am4;
467 //note: cartindex4 - i4 = INT_CARTINDEX(am4-1,i4,j4-1)
468 //note: cartindex3 - i3 = INT_CARTINDEX(am3+1,i3,j3+1)
469 int ci3 = cartindex3 + i3;
470 I0001[cartindex1234]
471 = I0010[ci2_I0010 + ci3 * size4m1]
472 + I0000[ci3_I0000]
473 * CmD1;
474 cartindex1234++;
475 //cartindex4++;
476 ///// CASE 2: i4 = 0 k4 > 0; shift on z
477 ci3++;
478 for (int ci4=0; ci4<am4; ci4++) {
479 //note: j4 = am4 - i4 - k4;
480 //note: cartindex4 - i4 - 1 = INT_CARTINDEX(am4-1,i4,j4)
481 //note: ci4 = cartindex4 - i4 - 1;
482 //note: cartindex3 - i3 - 1 = INT_CARTINDEX(am3+1,i3,j3)
483 I0001[cartindex1234]
484 = I0010[ci2_I0010 + ci3 * size4m1 + ci4 ]
485 + I0000[ci3_I0000 + ci4 ]
486 * CmD2;
487 cartindex1234++;
488 //cartindex4++;
489 }
490 ///// CASE 3: i4 > 0; shift on x
491 int ncart_remain = INT_NCART(am4) - (am4+1);
492 for (int ci4=0; ci4<ncart_remain; ci4++) {
493 //note: j4 = am4 - i4 - k4;
494 //note: cartindex4 - am4 - 1 = INT_CARTINDEX(am4-1,i4-1,j4)
495 //note: ci4 = cartindex4 - am4 - 1;
496 I0001[cartindex1234]
497 = I0010[ci3_I0010 + ci4]
498 + I0000[ci3_I0000 + ci4]
499 * CmD0;
500 cartindex1234++;
501 //cartindex4++;
502 }
503 cartindex3++;
504 }
505 }
506 cartindex2++;
507 }
508 }
509 cartindex1++;
510 }
511 }
512 }
513
514
515/////////////////////////////////////////////////////////////////////////////
516
517// Local Variables:
518// mode: c++
519// c-file-style: "CLJ-CONDENSED"
520// End:
Note: See TracBrowser for help on using the repository browser.