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_general.cc
31
32 @brief General functionality related to computation of integrals
33 involving Gaussian basis functions.
34
35 @author: Elias Rudberg <em>responsible</em>
36 */
37
38 /* Written by Elias Rudberg, KTH, Stockholm */
39 #include <stdlib.h>
40 #include <math.h>
41 #include <stdio.h>
42 #include <errno.h>
43 #include <memory.h>
44 #include <time.h>
45 #include <stdarg.h>
46 #include "memorymanag.h"
47 #include "pi.h"
48 #include "output.h"
49 #include "utilities.h"
50 #include "boysfunction.h"
51 #include "integral_info.h"
52 #include "integrals_general.h"
53
54
55
56 #define K_MAX_DIM 44
57
58
59 int
multiply_polynomials(ergo_real result[],polydeg1struct * polydeg1,int dim,ergo_real a[])60 multiply_polynomials(ergo_real result[],
61 polydeg1struct* polydeg1,
62 int dim,
63 ergo_real a[])
64 {
65 int i;
66 ergo_real p1[K_MAX_DIM + 1];
67 ergo_real p2[K_MAX_DIM + 1];
68 if(dim >= (K_MAX_DIM-1))
69 return -1;
70 for(i = 0; i <= dim; i++)
71 p1[i] = a[i]*polydeg1->a0;
72 p1[dim+1] = 0;
73 p2[0] = 0;
74 for(i = 0; i <= dim; i++)
75 p2[i+1] = a[i]*polydeg1->a1;
76 for(i = 0; i <= (dim+1); i++)
77 result[i] = p1[i] + p2[i];
78 return 0;
79 } /* END multiply_polynomials */
80
81
82
83
84
85 /*
86 get_product_simple_prims
87 This function calculates the product of two simple primitives.
88 The result is a list of simple primitives.
89 */
90 int
get_product_simple_prims(const DistributionSpecStruct & primA_in,const DistributionSpecStruct & primB_in,DistributionSpecStruct resultList[],int maxCount,ergo_real threshold)91 get_product_simple_prims(const DistributionSpecStruct& primA_in,
92 const DistributionSpecStruct& primB_in,
93 DistributionSpecStruct resultList[],
94 int maxCount,
95 ergo_real threshold)
96 {
97 // Use a coordinate system with primA at the origin.
98 // This solves the problem with extreme positions of the primitives.
99 DistributionSpecStruct primA_mod = primA_in;
100 DistributionSpecStruct primB_mod = primB_in;
101 int kk;
102 for(kk = 0; kk < 3; kk++)
103 {
104 primA_mod.centerCoords[kk] -= primA_in.centerCoords[kk];
105 primB_mod.centerCoords[kk] -= primA_in.centerCoords[kk];
106 }
107 DistributionSpecStruct* primA = &primA_mod;
108 DistributionSpecStruct* primB = &primB_mod;
109
110 ergo_real CxCyCz, AiAj, alphaNew;
111 ergo_real newCenter[3];
112 ergo_real poly0[K_MAX_DIM];
113 ergo_real poly1[K_MAX_DIM];
114 ergo_real poly2[K_MAX_DIM];
115 ergo_real tempPoly[K_MAX_DIM];
116 ergo_real tempPoly2[K_MAX_DIM];
117 ergo_real tempPoly3[K_MAX_DIM];
118 int tempPolyDegree, tempPoly2Degree;
119 int poly0degree, poly1degree, poly2degree, l, m, nn;
120 polydeg1struct polyDeg1;
121 ergo_real* poly;
122 int* degreePtr;
123 /* use the Gaussian product rule */
124 ergo_real sum = 0;
125 int k;
126 for(k = 0; k < 3; k++)
127 {
128 ergo_real temp = primA->centerCoords[k] - primB->centerCoords[k];
129 sum += temp * temp;
130 } /* END FOR k */
131 CxCyCz = template_blas_exp(-primA->exponent * primB->exponent *
132 sum / (primA->exponent + primB->exponent));
133
134 // FIXME: do this screening properly!
135 if(template_blas_fabs(CxCyCz) < threshold)
136 return 0;
137
138 AiAj = primA->coeff * primB->coeff;
139 alphaNew = primA->exponent + primB->exponent;
140 for(k = 0; k < 3; k++)
141 {
142 newCenter[k] =
143 (primA->exponent * primA->centerCoords[k] +
144 primB->exponent * primB->centerCoords[k]) /
145 (primA->exponent + primB->exponent);
146 } /* END FOR k */
147
148 /* do product of polynomials */
149 /* one coordinate at a time */
150 for(k = 0; k < 3; k++)
151 {
152 switch(k)
153 {
154 case 0: poly = poly0; degreePtr = &poly0degree; break;
155 case 1: poly = poly1; degreePtr = &poly1degree; break;
156 case 2: poly = poly2; degreePtr = &poly2degree; break;
157 default: return -1;
158 } /* END SWITCH k */
159 tempPoly[0] = 1;
160 tempPolyDegree = 0;
161 for(m = 0; m < primA->monomialInts[k]; m++)
162 {
163 polyDeg1.a0 = -primA->centerCoords[k];
164 polyDeg1.a1 = 1;
165 if(multiply_polynomials(tempPoly2, &polyDeg1,
166 tempPolyDegree, tempPoly) != 0)
167 return -1;
168 tempPolyDegree++;
169 memcpy(tempPoly,
170 tempPoly2,
171 (tempPolyDegree+1)*sizeof(ergo_real));
172 } /* END FOR m */
173 for(m = 0; m < primB->monomialInts[k]; m++)
174 {
175 polyDeg1.a0 = -primB->centerCoords[k];
176 polyDeg1.a1 = 1;
177 if(multiply_polynomials(tempPoly2, &polyDeg1,
178 tempPolyDegree, tempPoly) != 0)
179 return -1;
180 tempPolyDegree++;
181 memcpy(tempPoly,
182 tempPoly2,
183 (tempPolyDegree+1)*sizeof(ergo_real));
184 } /* END FOR m */
185
186 /* now do variable change */
187 for(m = 0; m < K_MAX_DIM; m++)
188 poly[m] = 0;
189 tempPoly2Degree = 0;
190 for(m = 0; m <= tempPolyDegree; m++)
191 {
192 tempPoly2[0] = tempPoly[m];
193 tempPoly2Degree = 0;
194 for(l = 0; l < m; l++)
195 {
196 polyDeg1.a0 = newCenter[k];
197 polyDeg1.a1 = 1;
198 if(multiply_polynomials(tempPoly3,
199 &polyDeg1,
200 tempPoly2Degree,
201 tempPoly2) != 0)
202 return -1;
203 tempPoly2Degree++;
204 memcpy(tempPoly2,
205 tempPoly3,
206 (tempPoly2Degree+1)*sizeof(ergo_real));
207 } /* END FOR l */
208 for(l = 0; l <= tempPoly2Degree; l++)
209 {
210 poly[l] += tempPoly2[l];
211 } /* END FOR l */
212 } /* END FOR m */
213 *degreePtr = tempPoly2Degree;
214 } /* END FOR k */
215
216 nn = 0;
217 for(k = 0; k <= poly0degree; k++)
218 {
219 int l;
220 for(l = 0; l <= poly1degree; l++)
221 {
222 int m;
223 for(m = 0; m <= poly2degree; m++)
224 {
225 ergo_real newCoeff = AiAj * CxCyCz * poly0[k] * poly1[l] * poly2[m];
226
227 ergo_real sqrtValue = template_blas_sqrt(pi / alphaNew);
228 ergo_real absvalue = newCoeff * sqrtValue * sqrtValue * sqrtValue;
229 if(absvalue < 0) absvalue *= -1;
230
231 /* add one function to final list */
232 resultList[nn].coeff = newCoeff;
233 resultList[nn].exponent = alphaNew;
234
235 memcpy(resultList[nn].centerCoords,
236 newCenter,
237 3 * sizeof(ergo_real));
238 resultList[nn].monomialInts[0] = k;
239 resultList[nn].monomialInts[1] = l;
240 resultList[nn].monomialInts[2] = m;
241
242 // Translate this term of result back to original coordinate system
243 for(kk = 0; kk < 3; kk++)
244 resultList[nn].centerCoords[kk] += primA_in.centerCoords[kk];
245
246 nn++;
247 if(nn >= maxCount)
248 {
249 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in get_product_simple_prims: "
250 "maxCount exceeded");
251 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "nn = %i, maxCount = %i",
252 nn, maxCount);
253 return -1;
254 }
255 } /* END FOR m */
256 } /* END FOR l */
257 } /* END FOR k */
258
259 return nn;
260 }
261
262
263
264
265
266
267 int
get_product_simple_primitives(const BasisInfoStruct & basisInfoA,int iA,const BasisInfoStruct & basisInfoB,int iB,DistributionSpecStruct resultList[],int maxCount,ergo_real threshold)268 get_product_simple_primitives(const BasisInfoStruct & basisInfoA, int iA,
269 const BasisInfoStruct & basisInfoB, int iB,
270 DistributionSpecStruct resultList[],
271 int maxCount,
272 ergo_real threshold)
273 {
274 BasisFuncStruct* basisFuncA = &basisInfoA.basisFuncList[iA];
275 int nPrimsA = basisFuncA->noOfSimplePrimitives;
276 int Astart = basisFuncA->simplePrimitiveIndex;
277 BasisFuncStruct* basisFuncB = &basisInfoB.basisFuncList[iB];
278 int nPrimsB = basisFuncB->noOfSimplePrimitives;
279 int Bstart = basisFuncB->simplePrimitiveIndex;
280 int n = 0;
281 int i;
282 if((nPrimsA <= 0) || (nPrimsB <= 0))
283 {
284 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in get_product_simple_primitives: "
285 "((nPrimsA <= 0) || (nPrimsB <= 0))\n");
286 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "nPrimsA = %i, nPrimsB = %i\n", nPrimsA, nPrimsB);
287 return -1;
288 }
289 for(i = 0; i < nPrimsA; i++)
290 {
291 const DistributionSpecStruct& primA =
292 basisInfoA.simplePrimitiveList[Astart + i];
293 int j;
294 for(j = 0; j < nPrimsB; j++)
295 {
296 const DistributionSpecStruct& primB =
297 basisInfoB.simplePrimitiveList[Bstart + j];
298
299 int nNewPrims = get_product_simple_prims(primA,
300 primB,
301 &resultList[n],
302 maxCount - n,
303 threshold);
304 if(nNewPrims < 0)
305 {
306 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in get_product_simple_prims");
307 return -1;
308 }
309
310 n += nNewPrims;
311 }
312 }
313 return n;
314 }
315
316
317 ergo_real
compute_integral_of_simple_prim(const DistributionSpecStruct & distr)318 compute_integral_of_simple_prim(const DistributionSpecStruct & distr) {
319 ergo_real result, alpha;
320 int i, n;
321 if( ((distr.monomialInts[0]|
322 distr.monomialInts[1]|
323 distr.monomialInts[2]) & 1) == 1) /* odd integrals disappear */
324 return 0;
325 alpha = distr.exponent;
326 result = distr.coeff * template_blas_pow((ergo_real)pi/alpha, (ergo_real)1.5);
327 ergo_real twoA = 2*alpha;
328 for(i = 0; i < 3; i++) {
329 n = distr.monomialInts[i];
330 for(int j=0; j<n; j+=2)
331 result *= (j+1)/twoA;
332 } /* END FOR i */
333 return result;
334 }
335
336
337 /**
338 Computes the largest integral of any primitive in the basis set,
339 when any x y z factors are ignored. This is useful for getting
340 rough estimates of basis function extents.
341 */
342 ergo_real
get_largest_simple_integral(const BasisInfoStruct & basisInfo)343 get_largest_simple_integral(const BasisInfoStruct & basisInfo)
344 {
345 int n = basisInfo.noOfBasisFuncs;
346 ergo_real A = 0;
347 int i;
348 for(i = 0; i < n; i++)
349 {
350 BasisFuncStruct* basisFunc = &basisInfo.basisFuncList[i];
351 int nPrims = basisFunc->noOfSimplePrimitives;
352 int start = basisFunc->simplePrimitiveIndex;
353 int j;
354 for(j = 0; j < nPrims; j++)
355 {
356 DistributionSpecStruct* prim = &basisInfo.simplePrimitiveList[start + j];
357 DistributionSpecStruct distr;
358 distr = *prim;
359 // Set monomialInts to zero to simplify things
360 distr.monomialInts[0] = 0;
361 distr.monomialInts[1] = 0;
362 distr.monomialInts[2] = 0;
363 ergo_real a = compute_integral_of_simple_prim(distr);
364 if(a > A)
365 A = a;
366 } // END FOR j
367 } // END FOR i
368 return A;
369 }
370
371
372
373 /**
374 Computes an estimate for the largest absolute value that any basis
375 function takes. Useful as "worst case" when you want to find out
376 the largest contribution to the density that a basis function can
377 be part of. */
get_max_basis_func_abs_value(const BasisInfoStruct & basisInfo)378 ergo_real get_max_basis_func_abs_value(const BasisInfoStruct & basisInfo) {
379 int n = basisInfo.noOfBasisFuncs;
380 ergo_real maxValue = 0;
381 for(int i = 0; i < n; i++) {
382 BasisFuncStruct* basisFunc = &basisInfo.basisFuncList[i];
383 int nPrims = basisFunc->noOfSimplePrimitives;
384 int start = basisFunc->simplePrimitiveIndex;
385 for(int j = 0; j < nPrims; j++) {
386 DistributionSpecStruct* prim = &basisInfo.simplePrimitiveList[start + j];
387 ergo_real valueAtCenter = template_blas_fabs(prim->coeff); // exp(0) = 1
388 if(valueAtCenter > maxValue)
389 maxValue = valueAtCenter;
390 } // END FOR j
391 } // END FOR i
392 return maxValue;
393 }
394
395
396 /**
397 Computes an "extent" for each basis function in the basis set.
398 The "extent" is such that the value of the function is smaller
399 than maxAbsValue at distances beyond the "extent".
400 */
401 int
get_basis_func_extent_list(const BasisInfoStruct & basisInfo,ergo_real * basisFuncExtentList,ergo_real maxAbsValue)402 get_basis_func_extent_list(const BasisInfoStruct & basisInfo, ergo_real* basisFuncExtentList, ergo_real maxAbsValue)
403 {
404 int n = basisInfo.noOfBasisFuncs;
405 for(int i = 0; i < n; i++)
406 {
407 BasisFuncStruct* basisFunc = &basisInfo.basisFuncList[i];
408 int nPrims = basisFunc->noOfSimplePrimitives;
409 int start = basisFunc->simplePrimitiveIndex;
410 ergo_real maxExtent = 0;
411 for(int j = 0; j < nPrims; j++)
412 {
413 DistributionSpecStruct* prim = &basisInfo.simplePrimitiveList[start + j];
414 ergo_real currExtent = template_blas_sqrt((1.0 / prim->exponent) * template_blas_log(template_blas_fabs(prim->coeff) / maxAbsValue));
415 if(currExtent > maxExtent)
416 maxExtent = currExtent;
417 } // END FOR j
418 basisFuncExtentList[i] = maxExtent;
419 } // END FOR i
420 return 0;
421 }
422