Merge branch 'deepmodeling:develop' into TDDFT_GPU_phase_1

Author: AsTonyshment

source/module_cell/klist.cpp

@@ -148,12 +148,12 @@ void K_Vectors::set(const UnitCell& ucell,
     // It's very important in parallel case,
     // firstly do the mpi_k() and then
     // do set_kup_and_kdw()
-    GlobalC::Pkpoints.kinfo(nkstot,
-                            GlobalV::KPAR,
-                            GlobalV::MY_POOL,
-                            GlobalV::RANK_IN_POOL,
-                            GlobalV::NPROC,
-                            nspin_in); // assign k points to several process pools
+    this->para_k.kinfo(nkstot,
+                       GlobalV::KPAR,
+                       GlobalV::MY_POOL,
+                       GlobalV::RANK_IN_POOL,
+                       GlobalV::NPROC,
+                       nspin_in); // assign k points to several process pools
 #ifdef __MPI
     // distribute K point data to the corresponding process
     this->mpi_k(); // 2008-4-29
@@ -1163,7 +1163,7 @@ void K_Vectors::mpi_k()
 
     Parallel_Common::bcast_double(koffset, 3);
 
-    this->nks = GlobalC::Pkpoints.nks_pool[GlobalV::MY_POOL];
+    this->nks = this->para_k.nks_pool[GlobalV::MY_POOL];
 
     GlobalV::ofs_running << std::endl;
     ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "k-point number in this process", nks);
@@ -1217,7 +1217,7 @@ void K_Vectors::mpi_k()
     for (int i = 0; i < nks; i++)
     {
         // 3 is because each k point has three value:kx, ky, kz
-        k_index = i + GlobalC::Pkpoints.startk_pool[GlobalV::MY_POOL];
+        k_index = i + this->para_k.startk_pool[GlobalV::MY_POOL];
         kvec_c[i].x = kvec_c_aux[k_index * 3];
         kvec_c[i].y = kvec_c_aux[k_index * 3 + 1];
         kvec_c[i].z = kvec_c_aux[k_index * 3 + 2];

source/module_cell/klist.h

@@ -5,7 +5,7 @@
 #include "module_base/global_variable.h"
 #include "module_base/matrix3.h"
 #include "module_cell/unitcell.h"
-
+#include "parallel_kpoints.h"
 #include <vector>
 
 class K_Vectors
@@ -31,6 +31,9 @@ class K_Vectors
     K_Vectors& operator=(const K_Vectors&) = default;
     K_Vectors& operator=(K_Vectors&& rhs) = default;
 
+    Parallel_Kpoints para_k; ///< parallel for kpoints
+
+
     /**
      * @brief Set up the k-points for the system.
      *

source/module_cell/parallel_kpoints.cpp

@@ -3,14 +3,6 @@
 #include "module_base/parallel_common.h"
 #include "module_base/parallel_global.h"
 
-Parallel_Kpoints::Parallel_Kpoints()
-{
-}
-
-Parallel_Kpoints::~Parallel_Kpoints()
-{
-}
-
 // the kpoints here are reduced after symmetry applied.
 void Parallel_Kpoints::kinfo(int& nkstot_in,
                              const int& kpar_in,
@@ -227,7 +219,7 @@ void Parallel_Kpoints::pool_collection(double* value_re,
     return;
 }
 
-void Parallel_Kpoints::pool_collection(std::complex<double>* value, const ModuleBase::ComplexArray& w, const int& ik)
+void Parallel_Kpoints::pool_collection(std::complex<double>* value, const ModuleBase::ComplexArray& w, const int& ik) const
 {
     const int dim2 = w.getBound2();
     const int dim3 = w.getBound3();
@@ -237,7 +229,7 @@ void Parallel_Kpoints::pool_collection(std::complex<double>* value, const Module
 }
 
 template <class T, class V>
-void Parallel_Kpoints::pool_collection_aux(T* value, const V& w, const int& dim, const int& ik)
+void Parallel_Kpoints::pool_collection_aux(T* value, const V& w, const int& dim, const int& ik) const
 {
 #ifdef __MPI
     const int ik_now = ik - this->startk_pool[this->my_pool];

source/module_cell/parallel_kpoints.h

@@ -9,8 +9,8 @@
 class Parallel_Kpoints
 {
   public:
-    Parallel_Kpoints();
-    ~Parallel_Kpoints();
+    Parallel_Kpoints(){};
+    ~Parallel_Kpoints(){};
 
     void kinfo(int& nkstot_in,
                const int& kpar_in,
@@ -28,9 +28,9 @@ class Parallel_Kpoints
                          const ModuleBase::realArray& a,
                          const ModuleBase::realArray& b,
                          const int& ik);
-    void pool_collection(std::complex<double>* value, const ModuleBase::ComplexArray& w, const int& ik);
+    void pool_collection(std::complex<double>* value, const ModuleBase::ComplexArray& w, const int& ik) const;
     template <class T, class V>
-    void pool_collection_aux(T* value, const V& w, const int& dim, const int& ik);
+    void pool_collection_aux(T* value, const V& w, const int& dim, const int& ik) const;
 #ifdef __MPI
     /**
      * @brief gather kpoints from all processors
@@ -46,8 +46,8 @@ class Parallel_Kpoints
     // int* nproc_pool = nullptr;    it is not used
 
     // inforamation about kpoints, dim: KPAR
-    std::vector<int> nks_pool;    // number of k-points in each pool
-    std::vector<int> startk_pool; // the first k-point in each pool
+    std::vector<int> nks_pool;    // number of k-points in each pool, here use k-points without spin
+    std::vector<int> startk_pool; // the first k-point in each pool, here use k-points without spin
 
     // information about which pool each k-point belongs to,
     std::vector<int> whichpool; // whichpool[k] : the pool which k belongs to, dim: nkstot_np

source/module_elecstate/elecstate.h

@@ -1,10 +1,9 @@
 #ifndef ELECSTATE_H
 #define ELECSTATE_H
-#include "module_parameter/parameter.h"
-
 #include "fp_energy.h"
 #include "module_cell/klist.h"
 #include "module_elecstate/module_charge/charge.h"
+#include "module_parameter/parameter.h"
 #include "module_psi/psi.h"
 #include "potentials/potential_new.h"
 
@@ -14,10 +13,10 @@ namespace elecstate
 class ElecState
 {
   public:
-    ElecState(){}
-    ElecState(Charge* charge_in,
-              ModulePW::PW_Basis* rhopw_in,
-              ModulePW::PW_Basis_Big* bigpw_in)
+    ElecState()
+    {
+    }
+    ElecState(Charge* charge_in, ModulePW::PW_Basis* rhopw_in, ModulePW::PW_Basis_Big* bigpw_in)
     {
         this->charge = charge_in;
         this->charge->set_rhopw(rhopw_in);
@@ -26,20 +25,20 @@ class ElecState
     }
     virtual ~ElecState()
     {
-        if(this->pot != nullptr) 
+        if (this->pot != nullptr)
         {
             delete this->pot;
             this->pot = nullptr;
         }
     }
-    void init_ks(Charge *chg_in, // pointer for class Charge
-                      const K_Vectors *klist_in,
-                      int nk_in, // number of k points
-                      ModulePW::PW_Basis* rhopw_in,
-                      const ModulePW::PW_Basis_Big* bigpw_in); 
+    void init_ks(Charge* chg_in, // pointer for class Charge
+                 const K_Vectors* klist_in,
+                 int nk_in, // number of k points
+                 ModulePW::PW_Basis* rhopw_in,
+                 const ModulePW::PW_Basis_Big* bigpw_in);
 
     // return current electronic density rho, as a input for constructing Hamiltonian
-    virtual const double *getRho(int spin) const;
+    virtual const double* getRho(int spin) const;
 
     // 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.
@@ -78,17 +77,14 @@ class ElecState
 
     // use occupied weights from INPUT and skip calculate_weights
     // mohan updated on 2024-06-08
-	void fixed_weights(
-			const std::vector<double>& ocp_kb,
-			const int &nbands,
-			const double &nelec);
+    void fixed_weights(const std::vector<double>& ocp_kb, const int& nbands, const double& nelec);
 
-    // if nupdown is not 0(TWO_EFERMI case), 
-    // nelec_spin will be fixed and weights will be constrained 
+    // if nupdown is not 0(TWO_EFERMI case),
+    // nelec_spin will be fixed and weights will be constrained
     void init_nelec_spin();
-    //used to record number of electrons per spin index
-    //for NSPIN=2, it will record number of spin up and number of spin down
-    //for NSPIN=4, it will record total number, magnetization for x, y, z direction  
+    // used to record number of electrons per spin index
+    // for NSPIN=2, it will record number of spin up and number of spin down
+    // for NSPIN=4, it will record total number, magnetization for x, y, z direction
     std::vector<double> nelec_spin;
 
     virtual void print_psi(const psi::Psi<double>& psi_in, const int istep = -1)
@@ -102,7 +98,7 @@ class ElecState
 
     /**
      * @brief Init rho_core, init rho, renormalize rho, init pot
-     * 
+     *
      * @param istep i-th step
      * @param ucell unit cell
      * @param strucfac structure factor
@@ -142,28 +138,24 @@ class ElecState
     void set_exx(const std::complex<double>& Eexx);
 #endif //__LCAO
 #endif //__EXX
- 
+
     double get_hartree_energy();
     double get_etot_efield();
     double get_etot_gatefield();
 
     double get_solvent_model_Ael();
     double get_solvent_model_Acav();
 
-    virtual double get_spin_constrain_energy() {
+    virtual double get_spin_constrain_energy()
+    {
         return 0.0;
     }
 
     double get_dftu_energy();
     double get_local_pp_energy();
 
-#ifdef __DEEPKS
-    double get_deepks_E_delta();
-    double get_deepks_E_delta_band();
-#endif
-
-    fenergy f_en;                                  ///< energies contribute to the total free energy
-    efermi eferm;                                  ///< fermi energies
+    fenergy f_en; ///< energies contribute to the total free energy
+    efermi eferm; ///< fermi energies
 
     // below defines the bandgap:
 

source/module_elecstate/elecstate_energy.cpp

@@ -1,10 +1,10 @@
-#include <cmath>
-
 #include "elecstate.h"
 #include "elecstate_getters.h"
 #include "module_base/global_variable.h"
 #include "module_base/parallel_reduce.h"
 #include "module_parameter/parameter.h"
+
+#include <cmath>
 #ifdef USE_PAW
 #include "module_hamilt_general/module_xc/xc_functional.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
@@ -103,10 +103,9 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const
     const double* v_eff = this->pot->get_effective_v(0);
     const double* v_fixed = this->pot->get_fixed_v();
     const double* v_ofk = nullptr;
-    const bool v_ofk_flag =(get_xc_func_type() == 3 
-                            || get_xc_func_type() == 5);
+    const bool v_ofk_flag = (get_xc_func_type() == 3 || get_xc_func_type() == 5);
 #ifdef USE_PAW
-    if(PARAM.inp.use_paw)
+    if (PARAM.inp.use_paw)
     {
         ModuleBase::matrix v_xc;
         const std::tuple<double, double, ModuleBase::matrix> etxc_vtxc_v
@@ -115,19 +114,19 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const
 
         for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
         {
-            deband_aux -= this->charge->rho[0][ir] * v_xc(0,ir);
+            deband_aux -= this->charge->rho[0][ir] * v_xc(0, ir);
         }
         if (PARAM.inp.nspin == 2)
         {
             for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
             {
-                deband_aux -= this->charge->rho[1][ir] * v_xc(1,ir);
+                deband_aux -= this->charge->rho[1][ir] * v_xc(1, ir);
             }
         }
     }
 #endif
 
-    if(!PARAM.inp.use_paw)
+    if (!PARAM.inp.use_paw)
     {
         for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
         {
@@ -137,16 +136,16 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const
         {
             v_ofk = this->pot->get_effective_vofk(0);
             // cause in the get_effective_vofk, the func will return nullptr
-            if(v_ofk==nullptr && this->charge->rhopw->nrxx>0)
+            if (v_ofk == nullptr && this->charge->rhopw->nrxx > 0)
             {
-                ModuleBase::WARNING_QUIT("ElecState::cal_delta_eband","v_ofk is nullptr");
+                ModuleBase::WARNING_QUIT("ElecState::cal_delta_eband", "v_ofk is nullptr");
             }
             for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
             {
                 deband_aux -= this->charge->kin_r[0][ir] * v_ofk[ir];
             }
         }
-        
+
         if (PARAM.inp.nspin == 2)
         {
             v_eff = this->pot->get_effective_v(1);
@@ -157,9 +156,9 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const
             if (v_ofk_flag)
             {
                 v_ofk = this->pot->get_effective_vofk(1);
-                if(v_ofk==nullptr && this->charge->rhopw->nrxx>0)
+                if (v_ofk == nullptr && this->charge->rhopw->nrxx > 0)
                 {
-                    ModuleBase::WARNING_QUIT("ElecState::cal_delta_eband","v_ofk is nullptr");
+                    ModuleBase::WARNING_QUIT("ElecState::cal_delta_eband", "v_ofk is nullptr");
                 }
                 for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
                 {
@@ -219,7 +218,7 @@ double ElecState::cal_delta_escf() const
         if (get_xc_func_type() == 3 || get_xc_func_type() == 5)
         {
             // cause in the get_effective_vofk, the func will return nullptr
-            assert(v_ofk!=nullptr);
+            assert(v_ofk != nullptr);
             descf -= (this->charge->kin_r[0][ir] - this->charge->kin_r_save[0][ir]) * v_ofk[ir];
         }
     }
@@ -278,7 +277,7 @@ void ElecState::cal_converged()
 /**
  * @brief calculate energies
  *
- * @param type: 1 means Harris-Foulkes functinoal; 
+ * @param type: 1 means Harris-Foulkes functinoal;
  * @param type: 2 means Kohn-Sham functional;
  */
 void ElecState::cal_energies(const int type)
@@ -292,7 +291,7 @@ void ElecState::cal_energies(const int type)
     //! energy from gate-field
     this->f_en.gatefield = get_etot_gatefield();
 
-    //! energy from implicit solvation model 
+    //! energy from implicit solvation model
     if (PARAM.inp.imp_sol)
     {
         this->f_en.esol_el = get_solvent_model_Ael();
@@ -305,27 +304,19 @@ void ElecState::cal_energies(const int type)
         this->f_en.escon = get_spin_constrain_energy();
     }
 
-     // energy from DFT+U
+    // energy from DFT+U
     if (PARAM.inp.dft_plus_u)
     {
         this->f_en.edftu = get_dftu_energy();
     }
 
     this->f_en.e_local_pp = get_local_pp_energy();
 
-#ifdef __DEEPKS
-    // energy from deepks
-    if (PARAM.inp.deepks_scf)
-    {
-        this->f_en.edeepks_scf = get_deepks_E_delta() - get_deepks_E_delta_band();
-    }
-#endif
-
     if (type == 1) // Harris-Foulkes functional
     {
         this->f_en.calculate_harris();
     }
-    else if (type == 2)// Kohn-Sham functional
+    else if (type == 2) // Kohn-Sham functional
     {
         this->f_en.calculate_etot();
     }