separate H from dF/dC or dK/dC

Author: maki49

source/module_lr/Grad/esolver_lr_grad.cpp

@@ -128,6 +128,7 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
         // seems because those two are classes having constructors 
         // but `cal_edm_from_XZ_istate` here is a functions
         std::weak_ptr<PotHxcLR> pot_weak = this->pot[ispin];
+        std::weak_ptr<PotHxcLR> pot_hxc_gs_weak = this->pot_hxc_gs;
 #ifdef __EXX
         std::weak_ptr<Exx_LRI<T>> exx_lri_weak = this->exx_lri;
 #endif
@@ -141,7 +142,8 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
 #ifdef __EXX
                 exx_lri_weak, this->exx_info.info_global.hybrid_alpha,
 #endif
-                this->gint_, pot_weak, this->kv, this->gd, this->paraX_, this->paraC_, this->paraMat_,
+                this->gint_, pot_weak, pot_hxc_gs_weak,
+                this->kv, this->gd, this->paraX_, this->paraC_, this->paraMat_,
                 has_local_xc(this->xc_kernel));
         elecstate::DensityMatrix<T, double> edm_real = LR_Util::build_dm_from_dmk<T, double>(edm_k,
             this->paraMat_, this->nk, this->kv.kvec_d, this->ucell, this->gd, this->orb_cutoff_);

source/module_lr/Grad/multipliers/cal_edm_from_multipliers.h

@@ -102,11 +102,12 @@ namespace LR
         // 2. edm of Z : $\sum_i \sum_a c_a epsilon_i Z_{ai} c_i
         std::vector<T> epsi_Z(px.get_local_size());
         multiply_eig_onto_vec(Z, eig_ks, px, epsi_Z.data());
-        const std::vector<ct::Tensor> cZc = cal_dm_trans_pblas(Z, px, c, pc, naos, nocc, nvirt, pmat);
+        std::vector<ct::Tensor> cZc = cal_dm_trans_pblas(Z, px, c, pc, naos, nocc, nvirt, pmat);
+        std::for_each(cZc.begin(), cZc.end(), [&](ct::Tensor& s) { LR_Util::matsym(s.data<T>(), naos, pmat); });
 
         //3. c * K_cvcx * c
-        const std::vector<ct::Tensor> cKc = cal_dm_trans_pblas(K_cvcx, px, c, pc, naos, nocc, nvirt, pmat, (T)2.0);
-        LR_Util::matsym(cKc[0].data<T>(), naos, pmat);
+        std::vector<ct::Tensor> cKc = cal_dm_trans_pblas(K_cvcx, px, c, pc, naos, nocc, nvirt, pmat, (T)2.0);
+        std::for_each(cKc.begin(), cKc.end(), [&](ct::Tensor& s) { LR_Util::matsym(s.data<T>(), naos, pmat); });
 
         // 4. $\sum_i (\Omega + \epsilon_i) \sum_{ab} C_{\mu a} X_{ia} C_{\nu b} X_{ib}$
         const std::vector<ct::Tensor> edm = cal_edm_term4(X, eig_ext_istate, eig_ks, c, px, pc, pmat);
@@ -145,6 +146,7 @@ namespace LR
 #endif 
         typename TGint<T>::type* gint,
         std::weak_ptr<PotHxcLR> pot,
+        std::weak_ptr<PotHxcLR> pot_hxc_gs,
         const K_Vectors& kv,
         const Grid_Driver& gd,
         const std::vector<Parallel_2D>& px,
@@ -162,7 +164,7 @@ namespace LR
 #ifdef __EXX
             exx_lri, exx_alpha,
 #endif
-            gint, pot, kv, px, pc, p_occ_occ, pmat, has_local_xc);
+            gint, pot_hxc_gs, kv, px, pc, p_occ_occ, pmat, has_local_xc);
         std::cout << "W: " << std::endl;
         LR_Util::print_value(W.data(), nk, p_occ_occ[0].get_col_size(), p_occ_occ[0].get_row_size());
 

source/module_lr/Grad/multipliers/cal_multiplier_w_from_z.h

@@ -36,6 +36,7 @@ namespace LR
             {
                 const int ga = px.local2global_row(la);
                 const double weight = eig_ext_istate - eig_ks[eigks_start_k + nocc + ga];
+                std::cout << "la=" << la << ", ks-eig=" << eig_ks[eigks_start_k + nocc + ga] << ", weight=" << weight << ", X2=" << X[x_start_k + la] << std::endl;
                 for (int li = 0;li < px.get_col_size();++li)
                 {
                     const int idx = li * px.get_row_size() + la;

source/module_lr/esolver_lrtd_lcao.cpp

@@ -669,11 +669,11 @@ void LR::ESolver_LR<T, TR>::set_X_initial_guess()
 template<typename T, typename TR>
 void LR::ESolver_LR<T, TR>::init_pot(const Charge& chg_gs)
 {
+    using ST = PotHxcLR::SpinType;
     this->pot.resize(nspin, nullptr);
     if (this->input.ri_hartree_benchmark != "none") { return; } //no need to initialize potential for Hxc kernel in the RI-benchmark routine
     switch (nspin)
     {
-        using ST = PotHxcLR::SpinType;
     case 1:
         this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, *this->pw_rho, ucell, chg_gs, Pgrid, ST::S1, input.lr_init_xc_kernel);
         break;
@@ -685,7 +685,12 @@ void LR::ESolver_LR<T, TR>::init_pot(const Charge& chg_gs)
         throw std::invalid_argument("ESolver_LR: nspin must be 1 or 2");
     }
     // ground-state potentials are needed for calculating the excited state force
-    if (PARAM.inp.cal_force) { this->init_pot_groundstate(chg_gs); }
+    if (PARAM.inp.cal_force)
+    {
+        this->init_pot_groundstate(chg_gs);
+        const std::string xc_kernel_gs = LR_Util::tolower(input.dft_functional);
+        this->pot_hxc_gs = std::make_shared<PotHxcLR>(xc_kernel_gs, *this->pw_rho, ucell, chg_gs, Pgrid, ST::S1, input.lr_init_xc_kernel);
+    }
 }
 
 template<typename T, typename TR>

source/module_lr/esolver_lrtd_lcao.h

@@ -72,6 +72,7 @@ namespace LR
         double etxc_gs = 0.;
         double vtxc_gs = 0.;
         ModuleBase::matrix wg_ks;   /// occupation number of ground state
+        std::shared_ptr<PotHxcLR> pot_hxc_gs; /// used in lr-grad, in the ground-state Hxc gradient term coming from dF/dC
 
         /// @brief Excited state wavefunction (locc, lvirt are local size of nocc and nvirt in each process)
         /// size of X: [neq][{nstate, nloc_per_band}], namely: