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 SCF_restricted.cc
31
32 @brief Class for self-consistent field (SCF) procedure;
33 spin-restricted case.
34
35 @author: Elias Rudberg <em>responsible</em>.
36 */
37
38 #include <sstream>
39 #include "SCF_restricted.h"
40 #include "output.h"
41 #include "scf_utils.h"
42 #include "utilities.h"
43 #include "diis_restricted.h"
44 #include "density_projection.h"
45 #include "density_description_file.h"
46 #include "matrix_utilities.h"
47 #include "machine_epsilon.h"
48 #include "units.h"
49 #include "atom_labels.h"
50 #include "integral_matrix_wrappers.h"
51 #include "dipole_moment.h"
52
53
SCF_restricted(const Molecule & molecule_,const Molecule & extraCharges_,const BasisInfoStruct & basisInfo_,const IntegralInfo & integralInfo_,const char * guessDmatFileNamePtr,const JK::Params & J_K_paramsPtr,const Dft::GridParams & gridParams_,const SCF::Options & scfopts,const SCF::MatOptions & matOpts_,ergo_real threshold_integrals_1el_input)54 SCF_restricted::SCF_restricted(const Molecule& molecule_,
55 const Molecule& extraCharges_,
56 const BasisInfoStruct& basisInfo_,
57 const IntegralInfo& integralInfo_,
58 const char *guessDmatFileNamePtr,
59 const JK::Params& J_K_paramsPtr,
60 const Dft::GridParams& gridParams_,
61 const SCF::Options& scfopts,
62 const SCF::MatOptions& matOpts_,
63 ergo_real threshold_integrals_1el_input)
64 : SCF_general(molecule_,
65 extraCharges_,
66 basisInfo_,
67 integralInfo_,
68 guessDmatFileNamePtr,
69 J_K_paramsPtr,
70 gridParams_,
71 scfopts,
72 matOpts_,
73 threshold_integrals_1el_input)
74 {
75 DIIS = new DIISManagerRestricted;
76
77
78 DensFromFock.do_restricted_calculations(); // set factor = 2
79 DensFromFock.set_no_occupied_orbs(noOfElectrons / 2);
80 }
81
82
~SCF_restricted()83 SCF_restricted::~SCF_restricted()
84 {
85 delete ((DIISManagerRestricted *)DIIS);
86 }
87
88
get_Fock_matrix(symmMatrix & FockMatrix_)89 void SCF_restricted::get_Fock_matrix(symmMatrix& FockMatrix_)
90 {
91 FockMatrix.readFromFile();
92 FockMatrix_ = FockMatrix;
93 FockMatrix.writeToFile();
94 }
95
96
get_density_matrix(symmMatrix & densityMatrix_)97 void SCF_restricted::get_density_matrix(symmMatrix& densityMatrix_)
98 {
99 densityMatrix.readFromFile();
100 densityMatrix_ = densityMatrix;
101 densityMatrix.writeToFile();
102 }
103
104
initialize_matrices()105 void SCF_restricted::initialize_matrices()
106 {
107 densityMatrix.resetSizesAndBlocks(matOpts.size_block_info,
108 matOpts.size_block_info);
109
110 densityMatrix_core.resetSizesAndBlocks(matOpts.size_block_info, matOpts.size_block_info);
111
112 twoel_matrix_core.resetSizesAndBlocks(matOpts.size_block_info, matOpts.size_block_info);
113
114 FockMatrix.resetSizesAndBlocks(matOpts.size_block_info,
115 matOpts.size_block_info);
116 Fprev.resetSizesAndBlocks(matOpts.size_block_info,
117 matOpts.size_block_info);
118 Dprev.resetSizesAndBlocks(matOpts.size_block_info,
119 matOpts.size_block_info);
120 F_ort_prev.resetSizesAndBlocks(matOpts.size_block_info,
121 matOpts.size_block_info);
122 // D_ort_prev.resetSizesAndBlocks(matOpts.size_block_info,
123 // matOpts.size_block_info);
124 bestFockMatrixSoFar.resetSizesAndBlocks(matOpts.size_block_info,
125 matOpts.size_block_info);
126 bestFockMatrixSoFar2.resetSizesAndBlocks(matOpts.size_block_info,
127 matOpts.size_block_info);
128 ErrorMatrix.resetSizesAndBlocks(matOpts.size_block_info,
129 matOpts.size_block_info);
130 }
131
132
check_params()133 void SCF_restricted::check_params()
134 {
135 if (noOfElectrons % 2 != 0)
136 {
137 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error: odd number of electrons in restricted calculation");
138 throw "error: odd number of electrons in restricted calculation";
139 }
140 }
141
142
get_starting_guess_density()143 void SCF_restricted::get_starting_guess_density()
144 {
145 // set up starting guess
146
147 int n = basisInfo.noOfBasisFuncs;
148
149 DensFromFock.set_SCF_step(SCF_step);
150
151 if (guessDmatFileName != NULL)
152 {
153 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "getting starting guess density from file '%s'", guessDmatFileName);
154 int noOfDensityMatrices = 1;
155 symmMatrix *matrixList[2];
156 matrixList[0] = &densityMatrix;
157
158 if (load_density_and_project_sparse(DensFromFock,
159 guessDmatFileName,
160 noOfDensityMatrices,
161 &integralInfo,
162 basisInfo,
163 S_symm,
164 matrixList,
165 &noOfElectrons,
166 matOpts.size_block_info,
167 matOpts.permutationHML,
168 matOpts.sparse_threshold) != 0)
169 {
170 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error in load_density_and_project_sparse");
171 throw "error in load_density_and_project_sparse";
172 }
173 }
174 else
175 {
176 if (scfopts.use_simple_starting_guess == 1)
177 {
178 if (get_simple_starting_guess_sparse(n, noOfElectrons, densityMatrix) != 0)
179 {
180 throw "error in get_simple_starting_guess_sparse";
181 }
182 densityMatrix.writeToFile();
183 }
184 else if (scfopts.use_diag_guess_from_file == 1)
185 {
186 if (get_diag_matrix_from_file(n, densityMatrix, "diagdens.txt",
187 matOpts.permutationHML) != 0)
188 {
189 throw "error in get_diag_matrix_from_file";
190 }
191 densityMatrix.writeToFile();
192 }
193 else
194 {
195 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
196 "calling get_dens_from_fock to diagonalize H_core for starting guess, n = %i, sparse_threshold = %g",
197 n, (double)matOpts.sparse_threshold);
198
199 symmMatrix F_ort_prev_dummy;
200 F_ort_prev_dummy.resetSizesAndBlocks(matOpts.size_block_info,
201 matOpts.size_block_info);
202 F_ort_prev_dummy.writeToFile();
203 densityMatrix.writeToFile();
204
205 DensFromFock.clean_eigs_intervals();
206
207 int use_diag_on_error = DensFromFock.get_use_diag_on_error();
208
209 if (DensFromFock.get_use_diag_on_error_guess() == 1)
210 {
211 DensFromFock.set_use_diag_on_error();
212 }
213
214
215 if (DensFromFock.get_dens_from_fock(H_core_Matrix,
216 densityMatrix,
217 F_ort_prev_dummy) != 0)
218 {
219 throw "SCF_restricted::get_starting_guess_density: Error in get_dens_from_fock_general";
220 }
221
222 if (use_diag_on_error != 1)
223 {
224 DensFromFock.unset_use_diag_on_error();
225 }
226 } // END ELSE use H_core
227 } // END ELSE no dmat given
228
229 densityMatrix.readFromFile();
230 output_sparsity_symm(n, densityMatrix, "starting guess density matrix");
231 densityMatrix.writeToFile();
232
233 densityMatrix_core.writeToFile(); // densityMatrix_core (if used) also needs to be written to file
234 }
235
236
add_random_disturbance_to_starting_guess()237 void SCF_restricted::add_random_disturbance_to_starting_guess()
238 {
239 if (scfopts.sg_disturb_specific_elements > SCF::DISTURB_ELEMENT_MAX_COUNT)
240 {
241 throw "Error in SCF_restricted::add_random_disturbance_to_starting_guess: (scfopts.sg_disturb_specific_elements > SCF::DISTURB_ELEMENT_MAX_COUNT)";
242 }
243 int n = basisInfo.noOfBasisFuncs;
244 densityMatrix.readFromFile();
245 add_disturbance_to_matrix(n,
246 densityMatrix,
247 scfopts.starting_guess_disturbance,
248 scfopts.sg_disturb_specific_elements,
249 scfopts.disturbedElementIndexVector,
250 matOpts.permutationHML);
251 densityMatrix.writeToFile();
252 }
253
254
initialize_homo_lumo_limits()255 void SCF_restricted::initialize_homo_lumo_limits()
256 {
257 intervalType hugeInterval(-1e22, 1e22);
258
259 homoInterval_F_ort_prev = hugeInterval;
260 lumoInterval_F_ort_prev = hugeInterval;
261 homoInterval_Fprev = hugeInterval;
262 lumoInterval_Fprev = hugeInterval;
263 }
264
265
write_matrices_to_file()266 void SCF_restricted::write_matrices_to_file()
267 {
268 FockMatrix.writeToFile();
269 Fprev.writeToFile();
270 Dprev.writeToFile();
271 bestFockMatrixSoFar.writeToFile();
272 bestFockMatrixSoFar2.writeToFile();
273 F_ort_prev.writeToFile();
274 }
275
276
output_diff_norm_values(symmMatrix const & F1,symmMatrix const & F2,ergo_real acc,const char * name)277 static void output_diff_norm_values(symmMatrix const& F1,
278 symmMatrix const& F2,
279 ergo_real acc,
280 const char *name)
281 {
282 ergo_real E_norm_frob = symmMatrix::frob_diff(F1, F2);
283
284 Util::TimeMeter timeMeterMixedDiff;
285 ergo_real E_norm_mixed = symmMatrix::mixed_diff(F1, F2, acc);
286 timeMeterMixedDiff.print(LOG_AREA_DENSFROMF, "symmMatrix::mixed_diff");
287 Util::TimeMeter timeMeterEuclDiff;
288 ergo_real E_norm_eucl = symmMatrix::eucl_diff(F1, F2, acc);
289 timeMeterEuclDiff.print(LOG_AREA_DENSFROMF, "symmMatrix::eucl_diff ");
290 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Nor of error matrix for '%s':", name);
291 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "frob : %22.15f", (double)E_norm_frob);
292 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "mixed: %22.15f", (double)E_norm_mixed);
293 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "eucl: %22.15f", (double)E_norm_eucl);
294 }
295
296
get_eucl_diff_with_adapted_accuracy(int n,const symmMatrix & F_w,const symmMatrix & F_ort_prev_w,ergo_real acc)297 static ergo_real get_eucl_diff_with_adapted_accuracy(int n,
298 const symmMatrix & F_w,
299 const symmMatrix & F_ort_prev_w,
300 ergo_real acc) {
301 // The symmMatrix::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.
302 int maxIterForEuclDiff = std::max(n / 10, 500);
303 ergo_real maxEigValMovement_eucl = -1; // Value will be set in try/catch code below.
304 try {
305 Util::TimeMeter timeMeterEuclDiff;
306 maxEigValMovement_eucl = symmMatrix::eucl_diff(F_w, F_ort_prev_w, acc, maxIterForEuclDiff) + acc;
307 timeMeterEuclDiff.print(LOG_AREA_SCF, "symmMatrix::eucl_diff for maxEigValMovement_eucl ");
308 }
309 catch(...) {
310 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "symmMatrix::eucl_diff() for maxEigValMovement_eucl failed for maxIterForEuclDiff=%d. Calling eucl_diff() again with lower accuracy requirement sqrt(acc).",
311 maxIterForEuclDiff);
312 ergo_real acc2 = template_blas_sqrt(acc);
313 Util::TimeMeter timeMeterEuclDiff;
314 maxEigValMovement_eucl = symmMatrix::eucl_diff(F_w, F_ort_prev_w, acc2) + acc2;
315 timeMeterEuclDiff.print(LOG_AREA_SCF, "symmMatrix::eucl_diff for maxEigValMovement_eucl ");
316 }
317 return maxEigValMovement_eucl;
318 }
319
320
get_2e_part_and_energy()321 void SCF_restricted::get_2e_part_and_energy()
322 {
323 densityMatrix.readFromFile();
324 densityMatrix_core.readFromFile();
325
326 bool scan_do_invcholfactor_transf = (bool)scfopts.scan_do_invcholfactor_transf;
327
328 if (scfopts.do_acc_scan_J)
329 {
330 do_acc_scan_J(densityMatrix,
331 integralInfo,
332 basisInfo,
333 invCholFactor,
334 scan_do_invcholfactor_transf,
335 J_K_params,
336 matOpts.size_block_info,
337 matOpts.permutationHML,
338 scfopts.scan_no_of_steps,
339 scfopts.scan_start_thresh,
340 scfopts.scan_step_factor);
341 }
342
343 if (scfopts.do_acc_scan_K)
344 {
345 do_acc_scan_K(densityMatrix,
346 integralInfo,
347 basisInfo,
348 invCholFactor,
349 scan_do_invcholfactor_transf,
350 CAM_params,
351 J_K_params,
352 matOpts.size_block_info,
353 matOpts.permutationHML,
354 matOpts.inversePermutationHML,
355 scfopts.scan_no_of_steps,
356 scfopts.scan_start_thresh,
357 scfopts.scan_step_factor);
358 }
359
360 if (scfopts.do_acc_scan_Vxc)
361 {
362 do_acc_scan_Vxc(densityMatrix,
363 integralInfo,
364 basisInfo,
365 molecule,
366 gridParams,
367 noOfElectrons,
368 invCholFactor,
369 scan_do_invcholfactor_transf,
370 matOpts.size_block_info,
371 matOpts.permutationHML,
372 matOpts.inversePermutationHML,
373 scfopts.scan_no_of_steps,
374 scfopts.scan_start_thresh,
375 scfopts.scan_step_factor);
376 }
377
378 symmMatrix G;
379 G.resetSizesAndBlocks(matOpts.size_block_info,
380 matOpts.size_block_info);
381
382 if (get_2e_matrix_and_energy_sparse(basisInfo,
383 molecule,
384 integralInfo,
385 G,
386 densityMatrix,
387 J_K_params,
388 CAM_params,
389 gridParams,
390 scfopts.use_dft,
391 &energy_2el, noOfElectrons,
392 matOpts.size_block_info,
393 matOpts.permutationHML,
394 matOpts.inversePermutationHML,
395 0,
396 J_matrix,
397 K_matrix,
398 Fxc_matrix,
399 *curr_cycle_stats) != 0)
400 {
401 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error in get_2e_matrix_and_energy_sparse");
402 throw "error in get_2e_matrix_and_energy_sparse";
403 }
404
405 if (scfopts.compute_core_density == 1)
406 {
407 // Get energy_2el_valence and H_matrix_core, but only if densityMatrix_core is nonzero (it is zero the first time, then we just skip).
408 if (densityMatrix_core.frob() != 0)
409 {
410 // Get energy_2el_valence
411 {
412 symmMatrix densityMatrix_valence(densityMatrix);
413 densityMatrix_valence += (ergo_real) - 1 * densityMatrix_core;
414 symmMatrix G_valence;
415 G_valence.resetSizesAndBlocks(matOpts.size_block_info, matOpts.size_block_info);
416 symmMatrix J_matrix_valence;
417 symmMatrix K_matrix_valence;
418 symmMatrix Fxc_matrix_valence;
419 SCF_statistics curr_cycle_stats_dummy;
420 int no_of_valence_electrons = noOfElectrons - scfopts.no_of_core_electrons;
421 if (get_2e_matrix_and_energy_sparse(basisInfo,
422 molecule,
423 integralInfo,
424 G_valence,
425 densityMatrix_valence,
426 J_K_params,
427 CAM_params,
428 gridParams,
429 scfopts.use_dft,
430 &energy_2el_valence, no_of_valence_electrons,
431 matOpts.size_block_info,
432 matOpts.permutationHML,
433 matOpts.inversePermutationHML,
434 0,
435 J_matrix_valence,
436 K_matrix_valence,
437 Fxc_matrix_valence,
438 curr_cycle_stats_dummy) != 0)
439 {
440 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "Error in get_2e_matrix_and_energy_sparse for densityMatrix_valence.");
441 throw "Error in get_2e_matrix_and_energy_sparse for densityMatrix_valence.";
442 }
443 }
444 // Compute J_matrix_core
445 {
446 symmMatrix J_matrix_core;
447 symmMatrix K_matrix_core;
448 symmMatrix Fxc_matrix_core;
449 SCF_statistics curr_cycle_stats_dummy;
450 energy_2el_core = 0;
451 if (get_2e_matrix_and_energy_sparse(basisInfo,
452 molecule,
453 integralInfo,
454 twoel_matrix_core,
455 densityMatrix_core,
456 J_K_params,
457 CAM_params,
458 gridParams,
459 scfopts.use_dft,
460 &energy_2el_core, scfopts.no_of_core_electrons,
461 matOpts.size_block_info,
462 matOpts.permutationHML,
463 matOpts.inversePermutationHML,
464 0,
465 J_matrix_core,
466 K_matrix_core,
467 Fxc_matrix_core,
468 curr_cycle_stats_dummy) != 0)
469 {
470 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "Error in get_2e_matrix_and_energy_sparse for densityMatrix_core.");
471 throw "Error in get_2e_matrix_and_energy_sparse for densityMatrix_core.";
472 }
473 }
474 }
475 }
476
477 densityMatrix.writeToFile();
478 densityMatrix_core.writeToFile();
479
480 // Check that G matrix is free from "inf", "nan" etc.
481 if (check_if_matrix_contains_strange_elements(G, matOpts.inversePermutationHML))
482 {
483 throw std::runtime_error("error in SCF_restricted::get_2e_part_and_energy(): G matrix contains inf or nan.");
484 }
485
486 // calculate Fock matrix F = H_core + G
487 H_core_Matrix.readFromFile();
488 FockMatrix.readFromFile();
489 FockMatrix = H_core_Matrix + G;
490 H_core_Matrix.writeToFile();
491 G.clear();
492
493 // Save a copy of FockMatrix before truncation if verification requested.
494 symmMatrix FockMatrixBeforeTruncation;
495 if (scfopts.do_f_thresh_verification == 1)
496 {
497 FockMatrixBeforeTruncation = FockMatrix;
498 }
499
500 // Do truncation of FockMatrix taking into account gap of Fprev and norm of Z.
501 ergo_real gapOfFprevMin = lumoInterval_Fprev.low() - homoInterval_Fprev.upp();
502 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "About to truncate FockMatrix, gap of Fprev >= %22.11f", (double)gapOfFprevMin);
503 {
504 // Compare FockMatrix to Fprev to check how far eigenvalues can have moved.
505 Fprev.readFromFile();
506 ergo_real maxEigValMovement_frob = symmMatrix::frob_diff(FockMatrix, Fprev);
507 ergo_real acc = template_blas_sqrt(get_machine_epsilon());
508 Util::TimeMeter timeMeterMixedDiff;
509 ergo_real maxEigValMovement_mixed = symmMatrix::mixed_diff(FockMatrix, Fprev, acc) + acc;
510 timeMeterMixedDiff.print(LOG_AREA_SCF, "symmMatrix::mixed_diff for F vs Fprev maxEigValMovement_mixed");
511 int n = basisInfo.noOfBasisFuncs;
512 ergo_real maxEigValMovement_eucl = get_eucl_diff_with_adapted_accuracy(n, FockMatrix, Fprev, acc);
513 Fprev.writeToFile();
514 // Increase HOMO/LUMO intervals so that they for sure contain the HOMO and LUMO eigenvalues of F_ort
515 intervalType homoInterval = homoInterval_Fprev;
516 intervalType lumoInterval = lumoInterval_Fprev;
517 homoInterval.increase(maxEigValMovement_eucl);
518 lumoInterval.increase(maxEigValMovement_eucl);
519 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "maxEigValMovement_frob = %22.11f", (double)maxEigValMovement_frob);
520 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "maxEigValMovement_mixed = %22.11f", (double)maxEigValMovement_mixed);
521 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "maxEigValMovement_eucl = %22.11f", (double)maxEigValMovement_eucl);
522 // Now we have homoInterval and lumoInterval valid for FockMatrix.
523 ergo_real gapMin = lumoInterval.low() - homoInterval.upp();
524 ergo_real gapMax = lumoInterval.upp() - homoInterval.low();
525 ergo_real threshold_1;
526 // Choose subspace error as for purification.
527 ergo_real subspaceThr_1 = 0.1 * scfopts.purification_subspace_err_limit;
528 // We consider the gap to be accurately known if the uncertainty is at most 10 %
529 if ((gapMin > 0) && ((gapMax - gapMin) / gapMin < 0.1))
530 {
531 // Gap is accurately known: we use gapMin
532 threshold_1 = subspaceThr_1 * gapMin / (1 + subspaceThr_1);
533 }
534 else
535 {
536 // Gap is not accurately known. To avoid choosing a very tight
537 // threshold value due to a small lower bound for the gap, we
538 // use the largest of 'gap_expected_lower_bound' and calculated
539 // 'gapMin':
540 threshold_1 = gapMin > scfopts.gap_expected_lower_bound ?
541 subspaceThr_1 * gapMin / (1 + subspaceThr_1) :
542 subspaceThr_1 * scfopts.gap_expected_lower_bound / (1 + subspaceThr_1);
543 }
544
545 /* Truncate matrix taking into account that we are in
546 * 'non-orthogonal basis', passing invCholFactor to thresh */
547 invCholFactor.readFromFile();
548 Util::TimeMeter timeMeterEuclThresh;
549 double nnzF_S_before_trunc_pc = (double)FockMatrix.nnz() * 100 / ((double)n * n);
550 ergo_real truncError_1 = FockMatrix.eucl_thresh(threshold_1, &invCholFactor);
551 double nnzF_S_after_trunc_pc = (double)FockMatrix.nnz() * 100 / ((double)n * n);
552 invCholFactor.writeToFile();
553 timeMeterEuclThresh.print(LOG_AREA_SCF, "FockMatrix.eucl_thresh() (with Z)");
554 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
555 "Truncated FockMatrix (eucl with Z), selected threshold = %10.6g, returned error = %10.6g, nnz before = %3.4f %%, nnz after = %3.4f %%",
556 (double)threshold_1, (double)truncError_1, nnzF_S_before_trunc_pc, nnzF_S_after_trunc_pc);
557 this->curr_cycle_stats->add_value("investigation_nnz_percentage_F_S", nnzF_S_after_trunc_pc);
558 }
559 // FockMatrix.frob_thresh(matOpts.sparse_threshold);
560
561 if (scfopts.do_f_thresh_verification == 1)
562 {
563 Util::TimeMeter timeMeterThreshVerification;
564 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "do_f_thresh_verification requested, computing error matrices...");
565 symmMatrix F1(FockMatrixBeforeTruncation);
566 symmMatrix F2(FockMatrix);
567 ergo_real acc = template_blas_sqrt(get_machine_epsilon());
568 // First get norm of error matrix in non-orthogonal basis.
569 output_diff_norm_values(F1, F2, acc, "F in in non-orthogonal basis");
570 // Now get norm of error matrix in non-orthogonal basis.
571 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "do_f_thresh_verification requested, doing Z multiplications...");
572 invCholFactor.readFromFile();
573 Util::TimeMeter timeMeterZFZ1;
574 F1 = transpose(invCholFactor) * F1 * invCholFactor;
575 timeMeterZFZ1.print(LOG_AREA_DENSFROMF, "transpose(invCholFactor) * F1 * invCholFactor");
576 Util::TimeMeter timeMeterZFZ2;
577 F2 = transpose(invCholFactor) * F2 * invCholFactor;
578 timeMeterZFZ2.print(LOG_AREA_DENSFROMF, "transpose(invCholFactor) * F2 * invCholFactor");
579 invCholFactor.writeToFile();
580 // Now we have the matrix in orthogonal basis, before and after truncation.
581 output_diff_norm_values(F1, F2, acc, "F in in orthogonal basis");
582 timeMeterThreshVerification.print(LOG_AREA_SCF, "do_f_thresh_verification stuff");
583 }
584
585 FockMatrix.writeToFile();
586 }
587
588
output_sparsity_S_F_D(SCF_statistics & stats)589 void SCF_restricted::output_sparsity_S_F_D(SCF_statistics& stats)
590 {
591 int n = basisInfo.noOfBasisFuncs;
592
593 S_symm.readFromFile();
594 output_sparsity_symm(n, S_symm, "S");
595 stats.add_value("nnz_S", S_symm.nnz());
596 S_symm.writeToFile();
597 FockMatrix.readFromFile();
598 output_sparsity_symm(n, FockMatrix, "F");
599 stats.add_value("nnz_F", FockMatrix.nnz());
600 FockMatrix.writeToFile();
601 densityMatrix.readFromFile();
602 output_sparsity_symm(n, densityMatrix, "D");
603 stats.add_value("nnz_D", densityMatrix.nnz());
604 densityMatrix.writeToFile();
605 }
606
607
calculate_energy()608 void SCF_restricted::calculate_energy()
609 {
610 // calculate energy
611 H_core_Matrix.readFromFile();
612 densityMatrix.readFromFile();
613 energy = symmMatrix::trace_ab(densityMatrix, H_core_Matrix) + energy_2el;
614 if (scfopts.compute_core_density == 1)
615 {
616 densityMatrix_core.readFromFile();
617 symmMatrix densityMatrix_valence(densityMatrix);
618 densityMatrix_valence += (ergo_real) - 1 * densityMatrix_core;
619 energy_of_valence = symmMatrix::trace_ab(densityMatrix_valence, H_core_Matrix) + symmMatrix::trace_ab(densityMatrix_valence, twoel_matrix_core) + energy_2el_valence;
620 energy_reference = symmMatrix::trace_ab(densityMatrix_core, H_core_Matrix) + energy_2el_core + nuclearEnergy;
621 densityMatrix_core.writeToFile();
622 }
623 densityMatrix.writeToFile();
624 H_core_Matrix.writeToFile();
625 energy += nuclearEnergy;
626 }
627
628
get_FDSminusSDF()629 void SCF_restricted::get_FDSminusSDF()
630 {
631 int n = basisInfo.noOfBasisFuncs;
632
633 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "calling compute_FDSminusSDF_sparse, n = %i", n);
634 densityMatrix.readFromFile();
635 S_symm.readFromFile();
636 FockMatrix.readFromFile();
637 compute_FDSminusSDF_sparse(n, FockMatrix, densityMatrix, S_symm,
638 ErrorMatrix, matOpts.sparse_threshold);
639 S_symm.writeToFile();
640 FockMatrix.writeToFile();
641 densityMatrix.writeToFile();
642 // write to file and read back again to reduce memory fragmentation.
643 ErrorMatrix.writeToFile();
644 ErrorMatrix.readFromFile();
645 output_sparsity(n, ErrorMatrix, "FDS-SDF");
646 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_restricted::get_FDSminusSDF finished.");
647 }
648
649
get_error_measure()650 void SCF_restricted::get_error_measure()
651 {
652 ergo_real error_maxabs = compute_maxabs_sparse(ErrorMatrix);
653 ergo_real error_frob = ErrorMatrix.frob();
654
655 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "maxabs FDS-SDF is %8.3g", (double)error_maxabs);
656 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "frob FDS-SDF is %8.3g", (double)error_frob);
657 errorMeasure = error_maxabs;
658 }
659
660
add_to_DIIS_list()661 void SCF_restricted::add_to_DIIS_list()
662 {
663 FockMatrix.readFromFile();
664 if (((DIISManagerRestricted *)DIIS)->AddIterationToList(FockMatrix, ErrorMatrix) != 0)
665 {
666 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error in DIIS AddIterationToList");
667 throw "error in DIIS AddIterationToList";
668 }
669 FockMatrix.writeToFile();
670 }
671
672
update_best_fock_so_far()673 void SCF_restricted::update_best_fock_so_far()
674 {
675 Fprev.readFromFile();
676 bestFockMatrixSoFar.readFromFile();
677 bestFockMatrixSoFar = Fprev;
678 Fprev.writeToFile();
679 bestFockMatrixSoFar.writeToFile();
680 FockMatrix.readFromFile();
681 bestFockMatrixSoFar2.readFromFile();
682 bestFockMatrixSoFar2 = FockMatrix;
683 FockMatrix.writeToFile();
684 bestFockMatrixSoFar2.writeToFile();
685 }
686
687
combine_old_fock_matrices(ergo_real stepLength)688 void SCF_restricted::combine_old_fock_matrices(ergo_real stepLength)
689 {
690 bestFockMatrixSoFar.readFromFile();
691 bestFockMatrixSoFar2.readFromFile();
692 FockMatrix.readFromFile();
693 FockMatrix = 0;
694 FockMatrix += stepLength * bestFockMatrixSoFar2;
695 FockMatrix += (1 - stepLength) * bestFockMatrixSoFar;
696 FockMatrix.writeToFile();
697 bestFockMatrixSoFar.writeToFile();
698 bestFockMatrixSoFar2.writeToFile();
699 }
700
701
use_diis_to_get_new_fock_matrix()702 void SCF_restricted::use_diis_to_get_new_fock_matrix()
703 {
704 symmMatrix newFsymm;
705
706 if (((DIISManagerRestricted *)DIIS)->GetCombinedFockMatrix(newFsymm) != 0)
707 {
708 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error in DIIS.GetCombinedFockMatrix");
709 throw "error in DIIS.GetCombinedFockMatrix";
710 }
711 FockMatrix.readFromFile();
712 FockMatrix = newFsymm;
713 FockMatrix.writeToFile();
714 }
715
716
clear_diis_list()717 void SCF_restricted::clear_diis_list()
718 {
719 ((DIISManagerRestricted *)DIIS)->ClearList();
720 }
721
722
clear_error_matrices()723 void SCF_restricted::clear_error_matrices()
724 {
725 ErrorMatrix.clear();
726 }
727
728
save_current_fock_as_fprev()729 void SCF_restricted::save_current_fock_as_fprev()
730 {
731 // save current Fock matrix as Fprev
732 FockMatrix.readFromFile();
733 Fprev.readFromFile();
734 Fprev = FockMatrix;
735 FockMatrix.writeToFile();
736 Fprev.writeToFile();
737 }
738
739
get_new_density_matrix()740 void SCF_restricted::get_new_density_matrix()
741 {
742 int n = basisInfo.noOfBasisFuncs;
743
744 DensFromFock.set_SCF_step(SCF_step);
745
746 // As input to the density matrix construction routine, the default
747 // is to use FockMatrix. However, if shift_using_prev_density_matrix
748 // we use a modified matrix instead.
749 symmMatrix *F_effective = &FockMatrix;
750 symmMatrix F_modified;
751 if (scfopts.shift_using_prev_density_matrix != 0)
752 {
753 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Using shift_using_prev_density_matrix, shifting by %9.5f a.u. = %9.5f eV",
754 (double)scfopts.shift_using_prev_density_matrix, (double)scfopts.shift_using_prev_density_matrix / UNIT_one_eV);
755 F_modified = FockMatrix;
756 F_modified.readFromFile();
757 // Get matrix SDS
758 symmMatrix SDS(densityMatrix);
759 SDS.readFromFile();
760 transform_with_S(SDS);
761 // Use factor 0.5 since this is restricted case.
762 F_modified += (ergo_real) - 0.5 * scfopts.shift_using_prev_density_matrix * SDS;
763 F_modified.writeToFile();
764 F_effective = &F_modified;
765 }
766
767 DensFromFock.set_eigs_F_ort_prev(homoInterval_F_ort_prev, lumoInterval_F_ort_prev);
768 DensFromFock.clean_puri_stats();
769
770 DensFromFock.set_generate_figures();
771
772 // if eigenvectors are needed, set params
773 int use_init_guess = 0;
774 if ((scfopts.use_prev_vector_as_initial_guess == 1) &&
775 (scfopts.min_number_of_iterations-1 <= SCF_step) &&
776 (SCF_step > 1) &&
777 !eigVecUNOCC.empty() && !eigVecOCC.empty()) // ensure that we computed vectors in previous cycle
778 {
779 use_init_guess = 1;
780 }
781
782 // if we use new purification
783 if ((SCF_step > 1) &&
784 (scfopts.min_number_of_iterations-1 <= SCF_step) &&
785 (DensFromFock.get_use_purification() == 1) &&
786 (DensFromFock.get_output_homo_and_lumo_eigenvectors() == 1))
787 {
788 DensFromFock.compute_eigenvectors(scfopts.eigenvectors_method,
789 scfopts.eigenvectors_iterative_method,
790 scfopts.eigensolver_accuracy,
791 scfopts.eigensolver_maxiter,
792 use_init_guess,
793 scfopts.try_eigv_on_next_iteration_if_fail);
794
795 DensFromFock.compute_eigenvectors_extra(scfopts.puri_compute_eigv_in_each_iteration, scfopts.run_shift_and_square_method_on_F);
796 }
797
798 if (scfopts.create_checkpoints == 1)
799 {
800 Util::TimeMeter timeMeter;
801 generalVector * eigVecLUMO = NULL;
802 generalVector * eigVecHOMO = NULL;
803 if(eigVecOCC.size() >= 1) eigVecHOMO = &eigVecOCC[0];
804 if(eigVecUNOCC.size() >= 1) eigVecLUMO = &eigVecUNOCC[0];
805 DensFromFock.create_checkpoint(*F_effective, // should normally be same as Fprev now
806 F_ort_prev,
807 eigVecLUMO,
808 eigVecHOMO,
809 scfopts.checkpoint_IDstr);
810 timeMeter.print(LOG_AREA_SCF, "in SCF_restricted DensFromFock::create_checkpoint took");
811 }
812
813 if (DensFromFock.get_dens_from_fock(*F_effective, // should normally be same as Fprev now
814 densityMatrix,
815 F_ort_prev) != 0)
816 {
817 throw "SCF_restricted::get_new_density_matrix: Error in get_dens_from_fock";
818 }
819
820 // return empty vectors in case nothing is computed
821 DensFromFock.get_computed_eigenpairs(eigVecUNOCC, eigVecOCC, eigValUNOCC, eigValOCC);
822
823 DensFromFock.get_eigs_F_ort_prev(homoInterval_F_ort_prev, lumoInterval_F_ort_prev);
824 DensFromFock.get_eigs_Fprev(homoInterval_Fprev, lumoInterval_Fprev);
825
826 std::map<std::string, double> puri_stats;
827 DensFromFock.get_puri_stats(puri_stats);
828 this->curr_cycle_stats->add_values(puri_stats);
829
830 DensFromFock.unset_generate_figures();
831 electronicEntropyTerm = DensFromFock.get_result_entropy_term();
832
833 if (scfopts.compute_core_density == 1)
834 {
835 DensFromFock.clean_eigs_intervals();
836 if (DensFromFock.get_dens_from_fock(*F_effective,
837 densityMatrix_core,
838 F_ort_prev) != 0)
839 {
840 throw "SCF_restricted::get_new_density_matrix: Error in get_dens_from_fock for core density matrix.";
841 }
842 }
843
844 // Report sparsity of D and trace(DS)
845 S_symm.readFromFile();
846 densityMatrix.readFromFile();
847 output_sparsity_symm(n, densityMatrix, "new density matrix");
848 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Tr( D * S ) = %22.11f", (double)symmMatrix::trace_ab(densityMatrix, S_symm));
849 densityMatrix.writeToFile();
850 S_symm.writeToFile();
851 }
852
853
write_matrix_to_file(symmMatrix & M,const std::vector<int> & inversePermutationHML,const BasisInfoStruct & basisInfo,const char * fileName)854 static int write_matrix_to_file(symmMatrix& M, const std::vector<int>& inversePermutationHML, const BasisInfoStruct& basisInfo, const char *fileName)
855 {
856 M.readFromFile();
857 matrix_description_struct matrixList[2];
858 int nvalues = M.nvalues();
859 std::vector<int> rowind;
860 rowind.reserve(nvalues);
861 std::vector<int> colind;
862 colind.reserve(nvalues);
863 std::vector<ergo_real> values;
864 values.reserve(nvalues);
865 M.get_all_values(rowind, colind, values, inversePermutationHML, inversePermutationHML);
866 M.writeToFile();
867 matrixList[0].nvalues = nvalues;
868 matrixList[0].rowind = &rowind[0];
869 matrixList[0].colind = &colind[0];
870 matrixList[0].values = &values[0];
871 if (ddf_writeShellListAndDensityMatricesToFile_sparse(basisInfo, 1, matrixList, fileName) != 0)
872 {
873 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "Error in ddf_writeShellListAndDensityMatricesToFile_sparse.");
874 return -1;
875 }
876 return 0;
877 }
878
879
write_density_to_file()880 void SCF_restricted::write_density_to_file()
881 {
882 if (write_matrix_to_file(densityMatrix, matOpts.inversePermutationHML, basisInfo, "density.bin") != 0)
883 {
884 throw "error in SCF_restricted::write_density_to_file(): write_matrix_to_file failed for densityMatrix.";
885 }
886 if (scfopts.compute_core_density == 1)
887 {
888 if (write_matrix_to_file(densityMatrix_core, matOpts.inversePermutationHML, basisInfo, "density_core.bin") != 0)
889 {
890 throw "error in SCF_restricted::write_density_to_file(): write_matrix_to_file failed for densityMatrix_core.";
891 }
892 }
893 }
894
895
save_final_potential()896 void SCF_restricted::save_final_potential()
897 {
898 if (save_symmetric_matrix(FockMatrix, basisInfo, "potential.bin",
899 matOpts.inversePermutationHML) != 0)
900 {
901 do_output(LOG_CAT_ERROR, LOG_AREA_SCF,
902 "error in ddf_writeShellListAndDensityMatricesToFile");
903 throw "error in ddf_writeShellListAndDensityMatricesToFile";
904 }
905 }
906
907
save_full_matrices_for_matlab()908 void SCF_restricted::save_full_matrices_for_matlab()
909 {
910 int n = basisInfo.noOfBasisFuncs;
911
912 FockMatrix.readFromFile();
913 write_full_matrix(n, FockMatrix, "matrix_F",
914 matOpts.inversePermutationHML);
915 FockMatrix.writeToFile();
916
917 S_symm.readFromFile();
918 write_full_matrix(n, S_symm, "matrix_S",
919 matOpts.inversePermutationHML);
920 S_symm.writeToFile();
921
922 densityMatrix.readFromFile();
923 write_full_matrix(n, densityMatrix, "matrix_D",
924 matOpts.inversePermutationHML);
925 densityMatrix.writeToFile();
926 }
927
928
output_expected_values_pos_operator()929 void SCF_restricted::output_expected_values_pos_operator()
930 {
931 if(eigVecOCC.size() == 1)
932 get_expected_values_pos_operator(eigVecOCC[0], "HOMO");
933 if(eigVecUNOCC.size() == 1)
934 get_expected_values_pos_operator(eigVecUNOCC[0], "LUMO");
935
936 if(eigVecOCC.size() > 1)
937 for (size_t i = 1; i < eigVecOCC.size(); i++) {
938 std::ostringstream name;
939 name << "OCCUPIED " << i;
940 get_expected_values_pos_operator(eigVecOCC[i], name.str().c_str());
941 }
942 if(eigVecUNOCC.size() > 1)
943 for (size_t i = 1; i < eigVecUNOCC.size(); i++) {
944 std::ostringstream name;
945 name << "UNOCCUPIED " << i;
946 get_expected_values_pos_operator(eigVecUNOCC[i], name.str().c_str());
947 }
948 }
949
950
951 // expected value of a measurement of the position of the particle
get_expected_values_pos_operator(generalVector & eigVec,const char * vector_name)952 void SCF_restricted::get_expected_values_pos_operator(generalVector& eigVec, const char *vector_name)
953 {
954 Util::TimeMeter timeMeter;
955
956 int n = basisInfo.noOfBasisFuncs;
957
958 if (eigVec.is_empty())
959 {
960 do_output(LOG_CAT_WARNING, LOG_AREA_SCF, "Failed to output expected value of a position operator for the %s eigenvector.", vector_name);
961 return;
962 }
963 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Computing expected value of a position operator for the %s eigenvector.", vector_name);
964
965 std::vector<ergo_real> vec(n);
966 eigVec.fullvector(vec);
967
968 // get density matrix corresponding to the eigenvector
969 std::vector<int> rows(n), cols(n);
970 std::vector<ergo_real> vals(n);
971 ergo_real tmp;
972 size_t count = 0;
973 for (int i = 0; i < n; ++i)
974 for (int j = i; j < n; ++j)
975 {
976 tmp = vec[i] * vec[j];
977 if(template_blas_fabs(tmp) < 1e-5)
978 continue;
979 rows[count] = i;
980 cols[count] = j;
981 vals[count] = tmp;
982 count++;
983 if(count % n == 0)
984 {
985 rows.resize(count+n);
986 cols.resize(count+n);
987 vals.resize(count+n);
988 }
989 }
990 rows.resize(count);
991 cols.resize(count);
992 vals.resize(count);
993
994 symmMatrix densityMatrix;
995 densityMatrix.resetSizesAndBlocks(matOpts.size_block_info, matOpts.size_block_info);
996 densityMatrix.assign_from_sparse(rows, cols, vals);
997
998 std::vector<ergo_real> mean;
999 std::vector<ergo_real> std;
1000 get_exp_value_pos_operator(basisInfo,
1001 molecule,
1002 densityMatrix,
1003 matOpts.size_block_info,
1004 matOpts.permutationHML,
1005 mean,
1006 std);
1007
1008
1009 if ((mean.size() != 3) || (std.size() != 3))
1010 {
1011 throw "Error in output_expected_values_pos_operator: wrong size of a vector.";
1012 }
1013
1014 ergo_real conv_const = UNIT_one_Angstrom; // convert a.u. to Angstrom
1015 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Expected value of position operator: \n (%.12lf, %.12lf, %.12lf) a.u. = (%.12lf, %.12lf, %.12lf) A",
1016 mean[0], mean[1], mean[2],
1017 mean[0] / conv_const, mean[1] / conv_const, mean[2] / conv_const);
1018 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Standard deviation of position operator: \n (%.12lf, %.12lf, %.12lf) a.u. = (%.12lf, %.12lf, %.12lf) A",
1019 std[0], std[1], std[2],
1020 std[0] / conv_const, std[1] / conv_const, std[2] / conv_const);
1021
1022 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_restricted::output_expected_values_pos_operator finished OK.");
1023 timeMeter.print(LOG_AREA_SCF, "SCF_restricted::output_expected_values_pos_operator");
1024 }
1025
1026
output_density_images()1027 void SCF_restricted::output_density_images()
1028 {
1029 Util::TimeMeter timeMeter;
1030
1031 int n = basisInfo.noOfBasisFuncs;
1032
1033 ergo_real *densityMatrixFull_tot = new ergo_real[n * n];
1034 ergo_real *densityMatrixFull_spin = new ergo_real[n * n];
1035
1036 // Get full matrix version of density matrix, and empty spin density matrix.
1037 {
1038 std::vector<ergo_real> densityMatrixFull(n *n);
1039
1040 densityMatrix.readFromFile();
1041 densityMatrix.fullMatrix(densityMatrixFull,
1042 matOpts.inversePermutationHML,
1043 matOpts.inversePermutationHML);
1044 densityMatrix.writeToFile();
1045
1046 for (int i = 0; i < n * n; i++)
1047 {
1048 densityMatrixFull_tot [i] = densityMatrixFull[i];
1049 densityMatrixFull_spin[i] = 0;
1050 }
1051 }
1052
1053 do_density_images(basisInfo,
1054 molecule,
1055 densityMatrixFull_tot,
1056 densityMatrixFull_spin,
1057 scfopts.output_density_images_boxwidth);
1058
1059 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_restricted::output_density_images finished OK.");
1060 timeMeter.print(LOG_AREA_SCF, "SCF_restricted::output_density_images");
1061 }
1062
1063
output_density_images_orbital(generalVector & eigVec,const std::string & filename_id)1064 void SCF_restricted::output_density_images_orbital(generalVector& eigVec, const std::string& filename_id)
1065 {
1066 Util::TimeMeter timeMeter;
1067
1068 int n = basisInfo.noOfBasisFuncs;
1069
1070 ergo_real *densityMatrixFull_tot = new ergo_real[n * n];
1071 ergo_real *densityMatrixFull_spin = new ergo_real[n * n];
1072
1073 if (eigVec.is_empty())
1074 {
1075 do_output(LOG_CAT_WARNING, LOG_AREA_SCF, "Failed to output density image for the eigenvector.");
1076 return;
1077 }
1078 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Creating density image for the eigenvector.");
1079 std::vector<ergo_real> vec_perm(n);
1080 eigVec.fullvector(vec_perm);
1081 // now we have the permuted vector
1082 std::vector<ergo_real> vec(n);
1083 for (int ind = 0; ind < n; ind++)
1084 {
1085 vec[ind] = vec_perm[matOpts.inversePermutationHML[ind]];
1086 }
1087
1088 // create density matrix corresponding to the given orbital
1089 std::vector<ergo_real> densityMatrixFull(n *n);
1090
1091 for (int j = 0; j < n; ++j)
1092 {
1093 for (int i = 0; i < n; ++i)
1094 {
1095 densityMatrixFull[i + j * n] = vec[i] * vec[j];
1096 }
1097 }
1098
1099 // Get full matrix version of density matrix, and empty spin density matrix.
1100 {
1101 for (int i = 0; i < n * n; i++)
1102 {
1103 densityMatrixFull_tot [i] = densityMatrixFull[i];
1104 densityMatrixFull_spin[i] = 0;
1105 }
1106 }
1107
1108 do_density_images(basisInfo,
1109 molecule,
1110 densityMatrixFull_tot,
1111 densityMatrixFull_spin,
1112 scfopts.output_density_images_boxwidth,
1113 filename_id); // add filename id
1114
1115 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_restricted::output_density_images_orbital finished OK.");
1116 timeMeter.print(LOG_AREA_SCF, "SCF_restricted::output_density_images_orbital");
1117 }
1118
1119
write_diag_dens_to_file()1120 void SCF_restricted::write_diag_dens_to_file()
1121 {
1122 int n = basisInfo.noOfBasisFuncs;
1123
1124 densityMatrix.readFromFile();
1125 write_diag_elements_to_file(n, densityMatrix, "diagdens.txt",
1126 matOpts.permutationHML);
1127 densityMatrix.writeToFile();
1128 }
1129
1130
report_final_results()1131 void SCF_restricted::report_final_results()
1132 {
1133 }
1134
1135
do_spin_flip(int atomCount)1136 void SCF_restricted::do_spin_flip(int atomCount)
1137 {
1138 throw "error: SCF_restricted::do_spin_flip does not make sense, should not have been called.";
1139 }
1140
1141
save_density_as_prevdens()1142 void SCF_restricted::save_density_as_prevdens()
1143 {
1144 densityMatrix.readFromFile();
1145 Dprev.readFromFile();
1146 Dprev = densityMatrix;
1147 densityMatrix.writeToFile();
1148 Dprev.writeToFile();
1149 }
1150
1151
report_density_difference()1152 void SCF_restricted::report_density_difference()
1153 {
1154 if (scfopts.do_report_density_diff == 0)
1155 {
1156 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
1157 "SCF_restricted::report_density_difference() skipping: (scfopts.do_report_density_diff == 0).");
1158 return;
1159 }
1160 Util::TimeMeter tm;
1161 densityMatrix.readFromFile();
1162 Dprev.readFromFile();
1163 symmMatrix diff(densityMatrix);
1164 diff += (ergo_real) - 1.0 * Dprev;
1165 ergo_real diff_eucl = GetEuclideanNormOfMatrix(diff);
1166 densityMatrix.writeToFile();
1167 Dprev.writeToFile();
1168 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
1169 "SCF_restricted::report_density_difference, diff_eucl = %22.11f", (double)diff_eucl);
1170 tm.print(LOG_AREA_SCF, "SCF_restricted::report_density_difference");
1171 }
1172
1173
compute_dipole_moment()1174 void SCF_restricted::compute_dipole_moment()
1175 {
1176 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_restricted::compute_dipole_moment");
1177 densityMatrix.readFromFile();
1178 get_dipole_moment(densityMatrix, basisInfo, matOpts.size_block_info, matOpts.permutationHML, molecule, LOG_AREA_SCF, "SCF");
1179 densityMatrix.writeToFile();
1180 }
1181
1182
do_mulliken_pop_stuff()1183 void SCF_restricted::do_mulliken_pop_stuff()
1184 {
1185 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_restricted::do_mulliken_pop_stuff");
1186 densityMatrix.readFromFile();
1187 S_symm.readFromFile();
1188 do_mulliken_atomic_charges(densityMatrix,
1189 S_symm,
1190 basisInfo,
1191 matOpts.size_block_info,
1192 matOpts.permutationHML,
1193 matOpts.inversePermutationHML,
1194 molecule);
1195 densityMatrix.writeToFile();
1196 S_symm.writeToFile();
1197 }
1198
1199
create_mtx_files_F(int const scfIter)1200 void SCF_restricted::create_mtx_files_F(int const scfIter)
1201 {
1202 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Creating mtx file for Fock matrix");
1203 std::stringstream ss_fileName;
1204 ss_fileName << "F_matrix_" << scfIter;
1205 std::stringstream ss_id;
1206 ss_id << scfopts.calculation_identifier << " - effective Hamiltonian matrix, SCF cycle " << scfIter;
1207 FockMatrix.readFromFile();
1208 write_matrix_in_matrix_market_format(FockMatrix, matOpts.inversePermutationHML, ss_fileName.str(),
1209 ss_id.str(), scfopts.method_and_basis_set);
1210 FockMatrix.writeToFile();
1211 }
1212
1213
create_mtx_files_D(int const scfIter)1214 void SCF_restricted::create_mtx_files_D(int const scfIter)
1215 {
1216 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Creating mtx file for density matrix");
1217 std::stringstream ss_fileName;
1218 ss_fileName << "D_matrix_" << scfIter;
1219 std::stringstream ss_id;
1220 ss_id << scfopts.calculation_identifier << " - density matrix, SCF cycle " << scfIter;
1221 densityMatrix.readFromFile();
1222 write_matrix_in_matrix_market_format(densityMatrix, matOpts.inversePermutationHML, ss_fileName.str(),
1223 ss_id.str(), scfopts.method_and_basis_set);
1224 densityMatrix.writeToFile();
1225 }
1226
1227
create_eigenvalues_files() const1228 void SCF_restricted::create_eigenvalues_files() const
1229 {
1230 if(! eigValOCC.empty() )
1231 {
1232 std::ofstream ff("occupied_spectrum.txt");
1233 for (size_t ind = 0; ind < eigValOCC.size(); ind++)
1234 {
1235 ff << (double)eigValOCC[ind] << std::endl;
1236 }
1237 ff.close();
1238 }
1239
1240 if(! eigValUNOCC.empty() )
1241 {
1242 std::ofstream ff("unoccupied_spectrum.txt");
1243 for (size_t ind = 0; ind < eigValUNOCC.size(); ind++)
1244 {
1245 ff << (double)eigValUNOCC[ind] << std::endl;
1246 }
1247 ff.close();
1248 }
1249 }
1250
1251
create_eigenvectors_files() const1252 void SCF_restricted::create_eigenvectors_files() const
1253 {
1254 create_eigenvalues_files();
1255
1256
1257 if(! eigVecOCC.empty() )
1258 {
1259 create_eigvec_file(eigVecOCC[0],
1260 "HOMO",
1261 "homo_coefficient_vec");
1262 for (size_t i = 1; i < eigVecOCC.size(); i++)
1263 {
1264 std::stringstream ss_name, ss_filename;
1265 ss_name << "HOMO-" << i;
1266 ss_filename << "occ_" << i << "_coefficient_vec";
1267 create_eigvec_file(eigVecOCC[i],
1268 ss_name.str().c_str(),
1269 ss_filename.str().c_str());
1270 }
1271 }
1272
1273 if(! eigVecUNOCC.empty())
1274 {
1275 create_eigvec_file(eigVecUNOCC[0],
1276 "LUMO",
1277 "lumo_coefficient_vec");
1278 for (size_t i = 1; i < eigVecUNOCC.size(); i++)
1279 {
1280 std::stringstream ss_name, ss_filename;
1281 ss_name << "LUMO-" << i;
1282 ss_filename << "unocc_" << i << "_coefficient_vec";
1283 create_eigvec_file(eigVecUNOCC[i],
1284 ss_name.str().c_str(),
1285 ss_filename.str().c_str());
1286 }
1287 }
1288 }
1289
1290
create_eigvec_file(const generalVector & eigVec,const char * vector_name,const char * filename_id) const1291 void SCF_restricted::create_eigvec_file(const generalVector& eigVec,
1292 const char *vector_name,
1293 const char *filename_id) const
1294 {
1295 if (eigVec.is_empty())
1296 {
1297 do_output(LOG_CAT_WARNING, LOG_AREA_SCF, "Failed to output %s to file. No %s eigenvector stored.", vector_name, vector_name);
1298 return;
1299 }
1300 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Storing %s eigenvector to file %s.", vector_name, filename_id);
1301 int n = basisInfo.noOfBasisFuncs;
1302 std::vector<ergo_real> vec_perm(n);
1303 eigVec.fullvector(vec_perm);
1304 // now we have the permuted vector
1305 std::vector<ergo_real> vec(n);
1306 for (int ind = 0; ind < n; ind++)
1307 {
1308 vec[matOpts.inversePermutationHML[ind]] = vec_perm[ind];
1309 }
1310 char ffname[888];
1311 sprintf(ffname, "%s.txt", filename_id);
1312 std::ofstream ff(ffname);
1313 for (int ind = 0; ind < n; ind++)
1314 {
1315 ff << (double)vec[ind] << std::endl;
1316 }
1317 ff.close();
1318 }
1319
1320
1321 static void
output_orbital_coeffs_in_gabedit_order(const BasisInfoStruct & basisInfo,std::vector<int> const & shellIdxList,std::ofstream & ff,std::vector<ergo_real> const & orbital_vec)1322 output_orbital_coeffs_in_gabedit_order(const BasisInfoStruct& basisInfo,
1323 std::vector<int> const& shellIdxList,
1324 std::ofstream& ff,
1325 std::vector<ergo_real> const& orbital_vec)
1326 {
1327 int count = 0;
1328
1329 for (int i = 0; i < basisInfo.noOfShells; i++)
1330 {
1331 int k = shellIdxList[i];
1332 int startIdx = basisInfo.shellList[k].startIndexInMatrix;
1333 switch (basisInfo.shellList[k].shellType)
1334 {
1335 case 0: // s-type shell
1336 ff << count + 1 << " " << (double)orbital_vec[startIdx] << std::endl;
1337 count++;
1338 break;
1339
1340 case 1: // p-type shell
1341 ff << count + 1 << " " << (double)orbital_vec[startIdx + 2] << std::endl;
1342 count++;
1343 ff << count + 1 << " " << (double)orbital_vec[startIdx + 0] << std::endl;
1344 count++;
1345 ff << count + 1 << " " << (double)orbital_vec[startIdx + 1] << std::endl;
1346 count++;
1347 break;
1348
1349 case 2: // d-type shell
1350 ff << count + 1 << " " << (double)orbital_vec[startIdx + 2] << std::endl;
1351 count++;
1352 ff << count + 1 << " " << (double)orbital_vec[startIdx + 3] << std::endl;
1353 count++;
1354 ff << count + 1 << " " << (double)orbital_vec[startIdx + 1] << std::endl;
1355 count++;
1356 ff << count + 1 << " " << (double)orbital_vec[startIdx + 4] << std::endl;
1357 count++;
1358 ff << count + 1 << " " << (double)orbital_vec[startIdx + 0] << std::endl;
1359 count++;
1360 break;
1361
1362 case 3: // f-type shell
1363 ff << count + 1 << " " << (double)orbital_vec[startIdx + 3] << std::endl;
1364 count++;
1365 ff << count + 1 << " " << (double)orbital_vec[startIdx + 4] << std::endl;
1366 count++;
1367 ff << count + 1 << " " << (double)orbital_vec[startIdx + 2] << std::endl;
1368 count++;
1369 ff << count + 1 << " " << (double)orbital_vec[startIdx + 5] << std::endl;
1370 count++;
1371 ff << count + 1 << " " << (double)orbital_vec[startIdx + 1] << std::endl;
1372 count++;
1373 ff << count + 1 << " " << (double)orbital_vec[startIdx + 6] << std::endl;
1374 count++;
1375 ff << count + 1 << " " << (double)orbital_vec[startIdx + 0] << std::endl;
1376 count++;
1377 break;
1378
1379 case 4: // g-type shell
1380 ff << count + 1 << " " << (double)orbital_vec[startIdx + 4] << std::endl;
1381 count++;
1382 ff << count + 1 << " " << (double)orbital_vec[startIdx + 5] << std::endl;
1383 count++;
1384 ff << count + 1 << " " << (double)orbital_vec[startIdx + 3] << std::endl;
1385 count++;
1386 ff << count + 1 << " " << (double)orbital_vec[startIdx + 6] << std::endl;
1387 count++;
1388 ff << count + 1 << " " << (double)orbital_vec[startIdx + 2] << std::endl;
1389 count++;
1390 ff << count + 1 << " " << (double)orbital_vec[startIdx + 7] << std::endl;
1391 count++;
1392 ff << count + 1 << " " << (double)orbital_vec[startIdx + 1] << std::endl;
1393 count++;
1394 ff << count + 1 << " " << (double)orbital_vec[startIdx + 8] << std::endl;
1395 count++;
1396 ff << count + 1 << " " << (double)orbital_vec[startIdx + 0] << std::endl;
1397 count++;
1398 break;
1399
1400 default:
1401 throw "error in output_orbital_coeffs_in_gabedit_order: shell types beyond g not implemented!";
1402 }
1403 }
1404 if (count != basisInfo.noOfBasisFuncs)
1405 {
1406 throw "error in output_orbital_coeffs_in_gabedit_order: (count != basisInfo.noOfBasisFuncs)";
1407 }
1408 }
1409
1410
create_gabedit_file() const1411 void SCF_restricted::create_gabedit_file() const
1412 {
1413 if (eigVecOCC.empty() || eigVecUNOCC.empty())
1414 {
1415 do_output(LOG_CAT_WARNING, LOG_AREA_SCF, "Failed to output eigenvectors to gabedit file; no eigenvectors info available.");
1416 return;
1417 }
1418 if (basisInfo.use_6_d_funcs == 1)
1419 {
1420 do_output(LOG_CAT_WARNING, LOG_AREA_SCF, "Failed to output eigenvectors to gabedit file; not implemented for use_6_d_funcs case.");
1421 return;
1422 }
1423 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Creating Gabedit file with eigenvector info.");
1424 int n = basisInfo.noOfBasisFuncs;
1425
1426 // Create Gabedit file.
1427 const char fileName [] = "gabeditfile.gab";
1428 std::ofstream ff(fileName);
1429
1430 /* FIXME: check if we should use "Cart" or "Sphe" here. That is,
1431 * should we use use_6_d_funcs? */
1432 int use_6_d_funcs = 0;
1433 ff << "[Gabedit Format] Sphe" << std::endl;
1434 ff << "[Atoms] Angs" << std::endl;
1435 for (int i = 0; i < molecule.getNoOfAtoms(); i++)
1436 {
1437 char atomLabelString[4];
1438 get_atom_label_from_charge_int(molecule.getAtom(i).charge, atomLabelString, 4);
1439 ff << atomLabelString << " " << i + 1 << " " << (double)molecule.getAtom(i).charge
1440 << " " << (double)(molecule.getAtom(i).coords[0] / UNIT_one_Angstrom)
1441 << " " << (double)(molecule.getAtom(i).coords[1] / UNIT_one_Angstrom)
1442 << " " << (double)(molecule.getAtom(i).coords[2] / UNIT_one_Angstrom)
1443 << std::endl;
1444 }
1445 std::vector<int> shellIdxList(basisInfo.noOfShells);
1446 int shellIdxCounter = 0;
1447 SquareFuncIntegrator sfi;
1448 ff << "[Basis]" << std::endl;
1449 for (int i = 0; i < molecule.getNoOfAtoms(); i++)
1450 {
1451 ff << i + 1 << " 0" << std::endl;
1452 // Now output info about shells for this atom.
1453 for (int k = 0; k < basisInfo.noOfShells; k++)
1454 {
1455 // Check if this shell belongs to the current atom.
1456 ergo_real absdx = template_blas_fabs(basisInfo.shellList[k].centerCoords[0] - molecule.getAtom(i).coords[0]);
1457 ergo_real absdy = template_blas_fabs(basisInfo.shellList[k].centerCoords[1] - molecule.getAtom(i).coords[1]);
1458 ergo_real absdz = template_blas_fabs(basisInfo.shellList[k].centerCoords[2] - molecule.getAtom(i).coords[2]);
1459 ergo_real distlimit = 0.01;
1460 if ((absdx > distlimit) || (absdy > distlimit) || (absdz > distlimit))
1461 {
1462 continue;
1463 }
1464 // OK, now we know this shell is at least very near the current atom.
1465 shellIdxList[shellIdxCounter] = k;
1466 shellIdxCounter++;
1467 char shellChar = 'x';
1468 int shellType = basisInfo.shellList[k].shellType;
1469 switch (shellType)
1470 {
1471 case 0:
1472 shellChar = 's';
1473 break;
1474
1475 case 1:
1476 shellChar = 'p';
1477 break;
1478
1479 case 2:
1480 shellChar = 'd';
1481 break;
1482
1483 case 3:
1484 shellChar = 'f';
1485 break;
1486
1487 case 4:
1488 shellChar = 'g';
1489 break;
1490
1491 default:
1492 throw "SCF_restricted::create_gabedit_file error: shell types beyond g not implemented!";
1493 }
1494 ff << shellChar << " " << basisInfo.shellList[k].noOfContr << " 1.00" << std::endl;
1495 for (int contridx = 0; contridx < basisInfo.shellList[k].noOfContr; contridx++)
1496 {
1497 ergo_real exponent = basisInfo.shellList[k].exponentList[contridx];
1498 ergo_real shellFactor = sfi.getShellFactor(integralInfo, exponent, shellType, use_6_d_funcs);
1499 ergo_real coeff = basisInfo.shellList[k].coeffList[contridx] / shellFactor;
1500 ff << (double)exponent << " " << (double)coeff << std::endl;
1501 }
1502 }
1503 // Blank line before shells for next atom.
1504 ff << std::endl;
1505 }
1506 if (shellIdxCounter != basisInfo.noOfShells)
1507 {
1508 throw "Error: (shellIdxCounter != basisInfo.noOfShells)";
1509 }
1510 // MO section.
1511 ff << "[MO]" << std::endl;
1512
1513
1514 // // HOMO
1515 // ff << "Spin=Alpha" << std::endl;
1516 // ff << "Occup= 2.000000" << std::endl;
1517 // output_orbital_coeffs_in_gabedit_order(basisInfo, shellIdxList, ff, homo_vec);
1518 // // LUMO
1519 // ff << "Spin=Alpha" << std::endl;
1520 // ff << "Occup= 0.000000" << std::endl;
1521 // output_orbital_coeffs_in_gabedit_order(basisInfo, shellIdxList, ff, lumo_vec);
1522
1523
1524
1525 // Occupied orbitals
1526 for (int orb = (int)eigVecOCC.size()-1; orb >= 0; orb--) // note that orb can be negative
1527 {
1528 std::vector<ergo_real> orb_vec_perm(n);
1529 eigVecOCC[orb].fullvector(orb_vec_perm);
1530 // now we have the permuted vector
1531 std::vector<ergo_real> orb_vec(n);
1532 for (int ind = 0; ind < n; ind++)
1533 {
1534 orb_vec[matOpts.inversePermutationHML[ind]] = orb_vec_perm[ind];
1535 }
1536 ff << "Spin=Alpha" << std::endl;
1537 if((int)eigValOCC.size() > orb) // in case we do not have eigenvalues
1538 ff << "Energy= " << (double)eigValOCC[orb] << std::endl;
1539 ff << "Occup= 2.000000" << std::endl;
1540 output_orbital_coeffs_in_gabedit_order(basisInfo, shellIdxList, ff, orb_vec);
1541 }
1542
1543
1544 // Unoccupied orbitals
1545 for (size_t orb = 0; orb < eigVecUNOCC.size(); orb++)
1546 {
1547 std::vector<ergo_real> orb_vec_perm(n);
1548 eigVecUNOCC[orb].fullvector(orb_vec_perm);
1549 // now we have the permuted vector
1550 std::vector<ergo_real> orb_vec(n);
1551 for (int ind = 0; ind < n; ind++)
1552 {
1553 orb_vec[matOpts.inversePermutationHML[ind]] = orb_vec_perm[ind];
1554 }
1555 ff << "Spin=Alpha" << std::endl;
1556 if(eigValUNOCC.size() > orb) // in case we do not have eigenvalues
1557 ff << "Energy= " << (double)eigValUNOCC[orb] << std::endl;
1558 ff << "Occup= 0.000000" << std::endl;
1559 output_orbital_coeffs_in_gabedit_order(basisInfo, shellIdxList, ff, orb_vec);
1560 }
1561
1562
1563 // Blank line before end of file.
1564 ff << std::endl;
1565 // Close file.
1566 ff.close();
1567 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Gabedit file '%s' with eigenvectors info created OK.", fileName);
1568 }
1569
1570
update_subspace_diff()1571 void SCF_restricted::update_subspace_diff()
1572 {
1573 densityMatrix.readFromFile();
1574 Dprev.readFromFile();
1575 ergo_real acc = template_blas_sqrt(get_machine_epsilon());
1576
1577 symmMatrix diff(densityMatrix);
1578 diff += (ergo_real) - 1.0 * Dprev;
1579
1580 transform_with_S(diff);
1581 transform_with_invChol(diff);
1582
1583 // Compensate for factor 2 (restricted case)
1584 diff *= (ergo_real)0.5;
1585
1586 ergo_real diff_eucl = diff.eucl(acc);
1587
1588 densityMatrix.writeToFile();
1589 Dprev.writeToFile();
1590
1591 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_restricted::update_subspace_diff, diff_eucl = %22.11f", (double)diff_eucl);
1592 curr_subspace_diff = diff_eucl;
1593 }
1594
1595
1596 struct RandomNumber
1597 {
accumulateRandomNumber1598 ergo_real accumulate(ergo_real& a, int const dummy1, int const dummy2)
1599 {
1600 a = rand() / (ergo_real)RAND_MAX;
1601 return 0;
1602 }
1603 };
1604
1605
1606 /** Transform matrix A to S*A*S */
transform_with_S(symmMatrix & A)1607 void SCF_restricted::transform_with_S(symmMatrix& A)
1608 {
1609 S_symm.readFromFile();
1610
1611 normalMatrix S_norm(S_symm);
1612 normalMatrix A_norm(A);
1613
1614 normalMatrix SA(S_symm);
1615 SA = (ergo_real)1.0 * S_norm * A_norm;
1616 normalMatrix SAS(S_symm);
1617 SAS = (ergo_real)1.0 * SA * S_norm;
1618
1619 A = SAS;
1620
1621 S_symm.writeToFile();
1622 }
1623
1624
1625 /** Transform matrix A to invCholT*A*invChol */
transform_with_invChol(symmMatrix & A)1626 void SCF_restricted::transform_with_invChol(symmMatrix& A)
1627 {
1628 invCholFactor.readFromFile();
1629 A = transpose(invCholFactor) * A * invCholFactor;
1630 invCholFactor.writeToFile();
1631 }
1632
1633
get_non_ort_err_mat_normalized_in_ort_basis(symmMatrix & randomMatrix,int transform_with_S_also)1634 void SCF_restricted::get_non_ort_err_mat_normalized_in_ort_basis(symmMatrix& randomMatrix, int transform_with_S_also)
1635 {
1636 symmMatrix randomMatrix1;
1637
1638 randomMatrix1.resetSizesAndBlocks(matOpts.size_block_info,
1639 matOpts.size_block_info);
1640 symmMatrix randomMatrix2;
1641 randomMatrix2.resetSizesAndBlocks(matOpts.size_block_info,
1642 matOpts.size_block_info);
1643 randomMatrix1.random();
1644 randomMatrix2.random();
1645 randomMatrix = 0;
1646 randomMatrix += (ergo_real)1.0 * randomMatrix1;
1647 randomMatrix += (ergo_real) - 1.0 * randomMatrix2;
1648
1649 symmMatrix randomMatrix_ort(randomMatrix);
1650
1651 if (transform_with_S_also)
1652 {
1653 transform_with_S(randomMatrix_ort);
1654 }
1655
1656 transform_with_invChol(randomMatrix_ort);
1657
1658 ergo_real acc = template_blas_sqrt(get_machine_epsilon());
1659 ergo_real randomMatrix_Norm = randomMatrix.eucl(acc);
1660 ergo_real randomMatrix_ort_Norm = randomMatrix_ort.eucl(acc);
1661
1662 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "norms of randomMatrix and randomMatrix_ort : %22.11f %22.11f",
1663 (double)randomMatrix_Norm, (double)randomMatrix_ort_Norm);
1664
1665 // Normalize randomMatrix so that randomMatrix_ort would have norm 1.
1666 randomMatrix *= (ergo_real)(1.0 / randomMatrix_ort_Norm);
1667
1668 ergo_real randomMatrixNormAfterNormalization = randomMatrix.eucl(acc);
1669 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "randomMatrixNormAfterNormalization = %22.11f", (double)randomMatrixNormAfterNormalization);
1670 }
1671
1672
disturb_dens_matrix(ergo_real subspaceError)1673 void SCF_restricted::disturb_dens_matrix(ergo_real subspaceError)
1674 {
1675 ergo_real gap = 2;
1676 ergo_real desiredErrorNorm = gap * subspaceError / (1 + subspaceError);
1677
1678 symmMatrix randomMatrix;
1679
1680 randomMatrix.resetSizesAndBlocks(matOpts.size_block_info,
1681 matOpts.size_block_info);
1682 get_non_ort_err_mat_normalized_in_ort_basis(randomMatrix, 1);
1683
1684 densityMatrix.readFromFile();
1685
1686 densityMatrix += desiredErrorNorm * randomMatrix;
1687
1688 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Disturbed density matrix with desiredErrorNorm = %22.11f (gap = %8.3f)",
1689 (double)desiredErrorNorm, (double)gap);
1690
1691 densityMatrix.writeToFile();
1692 }
1693
1694
disturb_dens_matrix_exact_try(const symmMatrix & randomMatrix,const symmMatrix & orgDensMatrix,ergo_real disturbanceFactor,ergo_real & resultSinTheta,symmMatrix & resultDensMatrix)1695 void SCF_restricted::disturb_dens_matrix_exact_try(const symmMatrix& randomMatrix,
1696 const symmMatrix& orgDensMatrix,
1697 ergo_real disturbanceFactor,
1698 ergo_real& resultSinTheta,
1699 symmMatrix& resultDensMatrix)
1700 {
1701 symmMatrix D(orgDensMatrix);
1702
1703 D += (ergo_real)1.0 * disturbanceFactor * randomMatrix;
1704
1705 symmMatrix SDS_symm(D);
1706 transform_with_S(SDS_symm);
1707 SDS_symm *= (ergo_real) - 1.0;
1708
1709 symmMatrix F_ort_prev_dummy;
1710 F_ort_prev_dummy.resetSizesAndBlocks(matOpts.size_block_info,
1711 matOpts.size_block_info);
1712 F_ort_prev_dummy.writeToFile();
1713
1714 resultDensMatrix.writeToFile();
1715 SDS_symm.writeToFile();
1716
1717
1718 int use_diag = DensFromFock.get_use_diagonalization();
1719 int use_diag_on_error = DensFromFock.get_use_diag_on_error();
1720
1721 DensFromFock.unset_use_diagonalization();
1722 DensFromFock.unset_use_diag_on_error();
1723
1724 DensFromFock.clean_eigs_intervals();
1725
1726 if (DensFromFock.get_dens_from_fock(SDS_symm,
1727 resultDensMatrix,
1728 F_ort_prev_dummy) != 0)
1729 {
1730 throw "SCF_restricted::disturb_dens_matrix_exact_try: Error in get_dens_from_fock";
1731 }
1732
1733 if (use_diag == 1)
1734 {
1735 DensFromFock.set_use_diagonalization();
1736 }
1737 if (use_diag_on_error == 1)
1738 {
1739 DensFromFock.set_use_diag_on_error();
1740 }
1741
1742
1743 // OK, now we have computed D_Pure which is the purified version of SDS_symm.
1744 // But D_Pure is not in orthogonal basis.
1745
1746 resultDensMatrix.readFromFile();
1747
1748 symmMatrix diff(resultDensMatrix);
1749 diff += (ergo_real) - 1.0 * orgDensMatrix;
1750
1751 transform_with_S(diff);
1752 transform_with_invChol(diff);
1753
1754 // Compensate for factor 2 (restricted case)
1755 diff *= (ergo_real)0.5;
1756
1757 ergo_real acc = template_blas_sqrt(get_machine_epsilon());
1758
1759 resultSinTheta = diff.eucl(acc);
1760 }
1761
1762
disturb_dens_matrix_exact(ergo_real subspaceError)1763 void SCF_restricted::disturb_dens_matrix_exact(ergo_real subspaceError)
1764 {
1765 //ergo_real gap = 2;
1766 //ergo_real desiredErrorNorm = gap * subspaceError / (1 + subspaceError);
1767
1768 symmMatrix randomMatrix;
1769
1770 randomMatrix.resetSizesAndBlocks(matOpts.size_block_info,
1771 matOpts.size_block_info);
1772 get_non_ort_err_mat_normalized_in_ort_basis(randomMatrix, 1);
1773
1774 symmMatrix newDensMatrix;
1775 newDensMatrix.resetSizesAndBlocks(matOpts.size_block_info,
1776 matOpts.size_block_info);
1777
1778 densityMatrix.readFromFile();
1779
1780 ergo_real currSinTheta;
1781 ergo_real disturbanceFactor_min = 0;
1782 ergo_real disturbanceFactor_max = 5;
1783
1784 int iterCount = 0;
1785 do
1786 {
1787 iterCount++;
1788 if (iterCount > 44)
1789 {
1790 throw "error in SCF_restricted::disturb_dens_matrix_exact, iterCount esceeded limit.";
1791 }
1792 ergo_real disturbanceFactor = 0.5 * (disturbanceFactor_min + disturbanceFactor_max);
1793 disturb_dens_matrix_exact_try(randomMatrix,
1794 densityMatrix,
1795 disturbanceFactor,
1796 currSinTheta,
1797 newDensMatrix);
1798 if (currSinTheta < subspaceError)
1799 {
1800 disturbanceFactor_min = disturbanceFactor;
1801 }
1802 else
1803 {
1804 disturbanceFactor_max = disturbanceFactor;
1805 }
1806 } while (template_blas_fabs(currSinTheta - subspaceError) > 0.001 * subspaceError);
1807
1808 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_restricted::disturb_dens_matrix_exact done, iterCount = %2i", iterCount);
1809
1810 densityMatrix = newDensMatrix;
1811
1812 densityMatrix.writeToFile();
1813 }
1814
1815
disturb_fock_matrix(ergo_real subspaceError)1816 void SCF_restricted::disturb_fock_matrix(ergo_real subspaceError)
1817 {
1818 symmMatrix F_ort_prev_dummy;
1819
1820 F_ort_prev_dummy.resetSizesAndBlocks(matOpts.size_block_info,
1821 matOpts.size_block_info);
1822 F_ort_prev_dummy.writeToFile();
1823
1824 symmMatrix densityMatrix_dummy;
1825 densityMatrix_dummy.resetSizesAndBlocks(matOpts.size_block_info,
1826 matOpts.size_block_info);
1827 densityMatrix_dummy.writeToFile();
1828
1829
1830
1831 int use_diag = DensFromFock.get_use_diagonalization();
1832 int use_diag_on_error = DensFromFock.get_use_diag_on_error();
1833
1834 DensFromFock.unset_use_diagonalization();
1835 DensFromFock.unset_use_diag_on_error();
1836
1837 DensFromFock.clean_eigs_intervals();
1838
1839 if (DensFromFock.get_dens_from_fock(FockMatrix,
1840 densityMatrix_dummy,
1841 F_ort_prev_dummy) != 0)
1842 {
1843 throw "SCF_restricted::disturb_fock_matrix: Error in get_dens_from_fock";
1844 }
1845
1846 if (use_diag == 1)
1847 {
1848 DensFromFock.set_use_diagonalization();
1849 }
1850 if (use_diag_on_error == 1)
1851 {
1852 DensFromFock.set_use_diag_on_error();
1853 }
1854
1855 intervalType homoInterval_tmp2;
1856 intervalType lumoInterval_tmp2;
1857 DensFromFock.get_eigs_F_ort_prev(homoInterval_tmp2, lumoInterval_tmp2);
1858
1859
1860 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_restricted::disturb_fock_matrix, interval sizes: %22.11f %22.11f",
1861 (double)(homoInterval_tmp2.upp() - homoInterval_tmp2.low()), (double)(lumoInterval_tmp2.upp() - lumoInterval_tmp2.low()));
1862
1863 ergo_real gap = lumoInterval_tmp2.low() - homoInterval_tmp2.upp();
1864
1865 ergo_real desiredErrorNorm = gap * subspaceError / (1 + subspaceError);
1866
1867 symmMatrix randomMatrix;
1868 randomMatrix.resetSizesAndBlocks(matOpts.size_block_info,
1869 matOpts.size_block_info);
1870 get_non_ort_err_mat_normalized_in_ort_basis(randomMatrix, 0);
1871
1872 FockMatrix.readFromFile();
1873
1874 FockMatrix += desiredErrorNorm * randomMatrix;
1875
1876 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Disturbed Fock matrix with desiredErrorNorm = %22.11f (gap = %8.3f)",
1877 (double)desiredErrorNorm, (double)gap);
1878
1879 FockMatrix.writeToFile();
1880 }
1881
1882
get_nucl_energy_for_given_mol_and_dens(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)1883 static ergo_real get_nucl_energy_for_given_mol_and_dens(const IntegralInfo& integralInfo,
1884 const Molecule& molecule,
1885 const BasisInfoStruct& basisInfo,
1886 const symmMatrix& D,
1887 ergo_real threshold_integrals_1el,
1888 mat::SizesAndBlocks const& matrix_size_block_info,
1889 std::vector<int> const& permutationHML)
1890 {
1891 ergo_real nuclearRepulsionEnergy = molecule.getNuclearRepulsionEnergyQuadratic();
1892 ergo_real elecNuclEnergy = get_electron_nuclear_attraction_energy(integralInfo,
1893 molecule,
1894 basisInfo,
1895 D,
1896 threshold_integrals_1el,
1897 matrix_size_block_info,
1898 permutationHML);
1899
1900 return nuclearRepulsionEnergy + elecNuclEnergy;
1901 }
1902
1903
1904 /* Compute gradient of energy with respect to nuclear positions, for
1905 * fixed electron density. */
1906 void
compute_gradient_fixeddens()1907 SCF_restricted::compute_gradient_fixeddens()
1908 {
1909 // Since we here regard the electron density as fixed, there are
1910 // only two terms in the energy which give nonzero contributions:
1911 // the nuclear-electron interaction term, and the nuclear-nuclear
1912 // interaction term.
1913 densityMatrix.readFromFile();
1914 int nAtoms = molecule.getNoOfAtoms();
1915 std::vector<ergo_real> gradient(nAtoms * 3);
1916 get_gradient_for_given_mol_and_dens(integralInfo,
1917 molecule,
1918 basisInfo,
1919 densityMatrix,
1920 threshold_integrals_1el,
1921 matOpts.size_block_info,
1922 matOpts.permutationHML,
1923 &gradient[0]);
1924 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Gradient of energy with respect to nuclear positions, "
1925 "for fixed electron density (no Pulay correction, only Hellmann-Feynman forces):");
1926 for (int i = 0; i < nAtoms; i++)
1927 {
1928 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Atom %6d: %22.11f %22.11f %22.11f",
1929 i,
1930 (double)gradient[i * 3 + 0],
1931 (double)gradient[i * 3 + 1],
1932 (double)gradient[i * 3 + 2]);
1933 }
1934 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "(End of gradient)");
1935
1936 if (scfopts.verify_gradient_fixeddens == 1)
1937 {
1938 std::vector<ergo_real> gradient_for_verification(nAtoms * 3);
1939 for (int i = 0; i < nAtoms; i++)
1940 {
1941 for (int coordIdx = 0; coordIdx < 3; coordIdx++)
1942 {
1943 const ergo_real h = 1e-3;
1944 Molecule moleculeTmp = molecule;
1945 Atom atomTmp = molecule.getAtom(i);
1946 atomTmp.coords[coordIdx] += h;
1947 moleculeTmp.replaceAtom(i, atomTmp);
1948 ergo_real E1 = get_nucl_energy_for_given_mol_and_dens(integralInfo,
1949 moleculeTmp,
1950 basisInfo,
1951 densityMatrix,
1952 threshold_integrals_1el,
1953 matOpts.size_block_info,
1954 matOpts.permutationHML);
1955 moleculeTmp = molecule;
1956 atomTmp = molecule.getAtom(i);
1957 atomTmp.coords[coordIdx] -= h;
1958 moleculeTmp.replaceAtom(i, atomTmp);
1959 ergo_real E2 = get_nucl_energy_for_given_mol_and_dens(integralInfo,
1960 moleculeTmp,
1961 basisInfo,
1962 densityMatrix,
1963 threshold_integrals_1el,
1964 matOpts.size_block_info,
1965 matOpts.permutationHML);
1966 ergo_real gradientComponent = (E1 - E2) / (2 * h);
1967 gradient_for_verification[i * 3 + coordIdx] = gradientComponent;
1968 } // END FOR coordIdx
1969 } // END FOR i
1970 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Gradient of energy with respect to nuclear positions, for fixed electron density (for verification, computed using finite differences):");
1971 for (int i = 0; i < nAtoms; i++)
1972 {
1973 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Atom %6d: %22.11f %22.11f %22.11f",
1974 i,
1975 (double)gradient_for_verification[i * 3 + 0],
1976 (double)gradient_for_verification[i * 3 + 1],
1977 (double)gradient_for_verification[i * 3 + 2]);
1978 }
1979 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "(End of gradient)");
1980 ergo_real maxAbsDiff = -1;
1981 int atomIdx_saved = -1;
1982 int coordIdx_saved = -1;
1983 for (int atomIdx = 0; atomIdx < nAtoms; atomIdx++)
1984 {
1985 for (int coordIdx = 0; coordIdx < 3; coordIdx++)
1986 {
1987 ergo_real absdiff = template_blas_fabs(gradient[atomIdx * 3 + coordIdx] - gradient_for_verification[atomIdx * 3 + coordIdx]);
1988 if (absdiff > maxAbsDiff)
1989 {
1990 maxAbsDiff = absdiff;
1991 atomIdx_saved = atomIdx;
1992 coordIdx_saved = coordIdx;
1993 }
1994 }
1995 }
1996 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "maxAbsDiff %22.11f = %9.4g, for atomIdx %d and coordIdx %d", maxAbsDiff, maxAbsDiff, atomIdx_saved, coordIdx_saved);
1997 }
1998
1999 densityMatrix.writeToFile();
2000 }
2001