diff --git a/source/module_hamilt_general/module_vdw/test/vdw_test.cpp b/source/module_hamilt_general/module_vdw/test/vdw_test.cpp index 8b4bf62aae..ed5483b40c 100644 --- a/source/module_hamilt_general/module_vdw/test/vdw_test.cpp +++ b/source/module_hamilt_general/module_vdw/test/vdw_test.cpp @@ -377,7 +377,7 @@ class vdwd3Test: public testing::Test TEST_F(vdwd3Test, D30Default) { vdw::Vdwd3 vdwd3_test(ucell); - vdwd3_test.parameter().initial_parameters(input); + vdwd3_test.parameter().initial_parameters("pbe", input); EXPECT_EQ(vdwd3_test.parameter().s6(), 1.0); EXPECT_EQ(vdwd3_test.parameter().s18(), 0.7875); @@ -396,7 +396,8 @@ TEST_F(vdwd3Test, D30UnitA) input.vdw_cn_thr_unit = "A"; vdw::Vdwd3 vdwd3_test(ucell); - vdwd3_test.parameter().initial_parameters(input); + const std::string xc = "pbe"; + vdwd3_test.parameter().initial_parameters(xc, input); EXPECT_EQ(vdwd3_test.parameter().rthr2(), std::pow(95/ModuleBase::BOHR_TO_A, 2)); EXPECT_EQ(vdwd3_test.parameter().cn_thr2(), std::pow(40/ModuleBase::BOHR_TO_A, 2)); @@ -407,7 +408,8 @@ TEST_F(vdwd3Test, D30Period) input.vdw_cutoff_type = "period"; vdw::Vdwd3 vdwd3_test(ucell); - vdwd3_test.parameter().initial_parameters(input); + const std::string xc = "pbe"; + vdwd3_test.parameter().initial_parameters(xc, input); vdwd3_test.init(); std::vector rep_vdw_ref = {input.vdw_cutoff_period.x, input.vdw_cutoff_period.y, input.vdw_cutoff_period.z}; diff --git a/source/module_hamilt_general/module_vdw/vdw.cpp b/source/module_hamilt_general/module_vdw/vdw.cpp index b070d527ae..4842dbb0f2 100644 --- a/source/module_hamilt_general/module_vdw/vdw.cpp +++ b/source/module_hamilt_general/module_vdw/vdw.cpp @@ -1,7 +1,41 @@ +#include +#include #include "vdw.h" #include "vdwd2.h" #include "vdwd3.h" +#include "module_base/tool_quit.h" + +std::string parse_xcname(const std::string &xc_input, + const std::vector &xc_psp) +{ + if (xc_input != "default") + { + return xc_input; + } + + if (xc_psp.size() <= 0) + { + ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::parse_xcname", + "XC name automatic inference failed: no pseudopotential files are found"); + } + std::vector xc_psp_uniq = xc_psp; + std::sort(xc_psp_uniq.begin(), xc_psp_uniq.end()); + auto last = std::unique(xc_psp_uniq.begin(), xc_psp_uniq.end()); + xc_psp_uniq.erase(last, xc_psp_uniq.end()); + + if (xc_psp_uniq.size() > 1) + { + ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::parse_xcname", + "XC name automatic inference failed: inconsistency in XC names is found" + " in the pseudopotential files"); + } + const std::string xc = xc_psp_uniq[0]; + std::cout << " ***WARNING*** ModuleHamiltGeneral::ModuleVDW::parse_xcname: " + << "XC name is automatically inferred from pseudopotential as `" + << xc << "`" << std::endl; + return xc; +} namespace vdw { @@ -24,8 +58,13 @@ std::unique_ptr make_vdw(const UnitCell &ucell, } else if (input.vdw_method == "d3_0" || input.vdw_method == "d3_bj") { + std::vector xc_psp(ucell.ntype); + for (int it = 0; it < ucell.ntype; it++) + { + xc_psp[it] = ucell.atoms[it].ncpp.xc_func; + } std::unique_ptr vdw_ptr = make_unique(ucell); - vdw_ptr->parameter().initial_parameters(input, plog); + vdw_ptr->parameter().initial_parameters(parse_xcname(input.dft_functional, xc_psp), input, plog); return vdw_ptr; } else if (input.vdw_method != "none") diff --git a/source/module_hamilt_general/module_vdw/vdw.h b/source/module_hamilt_general/module_vdw/vdw.h index 257b80cd98..a4b494c8ad 100644 --- a/source/module_hamilt_general/module_vdw/vdw.h +++ b/source/module_hamilt_general/module_vdw/vdw.h @@ -11,9 +11,10 @@ namespace vdw { -template -std::unique_ptr make_unique(Args &&... args) { - return std::unique_ptr(new T(std::forward(args)...)); + template + std::unique_ptr make_unique(Args &&... args) { + return std::unique_ptr(new T(std::forward(args)...)); + } class Vdw diff --git a/source/module_hamilt_general/module_vdw/vdwd3_autoset_xcparam.cpp b/source/module_hamilt_general/module_vdw/vdwd3_autoset_xcparam.cpp index 4f2f077dac..be02f156a1 100644 --- a/source/module_hamilt_general/module_vdw/vdwd3_autoset_xcparam.cpp +++ b/source/module_hamilt_general/module_vdw/vdwd3_autoset_xcparam.cpp @@ -70,6 +70,7 @@ // DFT-D3(BJ) const std::map> bj = { + {"__default__", {1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 14.0, 0.0}}, {"bp", {1.0, 0.3946, 0.3946, 3.2822, 4.8516, 4.8516, 1.0, 14.0, 0.0}}, {"blyp", {1.0, 0.4298, 0.4298, 2.6996, 4.2359, 4.2359, 1.0, 14.0, 0.0}}, {"revpbe", {1.0, 0.5238, 0.5238, 2.355, 3.5016, 3.5016, 1.0, 14.0, 0.0}}, @@ -231,6 +232,7 @@ const std::map> bj = { }; // DFT-D3(0) const std::map> zero = { + {"__default__", {1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 14.0, 0.0}}, {"slaterdirac", {1.0, 0.999, 0.999, -1.957, 0.697, 0.697, 1.0, 14.0, 0.0}}, {"bp", {1.0, 1.139, 1.139, 1.683, 1.0, 1.0, 1.0, 14.0, 0.0}}, {"blyp", {1.0, 1.094, 1.094, 1.682, 1.0, 1.0, 1.0, 14.0, 0.0}}, @@ -318,6 +320,7 @@ const std::map> zero = { }; // DFT-D3M(BJ): not implemented for beta const std::map> bjm = { + {"__default__", {1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 14.0, 0.0}}, {"bp", {1.0, 0.82185, 0.82185, 3.140281, 2.728151, 2.728151, 1.0, 14.0, 0.0}}, {"blyp", {1.0, 0.448486, 0.448486, 1.875007, 3.610679, 3.610679, 1.0, 14.0, 0.0}}, {"b97_d", {1.0, 0.240184, 0.240184, 1.206988, 3.864426, 3.864426, 1.0, 14.0, 0.0}}, @@ -329,6 +332,7 @@ const std::map> bjm = { }; // DFT-D3M(0): not implemented for beta const std::map> zerom = { + {"__default__", {1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 14.0, 0.0}}, {"bp", {1.0, 1.23346, 1.23346, 1.945174, 1.0, 1.0, 1.0, 14.0, 0.0}}, {"blyp", {1.0, 1.279637, 1.279637, 1.841686, 1.0, 1.0, 1.0, 14.0, 0.01437}}, {"b97_d", {1.0, 1.151808, 1.151808, 1.020078, 1.0, 1.0, 1.0, 14.0, 0.035964}}, @@ -340,6 +344,8 @@ const std::map> zerom = { }; // DFT-D3(OptimizedPower) const std::map> op = { + // {'s6', 'rs6', 'a1', 's8', 'rs8', 'a2', 's9', 'alp', 'bet'} + {"__default__", {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 14.0, 0.0}}, {"blyp", {1.0, 0.425, 0.425, 1.31867, 3.5, 3.5, 1.0, 14.0, 2.0}}, {"revpbe", {1.0, 0.6, 0.6, 1.44765, 2.5, 2.5, 1.0, 14.0, 0.0}}, {"b97_d", {1.0, 0.6, 0.6, 1.46861, 2.5, 2.5, 1.0, 14.0, 0.0}}, @@ -356,7 +362,18 @@ const std::map> op = { {"ms2h", {1.0, 0.65, 0.65, 1.69464, 4.75, 4.75, 1.0, 14.0, 0.0}}, }; - +std::vector _search_impl(const std::string& xc, + const std::map>& dict) +{ + if (dict.find(xc) != dict.end()) + { + return dict.at(xc); + } + else + { + return std::vector(); + } +} // 's6', 'rs6', 'a1', 's8', 'rs8', 'a2', 's9', 'alp', 'bet' /** * @brief Get the dftd3 params object. @@ -368,76 +385,49 @@ const std::map> op = { * @param param the dftd3 parameters, ALL_KEYS = {'s6', 'rs6', 'a1', 's8', 'rs8', 'a2', 's9', 'alp', 'bet'} */ void _search(const std::string& xc, - const std::string& method, - std::vector& param) + const std::string& method, + std::vector& param) { const std::string xc_lowercase = FmtCore::lower(xc); const std::vector allowed_ = { "bj", "zero", "bjm", "zerom", "op" }; - assert(std::find(allowed_.begin(), allowed_.end(), method) != allowed_.end()); - if (method == "op") - { - if (op.find(xc_lowercase) != op.end()) - { - param = op.at(xc_lowercase); - } - else - { - ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::_search", - "XC (`" + xc + "`)'s DFT-D3(OP) parameters not found"); - } - } - else if (method == "bjm") + const int i = std::find(allowed_.begin(), allowed_.end(), method) - allowed_.begin(); + std::map> const * pdict = nullptr; + switch (i) { - if (bjm.find(xc_lowercase) != bjm.end()) - { - param = bjm.at(xc_lowercase); - } - else - { - ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::_search", - "XC (`" + xc + "`)'s DFT-D3M(BJ) parameters not found"); - } + case 0: + pdict = &bj; + break; + case 1: + pdict = &zero; + break; + case 2: + pdict = &bjm; + break; + case 3: + pdict = &zerom; + break; + case 4: + pdict = &op; + break; + default: + pdict = nullptr; + break; } - else if (method == "bj") + if (pdict == nullptr) { - if (bj.find(xc_lowercase) != bj.end()) - { - param = bj.at(xc_lowercase); - } - else - { - ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::_search", - "XC (`" + xc + "`)'s DFT-D3(BJ) parameters not found"); - } - } - else if (method == "zerom") - { - if (zerom.find(xc_lowercase) != zerom.end()) - { - param = zerom.at(xc_lowercase); - } - else - { - ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::_search", - "XC (`" + xc + "`)'s DFT-D3M(0) parameters not found"); - } - } - else if (method == "zero") - { - if (zero.find(xc_lowercase) != zero.end()) - { - param = zero.at(xc_lowercase); - } - else - { - ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::_search", - "XC (`" + xc + "`)'s DFT-D3(0) parameters not found"); - } + ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::_search", + "Unknown DFT-D3 method: " + method); } - else // should not reach here + param = _search_impl(xc_lowercase, *pdict); + if (param.empty()) { ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::_search", - "Unknown DFT-D3 method: " + method); + "XC (`" + xc + "`)'s DFT-D3(" + method + ") parameters not found"); + // is it meaningful to return a so-called default value? + std::cout << " ***WARNING*** " + << "XC (`" << xc << "`)'s DFT-D3(" << method << ") parameters not found, " + << "using default values. Please use at your own risk!" << std::endl; + param = _search_impl("__default__", *pdict); } } diff --git a/source/module_hamilt_general/module_vdw/vdwd3_parameters.cpp b/source/module_hamilt_general/module_vdw/vdwd3_parameters.cpp index d9e18e33b8..fc2736529b 100644 --- a/source/module_hamilt_general/module_vdw/vdwd3_parameters.cpp +++ b/source/module_hamilt_general/module_vdw/vdwd3_parameters.cpp @@ -10,7 +10,9 @@ namespace vdw { -void Vdwd3Parameters::initial_parameters(const Input_para &input, std::ofstream* plog) +void Vdwd3Parameters::initial_parameters(const std::string& xc, + const Input_para& input, + std::ofstream* plog) { // initialize the dftd3 parameters mxc_.resize(max_elem_, 1); @@ -23,7 +25,7 @@ void Vdwd3Parameters::initial_parameters(const Input_para &input, std::ofstream* 5, std::vector>(max_elem_, std::vector(max_elem_, 0.0))))); - _vdwd3_autoset_xcparam(input.dft_functional, input.vdw_method, + _vdwd3_autoset_xcparam(xc, input.vdw_method, input.vdw_s6, input.vdw_s8, input.vdw_a1, input.vdw_a2, s6_, s18_, rs6_, rs18_, /* rs6: a1, rs18: a2 */ plog); diff --git a/source/module_hamilt_general/module_vdw/vdwd3_parameters.h b/source/module_hamilt_general/module_vdw/vdwd3_parameters.h index 75c5801eb3..ceaadd6307 100644 --- a/source/module_hamilt_general/module_vdw/vdwd3_parameters.h +++ b/source/module_hamilt_general/module_vdw/vdwd3_parameters.h @@ -27,7 +27,8 @@ class Vdwd3Parameters : public VdwParameters * @param input Parameter instance * @param plog optional, for logging the parameter setting process */ - void initial_parameters(const Input_para &input, + void initial_parameters(const std::string& xc, + const Input_para& input, std::ofstream* plog = nullptr); // for logging the parameter autoset inline const std::string &version() const { return version_; }