return ne_diff

Author: simonpintarelli

src/dft/dft_ground_state.cpp

@@ -69,7 +69,7 @@ DFT_ground_state::energy_kin_sum_pw() const
     for (auto it : kset_.spl_num_kpoints()) {
         auto kp = kset_.get<double>(it.i);
 
-        #pragma omp parallel for schedule(static) reduction(+:ekin)
+        #pragma omp parallel for schedule(static) reduction(+ : ekin)
         for (int igloc = 0; igloc < kp->num_gkvec_loc(); igloc++) {
             auto Gk = kp->gkvec().gkvec_cart(gvec_index_t::local(igloc));
 
@@ -201,6 +201,8 @@ DFT_ground_state::find(double density_tol__, double energy_tol__, double iter_so
       << "num_dft_iter              : " << num_dft_iter__;
     ctx_.message(1, __func__, s);
 
+    double ne_diff = 0;
+
     for (int iter = 0; iter < num_dft_iter__; iter++) {
         PROFILE("sirius::DFT_ground_state::scf_loop|iteration");
         std::stringstream s;
@@ -225,7 +227,7 @@ DFT_ground_state::find(double density_tol__, double energy_tol__, double iter_so
                                                             ctx_.cfg().iterative_solver().num_steps());
             }
             /* find band occupancies */
-            kset_.find_band_occupancies<float>();
+            ne_diff = kset_.find_band_occupancies<float>();
             /* generate new density from the occupied wave-functions */
             density_.generate<float>(kset_, ctx_.use_symmetry(), true, true);
 #else
@@ -237,7 +239,7 @@ DFT_ground_state::find(double density_tol__, double energy_tol__, double iter_so
             result = sirius::diagonalize<double, double>(H0, kset_, iter_solver_tol__,
                                                          ctx_.cfg().iterative_solver().num_steps());
             /* find band occupancies */
-            kset_.find_band_occupancies<double>();
+            ne_diff = kset_.find_band_occupancies<double>();
             /* generate new density from the occupied wave-functions */
             density_.generate<double>(kset_, ctx_.use_symmetry(), true, true);
         }
@@ -356,6 +358,12 @@ DFT_ground_state::find(double density_tol__, double energy_tol__, double iter_so
             converged = converged && (rms < density_tol__);
         }
         if (converged) {
+            if (std::abs(ne_diff) > 0) {
+                std::stringstream ss;
+                ss << "Newton minimization didn't respect correct number of electrons, ne_diff=" << ne_diff;
+                ss << "\nReduce smearing width!";
+                RTE_THROW(ss.str());
+            }
             std::stringstream out;
             out << std::endl;
             out << "converged after " << iter + 1 << " SCF iterations!" << std::endl;

src/k_point/k_point_set.cpp

@@ -199,7 +199,7 @@ bisection_search(F&& f, double a, double b, double tol, int maxstep = 1000)
 }
 
 template <typename T>
-void
+double
 K_point_set::find_band_occupancies()
 {
     PROFILE("sirius::K_point_set::find_band_occupancies");
@@ -209,7 +209,7 @@ K_point_set::find_band_occupancies()
     auto band_occ_callback = ctx_.band_occ_callback();
     if (band_occ_callback) {
         band_occ_callback();
-        return;
+        return 0;
     }
 
     /* target number of electrons */
@@ -239,7 +239,7 @@ K_point_set::find_band_occupancies()
         }
 
         this->sync_band<T, sync_band_t::occupancy>();
-        return;
+        return 0;
     }
 
     if (ctx_.smearing_width() == 0) {
@@ -298,6 +298,7 @@ K_point_set::find_band_occupancies()
     }
     this->energy_fermi_ = res_bisection.value();
 
+    double ne_diff = 0;
     /* for cold and Methfessel Paxton smearing start newton minimization  */
     if (ctx_.smearing() == smearing::smearing_t::cold || ctx_.smearing() == smearing::smearing_t::methfessel_paxton) {
         f               = smearing::occupancy(ctx_.smearing(), ctx_.smearing_width());
@@ -313,6 +314,7 @@ K_point_set::find_band_occupancies()
                 RTE_OUT(ctx_.out()) << "newton iteration converged after " << res_newton.value().iter << " steps\n";
                 RTE_OUT(ctx_.out()) << "newton iteration ne_diff " << res_newton.value().ne_diff << " steps\n";
             }
+            ne_diff = res_newton.value().ne_diff;
         } else {
             // Newton has failed, fallback to bisection
             if (ctx_.verbosity() >= 2) {
@@ -373,12 +375,14 @@ K_point_set::find_band_occupancies()
             band_gap_ = eband[ist].first - eband[ist - 1].second;
         }
     }
+    return ne_diff;
 }
 
-template void
+
+template double
 K_point_set::find_band_occupancies<double>();
 #if defined(SIRIUS_USE_FP32)
-template void
+template double
 K_point_set::find_band_occupancies<float>();
 #endif
 

src/k_point/k_point_set.hpp

@@ -109,9 +109,9 @@ class K_point_set
     void
     sync_band();
 
-    /// Find Fermi energy and band occupation numbers.
+    /// Find Fermi energy and band occupation numbers. Returns ne_diff (Newton protocol)
     template <typename T>
-    void
+    double
     find_band_occupancies();
 
     /// Print basic info to the standard output.