From 28ac9ebd474299774261c95ae9dc3aa53cecea1c Mon Sep 17 00:00:00 2001 From: Mathieu Taillefumier Date: Fri, 1 Dec 2023 15:56:09 +0100 Subject: [PATCH] [API] Add one function that parses QE hubbard.dat file and initialize the config structure accordingly --- src/CMakeLists.txt | 1 + src/api/sirius_api.cpp | 79 ++++++------- src/context/config.hpp | 15 +++ src/context/input_schema.json | 8 +- src/hubbard/hubbard.cpp | 1 - src/hubbard/hubbard.hpp | 5 + src/hubbard/hubbard_utils.cpp | 215 ++++++++++++++++++++++++++++++++++ 7 files changed, 277 insertions(+), 47 deletions(-) create mode 100644 src/hubbard/hubbard_utils.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 7dc6441e4..aa9fdb8de 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -42,6 +42,7 @@ set(_SOURCES "hubbard/hubbard_occupancies_derivatives.cpp" "hubbard/hubbard_potential_energy.cpp" "hubbard/hubbard_matrix.cpp" + "hubbard/hubbard_utils.cpp" "potential/generate_d_operator_matrix.cpp" "potential/generate_pw_coeffs.cpp" "potential/paw_potential.cpp" diff --git a/src/api/sirius_api.cpp b/src/api/sirius_api.cpp index c3b1f4aff..eee592fab 100644 --- a/src/api/sirius_api.cpp +++ b/src/api/sirius_api.cpp @@ -1961,7 +1961,6 @@ sirius_add_atom_type(void* const* handler__, char const* label__, char const* fn std::string label = std::string(label__); std::string fname = (fname__ == nullptr) ? std::string("") : std::string(fname__); sim_ctx.unit_cell().add_atom_type(label, fname); - auto& type = sim_ctx.unit_cell().atom_type(label); if (zn__ != nullptr) { type.set_zn(*zn__); @@ -5795,50 +5794,40 @@ sirius_add_hubbard_atom_pair(void* const* handler__, int* const atom_pair__, int { call_sirius( [&]() { - auto& sim_ctx = get_sim_ctx(handler__); - auto conf_dict = sim_ctx.cfg().hubbard(); + auto& sim_ctx = get_sim_ctx(handler__); + sirius::add_hubbard_atom_pair(sim_ctx, atom_pair__, translation__, n__, l__, *coupling__); + }, + error_code__); +} - json elem; - std::vector atom_pair(atom_pair__, atom_pair__ + 2); - /* Fortran indices start from 1 */ - atom_pair[0] -= 1; - atom_pair[1] -= 1; - std::vector n(n__, n__ + 2); - std::vector l(l__, l__ + 2); - std::vector translation(translation__, translation__ + 3); - - elem["atom_pair"] = atom_pair; - elem["T"] = translation; - elem["n"] = n; - elem["l"] = l; - elem["V"] = *coupling__; - - bool test{false}; - - for (int idx = 0; idx < conf_dict.nonlocal().size(); idx++) { - auto v = conf_dict.nonlocal(idx); - auto at_pr = v.atom_pair(); - /* search if the pair is already present */ - if ((at_pr[0] == atom_pair[0]) && (at_pr[1] == atom_pair[1])) { - auto tr = v.T(); - if ((tr[0] = translation[0]) && (tr[1] = translation[1]) && (tr[2] = translation[2])) { - auto lvl = v.n(); - if ((lvl[0] == n[0]) && (lvl[0] == n[1])) { - auto li = v.l(); - if ((li[0] == l[0]) && (li[1] == l[1])) { - test = true; - break; - } - } - } - } - } +/* +@api_begin +sirius_add_hubbard_file: + doc: Read all information for the hubbard correction from a QE file + arguments: + handler: + type: ctx_handler + attr: in, required + doc: Simulation context handler. + file_name: + type: string + attr: in, required + doc: name of the file containing the information + error_code: + type: int + attr: out, optional + doc: Error code. +@api end +*/ - if (!test) { - conf_dict.nonlocal().append(elem); - } else { - RTE_THROW("Atom pair for hubbard correction is already present"); - } +void +sirius_parse_hubbard_file(void* const* handler__, char* const file_name__, int* const error_code__) +{ + call_sirius( + [&]() { + auto& sim_ctx = get_sim_ctx(handler__); + std::string file_name(file_name__); + sirius::parse_hubbard_file(sim_ctx, file_name); }, error_code__); } @@ -5876,7 +5865,7 @@ sirius_add_hubbard_atom_pair(void* const* handler__, int* const atom_pair__, int type: int attr: out, optional doc: Error code. - @api end +@api end */ void sirius_set_hubbard_contrained_parameters(void* const* handler__, double const* hubbard_conv_thr__, @@ -5946,7 +5935,7 @@ sirius_set_hubbard_contrained_parameters(void* const* handler__, double const* h type: int attr: out, optional doc: Error code. - @api end +@api end */ void sirius_add_hubbard_atom_constraint(void* const* handler__, int* const atom_id__, int* const n__, int* const l__, diff --git a/src/context/config.hpp b/src/context/config.hpp index 55a6dcc63..f93d130f0 100644 --- a/src/context/config.hpp +++ b/src/context/config.hpp @@ -1587,6 +1587,21 @@ class config_t : dict_(dict__) { } + /// Name of the input file containing all information about hubbard orbital for U and V contributions. + /** + Name of the input file containing all information about hubbard orbital for U and V contribution. The schema should follow QE standard + */ + inline auto input_file() const + { + return dict_.at("/hubbard/input_file"_json_pointer).get(); + } + inline void input_file(std::string input_file__) + { + if (dict_.contains("locked")) { + throw std::runtime_error(locked_msg); + } + dict_["/hubbard/input_file"_json_pointer] = input_file__; + } /// Method to use for generating the hubbard subspace. [none] bare hubbard orbitals, [full_orthogonaliation] use all atomic wave functions to generate the hubbard subspace, [normalize] normalize the original hubbard wave functions, [orthogonalize] orthogonalize the hubbard wave functions inline auto hubbard_subspace_method() const { diff --git a/src/context/input_schema.json b/src/context/input_schema.json index 83018cf89..26d787d1f 100644 --- a/src/context/input_schema.json +++ b/src/context/input_schema.json @@ -854,6 +854,12 @@ "type": "object", "title": "Hubbard U correction", "properties": { + "input_file": { + "type": "string", + "default": "", + "title": "Name of the input file containing all information about hubbard orbital for U and V contributions.", + "description": "Name of the input file containing all information about hubbard orbital for U and V contribution. The schema should follow QE standard" + }, "hubbard_subspace_method": { "type": "string", "default": "none", @@ -923,7 +929,7 @@ "default": -1, "title": "angular momentum" }, - "lm_order" : { + "lm_order": { "type": "array", "default": [], "items": { diff --git a/src/hubbard/hubbard.cpp b/src/hubbard/hubbard.cpp index d2921a43a..73415d9b2 100644 --- a/src/hubbard/hubbard.cpp +++ b/src/hubbard/hubbard.cpp @@ -10,5 +10,4 @@ Hubbard::Hubbard(Simulation_context& ctx__) return; } } - } // namespace sirius diff --git a/src/hubbard/hubbard.hpp b/src/hubbard/hubbard.hpp index 5bd712036..53be633cf 100644 --- a/src/hubbard/hubbard.hpp +++ b/src/hubbard/hubbard.hpp @@ -196,6 +196,11 @@ class Hubbard } }; +extern void +add_hubbard_atom_pair(Simulation_context& ctx__, int* const atom_pair__, int* const translation__, int* const n__, + int* const l__, const double coupling__); +extern void +parse_hubbard_file(Simulation_context& ctx__, const std::string& data_file__); } // namespace sirius #endif // __HUBBARD_HPP__ diff --git a/src/hubbard/hubbard_utils.cpp b/src/hubbard/hubbard_utils.cpp new file mode 100644 index 000000000..77ebcf446 --- /dev/null +++ b/src/hubbard/hubbard_utils.cpp @@ -0,0 +1,215 @@ +#include +#include +#include +#include +#include +#include "context/simulation_context.hpp" + +namespace sirius { + +static std::vector +split(std::string str, std::string token) +{ + std::vector result; + while (str.size()) { + int index = str.find(token); + if (index != std::string::npos) { + result.push_back(str.substr(0, index)); + str = str.substr(index + token.size()); + if (str.size() == 0) + result.push_back(str); + } else { + result.push_back(str); + str = ""; + } + } + return result; +} + +static auto +parse_atom_string(std::string& at_str__) -> std::tuple +{ + std::vector split_str_ = split(at_str__, "_"); + int n_ = -1; + int l_ = -1; + + // a bit pedantic.... + n_ = std::stoi(split_str_[1]); + + switch (split_str_[1][1]) { + case 's': + case 'S': + l_ = 0; + break; + case 'p': + case 'P': + l_ = 1; + break; + case 'd': + case 'D': + l_ = 2; + break; + case 'f': + case 'F': + l_ = 3; + break; + case 'g': + case 'G': + l_ = 4; + break; + case 'h': + case 'H': + l_ = 5; + break; + case 'I': + case 'i': + l_ = 6; + break; + } + + return std::make_tuple(std::string{split_str_[0]}, n_, l_); +} + +void +add_hubbard_atom_pair(Simulation_context& ctx__, int* const atom_pair__, int* const translation__, int* const n__, + int* const l__, const double coupling__) +{ + json elem; + std::vector atom_pair(atom_pair__, atom_pair__ + 2); + /* Fortran indices start from 1 */ + atom_pair[0] -= 1; + atom_pair[1] -= 1; + std::vector n(n__, n__ + 2); + std::vector l(l__, l__ + 2); + std::vector translation(translation__, translation__ + 3); + + elem["atom_pair"] = atom_pair; + elem["T"] = translation; + elem["n"] = n; + elem["l"] = l; + elem["V"] = coupling__; + + bool test{false}; + + for (int idx = 0; idx < ctx__.cfg().hubbard().nonlocal().size(); idx++) { + auto v = ctx__.cfg().hubbard().nonlocal(idx); + auto at_pr = v.atom_pair(); + /* search if the pair is already present */ + if ((at_pr[0] == atom_pair[0]) && (at_pr[1] == atom_pair[1])) { + auto tr = v.T(); + if ((tr[0] == translation[0]) && (tr[1] == translation[1]) && (tr[2] == translation[2])) { + auto lvl = v.n(); + if ((lvl[0] == n[0]) && (lvl[0] == n[1])) { + auto li = v.l(); + if ((li[0] == l[0]) && (li[1] == l[1])) { + test = true; + break; + } + } + } + } + } + + if (!test) { + ctx__.cfg().hubbard().nonlocal().append(elem); + } else { + RTE_THROW("Atom pair for hubbard correction is already present"); + } +} + +void +parse_hubbard_file(Simulation_context& ctx__, const std::string& data_file__) +{ + std::ifstream infile(data_file__); + + if (!infile.is_open()) { + std::stringstream s; + s << "The file " << data_file__ << " is unreadable."; + RTE_THROW(s); + } + + std::string processed_str; + // std::vector> construct_neighbors_(ctx__.unit_cell().num_atoms() * 27); + + // for (int atom_id_ = 0 ; atom_id_ < ctx__.unit_cell().num_atoms(); atom_id_++) { + // construct_neighbors_[atom_id_] = std::make_tuple(0, 0, 0, atom_id_); + // } + + // index_ = ctx__.unit_cell().num_atoms(); + // we only need this list to initialize the input dictionary + // we can get the reverse from the index with this + + // for (int nx = -1; nx < 2; nx++) { + // for (int ny = -1; ny < 2; ny++) { + // for (int nz = -1; nz < 2; nz++) { + // for (int at = 0; at < ctx__.unit_cell().num_atoms(); at++) { + // construct_neighbors_[index_] = std::make_tuple(nz, + // ny, + // nx, + // ((3 * (nx + 1) + (ny + 1)) * 3 + nz + + // 1) * ctx__.unit_cell().num_atoms() + + // at); + // index_++; + // } + // } + // } + // } + + while (std::getline(infile, processed_str)) { + std::vector split_string = split(processed_str, std::string(" ")); + + if ((split_string.size() != 3) || (split_string.size() != 6) || (split_string[0] != "#")) { + std::stringstream s; + s << "The file " << data_file__ + << " seems to be corrupted.\nEach line should either be\nU atom_type-orbital real\n or \nV atom_type_1 " + "atom_type_2 int int real\n"; + RTE_THROW(s); + } + + if ((split_string[0] == "U") || (split_string[0] == "u")) { + const double U_ = std::stod(split_string[2]); + std::string atom_type_; + int n_; + int l_; + std::tie(atom_type_, n_, l_) = parse_atom_string(split_string[1]); + json elem; + elem["atom_type"] = atom_type_; + elem["n"] = n_; + elem["l"] = l_; + elem["total_initial_occupancy"] = 0; + elem["U"] = U_; + ctx__.cfg().hubbard().local().append(elem); + } + + if (((split_string[0] == "V") || (split_string[0] == "v"))) { + std::string atom_type1_; + std::string atom_type2_; + int n_pair_[2]; + int l_pair_[2]; + std::tie(atom_type1_, n_pair_[0], l_pair_[0]) = parse_atom_string(split_string[1]); + std::tie(atom_type2_, n_pair_[1], l_pair_[1]) = parse_atom_string(split_string[2]); + const int atom_index1_ = std::stoi(split_string[3]); + const int atom_index2_ = std::stoi(split_string[4]); + const double V_ = std::stoi(split_string[5]); + int translation1[3] = {0, 0, 0}; + int atom_pair_[2]; + int translation2[3] = {0, 0, 0}; + atom_pair_[0] = atom_index1_ % ctx__.unit_cell().num_atoms(); + translation1[0] = (atom_index1_ / ctx__.unit_cell().num_atoms()) % 3 - 1; + translation1[1] = (atom_index1_ / 3 * ctx__.unit_cell().num_atoms()) % 3 - 1; + translation1[2] = (atom_index1_ / 9 * ctx__.unit_cell().num_atoms()) % 3 - 1; + + atom_pair_[1] = atom_index1_ % ctx__.unit_cell().num_atoms(); + translation2[0] = (atom_index2_ / ctx__.unit_cell().num_atoms()) % 3 - 1; + translation2[1] = (atom_index2_ / 3 * ctx__.unit_cell().num_atoms()) % 3 - 1; + translation2[2] = (atom_index2_ / 9 * ctx__.unit_cell().num_atoms()) % 3 - 1; + + translation2[0] -= translation1[0]; + translation2[1] -= translation1[1]; + translation2[2] -= translation1[2]; + add_hubbard_atom_pair(ctx__, atom_pair_, translation2, n_pair_, l_pair_, V_); + } + } + // ignore all the other lines +} +} // namespace sirius