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_1el_potential.cc
31
32 @brief Code for 1-electron integrals, computation of
33 electron-nuclear potential energy matrix V.
34
35 @author: Elias Rudberg <em>responsible</em>
36 */
37
38 /* Written by Elias Rudberg. First version was written in 2006 or even
39 earlier, at KTH. Then the code has been updated many times after
40 that. */
41
42 /*
43 ELIAS NOTE 2014-05-31: gradient computation has been added now. The
44 code for computing the gradient uses that the nuclear coordinates
45 enter the coputation of the V matrix in two ways:
46
47 (1) Directly in integrals of the type (ij|A) where ij is a basis
48 function pair and A is an atom index. In this case we can get the
49 derivatives directly using get_related_integrals_hermite().
50
51 (2) Via the multipole tree. For this kind of contribution we can
52 compute the derivatives indirectly, by first computing derivatives
53 of th emultipole moments w.r.t. nuclear coordinates, and
54 derivatives of the energy w.r.t. the multipole moments. Then, we
55 multiply those values to get the final derivatives of the energy
56 w.r.t. nuclear coordinates. I.e., we compute the derivative df/dx
57 using the chain rule, df/dx = (df/dy)*(dy/dx).
58 */
59
60 #include <stdlib.h>
61 #include <math.h>
62 #include <stdio.h>
63 #include <errno.h>
64 #include <memory.h>
65 #include <time.h>
66 #include <stdarg.h>
67 #include <assert.h>
68 #include <vector>
69
70 #ifdef _OPENMP
71 #include <omp.h>
72 #endif
73
74 #include "integrals_1el.h"
75 #include "integrals_1el_potential.h"
76 #include "integrals_1el_potential_prep.h"
77 #include "memorymanag.h"
78 #include "pi.h"
79 #include "output.h"
80 #include "utilities.h"
81 #include "boysfunction.h"
82 #include "integral_info.h"
83 #include "integrals_general.h"
84 #include "box_system.h"
85 #include "multipole.h"
86 #include "integrals_2el_single.h"
87 #include "integrals_1el_single.h"
88 #include "integrals_hermite.h"
89 #include "matrix_norm.h"
90 #include "mm_limit_table.h"
91
92
93
94 typedef struct
95 {
96 box_struct_basic basicBox;
97 ergo_real centerOfChargeCoords[3];
98 multipole_struct_large multipole;
99 ergo_real* multipole_moment_derivatives; // If needed, list of derivatives of multipole moments w.r.t. atom coordinates for all atoms in this box.
100 ergo_real* derivatives_wrt_multipole_moments; // If needed, list of derivatives of the energy w.r.t. the multipole moments for this box.
101 } atom_box_struct;
102
103
104 static ergo_real
get_distance_3d(const ergo_real * x,const ergo_real * y)105 get_distance_3d(const ergo_real* x, const ergo_real* y)
106 {
107 ergo_real sum = 0;
108 for(int k = 0; k < 3; k++)
109 {
110 ergo_real dx = x[k] - y[k];
111 sum += dx * dx;
112 }
113 return template_blas_sqrt(sum);
114 }
115
116
get_multipole_contribs_for_atom(multipole_struct_large & boxMultipole,ergo_real * multipolePointCoords,const Atom & currAtom,const MMTranslator & translator)117 static void get_multipole_contribs_for_atom(multipole_struct_large & boxMultipole, ergo_real* multipolePointCoords, const Atom & currAtom, const MMTranslator & translator) {
118 multipole_struct_small multipoleCurrAtom;
119 memset(&multipoleCurrAtom, 0, sizeof(multipole_struct_small));
120 multipoleCurrAtom.degree = 0;
121 multipoleCurrAtom.noOfMoments = 1;
122 multipoleCurrAtom.momentList[0] = currAtom.charge;
123 ergo_real dx = currAtom.coords[0] - multipolePointCoords[0];
124 ergo_real dy = currAtom.coords[1] - multipolePointCoords[1];
125 ergo_real dz = currAtom.coords[2] - multipolePointCoords[2];
126
127 ergo_real W[MAX_NO_OF_MOMENTS_PER_MULTIPOLE*MAX_NO_OF_MOMENTS_PER_MULTIPOLE];
128 translator.getTranslationMatrix(dx, dy, dz, MAX_MULTIPOLE_DEGREE, 0, W);
129
130 multipole_struct_large translatedMultipole;
131 for(int A = 0; A < MAX_NO_OF_MOMENTS_PER_MULTIPOLE; A++)
132 {
133 ergo_real sum = 0;
134 for(int B = 0; B < 1; B++)
135 sum += W[A*1+B] * multipoleCurrAtom.momentList[B];
136 translatedMultipole.momentList[A] = sum;
137 } // END FOR A
138 for(int kk = 0; kk < 3; kk++)
139 translatedMultipole.centerCoords[kk] = multipolePointCoords[kk];
140 translatedMultipole.degree = MAX_MULTIPOLE_DEGREE;
141 translatedMultipole.noOfMoments = MAX_NO_OF_MOMENTS_PER_MULTIPOLE;
142
143 // add translated multipole to box multipole
144 for(int A = 0; A < MAX_NO_OF_MOMENTS_PER_MULTIPOLE; A++)
145 boxMultipole.momentList[A] += translatedMultipole.momentList[A];
146 }
147
148
init_multipole_struct_large(multipole_struct_large & boxMultipole,const ergo_real * multipolePointCoords)149 static void init_multipole_struct_large(multipole_struct_large & boxMultipole, const ergo_real* multipolePointCoords) {
150 memset(&boxMultipole, 0, sizeof(multipole_struct_large));
151 for(int kk = 0; kk < 3; kk++)
152 boxMultipole.centerCoords[kk] = multipolePointCoords[kk];
153 boxMultipole.degree = MAX_MULTIPOLE_DEGREE;
154 boxMultipole.noOfMoments = MAX_NO_OF_MOMENTS_PER_MULTIPOLE;
155 }
156
157
158 static int
create_nuclei_mm_tree(const IntegralInfo & integralInfo,int nAtoms,const Atom * atomList,ergo_real boxSize,BoxSystem & boxSystem,atom_box_struct ** return_boxList,int * return_numberOfLevels,Atom ** return_atomListReordered,int * return_atomPermutation,bool compute_gradient_also)159 create_nuclei_mm_tree(const IntegralInfo& integralInfo,
160 int nAtoms,
161 const Atom* atomList,
162 ergo_real boxSize,
163 BoxSystem & boxSystem,
164 atom_box_struct** return_boxList, // allocated by this routine, must be freed by caller.
165 int *return_numberOfLevels,
166 Atom **return_atomListReordered, // allocated by this routine, must be freed by caller.
167 int* return_atomPermutation, // list of int, length=nAtoms, saying how atoms have been reordered in return_atomListReordered.
168 bool compute_gradient_also
169 )
170 {
171 // Create box system based on atoms
172 box_item_struct* itemList = new box_item_struct[nAtoms];
173 for(int i = 0; i < nAtoms; i++)
174 {
175 for(int j = 0; j < 3; j++)
176 itemList[i].centerCoords[j] = atomList[i].coords[j];
177 itemList[i].originalIndex = i;
178 } // END FOR i
179
180 const ergo_real maxToplevelBoxSize = boxSize;
181
182 if(boxSystem.create_box_system(itemList,
183 nAtoms,
184 maxToplevelBoxSize) != 0)
185 {
186 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in create_box_system");
187 return -1;
188 }
189
190 // Now the itemList has been reordered. Store the ordering info in the return_atomPermutation list.
191 for(int i = 0; i < nAtoms; i++)
192 return_atomPermutation[i] = itemList[i].originalIndex;
193
194 int numberOfLevels = boxSystem.noOfLevels;
195
196 atom_box_struct* boxList = new atom_box_struct[boxSystem.totNoOfBoxes];
197 for(int i = 0; i < boxSystem.totNoOfBoxes; i++)
198 boxList[i].basicBox = boxSystem.boxList[i];
199
200 // create new list of atoms, where they are ordered box by box at the level of smallest boxes.
201 Atom* atomList2 = new Atom[nAtoms];
202 for(int i = 0; i < nAtoms; i++)
203 atomList2[i] = atomList[itemList[i].originalIndex];
204
205 delete [] itemList;
206 itemList = NULL;
207
208 // Find center-of-charge for each box at top level
209 atom_box_struct* boxListTopLevel = &boxList[boxSystem.levelList[numberOfLevels-1].startIndexInBoxList];
210 int noOfBoxesTopLevel = boxSystem.levelList[numberOfLevels-1].noOfBoxes;
211 for(int i = 0; i < noOfBoxesTopLevel; i++)
212 {
213 atom_box_struct* currBox = &boxListTopLevel[i];
214 // Find center-of-charge for this box.
215 ergo_real cocSumList[3];
216 for(int k = 0; k < 3; k++)
217 cocSumList[k] = 0;
218 ergo_real chargeSum = 0;
219 Atom* atomListCurrBox = &atomList2[currBox->basicBox.firstItemIndex];
220 int noOfAtomsCurrBox = currBox->basicBox.noOfItems;
221 int nPositiveCharges = 0;
222 int nNegativeCharges = 0;
223 for(int j = 0; j < noOfAtomsCurrBox; j++)
224 {
225 Atom* currAtom = &atomListCurrBox[j];
226 chargeSum += currAtom->charge;
227 if(currAtom->charge > 0)
228 nPositiveCharges++;
229 if(currAtom->charge < 0)
230 nNegativeCharges++;
231 for(int k = 0; k < 3; k++)
232 cocSumList[k] += currAtom->coords[k] * currAtom->charge;
233 } // END FOR j Find center-of-charge for this box.
234 // We only want to use averaging if all charges are of the same
235 // sign, otherwise it may happen that they cancel out and then
236 // the "average" has no meaning.
237 if(nPositiveCharges == noOfAtomsCurrBox || nNegativeCharges == noOfAtomsCurrBox) {
238 // OK, all charges are of the same sign.
239 for(int k = 0; k < 3; k++)
240 currBox->centerOfChargeCoords[k] = cocSumList[k] / chargeSum;
241 }
242 else {
243 // All charges are not of the same sign. In this case we just use the center of the box as centerOfChargeCoords.
244 for(int k = 0; k < 3; k++)
245 currBox->centerOfChargeCoords[k] = currBox->basicBox.centerCoords[k];
246 }
247
248 if(compute_gradient_also) {
249 // Always use center point in this case, to simplify computation of derivatives.
250 // FIXME: CHECK IF THIS IS REALLY NEEDED, MAYBE IT WAS REALLY BUGS IN OTHER PLACES THAT CAUSED THE PROBLEMS?
251 for(int k = 0; k < 3; k++)
252 currBox->centerOfChargeCoords[k] = currBox->basicBox.centerCoords[k];
253 }
254
255 } // END FOR i Find center-of-charge for each box at top level
256
257 // Create multipole for each box at top level (smallest boxes)
258 MMTranslator translator(integralInfo.GetMultipolePrep());
259 for(int i = 0; i < noOfBoxesTopLevel; i++)
260 {
261 atom_box_struct* currBox = &boxListTopLevel[i];
262 ergo_real* multipolePointCoords = currBox->centerOfChargeCoords;
263 int noOfAtomsCurrBox = currBox->basicBox.noOfItems;
264 if(compute_gradient_also) {
265 currBox->multipole_moment_derivatives = new ergo_real[noOfAtomsCurrBox*MAX_NO_OF_MOMENTS_PER_MULTIPOLE*3];
266 currBox->derivatives_wrt_multipole_moments = new ergo_real[MAX_NO_OF_MOMENTS_PER_MULTIPOLE];
267 memset(currBox->derivatives_wrt_multipole_moments, 0, MAX_NO_OF_MOMENTS_PER_MULTIPOLE*sizeof(ergo_real));
268 }
269 multipole_struct_large boxMultipole;
270 init_multipole_struct_large(boxMultipole, multipolePointCoords);
271 Atom* atomListCurrBox = &atomList2[currBox->basicBox.firstItemIndex];
272 // Go through all atoms in this box
273 for(int j = 0; j < noOfAtomsCurrBox; j++)
274 {
275 const Atom & currAtom = atomListCurrBox[j];
276 // take multipole for this atom, and translate it to center-of-charge point
277 // the "multipole" for this atom is of course only a monopole.
278 get_multipole_contribs_for_atom(boxMultipole, multipolePointCoords, currAtom, translator);
279 if(compute_gradient_also) {
280 for(int coordIdx = 0; coordIdx < 3; coordIdx++) {
281 multipole_struct_large boxMultipoleTmp1;
282 init_multipole_struct_large(boxMultipoleTmp1, multipolePointCoords);
283 multipole_struct_large boxMultipoleTmp2;
284 init_multipole_struct_large(boxMultipoleTmp2, multipolePointCoords);
285 const ergo_real eps = 1e-5;
286 Atom atomTmp1 = currAtom;
287 atomTmp1.coords[coordIdx] += eps;
288 Atom atomTmp2 = currAtom;
289 atomTmp2.coords[coordIdx] -= eps;
290 get_multipole_contribs_for_atom(boxMultipoleTmp1, multipolePointCoords, atomTmp1, translator);
291 get_multipole_contribs_for_atom(boxMultipoleTmp2, multipolePointCoords, atomTmp2, translator);
292 for(int ii = 0; ii < MAX_NO_OF_MOMENTS_PER_MULTIPOLE; ii++) {
293 ergo_real m1 = boxMultipoleTmp1.momentList[ii];
294 ergo_real m2 = boxMultipoleTmp2.momentList[ii];
295 ergo_real value = (m1 - m2) / (2*eps);
296 currBox->multipole_moment_derivatives[j*MAX_NO_OF_MOMENTS_PER_MULTIPOLE*3 + MAX_NO_OF_MOMENTS_PER_MULTIPOLE*coordIdx + ii] = value;
297 }
298 }
299 }
300 } // END FOR j Go through all atoms in this box
301
302 currBox->multipole = boxMultipole;
303 } // END FOR i Create multipole for each box at top level (smallest boxes)
304
305
306
307 // OK, multipoles created for top level.
308 // Now go through the other levels, joining multipoles from child boxes to a single multipole in parent box
309
310 for(int levelNumber = numberOfLevels-2; levelNumber >= 0; levelNumber--)
311 {
312 int noOfBoxesCurrLevel = boxSystem.levelList[levelNumber].noOfBoxes;
313 atom_box_struct* boxListCurrLevel = &boxList[boxSystem.levelList[levelNumber].startIndexInBoxList];
314
315 for(int boxIndex = 0; boxIndex < noOfBoxesCurrLevel; boxIndex++)
316 {
317 atom_box_struct* currBox = &boxListCurrLevel[boxIndex];
318 int noOfChildren = currBox->basicBox.noOfChildBoxes;
319
320 if(noOfChildren == 0)
321 {
322 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "ERROR: (noOfChildren == 0)");
323 return -1;
324 }
325
326 multipole_struct_large* newMultipole = &currBox->multipole;
327 memset(newMultipole, 0, sizeof(multipole_struct_large));
328
329 // get average position of child multipoles
330 ergo_real avgPosList[3];
331 int kk;
332 for(kk = 0; kk < 3; kk++)
333 avgPosList[kk] = 0;
334 int childIndex;
335 for(childIndex = 0; childIndex < noOfChildren; childIndex++)
336 {
337 int childIndexInBoxList = currBox->basicBox.firstChildBoxIndex + childIndex;
338 atom_box_struct* childBox = &boxList[childIndexInBoxList];
339 for(kk = 0; kk < 3; kk++)
340 avgPosList[kk] += childBox->multipole.centerCoords[kk];
341 } // END FOR childIndex
342 for(kk = 0; kk < 3; kk++)
343 newMultipole->centerCoords[kk] = avgPosList[kk] / noOfChildren;
344
345 if(compute_gradient_also) {
346 // Always use center point in this case, to simplify computation of derivatives.
347 // FIXME: CHECK IF THIS IS REALLY NEEDED, MAYBE IT WAS REALLY BUGS IN OTHER PLACES THAT CAUSED THE PROBLEMS?
348 for(kk = 0; kk < 3; kk++)
349 newMultipole->centerCoords[kk] = currBox->basicBox.centerCoords[kk];
350 }
351
352 newMultipole->degree = MAX_MULTIPOLE_DEGREE;
353 newMultipole->noOfMoments = MAX_NO_OF_MOMENTS_PER_MULTIPOLE;
354
355 // Now translate child multipoles and add to parent multipole
356 for(childIndex = 0; childIndex < noOfChildren; childIndex++)
357 {
358 int childIndexInBoxList = currBox->basicBox.firstChildBoxIndex + childIndex;
359 atom_box_struct* childBox = &boxList[childIndexInBoxList];
360 multipole_struct_large* childMultipole = &childBox->multipole;
361
362 ergo_real dx = childMultipole->centerCoords[0] - newMultipole->centerCoords[0];
363 ergo_real dy = childMultipole->centerCoords[1] - newMultipole->centerCoords[1];
364 ergo_real dz = childMultipole->centerCoords[2] - newMultipole->centerCoords[2];
365
366 ergo_real W[MAX_NO_OF_MOMENTS_PER_MULTIPOLE*MAX_NO_OF_MOMENTS_PER_MULTIPOLE];
367 translator.getTranslationMatrix(dx, dy, dz, MAX_MULTIPOLE_DEGREE, MAX_MULTIPOLE_DEGREE, W);
368
369 multipole_struct_large translatedMultipole;
370 int A, B;
371 for(A = 0; A < MAX_NO_OF_MOMENTS_PER_MULTIPOLE; A++)
372 {
373 ergo_real sum = 0;
374 for(B = 0; B < MAX_NO_OF_MOMENTS_PER_MULTIPOLE; B++)
375 sum += W[A*MAX_NO_OF_MOMENTS_PER_MULTIPOLE+B] * childMultipole->momentList[B];
376 translatedMultipole.momentList[A] = sum;
377 } // END FOR A
378 for(kk = 0; kk < 3; kk++)
379 translatedMultipole.centerCoords[kk] = newMultipole->centerCoords[kk];
380 translatedMultipole.degree = MAX_MULTIPOLE_DEGREE;
381 translatedMultipole.noOfMoments = MAX_NO_OF_MOMENTS_PER_MULTIPOLE;
382
383 // add translated multipole to parent multipole
384 for(A = 0; A < MAX_NO_OF_MOMENTS_PER_MULTIPOLE; A++)
385 newMultipole->momentList[A] += translatedMultipole.momentList[A];
386 } // END FOR childIndex
387
388 if(compute_gradient_also) {
389 // Get multipole for large box by going through all atoms instead of translating child multipoles
390 ergo_real* multipolePointCoords = newMultipole->centerCoords;
391 multipole_struct_large boxMultipole;
392 init_multipole_struct_large(boxMultipole, multipolePointCoords);
393 Atom* atomListCurrBox = &atomList2[currBox->basicBox.firstItemIndex];
394 int noOfAtomsCurrBox = currBox->basicBox.noOfItems;
395 currBox->multipole_moment_derivatives = new ergo_real[noOfAtomsCurrBox*MAX_NO_OF_MOMENTS_PER_MULTIPOLE*3];
396 currBox->derivatives_wrt_multipole_moments = new ergo_real[MAX_NO_OF_MOMENTS_PER_MULTIPOLE];
397 memset(currBox->derivatives_wrt_multipole_moments, 0, MAX_NO_OF_MOMENTS_PER_MULTIPOLE*sizeof(ergo_real));
398 // Go through all atoms in this box
399 for(int j = 0; j < noOfAtomsCurrBox; j++) {
400 const Atom & currAtom = atomListCurrBox[j];
401 // take multipole for this atom, and translate it to center-of-charge point
402 // the "multipole" for this atom is of course only a monopole.
403 get_multipole_contribs_for_atom(boxMultipole, multipolePointCoords, currAtom, translator);
404 for(int coordIdx = 0; coordIdx < 3; coordIdx++) {
405 multipole_struct_large boxMultipoleTmp1;
406 init_multipole_struct_large(boxMultipoleTmp1, multipolePointCoords);
407 multipole_struct_large boxMultipoleTmp2;
408 init_multipole_struct_large(boxMultipoleTmp2, multipolePointCoords);
409 const ergo_real eps = 1e-5;
410 Atom atomTmp1 = currAtom;
411 atomTmp1.coords[coordIdx] += eps;
412 Atom atomTmp2 = currAtom;
413 atomTmp2.coords[coordIdx] -= eps;
414 get_multipole_contribs_for_atom(boxMultipoleTmp1, multipolePointCoords, atomTmp1, translator);
415 get_multipole_contribs_for_atom(boxMultipoleTmp2, multipolePointCoords, atomTmp2, translator);
416 for(int ii = 0; ii < MAX_NO_OF_MOMENTS_PER_MULTIPOLE; ii++) {
417 ergo_real m1 = boxMultipoleTmp1.momentList[ii];
418 ergo_real m2 = boxMultipoleTmp2.momentList[ii];
419 ergo_real value = (m1 - m2) / (2*eps);
420 currBox->multipole_moment_derivatives[j*MAX_NO_OF_MOMENTS_PER_MULTIPOLE*3 + MAX_NO_OF_MOMENTS_PER_MULTIPOLE*coordIdx + ii] = value;
421 }
422 }
423 } // END FOR j Go through all atoms in this box
424 }
425
426 } // END FOR boxIndex
427 } // END FOR levelNumber
428
429
430 // Prepare info needed later to determine needed multipole degree
431 for(int levelNumber = 0; levelNumber < numberOfLevels; levelNumber++)
432 {
433 int noOfBoxesCurrLevel = boxSystem.levelList[levelNumber].noOfBoxes;
434 atom_box_struct* boxListCurrLevel = &boxList[boxSystem.levelList[levelNumber].startIndexInBoxList];
435 for(int boxIndex = 0; boxIndex < noOfBoxesCurrLevel; boxIndex++)
436 {
437 atom_box_struct* currBox = &boxListCurrLevel[boxIndex];
438 setup_multipole_maxAbsMomentList(&currBox->multipole);
439 }
440 }
441
442
443 *return_boxList = boxList;
444 *return_numberOfLevels = numberOfLevels;
445 *return_atomListReordered = atomList2;
446 return 0;
447 }
448
449
450
451 typedef struct {
452 DistributionSpecStruct distr;
453 int pairIdx;
454 int basisFuncIdx1;
455 int basisFuncIdx2;
456 } DistributionSpecStructWithIndexes;
457
458
459 /**
460 Take care of interaction between list of distrs and box.
461 */
462 static int
do_interaction_recursive(const IntegralInfo & integralInfo,ergo_real * V_list,int noOfBasisFuncIndexPairs,const basis_func_index_pair_struct_1el * basisFuncIndexPairList,const DistributionSpecStructWithIndexes * list,int nDistrs,const multipole_struct_small * multipoleList,const ergo_real * maxMomentVectorNormForDistrsList,int maxNoOfMomentsForDistrs,int maxDegreeForDistrs,ergo_real distrExtent,const Atom * atomListReordered,const int * atomPermutation,ergo_real threshold,const atom_box_struct * boxList,MMInteractor & interactor,int boxIndex,int currLevel,int numberOfLevels,bool compute_gradient_also,const ergo_real * D_list,ergo_real * result_gradient_list)463 do_interaction_recursive(const IntegralInfo & integralInfo,
464 ergo_real* V_list,
465 int noOfBasisFuncIndexPairs,
466 const basis_func_index_pair_struct_1el* basisFuncIndexPairList,
467 const DistributionSpecStructWithIndexes* list,
468 int nDistrs,
469 const multipole_struct_small* multipoleList,
470 const ergo_real* maxMomentVectorNormForDistrsList,
471 int maxNoOfMomentsForDistrs,
472 int maxDegreeForDistrs,
473 ergo_real distrExtent,
474 const Atom *atomListReordered,
475 const int* atomPermutation, // list of int, length=nAtoms, saying how atoms have been reordered in return_atomListReordered.
476 ergo_real threshold,
477 const atom_box_struct* boxList,
478 MMInteractor & interactor,
479 int boxIndex,
480 int currLevel,
481 int numberOfLevels,
482 bool compute_gradient_also,
483 const ergo_real* D_list, // used in compute_gradient_also case, NULL otherwise
484 ergo_real* result_gradient_list // used in compute_gradient_also case, NULL otherwise
485 )
486 {
487 const atom_box_struct* currBox = &boxList[boxIndex];
488 const multipole_struct_large* boxMultipole = &currBox->multipole;
489
490 // check if current box is far enough away so that we can use multipole description.
491 ergo_real distance = get_distance_3d(list[0].distr.centerCoords, currBox->basicBox.centerCoords);
492 ergo_real boxRadius = currBox->basicBox.width * 0.5 * template_blas_sqrt((ergo_real)3);
493 ergo_real requiredDistance = boxRadius + distrExtent;
494
495 // Note that the distance to the box multipole is different, since
496 // it is not necessarily placed at the box center.
497 ergo_real multipoleDistance = get_distance_3d(list[0].distr.centerCoords, boxMultipole->centerCoords);
498
499 int degreeNeeded = integralInfo.GetMMLimitTable().get_minimum_multipole_degree_needed(multipoleDistance,
500 boxMultipole,
501 maxDegreeForDistrs,
502 maxMomentVectorNormForDistrsList,
503 threshold);
504 if(degreeNeeded < 0)
505 return -1;
506 degreeNeeded+=2; // We need a couple of extra degrees to handle gradient computation. FIXME: IS THIS REALL NEEDED? WHY? MAYBE IT WAS BUGS IN OTHER PLACES THAT CAUSED THE PROBLEMS?
507
508 bool multipoleDegreeIsSafe = false;
509 // Demand at least two degrees margin compared to
510 // MAX_MULTIPOLE_DEGREE, otherwise this may fail due to alternating
511 // odd/even degrees where only one of them are significant. This has
512 // happened in some test cases.
513 if(degreeNeeded <= MAX_MULTIPOLE_DEGREE-2)
514 multipoleDegreeIsSafe = true;
515
516 if(distance > requiredDistance && multipoleDegreeIsSafe)
517 {
518 // OK, use multipole description of atom charges.
519 int boxNeededNoOfMoments = (degreeNeeded+1)*(degreeNeeded+1);
520 // create interaction matrix
521 ergo_real T[boxNeededNoOfMoments * maxNoOfMomentsForDistrs];
522 ergo_real dx = boxMultipole->centerCoords[0] - list[0].distr.centerCoords[0];
523 ergo_real dy = boxMultipole->centerCoords[1] - list[0].distr.centerCoords[1];
524 ergo_real dz = boxMultipole->centerCoords[2] - list[0].distr.centerCoords[2];
525
526 interactor.getInteractionMatrix(dx, dy, dz, maxDegreeForDistrs, degreeNeeded, T);
527
528 ergo_real tempVector[MAX_NO_OF_MOMENTS_PER_MULTIPOLE];
529 for(int A = 0; A < maxNoOfMomentsForDistrs; A++)
530 {
531 ergo_real sum = 0;
532 for(int B = 0; B < boxNeededNoOfMoments; B++)
533 {
534 ergo_real momB = boxMultipole->momentList[B];
535 ergo_real Telement = T[A*boxNeededNoOfMoments+B];
536 sum += momB * Telement;
537 }
538 tempVector[A] = sum;
539 }
540
541 for(int i = 0; i < nDistrs; i++)
542 {
543 ergo_real sum = 0;
544 for(int A = 0; A < multipoleList[i].noOfMoments; A++)
545 sum += tempVector[A] * multipoleList[i].momentList[A];
546 V_list[list[i].pairIdx] += -1 * sum;
547 }
548
549 if(compute_gradient_also) {
550 // we need to compute the derivatives of the energy with respect to the multipole moments.
551 // To reduce number of critical sections when OpenMP threading is used, use temporary derivatives_wrt_multipole_moments list first, and then update the real one.
552 ergo_real derivatives_wrt_multipole_moments_tmp[MAX_NO_OF_MOMENTS_PER_MULTIPOLE];
553 memset(derivatives_wrt_multipole_moments_tmp, 0, MAX_NO_OF_MOMENTS_PER_MULTIPOLE*sizeof(ergo_real));
554 for(int i = 0; i < nDistrs; i++) {
555 ergo_real tempVector2[boxNeededNoOfMoments];
556 for(int B = 0; B < boxNeededNoOfMoments; B++) {
557 ergo_real sum = 0;
558 for(int A = 0; A < multipoleList[i].noOfMoments; A++) {
559 ergo_real Telement = T[A*boxNeededNoOfMoments+B];
560 ergo_real mom = multipoleList[i].momentList[A];
561 sum += mom * Telement;
562 }
563 tempVector2[B] = sum;
564 } // end for B
565 int pairIdx = list[i].pairIdx;
566 ergo_real extraFactor = 1.0;
567 if(basisFuncIndexPairList[pairIdx].index_1 != basisFuncIndexPairList[pairIdx].index_2)
568 extraFactor = 2.0;
569 ergo_real densityMatrixElement = D_list[list[i].pairIdx];
570 for(int B = 0; B < boxNeededNoOfMoments; B++)
571 derivatives_wrt_multipole_moments_tmp[B] += -1 * extraFactor * densityMatrixElement * tempVector2[B];
572 } // end for i (loop over distrs)
573
574 #ifdef _OPENMP
575 #pragma omp critical
576 #endif
577 {
578 for(int B = 0; B < boxNeededNoOfMoments; B++)
579 currBox->derivatives_wrt_multipole_moments[B] += derivatives_wrt_multipole_moments_tmp[B];
580 }
581
582 } // end if compute_gradient_also
583
584 return 0;
585 }
586
587 // No, multipole description could not be used in this case.
588 if(currLevel == numberOfLevels-1)
589 {
590 // We are at top level, must compute explicit interactions
591
592 int Nmax = 0;
593 for(int i = 0; i < nDistrs; i++)
594 {
595 const DistributionSpecStruct* distr = &list[i].distr;
596 int N = distr->monomialInts[0] + distr->monomialInts[1] + distr->monomialInts[2];
597 if(N > Nmax)
598 Nmax = N;
599 }
600
601 const JK::ExchWeights CAM_params_not_used;
602
603 const Atom* currAtomList = &atomListReordered[currBox->basicBox.firstItemIndex];
604 int noOfAtomsCurrBox = currBox->basicBox.noOfItems;
605 for(int ia = 0; ia < noOfAtomsCurrBox; ia++)
606 {
607 const Atom* currAtom = &currAtomList[ia];
608
609 const DistributionSpecStruct* distr1 = &list[0].distr;
610 int noOfMonomials_1 = integralInfo.monomial_info.no_of_monomials_list[Nmax];
611 int noOfMonomials_2 = integralInfo.monomial_info.no_of_monomials_list[0];
612 ergo_real primitiveIntegralList_h[noOfMonomials_1*noOfMonomials_2];
613 ergo_real primitiveIntegralList_2[noOfMonomials_1*noOfMonomials_2];
614 // Let the distr be 1 and the pointcharge 2
615 ergo_real alpha1 = distr1->exponent;
616 ergo_real alpha0 = alpha1;
617 int n1 = Nmax;
618 int n2 = 0;
619 ergo_real dx0 = currAtom->coords[0] - distr1->centerCoords[0];
620 ergo_real dx1 = currAtom->coords[1] - distr1->centerCoords[1];
621 ergo_real dx2 = currAtom->coords[2] - distr1->centerCoords[2];
622 ergo_real resultPreFactor = 2 * pi / alpha1;
623 get_related_integrals_hermite(integralInfo,
624 CAM_params_not_used,
625 n1, noOfMonomials_1,
626 n2, noOfMonomials_2,
627 dx0,
628 dx1,
629 dx2,
630 alpha0,
631 resultPreFactor,
632 primitiveIntegralList_h);
633 integralInfo.multiply_by_hermite_conversion_matrix_from_right(n1,
634 n2,
635 1.0/alpha1,
636 primitiveIntegralList_h,
637 primitiveIntegralList_2);
638 for(int id = 0; id < nDistrs; id++)
639 {
640 const DistributionSpecStruct* distr = &list[id].distr;
641 int n1x = distr->monomialInts[0];
642 int n1y = distr->monomialInts[1];
643 int n1z = distr->monomialInts[2];
644 int monomialIndex = integralInfo.monomial_info.monomial_index_list[n1x][n1y][n1z];
645 ergo_real integralValue = currAtom->charge * distr->coeff * primitiveIntegralList_2[monomialIndex];
646 V_list[list[id].pairIdx] += -1 * integralValue;
647 }
648
649 } // END FOR ia
650
651 if(compute_gradient_also) {
652 for(int ia = 0; ia < noOfAtomsCurrBox; ia++) {
653 const Atom* currAtom = &currAtomList[ia];
654 const DistributionSpecStruct* distr1 = &list[0].distr;
655 int n1 = Nmax;
656 int n2 = 1;
657 int noOfMonomials_1 = integralInfo.monomial_info.no_of_monomials_list[n1];
658 int noOfMonomials_2 = integralInfo.monomial_info.no_of_monomials_list[n2];
659 ergo_real primitiveIntegralList_h[noOfMonomials_1*noOfMonomials_2];
660 // Let the distr be 1 and the pointcharge 2
661 ergo_real alpha1 = distr1->exponent;
662 ergo_real alpha0 = alpha1;
663 ergo_real dx0 = currAtom->coords[0] - distr1->centerCoords[0];
664 ergo_real dx1 = currAtom->coords[1] - distr1->centerCoords[1];
665 ergo_real dx2 = currAtom->coords[2] - distr1->centerCoords[2];
666 ergo_real resultPreFactor = 2 * pi / alpha1;
667 get_related_integrals_hermite(integralInfo,
668 CAM_params_not_used,
669 n1, noOfMonomials_1,
670 n2, noOfMonomials_2,
671 dx0,
672 dx1,
673 dx2,
674 alpha0,
675 resultPreFactor,
676 primitiveIntegralList_h);
677 for(int id = 0; id < nDistrs; id++) {
678 const DistributionSpecStruct* distr = &list[id].distr;
679 int n1b = n1;
680 int n2b = 0;
681 ergo_real primitiveIntegralList_h_components[3][noOfMonomials_1];
682 int monomialIndex_x = integralInfo.monomial_info.monomial_index_list[1][0][0];
683 int monomialIndex_y = integralInfo.monomial_info.monomial_index_list[0][1][0];
684 int monomialIndex_z = integralInfo.monomial_info.monomial_index_list[0][0][1];
685 for(int i = 0; i < noOfMonomials_1; i++) {
686 primitiveIntegralList_h_components[0][i] = primitiveIntegralList_h[i*noOfMonomials_2+monomialIndex_x];
687 primitiveIntegralList_h_components[1][i] = primitiveIntegralList_h[i*noOfMonomials_2+monomialIndex_y];
688 primitiveIntegralList_h_components[2][i] = primitiveIntegralList_h[i*noOfMonomials_2+monomialIndex_z];
689 }
690 ergo_real primitiveIntegralList_2_components[3][noOfMonomials_1];
691 for(int i = 0; i < 3; i++)
692 integralInfo.multiply_by_hermite_conversion_matrix_from_right(n1b, n2b, 1.0/alpha1, primitiveIntegralList_h_components[i], primitiveIntegralList_2_components[i]);
693 int n1x = distr->monomialInts[0];
694 int n1y = distr->monomialInts[1];
695 int n1z = distr->monomialInts[2];
696 int monomialIndex = integralInfo.monomial_info.monomial_index_list[n1x][n1y][n1z];
697 ergo_real value_x = currAtom->charge * distr->coeff * primitiveIntegralList_2_components[0][monomialIndex];
698 ergo_real value_y = currAtom->charge * distr->coeff * primitiveIntegralList_2_components[1][monomialIndex];
699 ergo_real value_z = currAtom->charge * distr->coeff * primitiveIntegralList_2_components[2][monomialIndex];
700 int pairIdx = list[id].pairIdx;
701 ergo_real extraFactor = 1.0;
702 if(basisFuncIndexPairList[pairIdx].index_1 != basisFuncIndexPairList[pairIdx].index_2)
703 extraFactor = 2.0;
704 ergo_real densityMatrixElement = D_list[list[id].pairIdx];
705 int atomIndex = atomPermutation[currBox->basicBox.firstItemIndex + ia];
706 #ifdef _OPENMP
707 #pragma omp critical
708 #endif
709 {
710 result_gradient_list[atomIndex*3+0] += -1 * value_x * densityMatrixElement * extraFactor;
711 result_gradient_list[atomIndex*3+1] += -1 * value_y * densityMatrixElement * extraFactor;
712 result_gradient_list[atomIndex*3+2] += -1 * value_z * densityMatrixElement * extraFactor;
713 }
714 }
715
716 }
717
718 }
719
720 return 0;
721 }
722
723 // Go to next level
724 int noOfChildren = currBox->basicBox.noOfChildBoxes;
725 for(int i = 0; i < noOfChildren; i++)
726 {
727 int childBoxIndex = currBox->basicBox.firstChildBoxIndex + i;
728 if(do_interaction_recursive(integralInfo,
729 V_list,
730 noOfBasisFuncIndexPairs,
731 basisFuncIndexPairList,
732 list,
733 nDistrs,
734 multipoleList,
735 maxMomentVectorNormForDistrsList,
736 maxNoOfMomentsForDistrs,
737 maxDegreeForDistrs,
738 distrExtent,
739 atomListReordered,
740 atomPermutation,
741 threshold,
742 boxList,
743 interactor,
744 childBoxIndex,
745 currLevel + 1,
746 numberOfLevels,
747 compute_gradient_also,
748 D_list,
749 result_gradient_list) != 0)
750 return -1;
751 } // END FOR i
752 return 0;
753 }
754
755
756
757 /**
758 Take care of interaction between list of distrs and box.
759 */
760 static int
do_interaction_recursive_2(const IntegralInfo & integralInfo,csr_matrix_struct * V_CSR,int noOfBasisFuncIndexPairs,const basis_func_index_pair_struct_1el * basisFuncIndexPairList,const DistributionSpecStructWithIndexes2 * list,int nDistrs,const multipole_struct_small * multipoleList,const ergo_real * maxMomentVectorNormForDistrsList,int maxNoOfMomentsForDistrs,int maxDegreeForDistrs,ergo_real distrExtent,const Atom * atomListReordered,const int * atomPermutation,ergo_real threshold,const atom_box_struct * boxList,MMInteractor & interactor,int boxIndex,int currLevel,int numberOfLevels)761 do_interaction_recursive_2(const IntegralInfo & integralInfo,
762 csr_matrix_struct* V_CSR,
763 int noOfBasisFuncIndexPairs,
764 const basis_func_index_pair_struct_1el* basisFuncIndexPairList,
765 const DistributionSpecStructWithIndexes2* list,
766 int nDistrs,
767 const multipole_struct_small* multipoleList,
768 const ergo_real* maxMomentVectorNormForDistrsList,
769 int maxNoOfMomentsForDistrs,
770 int maxDegreeForDistrs,
771 ergo_real distrExtent,
772 const Atom *atomListReordered,
773 const int* atomPermutation, // list of int, length=nAtoms, saying how atoms have been reordered in return_atomListReordered.
774 ergo_real threshold,
775 const atom_box_struct* boxList,
776 MMInteractor & interactor,
777 int boxIndex,
778 int currLevel,
779 int numberOfLevels
780 )
781 {
782 const atom_box_struct* currBox = &boxList[boxIndex];
783 const multipole_struct_large* boxMultipole = &currBox->multipole;
784
785 // check if current box is far enough away so that we can use multipole description.
786 ergo_real distance = get_distance_3d(list[0].distr.centerCoords, currBox->basicBox.centerCoords);
787 ergo_real boxRadius = currBox->basicBox.width * 0.5 * template_blas_sqrt((ergo_real)3);
788 ergo_real requiredDistance = boxRadius + distrExtent;
789
790 // Note that the distance to the box multipole is different, since
791 // it is not necessarily placed at the box center.
792 ergo_real multipoleDistance = get_distance_3d(list[0].distr.centerCoords, boxMultipole->centerCoords);
793 int degreeNeeded = integralInfo.GetMMLimitTable().get_minimum_multipole_degree_needed(multipoleDistance,
794 boxMultipole,
795 maxDegreeForDistrs,
796 maxMomentVectorNormForDistrsList,
797 threshold);
798 if(degreeNeeded < 0)
799 return -1;
800 degreeNeeded+=2; // We need a couple of extra degrees to handle gradient computation. FIXME: IS THIS REALL NEEDED? WHY? MAYBE IT WAS BUGS IN OTHER PLACES THAT CAUSED THE PROBLEMS?
801
802 bool multipoleDegreeIsSafe = false;
803 // Demand at least two degrees margin compared to
804 // MAX_MULTIPOLE_DEGREE, otherwise this may fail due to alternating
805 // odd/even degrees where only one of them are significant. This has
806 // happened in some test cases.
807 if(degreeNeeded <= MAX_MULTIPOLE_DEGREE-2)
808 multipoleDegreeIsSafe = true;
809
810 if(distance > requiredDistance && multipoleDegreeIsSafe)
811 {
812 // OK, use multipole description of atom charges.
813 int boxNeededNoOfMoments = (degreeNeeded+1)*(degreeNeeded+1);
814 // create interaction matrix
815 ergo_real T[boxNeededNoOfMoments * maxNoOfMomentsForDistrs];
816 ergo_real dx = boxMultipole->centerCoords[0] - list[0].distr.centerCoords[0];
817 ergo_real dy = boxMultipole->centerCoords[1] - list[0].distr.centerCoords[1];
818 ergo_real dz = boxMultipole->centerCoords[2] - list[0].distr.centerCoords[2];
819
820 interactor.getInteractionMatrix(dx, dy, dz, maxDegreeForDistrs, degreeNeeded, T);
821
822 ergo_real tempVector[MAX_NO_OF_MOMENTS_PER_MULTIPOLE];
823 for(int A = 0; A < maxNoOfMomentsForDistrs; A++)
824 {
825 ergo_real sum = 0;
826 for(int B = 0; B < boxNeededNoOfMoments; B++)
827 {
828 ergo_real momB = boxMultipole->momentList[B];
829 ergo_real Telement = T[A*boxNeededNoOfMoments+B];
830 sum += momB * Telement;
831 }
832 tempVector[A] = sum;
833 }
834
835 for(int i = 0; i < nDistrs; i++)
836 {
837 ergo_real sum = 0;
838 for(int A = 0; A < multipoleList[i].noOfMoments; A++)
839 sum += tempVector[A] * multipoleList[i].momentList[A];
840 int idx1 = list[i].basisFuncIdx1;
841 int idx2 = list[i].basisFuncIdx2;
842 if(ergo_CSR_add_to_element(V_CSR, idx1, idx2, -1 * sum) != 0)
843 return -1;
844 }
845
846 return 0;
847 }
848
849 // No, multipole description could not be used in this case.
850 if(currLevel == numberOfLevels-1)
851 {
852 // We are at top level, must compute explicit interactions
853
854 int Nmax = 0;
855 for(int i = 0; i < nDistrs; i++)
856 {
857 const DistributionSpecStruct* distr = &list[i].distr;
858 int N = distr->monomialInts[0] + distr->monomialInts[1] + distr->monomialInts[2];
859 if(N > Nmax)
860 Nmax = N;
861 }
862
863 const JK::ExchWeights CAM_params_not_used;
864
865 const Atom* currAtomList = &atomListReordered[currBox->basicBox.firstItemIndex];
866 int noOfAtomsCurrBox = currBox->basicBox.noOfItems;
867 for(int ia = 0; ia < noOfAtomsCurrBox; ia++)
868 {
869 const Atom* currAtom = &currAtomList[ia];
870
871 const DistributionSpecStruct* distr1 = &list[0].distr;
872 int noOfMonomials_1 = integralInfo.monomial_info.no_of_monomials_list[Nmax];
873 int noOfMonomials_2 = integralInfo.monomial_info.no_of_monomials_list[0];
874 ergo_real primitiveIntegralList_h[noOfMonomials_1*noOfMonomials_2];
875 ergo_real primitiveIntegralList_2[noOfMonomials_1*noOfMonomials_2];
876 // Let the distr be 1 and the pointcharge 2
877 ergo_real alpha1 = distr1->exponent;
878 ergo_real alpha0 = alpha1;
879 int n1 = Nmax;
880 int n2 = 0;
881 ergo_real dx0 = currAtom->coords[0] - distr1->centerCoords[0];
882 ergo_real dx1 = currAtom->coords[1] - distr1->centerCoords[1];
883 ergo_real dx2 = currAtom->coords[2] - distr1->centerCoords[2];
884 ergo_real resultPreFactor = 2 * pi / alpha1;
885 get_related_integrals_hermite(integralInfo,
886 CAM_params_not_used,
887 n1, noOfMonomials_1,
888 n2, noOfMonomials_2,
889 dx0,
890 dx1,
891 dx2,
892 alpha0,
893 resultPreFactor,
894 primitiveIntegralList_h);
895 integralInfo.multiply_by_hermite_conversion_matrix_from_right(n1,
896 n2,
897 1.0/alpha1,
898 primitiveIntegralList_h,
899 primitiveIntegralList_2);
900 for(int id = 0; id < nDistrs; id++)
901 {
902 const DistributionSpecStruct* distr = &list[id].distr;
903 int n1x = distr->monomialInts[0];
904 int n1y = distr->monomialInts[1];
905 int n1z = distr->monomialInts[2];
906 int monomialIndex = integralInfo.monomial_info.monomial_index_list[n1x][n1y][n1z];
907 ergo_real integralValue = currAtom->charge * distr->coeff * primitiveIntegralList_2[monomialIndex];
908 int idx1 = list[id].basisFuncIdx1;
909 int idx2 = list[id].basisFuncIdx2;
910 if(ergo_CSR_add_to_element(V_CSR, idx1, idx2, -1 * integralValue) != 0)
911 return -1;
912 }
913
914 } // END FOR ia
915
916 return 0;
917 }
918
919 // Go to next level
920 int noOfChildren = currBox->basicBox.noOfChildBoxes;
921 for(int i = 0; i < noOfChildren; i++)
922 {
923 int childBoxIndex = currBox->basicBox.firstChildBoxIndex + i;
924 if(do_interaction_recursive_2(integralInfo,
925 V_CSR,
926 noOfBasisFuncIndexPairs,
927 basisFuncIndexPairList,
928 list,
929 nDistrs,
930 multipoleList,
931 maxMomentVectorNormForDistrsList,
932 maxNoOfMomentsForDistrs,
933 maxDegreeForDistrs,
934 distrExtent,
935 atomListReordered,
936 atomPermutation,
937 threshold,
938 boxList,
939 interactor,
940 childBoxIndex,
941 currLevel + 1,
942 numberOfLevels) != 0)
943 return -1;
944 } // END FOR i
945 return 0;
946 }
947
948
949
950
951 static int
get_list_of_distrs_for_V(const BasisInfoStruct & basisInfo,const basis_func_index_pair_struct_1el * basisFuncIndexPairList,int noOfBasisFuncIndexPairs,ergo_real threshold,ergo_real maxCharge,DistributionSpecStructWithIndexes * resultList,int maxCountResult)952 get_list_of_distrs_for_V(const BasisInfoStruct& basisInfo,
953 const basis_func_index_pair_struct_1el* basisFuncIndexPairList,
954 int noOfBasisFuncIndexPairs,
955 ergo_real threshold,
956 ergo_real maxCharge,
957 DistributionSpecStructWithIndexes* resultList,
958 int maxCountResult)
959 {
960 int distrCount = 0;
961 for(int kk = 0; kk < noOfBasisFuncIndexPairs; kk++)
962 {
963 int i = basisFuncIndexPairList[kk].index_1;
964 int j = basisFuncIndexPairList[kk].index_2;
965 const int maxCountProduct = POLY_PRODUCT_MAX_DISTRS;
966 DistributionSpecStruct psi_list[maxCountProduct];
967 /* form product of basisfuncs i and j, store product in psi_list */
968 int n_psi = get_product_simple_primitives(basisInfo, i,
969 basisInfo, j,
970 psi_list, maxCountProduct, 0);
971 if(n_psi < 0)
972 {
973 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in get_product_simple_primitives");
974 return -1;
975 }
976 for(int k = 0; k < n_psi; k++)
977 {
978 // now take care of psi_list[k]
979 // Here, we estimate the largest possible contribution to V from this distr as (maxCharge * 2 * pi * coeff / exponent).
980 DistributionSpecStruct* prim = &psi_list[k];
981 ergo_real maxContrib = template_blas_fabs(maxCharge * 2 * pi * prim->coeff / prim->exponent);
982 if(maxContrib > threshold)
983 {
984 if(maxCountResult > 0 && distrCount >= maxCountResult)
985 {
986 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in get_list_of_distrs_for_V: (maxCountResult > 0 && distrCount >= maxCountResult).");
987 return -1;
988 }
989 if(resultList != NULL)
990 {
991 resultList[distrCount].distr = *prim;
992 resultList[distrCount].pairIdx = kk;
993 resultList[distrCount].basisFuncIdx1 = i;
994 resultList[distrCount].basisFuncIdx2 = j;
995 }
996 distrCount++;
997 } // END IF above threshold
998 } // END FOR k
999 } // END FOR kk
1000 return distrCount;
1001 }
1002
1003
1004
1005 static ergo_real
get_nucl_repulsion_energy_using_multipoles(const Atom * atomListReordered,const int * atomPermutation,ergo_real threshold,const atom_box_struct * boxList,MMInteractor & interactor,int boxIndex1,int boxIndex2,int currLevel,int numberOfLevels)1006 get_nucl_repulsion_energy_using_multipoles(const Atom *atomListReordered,
1007 const int* atomPermutation, // list of int, length=nAtoms, saying how atoms have been reordered in return_atomListReordered.
1008 ergo_real threshold,
1009 const atom_box_struct* boxList,
1010 MMInteractor & interactor,
1011 int boxIndex1,
1012 int boxIndex2,
1013 int currLevel,
1014 int numberOfLevels) {
1015 const atom_box_struct* currBox1 = &boxList[boxIndex1];
1016 const atom_box_struct* currBox2 = &boxList[boxIndex2];
1017 const multipole_struct_large* boxMultipole1 = &currBox1->multipole;
1018 const multipole_struct_large* boxMultipole2 = &currBox2->multipole;
1019 // Check if we are at level of smallest boxes
1020 if(currLevel == numberOfLevels-1) {
1021 // Do interaction explicitly
1022 const Atom* currAtomList1 = &atomListReordered[currBox1->basicBox.firstItemIndex];
1023 int noOfAtomsCurrBox1 = currBox1->basicBox.noOfItems;
1024 const Atom* currAtomList2 = &atomListReordered[currBox2->basicBox.firstItemIndex];
1025 int noOfAtomsCurrBox2 = currBox2->basicBox.noOfItems;
1026 ergo_real sum = 0;
1027 for(int i1 = 0; i1 < noOfAtomsCurrBox1; i1++)
1028 for(int i2 = 0; i2 < noOfAtomsCurrBox2; i2++) {
1029 if(boxIndex1 == boxIndex2 && i1 == i2)
1030 continue; // skip interaction of an atom with itself
1031 if(boxIndex1 == boxIndex2 && i1 > i2)
1032 continue; // do not double-count interactions
1033 const Atom* currAtom1 = &currAtomList1[i1];
1034 const Atom* currAtom2 = &currAtomList2[i2];
1035 ergo_real dx0 = currAtom1->coords[0] - currAtom2->coords[0];
1036 ergo_real dx1 = currAtom1->coords[1] - currAtom2->coords[1];
1037 ergo_real dx2 = currAtom1->coords[2] - currAtom2->coords[2];
1038 ergo_real distance = template_blas_sqrt(dx0*dx0+dx1*dx1+dx2*dx2);
1039 sum += currAtom1->charge * currAtom2->charge / distance;
1040 }
1041 return sum;
1042 }
1043 // Check if multipole representation can be used
1044 // check if boxes are far enough apart so that we can consider using multipole description.
1045 ergo_real distance_between_box_centers = get_distance_3d(currBox1->basicBox.centerCoords, currBox2->basicBox.centerCoords);
1046 ergo_real boxRadius = currBox1->basicBox.width * 0.5 * template_blas_sqrt((ergo_real)3);
1047 ergo_real requiredDistance = 2.2*boxRadius; // 2*boxRadius should be enough, use 2.2*boxRadius to have some margin
1048 if(distance_between_box_centers > requiredDistance) {
1049 // Try using multipole description of atom charges.
1050 // Try with different multipole degree and check the difference, if nearly same result is obtained we trust it, otherwise not.
1051 const int nDegreesToTry = 3;
1052 ergo_real resultMin = 0;
1053 ergo_real resultMax = 0;
1054 const int highest_degree_to_use = MAX_MULTIPOLE_DEGREE / 2; // different choices possible here, pick something that gives good performance
1055 ergo_real dx = boxMultipole2->centerCoords[0] - boxMultipole1->centerCoords[0];
1056 ergo_real dy = boxMultipole2->centerCoords[1] - boxMultipole1->centerCoords[1];
1057 ergo_real dz = boxMultipole2->centerCoords[2] - boxMultipole1->centerCoords[2];
1058 int boxNeededNoOfMomentsMax = (highest_degree_to_use+1)*(highest_degree_to_use+1);
1059 ergo_real T[boxNeededNoOfMomentsMax * boxNeededNoOfMomentsMax];
1060 interactor.getInteractionMatrix(dx, dy, dz, highest_degree_to_use, highest_degree_to_use, T);
1061 for(int i = 0; i < nDegreesToTry; i++) {
1062 int degree = highest_degree_to_use - i;
1063 if(degree < 0)
1064 throw std::runtime_error("Error in get_nucl_repulsion_energy_using_multipoles: (degree < 0).");
1065 int boxNeededNoOfMoments = (degree+1)*(degree+1);
1066 // create interaction matrix
1067 ergo_real sum = 0;
1068 for(int A = 0; A < boxNeededNoOfMoments; A++)
1069 for(int B = 0; B < boxNeededNoOfMoments; B++)
1070 sum += T[A*boxNeededNoOfMomentsMax+B] * boxMultipole1->momentList[A] * boxMultipole2->momentList[B];
1071 if(i == 0) {
1072 resultMin = sum;
1073 resultMax = sum;
1074 }
1075 if(sum < resultMin)
1076 resultMin = sum;
1077 if(sum > resultMax)
1078 resultMax = sum;
1079 } // end for i
1080 ergo_real diff = resultMax - resultMin;
1081 if(diff < threshold)
1082 return 0.5*(resultMin+resultMax);
1083 } // end if distance large enough to consider multipoles
1084 // Go to next level, smaller boxes
1085 ergo_real sum = 0;
1086 int noOfChildren1 = currBox1->basicBox.noOfChildBoxes;
1087 int noOfChildren2 = currBox2->basicBox.noOfChildBoxes;
1088 for(int i1 = 0; i1 < noOfChildren1; i1++)
1089 for(int i2 = 0; i2 < noOfChildren2; i2++) {
1090 if(boxIndex1 == boxIndex2 && i1 > i2)
1091 continue; // do not double-count interactions
1092 int childBoxIndex1 = currBox1->basicBox.firstChildBoxIndex + i1;
1093 int childBoxIndex2 = currBox2->basicBox.firstChildBoxIndex + i2;
1094 sum += get_nucl_repulsion_energy_using_multipoles(atomListReordered,
1095 atomPermutation, // list of int, length=nAtoms, saying how atoms have been reordered in return_atomListReordered.
1096 threshold,
1097 boxList,
1098 interactor,
1099 childBoxIndex1,
1100 childBoxIndex2,
1101 currLevel + 1,
1102 numberOfLevels);
1103 } // END FOR i1 i2
1104 return sum;
1105 }
1106
1107
compute_V_and_gradient_linear(const BasisInfoStruct & basisInfo,const IntegralInfo & integralInfo,const Molecule & molecule,ergo_real threshold,ergo_real boxSize,const basis_func_index_pair_struct_1el * basisFuncIndexPairList,ergo_real * V_list,int noOfBasisFuncIndexPairs,bool compute_gradient_also,const ergo_real * D_list,ergo_real * gradient_list,ergo_real & result_nuclearRepulsionEnergy)1108 int compute_V_and_gradient_linear(const BasisInfoStruct& basisInfo,
1109 const IntegralInfo& integralInfo,
1110 const Molecule& molecule,
1111 ergo_real threshold,
1112 ergo_real boxSize,
1113 const basis_func_index_pair_struct_1el* basisFuncIndexPairList,
1114 ergo_real* V_list,
1115 int noOfBasisFuncIndexPairs,
1116 bool compute_gradient_also,
1117 const ergo_real* D_list, // List of corresponding density matrix elemets; used for compute_gradient_also case, NULL otherwise
1118 ergo_real* gradient_list, // list of result gradient values; used for compute_gradient_also case, NULL otherwise
1119 ergo_real & result_nuclearRepulsionEnergy
1120 )
1121 {
1122 result_nuclearRepulsionEnergy = 0; // to be set later
1123 int errorCount = 0;
1124 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "compute_V_and_gradient_linear start, compute_gradient_also = %d.", compute_gradient_also);
1125
1126 Util::TimeMeter timeMeterTotal;
1127 Util::TimeMeter timeMeterInitPart;
1128
1129 // Get maxCharge
1130 ergo_real maxCharge = 0;
1131 for(int i = 0; i < molecule.getNoOfAtoms(); i++) {
1132 ergo_real currCharge = molecule.getAtom(i).charge;
1133 if(currCharge > maxCharge)
1134 maxCharge = currCharge;
1135 }
1136
1137 // Create list of distributions
1138 int nDistrs = get_list_of_distrs_for_V(basisInfo,
1139 basisFuncIndexPairList,
1140 noOfBasisFuncIndexPairs,
1141 threshold,
1142 maxCharge,
1143 NULL,
1144 0);
1145 if(nDistrs <= 0)
1146 {
1147 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in compute_V_and_gradient_linear: (nDistrs <= 0).");
1148 return -1;
1149 }
1150 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "compute_V_and_gradient_linear nDistrs = %9i", nDistrs);
1151 std::vector<DistributionSpecStructWithIndexes> list(nDistrs);
1152 if(get_list_of_distrs_for_V(basisInfo,
1153 basisFuncIndexPairList,
1154 noOfBasisFuncIndexPairs,
1155 threshold,
1156 maxCharge,
1157 &list[0],
1158 nDistrs) != nDistrs)
1159 {
1160 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in compute_V_and_gradient_linear, in get_list_of_distrs_for_V.");
1161 return -1;
1162 }
1163
1164 // Sort list of distrs by x, y, z, exponent.
1165 // The point of this is to group together distrs that have same center and same exponent.
1166 sort_distr_list(&list[0], nDistrs);
1167
1168 // identify groups of distrs that have same center and same exponent.
1169 // Allocate according to worst case, each distr being a separate group.
1170 std::vector<group_struct> groupList(nDistrs);
1171 int ind = 0;
1172 int currGroupInd = 0;
1173 int groupCount = 0;
1174 int maxNDistrsPerGroup = 0;
1175 while(ind < nDistrs)
1176 {
1177 ind++;
1178 if(ind < nDistrs)
1179 {
1180 if(compare_distrs<DistributionSpecStructWithIndexes>(&list[ind], &list[currGroupInd]) == 0)
1181 continue;
1182 }
1183 // define new group
1184 groupList[groupCount].startIndex = currGroupInd;
1185 groupList[groupCount].count = ind - currGroupInd;
1186 if (groupList[groupCount].count > maxNDistrsPerGroup)
1187 maxNDistrsPerGroup = groupList[groupCount].count;
1188 groupCount++;
1189 // start next group
1190 currGroupInd = ind;
1191 }
1192 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "compute_V_and_gradient_linear groupCount = %9i", groupCount);
1193
1194 // Note that boxList and atomListReordered are allocated by create_nuclei_mm_tree,
1195 // we must remember to free them in the end.
1196 BoxSystem boxSystem;
1197 std::vector<int> atomPermutation(molecule.getNoOfAtoms());
1198 atom_box_struct* boxList = NULL;
1199 Atom *atomListReordered = NULL;
1200 int numberOfLevels = -1;
1201 if(create_nuclei_mm_tree(integralInfo,
1202 molecule.getNoOfAtoms(), molecule.getAtomListPtr(), boxSize,
1203 boxSystem,
1204 &boxList, &numberOfLevels,
1205 &atomListReordered,
1206 &atomPermutation[0],
1207 compute_gradient_also) != 0)
1208 {
1209 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "create_nuclei_mm_tree failed");
1210 return -1;
1211 }
1212
1213 timeMeterInitPart.print(LOG_AREA_INTEGRALS, "compute_V_and_gradient_linear init part (including create_nuclei_mm_tree)");
1214
1215 // Compute nuclear repulsion energy using multipole tree.
1216 Util::TimeMeter timeMeterNuclRep;
1217 {
1218 MMInteractor interactor(integralInfo.GetMultipolePrep());
1219 int boxIndex1 = 0;
1220 int boxIndex2 = 0;
1221 int currLevel = 0;
1222 ergo_real nuclRepEnergy = get_nucl_repulsion_energy_using_multipoles(atomListReordered,
1223 &atomPermutation[0],
1224 threshold,
1225 boxList,
1226 interactor,
1227 boxIndex1,
1228 boxIndex2,
1229 currLevel,
1230 numberOfLevels);
1231 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "get_nucl_repulsion_energy_using_multipoles() gave nuclRepEnergy = %22.11f", nuclRepEnergy);
1232 result_nuclearRepulsionEnergy = nuclRepEnergy;
1233 }
1234 timeMeterNuclRep.print(LOG_AREA_INTEGRALS, "get_nucl_repulsion_energy_using_multipoles");
1235
1236 Util::TimeMeter timeMeterMainPart;
1237
1238 const JK::ExchWeights CAM_params_not_used;
1239
1240 #ifdef _OPENMP
1241 #pragma omp parallel
1242 #endif
1243 {
1244 MMInteractor interactor(integralInfo.GetMultipolePrep());
1245 ergo_real *private_V_list = V_list;
1246 ergo_real* private_gradient_list = gradient_list;
1247 multipole_struct_small* multipoleList =
1248 new multipole_struct_small[maxNDistrsPerGroup];
1249 #ifdef _OPENMP
1250 if (omp_get_thread_num() != 0) {
1251 private_V_list = new ergo_real[noOfBasisFuncIndexPairs];
1252 if(compute_gradient_also)
1253 private_gradient_list = new ergo_real[3*molecule.getNoOfAtoms()];
1254 }
1255 #endif
1256
1257 memset(private_V_list, 0, noOfBasisFuncIndexPairs*sizeof(ergo_real));
1258 if(compute_gradient_also)
1259 memset(private_gradient_list, 0, 3*molecule.getNoOfAtoms()*sizeof(ergo_real));
1260 #ifdef _OPENMP
1261 #pragma omp for
1262 #endif
1263 for(int groupIndex = 0; groupIndex < groupCount; groupIndex++)
1264 {
1265 DistributionSpecStructWithIndexes* currList = &list[groupList[groupIndex].startIndex];
1266 int nDistrsCurrGroup = groupList[groupIndex].count;
1267
1268 // Create multipoles for distrs in this group.
1269 memset(multipoleList, 0, nDistrsCurrGroup*sizeof(multipole_struct_small));
1270 for(int i = 0; i < nDistrsCurrGroup; i++)
1271 {
1272 compute_multipole_moments(integralInfo,
1273 &currList[i].distr,
1274 &multipoleList[i]);
1275 }
1276
1277 int maxNoOfMoments = 0;
1278 int maxDegree = 0;
1279 ergo_real maxMomentVectorNormForDistrsList[MAX_MULTIPOLE_DEGREE_BASIC+1];
1280 for(int l = 0; l <= MAX_MULTIPOLE_DEGREE_BASIC; l++)
1281 maxMomentVectorNormForDistrsList[l] = 0;
1282 for(int i = 0; i < nDistrsCurrGroup; i++)
1283 {
1284 if(multipoleList[i].noOfMoments > maxNoOfMoments)
1285 maxNoOfMoments = multipoleList[i].noOfMoments;
1286 if(multipoleList[i].degree > maxDegree)
1287 maxDegree = multipoleList[i].degree;
1288 const multipole_struct_small* distrMultipole = &multipoleList[i];
1289 for(int l = 0; l <= distrMultipole->degree; l++)
1290 {
1291 int startIndex = l*l;
1292 int endIndex = (l+1)*(l+1);
1293 ergo_real sum = 0;
1294 for(int A = startIndex; A < endIndex; A++)
1295 sum += distrMultipole->momentList[A]*distrMultipole->momentList[A];
1296 ergo_real subNorm = template_blas_sqrt(sum);
1297 if(subNorm > maxMomentVectorNormForDistrsList[l])
1298 maxMomentVectorNormForDistrsList[l] = subNorm;
1299 }
1300 }
1301
1302 // Here we use an extent such that beyond the extent the abs
1303 // value of any distr is smaller than threshold/maxCharge.
1304 ergo_real maxabscoeff = 0;
1305 for(int i = 0; i < nDistrsCurrGroup; i++)
1306 {
1307 ergo_real abscoeff = template_blas_fabs(currList[i].distr.coeff);
1308 if(abscoeff > maxabscoeff)
1309 maxabscoeff = abscoeff;
1310 }
1311 ergo_real exponent = currList[0].distr.exponent; // all exponents in group are equal anyway.
1312 ergo_real R2 = -1 * (1/exponent) * template_blas_log(threshold/(maxabscoeff*maxCharge));
1313 ergo_real extent = 0;
1314 if(R2 > 0) // R2 can become negative, e.g. if maxabscoeff is very small, in such cases we let extent be zero.
1315 extent = template_blas_sqrt(R2);
1316 // Take care of interaction of this group with MM tree
1317 if(do_interaction_recursive(integralInfo,
1318 private_V_list,
1319 noOfBasisFuncIndexPairs,
1320 basisFuncIndexPairList,
1321 currList,
1322 nDistrsCurrGroup,
1323 multipoleList,
1324 maxMomentVectorNormForDistrsList,
1325 maxNoOfMoments,
1326 maxDegree,
1327 extent,
1328 atomListReordered,
1329 &atomPermutation[0],
1330 threshold,
1331 boxList,
1332 interactor,
1333 0,
1334 0,
1335 numberOfLevels,
1336 compute_gradient_also,
1337 D_list,
1338 private_gradient_list) != 0)
1339 {
1340 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in do_interaction_recursive");
1341 #ifdef _OPENMP
1342 #pragma omp atomic
1343 #endif
1344 errorCount++;
1345 continue;
1346 }
1347 }
1348 #ifdef _OPENMP
1349 if(omp_get_thread_num() != 0) {
1350 /* collect data from different threads. */
1351 #pragma omp critical
1352 {
1353 for(int i=0; i<noOfBasisFuncIndexPairs; i++)
1354 V_list[i] += private_V_list[i];
1355 if(compute_gradient_also)
1356 for(int i = 0; i < 3*molecule.getNoOfAtoms(); i++)
1357 gradient_list[i] += private_gradient_list[i];
1358 }
1359 delete [] private_V_list;
1360 delete [] private_gradient_list;
1361 }
1362 #endif
1363 delete []multipoleList;
1364 }
1365
1366 timeMeterMainPart.print(LOG_AREA_INTEGRALS, "compute_V_and_gradient_linear main part");
1367
1368 if(compute_gradient_also) {
1369 Util::TimeMeter timeMeterGradientMultipolePart;
1370 // Update gradient with multipole contributions
1371 for(int boxIdx = 0; boxIdx < boxSystem.totNoOfBoxes; boxIdx++) {
1372 atom_box_struct* currBox = &boxList[boxIdx];
1373 int noOfAtomsCurrBox = currBox->basicBox.noOfItems;
1374 // Go through all atoms in this box
1375 for(int j = 0; j < noOfAtomsCurrBox; j++) {
1376 int atomIndex = atomPermutation[currBox->basicBox.firstItemIndex + j];
1377 for(int coordIdx = 0; coordIdx < 3; coordIdx++) {
1378 for(int k = 0; k < MAX_NO_OF_MOMENTS_PER_MULTIPOLE; k++) {
1379 ergo_real multipole_moment_derivative = currBox->multipole_moment_derivatives[j*MAX_NO_OF_MOMENTS_PER_MULTIPOLE*3 + MAX_NO_OF_MOMENTS_PER_MULTIPOLE*coordIdx + k];
1380 ergo_real derivative_wrt_multipole_moment = currBox->derivatives_wrt_multipole_moments[k];
1381 ergo_real gradient_contrib = multipole_moment_derivative * derivative_wrt_multipole_moment;
1382 gradient_list[atomIndex*3+coordIdx] += gradient_contrib;
1383 } // end for k
1384 } // end for coordIdx
1385 } // end for j
1386 } // end for boxIdx
1387 timeMeterGradientMultipolePart.print(LOG_AREA_INTEGRALS, "compute_V_and_gradient_linear gradient multipole contributions part");
1388
1389 // Update gradient with nuclear-nuclear energy contribution
1390 Util::TimeMeter timeMeterNuclearRepulsionEnergyGradientContrib;
1391 molecule.getNuclearRepulsionEnergyGradientContribQuadratic(gradient_list);
1392 timeMeterNuclearRepulsionEnergyGradientContrib.print(LOG_AREA_INTEGRALS, "compute_V_and_gradient_linear getNuclearRepulsionEnergyGradientContribQuadratic");
1393
1394 // Cleanup: delete buffers that were allocated during creation of multipole tree
1395 for(int boxIdx = 0; boxIdx < boxSystem.totNoOfBoxes; boxIdx++) {
1396 atom_box_struct* currBox = &boxList[boxIdx];
1397 delete [] currBox->multipole_moment_derivatives;
1398 currBox->multipole_moment_derivatives = NULL;
1399 delete [] currBox->derivatives_wrt_multipole_moments;
1400 currBox->derivatives_wrt_multipole_moments = NULL;
1401 }
1402
1403 } // end if compute_gradient_also
1404
1405 delete [] boxList;
1406 delete [] atomListReordered;
1407
1408 timeMeterTotal.print(LOG_AREA_INTEGRALS, "compute_V_and_gradient_linear total");
1409
1410 return -errorCount;
1411 }
1412
1413
1414
compute_V_hierarchical(const BasisInfoStruct & basisInfo,const IntegralInfo & integralInfo,const Molecule & molecule,ergo_real threshold,ergo_real boxSize,const basis_func_index_pair_struct_1el * basisFuncIndexPairList,int noOfBasisFuncIndexPairs,csr_matrix_struct * V_CSR,ergo_real & result_nuclearRepulsionEnergy)1415 int compute_V_hierarchical(const BasisInfoStruct& basisInfo,
1416 const IntegralInfo& integralInfo,
1417 const Molecule& molecule,
1418 ergo_real threshold,
1419 ergo_real boxSize,
1420 const basis_func_index_pair_struct_1el* basisFuncIndexPairList,
1421 int noOfBasisFuncIndexPairs,
1422 csr_matrix_struct* V_CSR,
1423 ergo_real & result_nuclearRepulsionEnergy
1424 ) {
1425 result_nuclearRepulsionEnergy = 0; // computed later
1426 bool compute_gradient_also = false;
1427 int errorCount = 0;
1428 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "compute_V_hierarchical start.");
1429
1430 Util::TimeMeter timeMeterTotal;
1431 Util::TimeMeter timeMeterInitPart;
1432
1433 // Get maxCharge
1434 ergo_real maxCharge = 0;
1435 for(int i = 0; i < molecule.getNoOfAtoms(); i++) {
1436 ergo_real currCharge = molecule.getAtom(i).charge;
1437 if(currCharge > maxCharge)
1438 maxCharge = currCharge;
1439 }
1440
1441 // Create list of distributions
1442 int nDistrs = get_list_of_distrs_for_V(basisInfo, basisFuncIndexPairList, noOfBasisFuncIndexPairs,
1443 threshold, maxCharge, NULL, 0);
1444 if(nDistrs <= 0) {
1445 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in compute_V_hierarchical: (nDistrs <= 0).");
1446 return -1;
1447 }
1448 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "compute_V_hierarchical nDistrs = %9i", nDistrs);
1449 std::vector<DistributionSpecStructWithIndexes> listTmp(nDistrs);
1450 if(get_list_of_distrs_for_V(basisInfo, basisFuncIndexPairList, noOfBasisFuncIndexPairs,
1451 threshold, maxCharge, &listTmp[0], nDistrs) != nDistrs) {
1452 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in compute_V_hierarchical, in get_list_of_distrs_for_V.");
1453 return -1;
1454 }
1455
1456 std::vector<DistributionSpecStructWithIndexes2> listTmp2(nDistrs);
1457 for(int i = 0; i < nDistrs; i++) {
1458 listTmp2[i].distr = listTmp[i].distr;
1459 listTmp2[i].basisFuncIdx1 = listTmp[i].basisFuncIdx1;
1460 listTmp2[i].basisFuncIdx2 = listTmp[i].basisFuncIdx2;
1461 }
1462 listTmp.clear();
1463
1464 SetOfDistrsForV setOfDistrsForV;
1465 organize_distrs_for_V(integralInfo, setOfDistrsForV, listTmp2, threshold, maxCharge);
1466 listTmp2.clear();
1467
1468 // Note that boxList and atomListReordered are allocated by create_nuclei_mm_tree,
1469 // we must remember to free them in the end.
1470 BoxSystem boxSystem;
1471 std::vector<int> atomPermutation(molecule.getNoOfAtoms());
1472 atom_box_struct* boxList = NULL;
1473 Atom *atomListReordered = NULL;
1474 int numberOfLevels = -1;
1475 if(create_nuclei_mm_tree(integralInfo,
1476 molecule.getNoOfAtoms(), molecule.getAtomListPtr(), boxSize,
1477 boxSystem, &boxList, &numberOfLevels, &atomListReordered,
1478 &atomPermutation[0], compute_gradient_also) != 0) {
1479 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "create_nuclei_mm_tree failed");
1480 return -1;
1481 }
1482
1483 timeMeterInitPart.print(LOG_AREA_INTEGRALS, "compute_V_hierarchical init part (including create_nuclei_mm_tree)");
1484
1485 // Compute nuclear repulsion energy using multipole tree.
1486 Util::TimeMeter timeMeterNuclRep;
1487 {
1488 MMInteractor interactor(integralInfo.GetMultipolePrep());
1489 int boxIndex1 = 0;
1490 int boxIndex2 = 0;
1491 int currLevel = 0;
1492 ergo_real nuclRepEnergy = get_nucl_repulsion_energy_using_multipoles(atomListReordered,
1493 &atomPermutation[0],
1494 threshold,
1495 boxList,
1496 interactor,
1497 boxIndex1,
1498 boxIndex2,
1499 currLevel,
1500 numberOfLevels);
1501 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "get_nucl_repulsion_energy_using_multipoles() gave nuclRepEnergy = %22.11f", nuclRepEnergy);
1502 result_nuclearRepulsionEnergy = nuclRepEnergy;
1503 }
1504 timeMeterNuclRep.print(LOG_AREA_INTEGRALS, "get_nucl_repulsion_energy_using_multipoles");
1505
1506 Util::TimeMeter timeMeterMainPart;
1507
1508 const JK::ExchWeights CAM_params_not_used;
1509 MMInteractor interactor(integralInfo.GetMultipolePrep());
1510 int groupCount = setOfDistrsForV.groupList.size();
1511 for(int groupIndex = 0; groupIndex < groupCount; groupIndex++) {
1512 int groupStartIdx = setOfDistrsForV.groupList[groupIndex].startIndex;
1513 DistributionSpecStructWithIndexes2* currList = &setOfDistrsForV.distrList[groupStartIdx];
1514 int nDistrsCurrGroup = setOfDistrsForV.groupList[groupIndex].count;
1515 ergo_real extent = setOfDistrsForV.groupList[groupIndex].maxExtent;
1516 const multipole_struct_small* multipoleList = &setOfDistrsForV.multipoleList[groupStartIdx];
1517 const ergo_real* maxMomentVectorNormForDistrsList = setOfDistrsForV.maxMomentVectorNormList[groupIndex].maxMomentVectorNormList;
1518 int maxNoOfMoments = setOfDistrsForV.groupList[groupIndex].maxNoOfMoments;
1519 int maxDegree = setOfDistrsForV.groupList[groupIndex].maxDegree;
1520 // Take care of interaction of this group with MM tree
1521 if(do_interaction_recursive_2(integralInfo, V_CSR, noOfBasisFuncIndexPairs,
1522 basisFuncIndexPairList, currList, nDistrsCurrGroup,
1523 multipoleList, maxMomentVectorNormForDistrsList, maxNoOfMoments,
1524 maxDegree, extent, atomListReordered,
1525 &atomPermutation[0], threshold, boxList,
1526 interactor, 0, 0, numberOfLevels) != 0) {
1527 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in do_interaction_recursive");
1528 errorCount++;
1529 continue;
1530 }
1531 }
1532
1533 timeMeterMainPart.print(LOG_AREA_INTEGRALS, "compute_V_hierarchical main part");
1534
1535 delete [] boxList;
1536 delete [] atomListReordered;
1537
1538 timeMeterTotal.print(LOG_AREA_INTEGRALS, "compute_V_hierarchical total");
1539
1540 return -errorCount;
1541 }
1542
1543
1544 static ergo_real
simplePrimVintegralSingle(const DistributionSpecStruct & prim,const Atom & atom,const IntegralInfo & integralInfo)1545 simplePrimVintegralSingle(const DistributionSpecStruct & prim,
1546 const Atom & atom,
1547 const IntegralInfo & integralInfo) {
1548 return do_1e_repulsion_integral_using_symb_info(prim,
1549 atom.charge,
1550 atom.coords,
1551 integralInfo);
1552 }
1553
1554
1555 ergo_real
simplePrimVintegral_list(const DistributionSpecStruct * list,int nPrims,const Atom & atom,const IntegralInfo & integralInfo)1556 simplePrimVintegral_list(const DistributionSpecStruct* list,
1557 int nPrims,
1558 const Atom & atom,
1559 const IntegralInfo & integralInfo) {
1560 ergo_real sum = 0;
1561 for(int k = 0; k < nPrims; k++) {
1562 const DistributionSpecStruct & currDistr = list[k];
1563 sum += simplePrimVintegralSingle(currDistr, atom, integralInfo);
1564 }
1565 return sum;
1566 }
1567
1568
1569 int
compute_V_matrix_full(const BasisInfoStruct & basisInfo,const IntegralInfo & integralInfo,int nAtoms,const Atom * atomList,ergo_real threshold,ergo_real * result)1570 compute_V_matrix_full(const BasisInfoStruct& basisInfo,
1571 const IntegralInfo& integralInfo,
1572 int nAtoms,
1573 const Atom* atomList,
1574 ergo_real threshold,
1575 ergo_real* result)
1576 {
1577 int mu, nu, A, j, k, nbast;
1578 nbast = basisInfo.noOfBasisFuncs;
1579
1580 for(mu = 0; mu < nbast; mu++)
1581 {
1582 BasisFuncStruct* basisFunc_mu = &basisInfo.basisFuncList[mu];
1583 int n_mu = basisFunc_mu->noOfSimplePrimitives;
1584 int start_prim_mu = basisFunc_mu->simplePrimitiveIndex;
1585 DistributionSpecStruct* list_mu =
1586 &basisInfo.simplePrimitiveList[start_prim_mu];
1587 for(nu = 0; nu <= mu; nu++)
1588 {
1589 BasisFuncStruct* basisFunc_nu = &basisInfo.basisFuncList[nu];
1590 int n_nu = basisFunc_nu->noOfSimplePrimitives;
1591 int start_prim_nu = basisFunc_nu->simplePrimitiveIndex;
1592 DistributionSpecStruct* list_nu =
1593 &basisInfo.simplePrimitiveList[start_prim_nu];
1594 /* compute matrix element [mu,nu] */
1595 ergo_real sum = 0;
1596 for(j = 0; j < n_mu; j++)
1597 {
1598 DistributionSpecStruct & prim_mu_j = list_mu[j];
1599 for(k = 0; k < n_nu; k++)
1600 {
1601 DistributionSpecStruct & prim_nu_k = list_nu[k];
1602 const int maxDistrsInTempList = 888;
1603 DistributionSpecStruct tempList[maxDistrsInTempList];
1604 int nNewPrims = get_product_simple_prims(prim_mu_j,
1605 prim_nu_k,
1606 tempList,
1607 maxDistrsInTempList,
1608 threshold);
1609 if(nNewPrims < 0)
1610 {
1611 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in get_product_simple_prims");
1612 exit(EXIT_FAILURE);
1613 }
1614 /* now loop over atoms */
1615 for(A = 0; A < nAtoms; A++)
1616 {
1617 sum += simplePrimVintegral_list(tempList,
1618 nNewPrims,
1619 atomList[A],
1620 integralInfo);
1621 } /* END FOR A */
1622 } /* END FOR k */
1623 } /* END FOR j */
1624 result[mu*nbast+nu] = -1 * sum;
1625 } /* END FOR nu */
1626 } /* END FOR mu */
1627
1628 // copy values to the other triangle
1629 for(mu = 0; mu < nbast; mu++)
1630 for(nu = mu+1; nu < nbast; nu++)
1631 result[mu*nbast+nu] = result[nu*nbast+mu];
1632
1633 return 0;
1634 }
1635
1636