Merge branch 'develop' into rdmft_PR

Author: WHUweiqingzhou

source/Makefile.Objects

@@ -247,13 +247,10 @@ OBJS_ESOLVER=esolver.o\
     esolver_of.o\
     esolver_of_tool.o\
     esolver_of_interface.o\
-    pw_init_globalc.o\
     pw_others.o\
 
 OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\
       esolver_ks_lcao_tddft.o\
-      dpks_cal_e_delta_band.o\
-      set_matrix_grid.o\
       lcao_before_scf.o\
       esolver_gets.o\
       lcao_others.o\

source/module_esolver/CMakeLists.txt

@@ -10,15 +10,12 @@ list(APPEND objects
     esolver_of.cpp
     esolver_of_interface.cpp
     esolver_of_tool.cpp
-    pw_init_globalc.cpp
     pw_others.cpp
 )
 if(ENABLE_LCAO)
   list(APPEND objects
       esolver_ks_lcao.cpp
       esolver_ks_lcao_tddft.cpp
-      dpks_cal_e_delta_band.cpp
-      set_matrix_grid.cpp
       lcao_before_scf.cpp
       esolver_gets.cpp
       lcao_others.cpp

source/module_esolver/dpks_cal_e_delta_band.cpp

@@ -67,14 +67,6 @@ ESolver_KS<T, Device>::ESolver_KS()
     ///----------------------------------------------------------
     p_chgmix = new Charge_Mixing();
     p_chgmix->set_rhopw(this->pw_rho, this->pw_rhod);
-
-    ///----------------------------------------------------------
-    /// wavefunc
-    ///----------------------------------------------------------
-    this->wf.init_wfc = PARAM.inp.init_wfc;
-    this->wf.mem_saver = PARAM.inp.mem_saver;
-    this->wf.out_wfc_pw = PARAM.inp.out_wfc_pw;
-    this->wf.out_wfc_r = PARAM.inp.out_wfc_r;
 }
 
 //------------------------------------------------------------------------------

source/module_esolver/esolver_ks.cpp

@@ -560,7 +560,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(const int istep, const int iter)
     // mohan move it outside 2011-01-13
     // first need to calculate the weight according to
     // electrons number.
-    if (istep == 0 && this->wf.init_wfc == "file")
+    if (istep == 0 && PARAM.inp.init_wfc == "file")
     {
         if (iter == 1)
         {
@@ -892,7 +892,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(const int istep, int& iter)
         const std::vector<std::vector<TK>>& dm
             = dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM()->get_DMK_vector();
 
-        this->dpks_cal_e_delta_band(dm);
+        GlobalC::ld.dpks_cal_e_delta_band(dm, this->kv.get_nks());
     }
 #endif
 

source/module_esolver/esolver_ks_lcao.cpp

@@ -97,13 +97,6 @@ class ESolver_KS_LCAO : public ESolver_KS<TK> {
     std::shared_ptr<Exx_LRI<std::complex<double>>> exx_lri_complex = nullptr;
 #endif
 
-  private:
-#ifdef __DEEPKS
-    void dpks_cal_e_delta_band(const std::vector<std::vector<TK>>& dm) const;
-
-    void dpks_cal_projected_DM(
-        const elecstate::DensityMatrix<TK, double>* dm) const;
-#endif
     friend class LR::ESolver_LR<double, double>;
     friend class LR::ESolver_LR<std::complex<double>, double>;
 };

source/module_esolver/esolver_ks_lcao.h

@@ -76,7 +76,7 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(const Input_para& inp, UnitCell&
 
 void ESolver_KS_LCAO_TDDFT::hamilt2density_single(const int istep, const int iter, const double ethr)
 {
-    if (wf.init_wfc == "file")
+    if (PARAM.inp.init_wfc == "file")
     {
         if (istep >= 1)
         {
@@ -256,7 +256,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter)
     const int nlocal = PARAM.globalv.nlocal;
 
     // store wfc and Hk laststep
-    if (istep >= (wf.init_wfc == "file" ? 0 : 1) && this->conv_esolver)
+    if (istep >= (PARAM.inp.init_wfc == "file" ? 0 : 1) && this->conv_esolver)
     {
         if (this->psi_laststep == nullptr)
         {
@@ -311,7 +311,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter)
         }
 
         // calculate energy density matrix for tddft
-        if (istep >= (wf.init_wfc == "file" ? 0 : 2) && module_tddft::Evolve_elec::td_edm == 0)
+        if (istep >= (PARAM.inp.init_wfc == "file" ? 0 : 2) && module_tddft::Evolve_elec::td_edm == 0)
         {
             elecstate::cal_edm_tddft(this->pv, this->pelec, this->kv, this->p_hamilt);
         }

source/module_esolver/esolver_ks_lcao_tddft.cpp

@@ -18,6 +18,7 @@
 #include "module_base/memory.h"
 #include "module_base/module_device/device.h"
 #include "module_elecstate/elecstate_pw.h"
+#include "module_elecstate/elecstate_pw_sdft.h"
 #include "module_hamilt_general/module_vdw/vdw.h"
 #include "module_hamilt_pw/hamilt_pwdft/elecond.h"
 #include "module_hamilt_pw/hamilt_pwdft/hamilt_pw.h"
@@ -136,14 +137,29 @@ void ESolver_KS_PW<T, Device>::before_all_runners(const Input_para& inp, UnitCel
     // 3) initialize ElecState,
     if (this->pelec == nullptr)
     {
-        this->pelec = new elecstate::ElecStatePW<T, Device>(this->pw_wfc,
-                                                            &(this->chr),
-                                                            &(this->kv),
-                                                            &ucell,
-                                                            &GlobalC::ppcell,
-                                                            this->pw_rhod,
-                                                            this->pw_rho,
-                                                            this->pw_big);
+        if (inp.esolver_type == "sdft")
+        {
+            //! SDFT only supports double precision currently
+            this->pelec = new elecstate::ElecStatePW_SDFT<std::complex<double>, Device>(this->pw_wfc,
+                                                                                        &(this->chr),
+                                                                                        &(this->kv),
+                                                                                        &ucell,
+                                                                                        &(GlobalC::ppcell),
+                                                                                        this->pw_rhod,
+                                                                                        this->pw_rho,
+                                                                                        this->pw_big);
+        }
+        else
+        {
+            this->pelec = new elecstate::ElecStatePW<T, Device>(this->pw_wfc,
+                                                                &(this->chr),
+                                                                &(this->kv),
+                                                                &ucell,
+                                                                &GlobalC::ppcell,
+                                                                this->pw_rhod,
+                                                                this->pw_rho,
+                                                                this->pw_big);
+        }
     }
 
     //! 4) inititlize the charge density.
@@ -165,12 +181,12 @@ void ESolver_KS_PW<T, Device>::before_all_runners(const Input_para& inp, UnitCel
     }
 
     //! 7) prepare some parameters for electronic wave functions initilization
-    this->p_wf_init = new psi::WFInit<T, Device>(PARAM.inp.init_wfc,
-                                                 PARAM.inp.ks_solver,
-                                                 PARAM.inp.basis_type,
-                                                 PARAM.inp.psi_initializer,
-                                                 &this->wf,
-                                                 this->pw_wfc);
+    this->p_wf_init = new psi::PSIInit<T, Device>(PARAM.inp.init_wfc,
+                                                  PARAM.inp.ks_solver,
+                                                  PARAM.inp.basis_type,
+                                                  PARAM.inp.psi_initializer,
+                                                  &this->wf,
+                                                  this->pw_wfc);
     this->p_wf_init->prepare_init(&(this->sf),
                                   &ucell,
                                   1,
@@ -180,8 +196,39 @@ void ESolver_KS_PW<T, Device>::before_all_runners(const Input_para& inp, UnitCel
 #endif
                                   &GlobalC::ppcell);
 
-    //! 8) setup global classes
-    this->Init_GlobalC(inp, ucell, GlobalC::ppcell);
+    if (this->psi != nullptr)
+    {
+        delete this->psi;
+    }
+
+    //! init pseudopotential
+    GlobalC::ppcell.init(ucell.ntype, &this->sf, this->pw_wfc);
+
+    //! initalize local pseudopotential
+    GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, this->pw_rhod);
+    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
+
+    //! Initalize non-local pseudopotential
+    GlobalC::ppcell.init_vnl(ucell, this->pw_rhod);
+    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "NON-LOCAL POTENTIAL");
+
+    //! Allocate psi
+    this->p_wf_init->allocate_psi(this->psi,
+                                  this->kv.get_nkstot(),
+                                  this->kv.get_nks(),
+                                  this->kv.ngk.data(),
+                                  this->pw_wfc->npwk_max,
+                                  &(this->sf));
+
+    this->kspw_psi = PARAM.inp.device == "gpu" || PARAM.inp.precision == "single"
+                         ? new psi::Psi<T, Device>(this->psi[0])
+                         : reinterpret_cast<psi::Psi<T, Device>*>(this->psi);
+
+    if (PARAM.inp.precision == "single")
+    {
+        ModuleBase::Memory::record("Psi_single", sizeof(T) * this->psi[0].size());
+    }
+    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT BASIS");
 
     //! 9) setup occupations
     if (PARAM.inp.ocp)
@@ -500,7 +547,7 @@ void ESolver_KS_PW<T, Device>::iter_finish(const int istep, int& iter)
         }
 
         // 4) Print out electronic wavefunctions
-        if (this->wf.out_wfc_pw == 1 || this->wf.out_wfc_pw == 2)
+        if (PARAM.inp.out_wfc_pw == 1 || PARAM.inp.out_wfc_pw == 2)
         {
             std::stringstream ssw;
             ssw << PARAM.globalv.global_out_dir << "WAVEFUNC";
@@ -526,7 +573,7 @@ void ESolver_KS_PW<T, Device>::after_scf(const int istep)
     ESolver_KS<T, Device>::after_scf(istep);
 
     // 3) output wavefunctions
-    if (this->wf.out_wfc_pw == 1 || this->wf.out_wfc_pw == 2)
+    if (PARAM.inp.out_wfc_pw == 1 || PARAM.inp.out_wfc_pw == 2)
     {
         std::stringstream ssw;
         ssw << PARAM.globalv.global_out_dir << "WAVEFUNC";
@@ -774,7 +821,7 @@ void ESolver_KS_PW<T, Device>::after_all_runners()
     }
 
     //! 6) Print out electronic wave functions in real space
-    if (this->wf.out_wfc_r == 1) // Peize Lin add 2021.11.21
+    if (PARAM.inp.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);
     }

source/module_esolver/esolver_ks_pw.cpp

@@ -46,18 +46,14 @@ class ESolver_KS_PW : public ESolver_KS<T, Device>
 
     virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override;
 
-    // temporary, this will be removed in the future;
-    // Init Global class
-    void Init_GlobalC(const Input_para& inp, UnitCell& ucell, pseudopot_cell_vnl& ppcell);
-
     virtual void allocate_hamilt();
     virtual void deallocate_hamilt();
 
     //! hide the psi in ESolver_KS for tmp use
     psi::Psi<std::complex<double>, base_device::DEVICE_CPU>* psi = nullptr;
 
     // psi_initializer controller
-    psi::WFInit<T, Device>* p_wf_init = nullptr;
+    psi::PSIInit<T, Device>* p_wf_init = nullptr;
 
     Device* ctx = {};
 

source/module_esolver/esolver_ks_pw.h

@@ -50,51 +50,9 @@ void ESolver_SDFT_PW<T, Device>::before_all_runners(const Input_para& inp, UnitC
     this->method_sto = inp.method_sto;
 
     // 2) run "before_all_runners" in ESolver_KS
-    ESolver_KS<T, Device>::before_all_runners(inp, ucell);
-
-    // 3) initialize the pointer for electronic states of SDFT
-    this->pelec = new elecstate::ElecStatePW_SDFT<T, Device>(this->pw_wfc,
-                                                             &(this->chr),
-                                                             &(this->kv),
-                                                             &ucell,
-                                                             &(GlobalC::ppcell),
-                                                             this->pw_rhod,
-                                                             this->pw_rho,
-                                                             this->pw_big);
-
-    // 4) inititlize the charge density.
-    this->pelec->charge->allocate(PARAM.inp.nspin);
-    this->pelec->omega = ucell.omega;
-
-    // 5) initialize the potential.
-    if (this->pelec->pot == nullptr)
-    {
-        this->pelec->pot = new elecstate::Potential(this->pw_rhod,
-                                                    this->pw_rho,
-                                                    &ucell,
-                                                    &(GlobalC::ppcell.vloc),
-                                                    &(this->sf),
-                                                    &(this->pelec->f_en.etxc),
-                                                    &(this->pelec->f_en.vtxc));
-    }
-
-    // 6) prepare some parameters for electronic wave functions initilization
-    this->p_wf_init = new psi::WFInit<T, Device>(PARAM.inp.init_wfc,
-                                                 PARAM.inp.ks_solver,
-                                                 PARAM.inp.basis_type,
-                                                 PARAM.inp.psi_initializer,
-                                                 &this->wf,
-                                                 this->pw_wfc);
-    // 7) set occupatio, redundant?
-    if (PARAM.inp.ocp)
-    {
-        this->pelec->fixed_weights(PARAM.inp.ocp_kb, PARAM.inp.nbands, PARAM.inp.nelec);
-    }
-
-    // 8) initialize the global classes
-    this->Init_GlobalC(inp, ucell, GlobalC::ppcell); // temporary
+    ESolver_KS_PW<T, Device>::before_all_runners(inp, ucell);
 
-    // 9) initialize the stochastic wave functions
+    // 3) initialize the stochastic wave functions
     this->stowf.init(&this->kv, this->pw_wfc->npwk_max);
     if (inp.nbands_sto != 0)
     {
@@ -117,7 +75,7 @@ void ESolver_SDFT_PW<T, Device>::before_all_runners(const Input_para& inp, UnitC
     }
     this->stowf.sync_chi0();
 
-    // 10) allocate spaces for \sqrt(f(H))|chi> and |\tilde{chi}>
+    // 4) allocate spaces for \sqrt(f(H))|chi> and |\tilde{chi}>
     size_t size = stowf.chi0->size();
     this->stowf.shchi
         = new psi::Psi<T, Device>(this->kv.get_nks(), this->stowf.nchip_max, this->wf.npwx, this->kv.ngk.data());