modify some factors and clarify in the comments

Author: maki49

source/module_lr/Grad/force/lr_force.cpp

@@ -89,16 +89,19 @@ namespace LR
         this->gint_->reset_DMRGint(dm_gs.get_DMR_vector().size());
         elecstate::Potential pot_hxc = this->dm_to_hxc_potential(dm_gs);
         this->gint_->reset_DMRGint(relax_diff_dm.get_DMR_vector().size());
+        // `cal_pulay_fs` calculates 1*Pulay-term. 
+        // For ground-state DFT, Pulay term = Hellmann-Feynman term, F = 1/2(Pulay + H-F) = Pulay, so directly call it once gives correct result.
         PulayForceStress::cal_pulay_fs(relax_diff_dm.get_DMR_vector().size()/*nspin*/, fhxc_dphi, stress_tmp,
             relax_diff_dm, this->ucell_, &pot_hxc, *this->gint_, true, false);
-        fhxc_dphi *= 0.5; // avoid double count
+        // fhxc_dphi *= 0.5; // avoid double count
 
         // 3.3 Hartree + xc (Hellmann-Feynman)
         ModuleBase::matrix fhxc_dvhxc(this->ucell_.nat, 3);
         elecstate::Potential pot_hxc_relaxed_diff = this->dm_to_hxc_potential(relax_diff_dm);
+        //`cal_pulay_fs` calculates only one spin channel because `relax_diff_dm` has only one.
         PulayForceStress::cal_pulay_fs(1/*nspin*/, fhxc_dvhxc, stress_tmp,
             dm_gs, this->ucell_, &pot_hxc_relaxed_diff, *this->gint_, true, false);
-        // fhxc_dvhxc *= 0.5; // avoid double count, but nspin=2 of ground-state dm cancels it here 
+        fhxc_dvhxc *= 2; // for the two channels of the ground-state dm. 
 
         // 4. kinetic (Pulay)
         std::vector<hamilt::HContainer<double>> dT = cal_hs_grad('T', this->ucell_, this->pv_, this->gd_, this->two_center_bundle_);
@@ -120,6 +123,9 @@ namespace LR
     template<typename TK>
     ModuleBase::matrix LR_Force<TK>::cal_force_hxc_dmtrans(const elecstate::DensityMatrix<TK, double>& dm_trans, const PotHxcLR& pot_hxc)
     {
+        // `dm_trans` (D^X) carries the singlet spin normalization (sqrt(2) per channel), 
+        // so D^X in pot_hxc and cal_pulay_fs together already contribute a factor 2.
+        // So cal_pulay_fs here returns 2*Pulay = Pulay + Hellmann-Feynman force. *2 is not needed here.
         return PulayForceStress::cal_pulay_fs(dm_trans, this->ucell_, &pot_hxc, *this->gint_);
     }
 
@@ -134,12 +140,13 @@ namespace LR
         auto& exx_lri_kernel = this->exx_lri_.lock()->get();
         exx_lri_kernel.set_Ds(dm_trans, this->exx_lri_.lock()->get_info().dm_threshold, spin_suffix);
         exx_lri_kernel.cal_Hs({ "", "", spin_suffix });
-        exx_lri_kernel.cal_force({ "", "", spin_suffix, "", "" });// using dm_trans
+        exx_lri_kernel.cal_force({ "", "", spin_suffix, "", "" });// using dm_trans，Pulay term only
         for (std::size_t idim = 0; idim < 3; ++idim)
             for (const auto& force_item : exx_lri_kernel.force[idim])
                 f_exx_dmtrans(force_item.first, idim) = std::real(force_item.second);
         const double fac = -2.0 * alpha; //-2 is the same as post_process_Hexx, Hartree to Ry (which didn't act on Hs)
-        return f_exx_dmtrans * fac; // dm_trans (DX) already contain the spin channel (sqrt(2) times of up/down channel DX)
+        const double pulay_to_total_sym = 2.0;  // Pulay -> Pulay + Hellmann-Feynman, only when Ds_left and Ds_right are equal
+        return f_exx_dmtrans * fac * pulay_to_total_sym; // dm_trans (DX) already contain the spin channel (sqrt(2) times of up/down channel DX)
         // return f_exx_dmtrans * fac * 2; // 2 is the same in post_process_Eexx at nspin=1 ( up->up + down->down, 2 spin-conserving transitions)
     }
 
@@ -199,8 +206,9 @@ namespace LR
 
         // -2 * 0.5 * alpha
         // -2 is the same as post_process_Hexx (a.u. to Ry, which didn't act on Hs)
-        // 0.5 is the 2-electron integral prefactor
-        const double fac = -alpha;
+        // 0.5 is the 2-electron integral prefactor，used in ground-state energy/force where two density matrix are identical
+        // But the LR-grad Lagrangian/force 2-e term here Tr[(T+Z)H[D]] does not have 1/2 factor.
+        const double fac = -2 * alpha;
         return f_exx_gs_diff * fac;
     }
 #endif

source/module_lr/Grad/force/lr_force_test.cpp

@@ -170,29 +170,29 @@ namespace LR
                             elecstate::DensityMatrix<TK, double> dm_kl_sym = init_dm_eff(k, l, true);
                             ModuleBase::matrix fhartree_pulay(this->ucell_.nat, 3), fhartree_h_f(this->ucell_.nat, 3);
                             ModuleBase::matrix stress_tmp;  // dummy
-                            PulayForceStress::cal_pulay_fs(1/*nspin*/, fhartree_pulay, stress_tmp, dm_ij_sym, this->ucell_, &pot_hxc_kl, *this->gint_, true, false);  // 2 * Pulay term
-                            PulayForceStress::cal_pulay_fs(1/*nspin*/, fhartree_h_f, stress_tmp, dm_kl_sym, this->ucell_, &pot_hxc_ij, *this->gint_, true, false);  // 2 * Hellmann-Feynman term
+                            PulayForceStress::cal_pulay_fs(1/*nspin*/, fhartree_pulay, stress_tmp, dm_ij_sym, this->ucell_, &pot_hxc_kl, *this->gint_, true, false);  // Pulay term
+                            PulayForceStress::cal_pulay_fs(1/*nspin*/, fhartree_h_f, stress_tmp, dm_kl_sym, this->ucell_, &pot_hxc_ij, *this->gint_, true, false);  // Hellmann-Feynman term
                             ModuleIO::print_force(GlobalV::ofs_running, this->ucell_,
                                 "H2_SZ_CENTER4_HXC_dtau(" + std::to_string(i) + std::to_string(j) + "|" + std::to_string(k) + std::to_string(l) + ") FORCE (Ry/au)",
-                                -(fhartree_pulay + fhartree_h_f), true);   // F_Hxc_ijkl = -1/2*dtau(ij|kl), dtau(ij|kl) = -2F = -(2F_pulay + 2F_H-F)
+                                -(fhartree_pulay + fhartree_h_f), true);   // F_Hxc_ijkl = -dtau(ij|kl)
                             ModuleIO::print_force(GlobalV::ofs_running, this->ucell_,
                                 "H2_SZ_CENTER4_HXC_Pulay_dtau(" + std::to_string(i) + std::to_string(j) + "|" + std::to_string(k) + std::to_string(l) + ") FORCE (Ry/au)",
-                                -fhartree_pulay, true);   // F_Hxc_ijkl = -1/2*dtau(ij|kl)
+                                -fhartree_pulay, true);   // F_Hxc_ijkl = -dtau(ij|kl)
                             ModuleIO::print_force(GlobalV::ofs_running, this->ucell_,
                                 "H2_SZ_CENTER4_HXC_H-F_dtau(" + std::to_string(i) + std::to_string(j) + "|" + std::to_string(k) + std::to_string(l) + ") FORCE (Ry/au)",
-                                -fhartree_h_f, true);   // F_Hxc_ijkl = -1/2*dtau(ij|kl)
+                                -fhartree_h_f, true);   // F_Hxc_ijkl = -dtau(ij|kl)
 #ifdef __EXX
                             if (!this->exx_lri_.expired())
                             {   // match the Gint result with LibRI
                                 auto ds_kl = LR_Util::get_exx_Ds_spin1(dm_kl, ucell_, kv, pv_); // returns ds_kl*0.5
                                 auto ds_ij = LR_Util::get_exx_Ds_spin1(dm_ij, ucell_, kv, pv_); // returns ds_ij*0.5
                                 // calulates F=0.5*d(ik|jl) (only one spin channel). 0.5 is the 2-electron integral prefactor.
                                 // 4 cancels the two 0.5s in Ds, induced by `split_m2D_ktoR`. 
-                                // No spin factor hard-coded in this function.
+                                // No spin factor or two-electron-energy factor (1/2) are hard-coded in this function.
                                 ModuleBase::matrix f_exx = this->cal_force_exx_gs_dm_relaxed_diff(ds_kl, ds_ij, alpha_ * 4.0, "");
                                 ModuleIO::print_force(GlobalV::ofs_running, ucell_,
                                     "H2_SZ_CENTER4_EXX_dtau(" + std::to_string(i) + std::to_string(j) + "|" + std::to_string(k) + std::to_string(l) + ") FORCE (Ry/au)",
-                                    f_exx * 2, true);   //d(ik|jl)=2F
+                                    f_exx , true);   //d(ik|jl)=F
                             }
 #endif
                         }