1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
30 /** @file integrals_2el_boxed.cc
31
32 @brief Code for 2-electron integrals, computation of Coulomb (J)
33 and HF exchange (K) matrices using a single box.
34
35 @author: Elias Rudberg <em>responsible</em>
36 */
37
38 #include <string.h>
39
40 #include "integrals_2el_boxed.h"
41 #include "integrals_2el_utils.h"
42 #include "organize_distrs.h"
43 #include "pi.h"
44 #include "utilities.h"
45
46
47 static const int HUGE_INTEGER_NUMBER = 2000000000;
48
49
50
51 typedef struct
52 {
53 int a, b, c, d;
54 int poly_ab_index;
55 int poly_cd_index;
56 int idx1;
57 int idx2;
58 ergo_real densValue;
59 } abcd_struct;
60
61 #define set_abcd_list_item_macro(i,A,B,C,D,v,i1,i2) \
62 list[i].a = A; list[i].b = B; list[i].c = C; list[i].d = D; list[i].densValue = v; list[i].idx1 = i1; list[i].idx2 = i2;
63
64
65
66 static int
get_JK_contribs_from_2_interacting_boxes(const BasisInfoStruct & basisInfo,const IntegralInfo & integralInfo,int maxNoOfMonomials,ergo_real * J,ergo_real * K,const ergo_real * dens,const minimal_distr_struct * minimalDistrList_1,int noOfGroups_1,const distr_group_struct * groupList_1,const minimal_distr_struct * minimalDistrList_2,int noOfGroups_2,const distr_group_struct * groupList_2,const cluster_struct * clusterList_1,int nClusters_1,const cluster_struct * clusterList_2,int nClusters_2,const batch_struct * batchList_1,int nBatchs_1,const batch_struct * batchList_2,int nBatchs_2,const basis_func_pair_struct * basisFuncPairList_1,const basis_func_pair_struct * basisFuncPairList_2,int interactionWithSelf,ergo_real threshold,JK_contribs_buffer_struct * bufferStructPtr)67 get_JK_contribs_from_2_interacting_boxes(const BasisInfoStruct & basisInfo,
68 const IntegralInfo & integralInfo,
69 int maxNoOfMonomials,
70 ergo_real* J,
71 ergo_real* K,
72 const ergo_real* dens,
73 const minimal_distr_struct* minimalDistrList_1,
74 int noOfGroups_1,
75 const distr_group_struct* groupList_1,
76 const minimal_distr_struct* minimalDistrList_2,
77 int noOfGroups_2,
78 const distr_group_struct* groupList_2,
79 const cluster_struct* clusterList_1,
80 int nClusters_1,
81 const cluster_struct* clusterList_2,
82 int nClusters_2,
83 const batch_struct* batchList_1,
84 int nBatchs_1,
85 const batch_struct* batchList_2,
86 int nBatchs_2,
87 const basis_func_pair_struct* basisFuncPairList_1,
88 const basis_func_pair_struct* basisFuncPairList_2,
89 int interactionWithSelf,
90 ergo_real threshold,
91 JK_contribs_buffer_struct* bufferStructPtr)
92 {
93 int n = basisInfo.noOfBasisFuncs;
94
95 const JK::ExchWeights CAM_params_not_used;
96
97 const ergo_real twoTimesPiToPow5half = 2 * pitopow52;// = 2 * pow(pi, 2.5);
98 ergo_real* summedIntegralList = bufferStructPtr->summedIntegralList;
99 ergo_real* primitiveIntegralList = bufferStructPtr->primitiveIntegralList;
100 ergo_real* primitiveIntegralList_work = bufferStructPtr->primitiveIntegralList_work;
101
102 for(int batch_i = 0; batch_i < nBatchs_1; batch_i++)
103 {
104 int batch_j_start = 0;
105 if(interactionWithSelf == 1)
106 batch_j_start = batch_i;
107 for(int batch_j = batch_j_start; batch_j < nBatchs_2; batch_j++)
108 {
109 int noOfBasisFuncPairs_1 = batchList_1[batch_i].noOfBasisFuncPairs;
110 int noOfBasisFuncPairs_2 = batchList_2[batch_j].noOfBasisFuncPairs;
111 // set integral list to zero
112 memset(summedIntegralList, 0, noOfBasisFuncPairs_1*noOfBasisFuncPairs_2*sizeof(ergo_real));
113
114 // get largest dmat element
115 ergo_real maxabsdmatelement = 0;
116 for(int i = 0; i < noOfBasisFuncPairs_1; i++)
117 for(int j = 0; j < noOfBasisFuncPairs_2; j++)
118 {
119 int a = basisFuncPairList_1[batchList_1[batch_i].basisFuncPairListIndex+i].index_1;
120 int b = basisFuncPairList_1[batchList_1[batch_i].basisFuncPairListIndex+i].index_2;
121 int c = basisFuncPairList_2[batchList_2[batch_j].basisFuncPairListIndex+j].index_1;
122 int d = basisFuncPairList_2[batchList_2[batch_j].basisFuncPairListIndex+j].index_2;
123 ergo_real absval;
124 if(J != NULL)
125 {
126 absval = template_blas_fabs(dens[a*n+b]);
127 if(absval > maxabsdmatelement)
128 maxabsdmatelement = absval;
129 absval = template_blas_fabs(dens[c*n+d]);
130 if(absval > maxabsdmatelement)
131 maxabsdmatelement = absval;
132 }
133 if(K != NULL)
134 {
135 absval = template_blas_fabs(dens[a*n+c]);
136 if(absval > maxabsdmatelement)
137 maxabsdmatelement = absval;
138 absval = template_blas_fabs(dens[a*n+d]);
139 if(absval > maxabsdmatelement)
140 maxabsdmatelement = absval;
141 absval = template_blas_fabs(dens[b*n+c]);
142 if(absval > maxabsdmatelement)
143 maxabsdmatelement = absval;
144 absval = template_blas_fabs(dens[b*n+d]);
145 if(absval > maxabsdmatelement)
146 maxabsdmatelement = absval;
147 }
148 } // END FOR i j get largest dmat element
149
150 int cluster_i_start = batchList_1[batch_i].clusterStartIndex;
151 int clusterCount1 = batchList_1[batch_i].noOfClusters;
152 for(int cluster_i = cluster_i_start; cluster_i < cluster_i_start + clusterCount1; cluster_i++)
153 {
154 int cluster_j_start = batchList_2[batch_j].clusterStartIndex;
155 int clusterCount2 = batchList_2[batch_j].noOfClusters;
156 int cluterIndexEnd2 = cluster_j_start + clusterCount2;
157 if(interactionWithSelf == 1 && batch_i == batch_j)
158 cluster_j_start = cluster_i;
159 for(int cluster_j = cluster_j_start; cluster_j < cluterIndexEnd2; cluster_j++)
160 {
161 // check if we can skip this combination of clusters
162 if(clusterList_1[cluster_i].maxLimitingFactorForCluster * clusterList_2[cluster_j].maxLimitingFactorForCluster * maxabsdmatelement < threshold)
163 continue;
164
165 int group_i_start = clusterList_1[cluster_i].groupStartIndex;
166 int group_i_end = group_i_start + clusterList_1[cluster_i].noOfGroups;
167 int group_j_start = clusterList_2[cluster_j].groupStartIndex;
168 int group_j_end = group_j_start + clusterList_2[cluster_j].noOfGroups;
169
170 int n1max = clusterList_1[cluster_i].nmax;
171 int n2max = clusterList_2[cluster_j].nmax;
172
173 // Now we can precompute things that depend only on exponents
174 ergo_real alpha_1 = groupList_1[group_i_start].exponent;
175 ergo_real alpha_2 = groupList_2[group_j_start].exponent;
176 ergo_real alphasum = alpha_1 + alpha_2;
177 ergo_real alphaproduct = alpha_1 * alpha_2;
178 ergo_real alpha_0 = alphaproduct / alphasum;
179
180 ergo_real resultPreFactor = twoTimesPiToPow5half / (alphaproduct*template_blas_sqrt(alphasum));
181
182 for(int group_i = group_i_start; group_i < group_i_end; group_i++)
183 {
184 if(interactionWithSelf == 1 && batch_i == batch_j && cluster_i == cluster_j)
185 group_j_start = group_i;
186 for(int group_j = group_j_start; group_j < group_j_end; group_j++)
187 {
188 if(K == NULL)
189 {
190 // Only J is considered; we can use maxAbsDmatElementGroup
191 ergo_real maxabs_1 = groupList_1[group_i].maxAbsDmatElementGroup;
192 ergo_real maxabs_2 = groupList_2[group_j].maxAbsDmatElementGroup;
193 if((groupList_1[group_i].maxLimitingFactorGroup * groupList_2[group_j].maxLimitingFactorGroup * maxabs_1 < threshold) &&
194 (groupList_1[group_i].maxLimitingFactorGroup * groupList_2[group_j].maxLimitingFactorGroup * maxabs_2 < threshold))
195 continue;
196 }
197 else
198 {
199 if(groupList_1[group_i].maxLimitingFactorGroup * groupList_2[group_j].maxLimitingFactorGroup * maxabsdmatelement < threshold)
200 continue;
201 }
202
203 // now we can do all integrals needed for this pair of groups
204 ergo_real dx = groupList_2[group_j].centerCoords[0] - groupList_1[group_i].centerCoords[0];
205 ergo_real dy = groupList_2[group_j].centerCoords[1] - groupList_1[group_i].centerCoords[1];
206 ergo_real dz = groupList_2[group_j].centerCoords[2] - groupList_1[group_i].centerCoords[2];
207
208 // now we have dx dy dz alpha0 alpha1 n1max n2max. Get all integrals for this case.
209 int noOfMonomials_1 = integralInfo.monomial_info.no_of_monomials_list[n1max];
210 int noOfMonomials_2 = integralInfo.monomial_info.no_of_monomials_list[n2max];
211
212 if(get_related_integrals_h(integralInfo,
213 CAM_params_not_used,
214 n1max, noOfMonomials_1,
215 n2max, noOfMonomials_2,
216 dx, dy, dz, alpha_1, alpha_2, alpha_0,
217 primitiveIntegralList,
218 primitiveIntegralList_work,
219 resultPreFactor
220 ) != 0)
221 {
222 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in get_related_integrals");
223 return -1;
224 }
225
226 int i_start = groupList_1[group_i].startIndex;
227 int i_end = i_start + groupList_1[group_i].distrCount;
228 for(int i = i_start; i < i_end; i++)
229 {
230 int idx_1 = minimalDistrList_1[i].basisFuncPairIndex;
231 int monomialIndex_1 = minimalDistrList_1[i].monomialIndex;
232 int j_start = groupList_2[group_j].startIndex;
233 int j_end = j_start + groupList_2[group_j].distrCount;
234 if(interactionWithSelf == 1 && group_j == group_i && batch_i == batch_j && cluster_i == cluster_j)
235 {
236 // take care of case i = j separately
237 ergo_real integralValue = primitiveIntegralList[monomialIndex_1*noOfMonomials_2+monomialIndex_1];
238 ergo_real integralValueCurr = minimalDistrList_1[i].coeff * minimalDistrList_1[i].coeff * integralValue;
239 integralValueCurr *= 0.5;
240 summedIntegralList[idx_1*noOfBasisFuncPairs_2 + idx_1] += integralValueCurr;
241 j_start = i+1;
242 }
243 for(int j = j_start; j < j_end; j++)
244 {
245 int idx_2 = minimalDistrList_2[j].basisFuncPairIndex;
246 int monomialIndex_2 = minimalDistrList_2[j].monomialIndex;
247 ergo_real integralValue = primitiveIntegralList[monomialIndex_1*noOfMonomials_2+monomialIndex_2];
248 ergo_real integralValueCurr = minimalDistrList_1[i].coeff * minimalDistrList_2[j].coeff * integralValue;
249 summedIntegralList[idx_1*noOfBasisFuncPairs_2 + idx_2] += integralValueCurr;
250 } // END FOR j
251 } // END FOR i
252
253 } // END FOR group_j
254 } // END FOR group_i
255 } // END FOR cluster_j
256 } // END FOR cluster_i
257
258 for(int idx_1 = 0; idx_1 < noOfBasisFuncPairs_1; idx_1++)
259 for(int idx_2 = 0; idx_2 < noOfBasisFuncPairs_2; idx_2++)
260 {
261 int a = basisFuncPairList_1[batchList_1[batch_i].basisFuncPairListIndex+idx_1].index_1;
262 int b = basisFuncPairList_1[batchList_1[batch_i].basisFuncPairListIndex+idx_1].index_2;
263 int c = basisFuncPairList_2[batchList_2[batch_j].basisFuncPairListIndex+idx_2].index_1;
264 int d = basisFuncPairList_2[batchList_2[batch_j].basisFuncPairListIndex+idx_2].index_2;
265 ergo_real integralValueCurr = summedIntegralList[idx_1*noOfBasisFuncPairs_2 + idx_2];
266
267 if(a == c && b == d)
268 integralValueCurr *= 2;
269
270 if(template_blas_fabs(integralValueCurr)*maxabsdmatelement < threshold)
271 continue;
272
273 if(a != b && c != d && a != c && a != d && b != c && b != d)
274 {
275 if(J != NULL)
276 {
277 J[a*n+b] += 2 * dens[c*n+d] * integralValueCurr;
278 J[c*n+d] += 2 * dens[a*n+b] * integralValueCurr;
279 }
280 if(K != NULL)
281 {
282 if(d >= a)
283 K[a*n+d] += -0.5 * dens[b*n+c] * integralValueCurr;
284 else
285 K[d*n+a] += -0.5 * dens[b*n+c] * integralValueCurr;
286 if(c >= a)
287 K[a*n+c] += -0.5 * dens[b*n+d] * integralValueCurr;
288 else
289 K[c*n+a] += -0.5 * dens[b*n+d] * integralValueCurr;
290 if(c >= b)
291 K[b*n+c] += -0.5 * dens[a*n+d] * integralValueCurr;
292 else
293 K[c*n+b] += -0.5 * dens[a*n+d] * integralValueCurr;
294 if(d >= b)
295 K[b*n+d] += -0.5 * dens[c*n+a] * integralValueCurr;
296 else
297 K[d*n+b] += -0.5 * dens[c*n+a] * integralValueCurr;
298 }
299 }
300 else
301 {
302 abcd_struct list[8];
303
304 /* determine unique configurations */
305 set_abcd_list_item_macro(0, a, b, c, d, 0, 0, 0);
306 set_abcd_list_item_macro(1, a, b, d, c, 0, 0, 0);
307 set_abcd_list_item_macro(2, b, a, c, d, 0, 0, 0);
308 set_abcd_list_item_macro(3, b, a, d, c, 0, 0, 0);
309 set_abcd_list_item_macro(4, c, d, a, b, 0, 0, 0);
310 set_abcd_list_item_macro(5, d, c, a, b, 0, 0, 0);
311 set_abcd_list_item_macro(6, c, d, b, a, 0, 0, 0);
312 set_abcd_list_item_macro(7, d, c, b, a, 0, 0, 0);
313
314 int ccc = 0;
315
316 for(int ii = 0; ii < 8; ii++)
317 {
318 abcd_struct* abcd = &list[ii];
319 int aa, bb, cc, dd;
320
321 /* check if this is a new unique configuration */
322 int unique = 1;
323 for(int jj = 0; jj < ii; jj++)
324 {
325 if(abcd->a == list[jj].a &&
326 abcd->b == list[jj].b &&
327 abcd->c == list[jj].c &&
328 abcd->d == list[jj].d)
329 unique = 0;
330 }
331 if(unique == 0)
332 continue;
333 /* now we know that this configuration is unique. */
334 aa = abcd->a;
335 bb = abcd->b;
336 cc = abcd->c;
337 dd = abcd->d;
338
339 ccc++;
340
341 /* add contribution to coulomb matrix */
342 if(bb >= aa && J != NULL)
343 J[aa*n+bb] += dens[cc*n+dd] * integralValueCurr;
344
345 if(dd >= aa && K != NULL)
346 K[aa*n+dd] += -0.5 * dens[bb*n+cc] * integralValueCurr;
347
348 } /* END FOR ii go through 8 configurations */
349 }
350
351 } // END FOR idx_1 idx_2
352 } // END FOR batch_j
353 } // END FOR batch_i
354
355 return 0;
356 }
357
358
359
360
361
362 typedef struct
363 {
364 int id;
365 ergo_real x[3];
366 } point_3d_struct;
367
368
369
370
371
372 int
compute_JK_single_box(const BasisInfoStruct & basisInfo,const IntegralInfo & integralInfo,ergo_real * J,ergo_real * K,const ergo_real * dens,ergo_real threshold)373 compute_JK_single_box(const BasisInfoStruct & basisInfo,
374 const IntegralInfo & integralInfo,
375 ergo_real* J,
376 ergo_real* K,
377 const ergo_real* dens,
378 ergo_real threshold)
379 {
380 Util::TimeMeter timeMeterTot;
381
382 Util::TimeMeter timeMeterDistrList;
383
384 int n = basisInfo.noOfBasisFuncs;
385 ergo_real maxDensityMatrixElement = get_max_abs_vector_element(n*n, dens);
386
387 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS,
388 "entering compute_JK_single_box, no of basis funcs = %5i, threshold = %7.3g",
389 n, (double)threshold);
390
391 // Require that threshold value is positive.
392 if(threshold <= 0)
393 {
394 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in compute_JK_single_box: (threshold <= 0)");
395 return -1;
396 }
397
398
399 // get largest limiting factor
400 Util::TimeMeter timeMeterTmp1;
401 ergo_real maxLimitingFactor = 0;
402 if(get_list_of_labeled_distrs_maxLimitingFactor(basisInfo,
403 integralInfo,
404 threshold,
405 &maxLimitingFactor,
406 maxDensityMatrixElement) != 0)
407 {
408 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in get_list_of_labeled_distrs_maxLimitingFactor");
409 return -1;
410 }
411 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS,
412 "get_list_of_labeled_distrs_maxLimitingFactor done, maxLimitingFactor = %22.11f",
413 (double)maxLimitingFactor);
414 timeMeterTmp1.print(LOG_AREA_INTEGRALS, "get_list_of_labeled_distrs_maxLimitingFactor");
415
416 // Get number of distributions
417 Util::TimeMeter timeMeterTmp2;
418 int distrCount = get_list_of_labeled_distrs(basisInfo,
419 integralInfo,
420 threshold,
421 NULL,
422 0,
423 maxLimitingFactor,
424 dens,
425 maxDensityMatrixElement);
426 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "getting distrCount done, distrCount = %12i", distrCount);
427 timeMeterTmp2.print(LOG_AREA_INTEGRALS, "get_list_of_labeled_distrs for getting distrCount");
428 if(distrCount == 0)
429 {
430 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "compute_JK_single_box: (distrCount == 0), skipping.");
431 memset(J, 0, n*n*sizeof(ergo_real));
432 memset(K, 0, n*n*sizeof(ergo_real));
433 return 0;
434 }
435 if(distrCount <= 0)
436 {
437 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in compute_JK_single_box: (distrCount <= 0)");
438 return -1;
439 }
440
441 std::vector<DistributionSpecStructLabeled> distrList(distrCount);
442
443 // create list of product primitives, with labels
444 Util::TimeMeter timeMeterTmp3;
445 int distrCountTemp = get_list_of_labeled_distrs(basisInfo,
446 integralInfo,
447 threshold,
448 &distrList[0],
449 distrCount,
450 maxLimitingFactor,
451 dens,
452 maxDensityMatrixElement);
453 if(distrCountTemp != distrCount)
454 {
455 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in compute_JK_single_box:(distrCountTemp != distrCount)");
456 return -1;
457 }
458 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "get_list_of_labeled_distrs done, distrCount = %12i", distrCount);
459 timeMeterTmp3.print(LOG_AREA_INTEGRALS, "get_list_of_labeled_distrs");
460
461
462 // compute extent for all distrs
463 Util::TimeMeter timeMeterComputeExtentForAllDistrs;
464 compute_extent_for_list_of_distributions(distrCount,
465 &distrList[0],
466 threshold,
467 maxLimitingFactor,
468 maxDensityMatrixElement);
469 timeMeterComputeExtentForAllDistrs.print(LOG_AREA_INTEGRALS, "Compute extent for all distrs");
470
471 // get maximum number of monomials
472 int maxNoOfMonomials = 0;
473 for(int i = 0; i < distrCount; i++)
474 {
475 int degree = 0;
476 for(int j = 0; j < 3; j++)
477 degree += distrList[i].distr.monomialInts[j];
478 int noOfMonomials = integralInfo.monomial_info.no_of_monomials_list[degree];
479 if(noOfMonomials > maxNoOfMonomials)
480 maxNoOfMonomials = noOfMonomials;
481 } // END FOR ABcount
482
483 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "distrCount = %i", distrCount);
484
485 std::vector<DistributionSpecStructLabeled> distrList2(distrCount);
486 int jcounter = 0;
487 for(int i = 0; i < distrCount; i++)
488 {
489 distrList2[jcounter] = distrList[i];
490 jcounter++;
491 }
492 distrCount = jcounter;
493 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "distrCount = %i (after removing negligible products)", distrCount);
494
495 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "Creating list of distributions done, distrCount = %9i", distrCount);
496 timeMeterDistrList.print(LOG_AREA_INTEGRALS, "Creating list of distributions");
497
498 #define NUMBER_OF_PARTS 1
499
500 int n_list[NUMBER_OF_PARTS];
501
502 distr_list_description_struct distr_list_description_list[NUMBER_OF_PARTS];
503
504 for(int i = 0; i < NUMBER_OF_PARTS; i++)
505 n_list[i] = 0;
506
507 for(int i = 0; i < distrCount; i++)
508 n_list[i % NUMBER_OF_PARTS]++;
509
510 ergo_real centerCoords[3];
511 memset(centerCoords, 0, 3*sizeof(ergo_real));
512 int currIndex = 0;
513 for(int i = 0; i < NUMBER_OF_PARTS; i++)
514 {
515 if(organize_distributions(integralInfo,
516 &distrList2[currIndex],
517 n_list[i],
518 &distr_list_description_list[i].org,
519 centerCoords,
520 HUGE_INTEGER_NUMBER) != 0)
521 {
522 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in organize_distributions");
523 return -1;
524 }
525 currIndex += n_list[i];
526 }
527
528
529 // Set J to zero
530 memset(J, 0, n*n*sizeof(ergo_real));
531
532 // Set K to zero
533 memset(K, 0, n*n*sizeof(ergo_real));
534
535
536 // Allocate buffers needed by integral code
537 JK_contribs_buffer_struct bufferStruct;
538 allocate_buffers_needed_by_integral_code(integralInfo, maxNoOfMonomials, 0, &bufferStruct);
539
540 for(int i = 0; i < NUMBER_OF_PARTS; i++)
541 for(int j = i; j < NUMBER_OF_PARTS; j++)
542 {
543 int self = 0;
544 if(i == j)
545 self = 1;
546 Util::TimeMeter timeMeterJKcontribs;
547 if(get_JK_contribs_from_2_interacting_boxes(basisInfo,
548 integralInfo,
549 maxNoOfMonomials,
550 J,
551 K,
552 dens,
553
554 &distr_list_description_list[i].org.minimalDistrList[0],
555 distr_list_description_list[i].org.groupList.size(),
556 &distr_list_description_list[i].org.groupList[0],
557
558 &distr_list_description_list[j].org.minimalDistrList[0],
559 distr_list_description_list[j].org.groupList.size(),
560 &distr_list_description_list[j].org.groupList[0],
561
562 &distr_list_description_list[i].org.clusterList[0],
563 distr_list_description_list[i].org.clusterList.size(),
564
565 &distr_list_description_list[j].org.clusterList[0],
566 distr_list_description_list[j].org.clusterList.size(),
567
568 &distr_list_description_list[i].org.batchList[0],
569 distr_list_description_list[i].org.batchList.size(),
570
571 &distr_list_description_list[j].org.batchList[0],
572 distr_list_description_list[j].org.batchList.size(),
573
574 &distr_list_description_list[i].org.basisFuncPairList[0],
575 &distr_list_description_list[j].org.basisFuncPairList[0],
576
577 self,
578 threshold,
579 &bufferStruct) != 0)
580 {
581 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in get_JK_contribs_from_2_interacting_boxes");
582 return -1;
583 }
584 timeMeterJKcontribs.print(LOG_AREA_INTEGRALS, "get_JK_contribs_from_2_interacting_boxes for both J and K together");
585 } // END FOR i j
586
587 // Fill the other triangle of K
588 for(int i = 0; i < n; i++)
589 for(int j = 0; j < i; j++)
590 K[i*n+j] = K[j*n+i];
591
592 // Fill the other triangle of J
593 for(int i = 0; i < n; i++)
594 for(int j = 0; j < i; j++)
595 J[i*n+j] = J[j*n+i];
596
597 free_buffers_needed_by_integral_code(&bufferStruct);
598
599 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "compute_JK_single_box ending OK.");
600 timeMeterTot.print(LOG_AREA_INTEGRALS, "compute_JK_single_box");
601
602 return 0;
603 }
604