diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
index 36700b8ff9..e61c1399d8 100644
--- a/docs/advanced/input_files/input-main.md
+++ b/docs/advanced/input_files/input-main.md
@@ -232,6 +232,7 @@
     - [of\_vw\_weight](#of_vw_weight)
     - [of\_wt\_alpha](#of_wt_alpha)
     - [of\_wt\_beta](#of_wt_beta)
+    - [of\_extwt\_kappa](#of_extwt_kappa)
     - [of\_wt\_rho0](#of_wt_rho0)
     - [of\_hold\_rho0](#of_hold_rho0)
     - [of\_lkt\_a](#of_lkt_a)
@@ -2395,7 +2396,8 @@
   - vw: von Weizsacker (vW) functional
   - tf+: TF + vW functional
   - wt: Wang-Teter (WT) functional
-  - xwm: XWM functional
+  - ext-wt: Extended Wang-Teter (ext-WT) functional
+  - xwm: Xu-Wang-Ma (XWM) functional
   - lkt: Luo-Karasiev-Trickey (LKT) functional
   - ml: Machine learning KEDF
   - mpn: MPN KEDF (automatically sets ml parameters)
@@ -2441,29 +2443,36 @@
 ### of_tf_weight
 
 - **Type**: Real
-- **Availability**: *OFDFT with of_kinetic=tf, tf+, wt, xwm*
+- **Availability**: *OFDFT with of_kinetic=tf, tf+, wt, ext-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, xwm*
+- **Availability**: *OFDFT with of_kinetic=vw, tf+, wt, ext-wt, lkt, xwm*
 - **Description**: Weight of vW KEDF (kinetic energy density functional).
 - **Default**: 1.0
 
 ### of_wt_alpha
 
 - **Type**: Real
-- **Availability**: *OFDFT with of_kinetic=wt*
+- **Availability**: *OFDFT with of_kinetic=wt, ext-wt*
 - **Description**: Parameter alpha of WT KEDF (kinetic energy density functional).
 
 ### of_wt_beta
 
 - **Type**: Real
-- **Availability**: *OFDFT with of_kinetic=wt*
+- **Availability**: *OFDFT with of_kinetic=wt, ext-wt*
 - **Description**: Parameter beta of WT KEDF (kinetic energy density functional).
 
+### of_extwt_kappa
+
+- **Type**: Real
+- **Availability**: *OFDFT with of_kinetic=ext-wt*
+- **Description**: Parameter kappa for EXT-WT KEDF.
+- **Default**: $\dfrac{1}{2(4/3)^{1/3}-1} \approx 0.832$
+
 ### of_wt_rho0
 
 - **Type**: Real
diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h
index 7ac2faf43e..cffcad7adc 100644
--- a/source/source_io/module_parameter/input_parameter.h
+++ b/source/source_io/module_parameter/input_parameter.h
@@ -196,6 +196,9 @@ struct Input_para
     double of_vw_weight = 1.0;                   ///< weight of vW KEDF
     double of_wt_alpha = 5. / 6.;                ///< parameter alpha of WT KEDF
     double of_wt_beta = 5. / 6.;                 ///< parameter beta of WT KEDF
+    double of_extwt_kappa = 1.0
+                            / (2.0 * std::pow(4. / 3., 1. / 3.) - 1.0);
+                                                 ///< parameter kappa of EXT-WT KEDF
     double of_wt_rho0 = 0.0;                     ///< set the average density of system, in Bohr^-3
     bool of_hold_rho0 = false;                   ///< If set to 1, the rho0 will be fixed even if the volume of
                                                  ///< system has changed, it will be set to 1 automatically if
diff --git a/source/source_io/module_parameter/read_input_item_ofdft.cpp b/source/source_io/module_parameter/read_input_item_ofdft.cpp
index 8254a1f1b0..31e3d2aeeb 100644
--- a/source/source_io/module_parameter/read_input_item_ofdft.cpp
+++ b/source/source_io/module_parameter/read_input_item_ofdft.cpp
@@ -20,6 +20,7 @@ void ReadInput::item_ofdft()
 * vw: von Weizsacker (vW) functional
 * tf+: TF + vW functional
 * wt: Wang-Teter (WT) functional
+* ext-wt: Extended Wang-Teter functional
 * xwm: XWM functional
 * lkt: Luo-Karasiev-Trickey (LKT) functional
 * ml: Machine learning KEDF
@@ -37,10 +38,11 @@ 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 != "ext-wt"
                 && 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, xwm, lkt, ml, mpn, or cpn5");
+                ModuleBase::WARNING_QUIT("ReadInput", "of_kinetic must be tf, vw, tf+, wt, ext-wt, xwm, lkt, ml, mpn, or cpn5");
             }
         };
         item.reset_value = [](const Input_Item& item, Parameter& para) {
@@ -175,7 +177,7 @@ void ReadInput::item_ofdft()
         item.description = "Weight of TF KEDF (kinetic energy density functional).";
         item.default_value = "1.0";
         item.unit = "";
-        item.availability = "OFDFT with of_kinetic=tf, tf+, wt, xwm";
+        item.availability = "OFDFT with of_kinetic=tf, tf+, wt, ext-wt, xwm";
         read_sync_double(input.of_tf_weight);
         this->add_item(item);
     }
@@ -187,7 +189,7 @@ void ReadInput::item_ofdft()
         item.description = "Weight of vW KEDF (kinetic energy density functional).";
         item.default_value = "1.0";
         item.unit = "";
-        item.availability = "OFDFT with of_kinetic=vw, tf+, wt, lkt, xwm";
+        item.availability = "OFDFT with of_kinetic=vw, tf+, wt, ext-wt, lkt, xwm";
         read_sync_double(input.of_vw_weight);
         this->add_item(item);
     }
@@ -199,7 +201,7 @@ void ReadInput::item_ofdft()
         item.description = "Parameter alpha of WT KEDF (kinetic energy density functional).";
         item.default_value = "";
         item.unit = "";
-        item.availability = "OFDFT with of_kinetic=wt";
+        item.availability = "OFDFT with of_kinetic=wt, ext-wt";
         read_sync_double(input.of_wt_alpha);
         this->add_item(item);
     }
@@ -211,10 +213,22 @@ void ReadInput::item_ofdft()
         item.description = "Parameter beta of WT KEDF (kinetic energy density functional).";
         item.default_value = "";
         item.unit = "";
-        item.availability = "OFDFT with of_kinetic=wt";
+        item.availability = "OFDFT with of_kinetic=wt, ext-wt";
         read_sync_double(input.of_wt_beta);
         this->add_item(item);
     }
+    {
+        Input_Item item("of_extwt_kappa");
+        item.annotation = "parameter kappa of EXT-WT KEDF";
+        item.category = "OFDFT: orbital free density functional theory";
+        item.type = "Real";
+        item.description = "Parameter kappa for EXT-WT KEDF.";
+        item.default_value = "1.0 / (2.0 * std::pow(4./3., 1./3.) - 1.0)";
+        item.unit = "";
+        item.availability = "OFDFT with of_kinetic=ext-wt";
+        read_sync_double(input.of_extwt_kappa);
+        this->add_item(item);
+    }
     {
         Input_Item item("of_wt_rho0");
         item.annotation = "the average density of system, used in WT KEDF, in Bohr^-3";
diff --git a/source/source_pw/module_ofdft/CMakeLists.txt b/source/source_pw/module_ofdft/CMakeLists.txt
index 8326faa0b4..404156457c 100644
--- a/source/source_pw/module_ofdft/CMakeLists.txt
+++ b/source/source_pw/module_ofdft/CMakeLists.txt
@@ -2,6 +2,7 @@ list(APPEND hamilt_ofdft_srcs
     kedf_tf.cpp
     kedf_vw.cpp
     kedf_wt.cpp
+    kedf_extwt.cpp
     kedf_xwm.cpp
     kedf_lkt.cpp
     kedf_manager.cpp
diff --git a/source/source_pw/module_ofdft/kedf_extwt.cpp b/source/source_pw/module_ofdft/kedf_extwt.cpp
new file mode 100644
index 0000000000..f336f38940
--- /dev/null
+++ b/source/source_pw/module_ofdft/kedf_extwt.cpp
@@ -0,0 +1,496 @@
+#include "./kedf_extwt.h"
+
+#include "source_io/module_parameter/parameter.h"
+#include <iostream>
+
+#include "source_base/global_variable.h"
+#include "source_base/parallel_reduce.h"
+#include "source_base/tool_quit.h"
+
+/**
+ * @brief Set the parameters of ext-WT KEDF, and initialize kernel
+ *
+ * @param dV the volume of one grid point in real space, omega/nxyz
+ * @param alpha
+ * @param beta
+ * @param nelec the number of electron
+ * @param tf_weight
+ * @param vw_weight
+ * @param of_extwt_kappa
+ * @param pw_rho pw_basis
+ */
+void KEDF_ExtWT::set_para(double dV,
+                       double alpha,
+                       double beta,
+                       double nelec,
+                       double tf_weight,
+                       double vw_weight,
+                       double of_extwt_kappa,
+                       ModulePW::PW_Basis* pw_rho)
+{
+    this->dV_ = dV;
+    // this->weightWT = weightWT;
+    this->alpha_ = alpha;
+    this->beta_ = beta;
+    this->kappa_ = of_extwt_kappa;
+    // std::cout << "kappa: " << this->kappa_ << std::endl;
+
+    this->rho0_ = 1. / (pw_rho->nxyz * dV) * nelec;
+
+    this->kf_ = std::pow(3. * std::pow(ModuleBase::PI, 2) * this->rho0_, 1. / 3.);
+    this->tkf_ = 2. * this->kf_;
+
+    this->wt_coef_
+        = 5. / (9. * this->alpha_ * this->beta_ * std::pow(this->rho0_, this->alpha_ + this->beta_ - 5. / 3.));
+
+    delete[] this->kernel_;
+    this->kernel_ = new double[pw_rho->npw];
+    this->fill_kernel(tf_weight, vw_weight, pw_rho);
+
+    delete[] this->dkernel_deta_;
+    this->dkernel_deta_ = new double[pw_rho->npw];
+
+    this->update_dkernel_deta(vw_weight, pw_rho);
+}
+
+void KEDF_ExtWT::update_rho0(const double* const* prho,
+                            ModulePW::PW_Basis* pw_rho)
+{
+    // Only for spin unpolarized case, need to be updated for spin polarized case
+    int nspin = 1;
+
+    this->sum_rho_kappa_ = 0.;
+    this->sum_rho_kappa_plus_one_ = 0.;
+    for (int is = 0; is < nspin; ++is)
+    {
+        for (int ir = 0; ir < pw_rho->nrxx; ++ir)
+        {
+            this->sum_rho_kappa_ += std::pow(prho[is][ir], this->kappa_);
+            this->sum_rho_kappa_plus_one_ += std::pow(prho[is][ir], this->kappa_ + 1.);
+        }
+    }
+    Parallel_Reduce::reduce_all(this->sum_rho_kappa_);
+    Parallel_Reduce::reduce_all(this->sum_rho_kappa_plus_one_);
+    this->sum_rho_kappa_ *= this->dV_;
+    this->sum_rho_kappa_plus_one_ *= this->dV_;
+
+    this->rho0_ = this->sum_rho_kappa_plus_one_ / this->sum_rho_kappa_;
+    // std::cout << "rho0: " << this->rho0_ << std::endl;
+}
+
+void KEDF_ExtWT::cal_kernel(double tf_weight,
+                double vw_weight,
+                double rho0,
+                ModulePW::PW_Basis* pw_rho)
+{
+    this->rho0_ = rho0;
+
+    this->kf_ = std::pow(3. * std::pow(ModuleBase::PI, 2) * this->rho0_, 1. / 3.);
+    this->tkf_ = 2. * this->kf_;
+
+    this->wt_coef_
+        = 5. / (9. * this->alpha_ * this->beta_ * std::pow(this->rho0_, this->alpha_ + this->beta_ - 5. / 3.));
+
+    this->fill_kernel(tf_weight, vw_weight, pw_rho);
+}
+
+/**
+ * @brief Fill the dkernel_deta (this->dkernel_deta_)
+ *
+ * @param alpha
+ * @param beta
+ * @param vw_weight
+ * @param pw_rho pw_basis
+ */
+void KEDF_ExtWT::update_dkernel_deta(const double &vw_weight, ModulePW::PW_Basis* pw_rho)
+{
+    double eta = 0.;
+    for (int ig = 0; ig < pw_rho->npw; ++ig)
+    {
+        eta = std::sqrt(pw_rho->gg[ig]) * pw_rho->tpiba / this->tkf_;
+        this->dkernel_deta_[ig] = - ((this->alpha_ + this->beta_ - 5./3.) * this->kernel_[ig] + 1./3. * eta * this->diff_linhard(eta, vw_weight))
+                                     * this->wt_coef_ / this->rho0_;
+    }
+}
+
+/**
+ * @brief Get the energy of ext-WT KEDF
+ *
+ * @param prho charge density
+ * @param pw_rho pw basis
+ * @return the energy of ext-WT KEDF
+ */
+double KEDF_ExtWT::get_energy(const double* const* prho, ModulePW::PW_Basis* pw_rho)
+{
+    double** kernelRhoBeta = new double*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is) {
+        kernelRhoBeta[is] = new double[pw_rho->nrxx];
+}
+    this->multi_kernel(prho, this->kernel_, kernelRhoBeta, this->beta_, 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->alpha_) * kernelRhoBeta[0][ir];
+        }
+        energy *= this->dV_ * this->c_tf_;
+    }
+    else if (PARAM.inp.nspin == 2)
+    {
+        // for (int is = 0; is < PARAM.inp.nspin; ++is)
+        // {
+        //     for (int ir = 0; ir < pw_rho->nrxx; ++ir)
+        //     {
+        //         energy += 2 * pphi[is][ir] * LapPhi[is][ir];
+        //     }
+        // }
+        // energy *= 0.5 * this->dV_ * 0.5;
+    }
+    this->extwt_energy = energy;
+    Parallel_Reduce::reduce_all(this->extwt_energy);
+
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        delete[] kernelRhoBeta[is];
+    }
+    delete[] kernelRhoBeta;
+
+    return this->extwt_energy;
+}
+
+/**
+ * @brief Get the energy density of ext-WT KEDF
+ *
+ * @param prho charge density
+ * @param is the index of spin
+ * @param ir the index of real space grid
+ * @param pw_rho pw basis
+ * @return the energy density of ext-WT KEDF
+ */
+double KEDF_ExtWT::get_energy_density(const double* const* prho, int is, int ir, ModulePW::PW_Basis* pw_rho)
+{
+    double** kernelRhoBeta = new double*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is) {
+        kernelRhoBeta[is] = new double[pw_rho->nrxx];
+}
+    this->multi_kernel(prho, this->kernel_, kernelRhoBeta, this->beta_, pw_rho);
+
+    double result = this->c_tf_ * std::pow(prho[is][ir], this->alpha_) * kernelRhoBeta[is][ir];
+
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        delete[] kernelRhoBeta[is];
+    }
+    delete[] kernelRhoBeta;
+    return result;
+}
+
+/**
+ * @brief Get the kinetic energy of ext-WT KEDF, and add it onto rtau_extwt
+ * 
+ * @param prho charge density
+ * @param pw_rho pw basis
+ * @param rtau_extwt rtau_extwt => rtau_extwt + tau_extwt
+ */
+void KEDF_ExtWT::tau_extwt(const double* const* prho, ModulePW::PW_Basis* pw_rho, double* rtau_extwt)
+{
+    double** kernelRhoBeta = new double*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is) {
+        kernelRhoBeta[is] = new double[pw_rho->nrxx];
+}
+    this->multi_kernel(prho, this->kernel_, kernelRhoBeta, this->beta_, pw_rho);
+
+    if (PARAM.inp.nspin == 1)
+    {
+        for (int ir = 0; ir < pw_rho->nrxx; ++ir)
+        {
+            rtau_extwt[ir] += std::pow(prho[0][ir], this->alpha_) * kernelRhoBeta[0][ir] * this->c_tf_;
+        }
+    }
+    else if (PARAM.inp.nspin == 2)
+    {
+        // Waiting for update
+    }
+
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        delete[] kernelRhoBeta[is];
+    }
+    delete[] kernelRhoBeta;
+}
+
+/**
+ * @brief Get the potential of ext-WT KEDF, and add it into rpotential,
+ * and the ext-WT energy will be calculated and stored in this->extwt_energy
+ * r)\rho^{\alpha}(r') dr'}] \f]
+ *
+ * @param prho charge density
+ * @param pw_rho pw basis
+ * @param rpotential rpotential => rpotential * 2 * phi + V_{extWT} * 2 * phi
+ */
+void KEDF_ExtWT::extwt_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpotential)
+{
+    ModuleBase::timer::start("KEDF_ExtWT", "extwt_potential");
+
+    // 1. WT potential
+    double** kernelRhoBeta = new double*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is) {
+        kernelRhoBeta[is] = new double[pw_rho->nrxx];
+    }
+    this->multi_kernel(prho, this->kernel_, kernelRhoBeta, this->beta_, pw_rho);
+
+    double** kernelRhoAlpha = new double*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is) {
+        kernelRhoAlpha[is] = new double[pw_rho->nrxx];
+    }
+    this->multi_kernel(prho, this->kernel_, kernelRhoAlpha, this->alpha_, pw_rho);
+
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        for (int ir = 0; ir < pw_rho->nrxx; ++ir)
+        {
+            rpotential(is, ir) += this->c_tf_
+                                  * (this->alpha_ * std::pow(prho[is][ir], this->alpha_ - 1.) * kernelRhoBeta[is][ir]
+                                     + this->beta_ * std::pow(prho[is][ir], this->beta_ - 1.) * kernelRhoAlpha[is][ir])
+                                     * 2. * std::sqrt(prho[is][ir]);
+        }
+    }
+
+    // 2. rho0 part
+    // 1) prepare variables
+    std::vector<std::vector<double>> rho_alpha(PARAM.inp.nspin, std::vector<double>(pw_rho->nrxx, 0.));
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        for (int ir = 0; ir < pw_rho->nrxx; ++ir)
+        {
+            rho_alpha[is][ir] = std::pow(prho[is][ir], this->alpha_);
+        }
+    }
+
+    // 2) calculate the constants
+    double coef = 0.;
+
+    double** dkernelRhoBeta = new double*[PARAM.inp.nspin];
+    for (int is = 0; is < PARAM.inp.nspin; ++is) {
+        dkernelRhoBeta[is] = new double[pw_rho->nrxx];
+    }
+    this->multi_kernel(prho, this->dkernel_deta_, dkernelRhoBeta, this->beta_, pw_rho);
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        for (int ir = 0; ir < pw_rho->nrxx; ++ir)
+        {
+            coef += rho_alpha[is][ir] * dkernelRhoBeta[is][ir];
+        }
+    }
+    Parallel_Reduce::reduce_all(coef);
+    coef *= this->dV_ * this->c_tf_;
+
+    // 3) calculate the total potential
+    int count = 0;
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        for (int ir = 0; ir < pw_rho->nrxx; ++ir)
+        {
+            double second_term = std::pow(prho[is][ir], this->kappa_ - 1. + 0.5); // 0.5 is for sqrt(rho)
+            if (second_term > 100 && count == 0)
+            {
+                count++;
+                std::cout << "Warning: second_term is too large: " << second_term << ", coef = " << coef << std::endl;
+            }
+            if (second_term >= 10000)
+            {
+                second_term = 10000;
+            }
+            rpotential(is, ir) += coef / this->sum_rho_kappa_ * ((this->kappa_ + 1) * std::pow(prho[is][ir], this->kappa_ + 0.5) - this->kappa_ * second_term * this->rho0_) * 2.; // 2 is for 2 * phi
+        }
+    }
+
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        delete[] dkernelRhoBeta[is];
+    }
+    delete[] dkernelRhoBeta;
+
+    // calculate energy
+    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->alpha_) * kernelRhoBeta[0][ir];
+        }
+        energy *= this->dV_ * this->c_tf_;
+    }
+    else if (PARAM.inp.nspin == 2)
+    {
+        // for (int is = 0; is < PARAM.inp.nspin; ++is)
+        // {
+        //     for (int ir = 0; ir < pw_rho->nrxx; ++ir)
+        //     {
+        //         energy += 2 * pphi[is][ir] * LapPhi[is][ir];
+        //     }
+        // }
+        // energy *= 0.5 * this->dV_ * 0.5;
+    }
+    this->extwt_energy = energy;
+    Parallel_Reduce::reduce_all(this->extwt_energy);
+
+    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    {
+        delete[] kernelRhoBeta[is];
+        delete[] kernelRhoAlpha[is];
+    }
+    delete[] kernelRhoBeta;
+    delete[] kernelRhoAlpha;
+    ModuleBase::timer::end("KEDF_ExtWT", "extwt_potential");
+}
+
+/**
+ * @brief Get the stress of ext-WT 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_ExtWT::get_stress(const double* const* prho, ModulePW::PW_Basis* pw_rho, double vw_weight)
+{
+    ModuleBase::WARNING_QUIT("KEDF_ExtWT", "ext-WT stress is not implemented yet!");
+}
+
+/**
+ * @brief Calculate the WT kernel according to Lindhard response function
+ *
+ * @param eta k / (2 * kF)
+ * @param tf_weight
+ * @param vw_weight
+ * @return W(eta)
+ */
+double KEDF_ExtWT::extwt_kernel(double eta, double tf_weight, double vw_weight)
+{
+    if (eta < 0.)
+    {
+        return 0.;
+    }
+    // limit for small eta
+    else if (eta < 1e-10)
+    {
+        return 1. - tf_weight + eta * eta * (1. / 3. - 3. * vw_weight);
+    }
+    // around the singularity
+    else if (std::abs(eta - 1.) < 1e-10)
+    {
+        return 2. - tf_weight - 3. * vw_weight + 20. * (eta - 1);
+    }
+    // Taylor expansion for high eta
+    else if (eta > 3.65)
+    {
+        double eta2 = eta * eta;
+        double invEta2 = 1. / eta2;
+        double LindG = 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)))))))))))))));
+        return LindG;
+    }
+    else
+    {
+        return 1. / (0.5 + 0.25 * (1. - eta * eta) / eta * std::log((1 + eta) / std::abs(1 - eta)))
+               - 3. * vw_weight * eta * eta - tf_weight;
+    }
+}
+
+/**
+ * @brief The derivative of the WT kernel
+ *
+ * @param eta k / (2 * kF)
+ * @param vw_weight
+ * @return d W(eta)/d eta
+ */
+double KEDF_ExtWT::diff_linhard(double eta, double vw_weight)
+{
+    if (eta < 0.)
+    {
+        return 0.;
+    }
+    else if (eta < 1e-10)
+    {
+        return 2. * eta * (1. / 3. - 3. * vw_weight);
+    }
+    else if (std::abs(eta - 1.) < 1e-10)
+    {
+        return 40.;
+    }
+    else
+    {
+        double eta2 = eta * eta;
+        return ((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;
+    }
+}
+
+/**
+ * @brief Calculate \int{W(r-r')rho^{exponent}(r') dr'}
+ *
+ * @param [in] prho charge density
+ * @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_ExtWT::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];
+            // recipkernelRho[is][ip] *= this->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_ExtWT::fill_kernel(double tf_weight, double vw_weight, ModulePW::PW_Basis* pw_rho)
+{
+    double eta = 0.;
+    for (int ig = 0; ig < pw_rho->npw; ++ig)
+    {
+        eta = std::sqrt(pw_rho->gg[ig]) * pw_rho->tpiba / this->tkf_;
+        this->kernel_[ig] = this->extwt_kernel(eta, tf_weight, vw_weight) * this->wt_coef_;
+    }
+}
\ No newline at end of file
diff --git a/source/source_pw/module_ofdft/kedf_extwt.h b/source/source_pw/module_ofdft/kedf_extwt.h
new file mode 100644
index 0000000000..539ee77f63
--- /dev/null
+++ b/source/source_pw/module_ofdft/kedf_extwt.h
@@ -0,0 +1,80 @@
+#ifndef KEDF_EXTWT_H
+#define KEDF_EXTWT_H
+#include <cmath>
+#include <cstdio>
+
+#include "source_base/global_function.h"
+#include "source_base/global_variable.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 extended Wang-Teter (ext-WT) KEDF.
+ * See Sun L, Chen M. arXiv preprint arXiv:2507.08442, 2025.
+ * @author sunliang on 2026-04
+ */
+class KEDF_ExtWT
+{
+  public:
+    KEDF_ExtWT()
+    {
+        this->stress.create(3, 3);
+    }
+    ~KEDF_ExtWT()
+    {
+        delete[] this->kernel_;
+        delete[] this->dkernel_deta_;
+    }
+
+    void set_para(double dV,
+                  double alpha,
+                  double beta,
+                  double nelec,
+                  double tf_weight,
+                  double vw_weight,
+                  double of_extwt_kappa,
+                  ModulePW::PW_Basis* pw_rho);
+
+    void update_rho0(const double* const* prho,
+                     ModulePW::PW_Basis* pw_rho);
+    void cal_kernel(double tf_weight,
+                    double vw_weight,
+                    double rho0,
+                    ModulePW::PW_Basis* pw_rho);
+
+    void update_dkernel_deta(const 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_extwt(const double* const* prho, ModulePW::PW_Basis* pw_rho, double* rtau_extwt);
+    void extwt_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 extwt_energy = 0.;
+    ModuleBase::matrix stress;
+
+    double rho0_ = 0.; // average rho
+
+  private:
+    double extwt_kernel(double eta, double tf_weight, double vw_weight);
+    double diff_linhard(double eta, double vw_weight);
+    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 kf_ = 0.;  // Fermi vector kF = (3 pi^2 rho)^(1/3)
+    double tkf_ = 0.; // 2 * kF
+    double alpha_ = 5. / 6.;
+    double beta_ = 5. / 6.;
+    // double weightWT = 1.;
+    const double c_tf_
+        = 3.0 / 10.0 * std::pow(3 * std::pow(M_PI, 2.0), 2.0 / 3.0)
+          * 2; // 10/3*(3*pi^2)^{2/3}, multiply by 2 to convert unit from Hartree to Ry, finally in Ry*Bohr^(-2)
+    double wt_coef_ = 0.; // coefficient of WT kernel
+    double* kernel_ = nullptr;
+    double* dkernel_deta_ = nullptr; // \partial w/ \partial rho0 = coef * ((alpha + beta - 5/3)F(eta) + 1/3 eta F'(eta))
+    double sum_rho_kappa_ = 0.;
+    double sum_rho_kappa_plus_one_ = 0.;
+    double kappa_ = 1.0 / (2.0 * std::pow(4./3., 1./3.) - 1.0);
+};
+#endif
\ No newline at end of file
diff --git a/source/source_pw/module_ofdft/kedf_manager.cpp b/source/source_pw/module_ofdft/kedf_manager.cpp
index 0916ebc25a..9c589fa9a9 100644
--- a/source/source_pw/module_ofdft/kedf_manager.cpp
+++ b/source/source_pw/module_ofdft/kedf_manager.cpp
@@ -22,6 +22,7 @@ void KEDF_Manager::init(
     if (this->of_kinetic_ == "tf"
      || this->of_kinetic_ == "tf+"
      || this->of_kinetic_ == "wt"
+        || this->of_kinetic_ == "ext-wt"
      || this->of_kinetic_ == "ml"
      || this->of_kinetic_ == "xwm")
     {
@@ -34,6 +35,7 @@ void KEDF_Manager::init(
 
     //! vW, TF+, WT, XWM, and LKT KEDFs
     if (this->of_kinetic_ == "vw" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt"
+        || this->of_kinetic_ == "ext-wt"
         || this->of_kinetic_ == "xwm" || this->of_kinetic_ == "lkt" || this->of_kinetic_ == "ml")
     {
         if (this->vw_ == nullptr)
@@ -63,6 +65,23 @@ void KEDF_Manager::init(
                             pw_rho);
     }
 
+    //! Extended Wang-Teter KEDF
+    if (this->of_kinetic_ == "ext-wt")
+    {
+        if (this->extwt_ == nullptr)
+        {
+            this->extwt_ = new KEDF_ExtWT();
+        }
+        this->extwt_->set_para(dV,
+                               inp.of_wt_alpha,
+                               inp.of_wt_beta,
+                               nelec,
+                               inp.of_tf_weight,
+                               inp.of_vw_weight,
+                               inp.of_extwt_kappa,
+                               pw_rho);
+    }
+
     //! Xu-Wang-Ma KEDF
     if (this->of_kinetic_ == "xwm")
     {
@@ -120,7 +139,8 @@ 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" || this->of_kinetic_ == "xwm")
+    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt"
+        || this->of_kinetic_ == "ext-wt" || this->of_kinetic_ == "xwm")
     {
         this->tf_->tf_potential(prho, rpot);
     }
@@ -154,10 +174,19 @@ void KEDF_Manager::get_potential(
     }
 
     if (this->of_kinetic_ == "vw" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt"
+        || this->of_kinetic_ == "ext-wt"
         || this->of_kinetic_ == "xwm" || this->of_kinetic_ == "lkt" || this->of_kinetic_ == "ml")
     {
         this->vw_->vw_potential(pphi, pw_rho, rpot);
     }
+
+    if (this->of_kinetic_ == "ext-wt")
+    {
+        this->extwt_->update_rho0(prho, pw_rho);
+        this->extwt_->cal_kernel(PARAM.inp.of_tf_weight, PARAM.inp.of_vw_weight, this->extwt_->rho0_, pw_rho);
+        this->extwt_->update_dkernel_deta(PARAM.inp.of_vw_weight, pw_rho);
+        this->extwt_->extwt_potential(prho, pw_rho, rpot);
+    }
 }
 
 /**
@@ -170,12 +199,14 @@ double KEDF_Manager::get_energy() const
 {
     double kinetic_energy = 0.0;
 
-    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt"|| this->of_kinetic_ == "xwm")
+    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt"
+        || this->of_kinetic_ == "ext-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_ == "ext-wt"
         || this->of_kinetic_ == "xwm" || this->of_kinetic_ == "lkt" || this->of_kinetic_ == "ml")
     {
         kinetic_energy += this->vw_->vw_energy;
@@ -186,6 +217,11 @@ double KEDF_Manager::get_energy() const
         kinetic_energy += this->wt_->wt_energy;
     }
 
+    if (this->of_kinetic_ == "ext-wt")
+    {
+        kinetic_energy += this->extwt_->extwt_energy;
+    }
+
     if (this->of_kinetic_ == "xwm")
     {
         kinetic_energy += this->xwm_->xwm_energy;
@@ -231,11 +267,13 @@ void KEDF_Manager::get_energy_density(
         rtau[0][ir] = 0.0;
     }
 
-    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt" || this->of_kinetic_ == "xwm")
+    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt"
+        || this->of_kinetic_ == "ext-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_ == "ext-wt"
         || this->of_kinetic_ == "xwm" || this->of_kinetic_ == "lkt")
     {
         this->vw_->tau_vw(pphi, pw_rho, rtau[0]);
@@ -244,6 +282,10 @@ void KEDF_Manager::get_energy_density(
     {
         this->wt_->tau_wt(prho, pw_rho, rtau[0]);
     }
+    if (this->of_kinetic_ == "ext-wt")
+    {
+        this->extwt_->tau_extwt(prho, pw_rho, rtau[0]);
+    }
     if (this->of_kinetic_ == "xwm")
     {
         this->xwm_->tau_xwm(prho, pw_rho, rtau[0]);
@@ -280,13 +322,15 @@ void KEDF_Manager::get_stress(
         }
     }
 
-    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt" || this->of_kinetic_ == "xwm")
+    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt"
+        || this->of_kinetic_ == "ext-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_ == "ext-wt"
         || this->of_kinetic_ == "xwm" || this->of_kinetic_ == "lkt")
     {
         this->vw_->get_stress(pphi, pw_rho);
@@ -299,6 +343,12 @@ void KEDF_Manager::get_stress(
         kinetic_stress_ += this->wt_->stress;
     }
 
+    if (this->of_kinetic_ == "ext-wt")
+    {
+        this->extwt_->get_stress(prho, pw_rho, PARAM.inp.of_vw_weight);
+        kinetic_stress_ += this->extwt_->stress;
+    }
+
     if (this->of_kinetic_ == "xwm")
     {
         this->xwm_->get_stress(prho, pw_rho, PARAM.inp.of_vw_weight);
@@ -321,12 +371,14 @@ void KEDF_Manager::record_energy(
     std::vector<double> &energies_Ry
 )
 {
-    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt" || this->of_kinetic_ == "xwm")
+    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt"
+        || this->of_kinetic_ == "ext-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_ == "ext-wt"
         || this->of_kinetic_ == "xwm" || this->of_kinetic_ == "lkt" || this->of_kinetic_ == "ml")
     {
         titles.push_back("vW KEDF");
@@ -337,6 +389,11 @@ void KEDF_Manager::record_energy(
         titles.push_back("WT KEDF");
         energies_Ry.push_back(this->wt_->wt_energy);
     }
+    if (this->of_kinetic_ == "ext-wt")
+    {
+        titles.push_back("EXT-WT KEDF");
+        energies_Ry.push_back(this->extwt_->extwt_energy);
+    }
     if (this->of_kinetic_ == "xwm")
     {
         titles.push_back("XWM KEDF");
diff --git a/source/source_pw/module_ofdft/kedf_manager.h b/source/source_pw/module_ofdft/kedf_manager.h
index af80646d6a..8cd556d786 100644
--- a/source/source_pw/module_ofdft/kedf_manager.h
+++ b/source/source_pw/module_ofdft/kedf_manager.h
@@ -8,6 +8,7 @@
 #include "kedf_tf.h"
 #include "kedf_vw.h"
 #include "kedf_wt.h"
+#include "kedf_extwt.h"
 #include "kedf_xwm.h"
 #include "kedf_ml.h"
 
@@ -22,6 +23,8 @@ class KEDF_Manager
         delete this->tf_;
         delete this->vw_;
         delete this->wt_;
+        delete this->extwt_;
+        delete this->xwm_;
 #ifdef __MLALGO
         delete this->ml_;
 #endif
@@ -75,6 +78,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_ExtWT* extwt_ = nullptr; // Extended Wang-Teter KEDF
     KEDF_XWM* xwm_ = nullptr; // Xu-Wang-Ma KEDF
 #ifdef __MLALGO
     KEDF_ML* ml_ = nullptr;   // Machine Learning KEDF
diff --git a/source/source_pw/module_ofdft/kedf_xwm.cpp b/source/source_pw/module_ofdft/kedf_xwm.cpp
index f80e5044a2..cacd3cdf23 100644
--- a/source/source_pw/module_ofdft/kedf_xwm.cpp
+++ b/source/source_pw/module_ofdft/kedf_xwm.cpp
@@ -265,7 +265,7 @@ void KEDF_XWM::xwm_potential(const double* const* prho, ModulePW::PW_Basis* pw_r
  */
 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;
+    ModuleBase::WARNING_QUIT("KEDF_XWM", "XWM stress is not implemented yet!");
 }
 
 /**
diff --git a/tests/07_OFDFT/34_OF_KE_extWT/INPUT b/tests/07_OFDFT/34_OF_KE_extWT/INPUT
new file mode 100644
index 0000000000..8ace44f827
--- /dev/null
+++ b/tests/07_OFDFT/34_OF_KE_extWT/INPUT
@@ -0,0 +1,23 @@
+INPUT_PARAMETERS 
+#Parameters (1.General)
+suffix            autotest
+calculation       scf
+esolver_type      ofdft
+
+symmetry          1
+pseudo_dir        ../../PP_ORB/
+pseudo_rcut       16
+nspin             1
+
+#Parameters (2.Iteration)
+ecutwfc           20
+scf_nmax          50
+
+#OFDFT
+of_kinetic        ext-wt
+of_method         tn
+of_conv           energy
+of_tole           2e-6
+
+#Parameters (3.Basis)
+basis_type        pw
diff --git a/tests/07_OFDFT/34_OF_KE_extWT/KPT b/tests/07_OFDFT/34_OF_KE_extWT/KPT
new file mode 100644
index 0000000000..c289c0158a
--- /dev/null
+++ b/tests/07_OFDFT/34_OF_KE_extWT/KPT
@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+1 1 1 0 0 0
diff --git a/tests/07_OFDFT/34_OF_KE_extWT/README b/tests/07_OFDFT/34_OF_KE_extWT/README
new file mode 100644
index 0000000000..fe2305ad46
--- /dev/null
+++ b/tests/07_OFDFT/34_OF_KE_extWT/README
@@ -0,0 +1 @@
+Test the energy of extended Wang-Teter (ext-WT) kinetic energy functional (of_method = ext-wt) in OFDFT, symmetry=on
diff --git a/tests/07_OFDFT/34_OF_KE_extWT/STRU b/tests/07_OFDFT/34_OF_KE_extWT/STRU
new file mode 100644
index 0000000000..e797c06432
--- /dev/null
+++ b/tests/07_OFDFT/34_OF_KE_extWT/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/34_OF_KE_extWT/result.ref b/tests/07_OFDFT/34_OF_KE_extWT/result.ref
new file mode 100644
index 0000000000..ecf4b063ca
--- /dev/null
+++ b/tests/07_OFDFT/34_OF_KE_extWT/result.ref
@@ -0,0 +1,6 @@
+etotref -57.93042539254519
+etotperatomref -57.9304253925
+pointgroupref O_h
+spacegroupref O_h
+nksibzref 1
+totaltimeref 0.39
diff --git a/tests/07_OFDFT/CASES_CPU.txt b/tests/07_OFDFT/CASES_CPU.txt
index a6e6110ca3..d58f16da56 100644
--- a/tests/07_OFDFT/CASES_CPU.txt
+++ b/tests/07_OFDFT/CASES_CPU.txt
@@ -31,3 +31,4 @@
 31_TDOFDFT_Al_CD
 32_TDOFDFT_Al_mCD
 33_OF_out_chg
+34_OF_KE_extWT
\ No newline at end of file