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 GetDensFromFock.h
31  *
32  *  @brief Routines for getting density matrix from a given Fock
33  *         matrix.
34  *
35  *  @author Anastasia Kruchinina <em>responsible</em>
36  *  @author Elias Rudberg
37  */
38 
39 #ifndef GETDENSFROMFOCKHEADER
40 #define GETDENSFROMFOCKHEADER
41 
42 #include "realtype.h"
43 #include "matrix_typedefs.h"
44 #include "matrix_typedefs_chtml.h"
45 #include "transform.h"
46 #include "output.h"
47 
48 
49 
50 /** GetDensFromFock class containing parameters and functions for computing density matrix.
51  *
52  * Flags are set to undefined value by default.  User should define
53  * them explicitly, otherwise exception is thrown if undefined flag is
54  * used.
55  */
56 class GetDensFromFock
57 {
58 public:
59 
60    static const int UNDEF_VALUE_UINT;
61    static const std::string NA_STRING;
62 
63    static const bool BOOL_TRUE;
64    static const bool BOOL_FALSE;
65 
66    void create_checkpoint(symmMatrix&   Finput,       /**< [in] Effective Hamiltonian matrix (written to file) */
67                           symmMatrix&   F_ort_prev,   /**< [in/out]
68                                                        *    Input: Previous F matrix in orthogonal basis. (written to file)
69                                                        *    Output: New F matrix in orthogonal basis ( ZT*Finput*Z ). (written to file) */
70                           generalVector *eigVecLUMO,  /**< [out] LUMO eigenvector */
71                           generalVector *eigVecHOMO,  /**< [out] HOMO eigenvector */
72                           std::string   IDstr         /**< [in] File identificator; added to the name of each file */
73                           );
74 
75 
76    static void restore_from_checkpoint(GetDensFromFock& DensFromFock,    /**< [out] Instance of GetDensFromFock class contatining all data for computing the density matrix */
77                                        symmMatrix&      Finput,          /**< [out] Effective Hamiltonian matrix (written to file) */
78                                        symmMatrix&      F_ort_prev,      /**< [out] F matrix in orthogonal basis ( ZT*Finput*Z ). (written to file) */
79                                        generalVector    *eigVecLUMO,     /**< [out] LUMO eigenvector */
80                                        generalVector    *eigVecHOMO,     /**< [out] HOMO eigenvector */
81                                        std::string      checkpoint_path, /**< [out] HOMO eigenvector */
82                                        std::string      IDstr,           /**< [in]  File identificator; added to the name of each file. */
83                                        int              SCF_step         /**< [in]  SCF step which should be restored; added to the name of each file in given SCF cycle. */
84                                        );
85 
86 
87    //constructor
GetDensFromFock()88    GetDensFromFock()
89    {
90       do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Create object from GetDensFromFock.");
91 
92       // set all variables and flags to default values
93       n = UNDEF_VALUE_UINT;
94       noOfOccupiedOrbs       = UNDEF_VALUE_UINT;
95       factor                 = UNDEF_VALUE_UINT;
96       invCholFactor_euclnorm = 0;
97       maxMul                 = 0;
98       plot_puri_results      = BOOL_FALSE;
99       SCF_step               = UNDEF_VALUE_UINT;
100 
101       use_diagonalization     = BOOL_FALSE;
102       use_purification        = BOOL_FALSE;
103 
104       electronicTemperature                = 0;
105       gap_expected_lower_bound             = 0;
106       eigvalueErrorLimit                   = 0;
107       subspaceErrorLimit                   = 0;
108       puri_eig_acc_factor_for_guess        = 0;
109       use_diag_on_error       = BOOL_FALSE;
110       use_diag_on_error_guess = BOOL_FALSE;
111 
112       create_m_files     = BOOL_FALSE;
113       output_homo_and_lumo_eigenvectors    = BOOL_FALSE;
114       number_of_occupied_eigenvectors = 0;
115       number_of_unoccupied_eigenvectors = 0;
116       go_back_X_iter_proj_method = 0;
117       jump_over_X_iter_proj_method = 0;
118       ignore_purification_failure          = BOOL_FALSE;
119       use_rand_perturbation_for_alleigsint = BOOL_FALSE;
120       use_acceleration                   = BOOL_FALSE;
121       use_new_stopping_criterion         = BOOL_FALSE;
122       store_all_eigenvalues_to_file      = BOOL_FALSE;
123       try_eigv_on_next_iteration_if_fail = BOOL_FALSE;
124 
125       leavesSizeMax         = UNDEF_VALUE_UINT;
126       blocksize             = UNDEF_VALUE_UINT;
127 
128       eigenvectors_method                 = NA_STRING;
129       eigenvectors_iterative_method       = NA_STRING;
130       use_prev_vector_as_initial_guess    = BOOL_FALSE;
131       puri_compute_eigv_in_each_iteration = BOOL_FALSE;
132       run_shift_and_square_method_on_F    = BOOL_FALSE;
133       save_permuted_F_matrix_in_bin       = BOOL_FALSE;
134 
135       eigensolver_accuracy = 0;
136       eigensolver_maxiter  = 0;
137 
138       std::string stats_prefix = ""; // default value
139 
140       truncationNormPurification = mat::euclNorm;  // default value
141       stopCriterionNormPurification = mat::euclNorm;  // default value
142 
143       clean_eigs_intervals();
144       clean_puri_stats();
145 
146 
147       filenameFinput        = "matrix_Finput";
148       filenameF_ort_prev    = "matrix_F_ort_prev";
149       filenameeigVecLUMO    = "vector_eigVecLUMO";
150       filenameeigVecHOMO    = "vector_eigVecHOMO";
151       filenameOverlap       = "matrix_Overlap";
152       filenameD_ort_prev    = "matrix_D_ort_prev";
153       filenameinvCholFactor = "matrix_invCholFactor";
154       file_for_basic_types  = "basic_types";
155    }
156 
157    /** Choose which method to use for computing the density matrix from Fock matrix.
158     */
159    int get_dens_from_fock(symmMatrix&   Finput,           /**< [in] Effective Hamiltonian matrix. (written to file) */
160                           symmMatrix&   resultDens,       /**< [out] Density matrix. (written to file) */
161                           symmMatrix&   F_ort_prev       /**< [in/out]
162                                                            *      Input: Previous F matrix in orthogonal basis. (written to file)
163                                                            *      Output: New F matrix in orthogonal basis ( ZT*Finput*Z ). (written to file) */
164                           );
165 
166 
167 
168    /** Use recursive expansion for computing the density matrix from Fock matrix.
169     */
170    int get_dens_from_fock_sparse(
171      symmMatrix&   F,                  /**< [in] Effective Hamiltonian matrix. (written to file) */
172      symmMatrix&   resultDens,         /**< [out] Density matrix. (written to file) */
173      symmMatrix&   F_ort_prev         /**< [in/out]
174      Input: Previous F matrix in orthogonal basis. (written to file)
175      Output: New F matrix in orthogonal basis ( ZT*Finput*Z ). (written to file)
176      */
177                                  );
178 
179 
get_computed_eigenpairs(std::vector<generalVector> & eigVecUNOCCref,std::vector<generalVector> & eigVecOCCref,std::vector<ergo_real> & eigValUNOCCref,std::vector<ergo_real> & eigValOCCref)180    void get_computed_eigenpairs(
181      std::vector<generalVector> &eigVecUNOCCref, /**< [out] Unoccupied eigenvectors */
182      std::vector<generalVector> &eigVecOCCref,  /**< [out] Occupied eigenvectors */
183      std::vector<ergo_real> &eigValUNOCCref,  /**< [out] Unoccupied eigenvalues */
184      std::vector<ergo_real> &eigValOCCref /**< [out] Occupied eigenvalues */
185    )
186    {
187      eigVecUNOCCref = eigVecUNOCC;
188      eigVecOCCref = eigVecOCC;
189      eigValOCCref = eigValOCC;
190      eigValUNOCCref = eigValUNOCC;
191    }
192 
193    /** Set bounds for HOMO and LUMO eigenvalues to -/+ inf, thus remove
194     *  any known bounds.
195     */
clean_eigs_intervals()196    inline void clean_eigs_intervals()
197    {
198       homoInterval_Finput = intervalType(-1e22, 1e22);
199       lumoInterval_Finput = intervalType(-1e22, 1e22);
200       homoInterval_Finput = intervalType(-1e22, 1e22);
201       lumoInterval_Finput = intervalType(-1e22, 1e22);
202 
203       homoInterval_F_ort_prev = intervalType(-1e22, 1e22);
204       lumoInterval_F_ort_prev = intervalType(-1e22, 1e22);
205       homoInterval_F_ort_prev = intervalType(-1e22, 1e22);
206       lumoInterval_F_ort_prev = intervalType(-1e22, 1e22);
207    }
208 
209 
210 
211 
set_SCF_step(int step)212    inline void set_SCF_step(int step /**< [in] Current SCF step */)
213    { SCF_step = step; }
unset_SCF_step()214    inline void unset_SCF_step()
215    { SCF_step = UNDEF_VALUE_UINT; }
216 
217 
218    /** Plot figures from the  recursive expansion.
219     */
220    inline void set_generate_figures(std::string str = "" /**< [in] String added to each generated file. */)
221    {
222       if (create_m_files)
223       {
224          assert(SCF_step >= 0);
225          plot_puri_results     = BOOL_TRUE;
226          plot_puri_results_str = str;
227       }
228    }
229 
230    /** Do not plot figures from the  recursive expansion.
231     */
unset_generate_figures()232    inline void unset_generate_figures()
233    {
234       plot_puri_results     = BOOL_FALSE;
235       plot_puri_results_str = "";
236    }
237 
set_general_params(const int n_,mat::SizesAndBlocks const & matrixSizesAndBlocks_)238    inline void set_general_params(const int n_,                   /**< [in] Number of basis functions. */
239                                   mat::SizesAndBlocks const& matrixSizesAndBlocks_ /**< [in] Matrix library parameters. */
240                                   )
241    {
242       assert(n_ >= 1);
243       n = n_;
244       matrixSizesAndBlocks = matrixSizesAndBlocks_;
245    }
246 
set_cht_matrix_params(const int leavesSizeMax_,const int blocksize_)247    inline void set_cht_matrix_params(const int leavesSizeMax_,  /**< [in] CHTMatrix library parameter leavesSizeMax. */
248                                      const int blocksize_       /**< [in] CHTMatrix library parameter blocksize. */
249                                      )
250    {
251       assert(leavesSizeMax_ >= 1);
252       assert(blocksize_ >= 1);
253       leavesSizeMax = leavesSizeMax_;
254       blocksize     = blocksize_;
255    }
256 
get_SizesAndBlocks(mat::SizesAndBlocks & matrixSizesAndBlocks_)257    inline void get_SizesAndBlocks(mat::SizesAndBlocks& matrixSizesAndBlocks_   /**< [out] Matrix library parameters. */
258                                   ) const
259    {
260       matrixSizesAndBlocks_ = matrixSizesAndBlocks;
261    }
262 
263    /** Set truncation norm used in the recursive expansion.
264     *  Possible norms: spectral, Frobenius or mixed.
265     */
set_truncationNormPurification(mat::normType const truncationNormPurification_)266    inline void set_truncationNormPurification(mat::normType const truncationNormPurification_ /**< [in]  Norm used in truncation. */)
267    { truncationNormPurification = truncationNormPurification_; }
268 
269    /** Set stopping criterion norm used in the recursive expansion.
270     *  Possible norms: spectral, Frobenius or mixed.
271     */
set_stopCriterionNormPurification(mat::normType const stopCriterionNormPurification_)272    inline void set_stopCriterionNormPurification(mat::normType const stopCriterionNormPurification_ /**< [in] Norm used in the stopping criterion. */)
273    { stopCriterionNormPurification = stopCriterionNormPurification_; }
274 
275 
do_restricted_calculations()276    inline void do_restricted_calculations()
277    { factor = 2; }
278 
do_unrestricted_calculations()279    inline void do_unrestricted_calculations()
280    { factor = 1; }
281 
set_no_occupied_orbs(int noOfOccupiedOrbs_)282    inline void set_no_occupied_orbs(int noOfOccupiedOrbs_)
283    {
284       assert(noOfOccupiedOrbs_ >= 0);
285       noOfOccupiedOrbs = noOfOccupiedOrbs_;
286    }
287 
clean_puri_stats()288    inline void clean_puri_stats()
289    { puri_stats.clear(); }
290 
291 
set_invCholFactor(triangMatrix const & invCholFactor_,ergo_real invCholFactor_euclnorm_)292    inline void set_invCholFactor(triangMatrix const& invCholFactor_,
293                                  ergo_real           invCholFactor_euclnorm_)
294    {
295       invCholFactor = invCholFactor_;
296       assert(invCholFactor_euclnorm_ >= 0);
297       invCholFactor_euclnorm = invCholFactor_euclnorm_;
298    }
299 
set_gap_expected_lower_bound(ergo_real gap_expected_lower_bound_)300    inline void set_gap_expected_lower_bound(ergo_real gap_expected_lower_bound_)
301    {
302       assert(gap_expected_lower_bound_ >= 0);
303       gap_expected_lower_bound = gap_expected_lower_bound_;
304    }
305 
306    /** Set maximum allowed number of iterations in recursive expansion.
307     */
set_purification_maxmul(ergo_real purification_maxmul_)308    inline void set_purification_maxmul(ergo_real purification_maxmul_)
309    {
310       assert(purification_maxmul_ > 0);
311       maxMul = purification_maxmul_;
312    }
313 
set_number_of_eigenvectors_to_compute(int occ,int unocc)314    inline void set_number_of_eigenvectors_to_compute(int occ, int unocc)
315    {
316      number_of_occupied_eigenvectors = occ;
317      number_of_unoccupied_eigenvectors = unocc;
318    }
319 
set_projection_method_params(int go_back,int step)320    inline void set_projection_method_params(int go_back, int step)
321    {
322      go_back_X_iter_proj_method = go_back;
323      jump_over_X_iter_proj_method = step;
324    }
325 
326    /****  SET/UNSET SECTION *****/
327 
get_purification_create_m_files()328    inline int get_purification_create_m_files() const
329    { return create_m_files == BOOL_TRUE; }
set_purification_create_m_files()330    inline void set_purification_create_m_files()
331    { create_m_files = BOOL_TRUE; }
unset_purification_create_m_files()332    inline void unset_purification_create_m_files()
333    { create_m_files = BOOL_FALSE; }
334 
335 
336 
get_output_homo_and_lumo_eigenvectors()337    inline int get_output_homo_and_lumo_eigenvectors() const
338    { return output_homo_and_lumo_eigenvectors == BOOL_TRUE; }
set_output_homo_and_lumo_eigenvectors()339    inline void set_output_homo_and_lumo_eigenvectors()
340    { output_homo_and_lumo_eigenvectors = BOOL_TRUE; }
unset_output_homo_and_lumo_eigenvectors()341    inline void unset_output_homo_and_lumo_eigenvectors()
342    { output_homo_and_lumo_eigenvectors = BOOL_FALSE; }
343 
344 
get_purification_ignore_failure()345    inline int get_purification_ignore_failure() const
346    { return ignore_purification_failure == BOOL_TRUE; }
set_purification_ignore_failure()347    inline void set_purification_ignore_failure()
348    { ignore_purification_failure = BOOL_TRUE; }
unset_purification_ignore_failure()349    inline void unset_purification_ignore_failure()
350    { ignore_purification_failure = BOOL_FALSE; }
351 
352 
get_use_rand_perturbation_for_alleigsint()353    inline int get_use_rand_perturbation_for_alleigsint() const
354    { return use_rand_perturbation_for_alleigsint == BOOL_TRUE; }
set_purification_use_rand_perturbation_for_alleigsint()355    inline void set_purification_use_rand_perturbation_for_alleigsint()
356    { use_rand_perturbation_for_alleigsint = BOOL_TRUE; }
unset_purification_use_rand_perturbation_for_alleigsint()357    inline void unset_purification_use_rand_perturbation_for_alleigsint()
358    { use_rand_perturbation_for_alleigsint = BOOL_FALSE; }
359 
get_use_diagonalization()360    inline int get_use_diagonalization() const
361    { return use_diagonalization == BOOL_TRUE; }
set_use_diagonalization()362    inline void set_use_diagonalization()
363    { use_diagonalization = BOOL_TRUE; }
unset_use_diagonalization()364    inline void unset_use_diagonalization()
365    { use_diagonalization = BOOL_FALSE; }
366 
367 
get_use_purification()368    inline int get_use_purification() const
369    { return use_purification == BOOL_TRUE; }
set_use_purification()370    inline void set_use_purification()
371    { use_purification = BOOL_TRUE; }
unset_use_purification()372    inline void unset_use_purification()
373    { use_purification = BOOL_FALSE; }
374 
375 
get_use_diag_on_error_guess()376    inline int get_use_diag_on_error_guess() const
377    { return use_diag_on_error_guess == BOOL_TRUE; }
set_use_diag_on_error_guess()378    inline void set_use_diag_on_error_guess()
379    { use_diag_on_error_guess = BOOL_TRUE; }
unset_use_diag_on_error_guess()380    inline void unset_use_diag_on_error_guess()
381    { use_diag_on_error_guess = BOOL_FALSE; }
382 
383 
get_use_diag_on_error()384    inline int get_use_diag_on_error() const
385    { return use_diag_on_error == BOOL_TRUE; }
set_use_diag_on_error()386    inline void set_use_diag_on_error()
387    { use_diag_on_error = BOOL_TRUE; }
unset_use_diag_on_error()388    inline void unset_use_diag_on_error()
389    { use_diag_on_error = BOOL_FALSE; }
390 
391 
get_stats_prefix()392    inline std::string get_stats_prefix() const
393    { return stats_prefix; }
set_stats_prefix(std::string stats_prefix_)394    inline void set_stats_prefix(std::string stats_prefix_)
395    { stats_prefix = stats_prefix_; }
unset_stats_prefix()396    inline void unset_stats_prefix()
397    { stats_prefix = ""; }
398 
399 
get_use_acceleration()400    inline int get_use_acceleration() const
401    { return use_acceleration == BOOL_TRUE; }
set_use_acceleration()402    inline void set_use_acceleration()
403    { use_acceleration = BOOL_TRUE; }
unset_use_acceleration()404    inline void unset_use_acceleration()
405    { use_acceleration = BOOL_FALSE; }
406 
get_use_new_stopping_criterion()407    inline int get_use_new_stopping_criterion() const
408    { return use_new_stopping_criterion == BOOL_TRUE; }
set_use_new_stopping_criterion()409    inline void set_use_new_stopping_criterion()
410    { use_new_stopping_criterion = BOOL_TRUE; }
unset_use_new_stopping_criterion()411    inline void unset_use_new_stopping_criterion()
412    { use_new_stopping_criterion = BOOL_FALSE; }
413 
get_store_all_eigenvalues_to_file()414    inline int get_store_all_eigenvalues_to_file() const
415    { return store_all_eigenvalues_to_file == BOOL_TRUE; }
set_store_all_eigenvalues_to_file()416    inline void set_store_all_eigenvalues_to_file()
417    { store_all_eigenvalues_to_file = BOOL_TRUE; }
unset_store_all_eigenvalues_to_file()418    inline void unset_store_all_eigenvalues_to_file()
419    { store_all_eigenvalues_to_file = BOOL_FALSE; }
420 
get_save_permuted_F_matrix_in_bin()421    inline int get_save_permuted_F_matrix_in_bin()
422    { return save_permuted_F_matrix_in_bin == BOOL_TRUE; }
set_save_permuted_F_matrix_in_bin()423    inline void set_save_permuted_F_matrix_in_bin()
424    { save_permuted_F_matrix_in_bin = BOOL_TRUE; }
unset_save_permuted_F_matrix_in_bin()425    inline void unset_save_permuted_F_matrix_in_bin()
426    { save_permuted_F_matrix_in_bin = BOOL_FALSE; }
427 
428 
429 
get_puri_compute_eigv_in_each_iteration()430    inline int get_puri_compute_eigv_in_each_iteration()
431    { return puri_compute_eigv_in_each_iteration == BOOL_TRUE; }
set_puri_compute_eigv_in_each_iteration()432    inline void set_puri_compute_eigv_in_each_iteration()
433    { puri_compute_eigv_in_each_iteration = BOOL_TRUE; }
unset_puri_compute_eigv_in_each_iteration()434    inline void unset_puri_compute_eigv_in_each_iteration()
435    { puri_compute_eigv_in_each_iteration = BOOL_FALSE; }
436 
437 
get_run_shift_and_square_method_on_F()438    inline int get_run_shift_and_square_method_on_F()
439    { return run_shift_and_square_method_on_F == BOOL_TRUE; }
set_run_shift_and_square_method_on_F()440    inline void set_run_shift_and_square_method_on_F()
441    { run_shift_and_square_method_on_F = BOOL_TRUE; }
unset_run_shift_and_square_method_on_F()442    inline void unset_run_shift_and_square_method_on_F()
443    { run_shift_and_square_method_on_F = BOOL_FALSE; }
444 
445 
get_try_eigv_on_next_iteration_if_fail()446    inline int get_try_eigv_on_next_iteration_if_fail()
447    { return try_eigv_on_next_iteration_if_fail == BOOL_TRUE; }
set_try_eigv_on_next_iteration_if_fail()448    inline void set_try_eigv_on_next_iteration_if_fail()
449    { try_eigv_on_next_iteration_if_fail = BOOL_TRUE; }
unset_try_eigv_on_next_iteration_if_fail()450    inline void unset_try_eigv_on_next_iteration_if_fail()
451    { try_eigv_on_next_iteration_if_fail = BOOL_FALSE; }
452 
453 
get_use_prev_vector_as_initial_guess()454    inline int get_use_prev_vector_as_initial_guess()
455    { return use_prev_vector_as_initial_guess == BOOL_TRUE; }
set_use_prev_vector_as_initial_guess()456    inline void set_use_prev_vector_as_initial_guess()
457    { use_prev_vector_as_initial_guess = BOOL_TRUE; }
unset_use_prev_vector_as_initial_guess()458    inline void unset_use_prev_vector_as_initial_guess()
459    { use_prev_vector_as_initial_guess = BOOL_FALSE; }
460 
461 
set_diagonalization_params(ergo_real electronicTemperature_,symmMatrix & overlapMatrix_)462    inline void set_diagonalization_params(ergo_real   electronicTemperature_,
463                                           symmMatrix& overlapMatrix_)
464    {
465       set_overlapMatrix(overlapMatrix_);
466       assert(electronicTemperature_ >= 0);
467       electronicTemperature = electronicTemperature_;
468    }
469 
set_overlapMatrix(symmMatrix & overlapMatrix_)470    inline void set_overlapMatrix(symmMatrix& overlapMatrix_)
471    { overlapMatrix = overlapMatrix_; }
472 
473 
474    inline void set_purification_limits(ergo_real subspaceErrorLimit_,
475                                        ergo_real eigvalueErrorLimit_ = 0,
476                                        ergo_real puri_eig_acc_factor_for_guess = 0)
477    {
478       set_eigvalueErrorLimit(eigvalueErrorLimit_);
479       set_subspaceErrorLimit(subspaceErrorLimit_);
480       set_puri_eig_acc_factor_for_guess(puri_eig_acc_factor_for_guess);
481    }
482 
483    /** Set maximum allowed error in eigenvalues of the density matrix.
484     */
set_eigvalueErrorLimit(ergo_real eigvalueErrorLimit_)485    inline void set_eigvalueErrorLimit(ergo_real eigvalueErrorLimit_)
486    { eigvalueErrorLimit = eigvalueErrorLimit_; }
487 
488    /** Set maximum allowed error in invariant subspaces of the density
489     * matrix.
490     */
set_subspaceErrorLimit(ergo_real subspaceErrorLimit_)491    inline void set_subspaceErrorLimit(ergo_real subspaceErrorLimit_)
492    { subspaceErrorLimit = subspaceErrorLimit_; }
493 
494    /** Set puri_eig_acc_factor_for_guess parameter.
495     *
496     * Obsolete parameter needed for the old stopping criterion for
497     * creating the initial guess.
498     */
set_puri_eig_acc_factor_for_guess(ergo_real puri_eig_acc_factor_for_guess_)499    inline void set_puri_eig_acc_factor_for_guess(ergo_real puri_eig_acc_factor_for_guess_)
500    { puri_eig_acc_factor_for_guess = puri_eig_acc_factor_for_guess_; }
501 
502 
503 
504    // get some results from the purification
505 
get_result_entropy_term()506    ergo_real get_result_entropy_term() const
507    { return resultEntropyTerm; }
508 
get_puri_stats(std::map<std::string,double> & puri_stats_)509    inline void get_puri_stats(std::map<std::string, double>& puri_stats_) const
510    { puri_stats_ = puri_stats; }
511 
512 
513 
514    // Fprev is effective Hamiltonian matrix (=Finput)
515    // Intervals contain the homo and lumo eigenvalues of Fprev
set_eigs_Fprev(intervalType & homoInterval_Finput_,intervalType & lumoInterval_Finput_)516    inline void set_eigs_Fprev(intervalType& homoInterval_Finput_,
517                               intervalType& lumoInterval_Finput_)
518    {
519       homoInterval_Finput = intervalType(homoInterval_Finput_);
520       lumoInterval_Finput = intervalType(lumoInterval_Finput_);
521    }
522 
get_eigs_Fprev(intervalType & homoInterval_Finput_,intervalType & lumoInterval_Finput_)523    inline void get_eigs_Fprev(intervalType& homoInterval_Finput_,
524                               intervalType& lumoInterval_Finput_) const
525    {
526       homoInterval_Finput_ = intervalType(homoInterval_Finput);
527       lumoInterval_Finput_ = intervalType(lumoInterval_Finput);
528    }
529 
530    // F_ort_prev is matrix in orthogonal basis
531    // Intervals contain the homo and lumo eigenvalues of F_ort_prev
set_eigs_F_ort_prev(intervalType & homoInterval_F_ort_prev_,intervalType & lumoInterval_F_ort_prev_)532    inline void set_eigs_F_ort_prev(intervalType& homoInterval_F_ort_prev_,
533                                    intervalType& lumoInterval_F_ort_prev_)
534    {
535       homoInterval_F_ort_prev = intervalType(homoInterval_F_ort_prev_);
536       lumoInterval_F_ort_prev = intervalType(lumoInterval_F_ort_prev_);
537    }
538 
get_eigs_F_ort_prev(intervalType & homoInterval_F_ort_prev_,intervalType & lumoInterval_F_ort_prev_)539    inline void get_eigs_F_ort_prev(intervalType& homoInterval_F_ort_prev_,
540                                    intervalType& lumoInterval_F_ort_prev_) const
541    {
542       homoInterval_F_ort_prev_       = intervalType(homoInterval_F_ort_prev);
543       lumoInterval_F_ort_prev_       = intervalType(lumoInterval_F_ort_prev);
544    }
545 
get_eigvalueErrorLimit()546    inline ergo_real get_eigvalueErrorLimit() const
547    { return eigvalueErrorLimit; }
get_subspaceErrorLimit()548    inline ergo_real get_subspaceErrorLimit() const
549    { return subspaceErrorLimit; }
get_puri_eig_acc_factor_for_guess()550    inline ergo_real get_puri_eig_acc_factor_for_guess() const
551    { return puri_eig_acc_factor_for_guess; }
552 
compute_eigenvectors(std::string eigenvectors_method_,std::string eigenvectors_iterative_method_,ergo_real eigensolver_accuracy_,int eigensolver_maxiter_,int use_prev_vector_as_initial_guess_,int try_eigv_on_next_iteration_if_fail_)553    inline void compute_eigenvectors(std::string eigenvectors_method_,
554                                     std::string eigenvectors_iterative_method_,
555                                     ergo_real eigensolver_accuracy_,
556                                     int eigensolver_maxiter_,
557                                     int use_prev_vector_as_initial_guess_,
558                                     int try_eigv_on_next_iteration_if_fail_)
559    {
560       assert(eigenvectors_method_ == "square" || eigenvectors_method_ == "projection");
561       eigenvectors_method = eigenvectors_method_;
562 
563       assert(eigenvectors_iterative_method_ == "power" || eigenvectors_iterative_method_ == "lanczos");
564       eigenvectors_iterative_method = eigenvectors_iterative_method_;
565 
566       eigensolver_accuracy = eigensolver_accuracy_;
567       eigensolver_maxiter  = eigensolver_maxiter_;
568 
569       if (use_prev_vector_as_initial_guess_ > 0)
570       {
571          set_use_prev_vector_as_initial_guess();
572       }
573       else
574       {
575          unset_use_prev_vector_as_initial_guess();
576       }
577       if (try_eigv_on_next_iteration_if_fail_ > 0)
578       {
579          set_try_eigv_on_next_iteration_if_fail();
580       }
581       else
582       {
583          unset_try_eigv_on_next_iteration_if_fail();
584       }
585 
586    }
587 
compute_eigenvectors_extra(int puri_compute_eigv_in_each_iteration_,int run_shift_and_square_method_on_F_)588    inline void compute_eigenvectors_extra(int puri_compute_eigv_in_each_iteration_, int run_shift_and_square_method_on_F_)
589    {
590       if (puri_compute_eigv_in_each_iteration_ > 0)
591       {
592          set_puri_compute_eigv_in_each_iteration();
593       }
594       else
595       {
596          unset_puri_compute_eigv_in_each_iteration();
597       }
598 
599       if (run_shift_and_square_method_on_F_ > 0)
600       {
601          set_run_shift_and_square_method_on_F();
602       }
603       else
604       {
605          unset_run_shift_and_square_method_on_F();
606       }
607    }
608 
609 private:
610 
611    int SCF_step;
612 
613 
614 
615    bool use_diagonalization;                  /**< Flag to turn on diagonalization. */
616 
617    bool use_purification;                     /**< Flag to turn on purification. */
618 
619    bool store_all_eigenvalues_to_file;        /**< Store eigenvalues to the file when doing diagonalization.
620                                               * NOTE: works just with diagonalization */
621 
622    bool try_eigv_on_next_iteration_if_fail;   /**< For square method: if eigenvector is not computed in iteration
623                                               * i, try to compute it in iteration i+1 */
624 
625    ergo_real electronicTemperature;          /**< Electronic temperature */
626 
627    ergo_real gap_expected_lower_bound;       /**< Expected lower bound for the gap to be used in early iterations.  */
628 
629    ergo_real eigvalueErrorLimit;             /**< Tolerated deviation of eigenvalues from 0 and 1 in the computed density matrix.  */
630 
631    ergo_real subspaceErrorLimit;             /**< Tolerated error in the occupied subspace as measured by the sinus of the largest canonical angle. */
632 
633    ergo_real puri_eig_acc_factor_for_guess;  /**< With this number will be multiplied the tolerated deviation of
634                                               * eigenvalues from 0 and 1 in the computed density matrix for the
635                                               * initial guess density matrix */
636 
637    bool use_diag_on_error;                    /**< Flag to fall back on diagonalization if purification fails. */
638    bool use_diag_on_error_guess;
639 
640    bool create_m_files;                       /**< Flag to create m-files with information about the purification process.  */
641 
642    bool output_homo_and_lumo_eigenvectors;    /**< Compute homo and lumo eigenvectors and write them to the file */
643 
644    int number_of_occupied_eigenvectors;    /**< Number of occupied eigenvectors to compute. */
645 
646    int number_of_unoccupied_eigenvectors;    /**< Number of unoccupied eigenvectors to compute. */
647 
648    int go_back_X_iter_proj_method; /**< Parameter used in the projection method for computing eigenvectors. Defines the iteration to start computation of the eigenvectors. */
649 
650    int jump_over_X_iter_proj_method;  /**< Parameter used in the projection method for computing eigenvectors. Defines how many iterations will be skipped before the next attempt if some eigenvectors are not computed. */
651 
652    bool use_prev_vector_as_initial_guess;     /**< Use eigenvector from the previous SCF cycle as an initial guess in this cycle */
653 
654    bool puri_compute_eigv_in_each_iteration;  /**< Compute eigenvectors in each iteration of the recursive expansion. */
655 
656    bool run_shift_and_square_method_on_F;     /**< (for comparison) Run shift_and_square method to
657                                               * get eigenvectors of the matrix F for various shifts. */
658 
659    bool save_permuted_F_matrix_in_bin;        /**< Save sparse matrix F into bin file in the current permutation of rows and columns. */
660 
661    bool ignore_purification_failure;          /**< Continue even if purification fails to converge.  */
662 
663    bool use_rand_perturbation_for_alleigsint; /**< Apply a random
664                                               * perturbation to (try
665                                               * to) improve the
666                                               * convergence speed of
667                                               * Lanczos calculation of
668                                               * extremal
669                                               * eigenvalues.  */
670 
671    std::string stats_prefix;                 /**< Prefix to be added to statistics files. */
672 
673    bool plot_puri_results;                    /**< Plot results of the purification from this function call */
674 
675    std::string plot_puri_results_str;
676 
677    bool use_acceleration;                        /**< Use acceleration in the purification */
678 
679    bool use_new_stopping_criterion;              /**< Use new parameterless stopping criterion */
680 
681    std::string eigenvectors_method;             /**< Method for computing eigenvectors: square or projection */
682 
683    std::string eigenvectors_iterative_method;   /**< Iterative method for computing eigenvectors: power or lanczos */
684 
685    ergo_real eigensolver_accuracy;              /**< The accuracy for the eigenvalue problem solver */
686 
687    int eigensolver_maxiter;                     /**< Maximum number of iterations for the eigenvalue problem solver */
688 
689    int n;                                       /**< System size. */
690 
691    int noOfOccupiedOrbs;                        /**< Number of occupied orbitals.  */
692 
693    ergo_real factor;                            /**< Factor to scale the resulting density matrix.
694                                                  * (for restricted vs unrestricted calc) */
695 
696    symmMatrix overlapMatrix;                    /**< Overlap matrix (written to file) */
697 
698    symmMatrix D_ort_prev;                       /**< Density matrix from previous SCF cycle (written to file) */
699 
700    triangMatrix invCholFactor;                  /**< Inverse Cholesky factor (written to file) */
701 
702    ergo_real invCholFactor_euclnorm;            /**< Euclidean norm of inverse Cholesky factor. */
703 
704    mat::normType truncationNormPurification;    /**< Norm to be used for truncation.  */
705 
706    mat::normType stopCriterionNormPurification; /**< Norm to be used for stopping criterion.  */
707 
708    int maxMul;                                  /**< Maximum allowed number of matrix multiplications in the purification */
709 
710    mat::SizesAndBlocks matrixSizesAndBlocks;    /**< Information about HML matrix block sizes etc. */
711 
712    int leavesSizeMax;                           /**< Information about leavesSizeMax and blocksize for CHTMatrix */
713    int blocksize;                               /**< Information about leavesSizeMax and blocksize for CHTMatrix */
714 
715    intervalType homoInterval_Finput;
716    intervalType lumoInterval_Finput;
717 
718    intervalType homoInterval_F_ort_prev;
719    intervalType lumoInterval_F_ort_prev;
720 
721    ergo_real resultEntropyTerm;
722    std::map<std::string, double> puri_stats;
723 
724    std::vector<generalVector> eigVecOCC;
725    std::vector<generalVector> eigVecUNOCC;
726    std::vector<ergo_real> eigValOCC;
727    std::vector<ergo_real> eigValUNOCC;
728 
729    // Names of files needed for checkpoints
730    const char *filenameFinput;
731    const char *filenameF_ort_prev;
732    const char *filenameeigVecLUMO;
733    const char *filenameeigVecHOMO;
734    const char *filenameOverlap;
735    const char *filenameD_ort_prev;
736    const char *filenameinvCholFactor;
737    const char *file_for_basic_types;
738 };
739 
740 
741 
742 #endif // GETDENSFROMFOCKHEADER
743