Refactor: remove nscf() (#5455)

* Refactor: remove nscf()

* move calculate_weight() out of psi2rho()

* update PW_DIAG_NMAX for nscf

Author: YuLiu98

source/Makefile.Objects

@@ -245,10 +245,8 @@ OBJS_ESOLVER=esolver.o\
     esolver_of.o\
     esolver_of_tool.o\
     esolver_of_interface.o\
-    pw_fun.o\
     pw_init_after_vc.o\
     pw_init_globalc.o\
-    pw_nscf.o\
     pw_others.o\
 
 OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\
@@ -259,7 +257,6 @@ OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\
       set_matrix_grid.o\
       lcao_before_scf.o\
       lcao_gets.o\
-      lcao_nscf.o\
       lcao_others.o\
       lcao_init_after_vc.o\
       lcao_fun.o\

source/driver_run.cpp

@@ -62,15 +62,14 @@ void Driver::driver_run() {
     {
         Run_MD::md_line(GlobalC::ucell, p_esolver, PARAM);
     }
-    else if (cal_type == "scf" || cal_type == "relax" || cal_type == "cell-relax")
+    else if (cal_type == "scf" || cal_type == "relax" || cal_type == "cell-relax" || cal_type == "nscf")
     {
         Relax_Driver rl_driver;
         rl_driver.relax_driver(p_esolver);
     }
     else
     {
         //! supported "other" functions:
-        //! nscf(PW,LCAO),
         //! get_pchg(LCAO),
         //! test_memory(PW,LCAO),
         //! test_neighbour(LCAO),

source/module_elecstate/elecstate_lcao.cpp

@@ -1,7 +1,5 @@
 #include "elecstate_lcao.h"
 
-#include <vector>
-
 #include "cal_dm.h"
 #include "module_base/timer.h"
 #include "module_elecstate/module_dm/cal_dm_psi.h"
@@ -11,6 +9,8 @@
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_parameter/parameter.h"
 
+#include <vector>
+
 namespace elecstate
 {
 
@@ -21,34 +21,31 @@ void ElecStateLCAO<std::complex<double>>::psiToRho(const psi::Psi<std::complex<d
     ModuleBase::TITLE("ElecStateLCAO", "psiToRho");
     ModuleBase::timer::tick("ElecStateLCAO", "psiToRho");
 
-    this->calculate_weights();
-
-    // the calculations of dm, and dm -> rho are, technically, two separate
-    // functionalities, as we cannot rule out the possibility that we may have a
-    // dm from other sources, such as read from file. However, since we are not
-    // separating them now, I opt to add a flag to control how dm is obtained as
-    // of now
-    if (!PARAM.inp.dm_to_rho)
-    {
-        this->calEBand();
-
-        ModuleBase::GlobalFunc::NOTE("Calculate the density matrix.");
-
-        // this part for calculating DMK in 2d-block format, not used for charge
-        // now
-        //    psi::Psi<std::complex<double>> dm_k_2d();
-
-        if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == "scalapack_gvx" || PARAM.inp.ks_solver == "lapack"
-            || PARAM.inp.ks_solver == "cusolver" || PARAM.inp.ks_solver == "cusolvermp"
-            || PARAM.inp.ks_solver == "cg_in_lcao") // Peize Lin test 2019-05-15
-        {
-            elecstate::cal_dm_psi(this->DM->get_paraV_pointer(),
-                                  this->wg,
-                                  psi,
-                                  *(this->DM));
-            this->DM->cal_DMR();
-        }
-    }
+    // // the calculations of dm, and dm -> rho are, technically, two separate
+    // // functionalities, as we cannot rule out the possibility that we may have a
+    // // dm from other sources, such as read from file. However, since we are not
+    // // separating them now, I opt to add a flag to control how dm is obtained as
+    // // of now
+    // if (!PARAM.inp.dm_to_rho)
+    // {
+    //     ModuleBase::GlobalFunc::NOTE("Calculate the density matrix.");
+
+    //     // this part for calculating DMK in 2d-block format, not used for charge
+    //     // now
+    //     //    psi::Psi<std::complex<double>> dm_k_2d();
+
+    //     if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver ==
+    //     "scalapack_gvx" || PARAM.inp.ks_solver == "lapack"
+    //         || PARAM.inp.ks_solver == "cusolver" || PARAM.inp.ks_solver == "cusolvermp"
+    //         || PARAM.inp.ks_solver == "cg_in_lcao") // Peize Lin test 2019-05-15
+    //     {
+    //         elecstate::cal_dm_psi(this->DM->get_paraV_pointer(),
+    //                               this->wg,
+    //                               psi,
+    //                               *(this->DM));
+    //         this->DM->cal_DMR();
+    //     }
+    // }
 
     for (int is = 0; is < PARAM.inp.nspin; is++)
     {
@@ -83,23 +80,6 @@ void ElecStateLCAO<double>::psiToRho(const psi::Psi<double>& psi)
     ModuleBase::TITLE("ElecStateLCAO", "psiToRho");
     ModuleBase::timer::tick("ElecStateLCAO", "psiToRho");
 
-    this->calculate_weights();
-    this->calEBand();
-
-    if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == "scalapack_gvx" || PARAM.inp.ks_solver == "lapack"
-        || PARAM.inp.ks_solver == "cusolver" || PARAM.inp.ks_solver == "cusolvermp" || PARAM.inp.ks_solver == "cg_in_lcao")
-    {
-        ModuleBase::timer::tick("ElecStateLCAO", "cal_dm_2d");
-
-        // get DMK in 2d-block format
-        elecstate::cal_dm_psi(this->DM->get_paraV_pointer(),
-                              this->wg,
-                              psi,
-                              *(this->DM));
-        this->DM->cal_DMR();
-        ModuleBase::timer::tick("ElecStateLCAO", "cal_dm_2d");
-    }
-
     for (int is = 0; is < PARAM.inp.nspin; is++)
     {
         ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is],

source/module_elecstate/elecstate_lcao.h

@@ -74,6 +74,8 @@ class ElecStateLCAO : public ElecState
     void dmToRho(std::vector<TK*> pexsi_DM, std::vector<TK*> pexsi_EDM);
 #endif
 
+    DensityMatrix<TK, double>* DM = nullptr;
+
   protected:
     // calculate electronic charge density on grid points or density matrix in real space
     // the consequence charge density rho saved into rho_out, preparing for charge mixing.
@@ -85,7 +87,6 @@ class ElecStateLCAO : public ElecState
 
     Gint_Gamma* gint_gamma = nullptr; // mohan add 2024-04-01
     Gint_k* gint_k = nullptr;         // mohan add 2024-04-01
-    DensityMatrix<TK, double>* DM = nullptr;
 };
 
 template <typename TK>

source/module_elecstate/elecstate_pw.cpp

@@ -89,9 +89,6 @@ void ElecStatePW<T, Device>::psiToRho(const psi::Psi<T, Device>& psi)
     ModuleBase::timer::tick("ElecStatePW", "psiToRho");
 
     this->init_rho_data();
-    this->calculate_weights();
-
-    this->calEBand();
 
     for(int is=0; is<PARAM.inp.nspin; is++)
 	{

source/module_elecstate/elecstate_pw_sdft.cpp

@@ -17,8 +17,6 @@ void ElecStatePW_SDFT<T, Device>::psiToRho(const psi::Psi<T, Device>& psi)
 
     if (GlobalV::MY_STOGROUP == 0)
     {
-        this->calEBand();
-
         for (int is = 0; is < nspin; is++)
         {
             setmem_var_op()(this->ctx, this->rho[is], 0, this->charge->nrxx);

source/module_esolver/CMakeLists.txt

@@ -10,10 +10,8 @@ list(APPEND objects
     esolver_of.cpp
     esolver_of_interface.cpp
     esolver_of_tool.cpp
-    pw_fun.cpp
     pw_init_after_vc.cpp
     pw_init_globalc.cpp
-    pw_nscf.cpp
     pw_others.cpp
 )
 if(ENABLE_LCAO)
@@ -26,7 +24,6 @@ if(ENABLE_LCAO)
       dftu_cal_occup_m.cpp
       lcao_before_scf.cpp
       lcao_gets.cpp
-      lcao_nscf.cpp
       lcao_others.cpp
       lcao_init_after_vc.cpp
       lcao_fun.cpp

source/module_esolver/esolver_fp.h

@@ -52,9 +52,6 @@ namespace ModuleESolver
         //! Electorn charge density
         Charge chr;
 
-        //! Non-Self-Consistant Filed (NSCF) calculations
-        virtual void nscf(){};
-
         //! Structure factors that used with plane-wave basis set
         Structure_Factor sf;
 

source/module_esolver/esolver_ks.cpp

@@ -402,7 +402,7 @@ void ESolver_KS<T, Device>::hamilt2density(const int istep, const int iter, cons
 
         drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec);
         hsolver_error = 0.0;
-        if (iter == 1)
+        if (iter == 1 && PARAM.inp.calculation != "nscf")
         {
             hsolver_error
                 = hsolver::cal_hsolve_error(PARAM.inp.basis_type, PARAM.inp.esolver_type, diag_ethr, PARAM.inp.nelec);
@@ -472,13 +472,10 @@ void ESolver_KS<T, Device>::runner(const int istep, UnitCell& ucell)
 
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF");
 
+    // 4) SCF iterations
     this->conv_esolver = false;
     this->niter = this->maxniter;
-
-    // 4) SCF iterations
     this->diag_ethr = PARAM.inp.pw_diag_thr;
-
-    std::cout << " * * * * * *\n << Start SCF iteration." << std::endl;
     for (int iter = 1; iter <= this->maxniter; ++iter)
     {
         // 6) initialization of SCF iterations
@@ -500,7 +497,6 @@ void ESolver_KS<T, Device>::runner(const int istep, UnitCell& ucell)
             break;
         }
     } // end scf iterations
-    std::cout << " >> Leave SCF iteration.\n * * * * * *" << std::endl;
 
     // 15) after scf
     ModuleBase::timer::tick(this->classname, "after_scf");
@@ -632,7 +628,7 @@ void ESolver_KS<T, Device>::iter_finish(const int istep, int& iter)
 
         // If drho < hsolver_error in the first iter or drho < scf_thr, we
         // do not change rho.
-        if (drho < hsolver_error || this->conv_esolver)
+        if (drho < hsolver_error || this->conv_esolver || PARAM.inp.calculation == "nscf")
         {
             if (drho < hsolver_error)
             {

source/module_esolver/esolver_ks.h

@@ -49,9 +49,6 @@ class ESolver_KS : public ESolver_FP
     // calculate electron density from a specific Hamiltonian with ethr
     virtual void hamilt2density_single(const int istep, const int iter, const double ethr);
 
-    // calculate electron states from a specific Hamiltonian
-    virtual void hamilt2estates(const double ethr) {};
-
     // calculate electron density from a specific Hamiltonian
     void hamilt2density(const int istep, const int iter, const double ethr);