Merge branch 'extended-dftd3' of https://github.com/kirk0830/abacus-develop into extended-dftd3

Author: kirk0830

docs/advanced/input_files/input-main.md

@@ -2923,27 +2923,27 @@ These variables are used to control vdW-corrected related parameters.
 
 - **Type**: Real
 - **Availability**: `vdw_method` is set to `d2`, `d3_0`, or `d3_bj`
-- **Description**: This scale factor is used to optimize the interaction energy deviations in van der Waals (vdW) corrected calculations. The recommended values of this parameter are dependent on the chosen vdW correction method and the DFT functional being used. For DFT-D2, the recommended values are 0.75 (PBE), 1.2 (BLYP), 1.05 (B-P86), 1.0 (TPSS), and 1.05 (B3LYP). If not set, will use values of PBE functional. For DFT-D3, recommended values with different DFT functionals can be found on the [here](https://www.chemiebn.uni-bonn.de/pctc/mulliken-center/software/dft-d3/dft-d3). If not set, will search in ABACUS built-in dataset based on the `dft_functional` keywords. User set value will overwrite the searched value.
+- **Description**: This scale factor is used to optimize the interaction energy deviations in van der Waals (vdW) corrected calculations. The recommended values of this parameter are dependent on the chosen vdW correction method and the DFT functional being used. For DFT-D2, the recommended values are 0.75 (PBE), 1.2 (BLYP), 1.05 (B-P86), 1.0 (TPSS), and 1.05 (B3LYP). If not set, will use values of PBE functional. For DFT-D3, recommended values with different DFT functionals can be found on the [here](https://github.com/dftd3/simple-dftd3/blob/main/assets/parameters.toml). If not set, will search in ABACUS built-in dataset based on the `dft_functional` keywords. User set value will overwrite the searched value.
 - **Default**:
   - 0.75: if `vdw_method` is set to `d2`
 
 ### vdw_s8
 
 - **Type**: Real
 - **Availability**: `vdw_method` is set to `d3_0` or `d3_bj`
-- **Description**: This scale factor is relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://www.chemiebn.uni-bonn.de/pctc/mulliken-center/software/dft-d3/dft-d3). If not set, will search in ABACUS built-in dataset based on the `dft_functional` keywords. User set value will overwrite the searched value.
+- **Description**: This scale factor is relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://github.com/dftd3/simple-dftd3/blob/main/assets/parameters.toml). If not set, will search in ABACUS built-in dataset based on the `dft_functional` keywords. User set value will overwrite the searched value.
 
 ### vdw_a1
 
 - **Type**: Real
 - **Availability**: `vdw_method` is set to `d3_0` or `d3_bj`
-- **Description**: This damping function parameter is relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://www.chemiebn.uni-bonn.de/pctc/mulliken-center/software/dft-d3/dft-d3). If not set, will search in ABACUS built-in dataset based on the `dft_functional` keywords. User set value will overwrite the searched value.
+- **Description**: This damping function parameter is relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://github.com/dftd3/simple-dftd3/blob/main/assets/parameters.toml). If not set, will search in ABACUS built-in dataset based on the `dft_functional` keywords. User set value will overwrite the searched value.
 
 ### vdw_a2
 
 - **Type**: Real
 - **Availability**: `vdw_method` is set to `d3_0` or `d3_bj`
-- **Description**: This damping function parameter is only relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://www.chemiebn.uni-bonn.de/pctc/mulliken-center/software/dft-d3/dft-d3). If not set, will search in ABACUS built-in dataset based on the `dft_functional` keywords. User set value will overwrite the searched value.
+- **Description**: This damping function parameter is only relevant for D3(0) and D3(BJ) van der Waals (vdW) correction methods. The recommended values of this parameter with different DFT functionals can be found on the [webpage](https://github.com/dftd3/simple-dftd3/blob/main/assets/parameters.toml). If not set, will search in ABACUS built-in dataset based on the `dft_functional` keywords. User set value will overwrite the searched value.
 
 ### vdw_d
 

source/module_esolver/esolver_ks_pw.cpp

@@ -210,7 +210,7 @@ void ESolver_KS_PW<T, Device>::before_scf(const int istep)
     //----------------------------------------------------------
     // about vdw, jiyy add vdwd3 and linpz add vdwd2
     //----------------------------------------------------------
-    auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp);
+    auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp, &(GlobalV::ofs_running));
     if (vdw_solver != nullptr)
     {
         this->pelec->f_en.evdw = vdw_solver->get_energy();

source/module_esolver/lcao_before_scf.cpp

@@ -198,7 +198,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(const int istep)
     //----------------------------------------------------------
     // about vdw, jiyy add vdwd3 and linpz add vdwd2
     //----------------------------------------------------------
-    auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp);
+    auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp, &(GlobalV::ofs_running));
     if (vdw_solver != nullptr)
     {
         this->pelec->f_en.evdw = vdw_solver->get_energy();

source/module_hamilt_general/module_vdw/dftd3_xc_name.h

@@ -344,11 +344,14 @@ namespace DFTD3 {
     void _xcname_libxc_xplusc(const std::string& xcpattern, std::string& xname)
     {
         std::vector<std::string> xc_words = FmtCore::split(xcpattern, "+");
+        std::for_each(xc_words.begin(), xc_words.end(), [](std::string& s) {
+             s = (FmtCore::startswith(s, "XC_")? s: "XC_" + s); }); // add XC_ if not present
         assert(xc_words.size() == 2);
 
         std::vector<std::string> words = FmtCore::split(xc_words[0], "_");
-        const std::string key = (words[2] == "X")? xcpattern: xc_words[1] + "+" + xc_words[0];
-        
+        const std::string key = (words[2] == "X")? 
+            xc_words[0] + "+" + xc_words[1]: xc_words[1] + "+" + xc_words[0];
+
         if (xcname_libxc_xplusc_.find(key) != xcname_libxc_xplusc_.end()) {
             xname = xcname_libxc_xplusc_.at(key);
         } else {
@@ -359,9 +362,11 @@ namespace DFTD3 {
 
     void _xcname_libxc_xc(const std::string& xcpattern, std::string& xname)
     {
-        std::vector<std::string> words = FmtCore::split(xcpattern, "_");
-        if (xcname_libxc_xc_.find(xcpattern) != xcname_libxc_xc_.end()) {
-            xname = xcname_libxc_xc_.at(xcpattern);
+        // add XC_ if not present
+        const std::string key = FmtCore::startswith(xcpattern, "XC_")? xcpattern: "XC_" + xcpattern;
+
+        if (xcname_libxc_xc_.find(key) != xcname_libxc_xc_.end()) {
+            xname = xcname_libxc_xc_.at(key);
         } else {
             ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::xcname_libxc_xc",
                                      "XC's LibXC-notation on `" + xcpattern + "` not recognized");
@@ -380,9 +385,9 @@ namespace DFTD3 {
     std::string _xcname(const std::string& xcpattern)
     {
         std::string xcname = xcpattern;
-        const std::regex pattern("XC_(LDA|GGA|MGGA|HYB|HYB_LDA|HYB_GGA|HYB_MGGA)_(X|C|XC|K)_(.*)");
+        const std::regex pattern("(LDA|GGA|MGGA|HYB|HYB_LDA|HYB_GGA|HYB_MGGA)_(X|C|XC|K)_(.*)");
         // as long as there is piece in xcpattern that can match, we can search for the corresponding name
-        if (std::regex_match(xcpattern, pattern)) {
+        if (std::regex_search(xcpattern, pattern)) {
             _xcname_libxc(xcpattern, xcname);
         }
         return xcname;

source/module_hamilt_general/module_vdw/dftd3_xc_param.h

@@ -421,7 +421,7 @@ namespace DFTD3 {
                                          "XC (`" + xc + "`)'s DFT-D3M(0) parameters not found");
             }
         }
-        else // zero
+        else if (method == "zero")
         {
             if (zero.find(xc_lowercase) != zero.end())
             {
@@ -433,6 +433,11 @@ namespace DFTD3 {
                                          "XC (`" + xc + "`)'s DFT-D3(0) parameters not found");
             }
         }
+        else // should not reach here
+        {
+            ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::DFTD3::_search",
+                                     "Unknown DFT-D3 method: " + method);
+        }
     }
 
     std::string _lowercase(const std::string& s)
@@ -474,7 +479,8 @@ namespace DFTD3 {
                       double& s6,
                       double& s8,
                       double& a1,
-                      double& a2)
+                      double& a2,
+                      std::ofstream* plog = nullptr)
     {
         const std::map<std::string, std::string> param_map = {
             {"d3_bj", "bj"}, {"d3_0", "zero"}, {"d3_bjm", "bjm"}, {"d3_0m", "zerom"},
@@ -498,12 +504,16 @@ namespace DFTD3 {
             s8 = (s8_in == "default") ? param[3] : std::stod(s8_in);
             a1 = (a1_in == "default") ? param[2] : std::stod(a1_in);
             a2 = (a2_in == "default") ? param[5] : std::stod(a2_in);
-            // param = {s6, s8, a1, a2};
-            // FmtTable vdwd3tab({"Parameters", "Original", "Autoset"}, 4, {"%10s", "%10s", "%10.4f"});
-            // const std::vector<std::string> items = {"s6", "s8", "a1", "a2"};
-            // vdwd3tab << items << flag << param;
-            // std::cout << "DFT-D3 Dispersion correction parameters autoset\n" << vdwd3tab.str()
-            //           << std::flush;
+            if (plog != nullptr) // logging the autoset
+            {
+                param = {s6, s8, a1, a2};
+                FmtTable vdwd3tab({"Parameters", "Original", "Autoset"}, 4, {"%10s", "%10s", "%10.4f"});
+                const std::vector<std::string> items = {"s6", "s8", "a1", "a2"};
+                vdwd3tab << items << flag << param;
+                (*plog) << "\nDFT-D3 Dispersion correction parameters autoset\n" << vdwd3tab.str()
+                        << "XC functional: " << xc_in << std::endl;
+            }
+
         }
     }
 }

source/module_hamilt_general/module_vdw/test/dftd3_xc_test.cpp

@@ -4,25 +4,31 @@
 
 TEST(DFTD3XCTest, SearchXcnameLibXCXplusC)
 {
-    std::string xcpattern = "XC_GGA_X_PBE+XC_GGA_C_OP_PBE";
     std::string xname;
-    DFTD3::_xcname_libxc_xplusc(xcpattern, xname);
+    DFTD3::_xcname_libxc_xplusc("XC_GGA_X_PBE+XC_GGA_C_OP_PBE", xname);
+    EXPECT_EQ(xname, "pbeop");
+    // then test the case with out prefix XC_
+    DFTD3::_xcname_libxc_xplusc("GGA_X_PBE+GGA_C_OP_PBE", xname);
     EXPECT_EQ(xname, "pbeop");
 }
 
 TEST(DFTD3XCTest, SearchXcnameLibXCXC)
 {
-    std::string xcpattern = "XC_LDA_XC_TETER93";
     std::string xname;
-    DFTD3::_xcname_libxc_xc(xcpattern, xname);
+    DFTD3::_xcname_libxc_xc("XC_LDA_XC_TETER93", xname);
+    EXPECT_EQ(xname, "teter93");
+    // then test the case with out prefix XC_
+    DFTD3::_xcname_libxc_xc("LDA_XC_TETER93", xname);
     EXPECT_EQ(xname, "teter93");
 }
 
 TEST(DFTD3XCTest, SearchXcnameLibXC)
 {
-    std::string xcpattern = "XC_GGA_X_PBE+XC_GGA_C_OP_PBE";
     std::string xname;
-    DFTD3::_xcname_libxc(xcpattern, xname);
+    DFTD3::_xcname_libxc("XC_GGA_X_PBE+XC_GGA_C_OP_PBE", xname);
+    EXPECT_EQ(xname, "pbeop");
+    // then test the case with out prefix XC_
+    DFTD3::_xcname_libxc("GGA_X_PBE+GGA_C_OP_PBE", xname);
     EXPECT_EQ(xname, "pbeop");
 }
 
@@ -45,6 +51,52 @@ TEST(DFTD3XCTest, SearchXcname)
     EXPECT_EQ(xname, "pbe");
 }
 
+TEST(DFTD3XCTest, SuccessfulSearch)
+{
+    std::string xc = "pbe";
+    std::string d3method = "d3_0";
+    std::string s6_in = "default";
+    std::string s8_in = "default";
+    std::string a1_in = "default";
+    std::string a2_in = "default";
+    double s6, s8, a1, a2;
+    DFTD3::dftd3_params(xc, d3method, s6_in, s8_in, a1_in, a2_in, s6, s8, a1, a2);
+    EXPECT_DOUBLE_EQ(s6, 1.0);
+    EXPECT_DOUBLE_EQ(s8, 0.722);
+    EXPECT_DOUBLE_EQ(a1, 1.217);
+    EXPECT_DOUBLE_EQ(a2, 1.0);
+
+    // a more complicated case: MGGA_X_SCAN+MGGA_C_SCAN
+    xc = "XC_MGGA_X_SCAN+XC_MGGA_C_SCAN";
+    DFTD3::dftd3_params(xc, d3method, s6_in, s8_in, a1_in, a2_in, s6, s8, a1, a2);
+    EXPECT_DOUBLE_EQ(s6, 1.0);
+    EXPECT_DOUBLE_EQ(s8, 0.0);
+    EXPECT_DOUBLE_EQ(a1, 1.324);
+    EXPECT_DOUBLE_EQ(a2, 1.0);
+
+    // user defines all parameters
+    s6_in = "1.1";
+    s8_in = "0.1";
+    a1_in = "1.325";
+    a2_in = "1.1";
+    DFTD3::dftd3_params(xc, d3method, s6_in, s8_in, a1_in, a2_in, s6, s8, a1, a2);
+    EXPECT_DOUBLE_EQ(s6, 1.1);
+    EXPECT_DOUBLE_EQ(s8, 0.1);
+    EXPECT_DOUBLE_EQ(a1, 1.325);
+    EXPECT_DOUBLE_EQ(a2, 1.1);
+
+    // user defines one parameter
+    s6_in = "1.1";
+    s8_in = "default";
+    a1_in = "default";
+    a2_in = "default";
+    DFTD3::dftd3_params(xc, d3method, s6_in, s8_in, a1_in, a2_in, s6, s8, a1, a2);
+    EXPECT_DOUBLE_EQ(s6, 1.1);
+    EXPECT_DOUBLE_EQ(s8, 0.0);
+    EXPECT_DOUBLE_EQ(a1, 1.324);
+    EXPECT_DOUBLE_EQ(a2, 1.0);
+}
+
 int main(int argc, char **argv)
 {
     ::testing::InitGoogleTest(&argc, argv);

source/module_hamilt_general/module_vdw/vdw.cpp

@@ -6,30 +6,35 @@
 namespace vdw
 {
 
-std::unique_ptr<Vdw> make_vdw(const UnitCell &ucell, const Input_para &input)
+std::unique_ptr<Vdw> make_vdw(const UnitCell &ucell, 
+                              const Input_para &input,
+                              std::ofstream* plog)
 {
-    if (ucell.nat < 2 && input.vdw_method != "none")
-    {
-        ModuleBase::WARNING("VDW", "Only one atom in this system, and will not do the calculation of VDW");
-        return nullptr;
-    }
-    else if (input.vdw_method == "d2")
+    // if (ucell.nat < 2 && input.vdw_method != "none")
+    // {
+    //     ModuleBase::WARNING("VDW", "Only one atom in this system, and will not do the calculation of VDW");
+    //     return nullptr;
+    // }
+    if (input.vdw_method == "d2")
     {
         std::unique_ptr<Vdwd2> vdw_ptr = make_unique<Vdwd2>(ucell);
-        vdw_ptr->parameter().initial_parameters(input);
+        vdw_ptr->parameter().initial_parameters(input, plog);
         vdw_ptr->parameter().initset(ucell);
         return vdw_ptr;
     }
     else if (input.vdw_method == "d3_0" || input.vdw_method == "d3_bj")
     {
         std::unique_ptr<Vdwd3> vdw_ptr = make_unique<Vdwd3>(ucell);
-        vdw_ptr->parameter().initial_parameters(input);
+        vdw_ptr->parameter().initial_parameters(input, plog);
         return vdw_ptr;
     }
-    else
+    else if (input.vdw_method != "none")
     {
+        ModuleBase::WARNING_QUIT("ModuleHamiltGeneral::ModuleVDW::make_vdw", 
+        "Unrecognized Van der Waals correction method: " + input.vdw_method);
         return nullptr;
     }
+    return nullptr; // "none" method
 }
 
 } // namespace vdw

source/module_hamilt_general/module_vdw/vdw.h

@@ -48,7 +48,17 @@ class Vdw
     virtual void cal_stress() = 0;
 };
 
-std::unique_ptr<Vdw> make_vdw(const UnitCell &ucell, const Input_para &input);
+/**
+ * @brief make vdw correction object
+ * 
+ * @param ucell UnitCell instance
+ * @param input Parameter instance
+ * @param plog optional, for logging the parameter setting process
+ * @return std::unique_ptr<Vdw> 
+ */
+std::unique_ptr<Vdw> make_vdw(const UnitCell &ucell, 
+                              const Input_para &input,
+                              std::ofstream* plog = nullptr);
 
 } // namespace vdw
 

source/module_hamilt_general/module_vdw/vdwd2_parameters.cpp

@@ -9,7 +9,7 @@
 namespace vdw
 {
 
-void Vdwd2Parameters::initial_parameters(const Input_para &input)
+void Vdwd2Parameters::initial_parameters(const Input_para &input, std::ofstream* plog)
 {
     scaling_ = std::stod(input.vdw_s6);
     damping_ = input.vdw_d;

source/module_hamilt_general/module_vdw/vdwd2_parameters.h

@@ -31,7 +31,15 @@ class Vdwd2Parameters : public VdwParameters
     void R0_input(const std::string &file, const std::string &unit);
 
     void initset(const UnitCell &ucell); // init sets of vdwd2 once this correction is called
-    void initial_parameters(const Input_para &input); // initial parameters of Vdwd2 with INPUT file
+
+    /**
+     * @brief initial parameters of Vdwd2 with INPUT file
+     * 
+     * @param input Parameter instance
+     * @param plog optional, for logging the parameter setting process (not implemented)
+     */
+    void initial_parameters(const Input_para &input, 
+                            std::ofstream* plog = nullptr); 
 
     inline const std::map<std::string, double> C6() const { return C6_; }
     inline const std::map<std::string, double> R0() const { return R0_; }