Fix: Enhance the computational stability of tau_vw.

Author: sunliang98

source/module_esolver/esolver_of_interface.cpp

@@ -193,7 +193,7 @@ void ESolver_OF::kinetic_energy_density(double** prho, double** pphi, double** r
     if (this->of_kinetic_ == "vw" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt"
         || this->of_kinetic_ == "lkt")
     {
-        this->vw_->tau_vw(prho, this->pw_rho, rtau[0]);
+        this->vw_->tau_vw(pphi, this->pw_rho, rtau[0]);
     }
     if (this->of_kinetic_ == "wt")
     {

source/module_hamilt_pw/hamilt_ofdft/kedf_vw.cpp

@@ -117,35 +117,42 @@ double KEDF_vW::get_energy_density(double** pphi, int is, int ir, ModulePW::PW_B
 
 /**
  * @brief Get the positive definite energy density of vW KEDF
- * \f[ \tau_{vW} = |\nabla \rho|^2 / (8 \rho) \f]
+ * \f[ \tau_{vW} = |\nabla \phi|^2 / 2 \f]
  * 
- * @param prho charge density
+ * @param pphi sqrt(rho)
  * @param pw_rho pw basis
  * @param rtau_vw rtau_vw => rtau_vw + tau_vw
  */
-void KEDF_vW::tau_vw(const double* const* prho, ModulePW::PW_Basis* pw_rho, double* rtau_vw)
+void KEDF_vW::tau_vw(const double* const* pphi, ModulePW::PW_Basis* pw_rho, double* rtau_vw)
 {
+    std::vector<double> abs_phi = std::vector<double>(pw_rho->nrxx, 0.);
+
     for (int is = 0; is < PARAM.inp.nspin; ++is)
     {
-        std::vector<std::vector<double>> nabla_rho(3, std::vector<double>(pw_rho->nrxx, 0.));
+        for (int ir = 0; ir < pw_rho->nrxx; ++ir)
+        {
+            abs_phi[ir] = std::abs(pphi[is][ir]);
+        }
+
+        std::vector<std::vector<double>> nabla_phi(3, std::vector<double>(pw_rho->nrxx, 0.));
+        std::vector<std::complex<double>> recip_phi(pw_rho->npw, 0.);
+        std::vector<std::complex<double>> recip_nabla_phi(pw_rho->npw, 0.);
 
-        std::vector<std::complex<double>> recip_rho(pw_rho->npw, 0.);
-        std::vector<std::complex<double>> recip_nabla_rho(pw_rho->npw, 0.);
-        pw_rho->real2recip(prho[is], recip_rho.data());
+        pw_rho->real2recip(abs_phi.data(), recip_phi.data());
 
         std::complex<double> img(0.0, 1.0);
         for (int j = 0; j < 3; ++j)
         {
             for (int ip = 0; ip < pw_rho->npw; ++ip)
             {
-                recip_nabla_rho[ip] = img * pw_rho->gcar[ip][j] * recip_rho[ip] * pw_rho->tpiba;
+                recip_nabla_phi[ip] = img * pw_rho->gcar[ip][j] * recip_phi[ip] * pw_rho->tpiba;
             }
 
-            pw_rho->recip2real(recip_nabla_rho.data(), nabla_rho[j].data());
+            pw_rho->recip2real(recip_nabla_phi.data(), nabla_phi[j].data());
 
             for (int ir = 0; ir < pw_rho->nrxx; ++ir)
             {
-                rtau_vw[ir] += nabla_rho[j][ir] * nabla_rho[j][ir] / (8. * prho[is][ir]) * 2.0; // convert Ha to Ry.
+                rtau_vw[ir] += nabla_phi[j][ir] * nabla_phi[j][ir] / 2. * this->vw_weight_ * 2.0; // convert Ha to Ry.
             }
         }
     }

source/module_hamilt_pw/hamilt_ofdft/kedf_vw.h

@@ -29,7 +29,7 @@ class KEDF_vW
 
     double get_energy(double** pphi, ModulePW::PW_Basis* pw_rho);
     double get_energy_density(double** pphi, int is, int ir, ModulePW::PW_Basis* pw_rho);
-    void tau_vw(const double* const* prho, ModulePW::PW_Basis* pw_rho, double* rtau_vw);
+    void tau_vw(const double* const* pphi, ModulePW::PW_Basis* pw_rho, double* rtau_vw);
     void vw_potential(const double* const* pphi, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpotential);
     void get_stress(const double* const* pphi, ModulePW::PW_Basis* pw_rho);