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_general.cc
31
32 @brief Class for self-consistent field (SCF) procedure; base class
33 that can be used for both restricted and unrestricted cases.
34
35 @author: Elias Rudberg <em>responsible</em>.
36 */
37
38 #include <sstream>
39 #include <sys/types.h>
40 #include <unistd.h>
41 #include "SCF_general.h"
42 #include "output.h"
43 #include "scf_utils.h"
44 #include "matrix_utilities.h"
45 #include "utilities.h"
46 #include "integral_matrix_wrappers.h"
47 #include "machine_epsilon.h"
48 #include "units.h"
49 #include "SCF_statistics.h"
50 #include "AllocatorManager.h"
51
52
get_eucl_norm_try_different_acc(const symmMatrix & A,ergo_real & acc)53 static ergo_real get_eucl_norm_try_different_acc(const symmMatrix & A, ergo_real & acc) {
54 // First try with acc = sqrt(epsilon) and if that is too difficult try with lower accuracy.
55 acc = template_blas_sqrt(mat::getMachineEpsilon<ergo_real>());
56 int maxIter = 1000;
57 ergo_real result_eucl_norm = -1;
58 try {
59 result_eucl_norm = A.eucl(acc, maxIter);
60 return result_eucl_norm;
61 }
62 catch(...) {
63 ergo_real accNew = template_blas_sqrt(acc);
64 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "get_eucl_norm_try_different_acc(): first attempt failed, now reducing accuracy from %g to %g", acc, accNew);
65 acc = accNew;
66 }
67 try {
68 result_eucl_norm = A.eucl(acc, maxIter);
69 return result_eucl_norm;
70 }
71 catch(...) {
72 ergo_real accNew = template_blas_sqrt(acc);
73 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "get_eucl_norm_try_different_acc(): second attempt failed, now reducing accuracy from %g to %g", acc, accNew);
74 acc = accNew;
75 }
76 // Now call eucl without maxIter, let it use as many iterations as it needs to reach the reduced accuracy acc.
77 return A.eucl(acc);
78 }
79
80
SCF_general(const Molecule & molecule_,const Molecule & extraCharges_,const BasisInfoStruct & basisInfo_,const IntegralInfo & integralInfo_,const char * guessDmatFileName_,const JK::Params & J_K_params_,const Dft::GridParams & gridParams_,const SCF::Options & scfoptsPtr,const SCF::MatOptions & matOpts_,ergo_real threshold_integrals_1el_input)81 SCF_general::SCF_general(const Molecule& molecule_,
82 const Molecule& extraCharges_,
83 const BasisInfoStruct & basisInfo_,
84 const IntegralInfo & integralInfo_,
85 const char* guessDmatFileName_,
86 const JK::Params& J_K_params_,
87 const Dft::GridParams& gridParams_,
88 const SCF::Options& scfoptsPtr,
89 const SCF::MatOptions& matOpts_,
90 ergo_real threshold_integrals_1el_input)
91 :
92 molecule(molecule_),
93 extraCharges(extraCharges_),
94 basisInfo(basisInfo_),
95 integralInfo(integralInfo_),
96 guessDmatFileName(guessDmatFileName_), // FIXME: copy this object properly
97 J_K_params(J_K_params_),
98 gridParams(gridParams_),
99 scfopts(scfoptsPtr),
100 matOpts(matOpts_),
101 threshold_integrals_1el(threshold_integrals_1el_input),
102 DIIS(NULL),
103 curr_cycle_stats(NULL)
104 {
105 int n = basisInfo.noOfBasisFuncs;
106 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_general constructor, number of basis functions: %i", n);
107
108 output_current_memory_usage(LOG_AREA_SCF, "beginning of SCF_general constructor");
109
110 // Report info about host name, process ID etc.
111 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "************** Some general info here **********************");
112 do_output_time(LOG_CAT_INFO, LOG_AREA_SCF, "VERSION: " VERSION " time : ");
113 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "machine_epsilon = %9.3g", (double)get_machine_epsilon());
114 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "sizeof(ergo_real) = %i", sizeof(ergo_real));
115 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "sizeof(size_t) = %i", sizeof(size_t));
116 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "sizeof(int) = %i", sizeof(int));
117 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "sizeof(long) = %i", sizeof(long));
118 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "sizeof(char*) = %i", sizeof(char*));
119 host_name_struct hostName;
120 get_host_name(&hostName);
121 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Host name: '%s'", hostName.s);
122 working_directory_struct workingDirectory;
123 get_working_directory(&workingDirectory);
124 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Working directory: '%s'", workingDirectory.s);
125 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Process ID (PID): %10i", getpid());
126 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "************************************************************");
127 ergo_real minDist, maxDist;
128 molecule.getExtremeInternuclearDistancesQuadratic(minDist, maxDist);
129 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Min internuclear distance: %12.5f a.u. = %12.5f Angstrom", (double)minDist, (double)(minDist/UNIT_one_Angstrom));
130 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Max internuclear distance: %12.5f a.u. = %12.5f Angstrom", (double)maxDist, (double)(maxDist/UNIT_one_Angstrom));
131
132 S_symm.resetSizesAndBlocks(matOpts.size_block_info,
133 matOpts.size_block_info);
134 if(compute_overlap_matrix_sparse(basisInfo, S_symm,
135 matOpts.permutationHML) != 0)
136 {
137 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error in compute_overlap_matrix_sparse");
138 throw "error in compute_overlap_matrix_sparse";
139 }
140
141 output_sparsity_symm(n, S_symm, "S_symm before trunc");
142 output_current_memory_usage(LOG_AREA_SCF, "after getting overlap matrix");
143 {
144 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
145 "truncating S using threshold value %6.2g", (double)scfopts.sparse_threshold_for_S);
146 double nnz_before_trunc_pc = (double)S_symm.nnz() * 100 / ((double)n*n);
147 ergo_real truncError = S_symm.eucl_thresh(scfopts.sparse_threshold_for_S);
148 double nnz_after_trunc_pc = (double)S_symm.nnz() * 100 / ((double)n*n);
149 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
150 "Truncated S (eucl), selected threshold = %10.6g, returned error = %10.6g, nnz before = %3.4f %%, nnz after = %3.4f %%",
151 scfopts.sparse_threshold_for_S, (double)truncError, nnz_before_trunc_pc, nnz_after_trunc_pc);
152 }
153
154
155 if ( scfopts.create_mtx_file_S == 1 ) {
156 // Write overlap matrix in matrix market format
157 std::stringstream ss_id;
158 ss_id << scfopts.calculation_identifier << " - overlap matrix";
159 write_matrix_in_matrix_market_format( S_symm, matOpts.inversePermutationHML, "S_matrix",
160 ss_id.str(), scfopts.method_and_basis_set );
161 }
162 if ( scfopts.create_basis_func_coord_file == 1 ) {
163 write_basis_func_coord_file(basisInfo);
164 }
165 if( scfopts.create_mtx_files_S_and_quit == 1 ) {
166 // Write overlap matrix in matrix market format using different basis function orderings, and then quit.
167 create_mtx_files_with_different_orderings(S_symm,
168 scfopts.calculation_identifier,
169 scfopts.method_and_basis_set,
170 matOpts.inversePermutationHML,
171 basisInfo);
172 throw "Breaking because create_mtx_files_S_and_quit was set.";
173 }
174 if ( scfopts.create_2el_integral_m_file == 1 )
175 write_2el_integral_m_file(basisInfo, integralInfo);
176
177 {
178 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Computing Euclidean norm of overlap matrix (just for fun)...");
179 Util::TimeMeter timeMeterEuclS;
180 ergo_real acc = -1; // value set by call below
181 ergo_real S_symm_euclnorm = get_eucl_norm_try_different_acc(S_symm, acc);
182 timeMeterEuclS.print(LOG_AREA_SCF, "S_symm.eucl()");
183 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Euclidean norm of overlap matrix = %22.11f (acc %7.3g)",
184 (double)S_symm_euclnorm, (double)acc);
185 }
186
187 invCholFactor.resetSizesAndBlocks(matOpts.size_block_info,
188 matOpts.size_block_info);
189 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Calling invCholFactor.inch with matOpts.threshold_inch = %g",
190 (double)matOpts.threshold_inch);
191 Util::TimeMeter timeMeterInch;
192 invCholFactor.inch(S_symm, matOpts.threshold_inch, mat::right);
193 timeMeterInch.print(LOG_AREA_SCF, "invCholFactor.inch");
194
195 output_sparsity_triang(n, invCholFactor, "invCholFactor before truncation");
196 {
197 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
198 "truncating Z using threshold value %6.2g", (double)scfopts.sparse_threshold_for_Z);
199 double nnz_before_trunc_pc = (double)invCholFactor.nnz() * 100 / ((double)n*n);
200 ergo_real truncError = invCholFactor.eucl_thresh(scfopts.sparse_threshold_for_Z);
201 double nnz_after_trunc_pc = (double)invCholFactor.nnz() * 100 / ((double)n*n);
202 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
203 "Truncated Z (eucl), selected threshold = %10.6g, returned error = %10.6g, nnz before = %3.4f %%, nnz after = %3.4f %%",
204 scfopts.sparse_threshold_for_Z, (double)truncError, nnz_before_trunc_pc, nnz_after_trunc_pc);
205 }
206 output_sparsity_triang(n, invCholFactor, "invCholFactor after truncation");
207 {
208 ergo_real invCholFactor_frobnorm = invCholFactor.frob();
209 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Frobenius norm of invCholFactor after truncation = %22.11f",
210 (double)invCholFactor_frobnorm);
211 Util::TimeMeter timeMeterEucl;
212 invCholFactor_euclnorm = invCholFactor.eucl( template_blas_sqrt(mat::getMachineEpsilon<ergo_real>()) );
213 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Euclidean norm of invCholFactor after truncation = %22.11f",
214 (double)invCholFactor_euclnorm);
215 timeMeterEucl.print(LOG_AREA_SCF, "invCholFactor.eucl()");
216 }
217 output_current_memory_usage(LOG_AREA_SCF, "after getting invCholFactor");
218
219 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "testing invCholFactor by computing ZT*S*Z");
220 {
221 Util::TimeMeter timeMeterZZTS;
222 normalMatrix noOver(S_symm);
223 normalMatrix inchcopy(invCholFactor);
224 normalMatrix tmpSZ, ID;
225 normalMatrix the_real_identity;
226 tmpSZ.resetSizesAndBlocks(matOpts.size_block_info,
227 matOpts.size_block_info);
228 ID.resetSizesAndBlocks(matOpts.size_block_info,
229 matOpts.size_block_info);
230 the_real_identity.resetSizesAndBlocks(matOpts.size_block_info,
231 matOpts.size_block_info);
232 tmpSZ = noOver * inchcopy;
233 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
234 "Truncating tmp matrix S*Z using eucl_thresh() with threshold value %6.3g",
235 (double)matOpts.sparse_threshold);
236 tmpSZ.eucl_thresh(matOpts.sparse_threshold);
237 ID = transpose(inchcopy) * tmpSZ;
238 the_real_identity = 1;
239 ergo_real frobenius_error = normalMatrix::frob_diff(ID, the_real_identity);
240 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "frobenius_error for ZT*S*Z is %10g", (double)frobenius_error);
241 timeMeterZZTS.print(LOG_AREA_SCF, "testing invCholFactor by computing ZT*S*Z");
242 }
243
244
245 S_symm.writeToFile();
246 if(scfopts.write_overlap_matrix) {
247 if(save_symmetric_matrix(S_symm, basisInfo, "overlap.bin",
248 matOpts.inversePermutationHML) != 0) {
249 do_output(LOG_CAT_ERROR, LOG_AREA_SCF,
250 "error in ddf_writeShellListAndDensityMatricesToFile");
251 throw "error in ddf_writeShellListAndDensityMatricesToFile";
252 }
253 }
254
255 invCholFactor.writeToFile();
256 output_current_memory_usage(LOG_AREA_SCF, "after writing invCholFactor and S_symm to file");
257
258 H_core_Matrix.resetSizesAndBlocks(matOpts.size_block_info, matOpts.size_block_info);
259 ergo_real nuclearRepulsionEnergyTmp = 0;
260 if(scfopts.skip_H_core == 1)
261 {
262 do_output(LOG_CAT_WARNING, LOG_AREA_SCF, "NOTE: skip_H_core parameter set, will skip construction of H_core matrix!");
263 do_output(LOG_CAT_WARNING, LOG_AREA_SCF, "NOTE: skip_H_core parameter set, results will be bogus!");
264 }
265 else
266 {
267 if(scfopts.use_simple_dense_H_core == 1) {
268 if(extraCharges.getNoOfAtoms() != 0) {
269 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error: (extraCharges.noOfAtoms != 0) not implemented for use_simple_dense_H_core case.");
270 throw "error: (extraCharges.noOfAtoms != 0) not implemented for use_simple_dense_H_core case.";
271 }
272 if(scfopts.electric_field.v[0] != 0 || scfopts.electric_field.v[1] != 0 || scfopts.electric_field.v[2] != 0) {
273 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error: electric field != 0 not implemented for use_simple_dense_H_core case.");
274 throw "error: electric field != 0 not implemented for use_simple_dense_H_core case.";
275 }
276 if(compute_h_core_matrix_simple_dense(integralInfo,
277 molecule,
278 basisInfo,
279 H_core_Matrix,
280 threshold_integrals_1el,
281 scfopts.no_of_threads_for_V,
282 matOpts.size_block_info,
283 matOpts.permutationHML,
284 nuclearRepulsionEnergyTmp) != 0)
285 {
286 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error in compute_h_core_matrix_simple_dense");
287 throw "error in compute_h_core_matrix_simple_dense";
288 }
289 }
290 else {
291 if(compute_h_core_matrix_sparse(integralInfo,
292 molecule,
293 extraCharges,
294 scfopts.electric_field.v[0],
295 scfopts.electric_field.v[1],
296 scfopts.electric_field.v[2],
297 basisInfo,
298 H_core_Matrix,
299 threshold_integrals_1el,
300 scfopts.no_of_threads_for_V,
301 scfopts.box_size_for_V_and_T,
302 nuclearRepulsionEnergyTmp,
303 matOpts.size_block_info,
304 matOpts.permutationHML,
305 scfopts.create_mtx_files_dipole,
306 &matOpts.inversePermutationHML,
307 &scfopts.calculation_identifier,
308 &scfopts.method_and_basis_set) != 0)
309 {
310 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error in compute_h_core_matrix_sparse");
311 throw "error in compute_h_core_matrix_sparse";
312 }
313 }
314 }
315 output_sparsity_symm(n, H_core_Matrix, "H_core_Matrix before trunc");
316 {
317 ergo_real subspaceThr = 0.001 * scfopts.purification_subspace_err_limit;
318 ergo_real threshold_Hcore = subspaceThr * scfopts.gap_expected_lower_bound / (1+subspaceThr);
319 double nnz_before_trunc_pc = (double)H_core_Matrix.nnz() * 100 / ((double)n*n);
320 invCholFactor.readFromFile();
321 ergo_real truncError = H_core_Matrix.eucl_thresh( threshold_Hcore, &invCholFactor );
322 invCholFactor.writeToFile();
323 double nnz_after_trunc_pc = (double)H_core_Matrix.nnz() * 100 / ((double)n*n);
324 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
325 "Truncated H_core_Matrix (eucl with Z), selected threshold = %10.6g, returned error = %10.6g, "
326 "nnz before = %3.4f %%, nnz after = %3.4f %%",
327 threshold_Hcore, truncError, nnz_before_trunc_pc, nnz_after_trunc_pc);
328 }
329 output_sparsity_symm(n, H_core_Matrix, "H_core_Matrix after trunc");
330
331 output_current_memory_usage(LOG_AREA_SCF, "after getting H_core_Matrix");
332
333 ergo_real nuclearElectricFieldEnergyTmp = molecule.getNuclearElectricFieldEnergy(scfopts.electric_field);
334 nuclearEnergy = nuclearRepulsionEnergyTmp + nuclearElectricFieldEnergyTmp;
335
336 if ( scfopts.create_mtx_file_H_core == 1 ) {
337 // Write H_core matrix in matrix market format
338 std::stringstream ss_id;
339 ss_id << scfopts.calculation_identifier << " - H_core matrix";
340 write_matrix_in_matrix_market_format( H_core_Matrix, matOpts.inversePermutationHML, "H_core_matrix",
341 ss_id.str(), scfopts.method_and_basis_set );
342 }
343
344 H_core_Matrix.writeToFile();
345 output_current_memory_usage(LOG_AREA_SCF, "after writing H_core_Matrix to file");
346
347 noOfElectrons = molecule.getNumberOfElectrons();
348
349 get_hf_weight_and_cam_params(scfopts.use_dft, &CAM_params.alpha,
350 &CAM_params.beta, &CAM_params.mu);
351
352 CAM_params.computeRangeSeparatedExchange =
353 CAM_params.beta != ergo_real(0.0);
354
355 energy_2el = 0;
356 energy = 0;
357 energy_2el_core = 0; // only used when "core density matrix" is used
358 energy_2el_valence = 0; // only used when "core density matrix" is used
359 energy_of_valence = 0; // only used when "core density matrix" is used
360 energy_reference = 0; // only used when "core density matrix" is used
361 electronicEntropyTerm = 0;
362
363
364 /*
365 Set parameters in the GetDensFromFock class.
366 All parameters must be set up, otherwise you will get exception.
367 */
368
369 // common parameters
370 DensFromFock.set_general_params(n,
371 matOpts.size_block_info);
372
373 #ifdef USE_CHUNKS_AND_TASKS
374 DensFromFock.set_cht_matrix_params(scfopts.cht_leavesSizeMax, scfopts.cht_blocksize);
375 #endif
376
377 DensFromFock.set_truncationNormPurification(scfopts.purification_truncation_norm);
378 DensFromFock.set_stopCriterionNormPurification(scfopts.purification_stop_crit_norm);
379 DensFromFock.set_invCholFactor(invCholFactor, invCholFactor_euclnorm);
380 DensFromFock.set_gap_expected_lower_bound(scfopts.gap_expected_lower_bound);
381 DensFromFock.set_purification_maxmul(scfopts.purification_maxmul);
382
383
384 if( scfopts.purification_create_m_files > 0 )
385 DensFromFock.set_purification_create_m_files();
386 else
387 DensFromFock.unset_purification_create_m_files();
388
389
390
391 if( scfopts.output_homo_and_lumo_eigenvectors > 0 )
392 DensFromFock.set_output_homo_and_lumo_eigenvectors();
393 else
394 DensFromFock.unset_output_homo_and_lumo_eigenvectors();
395
396
397 DensFromFock.set_number_of_eigenvectors_to_compute(
398 scfopts.number_of_occupied_eigenvectors, scfopts.number_of_unoccupied_eigenvectors);
399
400 DensFromFock.set_projection_method_params(
401 scfopts.go_back_X_iter_proj_method,
402 scfopts.jump_over_X_iter_proj_method);
403
404 /*
405 do not plot figures from the purification by default
406 set flag before calling function for purification
407 unset after calling function for purification
408 otherwise, it will be output from all its calls
409 */
410 DensFromFock.unset_generate_figures();
411
412
413 if( scfopts.purification_ignore_failure > 0 )
414 DensFromFock.set_purification_ignore_failure();
415 else
416 DensFromFock.unset_purification_ignore_failure();
417
418
419
420 if( scfopts.purification_use_rand_perturbation_for_alleigsint > 0 )
421 DensFromFock.set_purification_use_rand_perturbation_for_alleigsint();
422 else
423 DensFromFock.unset_purification_use_rand_perturbation_for_alleigsint();
424
425
426 // purification or something else?
427 if(scfopts.use_diagonalization == 0) // then use purification
428 DensFromFock.set_use_purification();
429 else
430 DensFromFock.unset_use_purification();
431
432 // purification settings
433 if(DensFromFock.get_use_purification() > 0)
434 {
435 DensFromFock.set_purification_limits(scfopts.purification_subspace_err_limit,
436 scfopts.purification_eigvalue_err_limit,
437 scfopts.puri_eig_acc_factor_for_guess);
438
439 if(scfopts.purification_with_acceleration > 0)
440 DensFromFock.set_use_acceleration();
441 else
442 DensFromFock.unset_use_acceleration();
443
444 if(scfopts.use_new_stopping_criterion > 0)
445 DensFromFock.set_use_new_stopping_criterion();
446 else
447 DensFromFock.unset_use_new_stopping_criterion();
448 }
449
450
451 if(scfopts.use_diag_on_error > 0)
452 DensFromFock.set_use_diag_on_error();
453 else
454 DensFromFock.unset_use_diag_on_error();
455
456 if(scfopts.use_diag_on_error_guess > 0)
457 DensFromFock.set_use_diag_on_error_guess();
458 else
459 DensFromFock.unset_use_diag_on_error_guess();
460
461
462
463
464 // diagonalization settings
465 if(scfopts.use_diagonalization > 0)
466 DensFromFock.set_use_diagonalization();
467 else
468 DensFromFock.unset_use_diagonalization();
469
470 // may be needed if scfopts.use_diag_on_error is set
471 DensFromFock.set_diagonalization_params(scfopts.electronic_temperature,
472 S_symm);
473
474 if(scfopts.use_diagonalization > 0 && scfopts.store_all_eigenvalues_to_file > 0)
475 DensFromFock.set_store_all_eigenvalues_to_file();
476 else
477 DensFromFock.unset_store_all_eigenvalues_to_file();
478
479
480 if(scfopts.save_permuted_F_matrix_in_bin > 0)
481 DensFromFock.set_save_permuted_F_matrix_in_bin();
482 else
483 DensFromFock.unset_save_permuted_F_matrix_in_bin();
484 }
485
486
~SCF_general()487 SCF_general::~SCF_general()
488 {
489 delete curr_cycle_stats;
490 }
491
492
GetEuclideanNormOfMatrix(const symmMatrix & A)493 ergo_real SCF_general::GetEuclideanNormOfMatrix(const symmMatrix & A)
494 {
495 ergo_real acc = template_blas_sqrt(mat::getMachineEpsilon<ergo_real>());
496 return A.eucl(acc);
497 }
498
499
get_overlap_matrix(symmMatrix & S)500 void SCF_general::get_overlap_matrix(symmMatrix & S)
501 {
502 S_symm.readFromFile();
503 S = S_symm;
504 S_symm.writeToFile();
505 }
506
507
get_invCholFactor_matrix(triangMatrix & invCholFactor_)508 void SCF_general::get_invCholFactor_matrix(triangMatrix & invCholFactor_)
509 {
510 invCholFactor.readFromFile();
511 invCholFactor_ = invCholFactor;
512 invCholFactor.writeToFile();
513 }
514
515
get_H_core_matrix(symmMatrix & H_core)516 void SCF_general::get_H_core_matrix(symmMatrix & H_core)
517 {
518 H_core_Matrix.readFromFile();
519 H_core = H_core_Matrix;
520 H_core_Matrix.writeToFile();
521 }
522
523
get_energy(ergo_real & E,ergo_real & E_nuclear)524 void SCF_general::get_energy(ergo_real & E, ergo_real & E_nuclear)
525 {
526 E = energy;
527 E_nuclear = nuclearEnergy;
528 }
529
530
do_SCF_iterations()531 void SCF_general::do_SCF_iterations()
532 {
533 Util::TimeMeter timeMeterTot;
534 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_general::do_SCF_iterations");
535
536 // Initialize DIIS
537 if(DIIS == NULL)
538 throw "ERROR: (DIIS == NULL)";
539 if(DIIS->Initialize(scfopts.max_no_of_diis_matrices) != 0)
540 throw "Error in DIIS->Initialize";
541
542 initialize_matrices();
543
544 if(J_K_params.threshold_J <= 0 || J_K_params.threshold_K <= 0)
545 throw "Error in SCF_general::do_SCF_iterations: (J_K_params.threshold_J <= 0 || J_K_params.threshold_K <= 0).";
546
547 // Check that parameters are reasonable, even number of electrons for restricted etc.
548 check_params();
549
550 // set up starting guess
551 get_starting_guess_density();
552
553 if(scfopts.spin_flip_atom_count > 0)
554 do_spin_flip(scfopts.spin_flip_atom_count);
555
556 if(scfopts.starting_guess_disturbance > 0)
557 add_random_disturbance_to_starting_guess();
558
559 if(scfopts.write_guess_density_only == 1) {
560 write_density_to_file();
561 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "throwing exception after write_density_to_file.");
562 throw "exiting after scfopts.write_guess_density_only";
563 }
564
565 if(scfopts.output_density_images_only == 1) {
566 output_density_images();
567 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "throwing exception after output_density_images.");
568 throw "exiting after scfopts.output_density_images_only";
569 }
570
571
572 // Use at least two iterations if disimple_starting_guess used.
573 int min_number_of_iterations = scfopts.min_number_of_iterations;
574 if(scfopts.use_simple_starting_guess && min_number_of_iterations < 2)
575 min_number_of_iterations = 2;
576
577 ergo_real stepLength = scfopts.step_length_start;
578 ergo_real best_energy_so_far = 0;
579 int step = 0;
580 SCF_step = 0;
581 int noOfFailuresInARow = 0;
582 int restartCount = 0;
583 int noOfImprovementsInARow = 0;
584 int diisUsedInLastIteration = 0;
585
586 curr_subspace_diff = 0;
587
588 initialize_homo_lumo_limits();
589
590 // The different Fockmatrix objects are empty at this point,
591 // But we write them to file because they are supposed to be on file
592 // in the beginning of each SCF cycle.
593 write_matrices_to_file();
594
595
596 output_current_memory_usage(LOG_AREA_SCF, "before main SCF loop");
597
598 Util::TimeMeter timeMeterScfMainLoop;
599
600 // main SCF loop
601 while(1)
602 {
603 curr_cycle_stats = new SCF_statistics;
604 curr_cycle_stats->start_timer("scf_cycle_time");
605 curr_cycle_stats->add_value("no_of_basis_func", basisInfo.noOfBasisFuncs);
606 curr_cycle_stats->add_value("sparse_matrix_block_size", matOpts.sparse_matrix_block_size);
607 step++;
608 SCF_step++;
609 curr_cycle_stats->add_value("scf_cycle", step);
610 char infoString[888];
611 sprintf(infoString, "Beginning of SCF cycle %i: ", step);
612 do_output_time(LOG_CAT_INFO, LOG_AREA_SCF, infoString);
613
614 output_current_memory_usage(LOG_AREA_SCF, infoString);
615
616 Util::TimeMeter timeMeterStep;
617
618 save_density_as_prevdens();
619
620 get_2e_part_and_energy();
621
622 //if(scfopts.use_artificial_subspace_disturbances == 1)
623 //disturb_fock_matrix(curr_subspace_diff * scfopts.subspace_factor_fock);
624
625 // Now we have created Fock matrix (or matrices), and written to file. Memory usage should be the same as before.
626 output_current_memory_usage(LOG_AREA_SCF, "After get_2e_part_and_energy");
627
628 output_sparsity_S_F_D(*curr_cycle_stats);
629
630 double virtMem = 0, resMem = 0, virtPeakMem = 0;
631 get_memory_usage_by_procfile(&virtMem, &resMem, &virtPeakMem);
632 curr_cycle_stats->add_value("peak_virt_mem_usage_GB", virtPeakMem);
633
634
635 calculate_energy();
636
637
638 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "nuclearEnergy, energy_2el, energy = %f, %f, %f",
639 (double)nuclearEnergy, (double)energy_2el, (double)energy);
640 // ELIAS NOTE 2014-01-01: FIXME: consider adding another output message here showing the "electronic energy" so that the difference between "electronic energy" and "total energy" becomes more clear.
641 if( step > 2)
642 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Energy %3i = %22.11f ( diff %22.15f )", step, (double)energy, (double)(energy-best_energy_so_far));
643 else
644 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Energy %3i = %22.11f", step, (double)energy);
645
646 if(scfopts.compute_core_density == 1 && step > 1) {
647 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Valence-only energy %3i = %22.11f", step, (double)energy_of_valence);
648 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Reference energy %3i = %22.11f", step, (double)energy_reference);
649 }
650
651 output_current_memory_usage(LOG_AREA_SCF, "After computing energy");
652
653 curr_cycle_stats->start_timer("get_FDSminusSDF");
654 Util::TimeMeter timeMeterFDSminusSDF;
655 get_FDSminusSDF();
656 timeMeterFDSminusSDF.print(LOG_AREA_SCF, "get_FDSminusSDF");
657 output_current_memory_usage(LOG_AREA_SCF, "After computing FDS-SDF");
658 curr_cycle_stats->stop_timer("get_FDSminusSDF");
659
660 Util::TimeMeter timeMeterGerErrorMeasure;
661 get_error_measure();
662 timeMeterGerErrorMeasure.print(LOG_AREA_SCF, "get_error_measure");
663 output_current_memory_usage(LOG_AREA_SCF, "After computing error measure");
664
665 // Check if converged
666 if(scfopts.use_artificial_subspace_disturbances == 1 && step > 1 && curr_subspace_diff < 1e-6)
667 {
668 do_output(LOG_CAT_RESULTS, LOG_AREA_SCF, "CONVERGED due to curr_subspace_diff after %3i iterations.", step);
669 do_output(LOG_CAT_RESULTS, LOG_AREA_SCF, "FINAL ENERGY: %22.11f", (double)energy);
670 report_final_results();
671 break;
672 }
673 if(errorMeasure < scfopts.convergence_threshold && step >= min_number_of_iterations)
674 {
675 // Write Fock matrix in matrix market format in the last SCF cycle
676 if ( scfopts.create_mtx_files_F == 1 )
677 create_mtx_files_F( step );
678
679 do_output(LOG_CAT_RESULTS, LOG_AREA_SCF, "CONVERGED after %3i iterations.", step);
680 do_output(LOG_CAT_RESULTS, LOG_AREA_SCF, "FINAL ENERGY: %22.11f", (double)energy);
681 report_final_results();
682 break;
683 }
684
685 // Check if max number of iterations reached
686 if(scfopts.max_number_of_iterations > 0 && step >= scfopts.max_number_of_iterations)
687 {
688 // Write Fock matrix in matrix market format in the last SCF cycle
689 if ( scfopts.create_mtx_files_F == 1 )
690 create_mtx_files_F( step );
691
692 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Maximum number of SCF iterations reached. Breaking SCF procedure.");
693 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Gave up after %i iterations.", step);
694 break;
695 }
696
697
698
699 if(step == 1)
700 {
701 // first time
702 add_to_DIIS_list();
703 }
704 else
705 {
706 // not first time
707 if((energy < best_energy_so_far) &&
708 (step > noOfImprovementsInARow+3) &&
709 (noOfImprovementsInARow < scfopts.no_of_impr_req_for_diis ||
710 errorMeasure > scfopts.error_maxabs_for_diis) &&
711 !(diisUsedInLastIteration == 1 && noOfImprovementsInARow > 1) &&
712 !(scfopts.use_diis_always == 1))
713 {
714 noOfImprovementsInARow++;
715 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "the energy got better, "
716 "but we do not dare to try DIIS yet, noOfImprovementsInARow = %i",
717 noOfImprovementsInARow);
718 best_energy_so_far = energy;
719 noOfFailuresInARow = 0;
720
721 update_best_fock_so_far();
722
723 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
724 "mixing best Fock matrix so far with next one, stepLength = %10.6f",
725 (double)stepLength);
726 combine_old_fock_matrices(stepLength);
727
728 diisUsedInLastIteration = 0;
729 stepLength *= 1.1;
730 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
731 "increased steplength by factor 1.1, stepLength = %10.6f", (double)stepLength);
732 }
733 else if(step <= scfopts.no_of_careful_first_scf_steps)
734 {
735 // "careful" option chosen: in this case we do not use
736 // DIIS for the early steps. Can be useful when the
737 // starting guess is very bad.
738 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Not considering DIIS now because 'careful' option "
739 "chosen, no_of_careful_first_scf_steps = %2d", scfopts.no_of_careful_first_scf_steps);
740 /* Elias note 2010-05-12: added || step == 2 in this if
741 statement condition to handle the case when
742 no_of_careful_first_scf_steps is used and the energy
743 increased in step 2, which heppened for the Umeda
744 protein molecule. */
745 if(energy < best_energy_so_far || step == 2) {
746 best_energy_so_far = energy;
747 update_best_fock_so_far();
748 stepLength *= 1.1;
749 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
750 "increased steplength by factor 1.1, stepLength = %10.6f", (double)stepLength);
751 }
752 else {
753 // Energy got worse.
754 noOfImprovementsInARow = 0;
755 noOfFailuresInARow++;
756 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "restarting DIIS because energy did not improve.");
757 clear_diis_list();
758 ergo_real newStepLength = stepLength * 0.5;
759 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "changing stepLength from %10.6f to %10.6f",
760 (double)stepLength, (double)newStepLength);
761 stepLength = newStepLength;
762 }
763 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
764 "mixing best Fock matrix so far with next one, stepLength = %10.6f",
765 (double)stepLength);
766 combine_old_fock_matrices(stepLength);
767 diisUsedInLastIteration = 0;
768 }
769 else if(step == 2 || energy < best_energy_so_far || scfopts.use_diis_always == 1)
770 {
771 if(step > 2 && energy < best_energy_so_far)
772 noOfImprovementsInARow++;
773 // energy got better
774 noOfFailuresInARow = 0;
775 best_energy_so_far = energy;
776
777 update_best_fock_so_far();
778
779 diisUsedInLastIteration = 1;
780
781 add_to_DIIS_list();
782
783 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "using DIIS to get combined Fock matrix, "
784 "number of iters used for DIIS: %2i", DIIS->GetNoOfIters());
785
786 use_diis_to_get_new_fock_matrix();
787 }
788 else
789 {
790 // energy got worse
791 if(scfopts.break_on_energy_increase == 1)
792 {
793 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Energy increased. Breaking SCF.");
794 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Stopped SCF after %i iterations.", step);
795 break;
796 }
797 if(stepLength < scfopts.step_length_giveup)
798 {
799 if(restartCount > scfopts.max_restart_count)
800 {
801 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
802 "The energy does not seem to get any lower. We give up!");
803 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Gave up after %i iterations.", step);
804 break;
805 }
806 restartCount++;
807 stepLength = scfopts.step_length_start;
808 noOfFailuresInARow = 0;
809 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "reset stepLength to %10.6f, restartCount = %i",
810 (double)scfopts.step_length_start, restartCount);
811 }
812 noOfImprovementsInARow = 0;
813 noOfFailuresInARow++;
814 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "restarting DIIS because energy did not improve.");
815 clear_diis_list();
816
817 if(diisUsedInLastIteration == 0)
818 {
819 ergo_real newStepLength = stepLength * 0.5;
820 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "changing stepLength from %10.6f to %10.6f",
821 (double)stepLength, (double)newStepLength);
822 stepLength = newStepLength;
823 }
824 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
825 "mixing best Fock matrix so far with next one, stepLength = %10.6f",
826 (double)stepLength);
827 combine_old_fock_matrices(stepLength);
828
829 diisUsedInLastIteration = 0;
830 }
831 }
832
833 output_current_memory_usage(LOG_AREA_SCF, "After creating lin comb F");
834
835 // Free memory used by Err_sparse
836 Util::TimeMeter timeMeterClearErrorMatrices;
837 clear_error_matrices();
838 timeMeterClearErrorMatrices.print(LOG_AREA_SCF, "clear_error_matrices");
839 output_current_memory_usage(LOG_AREA_SCF, "After clear_error_matrices");
840
841 // FIXME: Compare F and Fprev here to get info about gap of
842 // F? Such info will be needed as input to get_dens_from_fock?
843
844 Util::TimeMeter timeMeterSaveFockAsFprev;
845 save_current_fock_as_fprev();
846 timeMeterSaveFockAsFprev.print(LOG_AREA_SCF, "save_current_fock_as_fprev");
847
848 curr_cycle_stats->start_timer("get_new_density_matrix");
849 Util::TimeMeter timeMeterGetNewDensityMatrix;
850 get_new_density_matrix();
851 timeMeterGetNewDensityMatrix.print(LOG_AREA_SCF, "get_new_density_matrix");
852 curr_cycle_stats->stop_timer("get_new_density_matrix");
853
854 // At this point a new density matrix has just been computed, so the electronic entropy term has also been computed in the nonzero-temperature case.
855 if(scfopts.electronic_temperature > 0) {
856 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Electronic temperature: %22.11f a.u. = %22.2f Kelvin",
857 (double)scfopts.electronic_temperature, (double)scfopts.electronic_temperature/UNIT_one_Kelvin);
858 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Electronic entropy term: %22.11f, energy with entropy term included: %22.11f",
859 (double)electronicEntropyTerm, (double)(energy + electronicEntropyTerm));
860 }
861
862 if(step > 1)
863 report_density_difference();
864
865 if(scfopts.use_artificial_subspace_disturbances == 1 && step > 2)
866 disturb_dens_matrix_exact(curr_subspace_diff * scfopts.subspace_factor_dens);
867 //disturb_dens_matrix(curr_subspace_diff * scfopts.subspace_factor_dens);
868
869 if(scfopts.use_artificial_subspace_disturbances == 1)
870 update_subspace_diff();
871
872 if(scfopts.output_density_at_every_step == 1) {
873 Util::TimeMeter timeMeterWriteDensityToFile;
874 output_current_memory_usage(LOG_AREA_SCF, "before write_density_to_file()");
875 write_density_to_file();
876 output_current_memory_usage(LOG_AREA_SCF, "after write_density_to_file()");
877 timeMeterWriteDensityToFile.print(LOG_AREA_SCF, "write_density_to_file");
878 }
879
880 if ( scfopts.create_mtx_files_F == 1 )
881 // Write Fock matrix in matrix market format
882 create_mtx_files_F( step );
883 if ( scfopts.create_mtx_files_D == 1 )
884 // Write Fock matrix in matrix market format
885 create_mtx_files_D( step );
886
887 if ( scfopts.output_homo_and_lumo_eigenvectors ) {
888 // Write homo and lumo eigenvectors to file
889 create_eigenvectors_files();
890 create_gabedit_file();
891 }
892
893 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF cycle %3i finished.", step);
894 char tempString[88];
895 sprintf(tempString, "SCF cycle %3i", step);
896 timeMeterStep.print(LOG_AREA_SCF, tempString);
897 curr_cycle_stats->stop_timer("scf_cycle_time");
898 std::stringstream ss;
899 ss << "scf_cycle_" << step;
900 if(scfopts.output_statistics_mfiles)
901 curr_cycle_stats->output_mfile( ss.str() );
902 delete curr_cycle_stats;
903 curr_cycle_stats = NULL; // since curr_cycle_stats is deleted in destructor, we need to set to null here to avoid double-delete if we exit loop in some unusual way.
904 } // END WHILE main SCF loop
905
906 timeMeterScfMainLoop.print(LOG_AREA_SCF, "Main SCF loop");
907 output_current_memory_usage(LOG_AREA_SCF, "after main SCF loop");
908 std::string allocStatsStr = mat::AllocatorManager<ergo_real>::instance().getStatistics();
909 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "After main SCF loop: %s", allocStatsStr.c_str());
910
911 if(scfopts.output_density_at_every_step == 1 && step == 1) {
912 Util::TimeMeter timeMeterWriteDensityToFile;
913 output_current_memory_usage(LOG_AREA_SCF, "before write_density_to_file()");
914 write_density_to_file();
915 output_current_memory_usage(LOG_AREA_SCF, "after write_density_to_file()");
916 timeMeterWriteDensityToFile.print(LOG_AREA_SCF, "write_density_to_file");
917 }
918
919 compute_dipole_moment();
920
921 if(scfopts.output_mulliken_pop == 1)
922 do_mulliken_pop_stuff();
923
924 if(scfopts.compute_gradient_fixeddens == 1)
925 compute_gradient_fixeddens();
926
927 if(scfopts.save_full_matrices_for_matlab == 1)
928 save_full_matrices_for_matlab();
929
930 if(scfopts.save_final_potential == 1)
931 save_final_potential();
932
933 if(scfopts.output_expected_values_pos_operator == 1)
934 output_expected_values_pos_operator();
935
936 if(scfopts.output_density_images == 1)
937 output_density_images();
938
939 if(scfopts.write_diag_dens_to_file == 1)
940 write_diag_dens_to_file();
941
942 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_general::do_SCF_iterations finished.");
943 timeMeterTot.print(LOG_AREA_SCF, "SCF_general::do_SCF_iterations");
944 }
945