diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 422fcd3835..df92857bef 100644 --- a/docs/advanced/input_files/input-main.md +++ b/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 diff --git a/docs/parameters.yaml b/docs/parameters.yaml index ef311d26fc..157cd09bae 100644 --- a/docs/parameters.yaml +++ b/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 diff --git a/source/source_hamilt/module_vdw/test/vdw_test.cpp b/source/source_hamilt/module_vdw/test/vdw_test.cpp index 01aa74ec27..f08c92482f 100644 --- a/source/source_hamilt/module_vdw/test/vdw_test.cpp +++ b/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>& 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,7 +744,7 @@ 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) @@ -752,14 +752,14 @@ 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>& 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 diff --git a/source/source_hamilt/module_vdw/vdwd4.cpp b/source/source_hamilt/module_vdw/vdwd4.cpp index 8f66b68b2c..6acc5c76bb 100644 --- a/source/source_hamilt/module_vdw/vdwd4.cpp +++ b/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 method(xc_name_.begin(), xc_name_.end()); method.push_back('\0'); diff --git a/source/source_hamilt/module_vdw/vdwd4.h b/source/source_hamilt/module_vdw/vdwd4.h index f035db3d0b..994a5b9aac 100644 --- a/source/source_hamilt/module_vdw/vdwd4.h +++ b/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; diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 008414a9c3..c1a10113fe 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/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 diff --git a/source/source_io/module_parameter/read_inp_model.cpp b/source/source_io/module_parameter/read_inp_model.cpp index a0db9a9f4c..2ce78925dd 100644 --- a/source/source_io/module_parameter/read_inp_model.cpp +++ b/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"; diff --git a/source/source_io/test/read_input_ptest.cpp b/source/source_io/test/read_input_ptest.cpp index 5260238a3c..d5f81807b5 100644 --- a/source/source_io/test/read_input_ptest.cpp +++ b/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"); diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp index cb86a5d2d8..2dcbd8eba2 100644 --- a/source/source_io/test_serial/read_input_item_test.cpp +++ b/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;