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 boys_function_test.cc Tests the Boys function evaluation
31 accuracy. */
32
33 #include <stdio.h>
34 #include <unistd.h>
35 #include <memory>
36 #include <limits>
37 #include <vector>
38
39 #include "integral_info.h"
40 #include "matInclude.h"
41 #include "template_blas_common.h"
42 #include "config.h" // Needed to get the PRECISION_SINGLE macro
43
rand_0_to_1()44 static ergo_real rand_0_to_1() {
45 int randomint = rand();
46 ergo_real x = (ergo_real)randomint;
47 return x / RAND_MAX;
48 }
49
rand_1_to_10()50 static ergo_real rand_1_to_10() {
51 ergo_real rand_0_to_9 = rand_0_to_1() * 9;
52 ergo_real result = rand_0_to_9 + 1;
53 return result;
54 }
55
BoysFuncAccurate(int n,ergo_real x,const IntegralInfo & integralInfo)56 static ergo_real BoysFuncAccurate(int n, ergo_real x, const IntegralInfo & integralInfo) {
57 int noOfIntegrationIntervals = 40;
58 ergo_real machine_epsilon = mat::getMachineEpsilon<ergo_real>();
59 ergo_real diffLimit = template_blas_sqrt(template_blas_get_num_limit_min<ergo_real>());
60 ergo_real reldiffLimit = template_blas_pow(machine_epsilon, (ergo_real)0.75);
61 const int NMAX = 30;
62 ergo_real list[NMAX];
63 for(int i = 0; i < NMAX; i++) {
64 list[i] = integralInfo.BoysFunction_expensive(n, x, noOfIntegrationIntervals);
65 assert(list[i] >= 0);
66 if(i >= 2) {
67 // Check if difference between different results is small
68 ergo_real diff1 = template_blas_fabs(list[i-0]-list[i-1]);
69 ergo_real diff2 = template_blas_fabs(list[i-1]-list[i-2]);
70 if(diff1 < diffLimit && diff2 < diffLimit)
71 return list[i];
72 ergo_real reldiff1 = diff1 / list[i];
73 ergo_real reldiff2 = diff2 / list[i];
74 if(reldiff1 < reldiffLimit && reldiff2 < reldiffLimit)
75 return list[i];
76 }
77 noOfIntegrationIntervals *= 2;
78 }
79 printf("ERROR in BoysFuncAccurate(): failed to get accurate result.\n");
80 return 0;
81 }
82
main(int argc,char * argv[])83 int main(int argc, char *argv[]) {
84 IntegralInfo integralInfo(true);
85 integralInfo.init();
86 ergo_real machine_epsilon = mat::getMachineEpsilon<ergo_real>();
87 ergo_real num_limit_min = template_blas_get_num_limit_min<ergo_real>();
88 ergo_real min_value_for_meaningful_reldiff = template_blas_sqrt(num_limit_min);
89 printf("boys_function_test start, machine_epsilon = %g\n", (double)machine_epsilon);
90 // Different things to test: try different values of n and different values of x.
91 // Define constant POW_OF_10_LIMIT that determines the range of sizes of x-values that are tested.
92 #ifdef PRECISION_SINGLE
93 const int POW_OF_10_LIMIT = 4;
94 #else
95 const int POW_OF_10_LIMIT = 10;
96 #endif
97 int saved_n = -1;
98 ergo_real saved_x = 0;
99 ergo_real saved_absdiff = 0;
100 ergo_real saved_value = 0;
101 ergo_real saved_refvalue = 0;
102 ergo_real maxreldiff = 0;
103 // Try all values of n from 0 up to BOYS_N_MAX-1
104 printf("Testing Boys function accuracy for 0 <= n <= %d\n", BOYS_N_MAX-1);
105 for(int n = 0; n < BOYS_N_MAX; n++) {
106 // Try different order of magnitude of x-values
107 for(int i = -POW_OF_10_LIMIT; i < POW_OF_10_LIMIT; i++) {
108 ergo_real x_base = template_blas_pow((ergo_real)10, (ergo_real)i);
109 const int nTests = 10;
110 for(int k = 0; k < nTests; k++) {
111 // Get random x-value around x_base
112 ergo_real x = x_base * rand_1_to_10();
113 // Compare computed Boys func values for that x-value
114 ergo_real boysFuncResultRef = BoysFuncAccurate(n, x, integralInfo);
115 ergo_real boysFuncResult = integralInfo.BoysFunction(n, x);
116 ergo_real absdiff = template_blas_fabs(boysFuncResult - boysFuncResultRef);
117 ergo_real reldiff = 0;
118 if(template_blas_fabs(boysFuncResultRef) > min_value_for_meaningful_reldiff)
119 reldiff = absdiff / template_blas_fabs(boysFuncResultRef);
120 else {
121 // boysFuncResultRef is extremely small. In this case we just check that absdiff is extremely small also
122 if(absdiff > 2*min_value_for_meaningful_reldiff)
123 reldiff = 1;
124 }
125 if(reldiff > maxreldiff) {
126 maxreldiff = reldiff;
127 saved_n = n;
128 saved_x = x;
129 saved_absdiff = absdiff;
130 saved_refvalue = boysFuncResultRef;
131 saved_value = boysFuncResult;
132 }
133 }
134 }
135 }
136 printf("maxreldiff = %g\n", (double)maxreldiff);
137 printf("saved_n = %d\n", saved_n);
138 printf("saved_x = %g\n", (double)saved_x);
139 printf("saved_absdiff = %g\n", (double)saved_absdiff);
140 printf("saved_refvalue = %g\n", (double)saved_refvalue);
141 printf("saved_value = %g\n", (double)saved_value);
142 ergo_real tolerance = template_blas_pow(machine_epsilon, (ergo_real)0.75);
143 if(maxreldiff > tolerance) {
144 printf("Error in boys_function_test: (maxreldiff > tolerance). tolerance= %g\n", (double)tolerance);
145 return -1;
146 }
147 printf("boys_function_test finished OK (maxreldiff was below tolerance %g).\n", (double)tolerance);
148 return 0;
149 }
150