diff --git a/src/Control.cc b/src/Control.cc index 22de9b27..57e12336 100644 --- a/src/Control.cc +++ b/src/Control.cc @@ -224,6 +224,7 @@ void Control::print(std::ostream& os) << conv_tol << std::endl; os << std::fixed; os << " Density matrix mixing = " << dm_mix << std::endl; + os << " Density matrix tol = " << dm_tol << std::endl; if (DMEigensolver() == DMEigensolverType::Eigensolver) { os << " Density matrix computation algorithm = " @@ -439,7 +440,7 @@ void Control::sync(void) memset(&int_buffer[0], 0, size_int_buffer * sizeof(int)); } - const short size_float_buffer = 43; + const short size_float_buffer = 44; float* float_buffer = new float[size_float_buffer]; if (mype_ == 0) { @@ -485,6 +486,7 @@ void Control::sync(void) float_buffer[40] = threshold_eigenvalue_gram_quench_; float_buffer[41] = pair_mlwf_distance_threshold_; float_buffer[42] = e0_; + float_buffer[43] = dm_tol; } else { @@ -680,6 +682,7 @@ void Control::sync(void) threshold_eigenvalue_gram_quench_ = float_buffer[40]; pair_mlwf_distance_threshold_ = float_buffer[41]; e0_ = float_buffer[42]; + dm_tol = float_buffer[43]; max_electronic_steps_loose_ = max_electronic_steps; delete[] short_buffer; @@ -1720,6 +1723,7 @@ void Control::setOptions(const boost::program_options::variables_map& vm) lrs_extrapolation = 10; dm_mix = vm["DensityMatrix.mixing"].as(); + dm_tol = vm["DensityMatrix.tol"].as(); dm_inner_steps = vm["DensityMatrix.nb_inner_it"].as(); dm_use_old_ = vm["DensityMatrix.use_old"].as() ? 1 : 0; str = vm["DensityMatrix.algo"].as(); diff --git a/src/MVPSolver.cc b/src/MVPSolver.cc index 770acdd0..1c093e53 100644 --- a/src/MVPSolver.cc +++ b/src/MVPSolver.cc @@ -47,14 +47,16 @@ MVPSolver::MVPSolver(MPI_Comm comm, std::ostream& os, Electrostatic* electrostat, MGmol* mgmol_strategy, const int numst, const double kbT, const std::vector>& global_indexes, - const short n_inner_steps, const double mixing, const bool use_old_dm) + const short n_inner_steps, const double mixing, const double tol_de0, + const bool use_old_dm) : comm_(comm), os_(os), n_inner_steps_(n_inner_steps), use_old_dm_(use_old_dm), ions_(ions), numst_(numst), - mixing_(mixing) + mixing_(mixing), + tol_de0_(tol_de0) { Control& ct = *(Control::instance()); if (onpe0 && ct.verbose > 0) @@ -208,7 +210,6 @@ int MVPSolver::solve(OrbitalsType& orbitals) kbpsi.computeHvnlMatrix(&kbpsi, ions_, h11_nl); - const double tol_de0 = 1.e-12; for (int inner_it = 0; inner_it < n_inner_steps_; inner_it++) { if (onpe0 && ct.verbose > 1) @@ -268,13 +269,26 @@ int MVPSolver::solve(OrbitalsType& orbitals) MatrixType delta_dm("delta_dm", numst_, numst_); delta_dm = target; delta_dm -= dmInit; + + double de0 = evaluateDerivative(dmInit, delta_dm, ts0); + + // check for convergence + if (std::abs(de0) < tol_de0_ && inner_it > 0) + { + if (onpe0 && ct.verbose > 0) + std::cout << "MVP: de0 = " << de0 + << ", convergence achieved" << std::endl; + break; + } + double beta = 0.; if (mixing_ > 0.) { beta = mixing_; - if (onpe0 && ct.verbose > 1) + if (onpe0 && ct.verbose > 0) { - os_ << "MVP with beta = " << beta << std::endl; + if (ct.verbose > 1) + os_ << "MVP with beta = " << beta << std::endl; os_ << std::setprecision(12); os_ << std::fixed << "MVP inner iteration " << inner_it << ", E0=" << e0 << std::endl; @@ -282,16 +296,6 @@ int MVPSolver::solve(OrbitalsType& orbitals) } else { - double de0 = evaluateDerivative(dmInit, delta_dm, ts0); - - if (std::abs(de0) < tol_de0 && inner_it > 0) - { - if (onpe0 && ct.verbose > 0) - std::cout << "MVP: de0 = " << de0 - << ", convergence achieved" << std::endl; - break; - } - // // evaluate free energy at beta=1 // diff --git a/src/MVPSolver.h b/src/MVPSolver.h index 5558bcdd..77a6a345 100644 --- a/src/MVPSolver.h +++ b/src/MVPSolver.h @@ -35,6 +35,11 @@ class MVPSolver int numst_; double mixing_; + /*! + * tolerance on energy slope in inner iterations + */ + double tol_de0_; + Rho* rho_; Energy* energy_; Electrostatic* electrostat_; @@ -57,7 +62,8 @@ class MVPSolver Electrostatic* electrostat, MGmol* mgmol_strategy, const int numst, const double kbT, const std::vector>& global_indexes, - const short n_inner_steps, const double mixing, const bool use_old_dm); + const short n_inner_steps, const double mixing, const double tol_de0, + const bool use_old_dm); ~MVPSolver(); int solve(OrbitalsType& orbitals); diff --git a/src/MVP_DMStrategy.cc b/src/MVP_DMStrategy.cc index 899bca46..3fb53dee 100644 --- a/src/MVP_DMStrategy.cc +++ b/src/MVP_DMStrategy.cc @@ -54,7 +54,7 @@ int MVP_DMStrategy::update(OrbitalsType& orbitals) MVPSolver solver(comm_, os_, ions_, rho_, energy_, electrostat_, mgmol_strategy_, ct.numst, ct.occ_width, global_indexes_, - ct.dm_inner_steps, ct.dm_mix, use_old_dm_); + ct.dm_inner_steps, ct.dm_mix, ct.dm_tol, use_old_dm_); return solver.solve(orbitals); } diff --git a/src/read_config.cc b/src/read_config.cc index 642b1b99..ed0e66e1 100644 --- a/src/read_config.cc +++ b/src/read_config.cc @@ -297,7 +297,9 @@ int read_config(int argc, char** argv, po::variables_map& vm, po::value()->default_value(0), "Flag for computing new centers from extrapolated orbitals.")( "DensityMatrix.mixing", po::value()->default_value(-1.), - "Mixing coefficient for Density Matrix")("DensityMatrix.solver", + "Mixing coefficient for Density Matrix")("DensityMatrix.tol", + po::value()->default_value(1.e-12), + "Tolerance for Density Matrix convergence")("DensityMatrix.solver", po::value()->default_value("Mixing"), "Algorithm for updating Density Matrix: Mixing, MVP, HMVP")( "DensityMatrix.nb_inner_it", po::value()->default_value(3),