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 xcmat_test.cc Tests the DFT XC matrix construction.
31 This test computes the XC energy many times and checks that the
32 resulting energy is the same every time. If this fails, it is
33 probably because of some bug related to synchronization
34 of threads.
35 */
36
37 #include <stdio.h>
38 #include <unistd.h>
39 #include <memory>
40 #include <limits>
41
42 #include "integrals_1el_potential.h"
43 #include "integrals_2el.h"
44 #include "memorymanag.h"
45 #include "grid_reader.h"
46 #include "dft_common.h"
47 #include "xc_matrix.h"
48
49 static bool
compare_matrices(char mat_name,const real * computed,const long double * ref,int sz,ergo_real eps)50 compare_matrices(char mat_name,
51 const real *computed, const long double *ref, int sz,
52 ergo_real eps)
53 {
54 bool failed = false;
55 for(int row=0; row<sz; row++) {
56 for(int col=0; col<sz; col++) {
57 real currdiff = computed[row + col*sz]- ref[row+col*sz];
58 if (template_blas_fabs(currdiff)>eps) {
59 printf("%c (%d,%d): ref: %28.25Lf got: %28.25Lf diff: %12g eps: %g\n",
60 mat_name, row, col,
61 (long double)ref[row + col*sz],
62 (long double)computed[row + col*sz],
63 (double)currdiff,
64 (double)eps);
65 failed = true;
66 }
67 }
68 }
69 return failed;
70 }
71
72 static int
test_small(const IntegralInfo & ii,const char * functional,const Dft::GridParams::RadialScheme & gridScheme,const char * gridSchemeName,const int * charges,const real (* coords)[3],const long double (* XCRef)[2])73 test_small(const IntegralInfo& ii, const char *functional,
74 const Dft::GridParams::RadialScheme& gridScheme,
75 const char *gridSchemeName,
76 const int *charges, const real (*coords)[3],
77 const long double (*XCRef)[2])
78 {
79 BasisInfoStruct* bis = new BasisInfoStruct();
80 Molecule m;
81 /* The code later will change the order of atoms, this is why the
82 reference table may seem strange at the first sight. */
83 for(int i=0; i<2; i++) {
84 m.addAtom(charges[i], coords[i][0],coords[i][1],coords[i][2]);
85 }
86
87 if(bis->addBasisfuncsForMolecule(m, ERGO_SPREFIX "/basis/STO-3G",
88 0, NULL, ii, 0, 0, 0) != 0) {
89 printf("bis->addBasisfuncsForMolecule failed.\n");
90 return 1;
91 }
92
93 int n = bis->noOfBasisFuncs;
94
95 /* set up density matrix */
96 ergo_real *dmat= ergo_new(n*n, ergo_real);
97 dmat[0*n+0] = 1.1; dmat[0*n+1] = 0.2;
98 dmat[1*n+0] = 0.2; dmat[1*n+1] = 1.3;
99
100 dft_init();
101 if(dft_setfunc(functional) == 0)
102 {
103 printf("error in dft_setfunc\n");
104 return 1;
105 }
106 grid_set_tmpdir("/tmp");
107 static const ergo_real GRID_CELL_SIZE = 2.5;
108 Dft::GridParams gridParams(1e-5, 6, 7, GRID_CELL_SIZE);
109 gridParams.radialGridScheme = gridScheme;
110
111 ergo_real *xcmat= ergo_new(n*n, ergo_real);
112 ergo_real *xca = ergo_new(n*n, ergo_real);
113 ergo_real *xcb = ergo_new(n*n, ergo_real);
114 ergo_real *dmata = ergo_new(n*n, ergo_real);
115 for(int i=n*n-1; i>=0; --i) dmata[i] = 0.5*dmat[i];
116
117 int noOfElectrons = 2;
118 char mode;
119 ergo_real dftEnergy = 0;
120 dft_get_xc_mt(noOfElectrons, dmat, bis, &m, gridParams, xcmat, &dftEnergy);
121 /* We give some room to accumulation error. */
122 /* Since the reference values are computed using long double we
123 cannot check more accurately than that. */
124 ergo_real EPS = mat::getMachineEpsilon<ergo_real>();
125 if(EPS < mat::getMachineEpsilon<long double>())
126 EPS = mat::getMachineEpsilon<long double>();
127 ergo_real extraFactor = 20;
128 if(EPS <= mat::getMachineEpsilon<long double>())
129 extraFactor = 230;
130 /* ELIAS NOTE 2016-09-06: earlier the factor used here was selected like this:
131 (sizeof(ergo_real) == sizeof(ergo_long_real) ? 230 : 20)
132 unclear why that would be needed.
133 FIXME: find out why long double does not give as high accuracy as
134 expected here. */
135 EPS *= extraFactor;
136 int nrepeat = 2;
137 bool failed = false;
138 for(int i = 0; i < nrepeat; i++)
139 {
140 mode = 'R';
141 ergo_real dftEnergyAgain = 0, electronsR, electronsU, dftEnergyU;
142 memset(xcmat, 0, n*n*sizeof(ergo_real));
143 electronsR = dft_get_xc_mt(noOfElectrons, dmat, bis, &m, gridParams,
144 xcmat, &dftEnergyAgain);
145 failed = compare_matrices('R', xcmat, &XCRef[0][0], n, EPS);
146 if(template_blas_fabs(dftEnergyAgain - dftEnergy) > EPS)
147 {
148 printf("%s/%s energy repeatability test failed.\n",
149 selected_func->is_gga() ? "GGA" : "LDA", functional);
150 printf("i = %5i of %5i: computed: %20.19f diff: %g\n",
151 i, nrepeat,
152 (double)dftEnergyAgain, (double)(dftEnergy-dftEnergyAgain));
153 failed = true;
154 }
155 if(failed)
156 break;
157
158 mode = 'U';
159 memset(xca, 0, n*n*sizeof(ergo_real));
160 memset(xcb, 0, n*n*sizeof(ergo_real));
161 electronsU = dft_get_uxc_mt(noOfElectrons,
162 dmata, dmata,
163 bis, &m, gridParams,
164 xca, xcb, &dftEnergyU);
165 failed = compare_matrices('A', xca, &XCRef[0][0], n, EPS)
166 || compare_matrices('B', xcb, &XCRef[0][0], n, EPS);
167 if (template_blas_fabs(electronsU - electronsR) > EPS) {
168 printf("%s/%s Electrons restricted %28.25Lg unrestricted %28.25Lg\n",
169 selected_func->is_gga() ? "GGA" : "LDA", functional,
170 (long double)electronsR,
171 (long double)electronsU);
172 }
173 if(failed)
174 break;
175 }
176
177 ergo_free(dmat);
178 ergo_free(dmata);
179 ergo_free(xcmat);
180 ergo_free(xca);
181 ergo_free(xcb);
182 grid_free_files();
183 delete bis;
184 printf("%cXC %s %s/%s test %s\n", failed ? mode : ' ',
185 gridSchemeName,
186 selected_func->is_gga() ? "GGA" : "LDA",
187 functional, failed ? "failed" : "OK");
188 if(!failed)
189 unlink("ergoscf.out");
190 return failed ? 1 : 0;
191 }
192
193 static int
test_small_both()194 test_small_both()
195 {
196 int res = 0;
197 IntegralInfo ii(true);
198
199 static const int sys1Z[2] = { 2, 1 };
200 static const ergo_real sys1C[2][3] = { { 0, 0, 0 }, { 0, 0, 1.5 } };
201 static const int sys2Z[2] = { 2, 1 };
202 static const ergo_real sys2C[2][3] = { { 0, 0, 0 }, { 0, 0, 20.0 } };
203
204 /* ELIAS NOTE 2016-09-06: new reference values were computed (using
205 gcc 6.1.1) because the test failed after the previous fix adding
206 the L suffix for the pi constant in pi.h which is apparently used
207 by the XC computations. So the old reference values were wrong
208 because they were computed using an inaccurate pi value. Of
209 course it's a bad idea to have this kind of test, it just
210 verifies that we get the same result as earlier, it does not
211 really check if it is correct. But now it is like this, and it is
212 still better than no test at all. */
213 static const long double XCRefSys1Svwn5_GC2[2][2] = {
214 { -0.46846890237096027839L, -0.27741047198177963827L },
215 { -0.27741047198177963827L, -0.63052487439961731318L }
216 };
217 static const long double XCRefSys1Svwn5_Turbo[2][2] = {
218 { -0.46837701264321619278L, -0.27730481689558706597L },
219 { -0.27730481689558706597L, -0.63065863444345805105L }
220 };
221 static const long double XCRefSys1Svwn5_LMG[2][2] = {
222 { -0.46847089939399400585L, -0.27732916610776639027L },
223 { -0.27732916610776639027L, -0.63074970313439271847L }
224 };
225 static const long double XCRefSys1BP86_LMG[2][2] = {
226 { -0.48456321202299733965L, -0.28477354319527885414L },
227 { -0.28477354319527885414L, -0.65856888979121373437L }
228 };
229 static const long double XCRefSys1BP86_TURBO[2][2] = {
230 { -0.48447235314731952526L, -0.28476089225530221456L },
231 { -0.28476089225530221456L, -0.65847904553389222668L }
232 };
233 static const long double XCRefSys2Combine[2][2] = {
234 { -0.41581586869051081036L, 0.00000000000000000000L },
235 { 0.00000000000000000000L, -0.58371746046633447966L }
236 };
237 static const long double XCRefSys2Blyp[2][2] = {
238 { -0.43551596812630800643L, 0.00000000000000000000L },
239 { 0.00000000000000000000L, -0.61589810422545549677L }
240 };
241
242 res += test_small(ii, "SVWN5", Dft::GridParams::GC2, "GC2 ",
243 sys1Z, sys1C, &XCRefSys1Svwn5_GC2[0]);
244
245 res += test_small(ii, "SVWN5", Dft::GridParams::TURBO, "Turbo",
246 sys1Z, sys1C, &XCRefSys1Svwn5_Turbo[0]);
247
248 res += test_small(ii, "SVWN5", Dft::GridParams::LMG, "LMG ",
249 sys1Z, sys1C, &XCRefSys1Svwn5_LMG[0]);
250
251 res += test_small(ii, "BP86", Dft::GridParams::LMG, "LMG ",
252 sys1Z, sys1C, &XCRefSys1BP86_LMG[0]);
253
254 res += test_small(ii, "BP86", Dft::GridParams::TURBO, "Turbo",
255 sys1Z, sys1C, &XCRefSys1BP86_TURBO[0]);
256
257 res += test_small(ii, "Combine Slater=1 PZ81=1",
258 Dft::GridParams::LMG, "LMG ",
259 sys2Z, sys2C, &XCRefSys2Combine[0]);
260
261 res += test_small(ii, "BLYP", Dft::GridParams::LMG, "LMG ",
262 sys2Z, sys2C, &XCRefSys2Blyp[0]);
263 return res;
264 }
265
266 static int
test_mol(const char * mol_fname,const char * basisSet,const char * xcFunc)267 test_mol(const char *mol_fname, const char *basisSet, const char *xcFunc)
268 {
269 unlink("ergoscf.out");
270 dft_init();
271
272 Molecule m;
273
274 if(m.setFromMoleculeFile(mol_fname, 0, NULL) != 0) {
275 printf("Molecule::setFromMoleculeFile failed.\n");
276 return 1;
277 }
278
279 IntegralInfo integralInfo(true);
280 BasisInfoStruct *bis = new BasisInfoStruct();
281
282 if(bis->addBasisfuncsForMolecule(m, basisSet, 0, NULL,
283 integralInfo, 0, 1, 0) != 0) {
284 printf("bis->addBasisfuncsForMolecule failed.\n");
285 delete bis;
286 return 1;
287 }
288
289 int n = bis->noOfBasisFuncs;
290
291 /* Set up density matrix. Equivalent to use_simple_starting_guess. */
292 ergo_real *dmat= ergo_new(n*n, ergo_real);
293 int noOfElectrons = m.getNumberOfElectrons();
294 real diag = noOfElectrons/ergo_real(n);
295 for(int col=0; col<n; col++) {
296 for(int row=0; row<n; row++)
297 dmat[row+n*col] = 0.0;
298 dmat[col+n*col] = diag;
299 }
300
301 if(dft_setfunc(xcFunc) == 0) {
302 fprintf(stderr, "Error in dft_setfunc(%s)\n", xcFunc);
303 return 1;
304 }
305 grid_set_tmpdir("/tmp");
306 static const int ANGMIN = 6;
307 static const int ANGINT = 35;
308 static const real RADINT = 1e-7;
309
310 Dft::GridParams gridParams(RADINT, ANGMIN, ANGINT);
311
312 ergo_real *xcmat= ergo_new(n*n, ergo_real);
313
314 ergo_real dftEnergy = 0;
315 ergo_real integratedNoOfElectrons =
316 dft_get_xc_mt(noOfElectrons, dmat, bis, &m, gridParams, xcmat, &dftEnergy);
317
318 ergo_free(dmat);
319 ergo_free(xcmat);
320 grid_free_files();
321 delete bis;
322 printf("%s/%s benchmark executed. "
323 "Expected %d electrons. Integrated %f\n",
324 selected_func->is_gga() ? "GGA" : "LDA",
325 xcFunc, noOfElectrons, (double)integratedNoOfElectrons);
326 return 0;
327 }
328
main(int argc,char * argv[])329 int main(int argc, char *argv[])
330 {
331 const char *DEFAULT_XC_FUNC = "BP86";
332 if(argc<=1)
333 return test_small_both();
334 else if(argc==2)
335 return test_mol(argv[1], ERGO_SPREFIX "/basis/Turbomole-SVP",
336 DEFAULT_XC_FUNC);
337 else if(argc==3)
338 return test_mol(argv[1], argv[2], DEFAULT_XC_FUNC);
339 else if(argc==4)
340 return test_mol(argv[1], argv[2], argv[3]);
341 else {
342 fputs("Usage: xcmat_test [MOL_FILE [BASIS_SET]].\n", stderr);
343 return 1;
344 }
345 /* Not reached */
346 return 0;
347 }
348