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 organize_distrs_mm.cc
31
32 @brief Code for organizing a given set of primitive Gaussian
33 distributions (typically coming from basis function products)
34 regarding information related to multipole methods.
35
36 @author: Elias Rudberg <em>responsible</em>
37 */
38
39 #include "organize_distrs_mm.h"
40 #include "serialization_tools.h"
41 #include <stdexcept>
42
43 /* distr_org_mm_struct functions */
44
Data()45 distr_org_mm_struct::Data::Data() : chargeSum(0) {
46 multipole.degree = -1;
47 multipole.noOfMoments = 0;
48 memset(multipole.momentList, 0, MAX_NO_OF_MOMENTS_PER_MULTIPOLE*sizeof(ergo_real));
49 memset(&multipolePoint, 0x00, 3*sizeof(ergo_real));
50 memset(maxMomentVectorNormForDistrsList, 0, (MAX_MULTIPOLE_DEGREE_BASIC+1)*sizeof(ergo_real));
51 }
52
writeToBuffer(char * dataBuffer,size_t const bufferSize) const53 void distr_org_mm_struct::writeToBuffer(char* dataBuffer, size_t const bufferSize) const {
54 assert(bufferSize >= getSize());
55 char* p = dataBuffer;
56 memcpy(p, &data, sizeof(data));
57 p += sizeof(data);
58 std_vector_writeToBuffer_and_move_ptr(multipoleListForGroups, p);
59 std_vector_writeToBuffer_and_move_ptr(multipoleListForDistrs, p);
60 }
61
getSize() const62 size_t distr_org_mm_struct::getSize() const {
63 size_t size = sizeof(distr_org_mm_struct::Data);
64 size += std_vector_getSize(multipoleListForGroups);
65 size += std_vector_getSize(multipoleListForDistrs);
66 return size;
67 }
68
assignFromBuffer(char const * dataBuffer,size_t const bufferSize)69 void distr_org_mm_struct::assignFromBuffer(char const * dataBuffer, size_t const bufferSize) {
70 const char* p = dataBuffer;
71 size_t remainingBytes = bufferSize;
72 assert(remainingBytes >= sizeof(data));
73 memcpy(&data, p, sizeof(data));
74 p += sizeof(data);
75 const char* bufEndPtr = &dataBuffer[bufferSize];
76 std_vector_assignFromBuffer_and_move_ptr(multipoleListForGroups, p, bufEndPtr);
77 std_vector_assignFromBuffer_and_move_ptr(multipoleListForDistrs, p, bufEndPtr);
78 }
79
80 /* distr_list_description_struct functions */
81
writeToBuffer(char * dataBuffer,size_t const bufferSize) const82 void distr_list_description_struct::writeToBuffer(char* dataBuffer, size_t const bufferSize) const {
83 assert(bufferSize >= getSize());
84 char* p = dataBuffer;
85 org.writeToBuffer(p, org.getSize());
86 p += org.getSize();
87 org_mm.writeToBuffer(p, org_mm.getSize());
88 p += org_mm.getSize();
89 assert((size_t)(p - dataBuffer) <= bufferSize);
90 }
91
getSize() const92 size_t distr_list_description_struct::getSize() const {
93 return org.getSize() + org_mm.getSize();
94 }
95
assignFromBuffer(char const * dataBuffer,size_t const bufferSize)96 void distr_list_description_struct::assignFromBuffer(char const * dataBuffer, size_t const bufferSize) {
97 const char* p = dataBuffer;
98 org.assignFromBuffer(p, bufferSize);
99 p += org.getSize();
100 org_mm.assignFromBuffer(p, bufferSize - org.getSize());
101 p += org_mm.getSize();
102 assert((size_t)(p - dataBuffer) == bufferSize);
103 }
104
105 /* ************************************************** */
106
107 int
generate_multipoles_for_groups(const IntegralInfo & integralInfo,const distr_org_struct & org,distr_org_mm_struct & result_org_mm,ergo_real * averagePosList,int & avgPosCounter)108 generate_multipoles_for_groups(const IntegralInfo & integralInfo,
109 const distr_org_struct & org,
110 distr_org_mm_struct & result_org_mm,
111 ergo_real* averagePosList,
112 int & avgPosCounter
113 ) {
114 ergo_real chargeSum = 0;
115
116 const std::vector<batch_struct> & batchList = org.batchList;
117 const std::vector<cluster_struct> & clusterList = org.clusterList;
118 const std::vector<distr_group_struct> & groupList = org.groupList;
119 const std::vector<minimal_distr_struct> & minimalDistrList = org.minimalDistrList;
120 int batchCount = batchList.size();
121 const std::vector<basis_func_pair_struct> & basisFuncPairList = org.basisFuncPairList;
122
123 ergo_real* maxMomentVectorNormForDistrsList = result_org_mm.data.maxMomentVectorNormForDistrsList;
124 for(int l = 0; l <= MAX_MULTIPOLE_DEGREE_BASIC; l++)
125 maxMomentVectorNormForDistrsList[l] = 0;
126
127 int groupCount = org.groupList.size();
128 result_org_mm.multipoleListForGroups.resize(groupCount);
129
130 for(int batchIndex = 0; batchIndex < batchCount; batchIndex++) {
131 int clusterCount = batchList[batchIndex].noOfClusters;
132 int cluster_start = batchList[batchIndex].clusterStartIndex;
133 for(int clusterIndex = cluster_start; clusterIndex < cluster_start + clusterCount; clusterIndex++) {
134 int group_start = clusterList[clusterIndex].groupStartIndex;
135 int group_end = group_start + clusterList[clusterIndex].noOfGroups;
136 for(int groupIndex = group_start; groupIndex < group_end; groupIndex++) {
137 const distr_group_struct* currGroup = &groupList[groupIndex];
138
139 // Now create a single multipole description of the density of this group.
140 multipole_struct_small* multipoleCurrGroup = &result_org_mm.multipoleListForGroups[groupIndex];
141 multipoleCurrGroup->degree = -1;
142 multipoleCurrGroup->noOfMoments = 0;
143 multipoleCurrGroup->centerCoords[0] = currGroup->centerCoords[0];
144 multipoleCurrGroup->centerCoords[1] = currGroup->centerCoords[1];
145 multipoleCurrGroup->centerCoords[2] = currGroup->centerCoords[2];
146 memset(multipoleCurrGroup->momentList, 0, MAX_NO_OF_MOMENTS_PER_MULTIPOLE_BASIC*sizeof(ergo_real));
147
148 int distr_start = currGroup->startIndex;
149 int distr_end = distr_start + currGroup->distrCount;
150 for(int distrIndex = distr_start; distrIndex < distr_end; distrIndex++) {
151 int basisFuncPairIndex = minimalDistrList[distrIndex].basisFuncPairIndex;
152 int monomialIndex = minimalDistrList[distrIndex].monomialIndex;
153 ergo_real coeff = minimalDistrList[distrIndex].coeff;
154 // get monomialInts from monomialIndex
155 DistributionSpecStruct distr;
156 distr.monomialInts[0] = integralInfo.monomial_info.monomial_list[monomialIndex].ix;
157 distr.monomialInts[1] = integralInfo.monomial_info.monomial_list[monomialIndex].iy;
158 distr.monomialInts[2] = integralInfo.monomial_info.monomial_list[monomialIndex].iz;
159 distr.coeff = coeff;
160 distr.exponent = currGroup->exponent;
161 distr.centerCoords[0] = currGroup->centerCoords[0];
162 distr.centerCoords[1] = currGroup->centerCoords[1];
163 distr.centerCoords[2] = currGroup->centerCoords[2];
164
165 multipole_struct_small multipole;
166 if(compute_multipole_moments(integralInfo, &distr, &multipole) != 0) {
167 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in compute_multipole_moments");
168 return -1;
169 }
170
171 // add this multipole to multipole for group.
172 int a = basisFuncPairList[batchList[batchIndex].basisFuncPairListIndex+basisFuncPairIndex].index_1;
173 int b = basisFuncPairList[batchList[batchIndex].basisFuncPairListIndex+basisFuncPairIndex].index_2;
174 ergo_real factor = basisFuncPairList[batchList[batchIndex].basisFuncPairListIndex+basisFuncPairIndex].dmatElement;
175 if(a != b)
176 factor *= 2;
177
178 for(int l = 0; l <= multipole.degree; l++) {
179 int startIndex = l*l;
180 int endIndex = (l+1)*(l+1);
181 ergo_real sum = 0;
182 for(int A = startIndex; A < endIndex; A++)
183 sum += multipole.momentList[A]*multipole.momentList[A];
184 ergo_real subNorm = template_blas_sqrt(sum);
185 if(subNorm > maxMomentVectorNormForDistrsList[l])
186 maxMomentVectorNormForDistrsList[l] = subNorm;
187 }
188
189 for(int kk = 0; kk < multipole.noOfMoments; kk++)
190 multipoleCurrGroup->momentList[kk] += factor * multipole.momentList[kk];
191 if(multipole.degree > multipoleCurrGroup->degree)
192 multipoleCurrGroup->degree = multipole.degree;
193 if(multipole.noOfMoments > multipoleCurrGroup->noOfMoments)
194 multipoleCurrGroup->noOfMoments = multipole.noOfMoments;
195 } // END FOR distrIndex
196
197 // OK, multipoleCurrGroup is complete.
198 chargeSum += multipoleCurrGroup->momentList[0];
199 for(int kk = 0; kk < 3; kk++)
200 averagePosList[kk] += multipoleCurrGroup->centerCoords[kk];
201 avgPosCounter++;
202
203 } // END FOR groupIndex
204 } // END FOR clusterIndex
205 } // END FOR batchIndex
206
207 return 0;
208 }
209
210
211 int
get_multipole_pt_for_box(const ergo_real * boxCenterCoords,ergo_real boxWidth,const ergo_real * averagePosList,int avgPosCounter,ergo_real * resultMultipolePoint)212 get_multipole_pt_for_box(const ergo_real* boxCenterCoords,
213 ergo_real boxWidth,
214 const ergo_real* averagePosList,
215 int avgPosCounter,
216 ergo_real* resultMultipolePoint) {
217 // use average position instead of center-of-charge,
218 // because center-of-charge is ill-defined when some charges are negative.
219 if(avgPosCounter == 0) {
220 for(int kk = 0; kk < 3; kk++)
221 resultMultipolePoint[kk] = boxCenterCoords[kk];
222 }
223 else {
224 for(int kk = 0; kk < 3; kk++)
225 resultMultipolePoint[kk] = averagePosList[kk] / avgPosCounter;
226 }
227 // check that "resultMultipolePoint" is not too far from box center.
228 ergo_real sumofsquares = 0;
229 for(int kk = 0; kk < 3; kk++) {
230 ergo_real dx = resultMultipolePoint[kk] - boxCenterCoords[kk];
231 sumofsquares += dx*dx;
232 }
233 ergo_real distFromCenter = template_blas_sqrt(sumofsquares);
234 if(distFromCenter > boxWidth) {
235 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in get_multipole_pt_for_box: (distFromCenter > boxWidth).");
236 return -1;
237 }
238 return 0;
239 }
240
241
242 int
translate_multipoles_for_box(distr_org_mm_struct & result_org_mm,const distr_org_struct & org,const MMTranslator & translator)243 translate_multipoles_for_box(distr_org_mm_struct & result_org_mm,
244 const distr_org_struct & org,
245 const MMTranslator & translator
246 ) {
247
248 ergo_real* multipolePointCoords = result_org_mm.data.multipolePoint;
249 multipole_struct_large branchMultipole;
250 for(int A = 0; A < MAX_NO_OF_MOMENTS_PER_MULTIPOLE; A++)
251 branchMultipole.momentList[A] = 0;
252 for(int kk = 0; kk < 3; kk++)
253 branchMultipole.centerCoords[kk] = multipolePointCoords[kk];
254 branchMultipole.degree = MAX_MULTIPOLE_DEGREE;
255 branchMultipole.noOfMoments = MAX_NO_OF_MOMENTS_PER_MULTIPOLE;
256
257 const std::vector<batch_struct> & batchList = org.batchList;
258 const std::vector<cluster_struct> & clusterList = org.clusterList;
259 const std::vector<multipole_struct_small> & multipoleListForGroups = result_org_mm.multipoleListForGroups;
260
261 int batchCount = org.batchList.size();
262 for(int batchIndex = 0; batchIndex < batchCount; batchIndex++) {
263 int clusterCount = batchList[batchIndex].noOfClusters;
264 int cluster_start = batchList[batchIndex].clusterStartIndex;
265 for(int clusterIndex = cluster_start; clusterIndex < cluster_start + clusterCount; clusterIndex++) {
266 int group_start = clusterList[clusterIndex].groupStartIndex;
267 int group_end = group_start + clusterList[clusterIndex].noOfGroups;
268 for(int groupIndex = group_start; groupIndex < group_end; groupIndex++) {
269 const multipole_struct_small & multipoleCurrGroup = multipoleListForGroups[groupIndex];
270
271 // take multipole for this group, and translate it to center-of-charge point
272 ergo_real dx = multipoleCurrGroup.centerCoords[0] - multipolePointCoords[0];
273 ergo_real dy = multipoleCurrGroup.centerCoords[1] - multipolePointCoords[1];
274 ergo_real dz = multipoleCurrGroup.centerCoords[2] - multipolePointCoords[2];
275
276 ergo_real W[MAX_NO_OF_MOMENTS_PER_MULTIPOLE*MAX_NO_OF_MOMENTS_PER_MULTIPOLE];
277 translator.getTranslationMatrix
278 (dx, dy, dz, MAX_MULTIPOLE_DEGREE,
279 multipoleCurrGroup.degree, W);
280
281 multipole_struct_large translatedMultipole;
282 for(int A = 0; A < MAX_NO_OF_MOMENTS_PER_MULTIPOLE; A++) {
283 ergo_real sum = 0;
284 for(int B = 0; B < multipoleCurrGroup.noOfMoments; B++)
285 sum += W[A*multipoleCurrGroup.noOfMoments+B] * multipoleCurrGroup.momentList[B];
286 translatedMultipole.momentList[A] = sum;
287 } // END FOR A
288 for(int kk = 0; kk < 3; kk++)
289 translatedMultipole.centerCoords[kk] = multipolePointCoords[kk];
290 translatedMultipole.degree = MAX_MULTIPOLE_DEGREE;
291 translatedMultipole.noOfMoments = MAX_NO_OF_MOMENTS_PER_MULTIPOLE;
292
293 // add translated multipole to branch multipole
294 for(int A = 0; A < MAX_NO_OF_MOMENTS_PER_MULTIPOLE; A++)
295 branchMultipole.momentList[A] += translatedMultipole.momentList[A];
296 } // END FOR groupIndex
297 } // END FOR clusterIndex
298 } // END FOR batchIndex
299 setup_multipole_maxAbsMomentList(&branchMultipole);
300 result_org_mm.data.multipole = branchMultipole;
301
302 return 0;
303 }
304
305
306
307 int
combine_mm_info_for_child_boxes(distr_list_description_struct & result_box_branch,const distr_list_description_struct ** child_box_branches,int noOfChildren,const MMTranslator & translator)308 combine_mm_info_for_child_boxes(distr_list_description_struct & result_box_branch,
309 const distr_list_description_struct** child_box_branches,
310 int noOfChildren,
311 const MMTranslator & translator) {
312 multipole_struct_large & newMultipole = result_box_branch.org_mm.data.multipole;
313 for(int A = 0; A < MAX_NO_OF_MOMENTS_PER_MULTIPOLE; A++)
314 newMultipole.momentList[A] = 0;
315
316 // get average position of child multipoles
317 ergo_real avgPosList[3];
318 for(int kk = 0; kk < 3; kk++)
319 avgPosList[kk] = 0;
320
321 ergo_real* maxMomentVectorNormForDistrsList = result_box_branch.org_mm.data.maxMomentVectorNormForDistrsList;
322 for(int l = 0; l <= MAX_MULTIPOLE_DEGREE_BASIC; l++)
323 maxMomentVectorNormForDistrsList[l] = 0;
324
325 for(int childIndex = 0; childIndex < noOfChildren; childIndex++) {
326 for(int kk = 0; kk < 3; kk++)
327 avgPosList[kk] += child_box_branches[childIndex]->org_mm.data.multipole.centerCoords[kk];
328 } // END FOR childIndex
329
330 for(int kk = 0; kk < 3; kk++)
331 newMultipole.centerCoords[kk] = avgPosList[kk] / noOfChildren;
332 newMultipole.degree = MAX_MULTIPOLE_DEGREE;
333 newMultipole.noOfMoments = MAX_NO_OF_MOMENTS_PER_MULTIPOLE;
334
335 // We also want to get maxExtent and maxDistanceOutsideBox for parent box (use largest values found among the children).
336 ergo_real maxExtent = 0;
337 ergo_real maxDistanceOutsideBox = 0;
338
339 // Now translate child multipoles and add to parent multipole
340 for(int childIndex = 0; childIndex < noOfChildren; childIndex++) {
341 const multipole_struct_large* childMultipole = &child_box_branches[childIndex]->org_mm.data.multipole;
342
343 if(child_box_branches[childIndex]->org.data.maxExtent > maxExtent)
344 maxExtent = child_box_branches[childIndex]->org.data.maxExtent;
345
346 if(child_box_branches[childIndex]->org.data.maxDistanceOutsideBox > maxDistanceOutsideBox)
347 maxDistanceOutsideBox = child_box_branches[childIndex]->org.data.maxDistanceOutsideBox;
348
349 ergo_real dx = childMultipole->centerCoords[0] - newMultipole.centerCoords[0];
350 ergo_real dy = childMultipole->centerCoords[1] - newMultipole.centerCoords[1];
351 ergo_real dz = childMultipole->centerCoords[2] - newMultipole.centerCoords[2];
352
353 ergo_real W[MAX_NO_OF_MOMENTS_PER_MULTIPOLE*MAX_NO_OF_MOMENTS_PER_MULTIPOLE];
354 translator.getTranslationMatrix(dx, dy, dz,
355 MAX_MULTIPOLE_DEGREE,
356 MAX_MULTIPOLE_DEGREE, W);
357
358 multipole_struct_large translatedMultipole;
359 for(int A = 0; A < MAX_NO_OF_MOMENTS_PER_MULTIPOLE; A++) {
360 ergo_real sum = 0;
361 for(int B = 0; B < MAX_NO_OF_MOMENTS_PER_MULTIPOLE; B++)
362 sum += W[A*MAX_NO_OF_MOMENTS_PER_MULTIPOLE+B] * childMultipole->momentList[B];
363 translatedMultipole.momentList[A] = sum;
364 } // END FOR A
365 for(int kk = 0; kk < 3; kk++)
366 translatedMultipole.centerCoords[kk] = newMultipole.centerCoords[kk];
367 translatedMultipole.degree = MAX_MULTIPOLE_DEGREE;
368 translatedMultipole.noOfMoments = MAX_NO_OF_MOMENTS_PER_MULTIPOLE;
369
370 // add translated multipole to parent multipole
371 for(int A = 0; A < MAX_NO_OF_MOMENTS_PER_MULTIPOLE; A++)
372 newMultipole.momentList[A] += translatedMultipole.momentList[A];
373
374 for(int l = 0; l <= MAX_MULTIPOLE_DEGREE_BASIC; l++) {
375 ergo_real childValue = child_box_branches[childIndex]->org_mm.data.maxMomentVectorNormForDistrsList[l];
376 if(childValue > maxMomentVectorNormForDistrsList[l])
377 maxMomentVectorNormForDistrsList[l] = childValue;
378 }
379
380 } // END FOR childIndex
381
382 setup_multipole_maxAbsMomentList(&newMultipole);
383
384 result_box_branch.org.data.maxExtent = maxExtent;
385 result_box_branch.org.data.maxDistanceOutsideBox = maxDistanceOutsideBox;
386
387 return 0;
388 }
389