KEY fix: Hxc Hellmann-Feynman term with (T+DZ)

Author: maki49

source/module_lr/Grad/esolver_lr_grad.cpp

@@ -75,7 +75,7 @@ inline auto get_exx_Ds_spin1(const elecstate::DensityMatrix<T, T>& dm,
     return RI_2D_Comm::split_m2D_ktoR<T>(ucell, kv, DMk_trans_pointer, pmat, /*nspin=*/1)[0];
 }
 template <typename T>
-inline auto get_exx_Ds_gs(elecstate::DensityMatrix<T, double>& dm,
+inline auto get_exx_Ds_gs(const elecstate::DensityMatrix<T, double>& dm,
     const UnitCell& ucell, const K_Vectors& kv, const Parallel_Orbitals& pmat, const int nspin)
     -> std::map<int, std::map<std::pair<int, std::array<int, 3>>, RI::Tensor<T>>>
 {
@@ -245,7 +245,8 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
         std::cout << "Force (Hxc-DMTrans term) of state " << istate << ": " << std::endl;
         LR_Util::print_value(force_hxc_dmtrans.c, ucell.nat, 3);
 
-        ModuleBase::matrix force_hamiltgs_relaxed_diff = lr_force.cal_force_hamilt_gs_dm_relaxed_diff(relaxed_diff_dm_real, *pot_gs, /*with_ewald=*/false);
+        const elecstate::DensityMatrix<T, double>& dm_gs = this->cal_dm_gs();
+        ModuleBase::matrix force_hamiltgs_relaxed_diff = lr_force.cal_force_hamilt_gs_dm_relaxed_diff(relaxed_diff_dm_real, dm_gs, /*with_ewald=*/false);
         std::cout << "Force (GS-(T+Z) term) of state " << istate << ": " << std::endl;
         LR_Util::print_value(force_hamiltgs_relaxed_diff.c, ucell.nat, 3);
 
@@ -261,8 +262,6 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
         std::cout << "Force (EXX-DMTrans term) of state " << istate << ": " << std::endl;
         LR_Util::print_value(force_exx_dmtrans.c, ucell.nat, 3);
 
-        elecstate::DensityMatrix<T, double> dm_gs(&this->paraMat_, this->nspin, this->kv.kvec_d, this->nk);   //DX
-        elecstate::cal_dm_psi(&this->paraMat_all_, this->wg_ks_all, *this->psi_ks_all, dm_gs);   // nbands is important here
         const auto& Ds_gs = get_exx_Ds_gs(dm_gs, this->ucell, this->kv, this->paraMat_, this->nspin);
         const auto& Ds_relaxed_diff = get_exx_Ds_spin1(relaxed_diff_dm, this->ucell, this->kv, this->paraMat_);
         ModuleBase::matrix force_exx_gs_diff = lr_force.cal_force_exx_gs_dm_relaxed_diff(Ds_gs, Ds_relaxed_diff, alpha);
@@ -281,6 +280,16 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
     return forces;
 }
 
+template<typename T, typename TR>
+elecstate::DensityMatrix<T, double> LR::ESolver_LR<T, TR>::cal_dm_gs()
+{
+    elecstate::DensityMatrix<T, double> dm_gs(&this->paraMat_, this->nspin, this->kv.kvec_d, this->nk);
+    elecstate::cal_dm_psi(&this->paraMat_all_, this->wg_ks_all, *this->psi_ks_all, dm_gs);   // nbands is important here
+    LR_Util::initialize_DMR(dm_gs, this->paraMat_, this->ucell, this->gd, this->orb_cutoff_);   // nbands is not important here
+    dm_gs.cal_DMR();
+    return dm_gs;
+}
+
 template<typename T, typename TR>
 void LR::ESolver_LR<T, TR>::test_force()
 {
@@ -290,11 +299,8 @@ void LR::ESolver_LR<T, TR>::test_force()
 #endif
     );
 
-    elecstate::DensityMatrix<T, double> dm_gs(&this->paraMat_, this->nspin, this->kv.kvec_d, this->nk);   //DX
-    elecstate::cal_dm_psi(&this->paraMat_all_, this->wg_ks_all, *this->psi_ks_all, dm_gs);   // nbands is important here
-    LR_Util::initialize_DMR(dm_gs, this->paraMat_, this->ucell, this->gd, this->orb_cutoff_);   // nbands is not important here
-    dm_gs.cal_DMR();
-    LR_Util::print_DMR(dm_gs, this->ucell.nat, "DM(R) of ground state");
+    const elecstate::DensityMatrix<T, double>& dm_gs = this->cal_dm_gs();
+    // LR_Util::print_DMR(dm_gs, this->ucell.nat, "DM(R) of ground state");
     ///========================== test 1: reproduce the force of ground state =========================
     // energy density matrix of the ground state
     elecstate::DensityMatrix<T, double> edm_gs(&this->paraMat_, this->nspin, this->kv.kvec_d, this->nk);   //DX
@@ -305,7 +311,7 @@ void LR::ESolver_LR<T, TR>::test_force()
     LR_Util::initialize_DMR(edm_gs, this->paraMat_, this->ucell, this->gd, this->orb_cutoff_);
     edm_gs.cal_DMR();
     // ground-state force
-    ModuleBase::matrix force_gs = lr_force.reproduce_force_gs(dm_gs, edm_gs, *this->pot_gs);
+    ModuleBase::matrix force_gs = lr_force.reproduce_force_gs(dm_gs, edm_gs);
     ModuleIO::print_force(GlobalV::ofs_running, this->ucell, "Ground State FORCE (eV/Angstrom)", force_gs, false);
     /// ======================================= END test 1 =========================================
     ///========================== test 2: reproduce the DX Hartree term =========================

source/module_lr/Grad/force/lr_force.cpp

@@ -25,9 +25,34 @@ namespace LR
         return chr;
     }
 
+    template<typename TK>
+    elecstate::Potential LR_Force<TK>::dm_to_hxc_potential(const elecstate::DensityMatrix<TK, double>& dm)
+    {
+        elecstate::Potential pot(nullptr/*remaining for uspp*/,
+            &this->rhopw_, &this->ucell_, &this->locpp_.vloc, const_cast<Structure_Factor*>(&this->sf_),
+            nullptr/*surchem*/, nullptr/*etxc*/, nullptr/*vtxc*/);
+        PARAM.inp.vh_in_h ? pot.pot_register({ "hartree", "xc" }) : pot.pot_register({ "xc" });
+        const Charge& charge = this->dm_to_charge(dm);
+        pot.init_pot(0, &charge); // call update_from_charge inside
+        return pot;
+    }
+
+    template<typename TK>
+    elecstate::Potential LR_Force<TK>::local_potential()
+    {
+        elecstate::Potential pot(nullptr/*remaining for uspp*/,
+            &this->rhopw_, &this->ucell_, &this->locpp_.vloc, const_cast<Structure_Factor*>(&this->sf_),
+            nullptr/*surchem*/, nullptr/*etxc*/, nullptr/*vtxc*/);
+        pot.pot_register({ "local" });
+        pot.init_pot(0, nullptr);
+        return pot;
+    }
+
     template<typename TK>
     ModuleBase::matrix LR_Force<TK>::cal_force_hamilt_gs_dm_relaxed_diff(const elecstate::DensityMatrix<TK, double>& relax_diff_dm,
-        const elecstate::Potential& pot_gs, const bool with_ewald)
+        const elecstate::DensityMatrix<TK, double>& dm_gs,
+        // const elecstate::Potential& pot_gs,
+        const bool with_ewald)
     {
         const Charge chr_diff_relaxed = dm_to_charge(relax_diff_dm);
 
@@ -39,11 +64,38 @@ namespace LR
         // 2. nonlocal pp (Hellmann-Feynman + Pulay)
         ModuleBase::matrix fvnl = cal_force_nonlocal(this->ucell_, this->kvec_d_, this->gd_, this->two_center_bundle_, relax_diff_dm);
 
-        // 3. local pp (Pulay) + Hartree + xc (grid integration)
+        // // 3. local pp (Pulay) + Hartree + xc (grid integration)
+        // ModuleBase::matrix fvl_dphi(this->ucell_.nat, 3);
+        // ModuleBase::matrix stress_tmp;  // no use now, only for passing into interfaces
+        // PulayForceStress::cal_pulay_fs(relax_diff_dm.get_DMR_vector().size()/*nspin*/, fvl_dphi, stress_tmp,
+        //     relax_diff_dm, this->ucell_, &pot_gs, *this->gint_, true, false);
+
+        // 3.1. local pp (Pulay)
         ModuleBase::matrix fvl_dphi(this->ucell_.nat, 3);
         ModuleBase::matrix stress_tmp;  // no use now, only for passing into interfaces
+        elecstate::Potential pot_loc = this->local_potential();
         PulayForceStress::cal_pulay_fs(relax_diff_dm.get_DMR_vector().size()/*nspin*/, fvl_dphi, stress_tmp,
-            relax_diff_dm, this->ucell_, &pot_gs, *this->gint_, true, false);
+            relax_diff_dm, this->ucell_, &pot_loc, *this->gint_, true, false);
+
+        // 3.2. Hartree + xc (Pulay) 
+        //  method 1
+        // ModuleBase::matrix fgs_dphi(this->ucell_.nat, 3);
+        // PulayForceStress::cal_pulay_fs(relax_diff_dm.get_DMR_vector().size()/*nspin*/, fgs_dphi, stress_tmp,
+        //     relax_diff_dm, this->ucell_, &pot_gs, *this->gint_, true, false);
+        // ModuleBase::matrix fhxc_dphi = (fgs_dphi - fvl_dphi) * 0.5; // avoid double count of hxc Pulay term
+        // method 2 
+        ModuleBase::matrix fhxc_dphi(this->ucell_.nat, 3);
+        elecstate::Potential pot_hxc = this->dm_to_hxc_potential(dm_gs);
+        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
+
+        // 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);
+        PulayForceStress::cal_pulay_fs(dm_gs.get_DMR_vector().size()/*nspin*/, fhxc_dvhxc, stress_tmp,
+            relax_diff_dm, this->ucell_, &pot_hxc_relaxed_diff, *this->gint_, false, true);
+        fhxc_dvhxc *= 0.5; // avoid double count
 
         // 4. kinetic (Pulay)
         std::vector<hamilt::HContainer<double>> dT = cal_hs_grad('T', this->ucell_, this->pv_, this->gd_, this->two_center_bundle_);
@@ -54,11 +106,13 @@ namespace LR
             ModuleIO::print_force(GlobalV::ofs_running, this->ucell_, "PW      FORCE (eV/Angstrom)", f_pw, false);
             ModuleIO::print_force(GlobalV::ofs_running, this->ucell_, "NONLOCAL     FORCE (eV/Angstrom)", fvnl, false);
             ModuleIO::print_force(GlobalV::ofs_running, this->ucell_, "KINETIC     FORCE (eV/Angstrom)", ft_dphi, false);
-            ModuleIO::print_force(GlobalV::ofs_running, this->ucell_, "LOCAL Pulay FORCE (pp+Hartree+XC) (eV/Angstrom)", fvl_dphi, false);
+            ModuleIO::print_force(GlobalV::ofs_running, this->ucell_, "LOCAL-PP Pulay FORCE (eV/Angstrom)", fvl_dphi, false);
+            ModuleIO::print_force(GlobalV::ofs_running, this->ucell_, "HARTREE+XC Pulay FORCE (eV/Angstrom)", fhxc_dphi, false);
+            ModuleIO::print_force(GlobalV::ofs_running, this->ucell_, "HARTREE+XC Hellmann-Feynman FORCE (eV/Angstrom)", fhxc_dvhxc, false);
         }
 
         // from the formula, we do not need the non-ortho term (overlap*edm) here.
-        return f_pw + fvnl + ft_dphi + fvl_dphi;
+        return f_pw + fvnl + ft_dphi + fvl_dphi + fhxc_dphi + fhxc_dvhxc;
     }
     template<typename TK>
     ModuleBase::matrix LR_Force<TK>::cal_force_hxc_dmtrans(const elecstate::DensityMatrix<TK, double>& dm_trans, const PotHxcLR& pot_hxc)
@@ -85,12 +139,11 @@ namespace LR
     template<typename TK>
     ModuleBase::matrix LR_Force<TK>::reproduce_force_gs(
         const elecstate::DensityMatrix<TK, double>& dm_gs,
-        const elecstate::DensityMatrix<TK, double>& edm_gs,
-        const elecstate::Potential& pot_gs)
+        const elecstate::DensityMatrix<TK, double>& edm_gs)
     {
         this->gint_->reset_DMRGint(PARAM.inp.nspin);
-        // Hartree+xc term
-        ModuleBase::matrix f_gs_hf_pulay = cal_force_hamilt_gs_dm_relaxed_diff(dm_gs, pot_gs); // pw+vnl+t_dphi+vl_dphi
+        // local + Hartree + xc term, including Hellmann-Feynman and Pulay
+        ModuleBase::matrix f_gs_hf_pulay = cal_force_hamilt_gs_dm_relaxed_diff(dm_gs, dm_gs); // pw+vnl+t_dphi+vl_dphi
         // edm term
         ModuleBase::matrix f_nonortho = cal_force_overlap_edm(edm_gs); // overlap
         this->gint_->reset_DMRGint(1);

source/module_lr/Grad/force/lr_force.h

@@ -39,7 +39,8 @@ namespace LR
 
         /// 1. $Tr[H_{GS}^x * (T+D^Z)]$, where GS=groud state and $(T+D^Z)$ is the relaxed difference density matrix
         ModuleBase::matrix cal_force_hamilt_gs_dm_relaxed_diff(const elecstate::DensityMatrix<TK, double>& relaxed_diff_dm,
-            const elecstate::Potential& pot_gs, const bool with_ewald = true);
+            const elecstate::DensityMatrix<TK, double>& dm_gs, const bool with_ewald = true);
+        // const elecstate::Potential& pot_gs, const bool with_ewald = true);
 
         /// 2. $Tr[S^x * (EDM)]
         ModuleBase::matrix cal_force_overlap_edm(const elecstate::DensityMatrix<TK, double>& edm);
@@ -62,8 +63,8 @@ namespace LR
         // test functions
         /// reproduce the force of the ground state
         ModuleBase::matrix reproduce_force_gs(const elecstate::DensityMatrix<TK, double>& dm_gs,
-            const elecstate::DensityMatrix<TK, double>& edm_gs,
-            const elecstate::Potential& pot_gs);
+            const elecstate::DensityMatrix<TK, double>& edm_gs);
+
         /// repreduce the ground state local term
         ModuleBase::matrix reproduce_force_gs_loc(const elecstate::DensityMatrix<TK, double>& dm_gs,
             const elecstate::Potential& pot_gs);
@@ -84,7 +85,8 @@ namespace LR
 #endif
 
         Charge dm_to_charge(const elecstate::DensityMatrix<TK, double>& dm);
-
+        elecstate::Potential dm_to_hxc_potential(const elecstate::DensityMatrix<TK, double>& dm);
+        elecstate::Potential local_potential();
         // probably move frome the ground state?
         // void build_dHS()
     };

source/module_lr/esolver_lrtd_lcao.h

@@ -148,6 +148,7 @@ namespace LR
         ct::Tensor solve_zvector_eqation(const int ispin);
         std::vector<ModuleBase::matrix> cal_force(const int ispin);
         void test_force();   // test: reproduce the force of ground state
+        elecstate::DensityMatrix<T, double> cal_dm_gs();  ///< ground-state density matrix
 
 #ifdef __EXX
         /// Tdata of Exx_LRI is same as T, for the reason, see operator_lr_exx.h