Feature: Implement Xu-Wang-Ma KEDF

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

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\
 

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

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";

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);

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)

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

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,14 +120,18 @@ 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);
     }
     if (this->of_kinetic_ == "wt")
     {
         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,19 +231,23 @@ 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]);
     }
     if (this->of_kinetic_ == "wt")
     {
         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");

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