fix wrong nscf energy of sDFT (#5618)

Qianruipku · web-flow · commit 0f569ff1bb06 · 2024-11-26T22:11:36.000+08:00
add nscf energy tests
diff --git a/source/module_elecstate/elecstate_pw.cpp b/source/module_elecstate/elecstate_pw.cpp
@@ -135,11 +135,6 @@ void ElecStatePW<T, Device>::parallelK()
 {
 #ifdef __MPI
     this->charge->rho_mpi();
-    if(PARAM.inp.esolver_type == "sdft") //qinarui add it 2021-7-21
-	{
-        this->f_en.eband /= GlobalV::NPROC_IN_POOL;
-        MPI_Allreduce(MPI_IN_PLACE, &this->f_en.eband, 1, MPI_DOUBLE, MPI_SUM, STO_WORLD);
-    }
 #endif
 }
 
diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp
@@ -442,15 +442,16 @@ void ESolver_KS<T, Device>::runner(UnitCell& ucell, const int istep)
     this->diag_ethr = PARAM.inp.pw_diag_thr;
     for (int iter = 1; iter <= this->maxniter; ++iter)
     {
-        // 6) initialization of SCF iterations
+        // 5) initialization of SCF iterations
         this->iter_init(ucell, istep, iter);
 
+        // 6) use Hamiltonian to obtain charge density
         this->hamilt2density(ucell, istep, iter, diag_ethr);
 
-        // 10) finish scf iterations
+        // 7) finish scf iterations
         this->iter_finish(ucell, istep, iter);
 
-        // 13) check convergence
+        // 8) check convergence
         if (this->conv_esolver || this->oscillate_esolver)
         {
             this->niter = iter;
@@ -462,7 +463,7 @@ void ESolver_KS<T, Device>::runner(UnitCell& ucell, const int istep)
         }
     } // end scf iterations
 
-    // 15) after scf
+    // 9) after scf
     ModuleBase::timer::tick(this->classname, "after_scf");
     this->after_scf(ucell, istep);
     ModuleBase::timer::tick(this->classname, "after_scf");
diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_iter.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_iter.cpp
@@ -472,13 +472,13 @@ void Stochastic_Iter<T, Device>::calHsqrtchi(Stochastic_WF<T, Device>& stowf)
 }
 
 template <typename T, typename Device>
-void Stochastic_Iter<T, Device>::sum_stoband(Stochastic_WF<T, Device>& stowf,
+void Stochastic_Iter<T, Device>::sum_stoeband(Stochastic_WF<T, Device>& stowf,
                                              elecstate::ElecStatePW<T, Device>* pes,
                                              hamilt::Hamilt<T, Device>* pHamilt,
                                              ModulePW::PW_Basis_K* wfc_basis)
 {
-    ModuleBase::TITLE("Stochastic_Iter", "sum_stoband");
-    ModuleBase::timer::tick("Stochastic_Iter", "sum_stoband");
+    ModuleBase::TITLE("Stochastic_Iter", "sum_stoeband");
+    ModuleBase::timer::tick("Stochastic_Iter", "sum_stoeband");
     const int npwx = wfc_basis->npwk_max;
     const int norder = p_che->norder;
 
@@ -556,7 +556,7 @@ void Stochastic_Iter<T, Device>::sum_stoband(Stochastic_WF<T, Device>& stowf,
     MPI_Allreduce(MPI_IN_PLACE, &sto_eband, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 #endif
     pes->f_en.eband += sto_eband;
-    ModuleBase::timer::tick("Stochastic_Iter", "sum_stoband");
+    ModuleBase::timer::tick("Stochastic_Iter", "sum_stoeband");
 }
 
 template <typename T, typename Device>
diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_iter.h b/source/module_hamilt_pw/hamilt_stodft/sto_iter.h
@@ -52,7 +52,7 @@ class Stochastic_Iter
      * @param pHamilt hamiltonian
      * @param wfc_basis wfc pw basis
      */
-    void sum_stoband(Stochastic_WF<T, Device>& stowf,
+    void sum_stoeband(Stochastic_WF<T, Device>& stowf,
                      elecstate::ElecStatePW<T, Device>* pes,
                      hamilt::Hamilt<T, Device>* pHamilt,
                      ModulePW::PW_Basis_K* wfc_basis);
diff --git a/source/module_hsolver/hsolver_pw_sdft.cpp b/source/module_hsolver/hsolver_pw_sdft.cpp
@@ -108,29 +108,39 @@ void HSolverPW_SDFT<T, Device>::solve(hamilt::Hamilt<T, Device>* pHamilt,
     // prepare sqrt{f(\hat{H})}|\chi> to calculate density, force and stress
     stoiter.calHsqrtchi(stowf);
 
+    // calculate eband = \sum_{ik,ib} w(ik)f(ik,ib)e_{ikib}, demet = -TS
     elecstate::ElecStatePW<T, Device>* pes_pw = static_cast<elecstate::ElecStatePW<T, Device>*>(pes);
     if (GlobalV::MY_STOGROUP == 0)
     {
         pes_pw->calEBand();
     }
+    if (nbands > 0)
+    {
+#ifdef __MPI
+        pes->f_en.eband /= GlobalV::NPROC_IN_POOL;
+        MPI_Allreduce(MPI_IN_PLACE, &pes->f_en.eband, 1, MPI_DOUBLE, MPI_SUM, STO_WORLD);
+        MPI_Bcast(&pes->f_en.eband, 1, MPI_DOUBLE, 0, PARAPW_WORLD);
+#endif
+    }
+    stoiter.sum_stoeband(stowf, pes_pw, pHamilt, wfc_basis);
+    
+    
+
+    // for nscf, skip charge
     if (skip_charge)
     {
         ModuleBase::timer::tick("HSolverPW_SDFT", "solve");
         return;
     }
+
     //(5) calculate new charge density
     // calculate KS rho.
     pes_pw->init_rho_data();
     if (nbands > 0)
     {
         pes_pw->psiToRho(psi);
-#ifdef __MPI
-        MPI_Bcast(&pes->f_en.eband, 1, MPI_DOUBLE, 0, PARAPW_WORLD);
-#endif
     }
-    
     // calculate stochastic rho
-    stoiter.sum_stoband(stowf, pes_pw, pHamilt, wfc_basis);
     stoiter.cal_storho(stowf, pes_pw, wfc_basis);
 
     // will do rho symmetry and energy calculation in esolver
diff --git a/source/module_hsolver/test/test_hsolver_sdft.cpp b/source/module_hsolver/test/test_hsolver_sdft.cpp
@@ -141,7 +141,7 @@ void Stochastic_Iter<T, Device>::calHsqrtchi(Stochastic_WF<T, Device>& stowf)
 }
 
 template <typename T, typename Device>
-void Stochastic_Iter<T, Device>::sum_stoband(Stochastic_WF<T, Device>& stowf,
+void Stochastic_Iter<T, Device>::sum_stoeband(Stochastic_WF<T, Device>& stowf,
                                              elecstate::ElecStatePW<T, Device>* pes,
                                              hamilt::Hamilt<T, Device>* pHamilt,
                                              ModulePW::PW_Basis_K* wfc_basis)
diff --git a/tests/integrate/107_PW_OBOD_MemSaver/result.ref b/tests/integrate/107_PW_OBOD_MemSaver/result.ref
@@ -1,3 +1,5 @@
+etotref -206.9728084455294
+etotperatomref -103.4864042228
 totaldosref 1190.44
 CompareBand_pass 0
-totaltimeref 0.33861
+totaltimeref 0.38
diff --git a/tests/integrate/107_PW_OB_outputbands/result.ref b/tests/integrate/107_PW_OB_outputbands/result.ref
@@ -1,2 +1,4 @@
+etotref -206.9727001566756
+etotperatomref -103.4863500783
 CompareBand_pass 0
-totaltimeref 0.28457
+totaltimeref 0.23
diff --git a/tests/integrate/107_PW_OD_outputdos/result.ref b/tests/integrate/107_PW_OD_outputdos/result.ref
@@ -1,2 +1,4 @@
+etotref -206.9727320695148
+etotperatomref -103.4863660348
 totaldosref 1190.44
-totaltimeref 
+totaltimeref 0.26
diff --git a/tests/integrate/107_PW_W90/result.ref b/tests/integrate/107_PW_W90/result.ref
@@ -1,4 +1,6 @@
+etotref -308.0061503011041
+etotperatomref -154.0030751506
 CompareAMN_pass 0
 CompareMMN_pass 0
 CompareEIG_pass 0
-totaltimeref 9.3341
+totaltimeref 2.15
diff --git a/tests/integrate/117_PW_out_pot_nscf/result.ref b/tests/integrate/117_PW_out_pot_nscf/result.ref
@@ -1,5 +1,7 @@
+etotref -3526.200848639589
+etotperatomref -3526.2008486396
 ComparePot1_pass 0
 pointgroupref O_h
 spacegroupref O_h
 nksibzref 3
-totaltimeref 0.35
+totaltimeref 0.33
diff --git a/tests/integrate/186_PW_NSCF_KG_100/result.ref b/tests/integrate/186_PW_NSCF_KG_100/result.ref
@@ -1,5 +1,7 @@
+etotref -150.4802162021569
+etotperatomref -150.4802162022
 CompareH_Failed 0
 pointgroupref O_h
 spacegroupref O_h
 nksibzref 3
-totaltimeref 0.67841
+totaltimeref 1.46
diff --git a/tests/integrate/186_PW_SDOS_MALL/result.ref b/tests/integrate/186_PW_SDOS_MALL/result.ref
@@ -1,7 +1,7 @@
-etotref 0
-etotperatomref 0.0000000000
+etotref -150.4802139726653
+etotperatomref -150.4802139727
 totaldosref 96.4227
 pointgroupref O_h
 spacegroupref O_h
 nksibzref 3
-totaltimeref 5.1234
+totaltimeref 11.83
diff --git a/tests/integrate/186_PW_SKG_10D10S/result.ref b/tests/integrate/186_PW_SKG_10D10S/result.ref
@@ -1,6 +1,7 @@
+etotref -150.5555663868826
+etotperatomref -150.5555663869
 CompareH_Failed 0
 pointgroupref O_h
 spacegroupref O_h
 nksibzref 3
-totaltimeref 58.0934
-4.88
+totaltimeref 3.67
diff --git a/tests/integrate/207_NO_KP_OB/result.ref b/tests/integrate/207_NO_KP_OB/result.ref
@@ -1,3 +1,5 @@
+etotref -205.4532897097782
+etotperatomref -102.7266448549
 CompareOrb_pass 0
 CompareBand_pass 0
-totaltimeref 0.27705
+totaltimeref 0.31
diff --git a/tests/integrate/207_NO_KP_OD/result.ref b/tests/integrate/207_NO_KP_OD/result.ref
@@ -1,2 +1,4 @@
+etotref -205.4532897097782
+etotperatomref -102.7266448549
 totaldosref 1189.14
-totaltimeref 0.31661
+totaltimeref 0.31
diff --git a/tests/integrate/207_NO_KP_OH_nscf/result.ref b/tests/integrate/207_NO_KP_OH_nscf/result.ref
@@ -1,4 +1,6 @@
+etotref -205.4532897097782
+etotperatomref -102.7266448549
 CompareSR_pass 0
 ComparerR_pass 0
 ComparerTR_pass 0
-totaltimeref 1.88
+totaltimeref 2.21
diff --git a/tests/integrate/214_NO_mulliken_nscf/result.ref b/tests/integrate/214_NO_mulliken_nscf/result.ref
@@ -1,2 +1,4 @@
+etotref -50.22627762925340
+etotperatomref -25.1131388146
 Compare_mulliken_pass 0
-totaltimeref 0.20
+totaltimeref 0.76
diff --git a/tests/integrate/217_NO_out_pot_nscf/result.ref b/tests/integrate/217_NO_out_pot_nscf/result.ref
@@ -1,5 +1,7 @@
+etotref -196.6254071148604
+etotperatomref -98.3127035574
 ComparePot1_pass 0
 pointgroupref T_d
 spacegroupref O_h
 nksibzref 1
-totaltimeref 0.49
+totaltimeref 0.35
diff --git a/tests/integrate/tools/catch_properties.sh b/tests/integrate/tools/catch_properties.sh
@@ -83,7 +83,7 @@ if [ ! -z $esolver_type ] && ([ $esolver_type == "lr" ] || [ $esolver_type == "k
 	is_lr=1
 fi
 
-if [ $calculation != "nscf" ] && [ $calculation != "get_wf" ]\
+if [ $calculation != "get_wf" ]\
 && [ $calculation != "get_pchg" ] && [ $calculation != "get_S" ]\
 && [ $is_lr == 0 ]; then
 	#etot=`grep ETOT_ $running_path | awk '{print $2}'` 

Original file line number	Diff line number	Diff line change
`@@ -135,11 +135,6 @@ void ElecStatePW<T, Device>::parallelK()`
`135`	`135`	`{`
`136`	`136`	`#ifdef __MPI`
`137`	`137`	`this->charge->rho_mpi();`
`138`		`- if(PARAM.inp.esolver_type == "sdft") //qinarui add it 2021-7-21`
`139`		`- {`
`140`		`- this->f_en.eband /= GlobalV::NPROC_IN_POOL;`
`141`		`- MPI_Allreduce(MPI_IN_PLACE, &this->f_en.eband, 1, MPI_DOUBLE, MPI_SUM, STO_WORLD);`
`142`		`- }`
`143`	`138`	`#endif`
`144`	`139`	`}`
`145`	`140`
Original file line number	Diff line number	Diff line change
`@@ -141,7 +141,7 @@ void Stochastic_Iter<T, Device>::calHsqrtchi(Stochastic_WF<T, Device>& stowf)`
`141`	`141`	`}`
`142`	`142`
`143`	`143`	`template <typename T, typename Device>`
`144`		`-void Stochastic_Iter<T, Device>::sum_stoband(Stochastic_WF<T, Device>& stowf,`
	`144`	`+void Stochastic_Iter<T, Device>::sum_stoeband(Stochastic_WF<T, Device>& stowf,`
`145`	`145`	`elecstate::ElecStatePW<T, Device>* pes,`
`146`	`146`	`hamilt::Hamilt<T, Device>* pHamilt,`
`147`	`147`	`ModulePW::PW_Basis_K* wfc_basis)`