add bcast for energy

Author: Qianruipku

source/module_esolver/esolver_ks_pw.cpp

@@ -441,6 +441,7 @@ void ESolver_KS_PW<T, Device>::update_pot(const int istep, const int iter)
         }
         this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell);
         this->pelec->f_en.descf = this->pelec->cal_delta_escf();
+        MPI_Bcast(&(this->pelec->f_en.descf), 1, MPI_DOUBLE, 0, PARAPW_WORLD);
     }
     else
     {

source/module_esolver/esolver_sdft_pw.cpp

@@ -221,6 +221,7 @@ void ESolver_SDFT_PW::hamilt2density(int istep, int iter, double ethr)
             srho.begin(is, *(this->pelec->charge), pw_rho, GlobalC::ucell.symm);
         }
         this->pelec->f_en.deband = this->pelec->cal_delta_eband();
+        MPI_Bcast(&(this->pelec->f_en.deband), 1, MPI_DOUBLE, 0, PARAPW_WORLD);
     }
     else
     {

source/module_hamilt_pw/hamilt_stodft/sto_iter.cpp

@@ -414,7 +414,6 @@ void Stochastic_Iter::calHsqrtchi(Stochastic_WF& stowf)
     p_che->calcoef_real(nroot_fd);
     for (int ik = 0; ik < this->pkv->get_nks(); ++ik)
     {
-        p_hamilt_sto->updateHk(ik);
         this->calTnchi_ik(ik, stowf);
     }
 }
@@ -481,7 +480,7 @@ void Stochastic_Iter::sum_stoband(Stochastic_WF& stowf,
             const int nchip_ik = nchip[ik];
             if (this->pkv->get_nks() > 1)
             {
-                pHamilt->updateHk(ik);
+                pHamilt->updateHk(ik); //can be merged with calTnchi_ik, but it does not nearly cost time.
                 stowf.shchi->fix_k(ik);
             }
             const int npw = this->pkv->ngk[ik];
@@ -632,6 +631,10 @@ void Stochastic_Iter::calTnchi_ik(const int& ik, Stochastic_WF& stowf)
     }
     else
     {
+        if (this->pkv->get_nks() > 1)
+        {
+            p_hamilt_sto->updateHk(ik); //necessary, because itermu should be called before this function
+        }
         auto hchi_norm = std::bind(&hamilt::HamiltSdftPW<std::complex<double>>::hPsi_norm,
                                    p_hamilt_sto,
                                    std::placeholders::_1,

source/module_hsolver/hsolver_pw_sdft.cpp

@@ -68,12 +68,15 @@ void HSolverPW_SDFT::solve(hamilt::Hamilt<std::complex<double>>* pHamilt,
         // init k
         if (nks > 1)
         {
-            pHamilt->updateHk(ik);
+            pHamilt->updateHk(ik); // necessary , because emax and emin should be decided first
         }
         stoiter.calPn(ik, stowf);
     }
 
+    // iterate to get mu
     stoiter.itermu(iter, pes);
+
+    // prepare sqrt{f(\hat{H})}|\chi> to calculate density, force and stress
     stoiter.calHsqrtchi(stowf);
     if (skip_charge)
     {