|
11 | 11 |
|
12 | 12 | #include <xc.h> |
13 | 13 | #include <vector> |
| 14 | +#include <regex> |
| 15 | +#include <map> |
| 16 | +#include <algorithm> |
| 17 | +#include <cassert> |
| 18 | + |
14 | 19 | bool not_supported_xc_with_laplacian(const std::string& xc_func_in) |
15 | 20 | { |
16 | 21 | // see Pyscf: https://github.com/pyscf/pyscf/blob/master/pyscf/dft/libxc.py#L1062 |
@@ -55,118 +60,174 @@ bool not_supported_xc_with_nonlocal_vdw(const std::string& xc_func_in) |
55 | 60 | return false; |
56 | 61 | } |
57 | 62 |
|
58 | | -std::pair<int,std::vector<int>> XC_Functional_Libxc::set_xc_type_libxc(std::string xc_func_in) |
| 63 | +int xc_func_type_classifier(const std::string& xc_func, |
| 64 | + const std::map<std::string, int>& mymap = { |
| 65 | + {"LDA", 1}, |
| 66 | + {"GGA", 2}, |
| 67 | + {"MGGA", 3}, |
| 68 | + {"HYB_LDA", 4}, |
| 69 | + {"HYB_GGA", 4}, |
| 70 | + {"HYB_MGGA", 5} |
| 71 | + }) |
| 72 | +{ |
| 73 | + // the libxc standard functional pattern is like: |
| 74 | + // "(XC_)?(LDA|GGA|MGGA|HYB_GGA|HYB_MGGA|HYB_LDA)_(X|C|XC|K)(_(.*))?" |
| 75 | + std::regex pattern(R"((XC_)?(LDA|GGA|MGGA|HYB_GGA|HYB_MGGA|HYB_LDA)_(X|C|XC|K)(_(.*))?)"); |
| 76 | + std::smatch match; |
| 77 | + if (std::regex_match(xc_func, match, pattern)) { |
| 78 | + std::string type = match[2].str(); |
| 79 | + auto it = mymap.find(type); |
| 80 | + if (it != mymap.end()) { |
| 81 | + return it->second; |
| 82 | + } else { |
| 83 | + ModuleBase::WARNING_QUIT("XC_Functional_Libxc::xc_func_type_classifier", |
| 84 | + "Unrecognized functional type: " + type); |
| 85 | + } |
| 86 | + } else { |
| 87 | + ModuleBase::WARNING_QUIT("XC_Functional_Libxc::xc_func_type_classifier", |
| 88 | + "Unrecognized functional format: " + xc_func); |
| 89 | + } |
| 90 | +} |
| 91 | + |
| 92 | +std::pair<int,std::vector<int>> |
| 93 | +XC_Functional_Libxc::set_xc_type_libxc(const std::string& xc_func_in) |
59 | 94 | { |
60 | | - // determine the type (lda/gga/mgga) |
| 95 | + // check if the functional involves Laplacian of rho |
61 | 96 | if (not_supported_xc_with_laplacian(xc_func_in)) |
62 | 97 | { |
63 | 98 | ModuleBase::WARNING_QUIT("XC_Functional::set_xc_type_libxc", |
64 | 99 | "XC Functional involving Laplacian of rho is not implemented."); |
65 | 100 | } |
66 | | - int func_type; //0:none, 1:lda, 2:gga, 3:mgga, 4:hybrid lda/gga, 5:hybrid mgga |
| 101 | + |
| 102 | + // check if the functional involves non-local dispersion |
67 | 103 | if(not_supported_xc_with_nonlocal_vdw(xc_func_in)) |
68 | | - { ModuleBase::WARNING_QUIT("XC_Functional::set_xc_type_libxc","functionals with non-local dispersion are not supported."); } |
69 | | - func_type = 1; |
70 | | - if(xc_func_in.find("GGA") != std::string::npos) { func_type = 2; } |
71 | | - if(xc_func_in.find("MGGA") != std::string::npos) { func_type = 3; } |
72 | | - if(xc_func_in.find("HYB") != std::string::npos) { func_type =4; } |
73 | | - if(xc_func_in.find("HYB") != std::string::npos && xc_func_in.find("MGGA") != std::string::npos) { func_type =5; } |
74 | | - |
75 | | - // determine the id |
76 | | - std::vector<int> func_id; // libxc id of functional |
77 | | - size_t pos = 0; |
78 | | - std::string delimiter = "+"; |
79 | | - std::string token; |
80 | | - while ((pos = xc_func_in.find(delimiter)) != std::string::npos) |
81 | | - { |
82 | | - token = xc_func_in.substr(0, pos); |
83 | | - int id = xc_functional_get_number(token.c_str()); |
84 | | - std::cout << "func,id" << token << " " << id << std::endl; |
85 | | - if (id == -1) |
86 | | - { |
87 | | - std::string message = "Unrecognized exchange-correlation functional '"+ xc_func_in +"'.\n" |
88 | | - " Possible source: Pseudopotential file or dft_functional parameter.\n" |
89 | | - " Please explicitly set dft_functional in INPUT,\n" |
90 | | - " or verify the functional name is supported."; |
91 | | - ModuleBase::WARNING_QUIT("XC_Functional::set_xc_type_libxc",message); |
92 | | - } |
93 | | - func_id.push_back(id); |
94 | | - xc_func_in.erase(0, pos + delimiter.length()); |
95 | | - } |
96 | | - int id = xc_functional_get_number(xc_func_in.c_str()); |
97 | | - std::cout << "func,id" << xc_func_in << " " << id << std::endl; |
98 | | - if (id == -1) |
99 | 104 | { |
100 | | - std::string message = "Unrecognized exchange-correlation functional '"+ xc_func_in +"'.\n" |
| 105 | + ModuleBase::WARNING_QUIT("XC_Functional::set_xc_type_libxc", |
| 106 | + "functionals with non-local dispersion are not supported."); |
| 107 | + } |
| 108 | + |
| 109 | + // check the consistency of the functional type (LDA, GGA, MGGA, HYB_LDA, HYB_GGA, HYB_MGGA) |
| 110 | + const std::vector<std::string> xcfunc_words_ = FmtCore::split(xc_func_in, "+"); |
| 111 | + std::vector<int> xcfunc_type_(xcfunc_words_.size(), 0); // 0: None |
| 112 | + std::transform(xcfunc_words_.begin(), xcfunc_words_.end(), xcfunc_type_.begin(), |
| 113 | + [](const std::string& func) { return xc_func_type_classifier(func); }); |
| 114 | + if (std::adjacent_find(xcfunc_type_.begin(), xcfunc_type_.end(), |
| 115 | + [](int a, int b) { return a != b; }) != xcfunc_type_.end()) |
| 116 | + { |
| 117 | + ModuleBase::WARNING_QUIT("XC_Functional::set_xc_type_libxc", |
| 118 | + "All exchange-correlation functionals must be of the same type" |
| 119 | + "(LDA, GGA, MGGA, HYB_LDA, HYB_GGA, HYB_MGGA)."); |
| 120 | + } |
| 121 | + |
| 122 | + // check if there is None (no, we dont check it) |
| 123 | + int func_type = xcfunc_type_.front(); // all functionals are of the same type |
| 124 | + // if (func_type == 0) |
| 125 | + // { |
| 126 | + // ModuleBase::WARNING_QUIT("XC_Functional::set_xc_type_libxc", |
| 127 | + // "Unrecognized functional type in '" + xc_func_in + "'."); |
| 128 | + // } |
| 129 | + |
| 130 | + // determine the functional id |
| 131 | + std::vector<int> func_id(xcfunc_words_.size(), -1); |
| 132 | + std::transform(xcfunc_words_.begin(), xcfunc_words_.end(), func_id.begin(), |
| 133 | + [](const std::string& func) { return xc_functional_get_number(func.c_str()); }); |
| 134 | + // if there is any -1, it means the functional is not recognized |
| 135 | + const bool not_recognized_xc = std::any_of(func_id.begin(), func_id.end(), |
| 136 | + [](int id) { return id == -1; }); |
| 137 | + if (not_recognized_xc) |
| 138 | + { |
| 139 | + std::string message = "Unrecognized exchange-correlation functional '" + xc_func_in + "'.\n" |
101 | 140 | " Possible source: Pseudopotential file or dft_functional parameter.\n" |
102 | 141 | " Please explicitly set dft_functional in INPUT,\n" |
103 | 142 | " or verify the functional name is supported."; |
104 | | - ModuleBase::WARNING_QUIT("XC_Functional::set_xc_type_libxc",message); |
| 143 | + ModuleBase::WARNING_QUIT("XC_Functional::set_xc_type_libxc", message); |
105 | 144 | } |
106 | | - func_id.push_back(id); |
107 | 145 |
|
| 146 | + // return |
108 | 147 | return std::make_pair(func_type, func_id); |
109 | 148 | } |
110 | 149 |
|
111 | | -std::vector<xc_func_type> XC_Functional_Libxc::init_func(const std::vector<int> &func_id, const int xc_polarized) |
| 150 | +const std::vector<double> in_built_xc_func_ext_params(const int id) |
112 | 151 | { |
113 | | - // 'funcs' is the return value |
114 | | - std::vector<xc_func_type> funcs; |
115 | | - |
116 | | - //------------------------------------------- |
117 | | - // define a function named 'add_func', which |
118 | | - // initialize a functional according to its ID |
119 | | - //------------------------------------------- |
120 | | - auto add_func = [&]( const int func_id ) |
121 | | - { |
122 | | - funcs.push_back({}); |
123 | | - // 'xc_func_init' is defined in Libxc |
124 | | - xc_func_init( &funcs.back(), func_id, xc_polarized ); |
| 152 | + const std::map<int, std::vector<double>> mymap = { |
| 153 | + // finite temperature XC functionals |
| 154 | + {XC_LDA_XC_KSDT, {PARAM.inp.xc_temperature * 0.5}}, |
| 155 | + {XC_LDA_XC_CORRKSDT, {PARAM.inp.xc_temperature * 0.5}}, |
| 156 | + {XC_LDA_XC_GDSMFB, {PARAM.inp.xc_temperature * 0.5}}, |
| 157 | + // hybrid functionals |
| 158 | +#ifdef __EXX |
| 159 | + {XC_HYB_GGA_XC_PBEH, {GlobalC::exx_info.info_global.hybrid_alpha, |
| 160 | + GlobalC::exx_info.info_global.hse_omega, |
| 161 | + GlobalC::exx_info.info_global.hse_omega}}, |
| 162 | + {XC_HYB_GGA_XC_HSE06, {GlobalC::exx_info.info_global.hybrid_alpha, |
| 163 | + GlobalC::exx_info.info_global.hse_omega, |
| 164 | + GlobalC::exx_info.info_global.hse_omega}}, |
| 165 | + // short-range of B88_X |
| 166 | + {XC_GGA_X_ITYH, {PARAM.inp.exx_hse_omega}}, |
| 167 | + // short-range of LYP_C |
| 168 | + {XC_GGA_C_LYPR, {0.04918, 0.132, 0.2533, 0.349, |
| 169 | + 0.35/2.29, 2.0/2.29, PARAM.inp.exx_hse_omega}}, |
| 170 | +#endif |
125 | 171 | }; |
| 172 | + auto it = mymap.find(id); |
| 173 | + return (it != mymap.end()) ? it->second : std::vector<double>{}; |
| 174 | +} |
126 | 175 |
|
127 | | - for(int id : func_id) |
128 | | - { |
129 | | - if(id == XC_LDA_XC_KSDT || id == XC_LDA_XC_CORRKSDT || id == XC_LDA_XC_GDSMFB) //finite temperature XC functionals |
| 176 | +const std::vector<double> external_xc_func_ext_params(const int id) |
| 177 | +{ |
| 178 | + const std::map<int, std::vector<double>> mymap = { |
130 | 179 | { |
131 | | - add_func(id); |
132 | | - double parameter_finitet[1] = {PARAM.inp.xc_temperature * 0.5}; // converts to Hartree for libxc |
133 | | - xc_func_set_ext_params(&funcs.back(), parameter_finitet); |
134 | | - } |
135 | | -#ifdef __EXX |
136 | | - else if( id == XC_HYB_GGA_XC_PBEH ) // PBE0 |
137 | | - { |
138 | | - add_func( XC_HYB_GGA_XC_PBEH ); |
139 | | - double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha, |
140 | | - GlobalC::exx_info.info_global.hse_omega, |
141 | | - GlobalC::exx_info.info_global.hse_omega }; |
142 | | - xc_func_set_ext_params(&funcs.back(), parameter_hse); |
143 | | - } |
144 | | - else if( id == XC_HYB_GGA_XC_HSE06 ) // HSE06 hybrid functional |
145 | | - { |
146 | | - add_func( XC_HYB_GGA_XC_HSE06 ); |
147 | | - double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha, |
148 | | - GlobalC::exx_info.info_global.hse_omega, |
149 | | - GlobalC::exx_info.info_global.hse_omega }; |
150 | | - xc_func_set_ext_params(&funcs.back(), parameter_hse); |
151 | | - } |
152 | | - // added by jghan, 2024-07-06 |
153 | | - else if( id == XC_GGA_X_ITYH ) // short-range of B88_X |
154 | | - { |
155 | | - add_func( XC_GGA_X_ITYH ); |
156 | | - double parameter_omega[1] = {PARAM.inp.exx_hse_omega}; // GlobalC::exx_info.info_global.hse_omega |
157 | | - xc_func_set_ext_params(&funcs.back(), parameter_omega); |
158 | | - } |
159 | | - else if( id == XC_GGA_C_LYPR ) // short-range of LYP_C |
160 | | - { |
161 | | - add_func( XC_GGA_C_LYPR ); |
162 | | - // the first six parameters come from libxc, and may need to be modified in some cases |
163 | | - double parameter_lypr[7] = {0.04918, 0.132, 0.2533, 0.349, 0.35/2.29, 2.0/2.29, PARAM.inp.exx_hse_omega}; |
164 | | - xc_func_set_ext_params(&funcs.back(), parameter_lypr); |
| 180 | + PARAM.inp.xc_exch_ext[0], |
| 181 | + std::vector<double>(PARAM.inp.xc_exch_ext.begin()+1, |
| 182 | + PARAM.inp.xc_exch_ext.end()) |
| 183 | + }, |
| 184 | + { |
| 185 | + PARAM.inp.xc_corr_ext[0], |
| 186 | + std::vector<double>(PARAM.inp.xc_corr_ext.begin()+1, |
| 187 | + PARAM.inp.xc_corr_ext.end()) |
165 | 188 | } |
166 | | -#endif |
167 | | - else |
| 189 | + }; |
| 190 | + auto it = mymap.find(id); |
| 191 | + return (it != mymap.end()) ? it->second : std::vector<double>{}; |
| 192 | +} |
| 193 | + |
| 194 | +std::vector<xc_func_type> |
| 195 | +XC_Functional_Libxc::init_func(const std::vector<int> &func_id, |
| 196 | + const int xc_polarized) |
| 197 | +{ |
| 198 | + std::vector<xc_func_type> funcs; |
| 199 | + for (int id : func_id) |
| 200 | + { |
| 201 | + funcs.push_back({}); // create placeholder |
| 202 | + xc_func_init(&funcs.back(), id, xc_polarized); // instantiate the XC term |
| 203 | + |
| 204 | + // search for external parameters |
| 205 | + const std::vector<double> in_built_ext_params = in_built_xc_func_ext_params(id); |
| 206 | + const std::vector<double> external_ext_params = external_xc_func_ext_params(id); |
| 207 | + // for temporary use, I name their size as n1 and n2 |
| 208 | + const int n1 = in_built_ext_params.size(); |
| 209 | + const int n2 = external_ext_params.size(); |
| 210 | + |
| 211 | +// #ifdef __DEBUG // will the following assertion cause performance issue? |
| 212 | + // assert the number of parameters should be either zero or the value from |
| 213 | + // libxc function xc_func_info_get_n_ext_params, this is to avoid the undefined |
| 214 | + // behavior of illegal memory access |
| 215 | + const xc_func_info_type* info = xc_func_get_info(&funcs.back()); |
| 216 | + const int nref = xc_func_info_get_n_ext_params(info); |
| 217 | + assert ((n1 == 0) || (n1 == nref) || (n2 == 0) || (n2 == nref)); |
| 218 | +// #endif |
| 219 | + |
| 220 | + // external overwrites in-built if the same functional id is found in both maps |
| 221 | + const double* xc_func_ext_params = |
| 222 | + (n2 > 0) ? external_ext_params.data() : |
| 223 | + (n1 > 0) ? in_built_ext_params.data() : |
| 224 | + nullptr; // nullptr if no external parameters are found |
| 225 | + |
| 226 | + // if there are no external parameters, do nothing, otherwise we set |
| 227 | + if(xc_func_ext_params != nullptr) |
168 | 228 | { |
169 | | - add_func( id ); |
| 229 | + // set the external parameters |
| 230 | + xc_func_set_ext_params(&funcs.back(), const_cast<double*>(xc_func_ext_params)); |
170 | 231 | } |
171 | 232 | } |
172 | 233 | return funcs; |
|
0 commit comments