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 integral_matrix_wrappers.cc
31
32 @brief Wrapper routines for different parts of the integral code,
33 including conversion of matrices from/to the hierarchical matrix
34 library (HML) format.
35
36 @author: Elias Rudberg <em>responsible</em>
37 */
38
39 #include <string.h>
40
41 #include "integral_matrix_wrappers.h"
42 #include "output.h"
43 #include "basis_func_pair_list_1el.h"
44 #include "integrals_1el_kinetic.h"
45 #include "integrals_1el_potential.h"
46 #include "operator_matrix.h"
47 #include "basis_func_pair_list.h"
48 #include "integrals_2el_boxed.h"
49 #include "integrals_2el_J.h"
50 #include "integrals_2el_K.h"
51 #include "integrals_2el_K_prep.h"
52 #include "utilities.h"
53 #include "memorymanag.h"
54
55
get_max_charge(const Molecule & molecule)56 static ergo_real get_max_charge(const Molecule& molecule) {
57 ergo_real maxCharge = 0;
58 for(int i = 0; i < molecule.getNoOfAtoms(); i++) {
59 ergo_real currCharge = molecule.getAtom(i).charge;
60 if(currCharge > maxCharge)
61 maxCharge = currCharge;
62 }
63 return maxCharge;
64 }
65
66
convert_symm_CSR_to_HML_and_destroy_CSR(csr_matrix_struct & M_CSR,symmMatrix & M_HML,std::vector<int> const & permutationHML)67 static void convert_symm_CSR_to_HML_and_destroy_CSR(csr_matrix_struct & M_CSR,
68 symmMatrix & M_HML,
69 std::vector<int> const & permutationHML) {
70 // Convert matrix from CSR to HML format.
71 int nvalues = ergo_CSR_get_nvalues(&M_CSR);
72 std::vector<int> rowind(nvalues);
73 std::vector<int> colind(nvalues);
74 std::vector<ergo_real> values(nvalues);
75 if(ergo_CSR_get_values(&M_CSR, rowind, colind, values, nvalues) != 0) {
76 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in ergo_CSR_get_values.");
77 throw std::runtime_error("Error in convert_symm_CSR_to_HML_and_destroy_CSR(): ergo_CSR_get_values() failed.");
78 }
79 // now the information is in the vectors rowind colind values, we can free memory for M_CSR.
80 ergo_CSR_destroy(&M_CSR);
81 M_HML.assign_from_sparse(rowind, colind, values, permutationHML, permutationHML);
82 }
83
84
85 int
compute_V_sparse(const BasisInfoStruct & basisInfo,const IntegralInfo & integralInfo,const Molecule & molecule,ergo_real threshold,ergo_real boxSize,symmMatrix & V,std::vector<int> const & permutationHML,ergo_real & result_nuclearRepulsionEnergy)86 compute_V_sparse(const BasisInfoStruct& basisInfo,
87 const IntegralInfo& integralInfo,
88 const Molecule& molecule,
89 ergo_real threshold,
90 ergo_real boxSize,
91 symmMatrix & V,
92 std::vector<int> const & permutationHML,
93 ergo_real & result_nuclearRepulsionEnergy)
94 {
95 int n = basisInfo.noOfBasisFuncs;
96
97 ergo_real maxCharge = get_max_charge(molecule);
98
99 do_output(LOG_CAT_INFO, LOG_AREA_UNDEFINED,
100 "compute_V_sparse, no of atoms = %5i, maxCharge = %6.2f",
101 molecule.getNoOfAtoms(), (double)maxCharge);
102
103 int noOfBasisFuncIndexPairs = get_basis_func_pair_list_1el_for_V(basisInfo,
104 threshold,
105 boxSize,
106 maxCharge,
107 NULL,
108 2000000000);
109 if(noOfBasisFuncIndexPairs <= 0)
110 {
111 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in get_basis_func_pair_list_1el, noOfBasisFuncIndexPairs = %i", noOfBasisFuncIndexPairs);
112 return -1;
113 }
114 std::vector<basis_func_index_pair_struct_1el> basisFuncIndexPairList(noOfBasisFuncIndexPairs);
115 std::vector<ergo_real> V_list(noOfBasisFuncIndexPairs);
116 noOfBasisFuncIndexPairs = get_basis_func_pair_list_1el_for_V(basisInfo,
117 threshold,
118 boxSize,
119 maxCharge,
120 &basisFuncIndexPairList[0],
121 noOfBasisFuncIndexPairs);
122 if(noOfBasisFuncIndexPairs <= 0)
123 {
124 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in get_basis_func_pair_list_1el, noOfBasisFuncIndexPairs = %i", noOfBasisFuncIndexPairs);
125 return -1;
126 }
127
128 do_output(LOG_CAT_INFO, LOG_AREA_UNDEFINED,
129 "noOfBasisFuncIndexPairs = %i ==> storing %6.2f %% of a full matrix",
130 noOfBasisFuncIndexPairs, (double)100*noOfBasisFuncIndexPairs/((double)n*n));
131
132 if(compute_V_and_gradient_linear(basisInfo,
133 integralInfo,
134 molecule,
135 threshold,
136 boxSize,
137 &basisFuncIndexPairList[0],
138 &V_list[0],
139 noOfBasisFuncIndexPairs,
140 false, NULL, NULL,
141 result_nuclearRepulsionEnergy) != 0)
142 {
143 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in compute_V_and_gradient_linear");
144 return -1;
145 }
146
147 // Now transfer result from V_list to V
148 {
149 std::vector<int> rowind(noOfBasisFuncIndexPairs);
150 std::vector<int> colind(noOfBasisFuncIndexPairs);
151 for(int i = 0; i < noOfBasisFuncIndexPairs; i++)
152 {
153 rowind[i] = basisFuncIndexPairList[i].index_1;
154 colind[i] = basisFuncIndexPairList[i].index_2;
155 }
156 V.assign_from_sparse(rowind,
157 colind,
158 V_list,
159 permutationHML,
160 permutationHML);
161 }
162
163 return 0;
164 }
165
166
167 int
compute_V_sparse_hierarchical(const BasisInfoStruct & basisInfo,const IntegralInfo & integralInfo,const Molecule & molecule,ergo_real threshold,ergo_real boxSize,symmMatrix & V,std::vector<int> const & permutationHML,ergo_real & result_nuclearRepulsionEnergy)168 compute_V_sparse_hierarchical(const BasisInfoStruct& basisInfo,
169 const IntegralInfo& integralInfo,
170 const Molecule& molecule,
171 ergo_real threshold,
172 ergo_real boxSize,
173 symmMatrix & V,
174 std::vector<int> const & permutationHML,
175 ergo_real & result_nuclearRepulsionEnergy) {
176 //
177 // Do "get_basis_func_pair_list_1el_hierarchical" or similar
178 //
179
180 int n = basisInfo.noOfBasisFuncs;
181 ergo_real maxCharge = get_max_charge(molecule);
182 do_output(LOG_CAT_INFO, LOG_AREA_UNDEFINED, "compute_V_sparse_hierarchical, no of atoms = %5i, maxCharge = %6.2f",
183 molecule.getNoOfAtoms(), (double)maxCharge);
184 int noOfBasisFuncIndexPairs = get_basis_func_pair_list_1el_for_V(basisInfo, threshold, boxSize, maxCharge, NULL, 2000000000);
185 if(noOfBasisFuncIndexPairs <= 0) {
186 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in get_basis_func_pair_list_1el, noOfBasisFuncIndexPairs = %i", noOfBasisFuncIndexPairs);
187 return -1;
188 }
189 std::vector<basis_func_index_pair_struct_1el> basisFuncIndexPairList(noOfBasisFuncIndexPairs);
190 noOfBasisFuncIndexPairs = get_basis_func_pair_list_1el_for_V(basisInfo, threshold, boxSize, maxCharge, &basisFuncIndexPairList[0], noOfBasisFuncIndexPairs);
191 if(noOfBasisFuncIndexPairs <= 0) {
192 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in get_basis_func_pair_list_1el, noOfBasisFuncIndexPairs = %i", noOfBasisFuncIndexPairs);
193 return -1;
194 }
195 do_output(LOG_CAT_INFO, LOG_AREA_UNDEFINED, "noOfBasisFuncIndexPairs = %i ==> storing %6.2f %% of a full matrix",
196 noOfBasisFuncIndexPairs, (double)100*noOfBasisFuncIndexPairs/((double)n*n));
197
198 //
199 // Set up CSR matrix for storing resulting V matrix (similar to compute_K_by_boxes_sparse)
200 //
201 csr_matrix_struct V_CSR;
202 {
203 memset(&V_CSR, 0, sizeof(csr_matrix_struct));
204 int symmetryFlagForCSR = 1;
205 std::vector<int> rowind(noOfBasisFuncIndexPairs);
206 std::vector<int> colind(noOfBasisFuncIndexPairs);
207 for(int k = 0; k < noOfBasisFuncIndexPairs; k++) {
208 rowind[k] = basisFuncIndexPairList[k].index_1;
209 colind[k] = basisFuncIndexPairList[k].index_2;
210 }
211 if(ergo_CSR_create(&V_CSR, symmetryFlagForCSR, n, noOfBasisFuncIndexPairs, rowind, colind) != 0) {
212 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in ergo_CSR_create for V_CSR");
213 return -1;
214 }
215 }
216
217 //
218 // Call "compute_V_hierarchical" or similar
219 //
220 if(compute_V_hierarchical(basisInfo, integralInfo, molecule, threshold, boxSize,
221 &basisFuncIndexPairList[0], noOfBasisFuncIndexPairs, &V_CSR,
222 result_nuclearRepulsionEnergy) != 0) {
223 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in compute_V_hierarchical");
224 return -1;
225 }
226
227 //
228 // Convert result V matrix from CSR form to symmMatrix (similar to compute_K_by_boxes_sparse)
229 //
230 convert_symm_CSR_to_HML_and_destroy_CSR(V_CSR, V, permutationHML);
231
232 return 0;
233 }
234
235
236 int
compute_gradient_of_nucl_and_trDV(const BasisInfoStruct & basisInfo,const IntegralInfo & integralInfo,const Molecule & molecule,ergo_real threshold,ergo_real boxSize,const symmMatrix & densityMatrix_sparse,std::vector<int> const & permutationHML,ergo_real * result_gradient_list)237 compute_gradient_of_nucl_and_trDV(const BasisInfoStruct& basisInfo,
238 const IntegralInfo& integralInfo,
239 const Molecule& molecule,
240 ergo_real threshold,
241 ergo_real boxSize,
242 const symmMatrix & densityMatrix_sparse,
243 std::vector<int> const & permutationHML,
244 ergo_real* result_gradient_list) {
245 int n = basisInfo.noOfBasisFuncs;
246
247 ergo_real maxCharge = get_max_charge(molecule);
248
249 do_output(LOG_CAT_INFO, LOG_AREA_UNDEFINED,
250 "compute_gradient_of_tr_D_V, no of atoms = %5i, maxCharge = %6.2f",
251 molecule.getNoOfAtoms(), (double)maxCharge);
252
253 int noOfBasisFuncIndexPairs = get_basis_func_pair_list_1el_for_V(basisInfo, threshold, boxSize, maxCharge, NULL, 2000000000);
254 if(noOfBasisFuncIndexPairs <= 0) {
255 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in get_basis_func_pair_list_1el, noOfBasisFuncIndexPairs = %i", noOfBasisFuncIndexPairs);
256 return -1;
257 }
258 std::vector<basis_func_index_pair_struct_1el> basisFuncIndexPairList(noOfBasisFuncIndexPairs);
259 std::vector<ergo_real> V_list(noOfBasisFuncIndexPairs);
260 noOfBasisFuncIndexPairs = get_basis_func_pair_list_1el_for_V(basisInfo, threshold, boxSize, maxCharge, &basisFuncIndexPairList[0], noOfBasisFuncIndexPairs);
261 if(noOfBasisFuncIndexPairs <= 0) {
262 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in get_basis_func_pair_list_1el, noOfBasisFuncIndexPairs = %i", noOfBasisFuncIndexPairs);
263 return -1;
264 }
265
266 // Get density matrix elements
267 std::vector<ergo_real> D_list(noOfBasisFuncIndexPairs);
268 {
269 std::vector<int> rowind(noOfBasisFuncIndexPairs);
270 std::vector<int> colind(noOfBasisFuncIndexPairs);
271 for(int i = 0; i < noOfBasisFuncIndexPairs; i++) {
272 rowind[i] = basisFuncIndexPairList[i].index_1;
273 colind[i] = basisFuncIndexPairList[i].index_2;
274 }
275 densityMatrix_sparse.get_values(rowind, colind, D_list, permutationHML, permutationHML);
276 }
277
278 do_output(LOG_CAT_INFO, LOG_AREA_UNDEFINED,
279 "noOfBasisFuncIndexPairs = %i ==> storing %6.2f %% of a full matrix",
280 noOfBasisFuncIndexPairs, (double)100*noOfBasisFuncIndexPairs/((double)n*n));
281
282 ergo_real nuclRepEnergyDummy = 0;
283 if(compute_V_and_gradient_linear(basisInfo,
284 integralInfo,
285 molecule,
286 threshold,
287 boxSize,
288 &basisFuncIndexPairList[0],
289 &V_list[0],
290 noOfBasisFuncIndexPairs,
291 true, &D_list[0],
292 result_gradient_list,
293 nuclRepEnergyDummy) != 0)
294 {
295 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in compute_V_and_gradient_linear");
296 return -1;
297 }
298
299 return 0;
300 }
301
302
303 ergo_real
get_electron_nuclear_attraction_energy(const IntegralInfo & integralInfo,const Molecule & molecule,const BasisInfoStruct & basisInfo,const symmMatrix & D,ergo_real threshold_integrals_1el,mat::SizesAndBlocks const & matrix_size_block_info,std::vector<int> const & permutationHML)304 get_electron_nuclear_attraction_energy(const IntegralInfo& integralInfo,
305 const Molecule& molecule,
306 const BasisInfoStruct& basisInfo,
307 const symmMatrix & D,
308 ergo_real threshold_integrals_1el,
309 mat::SizesAndBlocks const & matrix_size_block_info,
310 std::vector<int> const & permutationHML)
311 {
312 // Get V
313 symmMatrix V;
314 V.resetSizesAndBlocks(matrix_size_block_info, matrix_size_block_info);
315 ergo_real boxSize = 10.0;
316 ergo_real nuclRepEnergyDummy = 0;
317 if(compute_V_sparse(basisInfo,
318 integralInfo,
319 molecule,
320 threshold_integrals_1el,
321 boxSize,
322 V,
323 permutationHML,
324 nuclRepEnergyDummy) != 0)
325 throw "error in compute_V_sparse";
326 return symmMatrix::trace_ab(D, V);
327 }
328
329
330 int
compute_T_sparse_linear(const BasisInfoStruct & basisInfo,const IntegralInfo & integralInfo,ergo_real threshold,ergo_real boxSize,symmMatrix & T,std::vector<int> const & permutationHML)331 compute_T_sparse_linear(const BasisInfoStruct& basisInfo,
332 const IntegralInfo& integralInfo,
333 ergo_real threshold,
334 ergo_real boxSize,
335 symmMatrix & T,
336 std::vector<int> const & permutationHML) {
337 int n = basisInfo.noOfBasisFuncs;
338 std::vector<int> nvaluesList(n);
339 std::vector<int*> colindList(n);
340 std::vector<ergo_real*> valuesList(n);
341 if(compute_T_matrix_sparse_linear(basisInfo,
342 threshold,
343 boxSize,
344 &nvaluesList[0],
345 &colindList[0],
346 &valuesList[0]) != 0) {
347 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in compute_T_matrix_sparse");
348 return -1;
349 }
350 // Now convert result to three vectors so prepare for HML storage.
351 int nvalues = 0;
352 for(int i = 0; i < n; i++)
353 nvalues += nvaluesList[i];
354 // allocate vectors
355 std::vector<int> rowind(nvalues);
356 std::vector<int> colind(nvalues);
357 std::vector<ergo_real> values(nvalues);
358 // populate vectors
359 int count = 0;
360 for(int i = 0; i < n; i++)
361 {
362 for(int j = 0; j < nvaluesList[i]; j++)
363 {
364 rowind[count] = i;
365 colind[count] = colindList[i][j];
366 values[count] = valuesList[i][j];
367 count++;
368 } // END FOR j
369 } // END FOR i
370 // Now the information is in rowind colind values.
371 // free memory allocated by compute_operator_matrix_sparse.
372 for(int i = 0; i < n; i++) {
373 ergo_free(colindList[i]);
374 ergo_free(valuesList[i]);
375 }
376 T.assign_from_sparse(rowind,
377 colind,
378 values,
379 permutationHML,
380 permutationHML);
381 return 0;
382 }
383
384
385
386 static int
check_diagonal_elements_of_overlap_matrix(int n,const symmMatrix & S_symm)387 check_diagonal_elements_of_overlap_matrix(int n, const symmMatrix & S_symm)
388 {
389 std::vector<int> rowind(n);
390 std::vector<int> colind(n);
391 std::vector<ergo_real> values(n);
392 for(int i = 0; i < n; i++)
393 {
394 rowind[i] = i;
395 colind[i] = i;
396 }
397 S_symm.get_values(rowind,
398 colind,
399 values);
400 ergo_real maxabsdiff = 0;
401 for(int i = 0; i < n; i++)
402 {
403 ergo_real absdiff = template_blas_fabs(values[i] - 1);
404 if(absdiff > maxabsdiff)
405 maxabsdiff = absdiff;
406 }
407 do_output(LOG_CAT_INFO, LOG_AREA_UNDEFINED, "diag elements of ovl matrix differ from 1 by at most %9.6g",
408 (double)maxabsdiff);
409 if(maxabsdiff > 1e-3)
410 {
411 do_output(LOG_CAT_WARNING, LOG_AREA_UNDEFINED, "WARNING: bad overlap matrix: diag elements of ovl matrix differ from 1 by more than safe limit.");
412 do_output(LOG_CAT_WARNING, LOG_AREA_UNDEFINED, "WARNING: bad overlap matrix: diag elements of ovl matrix differ from 1 by more than safe limit.");
413 do_output(LOG_CAT_WARNING, LOG_AREA_UNDEFINED, "WARNING: bad overlap matrix: diag elements of ovl matrix differ from 1 by more than safe limit.");
414 do_output(LOG_CAT_WARNING, LOG_AREA_UNDEFINED, "WARNING: this means that basis functions are not completely normlized.");
415 do_output(LOG_CAT_WARNING, LOG_AREA_UNDEFINED, "WARNING: this means that basis functions are not completely normlized.");
416 do_output(LOG_CAT_WARNING, LOG_AREA_UNDEFINED, "WARNING: this means that basis functions are not completely normlized.");
417 }
418 return 0;
419 }
420
421
422 int
compute_overlap_matrix_sparse(const BasisInfoStruct & basisInfo,symmMatrix & S_symm,std::vector<int> const & permutationHML)423 compute_overlap_matrix_sparse(const BasisInfoStruct& basisInfo,
424 symmMatrix & S_symm,
425 std::vector<int> const & permutationHML)
426 {
427 if(compute_operator_matrix_sparse_symm(basisInfo,
428 0,
429 0,
430 0,
431 S_symm,
432 permutationHML) != 0)
433 {
434 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in compute_operator_matrix_sparse_symm");
435 return -1;
436 }
437
438 // check diagonal elements of overlap matrix
439 int n = basisInfo.noOfBasisFuncs;
440 check_diagonal_elements_of_overlap_matrix(n, S_symm);
441 return 0;
442 }
443
444
445 int
compute_R_matrix_sparse(const BasisInfoStruct & basisInfo_A,const BasisInfoStruct & basisInfo_B,normalMatrix & result_R,ergo_real sparse_threshold,std::vector<int> const & matrixPermutationVec_A,std::vector<int> const & matrixPermutationVec_B)446 compute_R_matrix_sparse(const BasisInfoStruct & basisInfo_A,
447 const BasisInfoStruct & basisInfo_B,
448 normalMatrix & result_R,
449 ergo_real sparse_threshold,
450 std::vector<int> const & matrixPermutationVec_A,
451 std::vector<int> const & matrixPermutationVec_B)
452 {
453 int n_A = basisInfo_A.noOfBasisFuncs;
454 int n_B = basisInfo_B.noOfBasisFuncs;
455
456 std::vector<int> nvaluesList(n_A);
457 std::vector< std::vector<int> > colindList(n_A);
458 std::vector< std::vector<ergo_real> > valuesList(n_A);
459
460 if(compute_operator_matrix_sparse(basisInfo_A,
461 basisInfo_B,
462 0,
463 0,
464 0,
465 n_A,
466 n_B,
467 nvaluesList,
468 colindList,
469 valuesList) != 0)
470 {
471 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in compute_operator_matrix_sparse");
472 return -1;
473 }
474
475 // Now convert result to three vectors to prepare for HML storage.
476 long nvalues = 0;
477 for(int i = 0; i < n_A; i++)
478 nvalues += nvaluesList[i];
479 // allocate vectors
480 std::vector<int> rowind(nvalues);
481 std::vector<int> colind(nvalues);
482 std::vector<ergo_real> values(nvalues);
483 // populate vectors
484 long count = 0;
485 for(long i = 0; i < n_A; i++)
486 {
487 for(long j = 0; j < nvaluesList[i]; j++)
488 {
489 rowind[count] = i;
490 colind[count] = colindList[i][j];
491 values[count] = valuesList[i][j];
492 count++;
493 } // END FOR j
494 } // END FOR i
495
496 // Now the information is in rowind colind values.
497
498 /* FIXME: there is confusion about the meaning of "rows" and "columns".
499 Here it is reversed, otherwise it did not work.
500 We should use the same convention for rows and columns everywhere. */
501 result_R.assign_from_sparse(colind,
502 rowind,
503 values,
504 matrixPermutationVec_B,
505 matrixPermutationVec_A);
506
507 result_R.eucl_thresh(sparse_threshold);
508
509 // Write to file and read again to reduce memory fragmentation.
510 result_R.writeToFile();
511 result_R.readFromFile();
512
513 return 0;
514 }
515
516
517
518 typedef int* int_ptr;
519 typedef ergo_real* ergo_real_ptr;
520
521 int
compute_operator_matrix_sparse_symm(const BasisInfoStruct & basisInfo,int pow_x,int pow_y,int pow_z,symmMatrix & A_symm,std::vector<int> const & permutationHML)522 compute_operator_matrix_sparse_symm(const BasisInfoStruct& basisInfo,
523 int pow_x,
524 int pow_y,
525 int pow_z,
526 symmMatrix & A_symm,
527 std::vector<int> const & permutationHML)
528 {
529 int n = basisInfo.noOfBasisFuncs;
530
531 std::vector<int> nvaluesList(n);
532 std::vector< std::vector<int> > colindList(n);
533 std::vector< std::vector<ergo_real> > valuesList(n);
534 if(compute_operator_matrix_sparse(basisInfo,
535 basisInfo,
536 pow_x,
537 pow_y,
538 pow_z,
539 n,
540 n,
541 nvaluesList,
542 colindList,
543 valuesList) != 0)
544 {
545 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in compute_operator_matrix_sparse");
546 return -1;
547 }
548
549 // Now convert result to three vectors so prepare for HML storage.
550 int nvalues = 0;
551 for(int i = 0; i < n; i++)
552 nvalues += nvaluesList[i];
553
554 // allocate vectors
555 std::vector<int> rowind(nvalues);
556 std::vector<int> colind(nvalues);
557 std::vector<ergo_real> values(nvalues);
558 // populate vectors
559 int count = 0;
560 for(int i = 0; i < n; i++)
561 {
562 for(int j = 0; j < nvaluesList[i]; j++)
563 {
564 rowind[count] = i;
565 colind[count] = colindList[i][j];
566 values[count] = valuesList[i][j];
567 count++;
568 } // END FOR j
569 } // END FOR i
570
571 normalMatrix A_norm(A_symm);
572
573 A_norm.assign_from_sparse(rowind,
574 colind,
575 values,
576 permutationHML,
577 permutationHML);
578
579 // Convert to symmetric form.
580 A_symm = A_norm;
581 A_norm.clear();
582
583 // Earlier, the matrix was truncated here but we changed this so it is truncated outside this routine.
584
585 // Write to file and read again to reduce memory fragmentation.
586 A_symm.writeToFile();
587 A_symm.readFromFile();
588
589 return 0;
590 }
591
592
593
594
595 int
compute_J_by_boxes_sparse(const BasisInfoStruct & basisInfo,const IntegralInfo & integralInfo,const JK::Params & J_K_params,symmMatrix & J,const symmMatrix & densityMatrix_sparse,std::vector<int> const & permutationHML)596 compute_J_by_boxes_sparse(const BasisInfoStruct& basisInfo,
597 const IntegralInfo& integralInfo,
598 const JK::Params& J_K_params,
599 symmMatrix & J,
600 const symmMatrix & densityMatrix_sparse,
601 std::vector<int> const & permutationHML)
602 {
603 Util::TimeMeter timeMeter;
604
605 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "compute_J_by_boxes_sparse() start.");
606 Util::TimeMeter timeMeterWriteAndReadAll;
607 std::string sizesStr = mat::FileWritable::writeAndReadAll();
608 timeMeterWriteAndReadAll.print(LOG_AREA_DENSFROMF, "FileWritable::writeAndReadAll");
609 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, ((std::string)"writeAndReadAll: '" + sizesStr).c_str());
610
611 int n = basisInfo.noOfBasisFuncs;
612
613 ergo_real maxDensityMatrixElement = densityMatrix_sparse.maxAbsValue();
614 do_output(LOG_CAT_INFO, LOG_AREA_UNDEFINED, "compute_J_by_boxes_sparse, maxDensityMatrixElement = %22.11f",
615 (double)maxDensityMatrixElement);
616
617 std::vector<basis_func_index_pair_struct> basisFuncIndexPairList;
618 int noOfBasisFuncIndexPairs = get_basis_func_pair_list_2el(basisInfo,
619 integralInfo,
620 J_K_params.threshold_J,
621 maxDensityMatrixElement,
622 basisFuncIndexPairList);
623 if(noOfBasisFuncIndexPairs <= 0) {
624 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED,
625 "error in get_basis_func_pair_list, noOfBasisFuncIndexPairs = %i", noOfBasisFuncIndexPairs);
626 return -1;
627 }
628
629 do_output(LOG_CAT_INFO, LOG_AREA_UNDEFINED,
630 "compute_J_by_boxes_sparse: noOfBasisFuncIndexPairs = %i ==> storing %6.2f %% of a full matrix",
631 noOfBasisFuncIndexPairs, (double)100*noOfBasisFuncIndexPairs/((double)n*n));
632
633 // Setup D_list
634 std::vector<ergo_real> D_list(noOfBasisFuncIndexPairs);
635 {
636 std::vector<int> rowind(noOfBasisFuncIndexPairs);
637 std::vector<int> colind(noOfBasisFuncIndexPairs);
638 for(int i = 0; i < noOfBasisFuncIndexPairs; i++)
639 {
640 rowind[i] = basisFuncIndexPairList[i].index_1;
641 colind[i] = basisFuncIndexPairList[i].index_2;
642 }
643 densityMatrix_sparse.get_values(rowind,
644 colind,
645 D_list,
646 permutationHML,
647 permutationHML);
648 }
649
650 std::vector<ergo_real> J_list(noOfBasisFuncIndexPairs);
651
652 output_current_memory_usage(LOG_AREA_UNDEFINED, "Before calling compute_J_by_boxes_linear");
653 if(compute_J_by_boxes_linear(basisInfo,
654 integralInfo,
655 J_K_params,
656 &basisFuncIndexPairList[0],
657 noOfBasisFuncIndexPairs,
658 &D_list[0],
659 &J_list[0],
660 noOfBasisFuncIndexPairs) != 0)
661 {
662 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in compute_J_by_boxes_linear");
663 return -1;
664 }
665 output_current_memory_usage(LOG_AREA_UNDEFINED, "After calling compute_J_by_boxes_linear");
666
667 // Now transfer result from J_list to J
668 {
669 std::vector<int> rowind(noOfBasisFuncIndexPairs);
670 std::vector<int> colind(noOfBasisFuncIndexPairs);
671 for(int i = 0; i < noOfBasisFuncIndexPairs; i++)
672 {
673 rowind[i] = basisFuncIndexPairList[i].index_1;
674 colind[i] = basisFuncIndexPairList[i].index_2;
675 }
676 J.assign_from_sparse(rowind,
677 colind,
678 J_list,
679 permutationHML,
680 permutationHML);
681 output_current_memory_usage(LOG_AREA_UNDEFINED, "After J.assign_from_sparse");
682 }
683
684 timeMeter.print(LOG_AREA_UNDEFINED, "compute_J_by_boxes_sparse total");
685 return 0;
686 }
687
688
689
690 static int
get_CSR_from_symmMatrix(int n,const symmMatrix & A,std::vector<int> const & inversePermutationHML,csr_matrix_struct & CSR)691 get_CSR_from_symmMatrix(int n,
692 const symmMatrix & A,
693 std::vector<int> const & inversePermutationHML,
694 csr_matrix_struct & CSR)
695 {
696 std::vector<int> rowind;
697 std::vector<int> colind;
698 std::vector<ergo_real> values;
699 A.get_all_values(rowind,
700 colind,
701 values,
702 inversePermutationHML,
703 inversePermutationHML);
704 size_t nvalues = values.size();
705 // switch rows and columns if necessary
706 for(size_t i = 0; i < nvalues; i++)
707 {
708 int row = rowind[i];
709 int col = colind[i];
710 if(row > col)
711 {
712 rowind[i] = col;
713 colind[i] = row;
714 }
715 }
716 if(ergo_CSR_create(&CSR,
717 1,
718 n,
719 nvalues,
720 rowind,
721 colind) != 0)
722 {
723 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in ergo_CSR_create.");
724 return -1;
725 }
726 for(size_t i = 0; i < nvalues; i++)
727 {
728 ergo_CSR_add_to_element(&CSR,
729 rowind[i],
730 colind[i],
731 values[i]);
732 }
733 // check that CSR matrix is correct.
734 for(size_t i = 0; i < nvalues; i++)
735 {
736 if(ergo_CSR_get_element(&CSR, rowind[i], colind[i]) != values[i])
737 {
738 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error checking CSR matrix.");
739 return -1;
740 }
741 }
742 return 0;
743 }
744
745
746 static int
get_CSR_from_normalMatrix(int n,const normalMatrix & A,std::vector<int> const & inversePermutationHML,csr_matrix_struct & CSR)747 get_CSR_from_normalMatrix(int n,
748 const normalMatrix & A,
749 std::vector<int> const & inversePermutationHML,
750 csr_matrix_struct & CSR)
751 {
752 std::vector<int> rowind;
753 std::vector<int> colind;
754 std::vector<ergo_real> values;
755 A.get_all_values(rowind,
756 colind,
757 values,
758 inversePermutationHML,
759 inversePermutationHML);
760 size_t nvalues = values.size();
761 if(ergo_CSR_create(&CSR,
762 0,
763 n,
764 nvalues,
765 rowind,
766 colind) != 0)
767 {
768 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in ergo_CSR_create.");
769 return -1;
770 }
771 for(size_t i = 0; i < nvalues; i++) {
772 ergo_CSR_add_to_element(&CSR,
773 rowind[i],
774 colind[i],
775 values[i]);
776 }
777 // check that CSR matrix is correct.
778 for(size_t i = 0; i < nvalues; i++) {
779 if(ergo_CSR_get_element(&CSR, rowind[i], colind[i]) != values[i])
780 {
781 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "Error in get_CSR_from_normalMatrix, error checking CSR matrix.");
782 return -1;
783 }
784 }
785 return 0;
786 }
787
788
789
790 /** Returns the exchange matrix multiplied by 0.5.
791 * To get the correct value multiply K by 2.
792 */
793 int
compute_K_by_boxes_sparse(const BasisInfoStruct & basisInfo,const IntegralInfo & integralInfo,const JK::ExchWeights & CAM_params,const JK::Params & J_K_params,symmMatrix & K,symmMatrix & densityMatrix_sparse,std::vector<int> const & permutationHML,std::vector<int> const & inversePermutationHML)794 compute_K_by_boxes_sparse(const BasisInfoStruct& basisInfo,
795 const IntegralInfo& integralInfo,
796 const JK::ExchWeights & CAM_params,
797 const JK::Params& J_K_params,
798 symmMatrix & K,
799 symmMatrix & densityMatrix_sparse,
800 std::vector<int> const & permutationHML,
801 std::vector<int> const & inversePermutationHML)
802 {
803 Util::TimeMeter timeMeter;
804 int n = basisInfo.noOfBasisFuncs;
805 output_current_memory_usage(LOG_AREA_UNDEFINED, "Beginning of compute_K_by_boxes_sparse");
806
807 csr_matrix_struct dens_CSR;
808 memset(&dens_CSR, 0, sizeof(csr_matrix_struct));
809
810 csr_matrix_struct K_CSR;
811 memset(&K_CSR, 0, sizeof(csr_matrix_struct));
812
813 if(get_CSR_from_symmMatrix(n,
814 densityMatrix_sparse,
815 inversePermutationHML,
816 dens_CSR) != 0)
817 {
818 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in get_CSR_from_symmMatrix");
819 return -1;
820 }
821
822 // set up CSR K matrix for result
823 // note that this call will allocate memory, we must call ergo_CSR_destroy later.
824 if(create_CSR_for_K(basisInfo,
825 integralInfo,
826 J_K_params,
827 &dens_CSR,
828 &K_CSR,
829 1) != 0)
830 {
831 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in compute_K_by_boxes");
832 return -1;
833 }
834
835 densityMatrix_sparse.writeToFile();
836
837 output_current_memory_usage(LOG_AREA_UNDEFINED, "Before calling compute_K_by_boxes");
838 if(compute_K_by_boxes(basisInfo,
839 integralInfo,
840 CAM_params,
841 J_K_params,
842 &K_CSR,
843 &dens_CSR,
844 1) != 0)
845 {
846 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in compute_K_by_boxes");
847 return -1;
848 }
849 output_current_memory_usage(LOG_AREA_UNDEFINED, "After calling compute_K_by_boxes");
850
851 // collect result
852 convert_symm_CSR_to_HML_and_destroy_CSR(K_CSR, K, permutationHML);
853
854 ergo_CSR_destroy(&dens_CSR);
855
856 densityMatrix_sparse.readFromFile();
857
858 timeMeter.print(LOG_AREA_UNDEFINED, "compute_K_by_boxes_sparse total");
859 return 0;
860 }
861
862
863
864 #if 0
865 static void
866 get_fullmatrix_from_normalMatrix(int n,
867 const normalMatrix & A,
868 std::vector<int> const & inversePermutationHML,
869 ergo_real* full)
870 {
871 std::vector<int> rowind;
872 std::vector<int> colind;
873 std::vector<ergo_real> values;
874 A.get_all_values(rowind,
875 colind,
876 values,
877 inversePermutationHML,
878 inversePermutationHML);
879 size_t nvalues = values.size();
880 for(size_t i = 0; i < n*n; i++)
881 full[i] = 0;
882 for(size_t i = 0; i < nvalues; i++) {
883 int p = rowind[i];
884 int q = colind[i];
885 full[p*n+q] = values[i];
886 }
887 }
888 #endif
889
890
891 int
compute_K_by_boxes_sparse_nosymm(const BasisInfoStruct & basisInfo,const IntegralInfo & integralInfo,const JK::ExchWeights & CAM_params,const JK::Params & J_K_params,normalMatrix & K,normalMatrix & densityMatrix_sparse,std::vector<int> const & permutationHML,std::vector<int> const & inversePermutationHML)892 compute_K_by_boxes_sparse_nosymm(const BasisInfoStruct& basisInfo,
893 const IntegralInfo& integralInfo,
894 const JK::ExchWeights & CAM_params,
895 const JK::Params& J_K_params,
896 normalMatrix & K,
897 normalMatrix & densityMatrix_sparse,
898 std::vector<int> const & permutationHML,
899 std::vector<int> const & inversePermutationHML)
900 {
901 Util::TimeMeter timeMeter;
902 int n = basisInfo.noOfBasisFuncs;
903 output_current_memory_usage(LOG_AREA_UNDEFINED, "Beginning of compute_K_by_boxes_sparse");
904
905 csr_matrix_struct dens_CSR;
906 memset(&dens_CSR, 0, sizeof(csr_matrix_struct));
907
908 csr_matrix_struct K_CSR;
909 memset(&K_CSR, 0, sizeof(csr_matrix_struct));
910
911 if(get_CSR_from_normalMatrix(n,
912 densityMatrix_sparse,
913 inversePermutationHML,
914 dens_CSR) != 0)
915 {
916 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in get_CSR_from_normalMatrix");
917 return -1;
918 }
919
920 // set up CSR K matrix for result
921 // note that this call will allocate memory, we must call ergo_CSR_destroy later.
922 if(create_CSR_for_K(basisInfo,
923 integralInfo,
924 J_K_params,
925 &dens_CSR,
926 &K_CSR,
927 0) != 0)
928 {
929 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in compute_K_by_boxes");
930 return -1;
931 }
932
933 densityMatrix_sparse.writeToFile();
934
935 output_current_memory_usage(LOG_AREA_UNDEFINED, "Before calling compute_K_by_boxes");
936 if(compute_K_by_boxes(basisInfo,
937 integralInfo,
938 CAM_params,
939 J_K_params,
940 &K_CSR,
941 &dens_CSR,
942 0) != 0)
943 {
944 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in compute_K_by_boxes");
945 return -1;
946 }
947 output_current_memory_usage(LOG_AREA_UNDEFINED, "After calling compute_K_by_boxes");
948
949 // collect result
950 // Convert K matrix from CSR to HML format.
951 size_t nvalues = ergo_CSR_get_nvalues(&K_CSR);
952 std::vector<int> rowind(nvalues);
953 std::vector<int> colind(nvalues);
954 std::vector<ergo_real> values(nvalues);
955 if(ergo_CSR_get_values(&K_CSR,
956 rowind,
957 colind,
958 values,
959 nvalues) != 0)
960 {
961 do_output(LOG_CAT_ERROR, LOG_AREA_UNDEFINED, "error in ergo_CSR_get_values.");
962 return -1;
963 }
964 // now the information is in the vectors rowind colind values, we can free memory for K_CSR.
965
966 ergo_CSR_destroy(&K_CSR);
967 K.assign_from_sparse(rowind,
968 colind,
969 values,
970 permutationHML,
971 permutationHML);
972
973 ergo_CSR_destroy(&dens_CSR);
974
975 densityMatrix_sparse.readFromFile();
976
977 timeMeter.print(LOG_AREA_UNDEFINED, "compute_K_by_boxes_sparse total");
978 return 0;
979 }
980