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

Author: YuLiu98

source/Makefile.Objects

@@ -430,6 +430,8 @@ OBJS_RELAXATION=bfgs_basic.o\
     relax_old.o\
     relax.o\
     line_search.o\
+    bfgs.o\
+    
 
 OBJS_SURCHEM=surchem.o\
     H_correction_pw.o\

source/module_elecstate/elecstate_pw.h

@@ -17,6 +17,7 @@ class ElecStatePW : public ElecState
 {
   private:
     using Real = typename GetTypeReal<T>::type;
+
   public:
     ElecStatePW(ModulePW::PW_Basis_K* wfc_basis_in,
                 Charge* chg_in,
@@ -26,46 +27,59 @@ class ElecStatePW : public ElecState
                 ModulePW::PW_Basis* rhodpw_in,
                 ModulePW::PW_Basis* rhopw_in,
                 ModulePW::PW_Basis_Big* bigpw_in);
-    // void init(Charge* chg_in):charge(chg_in){} override;
 
     ~ElecStatePW();
-    // interface for HSolver to calculate rho from Psi
+
+    //! interface for HSolver to calculate rho from Psi
     virtual void psiToRho(const psi::Psi<T, Device>& psi);
-    // return current electronic density rho, as a input for constructing Hamiltonian
-    // const double* getRho(int spin) const override;
-    virtual void cal_tau(const psi::Psi<T, Device>& psi);
 
-    // update charge density for next scf step
-    // void getNewRho() override;
+    virtual void cal_tau(const psi::Psi<T, Device>& psi);
 
     Real* becsum = nullptr;
-    // init rho_data and kin_r_data
+
+    //! init rho_data and kin_r_data
     void init_rho_data();
-    Real ** rho = nullptr, ** kin_r = nullptr; //[Device] [spin][nrxx] rho and kin_r
+    Real** rho = nullptr; 
+    Real** kin_r = nullptr; //[Device] [spin][nrxx] rho and kin_r
 
   protected:
+
     ModulePW::PW_Basis* rhopw_smooth = nullptr;
+
     ModulePW::PW_Basis_K* basis = nullptr;
+
     UnitCell* ucell = nullptr;
+
     const pseudopot_cell_vnl* ppcell = nullptr;
 
-    // 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.
+    //! 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.
     void updateRhoK(const psi::Psi<T, Device>& psi); // override;
-    // sum over all pools for rho and ebands
+    
+    //! sum over all pools for rho and ebands
     void parallelK();
-    // calcualte rho for each k
+    
+    //! calcualte rho for each k
     void rhoBandK(const psi::Psi<T, Device>& psi);
-    // add to the charge density in reciprocal space the part which is due to the US augmentation.
+    
+    //! add to the charge density in reciprocal space the part which is due to the US augmentation.
     void add_usrho(const psi::Psi<T, Device>& psi);
-    // \sum_lm Q_lm(r) \sum_i <psi_i|beta_l><beta_m|psi_i> w_i
+
+    //! Non-local pseudopotentials
+    //! \sum_lm Q_lm(r) \sum_i <psi_i|beta_l><beta_m|psi_i> w_i
     void addusdens_g(const Real* becsum, T* rhog);
 
     Device * ctx = {};
+
     bool init_rho = false;
+
     mutable T* vkb = nullptr;
-    Real * rho_data = nullptr, * kin_r_data = nullptr;
-    T * wfcr = nullptr, * wfcr_another_spin = nullptr;
+
+    Real* rho_data = nullptr; 
+    Real* kin_r_data = nullptr;
+    T* wfcr = nullptr; 
+    T* wfcr_another_spin = nullptr;
+
   private:
     using meta_op = hamilt::meta_pw_op<Real, Device>;
     using elecstate_pw_op = elecstate::elecstate_pw_op<Real, Device>;
@@ -85,4 +99,4 @@ class ElecStatePW : public ElecState
 
 } // namespace elecstate
 
-#endif
+#endif

source/module_elecstate/test/charge_extra_test.cpp

@@ -87,7 +87,7 @@ Structure_Factor::Structure_Factor()
 Structure_Factor::~Structure_Factor()
 {
 }
-void Structure_Factor::setup_structure_factor(UnitCell* Ucell, const ModulePW::PW_Basis* rho_basis)
+void Structure_Factor::setup_structure_factor(const UnitCell* Ucell, const ModulePW::PW_Basis* rho_basis)
 {
 }
 

source/module_esolver/esolver_gets.cpp

@@ -39,7 +39,7 @@ void ESolver_GetS::before_all_runners(UnitCell& ucell, const Input_para& inp)
     }
 
     // 1.3) Setup k-points according to symmetry.
-    this->kv.set(ucell,ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec,GlobalV::ofs_running);
+    this->kv.set(ucell, ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
 
     ModuleIO::setup_parameters(ucell, this->kv);
@@ -101,14 +101,15 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep)
                          PARAM.inp.test_atom_input);
 
     Record_adj RA;
-    RA.for_2d(GlobalC::GridD, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
+    RA.for_2d(ucell, GlobalC::GridD, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
 
     if (this->p_hamilt == nullptr)
     {
         if (PARAM.inp.nspin == 4)
         {
             this->p_hamilt
-                = new hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>(GlobalC::GridD,
+                = new hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>(ucell,
+                                                                                     GlobalC::GridD,
                                                                                      &this->pv,
                                                                                      this->kv,
                                                                                      *(two_center_bundle_.overlap_orb),
@@ -118,7 +119,8 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep)
         }
         else
         {
-            this->p_hamilt = new hamilt::HamiltLCAO<std::complex<double>, double>(GlobalC::GridD,
+            this->p_hamilt = new hamilt::HamiltLCAO<std::complex<double>, double>(ucell,
+                                                                                  GlobalC::GridD,
                                                                                   &this->pv,
                                                                                   this->kv,
                                                                                   *(two_center_bundle_.overlap_orb),

source/module_esolver/esolver_ks_lcao.cpp

@@ -164,7 +164,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
     // 5) initialize Hamilt in LCAO
     // * allocate H and S matrices according to computational resources
     // * set the 'trace' between local H/S and global H/S
-    LCAO_domain::divide_HS_in_frag(PARAM.globalv.gamma_only_local, pv, this->kv.get_nks(), orb_);
+    LCAO_domain::divide_HS_in_frag(PARAM.globalv.gamma_only_local, ucell , pv, this->kv.get_nks(), orb_);
 
 #ifdef __EXX
     // 6) initialize exx
@@ -294,6 +294,7 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& for
                        PARAM.inp.cal_stress,
                        PARAM.inp.test_force,
                        PARAM.inp.test_stress,
+                       ucell,
                        this->pv,
                        this->pelec,
                        this->psi,
@@ -643,7 +644,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
             GlobalC::dftu.set_dmr(dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM());
         }
         // Calculate U and J if Yukawa potential is used
-        GlobalC::dftu.cal_slater_UJ(this->pelec->charge->rho, this->pw_rho->nrxx);
+        GlobalC::dftu.cal_slater_UJ(ucell,this->pelec->charge->rho, this->pw_rho->nrxx);
     }
 
 #ifdef __DEEPKS
@@ -874,11 +875,11 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
             {
                 const std::vector<std::vector<TK>>& tmp_dm
                     = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM()->get_DMK_vector();
-                ModuleDFTU::dftu_cal_occup_m(iter, tmp_dm, this->kv, this->p_chgmix->get_mixing_beta(), this->p_hamilt);
+                ModuleDFTU::dftu_cal_occup_m(iter, ucell,tmp_dm, this->kv, this->p_chgmix->get_mixing_beta(), this->p_hamilt);
             }
-            GlobalC::dftu.cal_energy_correction(istep);
+            GlobalC::dftu.cal_energy_correction(ucell,istep);
         }
-        GlobalC::dftu.output();
+        GlobalC::dftu.output(ucell);
     }
 
     // (7) for deepks, calculate delta_e

source/module_esolver/lcao_before_scf.cpp

@@ -1,9 +1,9 @@
+#include "module_elecstate/cal_ux.h"
 #include "module_elecstate/module_charge/symmetry_rho.h"
 #include "module_esolver/esolver_ks_lcao.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h"
 #include "module_hamilt_lcao/module_dftu/dftu.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
-#include "module_elecstate/cal_ux.h"
 //
 #include "module_base/timer.h"
 #include "module_cell/module_neighbor/sltk_atom_arrange.h"
@@ -127,7 +127,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     // (2)For each atom, calculate the adjacent atoms in different cells
     // and allocate the space for H(R) and S(R).
     // If k point is used here, allocate HlocR after atom_arrange.
-    this->RA.for_2d(GlobalC::GridD, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
+    this->RA.for_2d(ucell, GlobalC::GridD, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
 
     // 2. density matrix extrapolation
 
@@ -174,7 +174,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     }
 
     // prepare grid in Gint
-    LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, orb_, *this->pw_rho, *this->pw_big);
+    LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, ucell, orb_, *this->pw_rho, *this->pw_big);
 
     // init Hamiltonian
     if (this->p_hamilt != nullptr)
@@ -188,6 +188,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
         this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(
             PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr,
             PARAM.globalv.gamma_only_local ? nullptr : &(this->GK),
+            ucell,
             GlobalC::GridD,
             &this->pv,
             this->pelec->pot,
@@ -311,11 +312,9 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     // two cases are considered:
     // 1. DMK in DensityMatrix is not empty (istep > 0), then DMR is initialized by DMK
     // 2. DMK in DensityMatrix is empty (istep == 0), then DMR is initialized by zeros
-    if(istep > 0)
+    if (istep > 0)
     {
-        dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)
-            ->get_DM()
-            ->cal_DMR();
+        dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM()->cal_DMR();
     }
 
     if (PARAM.inp.dm_to_rho)
@@ -370,16 +369,18 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     this->p_hamilt->non_first_scf = istep;
 
     // update in ion-step
-    if( PARAM.inp.rdmft == true )
+    if (PARAM.inp.rdmft == true)
     {
         // necessary operation of these parameters have be done with p_esolver->Init() in source/driver_run.cpp
-        rdmft_solver.update_ion(ucell, *(this->pw_rho), this->ppcell.vloc, this->sf.strucFac);   // add by jghan, 2024-03-16/2024-10-08
+        rdmft_solver.update_ion(ucell,
+                                *(this->pw_rho),
+                                this->ppcell.vloc,
+                                this->sf.strucFac); // add by jghan, 2024-03-16/2024-10-08
     }
 
     return;
 }
 
-
 template class ESolver_KS_LCAO<double, double>;
 template class ESolver_KS_LCAO<std::complex<double>, double>;
 template class ESolver_KS_LCAO<std::complex<double>, std::complex<double>>;

source/module_esolver/lcao_others.cpp

@@ -1,9 +1,9 @@
+#include "module_elecstate/cal_ux.h"
 #include "module_elecstate/module_charge/symmetry_rho.h"
 #include "module_esolver/esolver_ks_lcao.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h"
 #include "module_hamilt_lcao/module_dftu/dftu.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
-#include "module_elecstate/cal_ux.h"
 //
 #include "module_base/timer.h"
 #include "module_cell/module_neighbor/sltk_atom_arrange.h"
@@ -128,7 +128,7 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
     // (2)For each atom, calculate the adjacent atoms in different cells
     // and allocate the space for H(R) and S(R).
     // If k point is used here, allocate HlocR after atom_arrange.
-    this->RA.for_2d(GlobalC::GridD, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
+    this->RA.for_2d(ucell, GlobalC::GridD, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
 
     // 2. density matrix extrapolation
 
@@ -175,7 +175,7 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
     }
 
     // prepare grid in Gint
-    LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, orb_, *this->pw_rho, *this->pw_big);
+    LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, ucell, orb_, *this->pw_rho, *this->pw_big);
 
     // init Hamiltonian
     if (this->p_hamilt != nullptr)
@@ -189,6 +189,7 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
         this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(
             PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr,
             PARAM.globalv.gamma_only_local ? nullptr : &(this->GK),
+            ucell,
             GlobalC::GridD,
             &this->pv,
             this->pelec->pot,

source/module_hamilt_general/module_surchem/cal_pseudo.cpp

@@ -8,7 +8,7 @@ void surchem::gauss_charge(const UnitCell& cell,
                            complex<double>* N,
                            Structure_Factor* sf)
 {
-    sf->setup_structure_factor(&GlobalC::ucell, rho_basis); // call strucFac(ntype,ngmc)
+    sf->setup_structure_factor(&cell, rho_basis); // call strucFac(ntype,ngmc)
     for (int it = 0; it < cell.ntype; it++)
     {
         double RCS = GetAtom.atom_RCS[cell.atoms[it].ncpp.psd];

source/module_hamilt_general/module_surchem/cal_vcav.cpp

@@ -2,7 +2,10 @@
 #include "module_hamilt_general/module_xc/xc_functional.h"
 #include "surchem.h"
 
-void lapl_rho(const std::complex<double>* rhog, double* lapn, const ModulePW::PW_Basis* rho_basis)
+void lapl_rho(const double& tpiba2,
+              const std::complex<double>* rhog, 
+              double* lapn, 
+              const ModulePW::PW_Basis* rho_basis)
 {
     std::complex<double> *gdrtmpg = new std::complex<double>[rho_basis->npw];
     ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis->nrxx);
@@ -22,9 +25,9 @@ void lapl_rho(const std::complex<double>* rhog, double* lapn, const ModulePW::PW
         // bring the gdr from G --> R
         rho_basis->recip2real(aux, aux);
         // remember to multily 2pi/a0, which belongs to G vectors.
-        for (int ir = 0; ir < rho_basis->nrxx; ir++) {
-            lapn[ir] -= aux[ir].real() * GlobalC::ucell.tpiba2;
-}
+        for (int ir = 0; ir < rho_basis->nrxx; ir++)
+            lapn[ir] -= aux[ir].real() * tpiba2;
+
     }
 
     delete[] gdrtmpg;
@@ -70,7 +73,7 @@ void surchem::createcavity(const UnitCell& ucell,
     ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis->nrxx);
 
     // nabla n
-    XC_Functional::grad_rho(PS_TOTN, nablan, rho_basis, GlobalC::ucell.tpiba);
+    XC_Functional::grad_rho(PS_TOTN, nablan, rho_basis, ucell.tpiba);
 
     //  |\nabla n |^2 = nablan_2
     for (int ir = 0; ir < rho_basis->nrxx; ir++)
@@ -79,7 +82,7 @@ void surchem::createcavity(const UnitCell& ucell,
     }
 
     // Laplacian of n
-    lapl_rho(PS_TOTN, lapn, rho_basis);
+    lapl_rho(ucell.tpiba2,PS_TOTN, lapn, rho_basis);
 
     //-------------------------------------------------------------
     // add -Lap(n)/|\nabla n| to vwork and copy \sqrt(|\nabla n|^2)
@@ -131,7 +134,7 @@ void surchem::createcavity(const UnitCell& ucell,
 
     // \nabla(1 / |\nabla n|), ggn in real space
     ModuleBase::Vector3<double> *ggn = new ModuleBase::Vector3<double>[rho_basis->nrxx];
-    XC_Functional::grad_rho(inv_gn, ggn, rho_basis, GlobalC::ucell.tpiba);
+    XC_Functional::grad_rho(inv_gn, ggn, rho_basis, ucell.tpiba);
 
     //-------------------------------------------------------------
     // add -(\nabla n . \nabla(1/ |\nabla n|)) to Vcav in real space

source/module_hamilt_general/module_surchem/cal_vel.cpp

@@ -17,6 +17,7 @@ void shape_gradn(const double* PS_TOTN_real, const ModulePW::PW_Basis* rho_basis
 }
 
 void eps_pot(const double* PS_TOTN_real,
+             const double& tpiba,
              const complex<double>* phi,
              const ModulePW::PW_Basis* rho_basis,
              double* d_eps,
@@ -36,7 +37,7 @@ void eps_pot(const double* PS_TOTN_real,
     double *phisq = new double[rho_basis->nrxx];
 
     // nabla phi
-    XC_Functional::grad_rho(phi, nabla_phi, rho_basis, GlobalC::ucell.tpiba);
+    XC_Functional::grad_rho(phi, nabla_phi, rho_basis, tpiba);
 
     for (int ir = 0; ir < rho_basis->nrxx; ir++)
     {
@@ -118,7 +119,7 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell& cell,
     this->Ael *= cell.omega / rho_basis->nxyz;
 
     // the 2nd item of tmp_Vel
-    eps_pot(PS_TOTN_real, Sol_phi, rho_basis, epsilon, epspot);
+    eps_pot(PS_TOTN_real, cell.tpiba, Sol_phi, rho_basis, epsilon, epspot);
 
     for (int i = 0; i < rho_basis->nrxx; i++)
     {