1 //
2 // comp_eri.cc
3 //
4 // Copyright (C) 2001 Edward Valeev
5 //
6 // Author: Edward Valeev <edward.valeev@chemistry.gatech.edu>
7 // Maintainer: EV
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27
28 #include <stdarg.h>
29
30 #include <util/misc/formio.h>
31 #include <chemistry/qc/cints/macros.h>
32 #include <chemistry/qc/cints/eri.h>
33 #include <chemistry/qc/cints/tform.h>
34 #ifdef DMALLOC
35 #include <dmalloc.h>
36 #endif
37
38 using namespace std;
39 using namespace sc;
40
41 static inline void
swtch(GaussianBasisSet * & i,GaussianBasisSet * & j)42 swtch(GaussianBasisSet* &i,GaussianBasisSet* &j)
43 {
44 GaussianBasisSet *tmp;
45 tmp = i;
46 i = j;
47 j = tmp;
48 }
49
50 static inline void
pswtch(void ** i,void ** j)51 pswtch(void**i,void**j)
52 {
53 void*tmp;
54 tmp = *i;
55 *i = *j;
56 *j = tmp;
57 }
58
59 static inline void
iswtch(int * i,int * j)60 iswtch(int *i,int *j)
61 {
62 int tmp;
63 tmp = *i;
64 *i = *j;
65 *j = tmp;
66 }
67
68 static void
fail()69 fail()
70 {
71 ExEnv::errn() << scprintf("failing module:\n%s",__FILE__) << endl;
72 abort();
73 }
74
75 void
compute_quartet(int * psh1,int * psh2,int * psh3,int * psh4)76 EriCints::compute_quartet(int *psh1, int *psh2, int *psh3, int *psh4)
77 {
78 #ifdef EREP_TIMING
79 char section[30];
80 #endif
81 GaussianBasisSet *pbs1=bs1_.pointer();
82 GaussianBasisSet *pbs2=bs2_.pointer();
83 GaussianBasisSet *pbs3=bs3_.pointer();
84 GaussianBasisSet *pbs4=bs4_.pointer();
85 int int_expweight1; // For exponent weighted contractions.
86 int int_expweight2; // For exponent weighted contractions.
87 int int_expweight3; // For exponent weighted contractions.
88 int int_expweight4; // For exponent weighted contractions.
89 int size;
90 int ii;
91 int size1, size2, size3, size4;
92 int tam1,tam2,tam3,tam4;
93 int i,j,k,l;
94 int pi, pj, pk, pl;
95 int gci, gcj, gck, gcl;
96 int sh1,sh2,sh3,sh4;
97 int osh1,osh2,osh3,osh4;
98 int am1,am2,am3,am4,am12,am34;
99 int minam1,minam2,minam3,minam4;
100 int redundant_index;
101 int e12,e13e24,e34;
102 int p12,p34,p13p24;
103 int eAB;
104
105 #ifdef DMALLOC
106 /*--- Test heap before ---*/
107 int
108 heapstate = dmalloc_verify(target_ints_buffer_);
109 if (heapstate == DMALLOC_VERIFY_ERROR)
110 fail();
111 heapstate = dmalloc_verify(cart_ints_);
112 if (heapstate == DMALLOC_VERIFY_ERROR)
113 fail();
114 heapstate = dmalloc_verify(sphharm_ints_);
115 if (heapstate == DMALLOC_VERIFY_ERROR)
116 fail();
117 heapstate = dmalloc_verify(perm_ints_);
118 if (heapstate == DMALLOC_VERIFY_ERROR)
119 fail();
120 heapstate = dmalloc_verify(tformbuf_);
121 if (heapstate == DMALLOC_VERIFY_ERROR)
122 fail();
123 #endif
124
125 osh1 = sh1 = *psh1;
126 osh2 = sh2 = *psh2;
127 osh3 = sh3 = *psh3;
128 osh4 = sh4 = *psh4;
129
130 /* Test the arguments to make sure that they are sensible. */
131 if ( sh1 < 0 || sh1 >= bs1_->nbasis()
132 || sh2 < 0 || sh2 >= bs2_->nbasis()
133 || sh3 < 0 || sh3 >= bs3_->nbasis()
134 || sh4 < 0 || sh4 >= bs4_->nbasis() ) {
135 ExEnv::errn() << scprintf("compute_erep has been incorrectly used\n");
136 ExEnv::errn() << scprintf("shells (bounds): %d (%d), %d (%d), %d (%d), %d (%d)\n",
137 sh1,bs1_->nbasis()-1,
138 sh2,bs2_->nbasis()-1,
139 sh3,bs3_->nbasis()-1,
140 sh4,bs4_->nbasis()-1);
141 fail();
142 }
143
144 /* Set up pointers to the current shells. */
145 int_shell1_ = &bs1_->shell(sh1);
146 int_shell2_ = &bs2_->shell(sh2);
147 int_shell3_ = &bs3_->shell(sh3);
148 int_shell4_ = &bs4_->shell(sh4);
149
150 /* Compute the maximum angular momentum on each centers to
151 * determine the most efficient way to invoke the building and shifting
152 * routines. The minimum angular momentum will be computed at the
153 * same time. */
154 minam1 = int_shell1_->min_am();
155 minam2 = int_shell2_->min_am();
156 minam3 = int_shell3_->min_am();
157 minam4 = int_shell4_->min_am();
158 am1 = int_shell1_->max_am();
159 am2 = int_shell2_->max_am();
160 am3 = int_shell3_->max_am();
161 am4 = int_shell4_->max_am();
162 am12 = am1 + am2;
163 am34 = am3 + am4;
164
165 // This condition being true is guaranteed by the constructor of IntegralCints
166 //if (minam1 != am1 ||
167 // minam2 != am2 ||
168 // minam3 != am3 ||
169 // minam4 != am4 ) {
170 // ExEnv::errn() << scprintf("Int2eCints::comp_eri() cannot yet handle fully general contractions") << endl;
171 // fail();
172 //}
173
174 /* See if need to transform to spherical harmonics */
175 bool need_cart2sph_transform = false;
176 if (int_shell1_->has_pure() ||
177 int_shell2_->has_pure() ||
178 int_shell3_->has_pure() ||
179 int_shell4_->has_pure())
180 need_cart2sph_transform = true;
181
182
183 /* See if contraction quartets need to be resorted into a shell quartet */
184 bool need_sort_to_shell_quartet = false;
185 int num_gen_shells = 0;
186 if (int_shell1_->ncontraction() > 1)
187 num_gen_shells++;
188 if (int_shell2_->ncontraction() > 1)
189 num_gen_shells++;
190 if (int_shell3_->ncontraction() > 1)
191 num_gen_shells++;
192 if (int_shell4_->ncontraction() > 1)
193 num_gen_shells++;
194 if (am12+am34 && num_gen_shells >= 1)
195 need_sort_to_shell_quartet = true;
196
197 /* Unique integrals are needed only ?*/
198 bool need_unique_ints_only = false;
199 if (!redundant_) {
200 e12 = 0;
201 if (int_shell1_ == int_shell2_ && int_shell1_->nfunction()>1)
202 e12 = 1;
203 e34 = 0;
204 if (int_shell3_ == int_shell4_ && int_shell3_->nfunction()>1)
205 e34 = 1;
206 e13e24 = 0;
207 if (int_shell1_ == int_shell3_ && int_shell2_ == int_shell4_ && int_shell1_->nfunction()*int_shell2_->nfunction()>1)
208 e13e24 = 1;
209
210 if ( e12 || e34 || e13e24 )
211 need_unique_ints_only = true;
212 }
213
214
215 #ifdef EREP_TIMING
216 sprintf(section,"erep am=%02d",am12+am34);
217 tim_enter(section);
218 tim_enter("setup");
219 #endif
220
221 /* Convert the integral to the most efficient form. */
222 p12 = 0;
223 p34 = 0;
224 p13p24 = 0;
225
226 if (am2 > am1) {
227 p12 = 1;
228 iswtch(&am1,&am2);iswtch(&sh1,&sh2);iswtch(psh1,psh2);
229 iswtch(&minam1,&minam2);
230 pswtch((void**)&int_shell1_,(void**)&int_shell2_);
231 swtch(pbs1,pbs2);
232 }
233 if (am4 > am3) {
234 p34 = 1;
235 iswtch(&am3,&am4);iswtch(&sh3,&sh4);iswtch(psh3,psh4);
236 iswtch(&minam3,&minam4);
237 pswtch((void**)&int_shell3_,(void**)&int_shell4_);
238 swtch(pbs3,pbs4);
239 }
240 if (am12 > am34) {
241 p13p24 = 1;
242 iswtch(&am1,&am3);iswtch(&sh1,&sh3);iswtch(psh1,psh3);
243 iswtch(&am2,&am4);iswtch(&sh2,&sh4);iswtch(psh2,psh4);
244 iswtch(&am12,&am34);
245 iswtch(&minam1,&minam3);
246 iswtch(&minam2,&minam4);
247 pswtch((void**)&int_shell1_,(void**)&int_shell3_);
248 swtch(pbs1,pbs3);
249 pswtch((void**)&int_shell2_,(void**)&int_shell4_);
250 swtch(pbs2,pbs4);
251 }
252 bool shells_were_permuted = (p12||p34||p13p24);
253
254 /* If the centers were permuted, then the int_expweighted variable may
255 * need to be changed. */
256 if (p12) {
257 iswtch(&int_expweight1,&int_expweight2);
258 }
259 if (p34) {
260 iswtch(&int_expweight3,&int_expweight4);
261 }
262 if (p13p24) {
263 iswtch(&int_expweight1,&int_expweight3);
264 iswtch(&int_expweight2,&int_expweight4);
265 }
266
267 /* Compute the shell sizes. */
268 size1 = int_shell1_->ncartesian();
269 size2 = int_shell2_->ncartesian();
270 size3 = int_shell3_->ncartesian();
271 size4 = int_shell4_->ncartesian();
272 size = size1*size2*size3*size4;
273
274 /* Compute center data for Libint */
275 int ctr1 = pbs1->shell_to_center(sh1);
276 int ctr2 = pbs2->shell_to_center(sh2);
277 int ctr3 = pbs3->shell_to_center(sh3);
278 int ctr4 = pbs4->shell_to_center(sh4);
279 for(i=0;i<3;i++) {
280 Libint_.AB[i] = pbs1->r(ctr1,i) - pbs2->r(ctr2,i);
281 Libint_.CD[i] = pbs3->r(ctr3,i) - pbs4->r(ctr4,i);
282 quartet_info_.A[i] = pbs1->r(ctr1,i);
283 quartet_info_.B[i] = pbs2->r(ctr2,i);
284 quartet_info_.C[i] = pbs3->r(ctr3,i);
285 quartet_info_.D[i] = pbs4->r(ctr4,i);
286 }
287 quartet_info_.AB2 = Libint_.AB[0]*Libint_.AB[0]+Libint_.AB[1]*Libint_.AB[1]+Libint_.AB[2]*Libint_.AB[2];
288 quartet_info_.CD2 = Libint_.CD[0]*Libint_.CD[0]+Libint_.CD[1]*Libint_.CD[1]+Libint_.CD[2]*Libint_.CD[2];
289
290 /* Set up pointers to the current shell pairs. */
291 quartet_info_.shell_pair12 = shell_pairs12_->shell_pair(osh1,osh2);
292 quartet_info_.shell_pair34 = shell_pairs34_->shell_pair(osh3,osh4);
293
294 /* Remember how permuted - will need to access shell pairs in grt_quartet_data_() using the original
295 primitive indices */
296 quartet_info_.p12 = p12;
297 quartet_info_.p34 = p34;
298 quartet_info_.p13p24 = p13p24;
299
300 /* Remember the original primitive indices to access shell pair data
301 Note the reverse order of switching, p13p24 first,
302 then p12 and p34 - because we need the inverse mapping! */
303 quartet_info_.op1 = &quartet_info_.p1;
304 quartet_info_.op2 = &quartet_info_.p2;
305 quartet_info_.op3 = &quartet_info_.p3;
306 quartet_info_.op4 = &quartet_info_.p4;
307 if (p13p24) {
308 pswtch((void **)&quartet_info_.op1,(void **)&quartet_info_.op3);
309 pswtch((void **)&quartet_info_.op2,(void **)&quartet_info_.op4);
310 }
311 if (p12)
312 pswtch((void **)&quartet_info_.op1,(void **)&quartet_info_.op2);
313 if (p34)
314 pswtch((void **)&quartet_info_.op3,(void **)&quartet_info_.op4);
315
316 /* Determine where integrals need to go at each stage */
317 if (shells_were_permuted)
318 if (need_sort_to_shell_quartet) {
319 prim_ints_ = cart_ints_;
320 if (need_cart2sph_transform)
321 contr_quartets_ = sphharm_ints_;
322 else
323 contr_quartets_ = cart_ints_;
324 shell_quartet_ = perm_ints_;
325 }
326 else {
327 prim_ints_ = cart_ints_;
328 if (need_cart2sph_transform) {
329 contr_quartets_ = sphharm_ints_;
330 shell_quartet_ = contr_quartets_;
331 }
332 else
333 shell_quartet_ = cart_ints_;
334 }
335 else
336 if (need_sort_to_shell_quartet) {
337 prim_ints_ = cart_ints_;
338 if (need_cart2sph_transform)
339 contr_quartets_ = sphharm_ints_;
340 else
341 contr_quartets_ = cart_ints_;
342 shell_quartet_ = target_ints_buffer_;
343 }
344 else {
345 if (need_cart2sph_transform) {
346 prim_ints_ = cart_ints_;
347 contr_quartets_ = target_ints_buffer_;
348 shell_quartet_ = target_ints_buffer_;
349 }
350 else {
351 prim_ints_ = target_ints_buffer_;
352 shell_quartet_ = target_ints_buffer_;
353 }
354 }
355
356 /* Begin loops over generalized contractions. */
357 int buffer_offset = 0;
358 for (gci=0; gci<int_shell1_->ncontraction(); gci++) {
359 tam1 = int_shell1_->am(gci);
360 int tsize1 = INT_NCART_NN(tam1);
361 quartet_info_.gc1 = gci;
362 for (gcj=0; gcj<int_shell2_->ncontraction(); gcj++) {
363 tam2 = int_shell2_->am(gcj);
364 int tsize2 = INT_NCART_NN(tam2);
365 quartet_info_.gc2 = gcj;
366 for (gck=0; gck<int_shell3_->ncontraction(); gck++) {
367 tam3 = int_shell3_->am(gck);
368 int tsize3 = INT_NCART_NN(tam3);
369 quartet_info_.gc3 = gck;
370 for (gcl=0; gcl<int_shell4_->ncontraction(); gcl++) {
371 tam4 = int_shell4_->am(gcl);
372 int tsize4 = INT_NCART_NN(tam4);
373 quartet_info_.gc4 = gcl;
374 quartet_info_.am = tam1 + tam2 + tam3 + tam4;
375 int size = tsize1*tsize2*tsize3*tsize4;
376
377 /* Begin loop over primitives. */
378 int num_prim_comb = 0;
379 for (pi=0; pi<int_shell1_->nprimitive(); pi++) {
380 quartet_info_.p1 = pi;
381 for (pj=0; pj<int_shell2_->nprimitive(); pj++) {
382 quartet_info_.p2 = pj;
383 for (pk=0; pk<int_shell3_->nprimitive(); pk++) {
384 quartet_info_.p3 = pk;
385 for (pl=0; pl<int_shell4_->nprimitive(); pl++) {
386 quartet_info_.p4 = pl;
387
388 /* Compute primitive data for Libint */
389 eri_quartet_data_(&(Libint_.PrimQuartet[num_prim_comb++]), 1.0);
390
391 }}}}
392 /* Compute the integrals */
393 if (quartet_info_.am) {
394 REALTYPE *raw_data = build_eri[tam1][tam2][tam3][tam4](&Libint_, num_prim_comb);
395 #ifdef NONDOUBLE_INTS
396 for(int ijkl=0; ijkl<size; ijkl++)
397 prim_ints_[buffer_offset + ijkl] = (double) raw_data[ijkl];
398 #else
399 memcpy(&(prim_ints_[buffer_offset]),raw_data,sizeof(REALTYPE)*size);
400 #endif
401 }
402 else {
403 REALTYPE temp = 0.0;
404 for(int i=0;i<num_prim_comb;i++)
405 temp += Libint_.PrimQuartet[i].F[0];
406 prim_ints_[buffer_offset] = (double) temp;
407 }
408 buffer_offset += size;
409 }}}}
410
411 /*-------------------------------------------
412 Transform to spherical harmonics if needed
413 -------------------------------------------*/
414 if (need_cart2sph_transform)
415 transform_contrquartets_(prim_ints_,contr_quartets_);
416
417 /*----------------------------------------------
418 Resort integrals from by-contraction-quartets
419 into shell-quartet order if needed
420 ----------------------------------------------*/
421 if (need_sort_to_shell_quartet)
422 sort_contrquartets_to_shellquartet_(contr_quartets_,shell_quartet_);
423
424 /*---------------------------------
425 Permute integrals back if needed
426 ---------------------------------*/
427 if ((!permute_)&&shells_were_permuted) {
428 // handle integrals first
429 permute_target_(shell_quartet_,target_ints_buffer_,p13p24,p12,p34);
430 // then indices
431 if (p13p24) {
432 iswtch(&sh1,&sh3);iswtch(psh1,psh3);
433 iswtch(&sh2,&sh4);iswtch(psh2,psh4);
434 iswtch(&am1,&am3);
435 iswtch(&am2,&am4);
436 iswtch(&am12,&am34);
437 pswtch((void**)&int_shell1_,(void**)&int_shell3_);
438 swtch(pbs1,pbs3);
439 pswtch((void**)&int_shell2_,(void**)&int_shell4_);
440 swtch(pbs2,pbs4);
441 iswtch(&int_expweight1,&int_expweight3);
442 iswtch(&int_expweight2,&int_expweight4);
443 }
444 if (p34) {
445 iswtch(&sh3,&sh4);iswtch(psh3,psh4);
446 iswtch(&am3,&am4);
447 pswtch((void**)&int_shell3_,(void**)&int_shell4_);
448 swtch(pbs3,pbs4);
449 iswtch(&int_expweight3,&int_expweight4);
450 }
451 if (p12) {
452 iswtch(&sh1,&sh2);iswtch(psh1,psh2);
453 iswtch(&am1,&am2);
454 pswtch((void**)&int_shell1_,(void**)&int_shell2_);
455 swtch(pbs1,pbs2);
456 iswtch(&int_expweight1,&int_expweight2);
457 }
458 }
459
460 /*--- Extract unique integrals if needed ---*/
461 if (need_unique_ints_only)
462 get_nonredundant_ints_(target_ints_buffer_,target_ints_buffer_,e13e24,e12,e34);
463
464 #ifdef DMALLOC
465 /*--- Test heap after ---*/
466 heapstate = dmalloc_verify(target_ints_buffer_);
467 if (heapstate == DMALLOC_VERIFY_ERROR)
468 fail();
469 heapstate = dmalloc_verify(cart_ints_);
470 if (heapstate == DMALLOC_VERIFY_ERROR)
471 fail();
472 heapstate = dmalloc_verify(sphharm_ints_);
473 if (heapstate == DMALLOC_VERIFY_ERROR)
474 fail();
475 heapstate = dmalloc_verify(perm_ints_);
476 if (heapstate == DMALLOC_VERIFY_ERROR)
477 fail();
478 heapstate = dmalloc_verify(tformbuf_);
479 if (heapstate == DMALLOC_VERIFY_ERROR)
480 fail();
481 #endif
482
483 return;
484 }
485
486
487 /////////////////////////////////////////////////////////////////////////////
488
489 // Local Variables:
490 // mode: c++
491 // c-file-style: "CLJ-CONDENSED"
492 // End:
493