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
31 /** @file GetDensFromFock.cc
32 *
33 * @brief Routines for getting density matrix from a given Fock
34 * matrix.
35 *
36 * @author Anastasia Kruchinina <em>responsible</em>
37 * @author Elias Rudberg
38 */
39
40
41
42 #include "output.h"
43 #include <memory.h>
44 #include <stdio.h>
45 #include <math.h>
46 #include <fstream>
47 #include <sstream>
48 #include "utilities.h"
49 #include "matrix_utilities.h"
50 #include "TC2.h"
51 #include "units.h"
52 #include "machine_epsilon.h"
53 #include "AllocatorManager.h"
54
55 #include "densfromf_full.h"
56
57 #include "purification_general.h"
58 #include "purification_sp2.h"
59 #include "purification_sp2acc.h"
60 #include "GetDensFromFock.h"
61
62 typedef generalVector VectorType;
63
64 const int GetDensFromFock::UNDEF_VALUE_UINT = -1;
65 const string GetDensFromFock::NA_STRING = "N/A";
66 const bool GetDensFromFock::BOOL_TRUE = true;
67 const bool GetDensFromFock::BOOL_FALSE = false;
68
69
70 /** Choose which method to use for computing the density matrix from Fock matrix.
71 *
72 * Possible alternatives:
73 * - use recursive expansion
74 * - use diagonalization
75 */
get_dens_from_fock(symmMatrix & Finput,symmMatrix & resultDens,symmMatrix & F_ort_prev)76 int GetDensFromFock::get_dens_from_fock(symmMatrix& Finput,
77 symmMatrix& resultDens,
78 symmMatrix& F_ort_prev)
79 {
80 Util::TimeMeter timeMeterTot;
81
82 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "get_dens_from_fock_general, n = %i, use_diagonalization = %i, use_diag_on_error = %i",
83 n, use_diagonalization, use_diag_on_error);
84 resultEntropyTerm = 0; // In nonzero temperature case, this will be set to nonzero value later.
85
86 std::string allocStatsStr1 = mat::AllocatorManager<ergo_real>::instance().getStatistics();
87 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Before writeAndReadAll(): %s", allocStatsStr1.c_str());
88
89 Util::TimeMeter timeMeterWriteAndReadAll;
90 std::string sizesStr = mat::FileWritable::writeAndReadAll();
91 timeMeterWriteAndReadAll.print(LOG_AREA_DENSFROMF, "FileWritable::writeAndReadAll");
92 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, ((std::string)"writeAndReadAll sizesStr: '" + sizesStr).c_str());
93
94 std::string allocStatsStr2 = mat::AllocatorManager<ergo_real>::instance().getStatistics();
95 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "After writeAndReadAll(): %s", allocStatsStr2.c_str());
96
97 if (noOfOccupiedOrbs == 0)
98 {
99 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "GetDensFromFock::get_dens_from_fock: (noOfOccupiedOrbs == 0), skipping.");
100 resultDens.readFromFile();
101 resultDens.clear();
102 resultDens.writeToFile();
103 return 0;
104 }
105
106
107
108 bool use_diag = false;
109 bool purification_has_failed = false;
110
111 if (use_diagonalization)
112 {
113 use_diag = true;
114 }
115 else
116 {
117 // Try purification
118
119 if (electronicTemperature != 0)
120 {
121 throw std::runtime_error("Error: (electronicTemperature != 0) not implemented for sparse case.");
122 }
123 resultDens.readFromFile();
124 resultDens.clear();
125
126 int puri_res;
127
128 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,
129 "calling get_dens_from_fock_sparse, n = %6i, subspaceErrorLimit = %g",
130 n, (double)subspaceErrorLimit);
131 puri_res = get_dens_from_fock_sparse(Finput,
132 resultDens,
133 F_ort_prev);
134
135
136 if (puri_res != 0)
137 {
138 if (use_diag_on_error)
139 {
140 do_output(LOG_CAT_ERROR, LOG_AREA_DENSFROMF, "get_dens_from_fock (GetDensFromFock class) : error while getting density matrix; trying with diagonalization instead.");
141 use_diag = true;
142 purification_has_failed = true;
143 }
144 else
145 {
146 do_output(LOG_CAT_ERROR, LOG_AREA_DENSFROMF, "get_dens_from_fock (GetDensFromFock class) : error while getting density matrix; aborting.");
147 return -1;
148 }
149 }
150 else
151 {
152 // Purification success!
153 do_output(LOG_CAT_ERROR, LOG_AREA_DENSFROMF, "get_dens_from_fock (GetDensFromFock class) : purification finished OK.");
154 }
155 resultDens.writeToFile();
156 }
157
158
159 if (use_diag)
160 {
161 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "calling get_dens_from_fock_full, n = %i", n);
162
163 std::vector<ergo_real> F_full(n *n);
164 std::vector<ergo_real> S_full(n *n);
165
166 {
167 // Create full matrix versions of F and S
168 normalMatrix *tmpMat;
169 Finput.readFromFile();
170 tmpMat = new normalMatrix(Finput);
171 Finput.writeToFile();
172 tmpMat->fullMatrix(F_full);
173 delete tmpMat;
174 overlapMatrix.readFromFile();
175 tmpMat = new normalMatrix(overlapMatrix);
176 overlapMatrix.writeToFile();
177 tmpMat->fullMatrix(S_full);
178 delete tmpMat;
179 }
180
181 std::vector<ergo_real> densityMatrixFull(n *n);
182
183
184 int number_of_unoccupied_eigenvectors_tmp = 0;
185 int number_of_occupied_eigenvectors_tmp = 0;
186 std::vector<std::vector<ergo_real> > eigVecOCC_tmp;
187 std::vector<std::vector<ergo_real> > eigVecUNOCC_tmp;
188 eigValOCC.clear();
189 eigValUNOCC.clear();
190
191 /* If it is requested to compute some eigenpairs, resize vectors for these eigenpairs. Otherwise, no eigenpairs is returned from get_dens_from_fock_full. */
192 if(output_homo_and_lumo_eigenvectors)
193 {
194 assert(number_of_unoccupied_eigenvectors >= 0);
195 assert(number_of_occupied_eigenvectors >= 0);
196 number_of_occupied_eigenvectors_tmp = number_of_occupied_eigenvectors;
197 number_of_unoccupied_eigenvectors_tmp = number_of_unoccupied_eigenvectors;
198 eigVecOCC_tmp.resize(number_of_occupied_eigenvectors_tmp);
199 eigVecUNOCC_tmp.resize(number_of_unoccupied_eigenvectors_tmp);
200 eigValOCC.resize(number_of_occupied_eigenvectors_tmp);
201 eigValUNOCC.resize(number_of_unoccupied_eigenvectors_tmp);
202 }
203
204 ergo_real gap = -1;
205 assert(factor == 1 || factor == 2);
206
207 if (get_dens_from_fock_full(n,
208 noOfOccupiedOrbs,
209 &densityMatrixFull[0],
210 &F_full[0],
211 &S_full[0],
212 factor,
213 electronicTemperature,
214 resultEntropyTerm,
215 gap,
216 store_all_eigenvalues_to_file,
217 number_of_occupied_eigenvectors_tmp,
218 number_of_unoccupied_eigenvectors_tmp,
219 eigVecOCC_tmp,
220 eigVecUNOCC_tmp,
221 eigValOCC,
222 eigValUNOCC) != 0)
223 {
224 throw std::runtime_error("error in get_dens_from_fock_full");
225 }
226
227 /* If the flag use_diag_on_error is set, we use diagonalization if the purification did not converge. Here we check the gap, and decide, if purification failure is acceptable or there might be some bug. */
228 if (purification_has_failed)
229 {
230 // Accept purification failure only in case gap is very small.
231 ergo_real gapLimit = 1e-4;
232 if (gap > gapLimit)
233 {
234 throw std::runtime_error("Error in GetDensFromFock::get_dens_from_fock: purification failed and (gap > gapLimit). Purification should not fail in such cases; something is wrong.");
235 }
236 }
237
238 resultDens.readFromFile();
239 resultDens.assignFromFull(densityMatrixFull);
240 resultDens.writeToFile();
241
242 if(output_homo_and_lumo_eigenvectors)
243 {
244 // save eigenvectors
245 eigVecUNOCC.clear(); eigVecOCC.clear();
246 eigVecUNOCC.reserve(number_of_unoccupied_eigenvectors_tmp); eigVecOCC.reserve(number_of_occupied_eigenvectors_tmp);
247 generalVector v;
248 for (int i = 0; i < number_of_occupied_eigenvectors_tmp; i++) {
249 v.assign_from_full(eigVecOCC_tmp[i], matrixSizesAndBlocks);
250 eigVecOCC.push_back(v);
251 }
252 for (int i = 0; i < number_of_unoccupied_eigenvectors_tmp; i++) {
253 v.assign_from_full(eigVecUNOCC_tmp[i], matrixSizesAndBlocks);
254 eigVecUNOCC.push_back(v);
255 }
256 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Computed %d occupied eigenpairs and %d unoccupied eigenpairs", number_of_occupied_eigenvectors_tmp, number_of_unoccupied_eigenvectors_tmp);
257 }
258
259 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "get_dens_from_fock_full finished");
260 }
261
262 timeMeterTot.print(LOG_AREA_DENSFROMF, "get_dens_from_fock");
263
264 return 0;
265 }
266
267
268 #ifndef USE_CHUNKS_AND_TASKS
get_eucl_diff_with_adapted_accuracy(int n,const symmMatrixWrap & F_w,const symmMatrixWrap & F_ort_prev_w,ergo_real acc)269 static ergo_real get_eucl_diff_with_adapted_accuracy(int n,
270 const symmMatrixWrap & F_w,
271 const symmMatrixWrap & F_ort_prev_w,
272 ergo_real acc) {
273 // The symmMatrixWrap::eucl_diff() call may be slow, we use a maxIter param to detect if it is a difficult case, and in such cases use a larger acc value.
274 int maxIterForEuclDiff = std::max(n / 10, 500);
275 ergo_real maxEigValMovement_eucl = -1; // Value will be set in try/catch code below.
276 try {
277 Util::TimeMeter timeMeterEuclDiff;
278 maxEigValMovement_eucl = symmMatrixWrap::eucl_diff(F_w, F_ort_prev_w, acc, maxIterForEuclDiff) + acc;
279 timeMeterEuclDiff.print(LOG_AREA_DENSFROMF, "symmMatrixWrap::eucl_diff for maxEigValMovement_eucl ");
280 }
281 catch(...) {
282 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "symmMatrixWrap::eucl_diff() for maxEigValMovement_eucl failed for maxIterForEuclDiff=%d. Calling eucl_diff() again with lower accuracy requirement sqrt(acc).",
283 maxIterForEuclDiff);
284 ergo_real acc2 = template_blas_sqrt(acc);
285 Util::TimeMeter timeMeterEuclDiff;
286 maxEigValMovement_eucl = symmMatrixWrap::eucl_diff(F_w, F_ort_prev_w, acc2) + acc2;
287 timeMeterEuclDiff.print(LOG_AREA_DENSFROMF, "symmMatrixWrap::eucl_diff for maxEigValMovement_eucl ");
288 }
289 return maxEigValMovement_eucl;
290 }
291 #endif
292
293
294 /** Use recursive expansion for computing the density matrix from Fock matrix.
295 *
296 * Construct approximation of the step function by recursive
297 * application of low order polynomials. Sparsity is preserved using
298 * truncation (see J. Chem. Phys. 128, 074106, 2008), which can be
299 * done using spectral, Frobenius or mixed norms (see
300 * J. Comput. Chem. 30.6 (2009): 974-977.).
301 *
302 * Possible alternatives (use_acceleration parameter):
303 * - SP2 recursive expansion
304 * - SP2 accelerated recursive expansion
305 */
get_dens_from_fock_sparse(symmMatrix & F,symmMatrix & resultDens,symmMatrix & F_ort_prev)306 int GetDensFromFock::get_dens_from_fock_sparse(symmMatrix& F,
307 symmMatrix& resultDens,
308 symmMatrix& F_ort_prev
309 )
310 {
311 #ifdef USE_CHUNKS_AND_TASKS
312 if (output_homo_and_lumo_eigenvectors)
313 {
314 throw std::invalid_argument("Error in get_dens_from_fock_sparse (GetDensFromFock class) : computation of eigenvectors is not implemented with Chunks and Tasks.");
315 }
316 if ((leavesSizeMax == UNDEF_VALUE_UINT) || (blocksize == UNDEF_VALUE_UINT))
317 {
318 throw std::invalid_argument("Error in get_dens_from_fock (GetDensFromFock class) : leavesSizeMax and/or blocksize flags are not specified");
319 }
320 size_t NRows = n; // defined in GetDensFromFock.h
321 size_t NCols = n;
322 #if defined(USE_CHUNKS_AND_TASKS_BSM) || defined(USE_CHUNKS_AND_TASKS_HBSM)
323 ParamsType params(leavesSizeMax, blocksize, NRows, NCols);
324 #else
325 ParamsType params(leavesSizeMax, NRows, NCols);
326 #endif
327 #else
328 ParamsType params;
329 #endif
330
331
332 Util::TimeMeter timeMeterTot;
333 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "get_dens_from_fock_sparse() start!");
334
335 if (use_new_stopping_criterion)
336 {
337 // If we use the new stopping criterion, then we should stop
338 // when the idempotency error cannot be decreased
339 // significantly anymore. The error in the subspace
340 // accumulates during the iterations, thus we can expect that
341 // it would be larger than the error in eigenvalues.
342 eigvalueErrorLimit = subspaceErrorLimit;
343 }
344
345 #ifdef USE_CHUNKS_AND_TASKS
346 if ((truncationNormPurification == mat::euclNorm) || (truncationNormPurification == mat::mixedNorm))
347 {
348 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "get_dens_from_fock_sparse(): set truncation "
349 "norm to mat::frobNorm since spectral and mixed norms are not implemented");
350 set_truncationNormPurification(mat::frobNorm);
351 }
352 if ((stopCriterionNormPurification == mat::euclNorm) || (stopCriterionNormPurification == mat::mixedNorm))
353 {
354 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "get_dens_from_fock_sparse(): set stopping criterion "
355 "norm to mat::frobNorm since spectral and mixed norms are is not implemented");
356 set_stopCriterionNormPurification(mat::frobNorm);
357 }
358 #endif
359
360 // Select tolerated errors in the occupied subspace for the three truncations
361 // and for purification.
362 ergo_real subspaceThr_1 = 0.1 * subspaceErrorLimit;
363 ergo_real subspaceThr_Puri = 0.7 * subspaceErrorLimit;
364 ergo_real subspaceThr_2 = 0.1 * subspaceErrorLimit;
365 ergo_real subspaceThr_3 = 0.1 * subspaceErrorLimit;
366
367 // Select tolerated errors in eigenvalues
368 ergo_real eigvalueThr_Puri = 0.7 * eigvalueErrorLimit;
369 ergo_real eigvalueThr_2 = 0.15 * eigvalueErrorLimit;
370 ergo_real eigvalueThr_3 = 0.15 * eigvalueErrorLimit;
371
372 symmMatrix F_tmp(F);
373 F_tmp.readFromFile(); // F stay written in the file
374 std::string allocStatsStr3 = mat::AllocatorManager<ergo_real>::instance().getStatistics();
375 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "After F.readFromFile(): %s", allocStatsStr3.c_str());
376
377
378 symmMatrixWrap F_w;
379 transform_matrix_from_to(F_tmp, F_w, params);
380
381
382
383 std::string allocStatsStr4 = mat::AllocatorManager<ergo_real>::instance().getStatistics();
384 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "After creating wrapper MatrixType F_w(F): %s", allocStatsStr4.c_str());
385 F_tmp.clear(); // we do not need this matrix anymore, it is saved in the wrapper
386
387 // transform to the orthogonal basis
388 {
389 triangMatrix invCholFactor_tmp(invCholFactor);
390 invCholFactor_tmp.readFromFile();
391 output_current_memory_usage(LOG_AREA_DENSFROMF, "In get_dens_from_fock_sparse, before F_w to orthogonal basis");
392
393 Util::TimeMeter timeMeterFortTransf;
394
395 triangMatrixWrap invCholFactor_tmp_w;
396 transform_matrix_from_to(invCholFactor_tmp, invCholFactor_tmp_w, params);
397 invCholFactor_tmp.clear();
398
399 F_w = transpose(invCholFactor_tmp_w) * F_w * invCholFactor_tmp_w;
400
401 timeMeterFortTransf.print(LOG_AREA_DENSFROMF, " F_w to orthogonal basis");
402 output_current_memory_usage(LOG_AREA_DENSFROMF, "In get_dens_from_fock_sparse, after "
403 "F_ort = tr(Z) * F_S * Z");
404 }
405
406 // Now F_w contains F_ort.
407
408 //Compare F to F_ort_prev to check how far eigenvalues move.
409 F_ort_prev.readFromFile();
410 output_current_memory_usage(LOG_AREA_DENSFROMF,
411 "In get_dens_from_fock_sparse, after F_ort_prev.readFromFile()");
412
413
414 symmMatrixWrap F_ort_prev_w;
415 transform_matrix_from_to(F_ort_prev, F_ort_prev_w, params);
416
417 output_current_memory_usage(LOG_AREA_DENSFROMF,
418 "In get_dens_from_fock_sparse, after symmMatrixWrap F_ort_prev_w(F_ort_prev)");
419 F_ort_prev.clear();
420
421 // intervals will be expanded to make sure that they contain the HOMO and LUMO eigenvalues of F_ort
422 intervalType lumoInterval_F_ort_prev_expanded;
423 intervalType homoInterval_F_ort_prev_expanded;
424
425 {
426 ergo_real maxEigValMovement_frob = symmMatrixWrap::frob_diff(F_w, F_ort_prev_w);
427 output_current_memory_usage(LOG_AREA_DENSFROMF,
428 "In get_dens_from_fock_sparse, after getting maxEigValMovement_frob ");
429
430 #ifndef USE_CHUNKS_AND_TASKS
431 ergo_real acc = template_blas_sqrt(get_machine_epsilon());
432 Util::TimeMeter timeMeterMixedDiff;
433 ergo_real maxEigValMovement_mixed = symmMatrixWrap::mixed_diff(F_w, F_ort_prev_w, acc) + acc;
434 timeMeterMixedDiff.print(LOG_AREA_DENSFROMF, "MatrixType::mixed_diff for maxEigValMovement_mixed");
435 output_current_memory_usage(LOG_AREA_DENSFROMF,
436 "In get_dens_from_fock_sparse, after getting maxEigValMovement_mixed");
437 ergo_real maxEigValMovement_eucl = get_eucl_diff_with_adapted_accuracy(n, F_w, F_ort_prev_w, acc);
438 output_current_memory_usage(LOG_AREA_DENSFROMF,
439 "In get_dens_from_fock_sparse, after getting maxEigValMovement_eucl ");
440 #endif
441
442 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "maxEigValMovement_frob = %22.11f", (double)maxEigValMovement_frob);
443 #ifndef USE_CHUNKS_AND_TASKS
444 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "maxEigValMovement_mixed = %22.11f", (double)maxEigValMovement_mixed);
445 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "maxEigValMovement_eucl = %22.11f", (double)maxEigValMovement_eucl);
446 #endif
447 // Increase HOMO/LUMO intervals so that they for sure contain the HOMO and LUMO eigenvalues of F_ort
448
449 // Anastasia comment: it may happen for very small cases that
450 // bounds for homo and lumo are very good and intervals are
451 // actually empty sets, then the increase() function will throw
452 // exception
453 if (homoInterval_F_ort_prev.low() >= homoInterval_F_ort_prev.upp())
454 {
455 homoInterval_F_ort_prev = intervalType(homoInterval_F_ort_prev.low(), homoInterval_F_ort_prev.low() + 1e-10);
456 }
457 if (lumoInterval_F_ort_prev.low() >= lumoInterval_F_ort_prev.upp())
458 {
459 lumoInterval_F_ort_prev = intervalType(lumoInterval_F_ort_prev.upp() - 1e-10, lumoInterval_F_ort_prev.upp());
460 }
461
462 #ifdef USE_CHUNKS_AND_TASKS
463 ergo_real maxEigValMovement = maxEigValMovement_frob;
464 #else
465 ergo_real maxEigValMovement = maxEigValMovement_eucl;
466 #endif
467
468 lumoInterval_F_ort_prev_expanded = lumoInterval_F_ort_prev;
469 homoInterval_F_ort_prev_expanded = homoInterval_F_ort_prev;
470
471 homoInterval_F_ort_prev_expanded.increase(maxEigValMovement);
472 lumoInterval_F_ort_prev_expanded.increase(maxEigValMovement);
473
474
475 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "lumo before truncation: [ %.12lf , %.12lf ]", (double)lumoInterval_F_ort_prev_expanded.low(), (double)lumoInterval_F_ort_prev_expanded.upp());
476 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "homo before truncation: [ %.12lf , %.12lf ]", (double)homoInterval_F_ort_prev_expanded.low(), (double)homoInterval_F_ort_prev_expanded.upp());
477 }
478
479
480
481 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Truncate matrix F and update intervals for homo and lumo.");
482
483 /* EMANUEL COMMENT:
484 * - truncationNorm can be set to any of
485 * mat::frobNorm, mat::mixedNorm, or mat::euclNorm
486 * The best choice depends on a trade-off between spending
487 * time in truncation and in matrix-matrix multiplication.
488 */
489 mat::normType truncationNorm = truncationNormPurification;
490
491
492 // Now, we will truncate F (and update eigenvalue intervals):
493 ergo_real truncError_1;
494 {
495 ergo_real gapMin = lumoInterval_F_ort_prev_expanded.low() - homoInterval_F_ort_prev_expanded.upp();
496 ergo_real gapMax = lumoInterval_F_ort_prev_expanded.upp() - homoInterval_F_ort_prev_expanded.low();
497 ergo_real threshold_1;
498 // We consider the gap to be accurately known if the uncertainty is at most 10 %
499 if ((gapMin > 0) && ((gapMax - gapMin) / gapMin < 0.1))
500 {
501 // Gap is accurately known: we use gapMin
502 threshold_1 = subspaceThr_1 * gapMin / (1 + subspaceThr_1);
503 }
504 else
505 {
506 // Gap is not accurately known. To avoid choosing a very tight
507 // threshold value due to a small lower bound for the gap, we
508 // use the largest of 'gap_expected_lower_bound' and calculated
509 // 'gapMin':
510 threshold_1 = gapMin > gap_expected_lower_bound ?
511 subspaceThr_1 * gapMin / (1 + subspaceThr_1) :
512 subspaceThr_1 * gap_expected_lower_bound / (1 + subspaceThr_1);
513 }
514
515 double nnzF_before_trunc_pc = (double)F_w.nnz() * 100 / ((double)n * n);
516 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Truncating F_ort ( %s ), selected threshold = %10.6g",
517 mat::getNormTypeString(truncationNorm).c_str(), (double)threshold_1);
518 Util::TimeMeter timeMeterFThresh;
519 #ifdef USE_CHUNKS_AND_TASKS
520 truncError_1 = F_w.thresh_frob(threshold_1);
521 #else
522 truncError_1 = F_w.thresh(threshold_1, truncationNorm);
523 #endif
524 double nnzF_after_trunc_pc = (double)F_w.nnz() * 100 / ((double)n * n);
525 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,
526 "Truncated F_ort ( %s ), selected threshold = %10.6g, returned error = %10.6g, nnz before = %3.4f %%, nnz after = %3.4f %%",
527 mat::getNormTypeString(truncationNorm).c_str(), (double)threshold_1, (double)truncError_1, nnzF_before_trunc_pc, nnzF_after_trunc_pc);
528 timeMeterFThresh.print(LOG_AREA_DENSFROMF, "Truncation of F_ort");
529 puri_stats[stats_prefix + "nnz_percentage_F_ort"] = nnzF_after_trunc_pc;
530
531 // Increase HOMO and LUMO intervals so that they contain the eigenvalues of the truncated matrix:
532 homoInterval_F_ort_prev_expanded.increase(truncError_1);
533 lumoInterval_F_ort_prev_expanded.increase(truncError_1);
534 }
535
536
537 F_ort_prev_w.clear();
538
539 transform_matrix_from_to(F_w, F_ort_prev, params);
540
541 F_ort_prev.writeToFile();
542
543
544
545 // The HOMO and LUMO intervals now contain the HOMO and LUMO
546 // eigenvalues of F_ort_prev but improved values will hopefully be
547 // calculated in purification.
548
549 intervalType homoIntervalSaved = homoInterval_F_ort_prev_expanded;
550 intervalType lumoIntervalSaved = lumoInterval_F_ort_prev_expanded;
551
552
553 #ifdef USE_CHUNKS_AND_TASKS
554 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen CHT Wrapper");
555 #else
556 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen Ergo Wrapper");
557 #endif
558
559
560 PurificationGeneral<symmMatrixWrap> *Puri; // abstract class
561
562 if (use_acceleration)
563 {
564 Puri = new Purification_sp2acc<symmMatrixWrap>();
565 }
566 else
567 {
568 Puri = new Purification_sp2<symmMatrixWrap>();
569 }
570
571 mat::Gblas::timekeeping = true;
572 mat::Gblas::time = 0;
573
574 Puri->initialize(F_w,
575 lumoIntervalSaved,
576 homoIntervalSaved,
577 maxMul,
578 subspaceThr_Puri,
579 eigvalueThr_Puri,
580 use_new_stopping_criterion, // 1 = new, 0 = old
581 truncationNorm,
582 stopCriterionNormPurification,
583 noOfOccupiedOrbs
584 );
585
586
587 if ((!output_homo_and_lumo_eigenvectors) || (eigenvectors_method == NA_STRING))
588 {
589 eigVecOCC.clear();
590 eigVecUNOCC.clear();
591 Puri->set_number_of_eigenvectors_to_compute(
592 0 /*number_of_occupied_eigenvectors*/,
593 0 /*number_of_unoccupied_eigenvectors*/);
594 }
595 else
596 {
597 Puri->set_number_of_eigenvectors_to_compute(
598 number_of_occupied_eigenvectors,
599 number_of_unoccupied_eigenvectors);
600
601 // used only of projection method is chosen
602 Puri->set_go_back_X_iter_proj_method(go_back_X_iter_proj_method);
603 Puri->set_jump_over_X_iter_proj_method(jump_over_X_iter_proj_method);
604
605
606 const std::vector<VectorType> & eigVecOCCref = eigVecOCC;
607 const std::vector<VectorType> & eigVecUNOCCref = eigVecUNOCC;
608
609 Puri->set_eigenvectors_params(eigenvectors_method,
610 eigenvectors_iterative_method,
611 eigensolver_accuracy,
612 eigensolver_maxiter,
613 SCF_step,
614 use_prev_vector_as_initial_guess,
615 try_eigv_on_next_iteration_if_fail,
616 eigVecOCCref, eigVecUNOCCref // used as initial guesses if use_prev_vector_as_initial_guess is set
617 );
618
619 if (puri_compute_eigv_in_each_iteration)
620 {
621 Puri->set_compute_eigenvectors_in_each_iteration();
622 }
623
624
625
626 if (run_shift_and_square_method_on_F)
627 {
628 // COMPUTE EIGENVALUES FOR F (FOR COMPARISON)
629 int eigsolver_maxiter = 5000;
630 Puri->compute_eigenvectors_without_diagonalization_on_F(F_w, eigsolver_maxiter);
631 }
632 }
633
634
635 // save Hamiltonian in this permutation if needed to compute eigenvectors
636 // for plotting error in matlab
637 if ((save_permuted_F_matrix_in_bin) && (SCF_step >= 1))
638 {
639 symmMatrix Y;
640 // get structure
641 mat::SizesAndBlocks rows;
642 F_ort_prev.readFromFile();
643 F_ort_prev.getRows(rows);
644 F_ort_prev.writeToFile();
645 Y.resetSizesAndBlocks(rows,rows);
646 transform_matrix_from_to(F_w, Y, params);
647 vector<int> Itmp, I, Jtmp, J;
648 vector<real> Vtmp, V;
649 Y.get_all_values(Itmp, Jtmp, Vtmp);
650
651 size_t nnz = 0;
652 // Count nonzeros
653 for (size_t i = 0; i < Itmp.size(); i++)
654 {
655 nnz += (Vtmp[i] != 0);
656 }
657
658 I.reserve(nnz);
659 J.reserve(nnz);
660 V.reserve(nnz);
661 // Extract nonzeros
662 for (size_t i = 0; i < Itmp.size(); i++)
663 {
664 if (Vtmp[i] != 0)
665 {
666 I.push_back(Itmp[i]);
667 J.push_back(Jtmp[i]);
668 V.push_back(Vtmp[i]);
669 }
670 }
671 ostringstream name;
672 name << "F.bin";
673 write_matrix_to_bin(name.str().c_str(), I, J, V, Y.get_nrows());
674 name.str("");
675 }
676
677
678 F_w.clear();
679
680 output_current_memory_usage(LOG_AREA_DENSFROMF, "Before Puri->PurificationStart()");
681 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,
682 "calling Puri->PurificationStart(), number of threads = %i, trunc norm '%s'",
683 mat::Params::getNProcs(), mat::getNormTypeString(truncationNorm).c_str());
684 mat::FileWritable::resetStats();
685 time_t puriStartWallTime;
686 time(&puriStartWallTime);
687
688 try
689 {
690 Util::TimeMeter timeMeterPurification;
691 Puri->PurificationStart();
692 timeMeterPurification.print(LOG_AREA_DENSFROMF, "Puri->PurificationStart()");
693 //Puri->info.print_collected_info();
694 }
695 catch (...)
696 {
697 if (Puri != NULL)
698 {
699 delete Puri;
700 }
701 throw;
702 }
703
704
705 {
706 std::stringstream ss;
707 ss << "Accumulated wall times for writeToFile in PurificationStart() : " << mat::FileWritable::getStatsTimeWrite();
708 do_output(LOG_CAT_TIMINGS, LOG_AREA_DENSFROMF, ss.str().c_str());
709 }
710 {
711 std::stringstream ss;
712 ss << "Accumulated wall times for readFromFile in PurificationStart() : " << mat::FileWritable::getStatsTimeRead();
713 do_output(LOG_CAT_TIMINGS, LOG_AREA_DENSFROMF, ss.str().c_str());
714 }
715 {
716 std::stringstream ss;
717 ss << "Accumulated wall times for copy and assign in PurificationStart() : " << mat::FileWritable::getStatsTimeCopyAndAssign();
718 do_output(LOG_CAT_TIMINGS, LOG_AREA_DENSFROMF, ss.str().c_str());
719 }
720
721
722 {
723 std::stringstream ss;
724 ss << "Number of calls to writeToFile in PurificationStart() : " << mat::FileWritable::getStatsCountWrite();
725 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, ss.str().c_str());
726 }
727 {
728 std::stringstream ss;
729 ss << "Number of calls to readFromFile in PurificationStart() : " << mat::FileWritable::getStatsCountRead();
730 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, ss.str().c_str());
731 }
732 {
733 std::stringstream ss;
734 ss << "Number of calls to FileWritable copy and assign in PurificationStart() : " << mat::FileWritable::getStatsCountCopyAndAssign();
735 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, ss.str().c_str());
736 }
737
738 do_output(LOG_CAT_TIMINGS, LOG_AREA_DENSFROMF,
739 "mat::Gblas::time after purification : %12.6f",
740 (double)mat::Gblas::time);
741
742 output_current_memory_usage(LOG_AREA_DENSFROMF, "After Puri->PurificationStart()");
743
744 symmMatrixWrap D_w(Puri->X);
745 Puri->clear(); // delete matrices from Puri
746 intervalType homoIntervalNew = intervalType(Puri->info.homo_estim_low_F, Puri->info.homo_estim_upp_F);
747 intervalType lumoIntervalNew = intervalType(Puri->info.lumo_estim_low_F, Puri->info.lumo_estim_upp_F);
748
749
750 if (plot_puri_results)
751 {
752 // plot results
753 ostringstream name;
754 name << "puri_out_error_" << SCF_step << plot_puri_results_str << ".m";
755 Puri->gen_matlab_file_norm_diff(name.str().c_str());
756 name.str("");
757
758 name << "puri_out_threshold_" << SCF_step << plot_puri_results_str << ".m";
759 Puri->gen_matlab_file_threshold(name.str().c_str());
760 name.str("");
761
762 name << "puri_out_nnz_" << SCF_step << plot_puri_results_str << ".m";
763 Puri->gen_matlab_file_nnz(name.str().c_str());
764 name.str("");
765
766 name << "puri_out_eigs_" << SCF_step << plot_puri_results_str << ".m";
767 Puri->gen_matlab_file_eigs(name.str().c_str());
768 name.str("");
769
770 name << "puri_out_time_" << SCF_step << plot_puri_results_str << ".m";
771 Puri->gen_matlab_file_time(name.str().c_str());
772 name.str("");
773 }
774 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Created .m files with results of the purification");
775
776 if (!Puri->info.converged)
777 {
778 if (! ignore_purification_failure)
779 {
780 do_output(LOG_CAT_ERROR, LOG_AREA_DENSFROMF,
781 "Error in purification: Puri->info.converged() "
782 "returned false.");
783 Puri->info.print_collected_info();
784 return -1;
785 }
786 else
787 {
788 do_output(LOG_CAT_WARNING, LOG_AREA_DENSFROMF, "Purification did NOT converge, ignoring.");
789 }
790 }
791 else
792 {
793 ergo_real acc_error = (ergo_real)Puri->info.accumulated_error_subspace;
794 if (acc_error == -1) // if we do not know that gap
795 {
796 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,
797 "Purification converged OK");
798 }
799 else
800 {
801 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,
802 "Purification converged OK, subspaceError <= %22.11f", acc_error);
803 }
804 }
805
806
807 Puri->extract_computed_eigenpairs(eigVecUNOCC, eigVecOCC, eigValUNOCC, eigValOCC);
808 assert(eigVecUNOCC.size() == eigValUNOCC.size());
809 assert(eigVecOCC.size() == eigValOCC.size());
810
811 // compute eigenvectors
812 if ((! eigVecOCC.empty()) || (! eigVecUNOCC.empty()))
813 {
814 triangMatrix invCholFactor_tmp(invCholFactor);
815 invCholFactor_tmp.readFromFile();
816 /* here: if eigenvector is not computed, it is empty, not NULL */
817 if (Puri->info.lumo_eigenvector_is_computed)
818 {
819 assert(! eigVecUNOCC.empty() );
820 // printf("1 %d %lf %d %lf\n",
821 // Puri->info.lumo_eigenvector_is_computed_in_iter,
822 // (double)Puri->info.eigValLUMO,
823 // Puri->info.lumo_eigensolver_iter,
824 // Puri->info.lumo_eigensolver_time);
825
826 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "LUMO eigenvector is computed.");
827 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Perform congruence transformation.");
828 // perform congruence transformation
829 eigVecUNOCC[0] = invCholFactor_tmp * eigVecUNOCC[0];
830
831
832 /* if we requested to compute more than one unoccupied eigenpair and if all of the requested eigenvectors are computed, then Puri->info.lumo_eigenvector_is_computed will be true */
833 if(number_of_unoccupied_eigenvectors > 1)
834 {
835 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "More unoccupied eigenvectors are computed. Perform congruence transformation.");
836 for (size_t i = 1; i < eigVecUNOCC.size(); i++) { // note, counter starts from 1 since we already performed the transformation on the lumo eigenvector
837 // perform congruence transformation
838 eigVecUNOCC[i] = invCholFactor_tmp * eigVecUNOCC[i];
839 }
840 }
841
842 }
843
844
845 if (Puri->info.homo_eigenvector_is_computed)
846 {
847 assert(! eigVecOCC.empty() );
848 // printf("2 %d %lf %d %lf\n",
849 // Puri->info.homo_eigenvector_is_computed_in_iter,
850 // (double)Puri->info.eigValHOMO,
851 // Puri->info.homo_eigensolver_iter,
852 // Puri->info.homo_eigensolver_time);
853
854 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "HOMO eigenvector is computed.");
855 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Perform congruence transformation.");
856 // perform congruence transformation
857 eigVecOCC[0] = invCholFactor_tmp * eigVecOCC[0];
858
859 /* if we requested to compute more than one occupied eigenpair and if all of the requested eigenvectors are computed, then Puri->info.homo_eigenvector_is_computed will be true */
860 if(number_of_occupied_eigenvectors > 1)
861 {
862 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "More occupied eigenvectors are computed. Perform congruence transformation.");
863 for (size_t i = 1; i < eigVecOCC.size(); i++) { // note, counter starts from 1 since we already performed the transformation on the homo eigenvector
864 // perform congruence transformation
865 eigVecOCC[i] = invCholFactor_tmp * eigVecOCC[i];
866 }
867 }
868
869 }
870 } // Note: invCholFactor_tmp goes out of scope
871
872
873 if (intervalType::intersect(lumoInterval_F_ort_prev_expanded, lumoIntervalNew).empty())
874 {
875 do_output(LOG_CAT_WARNING, LOG_AREA_DENSFROMF,
876 "The intersection of lumoInterval_F_ort_prev_expanded and lumoIntervalNew is empty set!");
877 }
878
879 if (intervalType::intersect(homoInterval_F_ort_prev_expanded, homoIntervalNew).empty())
880 {
881 do_output(LOG_CAT_WARNING, LOG_AREA_DENSFROMF,
882 "The intersection of homoInterval_F_ort_prev_expanded and homoIntervalNew is empty set!");
883 }
884
885 // Save the improved HOMO/LUMO intervals of F_ort:
886 homoInterval_F_ort_prev = homoIntervalNew;
887 lumoInterval_F_ort_prev = lumoIntervalNew;
888
889 // Calculate HOMO_LUMO intervals of Finput. We need to expand
890 // the F_ort intervals due to the truncation done earlier.
891 homoInterval_Finput = homoInterval_F_ort_prev;
892 lumoInterval_Finput = lumoInterval_F_ort_prev;
893
894 // Anastasia comment:
895 // it may happen that bounds for homo and lumo are very good and intervals are actually empty sets,
896 // then the increase() function will throw exception
897 if (homoInterval_Finput.low() >= homoInterval_Finput.upp())
898 {
899 homoInterval_Finput = intervalType(homoInterval_Finput.low(), homoInterval_Finput.low() + 1e-10);
900 }
901 if (lumoInterval_Finput.low() >= lumoInterval_Finput.upp())
902 {
903 lumoInterval_Finput = intervalType(lumoInterval_Finput.upp() - 1e-10, lumoInterval_Finput.upp());
904 }
905
906 homoInterval_Finput.increase(truncError_1);
907 lumoInterval_Finput.increase(truncError_1);
908
909 // Output info about gap.
910 ergo_real gapMin = lumoInterval_F_ort_prev.low() - homoInterval_F_ort_prev.upp();
911 ergo_real gapMax = lumoInterval_F_ort_prev.upp() - homoInterval_F_ort_prev.low();
912 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,
913 "E(LUMO) - E(HOMO) >= %22.11f = %22.11f eV",
914 (double)gapMin, (double)gapMin / UNIT_one_eV);
915 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,
916 "E(LUMO) - E(HOMO) <= %22.11f = %22.11f eV",
917 (double)gapMax, (double)gapMax / UNIT_one_eV);
918 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,
919 "HOMO interval : [ %17.12f %17.12f ]",
920 (double)homoInterval_F_ort_prev.low(), (double)homoInterval_F_ort_prev.upp());
921 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,
922 "LUMO interval : [ %17.12f %17.12f ]",
923 (double)lumoInterval_F_ort_prev.low(), (double)lumoInterval_F_ort_prev.upp());
924
925 puri_stats[stats_prefix + "HOMO_LUMO_gap_lo_eV"] = gapMin / UNIT_one_eV;
926 puri_stats[stats_prefix + "HOMO_LUMO_gap_hi_eV"] = gapMax / UNIT_one_eV;
927
928
929 // we do not need Puri anymore, then delete it
930 delete Puri;
931
932
933 // Check trace of resulting density matrix
934 ergo_real trace = D_w.trace();
935 ergo_real wantedTrace = noOfOccupiedOrbs;
936 ergo_real traceError = trace - wantedTrace;
937 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,
938 "Trace of resulting density matrix is %22.11f, error is %18.14f.",
939 (double)trace, (double)traceError);
940 // Check that relative error in trace is not unreasonably large
941 if (!ignore_purification_failure && (wantedTrace > 0) && (template_blas_fabs(traceError) > wantedTrace / 4))
942 {
943 throw std::runtime_error("Error in get_dens_from_fock_sparse (GetDensFromFock class): traceError is unreasonably large; seems like something went very wrong.");
944 }
945
946
947 // Do truncation to speed up following multiplication operation.
948 ergo_real threshold_2 = subspaceThr_2 * (1 - 2 * eigvalueThr_Puri) / (1 + subspaceThr_2);
949 // Make sure that eigenvalue movement is not too large:
950 threshold_2 = eigvalueThr_2 < threshold_2 ? eigvalueThr_2 : threshold_2;
951 size_t nnzD_before_trunc = D_w.nnz();
952 double nnzD_before_trunc_pc = (double)nnzD_before_trunc * 100 / ((double)n * n);
953 #ifdef USE_CHUNKS_AND_TASKS
954 ergo_real truncError_2 = D_w.thresh_frob(threshold_2);
955 #else
956 ergo_real truncError_2 = D_w.thresh(threshold_2, truncationNorm);
957 #endif
958 size_t nnzD_after_trunc = D_w.nnz();
959 double nnzD_after_trunc_pc = (double)nnzD_after_trunc * 100 / ((double)n * n);
960 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,
961 "Truncated D_ort ( %s ), selected threshold = %10.6g, "
962 "returned error = %10.6g, nnz before = %lu <-> %3.4f %%, nnz after = %lu <-> %3.4f %%",
963 mat::getNormTypeString(truncationNorm).c_str(), (double)threshold_2,
964 (double)truncError_2, nnzD_before_trunc, nnzD_before_trunc_pc, nnzD_after_trunc, nnzD_after_trunc_pc);
965 puri_stats[stats_prefix + "nnz_percentage_D_ort"] = nnzD_after_trunc_pc;
966
967 {
968 triangMatrix invCholFactor_tmp(invCholFactor);
969 invCholFactor_tmp.readFromFile();
970 output_current_memory_usage(LOG_AREA_DENSFROMF, "Before D_w.to_nonnorm_basis");
971 Util::TimeMeter timeMeterWriteAndReadAll;
972 std::string sizesStr = mat::FileWritable::writeAndReadAll();
973 timeMeterWriteAndReadAll.print(LOG_AREA_DENSFROMF, "FileWritable::writeAndReadAll");
974 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, ((std::string)"writeAndReadAll sizesStr: '" + sizesStr).c_str());
975 Util::TimeMeter timeMeterDortTransf;
976
977 triangMatrixWrap invCholFactor_tmp_w;
978 transform_matrix_from_to(invCholFactor_tmp, invCholFactor_tmp_w, params);
979 invCholFactor_tmp.clear();
980
981 D_w = invCholFactor_tmp_w * D_w * transpose(invCholFactor_tmp_w);
982
983 timeMeterDortTransf.print(LOG_AREA_DENSFROMF, "D_w to non-orthogonal basis");
984 output_current_memory_usage(LOG_AREA_DENSFROMF, "After D_w to non-orthogonal basis [D_S = Z * D_ort * ZT]");
985
986 #ifndef USE_CHUNKS_AND_TASKS // eucl_thresh is not implemented for CHT
987 // Do truncation again, to reduce memory usage.
988 ergo_real threshold_3 = subspaceThr_3 * (1 - 2 * eigvalueThr_Puri - 2 * truncError_2) / (1 + subspaceThr_3);
989
990 //Make sure that eigenvalue movement is not too large:
991 threshold_3 = eigvalueThr_3 < threshold_3 ? eigvalueThr_3 : threshold_3;
992
993 //Do truncation, taking into account that we are in 'non-orthogonal basis', passing invCholFactor to thresh
994 double nnzD_S_before_trunc_pc = (double)D_w.nnz() * 100 / ((double)n * n);
995 ergo_real truncError_3 = D_w.eucl_thresh(threshold_3, &invCholFactor_tmp_w);
996 double nnzD_S_after_trunc_pc = (double)D_w.nnz() * 100 / ((double)n * n);
997 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,
998 "Truncated D_S (eucl with Z), selected threshold = %10.6g, returned error = %10.6g, nnz before = %3.4f %%, nnz after = %3.4f %%",
999 (double)threshold_3, (double)truncError_3, nnzD_S_before_trunc_pc, nnzD_S_after_trunc_pc);
1000 puri_stats[stats_prefix + "nnz_percentage_D_S"] = nnzD_S_after_trunc_pc;
1001 #endif
1002 }
1003
1004 assert(factor == 1 || factor == 2);
1005 D_w *= factor;
1006 transform_matrix_from_to(D_w, resultDens, params);
1007 D_w.clear(); // clean wrapper
1008
1009 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "get_dens_from_fock_sparse ending OK");
1010 timeMeterTot.print(LOG_AREA_DENSFROMF, "get_dens_from_fock_sparse");
1011
1012 return 0;
1013 }
1014
1015
1016 /** Function save all needed data to files in order to repeat
1017 * recursive expansion in a desired SCF cycle later. The purpose of
1018 * the function is mainly testing.
1019 */
create_checkpoint(symmMatrix & Finput,symmMatrix & F_ort_prev,generalVector * eigVecLUMO,generalVector * eigVecHOMO,std::string IDstr)1020 void GetDensFromFock::create_checkpoint(symmMatrix& Finput,
1021 symmMatrix& F_ort_prev,
1022 generalVector *eigVecLUMO,
1023 generalVector *eigVecHOMO,
1024 std::string IDstr)
1025 {
1026 std::ostringstream checkpoint_ID;
1027 checkpoint_ID << IDstr << "_" << SCF_step;
1028 std::cout << "Create checkpoint with ID = " << checkpoint_ID.str() << std::endl;
1029
1030 /* Save input data */
1031
1032 std::ostringstream name;
1033 name << filenameFinput << "_" << checkpoint_ID.str() << ".bin";
1034 Finput.copyToFile(name.str().c_str());
1035
1036 name.clear();
1037 name.str("");
1038 name << filenameF_ort_prev << "_" << checkpoint_ID.str() << ".bin";
1039 F_ort_prev.copyToFile(name.str().c_str());
1040
1041
1042 /* Save class members */
1043
1044 // matrices
1045
1046 // Overlap matrix (written to file)
1047 name.clear();
1048 name.str("");
1049 name << filenameOverlap << "_" << checkpoint_ID.str() << ".bin";
1050 overlapMatrix.copyToFile(name.str().c_str());
1051
1052
1053 // Inverse Cholesky factor (written to file)
1054 name.clear();
1055 name.str("");
1056 name << filenameinvCholFactor << "_" << checkpoint_ID.str() << ".bin";
1057 invCholFactor.copyToFile(name.str().c_str());
1058
1059
1060
1061 // save all numbers of basic arithmetic types
1062 name.clear();
1063 name.str("");
1064 name << file_for_basic_types << "_" << checkpoint_ID.str() << ".data";
1065 ofstream file_basic(name.str().c_str(), ios::out | ios::app);
1066 if (!file_basic.is_open())
1067 {
1068 throw std::runtime_error("GetDensFromFock::create_checkpoint unable open file for writing.");
1069 }
1070
1071
1072 file_basic << SCF_step << "\n";
1073 file_basic << use_diagonalization << "\n";
1074 file_basic << use_purification << "\n";
1075 file_basic << store_all_eigenvalues_to_file << "\n";
1076 file_basic << try_eigv_on_next_iteration_if_fail << "\n";
1077 file_basic << (double)electronicTemperature << "\n";
1078 file_basic << std::setprecision(8) << (double)gap_expected_lower_bound << "\n";
1079 file_basic << std::setprecision(8) << (double)eigvalueErrorLimit << "\n";
1080 file_basic << std::setprecision(8) << (double)subspaceErrorLimit << "\n";
1081 file_basic << use_diag_on_error << "\n";
1082 file_basic << use_diag_on_error_guess << "\n";
1083 file_basic << create_m_files << "\n";
1084 file_basic << output_homo_and_lumo_eigenvectors << "\n";
1085 file_basic << use_prev_vector_as_initial_guess << "\n";
1086 file_basic << ignore_purification_failure << "\n";
1087 file_basic << use_rand_perturbation_for_alleigsint << "\n";
1088 file_basic << stats_prefix << "\n";
1089 file_basic << plot_puri_results << "\n";
1090 file_basic << plot_puri_results_str << "\n";
1091 file_basic << use_acceleration << "\n";
1092 file_basic << use_new_stopping_criterion << "\n";
1093 file_basic << eigenvectors_method << "\n";
1094 file_basic << eigenvectors_iterative_method << "\n";
1095 file_basic << 1 << "\n"; // old parameter number_of_eigenvalues, remains for the compatibility
1096 file_basic << std::setprecision(8) << (double)eigensolver_accuracy << "\n";
1097 file_basic << eigensolver_maxiter << "\n";
1098 file_basic << n << "\n";
1099 file_basic << noOfOccupiedOrbs << "\n";
1100 file_basic << (double)factor << "\n";
1101 file_basic << std::setprecision(8) << (double)invCholFactor_euclnorm << "\n";
1102 file_basic << maxMul << "\n";
1103 file_basic << leavesSizeMax << "\n";
1104 file_basic << blocksize << "\n";
1105
1106
1107 if ((output_homo_and_lumo_eigenvectors) && (eigenvectors_method != NA_STRING))
1108 {
1109 if (!eigVecLUMO->is_empty())
1110 {
1111 file_basic << 1 << "\n";
1112 eigVecLUMO->writeToFile();
1113
1114 name.clear();
1115 name.str("");
1116 name << filenameeigVecLUMO << "_" << checkpoint_ID.str() << ".bin";
1117 eigVecLUMO->copyToFile(name.str().c_str());
1118
1119 eigVecLUMO->readFromFile();
1120 }
1121 else
1122 {
1123 file_basic << 0 << "\n";
1124 }
1125
1126 if (!eigVecHOMO->is_empty())
1127 {
1128 file_basic << 1 << "\n";
1129 eigVecHOMO->writeToFile();
1130
1131 name.clear();
1132 name.str("");
1133 name << filenameeigVecHOMO << "_" << checkpoint_ID.str() << ".bin";
1134 eigVecHOMO->copyToFile(name.str().c_str());
1135
1136 eigVecHOMO->readFromFile();
1137 }
1138 else
1139 {
1140 file_basic << 0 << "\n";
1141 }
1142 }
1143
1144
1145 // all non-trivial
1146
1147 switch (truncationNormPurification)
1148 {
1149 case mat::frobNorm:
1150 file_basic << 1 << "\n";
1151 break;
1152
1153 case mat::mixedNorm:
1154 file_basic << 2 << "\n";
1155 break;
1156
1157 case mat::euclNorm:
1158 file_basic << 3 << "\n";
1159 break;
1160
1161 default:
1162 throw std::runtime_error("GetDensFromFock::create_checkpoint unknown truncation norm.");
1163 }
1164
1165 switch (stopCriterionNormPurification)
1166 {
1167 case mat::frobNorm:
1168 file_basic << 1 << "\n";
1169 break;
1170
1171 case mat::mixedNorm:
1172 file_basic << 2 << "\n";
1173 break;
1174
1175 case mat::euclNorm:
1176 file_basic << 3 << "\n";
1177 break;
1178
1179 default:
1180 throw std::runtime_error("GetDensFromFock::create_checkpoint unknown stopping criterion norm.");
1181 }
1182
1183 std::vector<int> blockSizesCopy;
1184 matrixSizesAndBlocks.getBlockSizeVector(blockSizesCopy);
1185 file_basic << blockSizesCopy.size() << "\n";
1186 for (unsigned int i = 0; i < blockSizesCopy.size(); ++i)
1187 {
1188 file_basic << blockSizesCopy[i] << "\n";
1189 }
1190
1191
1192 file_basic << std::setprecision(16) << (double)homoInterval_Finput.low() << "\n";
1193 file_basic << std::setprecision(16) << (double)homoInterval_Finput.upp() << "\n";
1194 file_basic << std::setprecision(16) << (double)lumoInterval_Finput.low() << "\n";
1195 file_basic << std::setprecision(16) << (double)lumoInterval_Finput.upp() << "\n";
1196
1197 file_basic << std::setprecision(16) << (double)homoInterval_F_ort_prev.low() << "\n";
1198 file_basic << std::setprecision(16) << (double)homoInterval_F_ort_prev.upp() << "\n";
1199 file_basic << std::setprecision(16) << (double)lumoInterval_F_ort_prev.low() << "\n";
1200 file_basic << std::setprecision(16) << (double)lumoInterval_F_ort_prev.upp() << "\n";
1201
1202 file_basic.close();
1203 }
1204
1205
1206
file_exist(const std::string & name)1207 inline bool file_exist(const std::string& name) {
1208 ifstream f(name.c_str());
1209 return f.good();
1210 }
1211
1212
1213 /** Function restores data from files in order to repeat recursive
1214 * expansion in a desired SCF cycle. The purpose of the function is
1215 * mainly testing.
1216 */
restore_from_checkpoint(GetDensFromFock & DensFromFock,symmMatrix & Finput,symmMatrix & F_ort_prev,generalVector * eigVecLUMO,generalVector * eigVecHOMO,std::string checkpoint_path,std::string IDstr,int SCF_step)1217 void GetDensFromFock::restore_from_checkpoint(GetDensFromFock& DensFromFock,
1218 symmMatrix& Finput,
1219 symmMatrix& F_ort_prev,
1220 generalVector *eigVecLUMO,
1221 generalVector *eigVecHOMO,
1222 std::string checkpoint_path,
1223 std::string IDstr,
1224 int SCF_step)
1225 {
1226 const char *filenameFinput = DensFromFock.filenameFinput;
1227 const char *filenameF_ort_prev = DensFromFock.filenameF_ort_prev;
1228 const char *filenameeigVecLUMO = DensFromFock.filenameeigVecLUMO;
1229 const char *filenameeigVecHOMO = DensFromFock.filenameeigVecHOMO;
1230 const char *filenameOverlap = DensFromFock.filenameOverlap;
1231 const char *filenameinvCholFactor = DensFromFock.filenameinvCholFactor;
1232 const char *file_for_basic_types = DensFromFock.file_for_basic_types;
1233
1234
1235 // read all basic data
1236 std::ostringstream name;
1237 name << checkpoint_path << "/" << file_for_basic_types << "_" << IDstr << "_" << SCF_step << ".data";
1238 ifstream file_basic(name.str().c_str(), ios::in);
1239 if (!file_basic.is_open())
1240 {
1241 throw std::runtime_error("GetDensFromFock::restore_from_checkpoint unable open file for reading.");
1242 }
1243
1244 string tmp;
1245 double tmp_double;
1246
1247 file_basic >> DensFromFock.SCF_step;
1248 file_basic >> DensFromFock.use_diagonalization;
1249 file_basic >> DensFromFock.use_purification;
1250 file_basic >> DensFromFock.store_all_eigenvalues_to_file;
1251 file_basic >> DensFromFock.try_eigv_on_next_iteration_if_fail;
1252 file_basic >> tmp_double;
1253 DensFromFock.electronicTemperature = tmp_double;
1254 file_basic >> tmp_double;
1255 DensFromFock.gap_expected_lower_bound = tmp_double;
1256 file_basic >> tmp_double;
1257 DensFromFock.eigvalueErrorLimit = tmp_double;
1258 file_basic >> tmp_double;
1259 DensFromFock.subspaceErrorLimit = tmp_double;
1260 file_basic >> DensFromFock.use_diag_on_error;
1261 file_basic >> DensFromFock.use_diag_on_error_guess;
1262 file_basic >> DensFromFock.create_m_files;
1263 file_basic >> DensFromFock.output_homo_and_lumo_eigenvectors;
1264 file_basic >> DensFromFock.use_prev_vector_as_initial_guess;
1265 file_basic >> DensFromFock.ignore_purification_failure;
1266 file_basic >> DensFromFock.use_rand_perturbation_for_alleigsint;
1267 getline(file_basic, tmp);
1268 getline(file_basic, DensFromFock.stats_prefix);
1269 //file_basic >> DensFromFock.stats_prefix;
1270 file_basic >> DensFromFock.plot_puri_results;
1271 getline(file_basic, tmp);
1272 getline(file_basic, DensFromFock.plot_puri_results_str);
1273 //file_basic >> DensFromFock.plot_puri_results_str;
1274 file_basic >> DensFromFock.use_acceleration;
1275 file_basic >> DensFromFock.use_new_stopping_criterion;
1276 getline(file_basic, tmp);
1277 getline(file_basic, DensFromFock.eigenvectors_method);
1278 //file_basic >> DensFromFock.eigenvectors_method;
1279 getline(file_basic, DensFromFock.eigenvectors_iterative_method);
1280 //file_basic >> DensFromFock.eigenvectors_iterative_method;
1281 int dummy_number_of_eigenvalues; // removed parameter, remains for the compatibility
1282 file_basic >> dummy_number_of_eigenvalues;
1283 file_basic >> tmp_double;
1284 DensFromFock.eigensolver_accuracy = tmp_double;
1285 file_basic >> DensFromFock.eigensolver_maxiter;
1286 file_basic >> DensFromFock.n;
1287 file_basic >> DensFromFock.noOfOccupiedOrbs;
1288 file_basic >> tmp_double;
1289 DensFromFock.factor = tmp_double;
1290 file_basic >> tmp_double;
1291 DensFromFock.invCholFactor_euclnorm = tmp_double;
1292 file_basic >> DensFromFock.maxMul;
1293 file_basic >> DensFromFock.leavesSizeMax;
1294 file_basic >> DensFromFock.blocksize;
1295
1296 int vector_lumo_not_null, vector_homo_not_null;
1297 file_basic >> vector_lumo_not_null;
1298 file_basic >> vector_homo_not_null;
1299
1300 int norm_trunc_ID;
1301 file_basic >> norm_trunc_ID;
1302 switch (norm_trunc_ID)
1303 {
1304 case 1:
1305 DensFromFock.truncationNormPurification = mat::frobNorm;
1306 break;
1307
1308 case 2:
1309 DensFromFock.truncationNormPurification = mat::mixedNorm;
1310 break;
1311
1312 case 3:
1313 DensFromFock.truncationNormPurification = mat::euclNorm;
1314 break;
1315
1316 default:
1317 throw std::runtime_error("GetDensFromFock::restore_from_checkpoint unknown truncation norm.");
1318 }
1319
1320
1321 int norm_st_crit_ID;
1322 file_basic >> norm_st_crit_ID;
1323 switch (norm_st_crit_ID)
1324 {
1325 case 1:
1326 DensFromFock.stopCriterionNormPurification = mat::frobNorm;
1327 break;
1328
1329 case 2:
1330 DensFromFock.stopCriterionNormPurification = mat::mixedNorm;
1331 break;
1332
1333 case 3:
1334 DensFromFock.stopCriterionNormPurification = mat::euclNorm;
1335 break;
1336
1337 default:
1338 throw std::runtime_error("GetDensFromFock::restore_from_checkpoint unknown stopping criterion norm.");
1339 }
1340
1341
1342 int blockSizes_size;
1343 file_basic >> blockSizes_size;
1344 std::vector<int> blockSizes(blockSizes_size);
1345 for (int i = 0; i < blockSizes_size; ++i)
1346 {
1347 file_basic >> blockSizes[i];
1348 }
1349
1350 DensFromFock.matrixSizesAndBlocks = mat::SizesAndBlocks(blockSizes, DensFromFock.n);
1351
1352 // eigenvalue bounds
1353 double lower, upper; // Elias note: changed from ergo_real to double here to make it work when ergo_real is __float128
1354 file_basic >> lower >> upper;
1355 DensFromFock.homoInterval_Finput = intervalType(lower, upper);
1356 file_basic >> lower >> upper;
1357 DensFromFock.lumoInterval_Finput = intervalType(lower, upper);
1358
1359 file_basic >> lower >> upper;
1360 DensFromFock.homoInterval_F_ort_prev = intervalType(lower, upper);
1361 file_basic >> lower >> upper;
1362 DensFromFock.lumoInterval_F_ort_prev = intervalType(lower, upper);
1363
1364
1365
1366 file_basic.close();
1367
1368 cout << "Finished to read text data. Starting with binary data." << endl;
1369
1370
1371 // initialize all matrices and vectors with corresponding block dimensions
1372
1373 name.clear();
1374 name.str("");
1375 name << checkpoint_path << "/" << filenameFinput << "_" << IDstr << "_" << SCF_step << ".bin";
1376
1377 if(!file_exist(name.str())) throw std::runtime_error("File " + name.str() + " does not exist!");
1378
1379 Finput.resetSizesAndBlocks(DensFromFock.matrixSizesAndBlocks, DensFromFock.matrixSizesAndBlocks); // set data structure
1380 Finput.copyFromFile(name.str().c_str());
1381
1382 name.clear();
1383 name.str("");
1384 name << checkpoint_path << "/" << filenameF_ort_prev << "_" << IDstr << "_" << SCF_step << ".bin";
1385
1386 if(!file_exist(name.str())) throw std::runtime_error("File " + name.str() + " does not exist!");
1387
1388 F_ort_prev.resetSizesAndBlocks(DensFromFock.matrixSizesAndBlocks, DensFromFock.matrixSizesAndBlocks); // set data structure
1389 F_ort_prev.copyFromFile(name.str().c_str());
1390
1391
1392 name.clear();
1393 name.str("");
1394 name << checkpoint_path << "/" << filenameOverlap << "_" << IDstr << "_" << SCF_step << ".bin";
1395
1396 if(!file_exist(name.str())) throw std::runtime_error("File " + name.str() + " does not exist!");
1397
1398 DensFromFock.overlapMatrix.resetSizesAndBlocks(DensFromFock.matrixSizesAndBlocks, DensFromFock.matrixSizesAndBlocks); // set data structure
1399 DensFromFock.overlapMatrix.copyFromFile(name.str().c_str());
1400
1401 name.clear();
1402 name.str("");
1403 name << checkpoint_path << "/" << filenameinvCholFactor << "_" << IDstr << "_" << SCF_step << ".bin";
1404
1405 if(!file_exist(name.str())) throw std::runtime_error("File " + name.str() + " does not exist!");
1406
1407 DensFromFock.invCholFactor.resetSizesAndBlocks(DensFromFock.matrixSizesAndBlocks, DensFromFock.matrixSizesAndBlocks); // set data structure
1408 DensFromFock.invCholFactor.copyFromFile(name.str().c_str());
1409
1410 // HOMO and LUMO eigenvectors
1411
1412 if (vector_lumo_not_null)
1413 {
1414 eigVecLUMO->resetSizesAndBlocks(DensFromFock.matrixSizesAndBlocks);// set data structure
1415 name.clear();
1416 name.str("");
1417 name << checkpoint_path << "/" << filenameeigVecLUMO << "_" << IDstr << "_" << SCF_step << ".bin";
1418
1419 if(!file_exist(name.str())) throw std::runtime_error("File " + name.str() + " does not exist!");
1420
1421 eigVecLUMO->copyFromFile(name.str().c_str());
1422 eigVecLUMO->readFromFile();
1423 }
1424
1425 if (vector_homo_not_null)
1426 {
1427 eigVecHOMO->resetSizesAndBlocks(DensFromFock.matrixSizesAndBlocks);// set data structure
1428 name.clear();
1429 name.str("");
1430 name << checkpoint_path << "/" << filenameeigVecHOMO << "_" << IDstr << "_" << SCF_step << ".bin";
1431
1432 if(!file_exist(name.str())) throw std::runtime_error("File " + name.str() + " does not exist!");
1433
1434 eigVecHOMO->copyFromFile(name.str().c_str());
1435 eigVecHOMO->readFromFile();
1436 }
1437
1438 }
1439