From d6c50472e0538ead7f6b083e58bb6efc882a1679 Mon Sep 17 00:00:00 2001
From: sunliang98 <1700011430@pku.edu.cn>
Date: Fri, 4 Jul 2025 17:04:37 +0800
Subject: [PATCH 1/4] Feature: Implement Xu-Wang-Ma KEDF

---
 source/Makefile.Objects                       |   1 +
 source/module_parameter/input_parameter.h     |   2 +
 source/source_io/read_input_item_ofdft.cpp    |  16 +-
 source/source_pw/hamilt_ofdft/CMakeLists.txt  |   1 +
 .../source_pw/hamilt_ofdft/kedf_manager.cpp   |  64 ++-
 source/source_pw/hamilt_ofdft/kedf_manager.h  |   2 +
 source/source_pw/hamilt_ofdft/kedf_xwm.cpp    | 381 ++++++++++++++++++
 source/source_pw/hamilt_ofdft/kedf_xwm.h      |  59 +++
 8 files changed, 510 insertions(+), 16 deletions(-)
 create mode 100644 source/source_pw/hamilt_ofdft/kedf_xwm.cpp
 create mode 100644 source/source_pw/hamilt_ofdft/kedf_xwm.h

diff --git a/source/Makefile.Objects b/source/Makefile.Objects
index f13d32f777..1412cd7d75 100644
--- a/source/Makefile.Objects
+++ b/source/Makefile.Objects
@@ -345,6 +345,7 @@ OBJS_HAMILT=hamilt_pw.o\
 OBJS_HAMILT_OF=kedf_tf.o\
     kedf_vw.o\
     kedf_wt.o\
+    kedf_xwm.o\
     kedf_lkt.o\
     kedf_manager.o\
 
diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h
index 16db768542..e137fbea9e 100644
--- a/source/module_parameter/input_parameter.h
+++ b/source/module_parameter/input_parameter.h
@@ -206,6 +206,8 @@ struct Input_para
                                                  ///< filled from file of_kernel_file, not from
                                                  ///< formula. Only usable for WT KEDF.
     std::string of_kernel_file = "WTkernel.txt"; ///< The name of WT kernel file.
+    double of_xwm_kappa = 0.0;                   ///< The parameter kappa of XWM KEDF
+    double of_xwm_rho_ref = 0.0;                 ///< The reference density of XWM KEDF
 
     // ML KEDF, sunliang added on 2022-11-07
     bool of_ml_gene_data = false; ///< Generate training data or not
diff --git a/source/source_io/read_input_item_ofdft.cpp b/source/source_io/read_input_item_ofdft.cpp
index 9f729d2793..794bbc4d30 100644
--- a/source/source_io/read_input_item_ofdft.cpp
+++ b/source/source_io/read_input_item_ofdft.cpp
@@ -18,10 +18,10 @@ void ReadInput::item_ofdft()
             }
 #endif
             if (para.input.of_kinetic != "tf" && para.input.of_kinetic != "vw" && para.input.of_kinetic != "wt"
-                && para.input.of_kinetic != "lkt" && para.input.of_kinetic != "tf+" 
+                && para.input.of_kinetic != "xwm" && para.input.of_kinetic != "lkt" && para.input.of_kinetic != "tf+" 
                 && para.input.of_kinetic != "ml" && para.input.of_kinetic != "mpn" && para.input.of_kinetic != "cpn5")
             {
-                ModuleBase::WARNING_QUIT("ReadInput", "of_kinetic must be tf, vw, tf+, wt, lkt, ml, mpn, or cpn5");
+                ModuleBase::WARNING_QUIT("ReadInput", "of_kinetic must be tf, vw, tf+, wt, xwm, lkt, ml, mpn, or cpn5");
             }
         };
         item.reset_value = [](const Input_Item& item, Parameter& para) {
@@ -209,6 +209,18 @@ void ReadInput::item_ofdft()
         read_sync_string(input.of_kernel_file);
         this->add_item(item);
     }
+    {
+        Input_Item item("of_xwm_kappa");
+        item.annotation = "The parameter kappa of XWM KEDF";
+        read_sync_double(input.of_xwm_kappa);
+        this->add_item(item);
+    }
+    {
+        Input_Item item("of_xwm_rho_ref");
+        item.annotation = "The reference density of XWM KEDF";
+        read_sync_double(input.of_xwm_rho_ref);
+        this->add_item(item);
+    }
     {
         Input_Item item("of_ml_gene_data");
         item.annotation = "Generate training data or not";
diff --git a/source/source_pw/hamilt_ofdft/CMakeLists.txt b/source/source_pw/hamilt_ofdft/CMakeLists.txt
index 7a98c732b4..fe6eec0d79 100644
--- a/source/source_pw/hamilt_ofdft/CMakeLists.txt
+++ b/source/source_pw/hamilt_ofdft/CMakeLists.txt
@@ -2,6 +2,7 @@ list(APPEND hamilt_ofdft_srcs
     kedf_tf.cpp
     kedf_vw.cpp
     kedf_wt.cpp
+    kedf_xwm.cpp
     kedf_lkt.cpp
     kedf_manager.cpp
     of_stress_pw.cpp
diff --git a/source/source_pw/hamilt_ofdft/kedf_manager.cpp b/source/source_pw/hamilt_ofdft/kedf_manager.cpp
index 3b1e7f6cc4..401d843baf 100644
--- a/source/source_pw/hamilt_ofdft/kedf_manager.cpp
+++ b/source/source_pw/hamilt_ofdft/kedf_manager.cpp
@@ -18,11 +18,12 @@ void KEDF_Manager::init(
 {
     this->of_kinetic_ = inp.of_kinetic;
 
-    //! Thomas-Fermi (TF) KEDF, TF+ KEDF, and Want-Teter (WT) KEDF
+    //! Thomas-Fermi (TF) KEDF, TF+ KEDF, Want-Teter (WT) KEDF, and XWM KEDF
     if (this->of_kinetic_ == "tf"
      || this->of_kinetic_ == "tf+"
      || this->of_kinetic_ == "wt"
-     || this->of_kinetic_ == "ml")
+     || this->of_kinetic_ == "ml"
+     || this->of_kinetic_ == "xwm")
     {
         if (this->tf_ == nullptr)
         {
@@ -31,9 +32,9 @@ void KEDF_Manager::init(
         this->tf_->set_para(pw_rho->nrxx, dV, inp.of_tf_weight);
     }
 
-    //! vW, TF+, WT, and LKT KEDFs
+    //! vW, TF+, WT, XWM, and LKT KEDFs
     if (this->of_kinetic_ == "vw" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt"
-        || this->of_kinetic_ == "lkt" || this->of_kinetic_ == "ml")
+        || this->of_kinetic_ == "xwm" || this->of_kinetic_ == "lkt" || this->of_kinetic_ == "ml")
     {
         if (this->vw_ == nullptr)
         {
@@ -62,6 +63,17 @@ void KEDF_Manager::init(
                             pw_rho);
     }
 
+    //! Xu-Wang-Ma KEDF
+    if (this->of_kinetic_ == "xwm")
+    {
+        if (this->xwm_ == nullptr)
+        {
+            this->xwm_ = new KEDF_XWM();
+        }
+        this->xwm_->set_para(dV, inp.of_xwm_rho_ref, inp.of_xwm_kappa, nelec,
+                            inp.of_tf_weight, inp.of_vw_weight, pw_rho);
+    }
+
     //! LKT KEDF
     if (this->of_kinetic_ == "lkt")
     {
@@ -108,7 +120,7 @@ void KEDF_Manager::get_potential(
     if (PARAM.inp.of_ml_local_test) this->ml_->localTest(prho, pw_rho);
 #endif
 
-    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt")
+    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt" || this->of_kinetic_ == "xwm")
     {
         this->tf_->tf_potential(prho, rpot);
     }
@@ -116,6 +128,10 @@ void KEDF_Manager::get_potential(
     {
         this->wt_->wt_potential(prho, pw_rho, rpot);
     }
+    if (this->of_kinetic_ == "xwm")
+    {
+        this->xwm_->xwm_potential(prho, pw_rho, rpot);
+    }
     if (this->of_kinetic_ == "lkt")
     {
         this->lkt_->lkt_potential(prho, pw_rho, rpot);
@@ -138,7 +154,7 @@ void KEDF_Manager::get_potential(
     }
 
     if (this->of_kinetic_ == "vw" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt"
-        || this->of_kinetic_ == "lkt" || this->of_kinetic_ == "ml")
+        || this->of_kinetic_ == "xwm" || this->of_kinetic_ == "lkt" || this->of_kinetic_ == "ml")
     {
         this->vw_->vw_potential(pphi, pw_rho, rpot);
     }
@@ -154,13 +170,13 @@ double KEDF_Manager::get_energy()
 {
     double kinetic_energy = 0.0;
 
-    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt")
+    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt"|| this->of_kinetic_ == "xwm")
     {
         kinetic_energy += this->tf_->tf_energy;
     }
 
     if (this->of_kinetic_ == "vw" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt"
-        || this->of_kinetic_ == "lkt" || this->of_kinetic_ == "ml")
+        || this->of_kinetic_ == "xwm" || this->of_kinetic_ == "lkt" || this->of_kinetic_ == "ml")
     {
         kinetic_energy += this->vw_->vw_energy;
     }
@@ -170,6 +186,11 @@ double KEDF_Manager::get_energy()
         kinetic_energy += this->wt_->wt_energy;
     }
 
+    if (this->of_kinetic_ == "xwm")
+    {
+        kinetic_energy += this->xwm_->xwm_energy;
+    }
+
     if (this->of_kinetic_ == "lkt")
     {
         kinetic_energy += this->lkt_->lkt_energy;
@@ -210,12 +231,12 @@ void KEDF_Manager::get_energy_density(
         rtau[0][ir] = 0.0;
     }
 
-    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt")
+    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt" || this->of_kinetic_ == "xwm")
     {
         this->tf_->tau_tf(prho, rtau[0]);
     }
     if (this->of_kinetic_ == "vw" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt"
-        || this->of_kinetic_ == "lkt")
+        || this->of_kinetic_ == "xwm" || this->of_kinetic_ == "lkt")
     {
         this->vw_->tau_vw(pphi, pw_rho, rtau[0]);
     }
@@ -223,6 +244,10 @@ void KEDF_Manager::get_energy_density(
     {
         this->wt_->tau_wt(prho, pw_rho, rtau[0]);
     }
+    if (this->of_kinetic_ == "xwm")
+    {
+        this->xwm_->tau_xwm(prho, pw_rho, rtau[0]);
+    }
     if (this->of_kinetic_ == "lkt")
     {
         this->lkt_->tau_lkt(prho, pw_rho, rtau[0]);
@@ -255,14 +280,14 @@ void KEDF_Manager::get_stress(
         }
     }
 
-    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt")
+    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt" || this->of_kinetic_ == "xwm")
     {
         this->tf_->get_stress(omega);
         kinetic_stress_ += this->tf_->stress;
     }
 
     if (this->of_kinetic_ == "vw" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt"
-        || this->of_kinetic_ == "lkt")
+        || this->of_kinetic_ == "xwm" || this->of_kinetic_ == "lkt")
     {
         this->vw_->get_stress(pphi, pw_rho);
         kinetic_stress_ += this->vw_->stress;
@@ -274,6 +299,12 @@ void KEDF_Manager::get_stress(
         kinetic_stress_ += this->wt_->stress;
     }
 
+    if (this->of_kinetic_ == "xwm")
+    {
+        this->xwm_->get_stress(prho, pw_rho, PARAM.inp.of_vw_weight);
+        kinetic_stress_ += this->xwm_->stress;
+    }
+
     if (this->of_kinetic_ == "lkt")
     {
         this->lkt_->get_stress(prho, pw_rho);
@@ -290,13 +321,13 @@ void KEDF_Manager::record_energy(
     std::vector<double> &energies_Ry
 )
 {
-    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt")
+    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt" || this->of_kinetic_ == "xwm")
     {
         titles.push_back("TF KEDF");
         energies_Ry.push_back(this->tf_->tf_energy);
     }
     if (this->of_kinetic_ == "vw" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt"
-        || this->of_kinetic_ == "lkt" || this->of_kinetic_ == "ml")
+        || this->of_kinetic_ == "xwm" || this->of_kinetic_ == "lkt" || this->of_kinetic_ == "ml")
     {
         titles.push_back("vW KEDF");
         energies_Ry.push_back(this->vw_->vw_energy);
@@ -306,6 +337,11 @@ void KEDF_Manager::record_energy(
         titles.push_back("WT KEDF");
         energies_Ry.push_back(this->wt_->wt_energy);
     }
+    if (this->of_kinetic_ == "xwm")
+    {
+        titles.push_back("XWM KEDF");
+        energies_Ry.push_back(this->xwm_->xwm_energy);
+    }
     if (this->of_kinetic_ == "lkt")
     {
         titles.push_back("LKT KEDF");
diff --git a/source/source_pw/hamilt_ofdft/kedf_manager.h b/source/source_pw/hamilt_ofdft/kedf_manager.h
index 92b2fb0114..b9b311832f 100644
--- a/source/source_pw/hamilt_ofdft/kedf_manager.h
+++ b/source/source_pw/hamilt_ofdft/kedf_manager.h
@@ -8,6 +8,7 @@
 #include "source_pw/hamilt_ofdft/kedf_tf.h"
 #include "source_pw/hamilt_ofdft/kedf_vw.h"
 #include "source_pw/hamilt_ofdft/kedf_wt.h"
+#include "source_pw/hamilt_ofdft/kedf_xwm.h"
 #include "source_pw/hamilt_ofdft/kedf_ml.h"
 
 class KEDF_Manager
@@ -74,6 +75,7 @@ class KEDF_Manager
     KEDF_TF* tf_ = nullptr;   // Thomas-Fermi KEDF
     KEDF_vW* vw_ = nullptr;   // von Weizsäcker KEDF
     KEDF_WT* wt_ = nullptr;   // Wang-Teter KEDF
+    KEDF_XWM* xwm_ = nullptr; // Xu-Wang-Ma KEDF
 #ifdef __MLALGO
     KEDF_ML* ml_ = nullptr;   // Machine Learning KEDF
 #endif
diff --git a/source/source_pw/hamilt_ofdft/kedf_xwm.cpp b/source/source_pw/hamilt_ofdft/kedf_xwm.cpp
new file mode 100644
index 0000000000..55ccfc92eb
--- /dev/null
+++ b/source/source_pw/hamilt_ofdft/kedf_xwm.cpp
@@ -0,0 +1,381 @@
+#include "./kedf_xwm.h"
+
+#include "module_parameter/parameter.h"
+#include "source_base/parallel_reduce.h"
+#include "source_base/tool_quit.h"
+
+/**
+ * @brief Set the parameters of XWM KEDF, and initialize kernel
+ *
+ * @param dV the volume of one grid point in real space, omega/nxyz
+ * @param rho_ref
+ * @param kappa
+ * @param nelec the number of electron
+ * @param tf_weight
+ * @param vw_weight
+ * @param pw_rho pw_basis
+ */
+void KEDF_XWM::set_para(double dV,
+                double rho_ref,
+                double kappa,
+                double nelec,
+                double tf_weight,
+                double vw_weight,
+                ModulePW::PW_Basis* pw_rho)
+{
+    this->dV_ = dV;
+    this->kappa_ = kappa;
+    this->kappa_5_6 = this->kappa_ + 5. / 6.;
+    this->kappa_11_6 = this->kappa_ + 11. / 6.;
+    this->kappa_1_6 = this->kappa_ - 1. / 6.;
+
+    this->rho0_ = 1. / (pw_rho->nxyz * dV) * nelec;
+    if (rho_ref != 0)
+    {
+        this->rho_ref_ = rho_ref;
+    }
+    else
+    {
+        this->rho_ref_ = this->rho0_;
+    }
+    this->c_kernel /= std::pow(this->rho0_, 2. * this->kappa_);
+
+    this->kf_ = std::pow(3. * std::pow(ModuleBase::PI, 2) * this->rho_ref_, 1. / 3.);
+    this->tkf_ = 2. * this->kf_;
+
+    double temp1 = 1. / (this->kappa_ + 5. / 6.);
+    double temp2 = 1. / (this->kappa_ + 11. / 6.);
+    this->c_0 = 0.5 * temp1 * temp1;
+    this->c_1 = temp1 * temp2;
+    this->c_2 = - temp1 * temp1 * this->rho_ref_;
+
+    this->kernel1_.resize(pw_rho->npw, 0.);
+    this->kernel2_.resize(pw_rho->npw, 0.);
+
+    this->fill_kernel(tf_weight, vw_weight, pw_rho);
+}
+
+/**
+ * @brief Get the energy of XWM KEDF
+ * 
+ * @param prho charge density
+ * @param pw_rho pw basis
+ * @return the energy of XWM KEDF
+ */
+double KEDF_XWM::get_energy(const double* const* prho, ModulePW::PW_Basis* pw_rho)
+{
+    double** w1Rho5_6 = new double*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        w1Rho5_6[is] = new double[pw_rho->nrxx];
+    }
+    this->multi_kernel(prho, this->kernel1_.data(), w1Rho5_6, this->kappa_5_6, pw_rho);
+
+    double** w2Rho5_6 = new double*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        w2Rho5_6[is] = new double[pw_rho->nrxx];
+    }
+    this->multi_kernel(prho, this->kernel2_.data(), w2Rho5_6, this->kappa_5_6, pw_rho);
+
+    double energy = 0.; // in Ry
+    if (PARAM.inp.nspin == 1)
+    {
+        for (int ir = 0; ir < pw_rho->nrxx; ++ir)
+        {
+            energy += std::pow(prho[0][ir], this->kappa_5_6) * w1Rho5_6[0][ir]
+                    + std::pow(prho[0][ir], this->kappa_11_6) * w2Rho5_6[0][ir];
+        }
+        energy += this->dV_;
+    }
+    else if (PARAM.inp.nspin == 2)
+    {
+        // TODO: spin polarized
+    }
+    this->xwm_energy = energy;
+    Parallel_Reduce::reduce_all(this->xwm_energy);
+
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        delete[] w1Rho5_6[is];
+        delete[] w2Rho5_6[is];
+    }
+    delete[] w1Rho5_6;
+    delete[] w2Rho5_6;
+
+    return energy;
+}
+
+/**
+ * @brief Get the energy density of XWM KEDF
+ * 
+ * @param prho charge density
+ * @param is spin index
+ * @param ir grid index
+ * @param pw_rho pw basis
+ * @return the energy density of XWM KEDF
+ */
+double KEDF_XWM::get_energy_density(const double* const* prho, int is, int ir, ModulePW::PW_Basis* pw_rho)
+{
+    double** w1Rho5_6 = new double*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        w1Rho5_6[is] = new double[pw_rho->nrxx];
+    }
+    this->multi_kernel(prho, this->kernel1_.data(), w1Rho5_6, this->kappa_5_6, pw_rho);
+
+    double** w2Rho5_6 = new double*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        w2Rho5_6[is] = new double[pw_rho->nrxx];
+    }
+    this->multi_kernel(prho, this->kernel2_.data(), w2Rho5_6, this->kappa_5_6, pw_rho);
+
+    double result = std::pow(prho[is][ir], this->kappa_5_6) * w1Rho5_6[is][ir]
+                  + std::pow(prho[is][ir], this->kappa_11_6) * w2Rho5_6[is][ir];
+    
+    result *= this->dV_;
+
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        delete[] w1Rho5_6[is];
+        delete[] w2Rho5_6[is];
+    }
+    delete[] w1Rho5_6;
+    delete[] w2Rho5_6;
+
+    return result;
+}
+
+/**
+ * @brief Get the kinetic energy of XWM KEDF, and add it onto rtau_xwm
+ * 
+ * @param prho charge density
+ * @param pw_rho pw basis
+ * @param rtau_xwm rtau_xwm => rtau_xwm + tau_xwm
+ */
+void KEDF_XWM::tau_xwm(const double* const* prho, ModulePW::PW_Basis* pw_rho, double* rtau_xwm)
+{
+    double** w1Rho5_6 = new double*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        w1Rho5_6[is] = new double[pw_rho->nrxx];
+    }
+    this->multi_kernel(prho, this->kernel1_.data(), w1Rho5_6, this->kappa_5_6, pw_rho);
+
+    double** w2Rho5_6 = new double*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        w2Rho5_6[is] = new double[pw_rho->nrxx];
+    }
+    this->multi_kernel(prho, this->kernel2_.data(), w2Rho5_6, this->kappa_5_6, pw_rho);
+
+    if (PARAM.inp.nspin == 1)
+    {
+        for (int ir = 0; ir < pw_rho->nrxx; ++ir)
+        {
+            rtau_xwm[ir] += std::pow(prho[0][ir], this->kappa_5_6) * w1Rho5_6[0][ir]
+                          + std::pow(prho[0][ir], this->kappa_11_6) * w2Rho5_6[0][ir];
+        }
+    }
+    else if (PARAM.inp.nspin == 2)
+    {
+        // TODO: spin polarized
+    }
+
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        delete[] w1Rho5_6[is];
+        delete[] w2Rho5_6[is];
+    }
+    delete[] w1Rho5_6;
+    delete[] w2Rho5_6;
+}
+
+/**
+ * @brief Get the potential of XWM KEDF, and add it into rpotential,
+ * and the XWM energy is stored in this->xwm_energy
+ * 
+ * @param prho charge density
+ * @param pw_rho pw basis
+ * @param rpotential rpotential => rpotential + V_{XWM}
+ */
+void KEDF_XWM::xwm_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpotential)
+{
+    ModuleBase::timer::tick("KEDF_XWM", "xwm_potential");
+    double** w1Rho5_6 = new double*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        w1Rho5_6[is] = new double[pw_rho->nrxx];
+    }
+    this->multi_kernel(prho, this->kernel1_.data(), w1Rho5_6, this->kappa_5_6, pw_rho);
+
+    double** w2Rho11_6 = new double*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        w2Rho11_6[is] = new double[pw_rho->nrxx];
+    }
+    this->multi_kernel(prho, this->kernel2_.data(), w2Rho11_6, this->kappa_11_6, pw_rho);
+
+    double** w2Rho5_6 = new double*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        w2Rho5_6[is] = new double[pw_rho->nrxx];
+    }
+    this->multi_kernel(prho, this->kernel2_.data(), w2Rho5_6, this->kappa_5_6, pw_rho);
+
+    double energy = 0.;
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        for (int ir = 0; ir < pw_rho->nrxx; ++ir)
+        {
+            double rho1_6 = std::pow(prho[is][ir], this->kappa_1_6);
+            double rho5_6 = std::pow(prho[is][ir], this->kappa_5_6);
+            double rho11_6 = std::pow(prho[is][ir], this->kappa_11_6);
+
+            rpotential(is, ir) += 2. * this->kappa_5_6 * rho1_6 * w1Rho5_6[is][ir]
+                                + this->kappa_5_6 * rho1_6 * w2Rho11_6[is][ir]
+                                + this->kappa_11_6 * rho5_6 * w2Rho5_6[is][ir];
+
+            energy += rho5_6 * w1Rho5_6[is][ir] + rho11_6 * w2Rho5_6[is][ir]; // FOR SPIN 1 !!!
+        }
+    }
+    energy *= this->dV_;
+    this->xwm_energy = energy;
+    Parallel_Reduce::reduce_all(this->xwm_energy);
+
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        delete[] w1Rho5_6[is];
+        delete[] w2Rho11_6[is];
+        delete[] w2Rho5_6[is];
+    }
+    delete[] w1Rho5_6;
+    delete[] w2Rho11_6;
+    delete[] w2Rho5_6;
+    ModuleBase::timer::tick("KEDF_XWM", "xwm_potential");
+}
+
+/**
+ * @brief Get the stress of XWM KEDF, and store it into this->stress
+ * 
+ * @param prho charge density
+ * @param pw_rho pw basis
+ * @param vw_weight the weight of vW KEDF
+ */
+void KEDF_XWM::get_stress(const double* const* prho, ModulePW::PW_Basis* pw_rho, double vw_weight)
+{
+    std::cout << "XWM stress is not implemented yet!" << std::endl;
+}
+
+/**
+ * @brief Calculate \int{W(r-r')rho^{exponent}(r') dr'}
+ *
+ * @param [in] prho charge density
+ * @param [in] kernel W(r-r')
+ * @param [out] rkernel_rho \int{W(r-r')rho^{exponent}(r') dr'}
+ * @param [in] exponent the exponent of rho
+ * @param [in] pw_rho pw_basis
+ */
+void KEDF_XWM::multi_kernel(const double* const* prho, const double* kernel, double** rkernel_rho, double exponent, ModulePW::PW_Basis* pw_rho)
+{
+    std::complex<double>** recipkernelRho = new std::complex<double>*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        recipkernelRho[is] = new std::complex<double>[pw_rho->npw];
+        for (int ir = 0; ir < pw_rho->nrxx; ++ir)
+        {
+            rkernel_rho[is][ir] = std::pow(prho[is][ir], exponent);
+        }
+        pw_rho->real2recip(rkernel_rho[is], recipkernelRho[is]);
+        for (int ip = 0; ip < pw_rho->npw; ++ip)
+        {
+            recipkernelRho[is][ip] *= kernel[ip];
+        }
+        pw_rho->recip2real(recipkernelRho[is], rkernel_rho[is]);
+    }
+
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        delete[] recipkernelRho[is];
+    }
+    delete[] recipkernelRho;
+}
+
+/**
+ * @brief Fill the kernel (this->kernel_)
+ *
+ * @param tf_weight
+ * @param vw_weight
+ * @param pw_rho pw_basis
+ */
+void KEDF_XWM::fill_kernel(double tf_weight, double vw_weight, ModulePW::PW_Basis* pw_rho)
+{
+    double eta = 0.;
+    double coef = - 1./6. / this->rho_ref_;
+    for (int ig = 0; ig < pw_rho->npw; ++ig)
+    {
+        eta = sqrt(pw_rho->gg[ig]) * pw_rho->tpiba / this->tkf_;
+
+        // calculate the lindhard response function and its derivative
+        double lindhard = 0.;
+        double diff_lindhard = 0.;
+        if (eta < 0.)
+        {
+            lindhard = 0.;
+            diff_lindhard = 0.;
+        }
+        // limit for small eta
+        else if (eta < 1e-10)
+        {
+            lindhard = 1. - tf_weight + eta * eta * (1. / 3. - 3. * vw_weight);
+            diff_lindhard = 2. * eta * (1. / 3. - 3. * vw_weight);
+        }
+        // around the singularity
+        else if (std::abs(eta - 1.) < 1e-10)
+        {
+            lindhard = 2. - tf_weight - 3. * vw_weight + 40. * (eta - 1);
+            diff_lindhard = 40.;
+        }
+        // Taylor expansion for high eta
+        else if (eta > 3.65)
+        {
+            double eta2 = eta * eta;
+            double invEta2 = 1. / eta2;
+            lindhard = 3. * (1. - vw_weight) * eta2
+                            -tf_weight-0.6
+                            + invEta2 * (-0.13714285714285712
+                            + invEta2 * (-6.39999999999999875E-2
+                            + invEta2 * (-3.77825602968460128E-2
+                            + invEta2 * (-2.51824061652633074E-2
+                            + invEta2 * (-1.80879839616166146E-2
+                            + invEta2 * (-1.36715733124818332E-2
+                            + invEta2 * (-1.07236045520990083E-2
+                            + invEta2 * (-8.65192783339199453E-3
+                            + invEta2 * (-7.1372762502456763E-3
+                            + invEta2 * (-5.9945117538835746E-3
+                            + invEta2 * (-5.10997527675418131E-3
+                            + invEta2 * (-4.41060829979912465E-3
+                            + invEta2 * (-3.84763737842981233E-3
+                            + invEta2 * (-3.38745061493813488E-3
+                            + invEta2 * (-3.00624946457977689E-3)))))))))))))));
+            diff_lindhard = ((eta2 + 1.) * 0.25 / eta2 * std::log(std::abs((1. + eta) / (1. - eta))) - 0.5 / eta)
+                            / std::pow((0.5 + 0.25 * (1. - eta2) * std::log((1. + eta) / std::abs(1. - eta)) / eta), 2)
+                            - 6. * eta * vw_weight;
+        }
+        else
+        {
+            double eta2 = eta * eta;
+            double log_term = std::log(std::abs((1. + eta) / (1. - eta)));
+            double term1 = 1. / (0.5 + 0.25 * (1. - eta2) / eta * log_term);
+            double term2 = 0.25 * (1. + eta2) / eta2 * log_term - 0.5 / eta;
+
+            lindhard = term1 - 3. * vw_weight * eta2 - tf_weight;
+            diff_lindhard = term1 * term1 * term2 - 6. * eta * vw_weight;
+        }
+
+        diff_lindhard = coef * eta * diff_lindhard;
+        this->kernel1_[ig] = (this->c_0 * lindhard + this->c_2 * diff_lindhard) * this->c_kernel;
+        this->kernel2_[ig] = this->c_1 * diff_lindhard * this->c_kernel;
+    }
+}
\ No newline at end of file
diff --git a/source/source_pw/hamilt_ofdft/kedf_xwm.h b/source/source_pw/hamilt_ofdft/kedf_xwm.h
new file mode 100644
index 0000000000..cc4ef2dbff
--- /dev/null
+++ b/source/source_pw/hamilt_ofdft/kedf_xwm.h
@@ -0,0 +1,59 @@
+#ifndef KEDF_XWM_H
+#define KEDF_XWM_H
+
+#include "source_base/matrix.h"
+#include "source_base/timer.h"
+#include "source_basis/module_pw/pw_basis.h"
+
+/**
+ * @brief A class which calculates the kinetic energy, potential, and stress with Xie-Wang-Morales (XWM) KEDF.
+ * See Xu Q, Wang Y, Ma Y. Physical Review B, 2019, 100(20): 205132.
+ * @author sunliang on 2025-01
+ */
+class KEDF_XWM
+{
+  public:
+    KEDF_XWM()
+    {
+        this->stress.create(3, 3);
+    }
+    ~KEDF_XWM(){}
+
+    void set_para(double dV,
+                  double rho_ref,
+                  double kappa,
+                  double nelec,
+                  double tf_weight,
+                  double vw_weight,
+                  ModulePW::PW_Basis* pw_rho);
+
+
+    double get_energy(const double* const* prho, ModulePW::PW_Basis* pw_rho);
+    double get_energy_density(const double* const* prho, int is, int ir, ModulePW::PW_Basis* pw_rho);
+    void tau_xwm(const double* const* prho, ModulePW::PW_Basis* pw_rho, double* rtau_xwm);
+    void xwm_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpotential);
+    void get_stress(const double* const* prho, ModulePW::PW_Basis* pw_rho, double vw_weight);
+    double xwm_energy = 0.;
+    ModuleBase::matrix stress;
+
+  private:
+    void multi_kernel(const double* const* prho, const double* kernel, double** rkernel_rho, double exponent, ModulePW::PW_Basis* pw_rho);
+    void fill_kernel(double tf_weight, double vw_weight, ModulePW::PW_Basis* pw_rho);
+
+    double dV_ = 0.;
+    double rho0_ = 0.; // average rho
+    double rho_ref_ = 0.; // reference rho
+    double kf_ = 0.;  // Fermi vector kF = (3 pi^2 rho_star_)^(1/3)
+    double tkf_ = 0.; // 2 * kF
+    double kappa_ = 0.;
+    double c_kernel = std::pow(ModuleBase::PI, 4./3.) / std::pow(3., 1. / 3.) * 2.; // multiply by 2 to convert unit from Hartree to Ry
+    double c_0 = 0.; // coef of T0, c_0 = 18 / (6 * kappa + 5) ^ 2
+    double c_1 = 0.; // coef of T1, c_1 = [(kappa + 5/6)(kappa + 11/6)] ^ -1
+    double c_2 = 0.; // coef of T1, c_2 = - (kappa + 5/6) ^ -2 * rho_ref_
+    double kappa_5_6 = 0.; // kappa + 5/6
+    double kappa_11_6 = 0.; // kappa + 11/6
+    double kappa_1_6 = 0.; // kappa - 1/6
+    std::vector<double> kernel1_ = {}; // w1 in PHYSICAL REVIEW B 110, 085113 (2024)
+    std::vector<double> kernel2_ = {}; // w2 in PHYSICAL REVIEW B 110, 085113 (2024)
+};
+#endif
\ No newline at end of file

From 6eb3d61c07e8a51bd6626677670029a22736e88b Mon Sep 17 00:00:00 2001
From: sunliang98 <1700011430@pku.edu.cn>
Date: Fri, 4 Jul 2025 17:07:40 +0800
Subject: [PATCH 2/4] Test: Add two integrate tests for XWM KEDF

---
 source/source_io/test/read_input_ptest.cpp |  2 ++
 source/source_io/test/support/INPUT        |  2 ++
 tests/07_OFDFT/28_OF_KE_XWM/INPUT          | 26 ++++++++++++++++++++
 tests/07_OFDFT/28_OF_KE_XWM/KPT            |  4 ++++
 tests/07_OFDFT/28_OF_KE_XWM/README         |  1 +
 tests/07_OFDFT/28_OF_KE_XWM/STRU           | 18 ++++++++++++++
 tests/07_OFDFT/28_OF_KE_XWM/result.ref     |  7 ++++++
 tests/07_OFDFT/29_OF_XWM_para/INPUT        | 28 ++++++++++++++++++++++
 tests/07_OFDFT/29_OF_XWM_para/KPT          |  4 ++++
 tests/07_OFDFT/29_OF_XWM_para/README       |  1 +
 tests/07_OFDFT/29_OF_XWM_para/STRU         | 18 ++++++++++++++
 tests/07_OFDFT/29_OF_XWM_para/result.ref   |  7 ++++++
 tests/07_OFDFT/CASES_CPU.txt               |  2 ++
 13 files changed, 120 insertions(+)
 create mode 100644 tests/07_OFDFT/28_OF_KE_XWM/INPUT
 create mode 100644 tests/07_OFDFT/28_OF_KE_XWM/KPT
 create mode 100644 tests/07_OFDFT/28_OF_KE_XWM/README
 create mode 100644 tests/07_OFDFT/28_OF_KE_XWM/STRU
 create mode 100644 tests/07_OFDFT/28_OF_KE_XWM/result.ref
 create mode 100644 tests/07_OFDFT/29_OF_XWM_para/INPUT
 create mode 100644 tests/07_OFDFT/29_OF_XWM_para/KPT
 create mode 100644 tests/07_OFDFT/29_OF_XWM_para/README
 create mode 100644 tests/07_OFDFT/29_OF_XWM_para/STRU
 create mode 100644 tests/07_OFDFT/29_OF_XWM_para/result.ref

diff --git a/source/source_io/test/read_input_ptest.cpp b/source/source_io/test/read_input_ptest.cpp
index 5755934032..35a6142c6c 100644
--- a/source/source_io/test/read_input_ptest.cpp
+++ b/source/source_io/test/read_input_ptest.cpp
@@ -365,6 +365,8 @@ TEST_F(InputParaTest, ParaRead)
     EXPECT_EQ(param.inp.of_full_pw_dim, 0);
     EXPECT_FALSE(param.inp.of_read_kernel);
     EXPECT_EQ(param.inp.of_kernel_file, "WTkernel.txt");
+    EXPECT_DOUBLE_EQ(param.inp.of_xwm_kappa, 1.);
+    EXPECT_DOUBLE_EQ(param.inp.of_xwm_rho_ref, 1.);
     EXPECT_EQ(param.inp.device, "cpu");
     EXPECT_NEAR(param.inp.force_thr_ev, 0.025711245953622324, 1e-8);
     EXPECT_DOUBLE_EQ(param.globalv.hubbard_u[0], 0);
diff --git a/source/source_io/test/support/INPUT b/source/source_io/test/support/INPUT
index f4814d34c2..db0146a4da 100644
--- a/source/source_io/test/support/INPUT
+++ b/source/source_io/test/support/INPUT
@@ -369,6 +369,8 @@ of_full_pw                     0 #If set to 1, ecut will be ignored when collect
 of_full_pw_dim                 0 #If of_full_pw = true, dimention of FFT is testricted to be (0) either odd or even; (1) odd only; (2) even only
 of_read_kernel                 0 #If set to 1, the kernel of WT KEDF will be filled from file of_kernel_file, not from formula. Only usable for WT KEDF
 of_kernel_file                 WTkernel.txt #The name of WT kernel file.
+of_xwm_kappa                   1 #The parameter kappa of XWM KEDF
+of_xwm_rho_ref                 1 #The reference density of XWM KEDF
 
 #Parameters (20.dft+u)
 dft_plus_u                     0 #true:DFT+U correction; false: standard DFT calcullation(default)
diff --git a/tests/07_OFDFT/28_OF_KE_XWM/INPUT b/tests/07_OFDFT/28_OF_KE_XWM/INPUT
new file mode 100644
index 0000000000..50408c9273
--- /dev/null
+++ b/tests/07_OFDFT/28_OF_KE_XWM/INPUT
@@ -0,0 +1,26 @@
+INPUT_PARAMETERS
+#Parameters (1.General)
+suffix		autotest
+calculation scf
+esolver_type    ofdft
+
+symmetry	1	
+pseudo_dir  ../../PP_ORB/
+pseudo_rcut 16
+nspin        1
+cal_force	1
+test_force	1
+
+#Parameters (2.Iteration)
+ecutwfc		20
+scf_nmax	50
+
+#OFDFT
+of_kinetic 	xwm
+of_method	tn
+of_conv 	energy
+of_tole 	2e-6
+
+#Parameters (3.Basis)
+basis_type	pw
+
diff --git a/tests/07_OFDFT/28_OF_KE_XWM/KPT b/tests/07_OFDFT/28_OF_KE_XWM/KPT
new file mode 100644
index 0000000000..c289c0158a
--- /dev/null
+++ b/tests/07_OFDFT/28_OF_KE_XWM/KPT
@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+1 1 1 0 0 0
diff --git a/tests/07_OFDFT/28_OF_KE_XWM/README b/tests/07_OFDFT/28_OF_KE_XWM/README
new file mode 100644
index 0000000000..12fc61c80d
--- /dev/null
+++ b/tests/07_OFDFT/28_OF_KE_XWM/README
@@ -0,0 +1 @@
+Test the energy and force of Xu-Wang-Ma (XWM) kinetic energy functional (of_method = xwm) in OFDFT, symmetry=on
diff --git a/tests/07_OFDFT/28_OF_KE_XWM/STRU b/tests/07_OFDFT/28_OF_KE_XWM/STRU
new file mode 100644
index 0000000000..e797c06432
--- /dev/null
+++ b/tests/07_OFDFT/28_OF_KE_XWM/STRU
@@ -0,0 +1,18 @@
+ATOMIC_SPECIES
+Al 26.98 al.lda.lps blps
+
+LATTICE_CONSTANT
+7.50241114482312  // add lattice constant
+
+LATTICE_VECTORS
+0.000000000000    0.500000000000    0.500000000000    
+0.500000000000    0.000000000000    0.500000000000    
+0.500000000000    0.500000000000    0.000000000000    
+
+ATOMIC_POSITIONS
+Direct
+
+Al
+0
+1
+    0.000000000000    0.000000000000    0.000000000000 1 1 1 
diff --git a/tests/07_OFDFT/28_OF_KE_XWM/result.ref b/tests/07_OFDFT/28_OF_KE_XWM/result.ref
new file mode 100644
index 0000000000..df43994366
--- /dev/null
+++ b/tests/07_OFDFT/28_OF_KE_XWM/result.ref
@@ -0,0 +1,7 @@
+etotref -57.9560215466872037
+etotperatomref -57.9560215467
+totalforceref 0.000000
+pointgroupref O_h
+spacegroupref O_h
+nksibzref 1
+totaltimeref 0.37
diff --git a/tests/07_OFDFT/29_OF_XWM_para/INPUT b/tests/07_OFDFT/29_OF_XWM_para/INPUT
new file mode 100644
index 0000000000..b8d23401f2
--- /dev/null
+++ b/tests/07_OFDFT/29_OF_XWM_para/INPUT
@@ -0,0 +1,28 @@
+INPUT_PARAMETERS
+#Parameters (1.General)
+suffix		autotest
+calculation scf
+esolver_type    ofdft
+
+symmetry	1	
+pseudo_dir  ../../PP_ORB/
+pseudo_rcut 16
+nspin        1
+cal_force	1
+test_force	1
+
+#Parameters (2.Iteration)
+ecutwfc		20
+scf_nmax	50
+
+#OFDFT
+of_kinetic 	xwm
+of_xwm_rho_ref  0.01
+of_xwm_kappa    0.1
+of_method	tn
+of_conv 	energy
+of_tole 	2e-6
+
+#Parameters (3.Basis)
+basis_type	pw
+
diff --git a/tests/07_OFDFT/29_OF_XWM_para/KPT b/tests/07_OFDFT/29_OF_XWM_para/KPT
new file mode 100644
index 0000000000..c289c0158a
--- /dev/null
+++ b/tests/07_OFDFT/29_OF_XWM_para/KPT
@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+1 1 1 0 0 0
diff --git a/tests/07_OFDFT/29_OF_XWM_para/README b/tests/07_OFDFT/29_OF_XWM_para/README
new file mode 100644
index 0000000000..ce9cbc1264
--- /dev/null
+++ b/tests/07_OFDFT/29_OF_XWM_para/README
@@ -0,0 +1 @@
+Test the energy and force of XWM, with of_xwm_rho_ref=0.01, of_xwm_kappa=0.1, symmetry=on
diff --git a/tests/07_OFDFT/29_OF_XWM_para/STRU b/tests/07_OFDFT/29_OF_XWM_para/STRU
new file mode 100644
index 0000000000..e797c06432
--- /dev/null
+++ b/tests/07_OFDFT/29_OF_XWM_para/STRU
@@ -0,0 +1,18 @@
+ATOMIC_SPECIES
+Al 26.98 al.lda.lps blps
+
+LATTICE_CONSTANT
+7.50241114482312  // add lattice constant
+
+LATTICE_VECTORS
+0.000000000000    0.500000000000    0.500000000000    
+0.500000000000    0.000000000000    0.500000000000    
+0.500000000000    0.500000000000    0.000000000000    
+
+ATOMIC_POSITIONS
+Direct
+
+Al
+0
+1
+    0.000000000000    0.000000000000    0.000000000000 1 1 1 
diff --git a/tests/07_OFDFT/29_OF_XWM_para/result.ref b/tests/07_OFDFT/29_OF_XWM_para/result.ref
new file mode 100644
index 0000000000..6ab2d9db1d
--- /dev/null
+++ b/tests/07_OFDFT/29_OF_XWM_para/result.ref
@@ -0,0 +1,7 @@
+etotref -57.9116843475440106
+etotperatomref -57.9116843475
+totalforceref 0.000000
+pointgroupref O_h
+spacegroupref O_h
+nksibzref 1
+totaltimeref 0.32
diff --git a/tests/07_OFDFT/CASES_CPU.txt b/tests/07_OFDFT/CASES_CPU.txt
index 7e7576b796..67e095d710 100644
--- a/tests/07_OFDFT/CASES_CPU.txt
+++ b/tests/07_OFDFT/CASES_CPU.txt
@@ -25,3 +25,5 @@
 25_OF_MD_LDA
 26_OF_MD_LibxcPBE
 27_OF_CR_LDA
+28_OF_KE_XWM
+29_OF_XWM_para
\ No newline at end of file

From bb5c65de1289fed64e3472ba5ebdcef478ad93da Mon Sep 17 00:00:00 2001
From: sunliang98 <1700011430@pku.edu.cn>
Date: Fri, 4 Jul 2025 17:17:24 +0800
Subject: [PATCH 3/4] Doc: Update document of XWM KEDF.

---
 docs/advanced/input_files/input-main.md | 21 +++++++++++++++++++--
 1 file changed, 19 insertions(+), 2 deletions(-)

diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
index b0df655cbf..d8dd5cd5ee 100644
--- a/docs/advanced/input_files/input-main.md
+++ b/docs/advanced/input_files/input-main.md
@@ -215,6 +215,8 @@
     - [of\_wt\_rho0](#of_wt_rho0)
     - [of\_hold\_rho0](#of_hold_rho0)
     - [of\_lkt\_a](#of_lkt_a)
+    - [of\_xwm\_rho\_ref](#of_xwm_rho_ref)
+    - [of\_xwm\_kappa](#of_xwm_kappa)
     - [of\_read\_kernel](#of_read_kernel)
     - [of\_kernel\_file](#of_kernel_file)
     - [of\_full\_pw](#of_full_pw)
@@ -2277,6 +2279,7 @@ Warning: this function is not robust enough for the current version. Please try
   - **vw**: von Weizsäcker.
   - **tf+**: TF $\rm{\lambda}$ vW, the parameter $\rm{\lambda}$ can be set by `of_vw_weight`.
   - **lkt**: Luo-Karasiev-Trickey.
+  - **xwm**: Xu-Wang-Ma
 
   Machine learning (ML) based functionals:
   - **ml**: ML-based KEDF allows for greater flexibility, enabling users to set related ML model parameters themselves. see [ML-KEDF: machine learning based kinetic energy density functional for OFDFT](#ml-kedf-machine-learning-based-kinetic-energy-density-functional-for-ofdft).
@@ -2323,14 +2326,14 @@ Warning: this function is not robust enough for the current version. Please try
 ### of_tf_weight
 
 - **Type**: Real
-- **Availability**: OFDFT with `of_kinetic=tf, tf+, wt`
+- **Availability**: OFDFT with `of_kinetic=tf, tf+, wt, xwm`
 - **Description**: Weight of TF KEDF (kinetic energy density functional).
 - **Default**: 1.0
 
 ### of_vw_weight
 
 - **Type**: Real
-- **Availability**: OFDFT with `of_kinetic=vw, tf+, wt, lkt`
+- **Availability**: OFDFT with `of_kinetic=vw, tf+, wt, lkt, xwm`
 - **Description**: Weight of vW KEDF (kinetic energy density functional).
 - **Default**: 1.0
 
@@ -2372,6 +2375,20 @@ Warning: this function is not robust enough for the current version. Please try
 - **Description**: Parameter a of LKT KEDF (kinetic energy density functional).
 - **Default**: 1.3
 
+### of_xwm_rho_ref
+
+- **Type**: Real
+- **Availability**: OFDFT with `of_kinetic=xwm`
+- **Description**: Reference charge density for XWM kinetic energy functional. If set to 0, the program will use average charge density.
+- **Default**: 0.0
+
+### of_xwm_kappa
+
+- **Type**: Real
+- **Availability**: OFDFT with `of_kinetic=xwm`
+- **Description**: Parameter $\kappa$ for XWM kinetic energy functional. See PHYSICAL REVIEW B 100, 205132 (2019) for optimal values.
+- **Default**: 0.0
+
 ### of_read_kernel
 
 - **Type**: Boolean

From 915165ac89bd0123d4e76006d82e4bdfd74d1dec Mon Sep 17 00:00:00 2001
From: sunliang98 <1700011430@pku.edu.cn>
Date: Fri, 4 Jul 2025 19:35:54 +0800
Subject: [PATCH 4/4] Fix: update kedf_manager.h

---
 source/source_pw/module_ofdft/kedf_manager.h | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/source/source_pw/module_ofdft/kedf_manager.h b/source/source_pw/module_ofdft/kedf_manager.h
index b9b311832f..fc1536b83c 100644
--- a/source/source_pw/module_ofdft/kedf_manager.h
+++ b/source/source_pw/module_ofdft/kedf_manager.h
@@ -4,12 +4,12 @@
 #include "module_parameter/parameter.h"
 #include "source_basis/module_pw/pw_basis.h"
 #include "source_estate/elecstate.h"
-#include "source_pw/hamilt_ofdft/kedf_lkt.h"
-#include "source_pw/hamilt_ofdft/kedf_tf.h"
-#include "source_pw/hamilt_ofdft/kedf_vw.h"
-#include "source_pw/hamilt_ofdft/kedf_wt.h"
-#include "source_pw/hamilt_ofdft/kedf_xwm.h"
-#include "source_pw/hamilt_ofdft/kedf_ml.h"
+#include "kedf_lkt.h"
+#include "kedf_tf.h"
+#include "kedf_vw.h"
+#include "kedf_wt.h"
+#include "kedf_xwm.h"
+#include "kedf_ml.h"
 
 class KEDF_Manager
 {