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 densitymanager.cc
31
32 @brief Functionality for working with the electron density as a
33 function of space, for a given basis set and density matrix.
34
35 @author: Elias Rudberg <em>responsible</em>
36 */
37
38 #include <stdlib.h>
39 #include <memory.h>
40 #include <math.h>
41 #include "memorymanag.h"
42 #include "output.h"
43 #include "densitymanager.h"
44 #include "pi.h"
45 #include "integrals_general.h"
46 #include "template_blas_common.h"
47
48
49 #define EXPONENT_DIFF_LIMIT 1e-22
50 #define DISTR_CENTER_DIST_LIMIT 1e-22
51
52
53 static ergo_real
compute_1d_gaussian_integral_recursive(ergo_real a,ergo_real b,int n,ergo_real alpha)54 compute_1d_gaussian_integral_recursive(ergo_real a, ergo_real b, int n, ergo_real alpha)
55 {
56 ergo_real result, sqrtalpha, term1, term2;
57 ergo_real aToPowerNminus1, bToPowerNminus1;
58 if(n == 0)
59 {
60 sqrtalpha = template_blas_sqrt(alpha);
61 result = template_blas_sqrt(pi/(4*alpha)) * (template_blas_erf(sqrtalpha*b) - template_blas_erf(sqrtalpha*a));
62 return result;
63 }
64 if(n == 1)
65 {
66 result = -(1 / (2*alpha)) * (template_blas_exp(-alpha*b*b) - template_blas_exp(-alpha*a*a));
67 return result;
68 }
69 if(n < 0)
70 {
71 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in 1dintegral: n < 0");
72 exit(0);
73 }
74 /* now we know that n >= 2 */
75 term1 = (n - 1) * compute_1d_gaussian_integral_recursive(a, b, n-2, alpha);
76 aToPowerNminus1 = template_blas_pow(a, (ergo_real)n-1);
77 bToPowerNminus1 = template_blas_pow(b, (ergo_real)n-1);
78 term2 =
79 bToPowerNminus1 * template_blas_exp(-alpha*b*b) -
80 aToPowerNminus1 * template_blas_exp(-alpha*a*a);
81 result = (term1 - term2) / (2 * alpha);
82 /* return 0; */
83 return result;
84 } /* END compute_1d_gaussian_integral_recursive */
85
86
87
88 static ergo_real
compute_integral_over_box(DistributionSpecStruct * distr,ergo_real * minVect,ergo_real * maxVect,std::vector<int> monomialIntsAdd=std::vector<int> (3,0))89 compute_integral_over_box(DistributionSpecStruct* distr,
90 ergo_real* minVect, ergo_real* maxVect,
91 std::vector<int> monomialIntsAdd = std::vector<int>(3, 0))
92 {
93 ergo_real result, a, b, alpha;
94 int i, n;
95
96 result = distr->coeff;
97 alpha = distr->exponent;
98 for(i = 0; i < 3; i++) // run over the coordinates
99 {
100 n = distr->monomialInts[i];
101 a = minVect[i] - distr->centerCoords[i];
102 b = maxVect[i] - distr->centerCoords[i];
103 result *= compute_1d_gaussian_integral_recursive(a, b, n, alpha);
104 } /* END FOR i */
105 return result;
106 } /* END compute_integral_over_box */
107
108
109
110
111
112
113
114 ergo_real
integrate_density_in_box(int nPrims,DistributionSpecStruct * rho,ergo_real mid_x,ergo_real mid_y,ergo_real mid_z,ergo_real box_width)115 integrate_density_in_box(int nPrims,
116 DistributionSpecStruct* rho,
117 ergo_real mid_x,
118 ergo_real mid_y,
119 ergo_real mid_z,
120 ergo_real box_width)
121 {
122 ergo_real minVect[3];
123 ergo_real maxVect[3];
124 minVect[0] = mid_x - 0.5 * box_width;
125 maxVect[0] = mid_x + 0.5 * box_width;
126 minVect[1] = mid_y - 0.5 * box_width;
127 maxVect[1] = mid_y + 0.5 * box_width;
128 minVect[2] = mid_z - 0.5 * box_width;
129 maxVect[2] = mid_z + 0.5 * box_width;
130 ergo_real sum = 0;
131 int i;
132 for(i = 0; i < nPrims; i++)
133 sum += compute_integral_over_box(&rho[i], minVect, maxVect);
134 return sum;
135 }
136
137
138
139 ergo_real
integrate_density_in_box_2(int nPrims,DistributionSpecStruct * rho,ergo_real * minVect,ergo_real * maxVect,std::vector<int> monomialIntsAdd)140 integrate_density_in_box_2(int nPrims,
141 DistributionSpecStruct* rho,
142 ergo_real* minVect,
143 ergo_real* maxVect,
144 std::vector<int> monomialIntsAdd)
145 {
146 ergo_real sum = 0;
147 int i;
148 for(i = 0; i < nPrims; i++)
149 sum += compute_integral_over_box(&rho[i], minVect, maxVect, monomialIntsAdd);
150 return sum;
151 }
152
153
154
155
156 int
get_no_of_primitives_for_density(ergo_real cutoff,const ergo_real * dmat,const BasisInfoStruct & basisInfo)157 get_no_of_primitives_for_density(ergo_real cutoff,
158 const ergo_real *dmat,
159 const BasisInfoStruct & basisInfo)
160 {
161 #define MAX_DISTR_IN_TEMP_LIST 888
162
163 int i, j;
164 int symmetryFactor;
165 int nBasisFuncs, nn;
166
167 nBasisFuncs = basisInfo.noOfBasisFuncs;
168 nn = 0;
169 for(i = 0; i < nBasisFuncs; i++)
170 {
171 for(j = 0; j < nBasisFuncs; j++)
172 {
173 DistributionSpecStruct tempList[MAX_DISTR_IN_TEMP_LIST];
174 int nPrimitives, k;
175 /* the matrix M is symmetric: include diagonal terms once, */
176 /* and include upper off-diagonal terms multiplied by 2 */
177 if(i == j)
178 symmetryFactor = 1;
179 else
180 symmetryFactor = 2;
181 if(i > j)
182 continue;
183 nPrimitives =
184 get_product_simple_primitives(basisInfo, i,
185 basisInfo, j,
186 tempList,
187 MAX_DISTR_IN_TEMP_LIST,
188 0);
189 if(nPrimitives <= 0)
190 {
191 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in get_product_simple_primitives");
192 return -1;
193 }
194 for(k = 0; k < nPrimitives; k++)
195 {
196 DistributionSpecStruct* currDistr = &tempList[k];
197 ergo_real Mij = dmat[i*nBasisFuncs+j];
198 ergo_real newCoeff = currDistr->coeff * Mij * symmetryFactor;
199 if(template_blas_fabs(newCoeff) > cutoff)
200 nn++;
201 }
202 }
203 }
204 return nn;
205 }
206
207
208
209
210
211
212
213
214
215
216 static int
do_merge_sort_distrs(int n,DistributionSpecStruct * list,DistributionSpecStruct * workList)217 do_merge_sort_distrs(int n,
218 DistributionSpecStruct* list,
219 DistributionSpecStruct* workList)
220 {
221 /* merge sort: */
222 /* first sort the first half, */
223 /* then sort the second half, */
224 /* then merge results to form final sorted list. */
225 int n1, n2, nn, decision, i1, i2, i;
226 DistributionSpecStruct* d1;
227 DistributionSpecStruct* d2;
228
229 if(n == 0)
230 return 0;
231
232 if(n < 1)
233 {
234 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "(n < 1)");
235 return -1;
236 }
237 if(n == 1)
238 return 0;
239
240 n1 = n / 2;
241 n2 = n - n1;
242
243 /* sort first half */
244 if(do_merge_sort_distrs(n1, list, workList) != 0)
245 return -1;
246
247 /* sort second half */
248 if(do_merge_sort_distrs(n2, &list[n1], workList) != 0)
249 return -1;
250
251 /* merge results */
252 nn = 0;
253 i1 = 0;
254 i2 = 0;
255 while(nn < n)
256 {
257 if((i1 < n1) && (i2 < n2))
258 {
259 /* compare */
260 d1 = &list[i1];
261 d2 = &list[n1+i2];
262 decision = 0;
263 for(i = 0; i < 3; i++)
264 {
265 if(decision == 0)
266 {
267 if(d1->monomialInts[i] != d2->monomialInts[i])
268 {
269 if(d1->monomialInts[i] > d2->monomialInts[i])
270 decision = 1;
271 else
272 decision = 2;
273 }
274 } /* END IF (decision == 0) */
275 } /* END FOR i */
276 if(decision == 0)
277 {
278 /* check exponents */
279 if(d1->exponent > d2->exponent)
280 decision = 1;
281 else
282 decision = 2;
283 }
284 }
285 else
286 {
287 if(i1 == n1)
288 decision = 2;
289 else
290 decision = 1;
291 }
292 if(decision <= 0)
293 {
294 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "(decision <= 0)");
295 return -1;
296 }
297 if(decision == 1)
298 {
299 memcpy(&workList[nn], &list[i1], sizeof(DistributionSpecStruct));
300 i1++;
301 }
302 else
303 {
304 memcpy(&workList[nn], &list[n1+i2], sizeof(DistributionSpecStruct));
305 i2++;
306 }
307 nn++;
308 } /* END WHILE (nn < n) */
309 if(i1 != n1)
310 {
311 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "(i1 != n1)");
312 return -1;
313 }
314 if(i2 != n2)
315 {
316 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "(i2 != n2)");
317 return -1;
318 }
319 if(nn != n)
320 {
321 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "(nn != n)");
322 return -1;
323 }
324 memcpy(list, workList, n * sizeof(DistributionSpecStruct));
325 return 0;
326 } /* END do_merge_sort_distrs */
327
328
329
330
331
332
333
334
get_density(const BasisInfoStruct & basisInfo,const ergo_real * dmat,ergo_real cutoff,int maxCountRho,DistributionSpecStruct * resultRho)335 int get_density(const BasisInfoStruct & basisInfo,
336 const ergo_real* dmat,
337 ergo_real cutoff,
338 int maxCountRho,
339 DistributionSpecStruct* resultRho)
340 {
341 #define MAX_DISTR_IN_TEMP_LIST 888
342
343 int i, j, k, kk;
344 DistributionSpecStruct* workList;
345 DistributionSpecStruct* rhoSaved;
346 ergo_real absvalue;
347 ergo_real absdiff;
348 ergo_real sqrtValue;
349 int sameYesNo, firstIndex, count, withinLimit, resultCount;
350 ergo_real coeffSum;
351 int* markList;
352 int symmetryFactor;
353 int nBasisFuncs, nn, nNeededForRho;
354 DistributionSpecStruct* rho;
355
356 nNeededForRho = maxCountRho;
357
358 /* allocate rho */
359 //rho = (DistributionSpecStruct*)ergo_malloc(nNeededForRho * sizeof(DistributionSpecStruct));
360 rho = resultRho;
361
362 nBasisFuncs = basisInfo.noOfBasisFuncs;
363 nn = 0;
364 for(i = 0; i < nBasisFuncs; i++)
365 {
366 for(j = 0; j < nBasisFuncs; j++)
367 {
368 DistributionSpecStruct tempList[MAX_DISTR_IN_TEMP_LIST];
369 int nPrimitives, k;
370 /* the matrix M is symmetric: include diagonal terms once, */
371 /* and include upper off-diagonal terms multiplied by 2 */
372 if(i == j)
373 symmetryFactor = 1;
374 else
375 symmetryFactor = 2;
376 if(i > j)
377 continue;
378 nPrimitives =
379 get_product_simple_primitives(basisInfo, i,
380 basisInfo, j,
381 tempList,
382 MAX_DISTR_IN_TEMP_LIST,
383 0);
384 if(nPrimitives <= 0)
385 {
386 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in get_product_simple_primitives");
387 return -1;
388 }
389 for(k = 0; k < nPrimitives; k++)
390 {
391 DistributionSpecStruct* currDistr = &tempList[k];
392 ergo_real Mij = dmat[i*nBasisFuncs+j];
393 ergo_real newCoeff = currDistr->coeff * Mij * symmetryFactor;
394 if(template_blas_fabs(newCoeff) > cutoff)
395 {
396 /* add to final list */
397 if(nn > nNeededForRho)
398 {
399 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error: (nn > nNeededForRho)");
400 return -1;
401 }
402 memcpy(&rho[nn], currDistr,
403 sizeof(DistributionSpecStruct));
404 rho[nn].coeff = newCoeff;
405 nn++;
406 }
407 }
408 }
409 }
410
411
412
413 /* Now all distributions are stored in the list 'rho'. */
414 /* The number of entries in the list is nn. */
415 /* It could happen that all entries are not unique. */
416 /* We want to join distributions that have the same center */
417 /* and the same exponent. */
418 /* To do this, start with sorting the list by nx, ny, nz, exponent. */
419 workList = (DistributionSpecStruct*)ergo_malloc(nn * sizeof(DistributionSpecStruct));
420 rhoSaved = (DistributionSpecStruct*)ergo_malloc(nn * sizeof(DistributionSpecStruct));
421 memcpy(rhoSaved, rho, nn * sizeof(DistributionSpecStruct));
422
423 if(do_merge_sort_distrs(nn, rho, workList) != 0)
424 {
425 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in do_merge_sort_distrs");
426 return -1;
427 }
428
429
430 /* check that list is sorted */
431 for(i = 0; i < (nn-1); i++)
432 {
433 if(rho[i].exponent < rho[i+1].exponent)
434 {
435 sameYesNo = 1;
436 for(j = 0; j < 3; j++)
437 {
438 if(rho[i].monomialInts[j] != rho[i+1].monomialInts[j])
439 sameYesNo = 0;
440 } /* END FOR j */
441 if(sameYesNo == 1)
442 {
443 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error: distr list NOT properly sorted");
444 return -1;
445 }
446 }
447 } /* END FOR i */
448
449
450 markList = (int*)ergo_malloc(nn * sizeof(int));
451 for(i = 0; i < nn; i++)
452 markList[i] = 0;
453
454 /* now go through sorted list, joining distributions where possible */
455 i = 0;
456 count = 0;
457 firstIndex = 0;
458 while(i < nn)
459 {
460 /* check if this entry has the same nx ny nz as current 'firstIndex' */
461 sameYesNo = 1;
462 for(j = 0; j < 3; j++)
463 {
464 if(rho[i].monomialInts[j] != rho[firstIndex].monomialInts[j])
465 sameYesNo = 0;
466 } /* END FOR j */
467 /* check exponent */
468 absdiff = template_blas_fabs(rho[i].exponent - rho[firstIndex].exponent);
469 if(absdiff > EXPONENT_DIFF_LIMIT)
470 sameYesNo = 0;
471 if(sameYesNo == 0)
472 {
473 for(j = firstIndex; j < i; j++)
474 {
475 if(markList[j] == 0)
476 {
477 markList[j] = 1;
478 /* join distrs that have centers within */
479 /* DISTR_CENTER_DIST_LIMIT of this one */
480 coeffSum = rho[j].coeff;
481 for(k = j+1; k < i; k++)
482 {
483 withinLimit = 1;
484 for(kk = 0; kk < 3; kk++)
485 {
486 absdiff = template_blas_fabs(rho[j].centerCoords[kk] -
487 rho[k].centerCoords[kk]);
488 if(absdiff > DISTR_CENTER_DIST_LIMIT)
489 withinLimit = 0;
490 } /* END FOR kk */
491 if(withinLimit == 1)
492 {
493 coeffSum += rho[k].coeff;
494 markList[k] = 1;
495 }
496 } /* END FOR k */
497 memcpy(&workList[count],
498 &rho[j],
499 sizeof(DistributionSpecStruct));
500 workList[count].coeff = coeffSum;
501 count++;
502 } /* END IF (markList[j] == 0) */
503 } /* END FOR j */
504 firstIndex = i;
505 }
506 else
507 {
508
509 }
510 i++;
511 } /* END WHILE (i < nn) */
512 /* take care of last part */
513 for(j = firstIndex; j < nn; j++)
514 {
515 if(markList[j] == 0)
516 {
517 markList[j] = 1;
518 /* join distrs that have centers within */
519 /* DISTR_CENTER_DIST_LIMIT of this one */
520 coeffSum = rho[j].coeff;
521 for(k = j+1; k < nn; k++)
522 {
523 withinLimit = 1;
524 for(kk = 0; kk < 3; kk++)
525 {
526 absdiff = template_blas_fabs(rho[j].centerCoords[kk] -
527 rho[k].centerCoords[kk]);
528 if(absdiff > DISTR_CENTER_DIST_LIMIT)
529 withinLimit = 0;
530 } /* END FOR kk */
531 if(withinLimit == 1)
532 {
533 coeffSum += rho[k].coeff;
534 markList[k] = 1;
535 }
536 } /* END FOR k */
537 memcpy(&workList[count], &rho[j], sizeof(DistributionSpecStruct));
538 workList[count].coeff = coeffSum;
539 count++;
540 } /* END IF (markList[j] == 0) */
541 } /* END FOR j */
542
543 for(j = 0; j < nn; j++)
544 {
545 if(markList[j] != 1)
546 {
547 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error: (markList[%i] != 1)", j);
548 return -1;
549 }
550 } /* END FOR j */
551
552
553 /* now move results back to list 'rho', */
554 /* skipping those that have too small coeff */
555 resultCount = 0;
556 for(i = 0; i < count; i++)
557 {
558 sqrtValue = template_blas_sqrt(pi / workList[i].exponent);
559 absvalue = workList[i].coeff * sqrtValue * sqrtValue * sqrtValue;
560 if(absvalue < 0) absvalue *= -1;
561 if(absvalue > cutoff)
562 {
563 memcpy(&rho[resultCount],
564 &workList[i],
565 sizeof(DistributionSpecStruct));
566 resultCount++;
567 }
568 } /* END FOR i */
569
570 ergo_free(workList);
571 ergo_free(markList);
572 ergo_free(rhoSaved);
573
574 return resultCount;
575
576 }
577
578
579