Merge branch 'develop' into mybfgs

Author: 19hello

source/Makefile.Objects

@@ -585,7 +585,6 @@ OBJS_LCAO=evolve_elec.o\
       LCAO_set_fs.o\
       LCAO_set_st.o\
       LCAO_nl_mu.o\
-      LCAO_nnr.o\
       LCAO_set_zero.o\
       LCAO_allocate.o\
       LCAO_set_mat2d.o\

source/module_base/formatter.h

@@ -145,6 +145,20 @@ class FmtCore
             [&delim](const std::string& acc, const std::string& s) { return acc + delim + s; });
     }
 
+    static std::string upper(const std::string& in)
+    {
+        std::string dst = in;
+        std::transform(dst.begin(), dst.end(), dst.begin(), ::toupper);
+        return dst;
+    }
+
+    static std::string lower(const std::string& in)
+    {
+        std::string dst = in;
+        std::transform(dst.begin(), dst.end(), dst.begin(), ::tolower);
+        return dst;
+    }
+
 private:
     std::string fmt_;
     template<typename T>

source/module_base/module_mixing/test/mixing_test.cpp

@@ -1,4 +1,6 @@
+#ifdef _OPENMP
 #include <omp.h>
+#endif
 
 #include "../broyden_mixing.h"
 #include "../plain_mixing.h"
@@ -151,7 +153,9 @@ class Mixing_Test : public testing::Test
 
 TEST_F(Mixing_Test, BroydenSolveLinearEq)
 {
+#ifdef _OPENMP
     omp_set_num_threads(1);
+#endif
     init_method("broyden");
     std::vector<double> x_in = xd_ref;
     std::vector<double> x_out(3);
@@ -196,7 +200,9 @@ TEST_F(Mixing_Test, BroydenSolveLinearEq)
 
 TEST_F(Mixing_Test, PulaySolveLinearEq)
 {
+#ifdef _OPENMP
     omp_set_num_threads(1);
+#endif
     init_method("pulay");
     std::vector<double> x_in = xd_ref;
     std::vector<double> x_out(3);
@@ -242,7 +248,9 @@ TEST_F(Mixing_Test, PulaySolveLinearEq)
 
 TEST_F(Mixing_Test, PlainSolveLinearEq)
 {
+#ifdef _OPENMP
     omp_set_num_threads(1);
+#endif
     init_method("plain");
     std::vector<double> x_in = xd_ref;
     std::vector<double> x_out(3);

source/module_basis/module_ao/parallel_orbitals.cpp

@@ -247,7 +247,10 @@ int Parallel_Orbitals::set_nloc_wfc_Eij(
         }
         else
         {
-            ModuleBase::WARNING_QUIT("Parallel_Orbitals::set_nloc_wfc_Eij", "some processor has no bands-row-blocks.");
+            ModuleBase::WARNING_QUIT("Parallel_Orbitals::set_nloc_wfc_Eij",
+                "The number of columns of the 2D process grid exceeds the number of bands. "
+                "Try launching the calculation with fewer MPI processes."
+            );
         }
     }
     int col_b_bands = block / dim1;

source/module_elecstate/test/charge_mixing_test.cpp

@@ -10,7 +10,10 @@
 #include "module_basis/module_pw/pw_basis.h"
 #include "module_hamilt_general/module_xc/xc_functional.h"
 #undef private
+
+#ifdef _OPENMP
 #include <omp.h>
+#endif
 
 int FUNC_TYPE = 1;
 
@@ -116,7 +119,9 @@ class ChargeMixingTest : public ::testing::Test
 
 TEST_F(ChargeMixingTest, SetMixingTest)
 {
+#ifdef _OPENMP
     omp_set_num_threads(1);
+#endif
     PARAM.input.nspin = 1;
     Charge_Mixing CMtest;
     CMtest.set_rhopw(&pw_basis, &pw_basis);
@@ -214,7 +219,9 @@ TEST_F(ChargeMixingTest, SetMixingTest)
 
 TEST_F(ChargeMixingTest, InitMixingTest)
 {
+#ifdef _OPENMP
     omp_set_num_threads(1);
+#endif
     PARAM.input.nspin = 1;
     FUNC_TYPE = 1;
     Charge_Mixing CMtest;

source/module_esolver/esolver_ks_lcao.cpp

@@ -634,19 +634,22 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(const int istep, const int iter)
 
 #ifdef __EXX
     // calculate exact-exchange
-    if (GlobalC::exx_info.info_ri.real_number)
-    {
-        this->exd->exx_eachiterinit(istep,
-                                    *dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
-                                    this->kv,
-                                    iter);
-    }
-    else
+    if (PARAM.inp.calculation != "nscf")
     {
-        this->exc->exx_eachiterinit(istep,
-                                    *dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
-                                    this->kv,
-                                    iter);
+        if (GlobalC::exx_info.info_ri.real_number)
+        {
+            this->exd->exx_eachiterinit(istep,
+                                        *dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
+                                        this->kv,
+                                        iter);
+        }
+        else
+        {
+            this->exc->exx_eachiterinit(istep,
+                                        *dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
+                                        this->kv,
+                                        iter);
+        }
     }
 #endif
 
@@ -728,13 +731,16 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2density_single(int istep, int iter, double
 
     // 5) what's the exd used for?
 #ifdef __EXX
-    if (GlobalC::exx_info.info_ri.real_number)
+    if (PARAM.inp.calculation != "nscf")
     {
-        this->exd->exx_hamilt2density(*this->pelec, this->pv, iter);
-    }
-    else
-    {
-        this->exc->exx_hamilt2density(*this->pelec, this->pv, iter);
+        if (GlobalC::exx_info.info_ri.real_number)
+        {
+            this->exd->exx_hamilt2density(*this->pelec, this->pv, iter);
+        }
+        else
+        {
+            this->exc->exx_hamilt2density(*this->pelec, this->pv, iter);
+        }
     }
 #endif
 
@@ -921,11 +927,29 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(const int istep, int& iter)
 
 #ifdef __EXX
     // 3) save exx matrix
-    if (GlobalC::exx_info.info_global.cal_exx)
+    if (PARAM.inp.calculation != "nscf")
     {
-        GlobalC::exx_info.info_ri.real_number ?
-            this->exd->exx_iter_finish(this->kv, GlobalC::ucell, *this->p_hamilt, *this->pelec, *this->p_chgmix, this->scf_ene_thr, iter, istep, this->conv_esolver) :
-            this->exc->exx_iter_finish(this->kv, GlobalC::ucell, *this->p_hamilt, *this->pelec, *this->p_chgmix, this->scf_ene_thr, iter, istep, this->conv_esolver);
+        if (GlobalC::exx_info.info_global.cal_exx)
+        {
+            GlobalC::exx_info.info_ri.real_number ? this->exd->exx_iter_finish(this->kv,
+                                                                               GlobalC::ucell,
+                                                                               *this->p_hamilt,
+                                                                               *this->pelec,
+                                                                               *this->p_chgmix,
+                                                                               this->scf_ene_thr,
+                                                                               iter,
+                                                                               istep,
+                                                                               this->conv_esolver)
+                                                  : this->exc->exx_iter_finish(this->kv,
+                                                                               GlobalC::ucell,
+                                                                               *this->p_hamilt,
+                                                                               *this->pelec,
+                                                                               *this->p_chgmix,
+                                                                               this->scf_ene_thr,
+                                                                               iter,
+                                                                               istep,
+                                                                               this->conv_esolver);
+        }
     }
 #endif
 
@@ -1040,17 +1064,20 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
 
 #ifdef __EXX
     // 5) write Hexx matrix for NSCF (see `out_chg` in docs/advanced/input_files/input-main.md)
-    if (GlobalC::exx_info.info_global.cal_exx && PARAM.inp.out_chg[0]
-        && istep % PARAM.inp.out_interval == 0) // Peize Lin add if 2022.11.14
+    if (PARAM.inp.calculation != "nscf")
     {
-        const std::string file_name_exx = PARAM.globalv.global_out_dir + "HexxR" + std::to_string(GlobalV::MY_RANK);
-        if (GlobalC::exx_info.info_ri.real_number)
+        if (GlobalC::exx_info.info_global.cal_exx && PARAM.inp.out_chg[0]
+            && istep % PARAM.inp.out_interval == 0) // Peize Lin add if 2022.11.14
         {
-            ModuleIO::write_Hexxs_csr(file_name_exx, GlobalC::ucell, this->exd->get_Hexxs());
-        }
-        else
-        {
-            ModuleIO::write_Hexxs_csr(file_name_exx, GlobalC::ucell, this->exc->get_Hexxs());
+            const std::string file_name_exx = PARAM.globalv.global_out_dir + "HexxR" + std::to_string(GlobalV::MY_RANK);
+            if (GlobalC::exx_info.info_ri.real_number)
+            {
+                ModuleIO::write_Hexxs_csr(file_name_exx, GlobalC::ucell, this->exd->get_Hexxs());
+            }
+            else
+            {
+                ModuleIO::write_Hexxs_csr(file_name_exx, GlobalC::ucell, this->exc->get_Hexxs());
+            }
         }
     }
 #endif

source/module_esolver/esolver_ks_pw.cpp

@@ -422,19 +422,21 @@ void ESolver_KS_PW<T, Device>::update_pot(const int istep, const int iter)
 template <typename T, typename Device>
 void ESolver_KS_PW<T, Device>::iter_finish(const int istep, int& iter)
 {
-    // call iter_finish() of ESolver_KS
+    // 1) Call iter_finish() of ESolver_KS
     ESolver_KS<T, Device>::iter_finish(istep, iter);
 
-    // liuyu 2023-10-24
+    // 2) Update USPP-related quantities
     // D in uspp need vloc, thus needs update when veff updated
     // calculate the effective coefficient matrix for non-local pseudopotential
-    // projectors
+    // projectors 
+    // liuyu 2023-10-24
     if (PARAM.globalv.use_uspp)
     {
         ModuleBase::matrix veff = this->pelec->pot->get_effective_v();
         GlobalC::ppcell.cal_effective_D(veff, this->pw_rhod, GlobalC::ucell);
     }
 
+    // 3) Print out charge density
     if (this->out_freq_elec && iter % this->out_freq_elec == 0)
     {
         if (PARAM.inp.out_chg[0] > 0)
@@ -475,7 +477,8 @@ void ESolver_KS_PW<T, Device>::iter_finish(const int istep, int& iter)
                 }
             }
         }
-        // output wavefunctions
+        
+        // 4) Print out electronic wavefunctions
         if (this->wf.out_wfc_pw == 1 || this->wf.out_wfc_pw == 2)
         {
             std::stringstream ssw;
@@ -509,6 +512,7 @@ void ESolver_KS_PW<T, Device>::after_scf(const int istep)
         ModuleIO::write_wfc_pw(ssw.str(), this->psi[0], this->kv, this->pw_wfc);
     }
 
+    // 4) Transfer data from GPU to CPU
     if (this->device == base_device::GpuDevice)
     {
         castmem_2d_d2h_op()(this->psi[0].get_device(),
@@ -518,7 +522,7 @@ void ESolver_KS_PW<T, Device>::after_scf(const int istep)
                             this->psi[0].size());
     }
 
-    // Calculate band-decomposed (partial) charge density
+    // 5) Calculate band-decomposed (partial) charge density
     const std::vector<int> bands_to_print = PARAM.inp.bands_to_print;
     if (bands_to_print.size() > 0)
     {
@@ -545,7 +549,7 @@ void ESolver_KS_PW<T, Device>::after_scf(const int istep)
                               PARAM.inp.if_separate_k);
     }
 
-    //! 6) calculate Wannier functions
+    //! 6) Calculate Wannier functions
     if (PARAM.inp.calculation == "nscf" && PARAM.inp.towannier90)
     {
         std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Wannier functions calculation");
@@ -638,6 +642,7 @@ void ESolver_KS_PW<T, Device>::cal_stress(ModuleBase::matrix& stress)
 template <typename T, typename Device>
 void ESolver_KS_PW<T, Device>::after_all_runners()
 {
+    //! 1) Output information to screen
     GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl;
     GlobalV::ofs_running << std::setprecision(16);
     GlobalV::ofs_running << " !FINAL_ETOT_IS " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl;
@@ -681,10 +686,11 @@ void ESolver_KS_PW<T, Device>::after_all_runners()
     {
         nspin0 = 2;
     }
-    //! print occupation in istate.info
+    
+    //! 2) Print occupation numbers into istate.info
     ModuleIO::write_istate_info(this->pelec->ekb, this->pelec->wg, this->kv, &(GlobalC::Pkpoints));
 
-    //! compute density of states
+    //! 3) Compute density of states (DOS)
     if (PARAM.inp.out_dos)
     {
         ModuleIO::write_dos_pw(this->pelec->ekb,
@@ -707,7 +713,8 @@ void ESolver_KS_PW<T, Device>::after_all_runners()
         }
     }
 
-    if (PARAM.inp.out_band[0]) // pengfei 2014-10-13
+    //! 4) Print out band structure information
+    if (PARAM.inp.out_band[0])
     {
         for (int is = 0; is < nspin0; is++)
         {
@@ -725,12 +732,10 @@ void ESolver_KS_PW<T, Device>::after_all_runners()
         }
     }
 
-    if (PARAM.inp.basis_type == "pw" && winput::out_spillage) // xiaohui add 2013-09-01
+    //! 5) Calculate the spillage value, used to generate numerical atomic orbitals
+    if (PARAM.inp.basis_type == "pw" && winput::out_spillage)
     {
-        //  calculate spillage value.
-
-        // ! Print out overlap before spillage optimization to generate atomic
-        // orbitals
+        // ! Print out overlap matrices
         if (winput::out_spillage <= 2)
         {
             for (int i = 0; i < PARAM.inp.bessel_nao_rcuts.size(); i++)
@@ -747,13 +752,13 @@ void ESolver_KS_PW<T, Device>::after_all_runners()
         }
     }
 
-    //! Print out wave functions in real space
+    //! 6) Print out electronic wave functions in real space
     if (this->wf.out_wfc_r == 1) // Peize Lin add 2021.11.21
     {
         ModuleIO::write_psi_r_1(this->psi[0], this->pw_wfc, "wfc_realspace", true, this->kv);
     }
 
-    //! Use Kubo-Greenwood method to compute conductivities
+    //! 7) Use Kubo-Greenwood method to compute conductivities
     if (PARAM.inp.cal_cond)
     {
         EleCond elec_cond(&GlobalC::ucell, &this->kv, this->pelec, this->pw_wfc, this->psi, &GlobalC::ppcell);

source/module_esolver/lcao_before_scf.cpp

@@ -210,13 +210,16 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(const int istep)
 
     // Peize Lin add 2016-12-03
 #ifdef __EXX // set xc type before the first cal of xc in pelec->init_scf
-    if (GlobalC::exx_info.info_ri.real_number)
+    if (PARAM.inp.calculation != "nscf")
     {
-        this->exd->exx_beforescf(istep, this->kv, *this->p_chgmix, GlobalC::ucell, orb_);
-    }
-    else
-    {
-        this->exc->exx_beforescf(istep, this->kv, *this->p_chgmix, GlobalC::ucell, orb_);
+        if (GlobalC::exx_info.info_ri.real_number)
+        {
+            this->exd->exx_beforescf(istep, this->kv, *this->p_chgmix, GlobalC::ucell, orb_);
+        }
+        else
+        {
+            this->exc->exx_beforescf(istep, this->kv, *this->p_chgmix, GlobalC::ucell, orb_);
+        }
     }
 #endif // __EXX