From 47d940f13cdc57f0e5ccbeae9e25441a96df9661 Mon Sep 17 00:00:00 2001
From: Growl <uwsy1059@qq.com>
Date: Sun, 24 May 2026 20:03:01 +0800
Subject: [PATCH 1/7] Add optional DFT-D4 support

---
 CMakeLists.txt                                |  13 +
 source/source_estate/elecstate_print.cpp      |   5 +
 .../source_hamilt/module_vdw/CMakeLists.txt   |  33 +-
 source/source_hamilt/module_vdw/vdw.cpp       |  26 ++
 source/source_hamilt/module_vdw/vdwd4.cpp     | 286 ++++++++++++++++++
 source/source_hamilt/module_vdw/vdwd4.h       |  50 +++
 .../module_parameter/input_parameter.h        |   3 +-
 .../read_input_item_model.cpp                 |  32 +-
 source/source_lcao/FORCE_STRESS.cpp           |   2 +-
 .../source_pw/module_ofdft/of_print_info.cpp  |   5 +
 toolchain/build_abacus_gcc-mkl.sh             |   1 +
 toolchain/build_abacus_gnu.sh                 |   1 +
 toolchain/build_abacus_intel.sh               |   1 +
 toolchain/scripts/lib/config_manager.sh       |   3 +-
 toolchain/scripts/lib/package_manager.sh      |   3 +-
 toolchain/scripts/lib/version_helper.sh       |   2 +-
 toolchain/scripts/package_versions.sh         |   8 +
 toolchain/scripts/stage4/install_dftd4.sh     | 107 +++++++
 toolchain/scripts/stage4/install_stage4.sh    |   1 +
 toolchain/toolchain_gcc-mkl.sh                |   4 +
 toolchain/toolchain_gnu.sh                    |   4 +
 toolchain/toolchain_intel.sh                  |   4 +
 22 files changed, 575 insertions(+), 19 deletions(-)
 create mode 100644 source/source_hamilt/module_vdw/vdwd4.cpp
 create mode 100644 source/source_hamilt/module_vdw/vdwd4.h
 create mode 100755 toolchain/scripts/stage4/install_dftd4.sh

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5c6ece90f9c..9df21ca5c9e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -34,6 +34,7 @@ option(USE_ELPA "Enable ELPA for LCAO" ON)
 option(ENABLE_LIBRI "Enable LibRI for hybrid functional" OFF)
 option(ENABLE_LIBCOMM "Enable LibComm" OFF)
 option(ENABLE_PEXSI "Enable PEXSI for LCAO" OFF)
+option(ENABLE_DFTD4 "Enable DFT-D4 dispersion correction" OFF)
 
 option(BUILD_TESTING "Build unittests" OFF)
 option(DEBUG_INFO "Print message to debug" OFF)
@@ -236,6 +237,18 @@ if(ENABLE_COVERAGE)
   add_coverage(${ABACUS_BIN_NAME})
 endif()
 
+if(ENABLE_DFTD4)
+  # DFTD4 requires enabling C and Fortran to work
+  enable_language(C)
+  enable_language(Fortran)
+  # Avoid custom-lapack fallback when resolving DFT-D4 dependencies
+  find_package(BLAS REQUIRED)
+  find_package(LAPACK REQUIRED)
+  find_package(dftd4 4.2.0 REQUIRED)
+  if(NOT TARGET dftd4::dftd4)
+    message(FATAL_ERROR "DFT-D4 was found, but target dftd4::dftd4 is missing.")
+  endif()
+endif()
 
 macro(set_if_higher VARIABLE VALUE)
   if(${VARIABLE} LESS ${VALUE})
diff --git a/source/source_estate/elecstate_print.cpp b/source/source_estate/elecstate_print.cpp
index 517c002dd6d..5c61755d6ff 100644
--- a/source/source_estate/elecstate_print.cpp
+++ b/source/source_estate/elecstate_print.cpp
@@ -250,6 +250,11 @@ void print_etot(const Magnetism& magnet,
             titles.push_back("E_vdwD3");
             energies_Ry.push_back(elec.f_en.evdw);
         }
+        else if (vdw_method == "d4")
+        {
+            titles.push_back("E_vdwD4");
+            energies_Ry.push_back(elec.f_en.evdw);
+        }
 
         // mohan add 20251108
 		if (PARAM.inp.dft_plus_u)
diff --git a/source/source_hamilt/module_vdw/CMakeLists.txt b/source/source_hamilt/module_vdw/CMakeLists.txt
index 6f2c1800a9f..f59434bde71 100644
--- a/source/source_hamilt/module_vdw/CMakeLists.txt
+++ b/source/source_hamilt/module_vdw/CMakeLists.txt
@@ -1,16 +1,25 @@
-add_library(
-        vdw
-        OBJECT
-        vdwd2_parameters.cpp
-        vdwd3_parameters_tab.cpp
-        vdwd3_parameters.cpp
-        vdwd2.cpp
-        vdwd3.cpp
-        vdwd3_autoset_xcname.cpp
-        vdwd3_autoset_xcparam.cpp
-        vdw.cpp
+set(vdw_sources
+    vdwd2_parameters.cpp
+    vdwd3_parameters_tab.cpp
+    vdwd3_parameters.cpp
+    vdwd2.cpp
+    vdwd3.cpp
+    vdwd3_autoset_xcname.cpp
+    vdwd3_autoset_xcparam.cpp
+    vdw.cpp
 )
 
+if(ENABLE_DFTD4)
+  list(APPEND vdw_sources vdwd4.cpp)
+endif()
+
+add_library(vdw OBJECT ${vdw_sources})
+
+if(ENABLE_DFTD4)
+  target_link_libraries(vdw PUBLIC dftd4::dftd4)
+  target_compile_definitions(vdw PUBLIC __DFTD4)
+endif()
+
 if(ENABLE_COVERAGE)
   add_coverage(vdw)
 endif()
@@ -19,4 +28,4 @@ if(BUILD_TESTING)
   if(ENABLE_MPI)
     add_subdirectory(test)
   endif()
-endif()
\ No newline at end of file
+endif()
diff --git a/source/source_hamilt/module_vdw/vdw.cpp b/source/source_hamilt/module_vdw/vdw.cpp
index b9c61858651..4d19d8a0bf3 100644
--- a/source/source_hamilt/module_vdw/vdw.cpp
+++ b/source/source_hamilt/module_vdw/vdw.cpp
@@ -6,6 +6,10 @@
 #include "vdwd3.h"
 #include "source_base/tool_quit.h"
 
+#ifdef __DFTD4
+#include "vdwd4.h"
+#endif
+
 std::string parse_xcname(const std::string &xc_input,
                          const std::vector<std::string> &xc_psp)
 {
@@ -70,6 +74,28 @@ std::unique_ptr<Vdw> make_vdw(const UnitCell &ucell,
         vdw_ptr->parameter().initial_parameters(parse_xcname(input.dft_functional, xc_psp), input, plog);
         return vdw_ptr;
     }
+    else if (input.vdw_method == "d4")
+    {
+#ifdef __DFTD4
+        std::vector<std::string> xc_psp(ucell.ntype);
+        for (int it = 0; it < ucell.ntype; it++)
+        {
+            xc_psp[it] = ucell.atoms[it].ncpp.xc_func;
+        }
+
+        std::string xc_name = input.vdw_d4_xc;
+        if (xc_name == "default")
+        {
+            xc_name = parse_xcname(input.dft_functional, xc_psp);
+        }
+
+        return vdw::make_unique<Vdwd4>(ucell, xc_name, input);
+#else
+        ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::make_vdw",
+                                 "DFT-D4 support was not enabled at build time. "
+                                 "Rebuild ABACUS with -DENABLE_DFTD4=ON.");
+#endif
+    }
     else if (input.vdw_method != "none")
     {
         ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::make_vdw", 
diff --git a/source/source_hamilt/module_vdw/vdwd4.cpp b/source/source_hamilt/module_vdw/vdwd4.cpp
new file mode 100644
index 00000000000..7e3516c2026
--- /dev/null
+++ b/source/source_hamilt/module_vdw/vdwd4.cpp
@@ -0,0 +1,286 @@
+#include "vdwd4.h"
+
+#include "source_base/constants.h"
+#include "source_base/element_name.h"
+#include "source_base/timer.h"
+#include "source_base/tool_quit.h"
+
+#include <algorithm>
+#include <array>
+#include <cctype>
+#include <cstddef>
+#include <cstring>
+#include <iostream>
+#include <string>
+#include <vector>
+
+#ifdef __DFTD4
+#include <dftd4.h>
+#endif
+
+namespace vdw
+{
+
+namespace
+{
+
+std::string normalize_xc_name(std::string xc)
+{
+    std::transform(xc.begin(), xc.end(), xc.begin(), [](unsigned char c) {
+        return static_cast<char>(std::tolower(c));
+    });
+    return xc;
+}
+
+void check_dftd4_error(dftd4_error error, const std::string& where)
+{
+    if (dftd4_check_error(error))
+    {
+        char buffer[1024];
+        std::memset(buffer, 0, sizeof(buffer));
+        int buffersize = static_cast<int>(sizeof(buffer));
+        dftd4_get_error(error, buffer, &buffersize);
+
+        ModuleBase::WARNING_QUIT("Vdwd4::" + where, std::string(buffer));
+    }
+}
+
+int atomic_number_from_symbol(const std::string& symbol)
+{
+    for (int i = 0; i < static_cast<int>(ModuleBase::element_name.size()); ++i)
+    {
+        if (symbol == ModuleBase::element_name[i])
+        {
+            return i + 1; // DFT-D4 expects the true atomic number.
+        }
+    }
+
+    ModuleBase::WARNING_QUIT("Vdwd4::atomic_number_from_symbol", "Unknown element symbol: " + symbol);
+    return 0;
+}
+
+double length_to_bohr(double value, const std::string& unit)
+{
+    if (unit == "A")
+    {
+        return value / ModuleBase::BOHR_TO_A;
+    }
+    if (unit == "Bohr")
+    {
+        return value;
+    }
+
+    ModuleBase::WARNING_QUIT("Vdwd4::length_to_bohr", "Unsupported length unit: " + unit);
+    return value;
+}
+
+double cutoff_to_bohr(const std::string& value, const std::string& unit)
+{
+    return length_to_bohr(std::stod(value), unit);
+}
+
+} // namespace
+
+Vdwd4::Vdwd4(const UnitCell& unit_in, const std::string& xc_name, const Input_para& input)
+    : Vdw(unit_in), xc_name_(normalize_xc_name(xc_name))
+{
+    cutoff_disp2_ = cutoff_to_bohr(input.vdw_cutoff_radius, input.vdw_radius_unit);
+    cutoff_disp3_ = std::min(40.0, cutoff_disp2_);
+    cutoff_cn_ = length_to_bohr(input.vdw_cn_thr, input.vdw_cn_thr_unit);
+}
+
+void Vdwd4::build_structure(std::vector<int>& numbers,
+                            std::vector<double>& positions,
+                            std::vector<double>& lattice,
+                            std::array<bool, 3>& periodic) const
+{
+    numbers.clear();
+    positions.clear();
+    lattice.clear();
+
+    numbers.reserve(ucell_.nat);
+    positions.reserve(3 * ucell_.nat);
+    lattice.reserve(9);
+
+    for (int it = 0; it < ucell_.ntype; ++it)
+    {
+        const int atomic_number = atomic_number_from_symbol(ucell_.atoms[it].ncpp.psd);
+
+        for (int ia = 0; ia < ucell_.atoms[it].na; ++ia)
+        {
+            const ModuleBase::Vector3<double> position = ucell_.atoms[it].tau[ia] * ucell_.lat0;
+
+            numbers.push_back(atomic_number);
+            positions.push_back(position.x);
+            positions.push_back(position.y);
+            positions.push_back(position.z);
+        }
+    }
+
+    const ModuleBase::Vector3<double> a1 = ucell_.a1 * ucell_.lat0;
+    const ModuleBase::Vector3<double> a2 = ucell_.a2 * ucell_.lat0;
+    const ModuleBase::Vector3<double> a3 = ucell_.a3 * ucell_.lat0;
+
+    // DFT-D4 C API documents lattice as lattice[3][3] in Bohr.
+    // Keep this order aligned with positions[natoms][3], and verify it with
+    // a non-orthogonal-cell regression test before relying on stress values.
+    lattice.push_back(a1.x);
+    lattice.push_back(a1.y);
+    lattice.push_back(a1.z);
+    lattice.push_back(a2.x);
+    lattice.push_back(a2.y);
+    lattice.push_back(a2.z);
+    lattice.push_back(a3.x);
+    lattice.push_back(a3.y);
+    lattice.push_back(a3.z);
+
+    periodic = {{true, true, true}};
+}
+
+void Vdwd4::compute(double& energy_ha,
+                    std::vector<double>* gradient_ha_bohr,
+                    std::array<double, 9>* sigma_ha)
+{
+#ifdef __DFTD4
+    std::vector<int> numbers;
+    std::vector<double> positions;
+    std::vector<double> lattice;
+    std::array<bool, 3> periodic;
+
+    build_structure(numbers, positions, lattice, periodic);
+
+    if (gradient_ha_bohr != nullptr
+        && gradient_ha_bohr->size() != static_cast<std::size_t>(3 * ucell_.nat))
+    {
+        ModuleBase::WARNING_QUIT("Vdwd4::compute",
+                                 "gradient_ha_bohr must have size 3 * nat when requested.");
+    }
+
+    // These vectors own all arrays passed to DFT-D4. Their data() pointers
+    // remain valid until all DFT-D4 handles created below are deleted.
+    dftd4_error error = dftd4_new_error();
+
+    dftd4_structure mol = dftd4_new_structure(error,
+                                              ucell_.nat,
+                                              numbers.data(),
+                                              positions.data(),
+                                              nullptr,
+                                              lattice.data(),
+                                              periodic.data());
+    check_dftd4_error(error, "dftd4_new_structure");
+
+    dftd4_model model = dftd4_new_d4_model(error, mol);
+    check_dftd4_error(error, "dftd4_new_d4_model");
+
+    dftd4_set_model_realspace_cutoff(error, model, cutoff_disp2_, cutoff_disp3_, cutoff_cn_);
+    check_dftd4_error(error, "dftd4_set_model_realspace_cutoff");
+
+    std::vector<char> method(xc_name_.begin(), xc_name_.end());
+    method.push_back('\0');
+
+    // Use the DFT-D4 library's internal rational damping parameters.
+    // The final boolean selects the three-body-specific parameterization in the C API.
+    const bool atm = true;
+    dftd4_param param = dftd4_load_rational_damping(error, method.data(), atm);
+    check_dftd4_error(error, "dftd4_load_rational_damping");
+
+    // These are borrowed output buffers for this synchronous C API call only.
+    double* gradient = gradient_ha_bohr ? gradient_ha_bohr->data() : nullptr;
+    double* sigma = sigma_ha ? sigma_ha->data() : nullptr;
+
+    dftd4_get_dispersion(error, mol, model, param, &energy_ha, gradient, sigma);
+    check_dftd4_error(error, "dftd4_get_dispersion");
+
+    dftd4_delete_param(&param);
+    dftd4_delete_model(&model);
+    dftd4_delete_structure(&mol);
+    dftd4_delete_error(&error);
+#else
+    ModuleBase::WARNING_QUIT("Vdwd4::compute", "DFT-D4 correction requires ABACUS to be built with DFTD4.");
+#endif
+}
+
+void Vdwd4::cal_energy()
+{
+    ModuleBase::TITLE("Vdwd4", "cal_energy");
+    ModuleBase::timer::start("Vdwd4", "cal_energy");
+
+    double energy_ha = 0.0;
+    compute(energy_ha, nullptr, nullptr);
+
+    // DFT-D4 returns Hartree; ABACUS vdW energies are stored in Ry.
+    energy_ = 2.0 * energy_ha;
+
+    ModuleBase::timer::end("Vdwd4", "cal_energy");
+}
+
+void Vdwd4::set_force_from_gradient(const std::vector<double>& gradient_ha_bohr)
+{
+    force_.clear();
+    force_.resize(ucell_.nat);
+
+    for (int iat = 0; iat < ucell_.nat; ++iat)
+    {
+        // DFT-D4 returns dE/dR in Ha/Bohr; ABACUS forces are -dE/dR in Ry/Bohr.
+        force_[iat].x = -2.0 * gradient_ha_bohr[3 * iat + 0];
+        force_[iat].y = -2.0 * gradient_ha_bohr[3 * iat + 1];
+        force_[iat].z = -2.0 * gradient_ha_bohr[3 * iat + 2];
+    }
+
+    has_force_cache_ = true;
+}
+
+void Vdwd4::set_stress_from_sigma(const std::array<double, 9>& sigma_ha)
+{
+    // Tentative mapping consistent with the current D3 convention.
+    // Confirm sign, transposition and volume normalization by finite-strain tests.
+    stress_ = ModuleBase::Matrix3(2.0 * sigma_ha[0], 2.0 * sigma_ha[1], 2.0 * sigma_ha[2],
+                                  2.0 * sigma_ha[3], 2.0 * sigma_ha[4], 2.0 * sigma_ha[5],
+                                  2.0 * sigma_ha[6], 2.0 * sigma_ha[7], 2.0 * sigma_ha[8])
+              / ucell_.omega;
+
+    has_stress_cache_ = true;
+}
+
+void Vdwd4::cal_force()
+{
+    ModuleBase::TITLE("Vdwd4", "cal_force");
+    ModuleBase::timer::start("Vdwd4", "cal_force");
+
+    if (!has_force_cache_ || !has_stress_cache_)
+    {
+        double energy_ha = 0.0;
+        std::vector<double> gradient(3 * ucell_.nat, 0.0);
+        std::array<double, 9> sigma;
+        sigma.fill(0.0);
+
+        // Request sigma together with the gradient.  The DFT-D4 C API computes
+        // sigma internally for gradient calculations anyway, so keep it and
+        // avoid a second expensive D4 call when ABACUS subsequently requests stress.
+        compute(energy_ha, &gradient, &sigma);
+        set_force_from_gradient(gradient);
+        set_stress_from_sigma(sigma);
+    }
+
+    ModuleBase::timer::end("Vdwd4", "cal_force");
+}
+
+void Vdwd4::cal_stress()
+{
+    ModuleBase::TITLE("Vdwd4", "cal_stress");
+    ModuleBase::timer::start("Vdwd4", "cal_stress");
+
+    if (!has_stress_cache_)
+    {
+        double energy_ha = 0.0;
+        std::array<double, 9> sigma;
+        sigma.fill(0.0);
+
+        compute(energy_ha, nullptr, &sigma);
+        set_stress_from_sigma(sigma);
+    }
+
+    ModuleBase::timer::end("Vdwd4", "cal_stress");
+}
+
+} // namespace vdw
diff --git a/source/source_hamilt/module_vdw/vdwd4.h b/source/source_hamilt/module_vdw/vdwd4.h
new file mode 100644
index 00000000000..3772c8646d2
--- /dev/null
+++ b/source/source_hamilt/module_vdw/vdwd4.h
@@ -0,0 +1,50 @@
+#ifndef VDWD4_H
+#define VDWD4_H
+
+#include "vdw.h"
+
+#include <array>
+#include <string>
+#include <vector>
+
+namespace vdw
+{
+
+class Vdwd4 : public Vdw
+{
+  public:
+    Vdwd4(const UnitCell& unit_in, const std::string& xc_name, const Input_para& input);
+
+    ~Vdwd4() override = default;
+
+  private:
+    std::string xc_name_;
+    double cutoff_disp2_ = 0.0; // Bohr, two-body dispersion cutoff
+    double cutoff_disp3_ = 0.0; // Bohr, three-body ATM cutoff
+    double cutoff_cn_ = 0.0;    // Bohr, coordination-number cutoff
+
+    bool has_force_cache_ = false;
+    bool has_stress_cache_ = false;
+
+    void set_force_from_gradient(const std::vector<double>& gradient_ha_bohr);
+    void set_stress_from_sigma(const std::array<double, 9>& sigma_ha);
+
+    void cal_energy() override;
+    void cal_force() override;
+    void cal_stress() override;
+
+    void build_structure(std::vector<int>& numbers,
+                         std::vector<double>& positions,
+                         std::vector<double>& lattice,
+                         std::array<bool, 3>& periodic) const;
+
+    // Optional output buffers are caller-owned and non-owning here.
+    // Vdwd4 writes to them during compute() and never stores their pointers.
+    void compute(double& energy_ha,
+                 std::vector<double>* gradient_ha_bohr,
+                 std::array<double, 9>* sigma_ha);
+};
+
+} // namespace vdw
+
+#endif // VDWD4_H
diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h
index 029ad364eb5..d953b542da5 100644
--- a/source/source_io/module_parameter/input_parameter.h
+++ b/source/source_io/module_parameter/input_parameter.h
@@ -509,7 +509,7 @@ struct Input_para
     //  vdw
     //  Peize Lin add 2014-03-31, jiyy update 2019-08-01
     // ==========================================================
-    std::string vdw_method = "none";                        ///< the method of calculating vdw (none; d2; d3_0; d3_bj)
+    std::string vdw_method = "none";                        ///< the method of calculating vdw (none; d2; d3_0; d3_bj; d4)
     std::string vdw_s6 = "default";                         ///< scale parameter of d2/d3_0/d3_bj
     std::string vdw_s8 = "default";                         ///< scale parameter of d3_0/d3_bj
     std::string vdw_a1 = "default";                         ///< damping parameter of d3_0/d3_bj
@@ -526,6 +526,7 @@ struct Input_para
     std::string vdw_radius_unit = "Bohr";                   ///< unit of radius cutoff for periodic structure
     double vdw_cn_thr = 40.0;                               ///< radius cutoff for cn
     std::string vdw_cn_thr_unit = "Bohr";                   ///< unit of cn_thr, Bohr or Angstrom
+    std::string vdw_d4_xc = "default";                     ///< functional name passed to DFT-D4
     ModuleBase::Vector3<int> vdw_cutoff_period = {3, 3, 3}; ///< periods of periodic structure
 
     // ==============   #Parameters (15.exx) ====================
diff --git a/source/source_io/module_parameter/read_input_item_model.cpp b/source/source_io/module_parameter/read_input_item_model.cpp
index e6c69139d9b..a3bd0b0a9f8 100644
--- a/source/source_io/module_parameter/read_input_item_model.cpp
+++ b/source/source_io/module_parameter/read_input_item_model.cpp
@@ -253,13 +253,14 @@ void ReadInput::item_model()
     // vdW Correction
     {
         Input_Item item("vdw_method");
-        item.annotation = "the method of calculating vdw (none ; d2 ; d3_0 ; d3_bj";
+        item.annotation = "the method of calculating vdw (none ; d2 ; d3_0 ; d3_bj ; d4)";
         item.category = "vdW correction";
         item.type = "String";
         item.description = R"(Specifies the method used for Van der Waals (VdW) correction. Available options are:
 * d2: Grimme's D2 dispersion correction method
 * d3_0: Grimme's DFT-D3(0) dispersion correction method (zero-damping)
 * d3_bj: Grimme's DFTD3(BJ) dispersion correction method (BJ-damping)
+* d4: Grimme's DFT-D4 dispersion correction method using the external DFT-D4 library
 * none: no vdW correction
 
 [NOTE] ABACUS supports automatic setting of DFT-D3 parameters for common functionals. To benefit from this feature, please specify the parameter dft_functional explicitly, otherwise the autoset procedure will crash. If not satisfied with the built-in parameters, any manual setting on vdw_s6, vdw_s8, vdw_a1 and vdw_a2 will overwrite the automatic values.)";
@@ -269,6 +270,19 @@ void ReadInput::item_model()
         read_sync_string(input.vdw_method);
         this->add_item(item);
     }
+    {
+        Input_Item item("vdw_d4_xc");
+        item.annotation = "functional name passed to DFT-D4";
+        item.category = "vdW correction";
+        item.type = "String";
+        item.description = R"(Functional name used to load DFT-D4 damping parameters from the DFT-D4 library.
+If set to default, ABACUS infers the functional name from dft_functional or pseudopotential metadata.)";
+        item.default_value = "default";
+        item.unit = "";
+        item.availability = "vdw_method is set to d4";
+        read_sync_string(input.vdw_d4_xc);
+        this->add_item(item);
+    }
     {
         Input_Item item("vdw_s6");
         item.annotation = "scale parameter of d2/d3_0/d3_bj";
@@ -507,6 +521,13 @@ Namely, each line contains the element name and the corresponding parameter.)";
                 {
                     para.input.vdw_cutoff_radius = "95";
                 }
+                else if (para.input.vdw_method == "d4")
+                {
+                    // DFT-D4 uses separate real-space cutoffs internally.
+                    // This input controls the two-body cutoff; the wrapper keeps
+                    // the three-body cutoff at min(40 Bohr, vdw_cutoff_radius).
+                    para.input.vdw_cutoff_radius = "60";
+                }
                 else
                 {
                     para.input.vdw_cutoff_radius = "0";
@@ -589,7 +610,14 @@ Namely, each line contains the element name and the corresponding parameter.)";
         item.description = "The cutoff radius when calculating coordination numbers.";
         item.default_value = "40";
         item.unit = "defined by vdw_cn_thr_unit (default: Bohr)";
-        item.availability = "vdw_method is set to d3_0 or d3_bj";
+        item.availability = "vdw_method is set to d3_0, d3_bj, or d4";
+        item.reset_value = [](const Input_Item& item, Parameter& para) {
+            if (!item.is_read() && para.input.vdw_method == "d4")
+            {
+                // DFT-D4 library default for coordination numbers.
+                para.input.vdw_cn_thr = 30.0;
+            }
+        };
         read_sync_double(input.vdw_cn_thr);
         item.check_value = [](const Input_Item& item, const Parameter& para) {
             if (para.input.vdw_cn_thr <= 0)
diff --git a/source/source_lcao/FORCE_STRESS.cpp b/source/source_lcao/FORCE_STRESS.cpp
index f9f493eaf1d..83d09c4f78c 100644
--- a/source/source_lcao/FORCE_STRESS.cpp
+++ b/source/source_lcao/FORCE_STRESS.cpp
@@ -768,7 +768,7 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
             ModuleIO::print_stress("XC       STRESS", sigmaxc, screen, ry, GlobalV::ofs_running);
             if (vdw_solver != nullptr)
             {
-                ModuleIO::print_stress("VDW      STRESS", sigmaxc, screen, ry, GlobalV::ofs_running);
+                ModuleIO::print_stress("VDW      STRESS", stress_vdw, screen, ry, GlobalV::ofs_running);
             }
             if (PARAM.inp.dft_plus_u)
             {
diff --git a/source/source_pw/module_ofdft/of_print_info.cpp b/source/source_pw/module_ofdft/of_print_info.cpp
index fa19083dcdc..1f3a14b7910 100644
--- a/source/source_pw/module_ofdft/of_print_info.cpp
+++ b/source/source_pw/module_ofdft/of_print_info.cpp
@@ -79,6 +79,11 @@ void OFDFT::print_info(const int iter,
             titles.push_back("E_vdwD3");
             energies_Ry.push_back(pelec->f_en.evdw);
         }
+        else if (vdw_method == "d4")
+        {
+            titles.push_back("E_vdwD4");
+            energies_Ry.push_back(pelec->f_en.evdw);
+        }
         if (PARAM.inp.imp_sol)
         {
             titles.push_back("E_sol_el");
diff --git a/toolchain/build_abacus_gcc-mkl.sh b/toolchain/build_abacus_gcc-mkl.sh
index 6bd0123bcb8..df71a0d65ff 100755
--- a/toolchain/build_abacus_gcc-mkl.sh
+++ b/toolchain/build_abacus_gcc-mkl.sh
@@ -70,6 +70,7 @@ cmake -B $BUILD_DIR -DCMAKE_INSTALL_PREFIX=$PREFIX \
         -DENABLE_LIBXC=ON \
         -DUSE_OPENMP=ON \
         -DUSE_ELPA=ON \
+        -DENABLE_DFTD4=ON \
         -DENABLE_RAPIDJSON=ON \
         -DRapidJSON_DIR=$RAPIDJSON \
         -DENABLE_LIBRI=ON \
diff --git a/toolchain/build_abacus_gnu.sh b/toolchain/build_abacus_gnu.sh
index 09e6ede052b..f621fbe2cc0 100755
--- a/toolchain/build_abacus_gnu.sh
+++ b/toolchain/build_abacus_gnu.sh
@@ -72,6 +72,7 @@ cmake -B $BUILD_DIR -DCMAKE_INSTALL_PREFIX=$PREFIX \
         -DENABLE_LIBXC=ON \
         -DUSE_OPENMP=ON \
         -DUSE_ELPA=ON \
+        -DENABLE_DFTD4=ON \
         -DENABLE_RAPIDJSON=ON \
         -DRapidJSON_DIR=$RAPIDJSON \
         -DENABLE_LIBRI=ON \
diff --git a/toolchain/build_abacus_intel.sh b/toolchain/build_abacus_intel.sh
index f4b282964be..f034988522a 100755
--- a/toolchain/build_abacus_intel.sh
+++ b/toolchain/build_abacus_intel.sh
@@ -71,6 +71,7 @@ cmake -B $BUILD_DIR -DCMAKE_INSTALL_PREFIX=$PREFIX \
         -DENABLE_LIBXC=ON \
         -DUSE_OPENMP=ON \
         -DUSE_ELPA=ON \
+        -DENABLE_DFTD4=ON \
         -DENABLE_RAPIDJSON=ON \
         -DRapidJSON_DIR=$RAPIDJSON \
 	    -DENABLE_LIBRI=ON \
diff --git a/toolchain/scripts/lib/config_manager.sh b/toolchain/scripts/lib/config_manager.sh
index e0228e81262..65ea5740d26 100644
--- a/toolchain/scripts/lib/config_manager.sh
+++ b/toolchain/scripts/lib/config_manager.sh
@@ -19,7 +19,7 @@ CONFIG_FILE_LOADED=false
 tool_list="gcc intel amd cmake"
 mpi_list="mpich openmpi intelmpi"
 math_list="mkl aocl openblas"
-lib_list="fftw libxc scalapack elpa cereal rapidjson libtorch libnpy libri libcomm nep"
+lib_list="fftw libxc scalapack elpa cereal rapidjson libtorch libnpy libri libcomm nep dftd4"
 package_list="${tool_list} ${mpi_list} ${math_list} ${lib_list}"
 
 # Configuration file paths for loading in advance
@@ -371,6 +371,7 @@ config_set_defaults() {
     CONFIG_CACHE["with_libri"]="__INSTALL__"
     CONFIG_CACHE["with_libcomm"]="__INSTALL__"
     CONFIG_CACHE["with_nep"]="__DONTUSE__"
+    CONFIG_CACHE["with_dftd4"]="__DONTUSE__"
     
     # Default enable options (following original script logic)
     CONFIG_CACHE["dry_run"]="__FALSE__"
diff --git a/toolchain/scripts/lib/package_manager.sh b/toolchain/scripts/lib/package_manager.sh
index 94a7c0b3079..71432578d81 100644
--- a/toolchain/scripts/lib/package_manager.sh
+++ b/toolchain/scripts/lib/package_manager.sh
@@ -60,6 +60,7 @@ package_manager_define_dependencies() {
     PACKAGE_DEPENDENCIES["elpa"]="scalapack"
     
     # Stage 4: Advanced Libraries
+    PACKAGE_DEPENDENCIES["dftd4"]="cmake"
     PACKAGE_DEPENDENCIES["cereal"]="gcc cmake"
     PACKAGE_DEPENDENCIES["rapidjson"]="gcc cmake"
     PACKAGE_DEPENDENCIES["libtorch"]="gcc cmake"
@@ -466,4 +467,4 @@ package_install_all_stages() {
     
     echo "All stages completed successfully."
     return 0
-}
\ No newline at end of file
+}
diff --git a/toolchain/scripts/lib/version_helper.sh b/toolchain/scripts/lib/version_helper.sh
index 46e7436cad0..d5bde1cfb69 100644
--- a/toolchain/scripts/lib/version_helper.sh
+++ b/toolchain/scripts/lib/version_helper.sh
@@ -335,4 +335,4 @@ export -f version_get_effective
 export -f version_load_package_vars
 export -f version_show_current
 export -f version_validate_config
-export -f version_helper_init
\ No newline at end of file
+export -f version_helper_init
diff --git a/toolchain/scripts/package_versions.sh b/toolchain/scripts/package_versions.sh
index d758c8ea0f3..6379ee8b2ba 100644
--- a/toolchain/scripts/package_versions.sh
+++ b/toolchain/scripts/package_versions.sh
@@ -84,6 +84,10 @@ scalapack_alt_sha256="a2f0c9180a210bf7ffe126c9cb81099cf337da1a7120ddb4cbe4894eb7
 # STAGE 4: Advanced Feature Libraries
 # =============================================================================
 
+# DFT-D4 dispersion correction
+dftd4_ver="4.2.0"
+dftd4_sha256="467e024071510ad82b862c66c383c2ebc164fc1140e15dfc79f48d2f999fd184"
+
 # LibTorch (supports dual versions) - main=2.1.2, alt=1.12.1
 libtorch_main_ver="2.1.2"
 libtorch_main_sha256="904b764df6106a8a35bef64c4b55b8c1590ad9d071eb276e680cf42abafe79e9"
@@ -222,6 +226,10 @@ load_package_vars() {
                 scalapack_sha256="${scalapack_main_sha256}"
             fi
             ;;
+        "dftd4")
+            dftd4_ver="${dftd4_ver}"
+            dftd4_sha256="${dftd4_sha256}"
+            ;;
         "libtorch")
             if [ "${version_suffix}" = "alt" ]; then
                 libtorch_ver="${libtorch_alt_ver}"
diff --git a/toolchain/scripts/stage4/install_dftd4.sh b/toolchain/scripts/stage4/install_dftd4.sh
new file mode 100755
index 00000000000..b9595226b84
--- /dev/null
+++ b/toolchain/scripts/stage4/install_dftd4.sh
@@ -0,0 +1,107 @@
+#!/bin/bash -e
+
+# TODO: Review and if possible fix shellcheck errors.
+# shellcheck disable=all
+
+[ "${BASH_SOURCE[0]}" ] && SCRIPT_NAME="${BASH_SOURCE[0]}" || SCRIPT_NAME=$0
+SCRIPT_DIR="$(cd "$(dirname "$SCRIPT_NAME")/.." && pwd -P)"
+
+source "${SCRIPT_DIR}"/common_vars.sh
+source "${SCRIPT_DIR}"/tool_kit.sh
+source "${SCRIPT_DIR}"/signal_trap.sh
+source "${SCRIPT_DIR}"/package_versions.sh
+
+# Load package variables with appropriate version
+load_package_vars "dftd4"
+dftd4_pkg="dftd4-${dftd4_ver}.tar.xz"
+source "${INSTALLDIR}"/toolchain.conf
+source "${INSTALLDIR}"/toolchain.env
+
+[ -f "${BUILDDIR}/setup_dftd4" ] && rm "${BUILDDIR}/setup_dftd4"
+
+! [ -d "${BUILDDIR}" ] && mkdir -p "${BUILDDIR}"
+cd "${BUILDDIR}"
+
+case "${with_dftd4}" in
+    __INSTALL__)
+        echo "==================== Installing DFT-D4 ===================="
+        url="https://github.com/dftd4/dftd4/releases/download/v${dftd4_ver}/${dftd4_pkg}"
+        pkg_install_dir="${INSTALLDIR}/dftd4-${dftd4_ver}"
+        install_lock_file="${pkg_install_dir}/install_successful"
+
+        if verify_checksums "${install_lock_file}"; then
+            echo "dftd4-${dftd4_ver} is already installed, skipping it."
+        else
+            retrieve_package "${dftd4_sha256}" "${dftd4_pkg}" "${url}"
+            if [ "${PACK_RUN}" = "__TRUE__" ]; then
+                echo "--pack-run mode specified, skip installation"
+                exit 0
+            fi
+
+            echo "Installing from scratch into ${pkg_install_dir}"
+            [ -d "dftd4-${dftd4_ver}" ] && rm -rf "dftd4-${dftd4_ver}"
+            tar -xJf "${dftd4_pkg}"
+            cd dftd4-${dftd4_ver}
+
+            mkdir build && cd build
+            CMAKE_PREFIX_PATH="${CMAKE_PREFIX_PATH}:${OPENBLAS_ROOT}" cmake \
+                -DCMAKE_INSTALL_PREFIX="${pkg_install_dir}" \
+                -DCMAKE_INSTALL_LIBDIR=lib \
+                -DCMAKE_BUILD_TYPE=Release \
+                -DCMAKE_VERBOSE_MAKEFILE=ON \
+                .. > cmake.log 2>&1 || tail -n ${LOG_LINES} cmake.log
+            make install -j $(get_nprocs) > make.log 2>&1 || tail -n ${LOG_LINES} make.log
+            cd ..
+            write_checksums "${install_lock_file}" "${SCRIPT_DIR}/stage4/$(basename ${SCRIPT_NAME})"
+        fi
+        ;;
+
+    __SYSTEM__)
+        echo "==================== Finding DFT-D4 from system paths ===================="
+        if [ "${PACK_RUN}" = "__TRUE__" ]; then
+            echo "--pack-run mode specified, skip system check"
+            exit 0
+        fi
+        check_command pkg-config --modversion dftd4
+        pkg_install_dir="$(dirname $(dirname $(find_in_paths "libdftd4.*" $LIB_PATHS)))"
+        if [ -d "${pkg_install_dir}/lib/cmake/dftd4" ] ||
+          [ -d "${pkg_install_dir}/lib64/cmake/dftd4" ]; then
+          echo "Package dftd4 is found and confirmed to be built with CMake."
+        else
+          echo "ERROR: ABACUS requires dftd4 to be built with CMake."
+          exit 1
+        fi
+        ;;
+
+    __DONTUSE__) ;;
+
+    *)
+        echo "==================== Linking DFT-D4 to user paths ===================="
+        pkg_install_dir="$with_dftd4"
+        check_dir "${pkg_install_dir}/include"
+        ;;
+esac
+
+if [ "${with_dftd4}" != "__DONTUSE__" ]; then
+    cat << EOF > "${BUILDDIR}/setup_dftd4"
+export DFTD4_VER="${dftd4_ver}"
+export DFTD4_ROOT="${pkg_install_dir}"
+EOF
+    if [ "${with_dftd4}" != "__SYSTEM__" ]; then
+        cat << EOF >> "${BUILDDIR}/setup_dftd4"
+prepend_path PATH "${pkg_install_dir}/bin"
+prepend_path LD_LIBRARY_PATH "${pkg_install_dir}/lib"
+prepend_path LD_RUN_PATH "${pkg_install_dir}/lib"
+prepend_path LIBRARY_PATH "${pkg_install_dir}/lib"
+prepend_path CPATH "${pkg_install_dir}/include"
+prepend_path CMAKE_PREFIX_PATH "${pkg_install_dir}"
+EOF
+    fi
+    filter_setup "${BUILDDIR}/setup_dftd4" "${SETUPFILE}"
+fi
+
+load "${BUILDDIR}/setup_dftd4"
+write_toolchain_env "${INSTALLDIR}"
+
+cd "${ROOTDIR}"
+report_timing "dftd4"
diff --git a/toolchain/scripts/stage4/install_stage4.sh b/toolchain/scripts/stage4/install_stage4.sh
index d2e3fa482a6..0980d809e5e 100755
--- a/toolchain/scripts/stage4/install_stage4.sh
+++ b/toolchain/scripts/stage4/install_stage4.sh
@@ -3,6 +3,7 @@
 # TODO: Review and if possible fix shellcheck errors.
 # shellcheck disable=all
 
+./scripts/stage4/install_dftd4.sh
 ./scripts/stage4/install_cereal.sh
 ./scripts/stage4/install_rapidjson.sh
 ./scripts/stage4/install_libtorch.sh
diff --git a/toolchain/toolchain_gcc-mkl.sh b/toolchain/toolchain_gcc-mkl.sh
index 7f7b9da8e17..bac479c38f7 100755
--- a/toolchain/toolchain_gcc-mkl.sh
+++ b/toolchain/toolchain_gcc-mkl.sh
@@ -44,6 +44,9 @@ WITH_LIBTORCH="no"
 WITH_LIBNPY="no"
 WITH_NEP="no"
 
+# Optional Features (DFT-D4 dispersion correction)
+WITH_DFTD4="install"
+
 # ELPA-GPU Support (uncomment and modify as needed)
 # ENABLE_CUDA="yes"
 # GPU_VERSION="75"  # Check your GPU compute capability
@@ -93,6 +96,7 @@ exec ./install_abacus_toolchain_new.sh \
   --with-libxc="$WITH_LIBXC" \
   --with-fftw="$WITH_FFTW" \
   --with-elpa="$WITH_ELPA" \
+  --with-dftd4="$WITH_DFTD4" \
   --with-cereal="$WITH_CEREAL" \
   --with-rapidjson="$WITH_RAPIDJSON" \
   --with-libtorch="$WITH_LIBTORCH" \
diff --git a/toolchain/toolchain_gnu.sh b/toolchain/toolchain_gnu.sh
index 456d81dacc5..a934fed0181 100755
--- a/toolchain/toolchain_gnu.sh
+++ b/toolchain/toolchain_gnu.sh
@@ -44,6 +44,9 @@ WITH_LIBTORCH="no"
 WITH_LIBNPY="no"
 WITH_NEP="no"
 
+# Optional Features (DFT-D4 dispersion correction)
+WITH_DFTD4="install"
+
 # ELPA-GPU Support (uncomment and modify as needed)
 # ENABLE_CUDA="yes"
 # GPU_VERSION="75"  # Check your GPU compute capability
@@ -95,6 +98,7 @@ exec ./install_abacus_toolchain_new.sh \
   --with-libxc="$WITH_LIBXC" \
   --with-fftw="$WITH_FFTW" \
   --with-elpa="$WITH_ELPA" \
+  --with-dftd4="$WITH_DFTD4" \
   --with-cereal="$WITH_CEREAL" \
   --with-rapidjson="$WITH_RAPIDJSON" \
   --with-libtorch="$WITH_LIBTORCH" \
diff --git a/toolchain/toolchain_intel.sh b/toolchain/toolchain_intel.sh
index 1f4dedbc215..ab40cf55b2b 100755
--- a/toolchain/toolchain_intel.sh
+++ b/toolchain/toolchain_intel.sh
@@ -45,6 +45,9 @@ WITH_LIBTORCH="no"
 WITH_LIBNPY="no"
 WITH_NEP="no"
 
+# Optional Features (DFT-D4 dispersion correction)
+WITH_DFTD4="install"
+
 # ELPA-GPU Support (uncomment and modify as needed)
 # ENABLE_CUDA="yes"
 # GPU_VERSION="75"  # Check your GPU compute capability
@@ -114,6 +117,7 @@ exec ./install_abacus_toolchain_new.sh \
   --with-libxc="$WITH_LIBXC" \
   --with-fftw="$WITH_FFTW" \
   --with-elpa="$WITH_ELPA" \
+  --with-dftd4="$WITH_DFTD4" \
   --with-cereal="$WITH_CEREAL" \
   --with-rapidjson="$WITH_RAPIDJSON" \
   --with-libtorch="$WITH_LIBTORCH" \

From f0c81210dc3841e7526eecd2a34a6d1d116db30a Mon Sep 17 00:00:00 2001
From: Growl <uwsy1059@qq.com>
Date: Tue, 26 May 2026 09:47:38 +0800
Subject: [PATCH 2/7] Docs and tests

---
 docs/advanced/input_files/input-main.md       | 16 ++++-
 docs/advanced/install.md                      | 18 +++++
 docs/parameters.yaml                          | 17 ++++-
 .../module_vdw/test/CMakeLists.txt            |  5 ++
 .../module_vdw/test/vdw_test.cpp              | 70 +++++++++++++++++++
 5 files changed, 120 insertions(+), 6 deletions(-)

diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
index b5f458c21ea..740b4f56a96 100644
--- a/docs/advanced/input_files/input-main.md
+++ b/docs/advanced/input_files/input-main.md
@@ -3640,11 +3640,21 @@
   - d2: Grimme's D2 dispersion correction method
   - d3_0: Grimme's DFT-D3(0) dispersion correction method (zero-damping)
   - d3_bj: Grimme's DFTD3(BJ) dispersion correction method (BJ-damping)
+  - d4: Grimme's DFT-D4 dispersion correction method using the external DFT-D4 library
   - none: no vdW correction
 
   > Note: ABACUS supports automatic setting of DFT-D3 parameters for common functionals. To benefit from this feature, please specify the parameter dft_functional explicitly, otherwise the autoset procedure will crash. If not satisfied with the built-in parameters, any manual setting on vdw_s6, vdw_s8, vdw_a1 and vdw_a2 will overwrite the automatic values.
+
+  > Note: DFT-D4 support requires ABACUS to be configured with ENABLE_DFTD4=ON and a CMake-installed dftd4 library exporting dftd4-config.cmake. DFT-D4 damping parameters are loaded from the external library.
 - **Default**: none
 
+### vdw_d4_xc
+
+- **Type**: String
+- **Availability**: *vdw_method is set to d4*
+- **Description**: Functional name passed to the DFT-D4 library to load its internal damping parameters. If set to default, ABACUS infers the functional name from dft_functional or pseudopotential metadata.
+- **Default**: default
+
 ### vdw_s6
 
 - **Type**: String
@@ -3737,7 +3747,7 @@
 
 - **Type**: String
 - **Availability**: *vdw_cutoff_type is set to radius*
-- **Description**: Defines the radius of the cutoff sphere when vdw_cutoff_type is set to radius. The default values depend on the chosen vdw_method.
+- **Description**: Defines the cutoff radius when vdw_cutoff_type is set to radius. The default values depend on the chosen vdw_method. For DFT-D4, this controls the two-body dispersion cutoff, while the three-body cutoff is internally limited to the DFT-D4 default value of 40 Bohr.
 - **Unit**: defined by vdw_radius_unit (default Bohr)
 
 ### vdw_radius_unit
@@ -3759,8 +3769,8 @@
 ### vdw_cn_thr
 
 - **Type**: Real
-- **Availability**: *vdw_method is set to d3_0 or d3_bj*
-- **Description**: The cutoff radius when calculating coordination numbers.
+- **Availability**: *vdw_method is set to d3_0, d3_bj, or d4*
+- **Description**: The cutoff radius when calculating coordination numbers. The default is 40 Bohr for DFT-D3 and 30 Bohr for DFT-D4.
 - **Default**: 40
 - **Unit**: defined by vdw_cn_thr_unit (default: Bohr)
 
diff --git a/docs/advanced/install.md b/docs/advanced/install.md
index 0df93f602a5..654e853e204 100644
--- a/docs/advanced/install.md
+++ b/docs/advanced/install.md
@@ -77,6 +77,24 @@ These two libraries are added as submodules in the [deps](https://github.com/dee
 
 If you prefer using manually downloaded libraries, provide `-DLIBRI_DIR=${path to your LibRI folder} -DLIBCOMM_DIR=${path to your LibComm folder}`.
 
+
+## Build with DFT-D4 support
+
+ABACUS can use the external [DFT-D4](https://github.com/dftd4/dftd4) library for Grimme's DFT-D4 dispersion correction. DFT-D4 support is optional and disabled by default.
+
+DFT-D4 must be built and installed with CMake, so that ABACUS can locate it through the exported `dftd4-config.cmake` package. Meson-built installations, including the Conda package of `dftd4`, are not supported unless they provide a compatible CMake package file.
+
+To build ABACUS with DFT-D4 support, pass `-DENABLE_DFTD4=ON` to CMake and provide `CMAKE_PREFIX_PATH` environment variable.
+
+In the input file, enable DFT-D4 with:
+
+```text
+vdw_method d4
+vdw_d4_xc pbe
+```
+
+If `vdw_d4_xc` is set to `default`, ABACUS will infer the functional name from `dft_functional` or pseudopotential metadata and pass it to the DFT-D4 library.
+
 ## Build Unit Tests
 
 To build tests for ABACUS, define `BUILD_TESTING` flag. You can also specify path to local installation of [Googletest](https://github.com/google/googletest) by setting `GTEST_DIR` flags. If not found in local, the configuration process will try to download it automatically.
diff --git a/docs/parameters.yaml b/docs/parameters.yaml
index 63ee83376c5..bb9d365b5ee 100644
--- a/docs/parameters.yaml
+++ b/docs/parameters.yaml
@@ -3738,12 +3738,23 @@ parameters:
       * d2: Grimme's D2 dispersion correction method
       * d3_0: Grimme's DFT-D3(0) dispersion correction method (zero-damping)
       * d3_bj: Grimme's DFTD3(BJ) dispersion correction method (BJ-damping)
+      * d4: Grimme's DFT-D4 dispersion correction method using the external DFT-D4 library
       * none: no vdW correction
 
       [NOTE] ABACUS supports automatic setting of DFT-D3 parameters for common functionals. To benefit from this feature, please specify the parameter dft_functional explicitly, otherwise the autoset procedure will crash. If not satisfied with the built-in parameters, any manual setting on vdw_s6, vdw_s8, vdw_a1 and vdw_a2 will overwrite the automatic values.
+
+      [NOTE] DFT-D4 support requires ABACUS to be configured with ENABLE_DFTD4=ON and a CMake-installed dftd4 library exporting `dftd4-config.cmake`. DFT-D4 damping parameters are loaded from the external library.
     default_value: none
     unit: ""
     availability: ""
+  - name: vdw_d4_xc
+    category: vdW correction
+    type: String
+    description: |
+      Functional name passed to the DFT-D4 library to load its internal damping parameters. If set to default, ABACUS infers the functional name from dft_functional or pseudopotential metadata.
+    default_value: default
+    unit: ""
+    availability: vdw_method is set to d4
   - name: vdw_s6
     category: vdW correction
     type: String
@@ -3852,7 +3863,7 @@ parameters:
     category: vdW correction
     type: String
     description: |
-      Defines the radius of the cutoff sphere when vdw_cutoff_type is set to radius. The default values depend on the chosen vdw_method.
+      Defines the cutoff radius when vdw_cutoff_type is set to radius. The default values depend on the chosen vdw_method. For DFT-D4, this controls the two-body dispersion cutoff, while the three-body cutoff is internally limited to the DFT-D4 default value of 40 Bohr.
     default_value: ""
     unit: defined by vdw_radius_unit (default Bohr)
     availability: vdw_cutoff_type is set to radius
@@ -3878,10 +3889,10 @@ parameters:
     category: vdW correction
     type: Real
     description: |
-      The cutoff radius when calculating coordination numbers.
+      The cutoff radius when calculating coordination numbers. The default is 40 Bohr for DFT-D3 and 30 Bohr for DFT-D4.
     default_value: "40"
     unit: "defined by vdw_cn_thr_unit (default: Bohr)"
-    availability: vdw_method is set to d3_0 or d3_bj
+    availability: vdw_method is set to d3_0, d3_bj, or d4
   - name: vdw_cn_thr_unit
     category: vdW correction
     type: String
diff --git a/source/source_hamilt/module_vdw/test/CMakeLists.txt b/source/source_hamilt/module_vdw/test/CMakeLists.txt
index 7e5c825128b..4b61f7f3000 100644
--- a/source/source_hamilt/module_vdw/test/CMakeLists.txt
+++ b/source/source_hamilt/module_vdw/test/CMakeLists.txt
@@ -10,3 +10,8 @@ AddTest(
   LIBS parameter ${math_libs} base device vdw 
   SOURCES vdw_test.cpp 
 )
+
+if(ENABLE_DFTD4)
+  target_compile_definitions(MODULE_HAMILT_vdwTest PRIVATE __DFTD4)
+  target_link_libraries(MODULE_HAMILT_vdwTest dftd4::dftd4)
+endif()
diff --git a/source/source_hamilt/module_vdw/test/vdw_test.cpp b/source/source_hamilt/module_vdw/test/vdw_test.cpp
index 5d66cb53a5c..72f76197d63 100644
--- a/source/source_hamilt/module_vdw/test/vdw_test.cpp
+++ b/source/source_hamilt/module_vdw/test/vdw_test.cpp
@@ -10,6 +10,9 @@
 #include "source_hamilt/module_vdw/vdwd3_parameters.h"
 #include "source_hamilt/module_vdw/vdwd2.h"
 #include "source_hamilt/module_vdw/vdwd3.h"
+#ifdef __DFTD4
+#include "source_hamilt/module_vdw/vdwd4.h"
+#endif
 #include "source_hamilt/module_vdw/vdw.h"
 #undef private
 
@@ -604,6 +607,73 @@ TEST_F(vdwd3abcTest, D3bjGetStress)
     EXPECT_NEAR(stress.e33, -3.4278442125590892e-05,1e-12);
 }
 
+#ifdef __DFTD4
+
+class vdwd4Test: public testing::Test
+{
+    protected:
+    UnitCell ucell;
+    Input_para input;
+
+    void SetUp(){
+        stru_ structure{std::vector<double>{0.5,0.5,0.0,0.5,0.0,0.5,0.0,0.5,0.5},
+                        std::vector<atomtype_>{atomtype_{"Si",
+                                                         std::vector<std::vector<double>>{
+                                                             {0., 0., 0.},
+                                                             {0.3, 0.25, 0.25}
+                                                         }}}};
+        construct_ucell(structure,ucell);
+
+        input.vdw_method = "d4";
+        input.vdw_d4_xc = "pbe";
+        input.vdw_cutoff_type = "radius";
+        input.vdw_radius_unit = "Bohr";
+        input.vdw_cutoff_radius = "60";
+        input.vdw_cn_thr_unit = "Bohr";
+        input.vdw_cn_thr = 30;
+    }
+
+    void TearDown(){
+        ClearUcell(ucell);
+    }
+};
+
+TEST_F(vdwd4Test, D4GetEnergy)
+{
+    auto vdw_solver = vdw::make_vdw(ucell, input);
+    double ene = vdw_solver->get_energy();
+    EXPECT_NEAR(ene, -0.04998837990336073, 1E-10);
+}
+
+TEST_F(vdwd4Test, D4GetForce)
+{
+    auto vdw_solver = vdw::make_vdw(ucell, input);
+    std::vector<ModuleBase::Vector3<double>> force = vdw_solver->get_force();
+    EXPECT_NEAR(force[0].x, -0.0023357259921368717, 1e-12);
+    EXPECT_NEAR(force[0].y, 0.0, 1e-12);
+    EXPECT_NEAR(force[0].z, 0.0, 1e-12);
+    EXPECT_NEAR(force[1].x, 0.0023357259921368730, 1e-12);
+    EXPECT_NEAR(force[1].y, 0.0, 1e-12);
+    EXPECT_NEAR(force[1].z, 0.0, 1e-12);
+}
+
+TEST_F(vdwd4Test, D4GetStress)
+{
+    auto vdw_solver = vdw::make_vdw(ucell, input);
+    ModuleBase::Matrix3 stress = vdw_solver->get_stress();
+    EXPECT_NEAR(stress.e11, 0.00015830384474877792, 1e-12);
+    EXPECT_NEAR(stress.e12, 0.0, 1e-12);
+    EXPECT_NEAR(stress.e13, 0.0, 1e-12);
+    EXPECT_NEAR(stress.e21, 0.0, 1e-12);
+    EXPECT_NEAR(stress.e22, 0.00016694998515968720, 1e-12);
+    EXPECT_NEAR(stress.e23, -1.5500973166318808e-05, 1e-12);
+    EXPECT_NEAR(stress.e31, 0.0, 1e-12);
+    EXPECT_NEAR(stress.e32, -1.5500973166318808e-05, 1e-12);
+    EXPECT_NEAR(stress.e33, 0.00016694998515968726, 1e-12);
+}
+
+#endif // __DFTD4
+
 int main(int argc, char **argv)
 {
     MPI_Init(&argc, &argv);

From 8c36a907410ffb77b28ff672a6f1196a98575185 Mon Sep 17 00:00:00 2001
From: Growl <uwsy1059@qq.com>
Date: Tue, 26 May 2026 10:16:18 +0800
Subject: [PATCH 3/7] Install dftd4 from toolchain in GitHub test

---
 .github/workflows/test.yml | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index a21d54b2a82..3c71534aa81 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -7,6 +7,10 @@ concurrency:
   group: ${{ github.workflow }}-${{ github.ref }}
   cancel-in-progress: true
 
+defaults:
+  run:
+    shell: bash
+
 jobs:
   test:
     name: Test
@@ -32,12 +36,20 @@ jobs:
       - name: Install CI tools
         run: |
           sudo apt-get update
-          sudo apt-get install -y ccache ca-certificates python-is-python3 python3-pip
+          sudo apt-get install -y gfortran ccache ca-certificates python-is-python3 python3-pip
           sudo pip install clang-format clang-tidy
 
+      - name: Install dftd4 from toolchain
+        run: |
+          cd toolchain
+          ./install_abacus_toolchain_new.sh --with-dftd4=install --dry-run -j8
+          ./scripts/stage4/install_dftd4.sh
+          cd ..
+
       - name: Configure
         run: |
-          cmake -B build -DBUILD_TESTING=ON -DENABLE_MLALGO=ON -DENABLE_LIBXC=ON -DENABLE_LIBRI=ON -DENABLE_GOOGLEBENCH=ON -DENABLE_RAPIDJSON=ON  -DCMAKE_EXPORT_COMPILE_COMMANDS=1 -DENABLE_FLOAT_FFTW=ON
+          source toolchain/install/setup
+          cmake -B build -DBUILD_TESTING=ON -DENABLE_MLALGO=ON -DENABLE_LIBXC=ON -DENABLE_LIBRI=ON -DENABLE_GOOGLEBENCH=ON -DENABLE_RAPIDJSON=ON  -DCMAKE_EXPORT_COMPILE_COMMANDS=1 -DENABLE_FLOAT_FFTW=ON -DENABLE_DFTD4=ON
 
 # Temporarily removed because no one maintains this now.
 # And it will break the CI test workflow.

From f93dc60fce95d5644044885dcbeed14802a13ebd Mon Sep 17 00:00:00 2001
From: Growl <uwsy1059@qq.com>
Date: Tue, 26 May 2026 10:31:26 +0800
Subject: [PATCH 4/7] Fix stress calculation

---
 source/source_hamilt/module_vdw/vdwd4.cpp | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/source/source_hamilt/module_vdw/vdwd4.cpp b/source/source_hamilt/module_vdw/vdwd4.cpp
index 7e3516c2026..03cfbfeb447 100644
--- a/source/source_hamilt/module_vdw/vdwd4.cpp
+++ b/source/source_hamilt/module_vdw/vdwd4.cpp
@@ -273,10 +273,13 @@ void Vdwd4::cal_stress()
     if (!has_stress_cache_)
     {
         double energy_ha = 0.0;
+        std::vector<double> gradient(3 * ucell_.nat, 0.0);
         std::array<double, 9> sigma;
         sigma.fill(0.0);
 
-        compute(energy_ha, nullptr, &sigma);
+        // DFT-D4 may require a valid gradient buffer when sigma is requested.
+        compute(energy_ha, &gradient, &sigma);
+        set_force_from_gradient(gradient);
         set_stress_from_sigma(sigma);
     }
 

From 61611376f69ceb2ad805f80623fda325bcdf005a Mon Sep 17 00:00:00 2001
From: Growl <uwsy1059@qq.com>
Date: Tue, 26 May 2026 10:57:40 +0800
Subject: [PATCH 5/7] Add regtest

---
 tests/03_NAO_multik/CASES_CPU.txt             |  1 +
 tests/03_NAO_multik/relax_cell_vdw4/INPUT     | 23 +++++++++++++++
 tests/03_NAO_multik/relax_cell_vdw4/KPT       |  4 +++
 tests/03_NAO_multik/relax_cell_vdw4/README    |  1 +
 tests/03_NAO_multik/relax_cell_vdw4/STRU      | 29 +++++++++++++++++++
 .../03_NAO_multik/relax_cell_vdw4/result.ref  |  5 ++++
 6 files changed, 63 insertions(+)
 create mode 100644 tests/03_NAO_multik/relax_cell_vdw4/INPUT
 create mode 100644 tests/03_NAO_multik/relax_cell_vdw4/KPT
 create mode 100644 tests/03_NAO_multik/relax_cell_vdw4/README
 create mode 100644 tests/03_NAO_multik/relax_cell_vdw4/STRU
 create mode 100644 tests/03_NAO_multik/relax_cell_vdw4/result.ref

diff --git a/tests/03_NAO_multik/CASES_CPU.txt b/tests/03_NAO_multik/CASES_CPU.txt
index c861b3dfe5a..43623591a26 100644
--- a/tests/03_NAO_multik/CASES_CPU.txt
+++ b/tests/03_NAO_multik/CASES_CPU.txt
@@ -51,6 +51,7 @@ relax_cell
 relax_cell_vdw2
 relax_cell_vdw3
 relax_cell_vdw3bj
+relax_cell_vdw4
 md_nvt
 md_msst
 md_out_wf
diff --git a/tests/03_NAO_multik/relax_cell_vdw4/INPUT b/tests/03_NAO_multik/relax_cell_vdw4/INPUT
new file mode 100644
index 00000000000..24acb742c38
--- /dev/null
+++ b/tests/03_NAO_multik/relax_cell_vdw4/INPUT
@@ -0,0 +1,23 @@
+INPUT_PARAMETERS 
+
+ecutwfc           10
+scf_thr           1e-06
+scf_nmax          400
+basis_type        lcao
+ks_solver         scalapack_gvx
+smearing_method   gaussian
+relax_nmax        1
+smearing_sigma    0.02
+mixing_type       broyden
+mixing_beta       0.7
+vdw_method        d4
+fixed_axes        c
+force_thr_ev      0.005
+stress_thr        0.1
+suffix            autotest
+calculation       cell-relax
+cal_force         1
+cal_stress        1
+pseudo_dir        ../../PP_ORB
+orbital_dir       ../../PP_ORB
+dft_functional    pbe
diff --git a/tests/03_NAO_multik/relax_cell_vdw4/KPT b/tests/03_NAO_multik/relax_cell_vdw4/KPT
new file mode 100644
index 00000000000..6c7f2228767
--- /dev/null
+++ b/tests/03_NAO_multik/relax_cell_vdw4/KPT
@@ -0,0 +1,4 @@
+K_POINTS
+0
+Monkhorst-Pack
+2 2 1 0 0 0
diff --git a/tests/03_NAO_multik/relax_cell_vdw4/README b/tests/03_NAO_multik/relax_cell_vdw4/README
new file mode 100644
index 00000000000..9126b6d767b
--- /dev/null
+++ b/tests/03_NAO_multik/relax_cell_vdw4/README
@@ -0,0 +1 @@
+test for SnTe system , VDWd4 effect on, LCAO base
diff --git a/tests/03_NAO_multik/relax_cell_vdw4/STRU b/tests/03_NAO_multik/relax_cell_vdw4/STRU
new file mode 100644
index 00000000000..accb0624488
--- /dev/null
+++ b/tests/03_NAO_multik/relax_cell_vdw4/STRU
@@ -0,0 +1,29 @@
+ATOMIC_SPECIES
+Sn  118.710  Sn_ONCV_PBE-1.0.upf
+Te  127.603  Te_ONCV_PBE-1.0.upf
+
+NUMERICAL_ORBITAL
+Sn_pbe_9.0au_100Ry_2s2p2d
+Te_pbe_9.0au_100Ry_2s2p2d
+
+LATTICE_CONSTANT
+1.89035917
+
+LATTICE_VECTORS
+4.646 0 0 #latvec1
+0 4.561 0 #latvec2
+0 0 16.32 #latvec3
+
+ATOMIC_POSITIONS
+Direct
+
+Sn #label
+0 #magnetism
+1 #number of atoms
+0.5493  0.5000000000000000  0.2 1 1 1
+
+Te #label
+0 #magnetism
+1 #number of atoms
+0.0000000000000000  0.0000000000000000  0.2 1 1 1
+
diff --git a/tests/03_NAO_multik/relax_cell_vdw4/result.ref b/tests/03_NAO_multik/relax_cell_vdw4/result.ref
new file mode 100644
index 00000000000..2dffcc4b937
--- /dev/null
+++ b/tests/03_NAO_multik/relax_cell_vdw4/result.ref
@@ -0,0 +1,5 @@
+etotref -4136.3779559087188318
+etotperatomref -2068.1889779544
+totalforceref 1.457000
+totalstressref 16.894570
+totaltimeref 1.50

From 57e6d04161530b38c9635eed247307fd378a11f8 Mon Sep 17 00:00:00 2001
From: Growl <uwsy1059@qq.com>
Date: Thu, 28 May 2026 01:12:17 +0800
Subject: [PATCH 6/7] Add D4S model

---
 docs/advanced/input_files/input-main.md       |  9 +++++
 docs/advanced/install.md                      |  3 +-
 docs/parameters.yaml                          |  8 ++++
 .../module_vdw/test/vdw_test.cpp              | 38 +++++++++++++++++++
 source/source_hamilt/module_vdw/vdwd4.cpp     | 26 +++++++++----
 source/source_hamilt/module_vdw/vdwd4.h       |  1 +
 .../module_parameter/input_parameter.h        |  1 +
 .../read_input_item_model.cpp                 | 22 +++++++++++
 tests/03_NAO_multik/CASES_CPU.txt             |  3 +-
 tests/03_NAO_multik/relax_cell_vdw4/README    |  1 -
 .../INPUT                                     |  1 +
 .../KPT                                       |  0
 tests/03_NAO_multik/relax_cell_vdw4_d4/README |  1 +
 .../STRU                                      |  0
 .../result.ref                                |  0
 tests/03_NAO_multik/relax_cell_vdw4_d4s/INPUT | 24 ++++++++++++
 tests/03_NAO_multik/relax_cell_vdw4_d4s/KPT   |  4 ++
 .../03_NAO_multik/relax_cell_vdw4_d4s/README  |  1 +
 tests/03_NAO_multik/relax_cell_vdw4_d4s/STRU  | 29 ++++++++++++++
 .../relax_cell_vdw4_d4s/result.ref            |  5 +++
 20 files changed, 167 insertions(+), 10 deletions(-)
 delete mode 100644 tests/03_NAO_multik/relax_cell_vdw4/README
 rename tests/03_NAO_multik/{relax_cell_vdw4 => relax_cell_vdw4_d4}/INPUT (96%)
 rename tests/03_NAO_multik/{relax_cell_vdw4 => relax_cell_vdw4_d4}/KPT (100%)
 create mode 100644 tests/03_NAO_multik/relax_cell_vdw4_d4/README
 rename tests/03_NAO_multik/{relax_cell_vdw4 => relax_cell_vdw4_d4}/STRU (100%)
 rename tests/03_NAO_multik/{relax_cell_vdw4 => relax_cell_vdw4_d4}/result.ref (100%)
 create mode 100644 tests/03_NAO_multik/relax_cell_vdw4_d4s/INPUT
 create mode 100644 tests/03_NAO_multik/relax_cell_vdw4_d4s/KPT
 create mode 100644 tests/03_NAO_multik/relax_cell_vdw4_d4s/README
 create mode 100644 tests/03_NAO_multik/relax_cell_vdw4_d4s/STRU
 create mode 100644 tests/03_NAO_multik/relax_cell_vdw4_d4s/result.ref

diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
index 740b4f56a96..a91c6dd7530 100644
--- a/docs/advanced/input_files/input-main.md
+++ b/docs/advanced/input_files/input-main.md
@@ -396,6 +396,8 @@
     - [sc\_scf\_thr](#sc_scf_thr)
   - [vdW correction](#vdw-correction)
     - [vdw\_method](#vdw_method)
+    - [vdw\_d4\_xc](#vdw_d4_xc)
+    - [vdw\_d4\_model](#vdw_d4_model)
     - [vdw\_s6](#vdw_s6)
     - [vdw\_s8](#vdw_s8)
     - [vdw\_a1](#vdw_a1)
@@ -3655,6 +3657,13 @@
 - **Description**: Functional name passed to the DFT-D4 library to load its internal damping parameters. If set to default, ABACUS infers the functional name from dft_functional or pseudopotential metadata.
 - **Default**: default
 
+### vdw_d4_model
+
+- **Type**: String
+- **Availability**: *vdw_method is set to d4*
+- **Description**: DFT-D4 dispersion model used by the external DFT-D4 library. Available options are d4 for the standard D4 model and d4s for the smooth D4S model.
+- **Default**: d4
+
 ### vdw_s6
 
 - **Type**: String
diff --git a/docs/advanced/install.md b/docs/advanced/install.md
index 654e853e204..520d1c386fc 100644
--- a/docs/advanced/install.md
+++ b/docs/advanced/install.md
@@ -91,9 +91,10 @@ In the input file, enable DFT-D4 with:
 ```text
 vdw_method d4
 vdw_d4_xc pbe
+vdw_d4_model d4    # or d4s for the smooth D4S model
 ```
 
-If `vdw_d4_xc` is set to `default`, ABACUS will infer the functional name from `dft_functional` or pseudopotential metadata and pass it to the DFT-D4 library.
+If `vdw_d4_xc` is set to `default`, ABACUS will infer the functional name from `dft_functional` or pseudopotential metadata and pass it to the DFT-D4 library. The `vdw_d4_model` keyword selects the dispersion model inside the DFT-D4 library; the default is `d4`, while `d4s` enables the smooth D4S model.
 
 ## Build Unit Tests
 
diff --git a/docs/parameters.yaml b/docs/parameters.yaml
index bb9d365b5ee..d8287e067c9 100644
--- a/docs/parameters.yaml
+++ b/docs/parameters.yaml
@@ -3755,6 +3755,14 @@ parameters:
     default_value: default
     unit: ""
     availability: vdw_method is set to d4
+  - name: vdw_d4_model
+    category: vdW correction
+    type: String
+    description: |
+      DFT-D4 dispersion model used by the external DFT-D4 library. Available options are d4 for the standard D4 model and d4s for the smooth D4S model.
+    default_value: d4
+    unit: ""
+    availability: vdw_method is set to d4
   - name: vdw_s6
     category: vdW correction
     type: String
diff --git a/source/source_hamilt/module_vdw/test/vdw_test.cpp b/source/source_hamilt/module_vdw/test/vdw_test.cpp
index 72f76197d63..faf62ba0139 100644
--- a/source/source_hamilt/module_vdw/test/vdw_test.cpp
+++ b/source/source_hamilt/module_vdw/test/vdw_test.cpp
@@ -626,6 +626,7 @@ class vdwd4Test: public testing::Test
 
         input.vdw_method = "d4";
         input.vdw_d4_xc = "pbe";
+        input.vdw_d4_model = "d4";
         input.vdw_cutoff_type = "radius";
         input.vdw_radius_unit = "Bohr";
         input.vdw_cutoff_radius = "60";
@@ -672,6 +673,43 @@ TEST_F(vdwd4Test, D4GetStress)
     EXPECT_NEAR(stress.e33, 0.00016694998515968726, 1e-12);
 }
 
+TEST_F(vdwd4Test, D4SGetEnergy)
+{
+    input.vdw_d4_model = "d4s";
+    auto vdw_solver = vdw::make_vdw(ucell, input);
+    double ene = vdw_solver->get_energy();
+    EXPECT_NEAR(ene, -0.05638517144755526, 1E-10);
+}
+
+TEST_F(vdwd4Test, D4SGetForce)
+{
+    input.vdw_d4_model = "d4s";
+    auto vdw_solver = vdw::make_vdw(ucell, input);
+    std::vector<ModuleBase::Vector3<double>> force = vdw_solver->get_force();
+    EXPECT_NEAR(force[0].x, -0.005448661796788402, 1e-12);
+    EXPECT_NEAR(force[0].y, 0.0, 1e-12);
+    EXPECT_NEAR(force[0].z, 0.0, 1e-12);
+    EXPECT_NEAR(force[1].x, 0.005448661796788397, 1e-12);
+    EXPECT_NEAR(force[1].y, 0.0, 1e-12);
+    EXPECT_NEAR(force[1].z, 0.0, 1e-12);
+}
+
+TEST_F(vdwd4Test, D4SGetStress)
+{
+    input.vdw_d4_model = "d4s";
+    auto vdw_solver = vdw::make_vdw(ucell, input);
+    ModuleBase::Matrix3 stress = vdw_solver->get_stress();
+    EXPECT_NEAR(stress.e11, 0.00013831119855416262, 1e-12);
+    EXPECT_NEAR(stress.e12, 0.0, 1e-12);
+    EXPECT_NEAR(stress.e13, 0.0, 1e-12);
+    EXPECT_NEAR(stress.e21, 0.0, 1e-12);
+    EXPECT_NEAR(stress.e22, 0.00015770515797834415, 1e-12);
+    EXPECT_NEAR(stress.e23, -3.862972112000666e-05, 1e-12);
+    EXPECT_NEAR(stress.e31, 0.0, 1e-12);
+    EXPECT_NEAR(stress.e32, -3.862972112000666e-05, 1e-12);
+    EXPECT_NEAR(stress.e33, 0.00015770515797834423, 1e-12);
+}
+
 #endif // __DFTD4
 
 int main(int argc, char **argv)
diff --git a/source/source_hamilt/module_vdw/vdwd4.cpp b/source/source_hamilt/module_vdw/vdwd4.cpp
index 03cfbfeb447..5eca0c63d32 100644
--- a/source/source_hamilt/module_vdw/vdwd4.cpp
+++ b/source/source_hamilt/module_vdw/vdwd4.cpp
@@ -10,7 +10,6 @@
 #include <cctype>
 #include <cstddef>
 #include <cstring>
-#include <iostream>
 #include <string>
 #include <vector>
 
@@ -24,12 +23,12 @@ namespace vdw
 namespace
 {
 
-std::string normalize_xc_name(std::string xc)
+std::string to_lower(std::string value)
 {
-    std::transform(xc.begin(), xc.end(), xc.begin(), [](unsigned char c) {
+    std::transform(value.begin(), value.end(), value.begin(), [](unsigned char c) {
         return static_cast<char>(std::tolower(c));
     });
-    return xc;
+    return value;
 }
 
 void check_dftd4_error(dftd4_error error, const std::string& where)
@@ -82,7 +81,7 @@ double cutoff_to_bohr(const std::string& value, const std::string& unit)
 } // namespace
 
 Vdwd4::Vdwd4(const UnitCell& unit_in, const std::string& xc_name, const Input_para& input)
-    : Vdw(unit_in), xc_name_(normalize_xc_name(xc_name))
+    : Vdw(unit_in), xc_name_(to_lower(xc_name)), model_name_(to_lower(input.vdw_d4_model))
 {
     cutoff_disp2_ = cutoff_to_bohr(input.vdw_cutoff_radius, input.vdw_radius_unit);
     cutoff_disp3_ = std::min(40.0, cutoff_disp2_);
@@ -169,8 +168,21 @@ void Vdwd4::compute(double& energy_ha,
                                               periodic.data());
     check_dftd4_error(error, "dftd4_new_structure");
 
-    dftd4_model model = dftd4_new_d4_model(error, mol);
-    check_dftd4_error(error, "dftd4_new_d4_model");
+    dftd4_model model = nullptr;
+    if (model_name_ == "d4")
+    {
+        model = dftd4_new_d4_model(error, mol);
+        check_dftd4_error(error, "dftd4_new_d4_model");
+    }
+    else if (model_name_ == "d4s")
+    {
+        model = dftd4_new_d4s_model(error, mol);
+        check_dftd4_error(error, "dftd4_new_d4s_model");
+    }
+    else
+    {
+        ModuleBase::WARNING_QUIT("Vdwd4::compute", "Unsupported DFT-D4 model: " + model_name_);
+    }
 
     dftd4_set_model_realspace_cutoff(error, model, cutoff_disp2_, cutoff_disp3_, cutoff_cn_);
     check_dftd4_error(error, "dftd4_set_model_realspace_cutoff");
diff --git a/source/source_hamilt/module_vdw/vdwd4.h b/source/source_hamilt/module_vdw/vdwd4.h
index 3772c8646d2..0cb46c2ca8c 100644
--- a/source/source_hamilt/module_vdw/vdwd4.h
+++ b/source/source_hamilt/module_vdw/vdwd4.h
@@ -19,6 +19,7 @@ class Vdwd4 : public Vdw
 
   private:
     std::string xc_name_;
+    std::string model_name_;
     double cutoff_disp2_ = 0.0; // Bohr, two-body dispersion cutoff
     double cutoff_disp3_ = 0.0; // Bohr, three-body ATM cutoff
     double cutoff_cn_ = 0.0;    // Bohr, coordination-number cutoff
diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h
index d953b542da5..6fac46679ec 100644
--- a/source/source_io/module_parameter/input_parameter.h
+++ b/source/source_io/module_parameter/input_parameter.h
@@ -527,6 +527,7 @@ struct Input_para
     double vdw_cn_thr = 40.0;                               ///< radius cutoff for cn
     std::string vdw_cn_thr_unit = "Bohr";                   ///< unit of cn_thr, Bohr or Angstrom
     std::string vdw_d4_xc = "default";                     ///< functional name passed to DFT-D4
+    std::string vdw_d4_model = "d4";                       ///< DFT-D4 dispersion model (d4 or d4s)
     ModuleBase::Vector3<int> vdw_cutoff_period = {3, 3, 3}; ///< periods of periodic structure
 
     // ==============   #Parameters (15.exx) ====================
diff --git a/source/source_io/module_parameter/read_input_item_model.cpp b/source/source_io/module_parameter/read_input_item_model.cpp
index a3bd0b0a9f8..2b29a9b5ec7 100644
--- a/source/source_io/module_parameter/read_input_item_model.cpp
+++ b/source/source_io/module_parameter/read_input_item_model.cpp
@@ -283,6 +283,28 @@ If set to default, ABACUS infers the functional name from dft_functional or pseu
         read_sync_string(input.vdw_d4_xc);
         this->add_item(item);
     }
+    {
+        Input_Item item("vdw_d4_model");
+        item.annotation = "DFT-D4 dispersion model";
+        item.category = "vdW correction";
+        item.type = "String";
+        item.description = R"(DFT-D4 dispersion model used by the external DFT-D4 library.
+Available options are:
+* d4: standard D4 model
+* d4s: smooth D4S model)";
+        item.default_value = "d4";
+        item.unit = "";
+        item.availability = "vdw_method is set to d4";
+        read_sync_string(input.vdw_d4_model);
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            if (para.input.vdw_d4_model != "d4" && para.input.vdw_d4_model != "d4s"
+                && para.input.vdw_d4_model != "D4" && para.input.vdw_d4_model != "D4S")
+            {
+                ModuleBase::WARNING_QUIT("ReadInput", "vdw_d4_model must be d4 or d4s");
+            }
+        };
+        this->add_item(item);
+    }
     {
         Input_Item item("vdw_s6");
         item.annotation = "scale parameter of d2/d3_0/d3_bj";
diff --git a/tests/03_NAO_multik/CASES_CPU.txt b/tests/03_NAO_multik/CASES_CPU.txt
index 43623591a26..beb5b151631 100644
--- a/tests/03_NAO_multik/CASES_CPU.txt
+++ b/tests/03_NAO_multik/CASES_CPU.txt
@@ -51,7 +51,8 @@ relax_cell
 relax_cell_vdw2
 relax_cell_vdw3
 relax_cell_vdw3bj
-relax_cell_vdw4
+relax_cell_vdw4_d4
+relax_cell_vdw4_d4s
 md_nvt
 md_msst
 md_out_wf
diff --git a/tests/03_NAO_multik/relax_cell_vdw4/README b/tests/03_NAO_multik/relax_cell_vdw4/README
deleted file mode 100644
index 9126b6d767b..00000000000
--- a/tests/03_NAO_multik/relax_cell_vdw4/README
+++ /dev/null
@@ -1 +0,0 @@
-test for SnTe system , VDWd4 effect on, LCAO base
diff --git a/tests/03_NAO_multik/relax_cell_vdw4/INPUT b/tests/03_NAO_multik/relax_cell_vdw4_d4/INPUT
similarity index 96%
rename from tests/03_NAO_multik/relax_cell_vdw4/INPUT
rename to tests/03_NAO_multik/relax_cell_vdw4_d4/INPUT
index 24acb742c38..f07f6fcbf69 100644
--- a/tests/03_NAO_multik/relax_cell_vdw4/INPUT
+++ b/tests/03_NAO_multik/relax_cell_vdw4_d4/INPUT
@@ -11,6 +11,7 @@ smearing_sigma    0.02
 mixing_type       broyden
 mixing_beta       0.7
 vdw_method        d4
+vdw_d4_model      d4
 fixed_axes        c
 force_thr_ev      0.005
 stress_thr        0.1
diff --git a/tests/03_NAO_multik/relax_cell_vdw4/KPT b/tests/03_NAO_multik/relax_cell_vdw4_d4/KPT
similarity index 100%
rename from tests/03_NAO_multik/relax_cell_vdw4/KPT
rename to tests/03_NAO_multik/relax_cell_vdw4_d4/KPT
diff --git a/tests/03_NAO_multik/relax_cell_vdw4_d4/README b/tests/03_NAO_multik/relax_cell_vdw4_d4/README
new file mode 100644
index 00000000000..b7b9258d5e1
--- /dev/null
+++ b/tests/03_NAO_multik/relax_cell_vdw4_d4/README
@@ -0,0 +1 @@
+test for SnTe system , VDWd4 effect on with model D4, LCAO base
diff --git a/tests/03_NAO_multik/relax_cell_vdw4/STRU b/tests/03_NAO_multik/relax_cell_vdw4_d4/STRU
similarity index 100%
rename from tests/03_NAO_multik/relax_cell_vdw4/STRU
rename to tests/03_NAO_multik/relax_cell_vdw4_d4/STRU
diff --git a/tests/03_NAO_multik/relax_cell_vdw4/result.ref b/tests/03_NAO_multik/relax_cell_vdw4_d4/result.ref
similarity index 100%
rename from tests/03_NAO_multik/relax_cell_vdw4/result.ref
rename to tests/03_NAO_multik/relax_cell_vdw4_d4/result.ref
diff --git a/tests/03_NAO_multik/relax_cell_vdw4_d4s/INPUT b/tests/03_NAO_multik/relax_cell_vdw4_d4s/INPUT
new file mode 100644
index 00000000000..30409c02a5b
--- /dev/null
+++ b/tests/03_NAO_multik/relax_cell_vdw4_d4s/INPUT
@@ -0,0 +1,24 @@
+INPUT_PARAMETERS 
+
+ecutwfc           10
+scf_thr           1e-06
+scf_nmax          400
+basis_type        lcao
+ks_solver         scalapack_gvx
+smearing_method   gaussian
+relax_nmax        1
+smearing_sigma    0.02
+mixing_type       broyden
+mixing_beta       0.7
+vdw_method        d4
+vdw_d4_model      d4s
+fixed_axes        c
+force_thr_ev      0.005
+stress_thr        0.1
+suffix            autotest
+calculation       cell-relax
+cal_force         1
+cal_stress        1
+pseudo_dir        ../../PP_ORB
+orbital_dir       ../../PP_ORB
+dft_functional    pbe
diff --git a/tests/03_NAO_multik/relax_cell_vdw4_d4s/KPT b/tests/03_NAO_multik/relax_cell_vdw4_d4s/KPT
new file mode 100644
index 00000000000..6c7f2228767
--- /dev/null
+++ b/tests/03_NAO_multik/relax_cell_vdw4_d4s/KPT
@@ -0,0 +1,4 @@
+K_POINTS
+0
+Monkhorst-Pack
+2 2 1 0 0 0
diff --git a/tests/03_NAO_multik/relax_cell_vdw4_d4s/README b/tests/03_NAO_multik/relax_cell_vdw4_d4s/README
new file mode 100644
index 00000000000..6b0b7c0b71b
--- /dev/null
+++ b/tests/03_NAO_multik/relax_cell_vdw4_d4s/README
@@ -0,0 +1 @@
+test for SnTe system , VDWd4 effect on with model D4S, LCAO base
diff --git a/tests/03_NAO_multik/relax_cell_vdw4_d4s/STRU b/tests/03_NAO_multik/relax_cell_vdw4_d4s/STRU
new file mode 100644
index 00000000000..accb0624488
--- /dev/null
+++ b/tests/03_NAO_multik/relax_cell_vdw4_d4s/STRU
@@ -0,0 +1,29 @@
+ATOMIC_SPECIES
+Sn  118.710  Sn_ONCV_PBE-1.0.upf
+Te  127.603  Te_ONCV_PBE-1.0.upf
+
+NUMERICAL_ORBITAL
+Sn_pbe_9.0au_100Ry_2s2p2d
+Te_pbe_9.0au_100Ry_2s2p2d
+
+LATTICE_CONSTANT
+1.89035917
+
+LATTICE_VECTORS
+4.646 0 0 #latvec1
+0 4.561 0 #latvec2
+0 0 16.32 #latvec3
+
+ATOMIC_POSITIONS
+Direct
+
+Sn #label
+0 #magnetism
+1 #number of atoms
+0.5493  0.5000000000000000  0.2 1 1 1
+
+Te #label
+0 #magnetism
+1 #number of atoms
+0.0000000000000000  0.0000000000000000  0.2 1 1 1
+
diff --git a/tests/03_NAO_multik/relax_cell_vdw4_d4s/result.ref b/tests/03_NAO_multik/relax_cell_vdw4_d4s/result.ref
new file mode 100644
index 00000000000..35bbc88f527
--- /dev/null
+++ b/tests/03_NAO_multik/relax_cell_vdw4_d4s/result.ref
@@ -0,0 +1,5 @@
+etotref -4136.3692216877407191
+etotperatomref -2068.1846108439
+totalforceref 1.428828
+totalstressref 16.950691
+totaltimeref 0.41

From 7d1148f1bae8154368de22e4b742f10a94777223 Mon Sep 17 00:00:00 2001
From: Growl <uwsy1059@qq.com>
Date: Sat, 30 May 2026 11:36:01 +0800
Subject: [PATCH 7/7] Add citations

---
 docs/CITATIONS.md | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/docs/CITATIONS.md b/docs/CITATIONS.md
index a25cb3e596c..43176829a1d 100644
--- a/docs/CITATIONS.md
+++ b/docs/CITATIONS.md
@@ -47,3 +47,11 @@ The following references are required to be cited when using ABACUS. Specificall
     Sun, Liang, and Mohan Chen. "Machine learning based nonlocal kinetic energy density functional for simple metals and alloys." Physical Review B 109.11 (2024): 115135.
 
     Sun, Liang, and Mohan Chen. "Multi-channel machine learning based nonlocal kinetic energy density functional for semiconductors." Electronic Structure 6.4 (2024): 045006.
+
+- **If DFT-D4 dispersion correction is used:**
+
+    Eike Caldeweyher, Sebastian Ehlert, Andreas Hansen, Hagen Neugebauer, Sebastian Spicher, Christoph Bannwarth, and Stefan Grimme. "A generally applicable atomic-charge dependent London dispersion correction." The Journal of Chemical Physics 150, 154122 (2019). DOI: 10.1063/1.5090222.
+
+    Eike Caldeweyher, Jan-Michael Mewes, Sebastian Ehlert, and Stefan Grimme. "Extension and evaluation of the D4 London-dispersion model for periodic systems." Physical Chemistry Chemical Physics 22, 8499-8512 (2020). DOI: 10.1039/D0CP00502A.
+
+    Nikolay V. Tkachenko, Linus B. Dittmer, Rebecca Tomann, and Martin Head-Gordon. "Smooth D4S model." The Journal of Physical Chemistry Letters 15, 10629-10637 (2024). DOI: 10.1021/acs.jpclett.4c02653.