@@ -51,7 +51,8 @@ namespace LR
5151 template <typename TK >
5252 ModuleBase::matrix LR_Force<TK >::cal_force_hamilt_gs_dm_relaxed_diff(const elecstate::DensityMatrix<TK , double >& relax_diff_dm,
5353 const elecstate::DensityMatrix<TK , double >& dm_gs,
54- const bool reproduce_gs)
54+ const bool reproduce_gs,
55+ const PotHxcLR* pot_hxc_gs)
5556 {
5657 const bool with_ewald = reproduce_gs;
5758 const Charge chr_diff_relaxed = dm_to_charge (relax_diff_dm);
@@ -94,10 +95,39 @@ namespace LR
9495
9596 // 3.3 Hartree + xc (Hellmann-Feynman)
9697 ModuleBase::matrix fhxc_dvhxc (this ->ucell_ .nat , 3 );
97- elecstate::Potential pot_hxc_relaxed_diff = this ->dm_to_hxc_potential (relax_diff_dm);
98- // `cal_pulay_fs` calculates only one spin channel because `relax_diff_dm` has only one.
99- PulayForceStress::cal_pulay_fs (1 /* nspin*/
100- dm_gs, this ->ucell_ , &pot_hxc_relaxed_diff, true , false );
98+ // The potential here must be the *linear response* of $V^\text{Hxc}$ to the difference
99+ // density, i.e. $v_H[\rho^{T+Z}] + f_{xc}[\rho^\text{gs}]\,\rho^{T+Z}$ -- NOT
100+ // $v_\text{Hxc}[\rho^{T+Z}]$. This term is
101+ // $\sum_{\kappa\lambda}(T{+}D^Z)_{\kappa\lambda}\int\phi_\kappa\phi_\lambda\,
102+ // f_{xc}\sum_{\alpha\beta}D^\text{gs}_{\alpha\beta}(\phi_\alpha\phi_\beta)^x$,
103+ // the half of $\partial_x V^\text{Hxc}$ whose basis derivative falls on the *ground-state*
104+ // pair. Hartree is linear in the density so feeding it $\rho^{T+Z}$ happens to be right;
105+ // xc is not -- $\rho^{T+Z}$ is not even positive everywhere, while LDA has
106+ // $v_{xc}\propto-\rho^{1/3}$.
107+ //
108+ // `pot_hxc_gs` supplies exactly this object (Hartree weight 1, xc = $(f_{uu}+f_{ud})/2$ at
109+ // nspin=2), with no extra factor. Verified on H2/SZ TDRPA@LDA, where $K^T\equiv0$ makes the
110+ // triplet gradient identical to $d(\varepsilon_a-\varepsilon_i)/dx$: analytic 28.5982 vs the
111+ // KS-gap finite difference 28.598156. It used to be off by -3.5 eV/Ang.
112+ //
113+ // `reproduce_gs` is the exception: there `relax_diff_dm` *is* the ground-state density
114+ // matrix and the term being checked is the true ground-state force, for which
115+ // $v_\text{Hxc}[\rho^\text{gs}]$ is the correct potential.
116+ if (reproduce_gs || pot_hxc_gs == nullptr )
117+ {
118+ elecstate::Potential pot_hxc_relaxed_diff = this ->dm_to_hxc_potential (relax_diff_dm);
119+ // `cal_pulay_fs` calculates only one spin channel because `relax_diff_dm` has only one.
120+ PulayForceStress::cal_pulay_fs (1 /* nspin*/
121+ dm_gs, this ->ucell_ , &pot_hxc_relaxed_diff, true , false );
122+ }
123+ else
124+ {
125+ ModuleBase::matrix v_lin (1 , this ->rhopw_ .nrxx ); // zero-initialized
126+ double * rho_in[1 ] = { const_cast <double *>(chr_diff_relaxed.rho [0 ]) };
127+ pot_hxc_gs->cal_v_eff (rho_in, this ->ucell_ , v_lin);
128+ std::vector<const double *> vr_eff = { v_lin.c };
129+ ModuleGint::cal_gint_fvl (1 , vr_eff, dm_gs.get_DMR_vector (), true , false , &fhxc_dvhxc, &stress_tmp);
130+ }
101131 if (!reproduce_gs) {fhxc_dvhxc *= 2 ;} // for the two channels of the ground-state dm.
102132
103133 // 4. kinetic (Pulay)
0 commit comments