vdW-D4: Switch to smooth cutoff

docs/advanced/input_files/input-main.md

@@ -431,6 +431,8 @@
     - [vdw\_cutoff\_type](#vdw_cutoff_type)
     - [vdw\_cutoff\_radius](#vdw_cutoff_radius)
     - [vdw\_radius\_unit](#vdw_radius_unit)
+    - [vdw\_cutoff\_smooth\_width\_2b](#vdw_cutoff_smooth_width_2b)
+    - [vdw\_cutoff\_smooth\_width\_3b](#vdw_cutoff_smooth_width_3b)
     - [vdw\_cutoff\_period](#vdw_cutoff_period)
     - [vdw\_cn\_thr](#vdw_cn_thr)
     - [vdw\_cn\_thr\_unit](#vdw_cn_thr_unit)
@@ -4062,6 +4064,24 @@
   - Bohr
 - **Default**: Bohr
 
+### vdw_cutoff_smooth_width_2b
+
+- **Type**: Real
+- **Availability**: *[`vdw_method`](#vdw_method)==d4*
+- **Description**: Width of the smooth switching region for the two-body pairwise dispersion real-space cutoff.
+  A value of zero disables smoothing for the two-body contribution.
+- **Default**: 0.05
+- **Unit**: Bohr
+
+### vdw_cutoff_smooth_width_3b
+
+- **Type**: Real
+- **Availability**: *[`vdw_method`](#vdw_method)==d4*
+- **Description**: Width of the smooth switching region for the three-body Axilrod-Teller-Muto (ATM) dispersion real-space cutoff.
+  A value of zero disables smoothing for the three-body contribution.
+- **Default**: 0.05
+- **Unit**: Bohr
+
 ### vdw_cutoff_period
 
 - **Type**: Integer Integer Integer

docs/parameters.yaml

@@ -4312,6 +4312,24 @@ parameters:
     default_value: Bohr
     unit: ""
     availability: vdw_cutoff_type==radius
+  - name: vdw_cutoff_smooth_width_2b
+    category: vdW correction
+    type: Real
+    description: |
+      Width of the smooth switching region for the two-body pairwise dispersion real-space cutoff.
+      A value of zero disables smoothing for the two-body contribution.
+    default_value: "0.05"
+    unit: Bohr
+    availability: vdw_method==d4
+  - name: vdw_cutoff_smooth_width_3b
+    category: vdW correction
+    type: Real
+    description: |
+      Width of the smooth switching region for the three-body Axilrod-Teller-Muto (ATM) dispersion real-space cutoff.
+      A value of zero disables smoothing for the three-body contribution.
+    default_value: "0.05"
+    unit: Bohr
+    availability: vdw_method==d4
   - name: vdw_cutoff_period
     category: vdW correction
     type: Integer Integer Integer

source/source_hamilt/module_vdw/test/vdw_test.cpp

@@ -690,7 +690,7 @@ TEST_F(vdwd4Test, D4GetEnergy)
     auto vdw_solver = vdw::make_vdw(ucell, input);
     const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(false, false));
     const double ene = result.energy;
-    EXPECT_NEAR(ene, -0.04998837990336073, 1E-10);
+    EXPECT_NEAR(ene, -0.049988405722573105, 1E-10);
 }
 
 TEST_F(vdwd4Test, D4GetEnergyForChargedSystem)
@@ -700,21 +700,21 @@ TEST_F(vdwd4Test, D4GetEnergyForChargedSystem)
     auto vdw_solver = vdw::make_vdw(ucell, input);
     const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(false, false));
     const double ene = result.energy;
-    EXPECT_NEAR(ene, -0.04359451765256733, 1E-10);
+    EXPECT_NEAR(ene, -0.04359454509118302, 1E-10);
 }
 
 TEST_F(vdwd4Test, D4GetForce)
 {
     auto vdw_solver = vdw::make_vdw(ucell, input);
     const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, false));
-    EXPECT_NEAR(result.energy, -0.04998837990336073, 1E-10);
+    EXPECT_NEAR(result.energy, -0.049988405722573105, 1E-10);
     ASSERT_TRUE(result.has_force);
     EXPECT_FALSE(result.has_stress);
     const std::vector<ModuleBase::Vector3<double>>& force = result.force;
-    EXPECT_NEAR(force[0].x, -0.0023357259921368717, 1e-12);
+    EXPECT_NEAR(force[0].x, -0.002339156758188389, 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].x, 0.0023391567581883886, 1e-12);
     EXPECT_NEAR(force[1].y, 0.0, 1e-12);
     EXPECT_NEAR(force[1].z, 0.0, 1e-12);
 }
@@ -723,19 +723,19 @@ TEST_F(vdwd4Test, D4GetStress)
 {
     auto vdw_solver = vdw::make_vdw(ucell, input);
     const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, true));
-    EXPECT_NEAR(result.energy, -0.04998837990336073, 1E-10);
+    EXPECT_NEAR(result.energy, -0.049988405722573105, 1E-10);
     ASSERT_TRUE(result.has_force);
     ASSERT_TRUE(result.has_stress);
     const ModuleBase::Matrix3& stress = result.stress;
-    EXPECT_NEAR(stress.e11, 0.00015830384474877792, 1e-12);
+    EXPECT_NEAR(stress.e11, 0.0001583939298091549, 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.e22, 0.00016697881796423088, 1e-12);
+    EXPECT_NEAR(stress.e23, -1.527806618822572e-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);
+    EXPECT_NEAR(stress.e32, -1.527806618822572e-05, 1e-12);
+    EXPECT_NEAR(stress.e33, 0.0001669788179642309, 1e-12);
 }
 
 TEST_F(vdwd4Test, D4SGetEnergy)
@@ -744,22 +744,22 @@ TEST_F(vdwd4Test, D4SGetEnergy)
     auto vdw_solver = vdw::make_vdw(ucell, input);
     const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(false, false));
     const double ene = result.energy;
-    EXPECT_NEAR(ene, -0.05638517144755526, 1E-10);
+    EXPECT_NEAR(ene, -0.05638520357171156, 1E-10);
 }
 
 TEST_F(vdwd4Test, D4SGetForce)
 {
     input.vdw_d4_model = "d4s";
     auto vdw_solver = vdw::make_vdw(ucell, input);
     const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, false));
-    EXPECT_NEAR(result.energy, -0.05638517144755526, 1E-10);
+    EXPECT_NEAR(result.energy, -0.05638520357171156, 1E-10);
     ASSERT_TRUE(result.has_force);
     EXPECT_FALSE(result.has_stress);
     const std::vector<ModuleBase::Vector3<double>>& force = result.force;
-    EXPECT_NEAR(force[0].x, -0.005448661796788402, 1e-12);
+    EXPECT_NEAR(force[0].x, -0.005452776236973487, 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].x, 0.005452776236973491, 1e-12);
     EXPECT_NEAR(force[1].y, 0.0, 1e-12);
     EXPECT_NEAR(force[1].z, 0.0, 1e-12);
 }
@@ -769,19 +769,19 @@ TEST_F(vdwd4Test, D4SGetStress)
     input.vdw_d4_model = "d4s";
     auto vdw_solver = vdw::make_vdw(ucell, input);
     const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, true));
-    EXPECT_NEAR(result.energy, -0.05638517144755526, 1E-10);
+    EXPECT_NEAR(result.energy, -0.05638520357171156, 1E-10);
     ASSERT_TRUE(result.has_force);
     ASSERT_TRUE(result.has_stress);
     const ModuleBase::Matrix3& stress = result.stress;
-    EXPECT_NEAR(stress.e11, 0.00013831119855416262, 1e-12);
+    EXPECT_NEAR(stress.e11, 0.0001384186027460731, 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.e22, 0.00015772616666498505, 1e-12);
+    EXPECT_NEAR(stress.e23, -3.836792114896563e-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);
+    EXPECT_NEAR(stress.e32, -3.836792114896563e-05, 1e-12);
+    EXPECT_NEAR(stress.e33, 0.0001577261666649851, 1e-12);
 }
 
 #endif // __DFTD4

source/source_hamilt/module_vdw/vdwd4.cpp

@@ -86,6 +86,19 @@ Vdwd4::Vdwd4(const UnitCell& unit_in, const std::string& xc_name, const Input_pa
     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);
+    smooth_width_2b_ = input.vdw_cutoff_smooth_width_2b;
+    smooth_width_3b_ = input.vdw_cutoff_smooth_width_3b;
+
+    if (smooth_width_2b_ < 0.0 || smooth_width_2b_ > cutoff_disp2_)
+    {
+        ModuleBase::WARNING_QUIT("Vdwd4::Vdwd4",
+                                 "vdw_cutoff_smooth_width_2b must satisfy 0 <= width <= two-body cutoff");
+    }
+    if (smooth_width_3b_ < 0.0 || smooth_width_3b_ > cutoff_disp3_)
+    {
+        ModuleBase::WARNING_QUIT("Vdwd4::Vdwd4",
+                                 "vdw_cutoff_smooth_width_3b must satisfy 0 <= width <= three-body cutoff");
+    }
 
     double valence_charge = 0.0;
     for (int it = 0; it < ucell_.ntype; ++it)
@@ -191,8 +204,14 @@ void Vdwd4::compute(double& energy_ha,
         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");
+    dftd4_set_model_realspace_cutoff_smooth(error,
+                                            model,
+                                            cutoff_disp2_,
+                                            cutoff_disp3_,
+                                            cutoff_cn_,
+                                            smooth_width_2b_,
+                                            smooth_width_3b_);
+    check_dftd4_error(error, "dftd4_set_model_realspace_cutoff_smooth");
 
     std::vector<char> method(xc_name_.begin(), xc_name_.end());
     method.push_back('\0');

source/source_hamilt/module_vdw/vdwd4.h

@@ -23,6 +23,8 @@ class Vdwd4 : public Vdw
     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
+    double smooth_width_2b_ = 0.0; // Bohr, two-body cutoff smoothing width
+    double smooth_width_3b_ = 0.0; // Bohr, three-body cutoff smoothing width
     double total_charge_ = 0.0; // e, total system charge (sum zv*na - nelec)
 
     void evaluate_impl(const VdwRequest& request, VdwResult& result) override;

source/source_io/module_parameter/input_parameter.h

@@ -569,6 +569,8 @@ struct Input_para
                                                             ///< structure, radius or period
     std::string vdw_cutoff_radius = "default";              ///< radius cutoff for periodic structure
     std::string vdw_radius_unit = "Bohr";                   ///< unit of radius cutoff for periodic structure
+    double vdw_cutoff_smooth_width_2b = 0.05;               ///< smooth cutoff width for two-body dispersion, Bohr
+    double vdw_cutoff_smooth_width_3b = 0.05;               ///< smooth cutoff width for three-body dispersion, Bohr
     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

source/source_io/module_parameter/read_inp_model.cpp

@@ -572,6 +572,44 @@ Namely, each line contains the element name and the corresponding parameter.)";
         };
         this->add_item(item);
     }
+    {
+        Input_Item item("vdw_cutoff_smooth_width_2b");
+        item.annotation = "smooth cutoff width for two-body dispersion";
+        item.category = "vdW correction";
+        item.type = "Real";
+        item.description = R"(Width of the smooth switching region for the two-body pairwise dispersion real-space cutoff.
+A value of zero disables smoothing for the two-body contribution.)";
+        item.default_value = "0.05";
+        item.unit = "Bohr";
+        item.set_availability("vdw_method==d4");
+        read_sync_double(input.vdw_cutoff_smooth_width_2b);
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            if (para.input.vdw_cutoff_smooth_width_2b < 0.0)
+            {
+                ModuleBase::WARNING_QUIT("ReadInput", "vdw_cutoff_smooth_width_2b must be non-negative");
+            }
+        };
+        this->add_item(item);
+    }
+    {
+        Input_Item item("vdw_cutoff_smooth_width_3b");
+        item.annotation = "smooth cutoff width for three-body dispersion";
+        item.category = "vdW correction";
+        item.type = "Real";
+        item.description = R"(Width of the smooth switching region for the three-body Axilrod-Teller-Muto (ATM) dispersion real-space cutoff.
+A value of zero disables smoothing for the three-body contribution.)";
+        item.default_value = "0.05";
+        item.unit = "Bohr";
+        item.set_availability("vdw_method==d4");
+        read_sync_double(input.vdw_cutoff_smooth_width_3b);
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            if (para.input.vdw_cutoff_smooth_width_3b < 0.0)
+            {
+                ModuleBase::WARNING_QUIT("ReadInput", "vdw_cutoff_smooth_width_3b must be non-negative");
+            }
+        };
+        this->add_item(item);
+    }
     {
         Input_Item item("vdw_cutoff_period");
         item.annotation = "periods of periodic structure";

source/source_io/test/read_input_ptest.cpp

@@ -274,6 +274,8 @@ TEST_F(InputParaTest, ParaRead)
     EXPECT_FALSE(param.inp.vdw_abc);
     EXPECT_EQ(std::stod(param.inp.vdw_cutoff_radius), 56.6918);
     EXPECT_EQ(param.inp.vdw_radius_unit, "Bohr");
+    EXPECT_DOUBLE_EQ(param.inp.vdw_cutoff_smooth_width_2b, 0.05);
+    EXPECT_DOUBLE_EQ(param.inp.vdw_cutoff_smooth_width_3b, 0.05);
     EXPECT_DOUBLE_EQ(param.inp.vdw_cn_thr, 40.0);
     EXPECT_EQ(param.inp.vdw_cn_thr_unit, "Bohr");
     EXPECT_EQ(param.inp.vdw_C6_file, "default");

source/source_io/test_serial/read_input_item_test.cpp

@@ -1440,6 +1440,22 @@ TEST_F(InputTest, Item_test2)
         output = testing::internal::GetCapturedStdout();
         EXPECT_THAT(output, testing::HasSubstr("NOTICE"));
     }
+    { // vdw_cutoff_smooth_width_2b
+        auto it = find_label("vdw_cutoff_smooth_width_2b", readinput.input_lists);
+        param.input.vdw_cutoff_smooth_width_2b = -1.0;
+        testing::internal::CaptureStdout();
+        EXPECT_EXIT(it->second.check_value(it->second, param), ::testing::ExitedWithCode(1), "");
+        output = testing::internal::GetCapturedStdout();
+        EXPECT_THAT(output, testing::HasSubstr("NOTICE"));
+    }
+    { // vdw_cutoff_smooth_width_3b
+        auto it = find_label("vdw_cutoff_smooth_width_3b", readinput.input_lists);
+        param.input.vdw_cutoff_smooth_width_3b = -1.0;
+        testing::internal::CaptureStdout();
+        EXPECT_EXIT(it->second.check_value(it->second, param), ::testing::ExitedWithCode(1), "");
+        output = testing::internal::GetCapturedStdout();
+        EXPECT_THAT(output, testing::HasSubstr("NOTICE"));
+    }
     { // vdw_cn_thr
         auto it = find_label("vdw_cn_thr", readinput.input_lists);
         param.input.vdw_cn_thr = -1;