Feature: support the default as the value of dft_functional when initialize vdw

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<int> rep_vdw_ref = {input.vdw_cutoff_period.x, input.vdw_cutoff_period.y, input.vdw_cutoff_period.z};
 

source/module_hamilt_general/module_vdw/vdw.cpp

@@ -1,7 +1,41 @@
+#include <algorithm>
+#include <cassert>
 
 #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<std::string> &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<std::string> 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<Vdw> make_vdw(const UnitCell &ucell,
     }
     else if (input.vdw_method == "d3_0" || input.vdw_method == "d3_bj")
     {
+        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::unique_ptr<Vdwd3> vdw_ptr = make_unique<Vdwd3>(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")

source/module_hamilt_general/module_vdw/vdw.h

@@ -11,9 +11,10 @@
 namespace vdw
 {
 
-template<typename T, typename... Args>
-std::unique_ptr<T> make_unique(Args &&... args) {
-    return std::unique_ptr<T>(new T(std::forward<Args>(args)...));
+    template<typename T, typename... Args>
+    std::unique_ptr<T> make_unique(Args &&... args) {
+        return std::unique_ptr<T>(new T(std::forward<Args>(args)...));
+
 }
 
 class Vdw

source/module_hamilt_general/module_vdw/vdwd3_autoset_xcparam.cpp

@@ -70,6 +70,7 @@
 
 // DFT-D3(BJ)
 const std::map<std::string, std::vector<double>> 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<std::string, std::vector<double>> bj = {
 };
 // DFT-D3(0)
 const std::map<std::string, std::vector<double>> 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<std::string, std::vector<double>> zero = {
 };
 // DFT-D3M(BJ): not implemented for beta
 const std::map<std::string, std::vector<double>> 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<std::string, std::vector<double>> bjm = {
 };
 // DFT-D3M(0): not implemented for beta
 const std::map<std::string, std::vector<double>> 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<std::string, std::vector<double>> zerom = {
 };
 // DFT-D3(OptimizedPower)
 const std::map<std::string, std::vector<double>> 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<std::string, std::vector<double>> op = {
     {"ms2h", {1.0, 0.65, 0.65, 1.69464, 4.75, 4.75, 1.0, 14.0, 0.0}},
 };
 
-
+std::vector<double> _search_impl(const std::string& xc,
+                                 const std::map<std::string, std::vector<double>>& dict)
+{
+    if (dict.find(xc) != dict.end())
+    {
+        return dict.at(xc);
+    }
+    else
+    {
+        return std::vector<double>();
+    }
+}
 // 's6', 'rs6', 'a1', 's8', 'rs8', 'a2', 's9', 'alp', 'bet'
 /**
  * @brief Get the dftd3 params object. 
@@ -368,76 +385,49 @@ const std::map<std::string, std::vector<double>> 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<double>& param)
+             const std::string& method, 
+             std::vector<double>& param)
 {
     const std::string xc_lowercase = FmtCore::lower(xc);
     const std::vector<std::string> 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<std::string, std::vector<double>> 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);
     }
 }
 

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<std::vector<double>>(max_elem_, std::vector<double>(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);

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_; }