Merge branch 'HSolver' of https://github.com/deepmodeling/abacus-develop into SDFT

Author: Qianruipku

CMakeLists.txt

@@ -295,6 +295,10 @@ target_link_libraries(${ABACUS_BIN_NAME}
     ri
     driver
     xc
+    hsolver
+    elecstate
+    hamilt
+    psi
     esolver
     -lm
 )

source/Makefile.Objects

@@ -61,6 +61,8 @@ sto_iter.o\
 sto_wf.o\
 sto_hchi.o\
 sto_che.o\
+sto_forces.o\
+sto_stress_pw.o
 
 OBJS_TOOLS=complexarray.o\
 complexmatrix.o \

source/input.cpp

@@ -383,7 +383,7 @@ void Input::Default(void)
 
     cell_factor = 1.2; // LiuXh add 20180619
 
-    GlobalV::out_mul = 0; // qi feng add 2019/9/10
+    out_mul = 0; // qi feng add 2019/9/10
 
     //----------------------------------------------------------			//Peize Lin add 2020-04-04
     // restart
@@ -1342,7 +1342,7 @@ bool Input::Read(const std::string &fn)
         }
         else if (strcmp("out_mul", word) == 0)
         {
-            read_value(ifs, GlobalV::out_mul);
+            read_value(ifs, out_mul);
         } // qifeng add 2019/9/10
         //----------------------------------------------------------
         // exx
@@ -1841,6 +1841,7 @@ void Input::Default_2(void) // jiyy add 2019-08-04
         }
     }
     if(calculation.substr(0,3) != "sto")    bndpar = 1;
+    if(bndpar > GlobalV::NPROC) bndpar = GlobalV::NPROC;
 }
 #ifdef __MPI
 void Input::Bcast()
@@ -2092,7 +2093,7 @@ void Input::Bcast()
     Parallel_Common::bcast_bool(test_just_neighbor);
     Parallel_Common::bcast_int(GlobalV::ocp);
     Parallel_Common::bcast_string(GlobalV::ocp_set);
-    Parallel_Common::bcast_int(GlobalV::out_mul); // qifeng add 2019/9/10
+    Parallel_Common::bcast_int(out_mul); // qifeng add 2019/9/10
 
     // Peize Lin add 2018-06-20
     Parallel_Common::bcast_string(dft_functional);
@@ -2314,9 +2315,9 @@ void Input::Check(void)
             ModuleBase::WARNING_QUIT("Input::Check", "calculate = istate is only availble for LCAO.");
         }
     }
-    else if (calculation == "md") // mohan add 2011-11-04
+    else if (calculation == "md" || calculation == "sto-md") // mohan add 2011-11-04
     {
-        GlobalV::CALCULATION = "md";
+        GlobalV::CALCULATION = calculation;
         symmetry = false;
         cal_force = 1;
         if (mdp.md_nstep == 0)

source/module_deepks/test/CMakeLists.txt

@@ -7,7 +7,7 @@ get_target_property(ABACUS_LINK_LIBRARIES ${ABACUS_BIN_NAME} LINK_LIBRARIES)
 target_link_libraries(
     test_deepks
     base cell symmetry md surchem xc
-    neighbor orb io ions lcao parallel mrrr pdiag pw ri driver esolver
+    neighbor orb io ions lcao parallel mrrr pdiag pw ri driver esolver hsolver psi elecstate hamilt
     pthread
     deepks
     ${ABACUS_LINK_LIBRARIES}

source/module_elecstate/CMakeLists.txt

@@ -6,3 +6,7 @@ add_library(
     elecstate_lcao.cpp
     dm2d_to_grid.cpp
 )
+
+IF (BUILD_TESTING)
+  add_subdirectory(test)
+endif()

source/module_elecstate/elecstate.cpp

@@ -16,7 +16,7 @@ const double* ElecState::getRho(int spin) const
     return &(this->charge->rho[spin][0]);
 }
 
-void ElecState::calculate_weights(void)
+void ElecState::calculate_weights()
 {
     ModuleBase::TITLE("ElecState", "calculate_weights");
 
@@ -30,6 +30,7 @@ void ElecState::calculate_weights(void)
         ekb_tmp[i] = &(this->ekb(i, 0));
     }
     int nbands = this->ekb.nc;
+    int nks = this->ekb.nr;
 
     if (GlobalV::KS_SOLVER == "selinv")
     {
@@ -41,46 +42,46 @@ void ElecState::calculate_weights(void)
     {
         if (GlobalV::TWO_EFERMI)
         {
-            Occupy::iweights(GlobalC::kv.nks,
-                             GlobalC::kv.wk,
+            Occupy::iweights(nks,
+                             this->klist->wk,
                              nbands,
                              GlobalC::ucell.magnet.get_nelup(),
                              ekb_tmp,
                              GlobalC::en.ef_up,
                              this->wg,
                              0,
-                             GlobalC::kv.isk);
-            Occupy::iweights(GlobalC::kv.nks,
-                             GlobalC::kv.wk,
+                             this->klist->isk);
+            Occupy::iweights(nks,
+                             this->klist->wk,
                              nbands,
                              GlobalC::ucell.magnet.get_neldw(),
                              ekb_tmp,
                              GlobalC::en.ef_dw,
                              this->wg,
                              1,
-                             GlobalC::kv.isk);
+                             this->klist->isk);
             // ef = ( ef_up + ef_dw ) / 2.0_dp need??? mohan add 2012-04-16
         }
         else
         {
             // -1 means don't need to consider spin.
-            Occupy::iweights(GlobalC::kv.nks,
-                             GlobalC::kv.wk,
+            Occupy::iweights(nks,
+                             this->klist->wk,
                              nbands,
-                             GlobalC::CHR.nelec,
+                             this->charge->nelec,
                              ekb_tmp,
                              this->ef,
                              this->wg,
                              -1,
-                             GlobalC::kv.isk);
+                             this->klist->isk);
         }
     }
     else if (Occupy::use_tetrahedron_method)
     {
         ModuleBase::WARNING_QUIT("calculate_weights", "not implemented yet,coming soon!");
         //		if(my_rank == 0)
         //		{
-        //			tweights(GlobalC::kv.nkstot, nspin, nbands, GlobalC::CHR.nelec, ntetra,tetra, GlobalC::wf.et,
+        //			tweights(GlobalC::kv.nkstot, nspin, nbands, this->charge->nelec, ntetra,tetra, GlobalC::wf.et,
         // this->ef, this->wg);
         //		}
     }
@@ -90,8 +91,8 @@ void ElecState::calculate_weights(void)
         {
             double demet_up = 0.0;
             double demet_dw = 0.0;
-            Occupy::gweights(GlobalC::kv.nks,
-                             GlobalC::kv.wk,
+            Occupy::gweights(nks,
+                             this->klist->wk,
                              nbands,
                              GlobalC::ucell.magnet.get_nelup(),
                              Occupy::gaussian_parameter,
@@ -101,9 +102,9 @@ void ElecState::calculate_weights(void)
                              demet_up,
                              this->wg,
                              0,
-                             GlobalC::kv.isk);
-            Occupy::gweights(GlobalC::kv.nks,
-                             GlobalC::kv.wk,
+                             this->klist->isk);
+            Occupy::gweights(nks,
+                             this->klist->wk,
                              nbands,
                              GlobalC::ucell.magnet.get_neldw(),
                              Occupy::gaussian_parameter,
@@ -113,35 +114,35 @@ void ElecState::calculate_weights(void)
                              demet_dw,
                              this->wg,
                              1,
-                             GlobalC::kv.isk);
-            GlobalC::en.demet = demet_up + demet_dw;
+                             this->klist->isk);
+            this->demet = demet_up + demet_dw;
         }
         else
         {
             // -1 means is no related to spin.
-            Occupy::gweights(GlobalC::kv.nks,
-                             GlobalC::kv.wk,
+            Occupy::gweights(nks,
+                             this->klist->wk,
                              nbands,
-                             GlobalC::CHR.nelec,
+                             this->charge->nelec,
                              Occupy::gaussian_parameter,
                              Occupy::gaussian_type,
                              ekb_tmp,
                              this->ef,
-                             GlobalC::en.demet,
+                             this->demet,
                              this->wg,
                              -1,
-                             GlobalC::kv.isk);
+                             this->klist->isk);
         }
 
         // qianrui fix a bug on 2021-7-21
-        Parallel_Reduce::reduce_double_allpool(GlobalC::en.demet);
+        Parallel_Reduce::reduce_double_allpool(this->demet);
     }
     else if (Occupy::fixed_occupations)
     {
         // fix occupations need nelup and neldw.
         // mohan add 2011-04-03
         this->ef = -1.0e+20;
-        for (int ik = 0; ik < GlobalC::kv.nks; ik++)
+        for (int ik = 0; ik < nks; ik++)
         {
             for (int ibnd = 0; ibnd < nbands; ibnd++)
             {
@@ -162,7 +163,7 @@ void ElecState::calculate_weights(void)
     {
         double ebotom = ekb_tmp[0][0];
         double etop = ekb_tmp[0][0];
-        for (int ik = 0; ik < GlobalC::kv.nks; ik++)
+        for (int ik = 0; ik < nks; ik++)
         {
             for (int ib = 0; ib < nbands; ib++)
             {
@@ -188,15 +189,27 @@ void ElecState::calculate_weights(void)
     return;
 }
 
-double ElecState::eBandK(const int& ik)
+void ElecState::calEBand()
 {
-    ModuleBase::TITLE("ElecStatePW", "eBandK");
-    double eband_k = 0.0;
-    for (int ibnd = 0; ibnd < this->ekb.nc; ibnd++)
+    ModuleBase::TITLE("ElecStatePW", "calEBand");
+    //calculate ebands using wg and ekb
+    this->eband = 0.0;
+    for (int ik = 0; ik < this->ekb.nr; ++ik)
     {
-        eband_k += this->ekb(ik, ibnd) * this->wg(ik, ibnd);
+        for (int ibnd = 0; ibnd < this->ekb.nc; ibnd++)
+        {
+            this->eband += this->ekb(ik, ibnd) * this->wg(ik, ibnd);
+        }
     }
-    return eband_k;
+    if(GlobalV::KPAR != 1)
+    {
+        //==================================
+        // Reduce all the Energy in each cpu
+        //==================================
+        this->eband /= GlobalV::NPROC_IN_POOL;
+        Parallel_Reduce::reduce_double_all(this->eband);
+    }
+    return;
 }
 
 void ElecState::print_band(const int& ik, const int& printe, const int& iter)

source/module_elecstate/elecstate.h

@@ -3,18 +3,22 @@
 
 #include "module_psi/psi.h"
 #include "src_pw/charge.h"
+#include "src_pw/klist.h"
 
 namespace elecstate
 {
 
 class ElecState
 {
   public:
+    ElecState(){};
     virtual void init(Charge *chg_in, // pointer for class Charge
+                      const K_Vectors *klist_in,
                       int nk_in, // number of k points
                       int nb_in) // number of bands
     {
         this->charge = chg_in;
+        this->klist = klist_in;
         this->ekb.create(nk_in, nb_in);
         this->wg.create(nk_in, nb_in);
     }
@@ -47,26 +51,33 @@ class ElecState
     }
 
     // calculate wg from ekb
-    void calculate_weights(void);
+    virtual void calculate_weights(void);
 
-    Charge *charge;
+    // pointer to charge density
+    Charge *charge = nullptr;
+    // pointer to k points lists
+    const K_Vectors* klist = nullptr;
     // energy for sum of electrons
     double eband = 0.0;
     // Fermi energy
     double ef = 0.0;
+    // correction energy for metals
+    double demet = 0.0;
 
     // band energy at each k point, each band.
     ModuleBase::matrix ekb;
     // occupation weight for each k-point and band
     ModuleBase::matrix wg;
 
+    std::string classname = "none";
+
   protected:
-    // calculate ebands for each k point
-    double eBandK(const int &ik);
+    // calculate ebands for all k points and all occupied bands
+    void calEBand();
 
     // print and check for band energy and occupations
     void print_band(const int &ik, const int &printe, const int &iter);
 };
 
 } // namespace elecstate
-#endif
+#endif

source/module_elecstate/elecstate_lcao.cpp

@@ -59,6 +59,9 @@ void ElecStateLCAO::psiToRho(const psi::Psi<std::complex<double>>& psi)
     ModuleBase::TITLE("ElecStateLCAO", "psiToRho");
     ModuleBase::timer::tick("ElecStateLCAO", "psiToRho");
 
+    this->calculate_weights();
+    this->calEBand();
+
     ModuleBase::GlobalFunc::NOTE("Calculate the density matrix.");
 
     // this part for calculating dm_k in 2d-block format, not used for charge now
@@ -109,6 +112,7 @@ void ElecStateLCAO::psiToRho(const psi::Psi<double>& psi)
     ModuleBase::timer::tick("ElecStateLCAO", "psiToRho");
 
     this->calculate_weights();
+    this->calEBand();
 
     if (GlobalV::KS_SOLVER == "genelpa" || GlobalV::KS_SOLVER == "scalapack_gvx" || GlobalV::KS_SOLVER == "lapack")
     {

source/module_elecstate/elecstate_lcao.h

@@ -13,16 +13,18 @@ class ElecStateLCAO : public ElecState
 {
   public:
     ElecStateLCAO(Charge* chg_in,
+                  const K_Vectors* klist_in, 
                   int nks_in,
                   int nbands_in,
                   Local_Orbital_Charge* loc_in,
                   LCAO_Hamilt* uhm_in,
                   Local_Orbital_wfc* lowf_in)
     {
-        init(chg_in, nks_in, nbands_in);
+        init(chg_in, klist_in, nks_in, nbands_in);
         this->loc = loc_in;
         this->uhm = uhm_in;
         this->lowf = lowf_in;
+        this->classname = "ElecStateLCAO";
     }
     // void init(Charge* chg_in):charge(chg_in){} override;