Add optional DFT-D4 support (#7380)

* Add optional DFT-D4 support

* Docs and tests

* Install dftd4 from toolchain in GitHub test

* Fix stress calculation

* Add regtest

* Add D4S model

* Add citations

Author: Growl1234

.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.

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

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.

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)
@@ -3640,11 +3642,28 @@
   - 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_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
@@ -3737,7 +3756,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 +3778,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)
 

docs/advanced/install.md

@@ -77,6 +77,25 @@ 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
+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. 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
 
 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.

docs/parameters.yaml

@@ -3738,12 +3738,31 @@ 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_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
@@ -3852,7 +3871,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 +3897,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

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)

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()
+endif()

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

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,111 @@ 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_d4_model = "d4";
+        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);
+}
+
+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)
 {
     MPI_Init(&argc, &argv);