Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[API] Add function for parsing hubbard.dat style files #933

Draft
wants to merge 1 commit into
base: develop
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
79 changes: 34 additions & 45 deletions src/api/sirius_api.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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__);
Expand Down Expand Up @@ -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<int> atom_pair(atom_pair__, atom_pair__ + 2);
/* Fortran indices start from 1 */
atom_pair[0] -= 1;
atom_pair[1] -= 1;
std::vector<int> n(n__, n__ + 2);
std::vector<int> l(l__, l__ + 2);
std::vector<int> 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__);
}
Expand Down Expand Up @@ -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__,
Expand Down Expand Up @@ -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__,
Expand Down
15 changes: 15 additions & 0 deletions src/context/config.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<std::string>();
}
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
{
Expand Down
8 changes: 7 additions & 1 deletion src/context/input_schema.json
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -923,7 +929,7 @@
"default": -1,
"title": "angular momentum"
},
"lm_order" : {
"lm_order": {
"type": "array",
"default": [],
"items": {
Expand Down
1 change: 0 additions & 1 deletion src/hubbard/hubbard.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -10,5 +10,4 @@ Hubbard::Hubbard(Simulation_context& ctx__)
return;
}
}

} // namespace sirius
5 changes: 5 additions & 0 deletions src/hubbard/hubbard.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -196,6 +196,11 @@ class Hubbard
}
};

extern void
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this are normal functions on the c++ side. No raw pointers as arguments

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__
215 changes: 215 additions & 0 deletions src/hubbard/hubbard_utils.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,215 @@
#include <cstdio>
#include <cstdlib>
#include <string>
#include <vector>
#include <tuple>
#include "context/simulation_context.hpp"

namespace sirius {

static std::vector<std::string>
split(std::string str, std::string token)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes it will. I will use it.

{
std::vector<std::string> 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::string, int, int>
{
std::vector<std::string> split_str_ = split(at_str__, "_");
int n_ = -1;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why underscore in local variables? it is not needed

int l_ = -1;

// a bit pedantic....
n_ = std::stoi(split_str_[1]);

switch (split_str_[1][1]) {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

there is a code repetition with https://github.com/electronic-structure/SIRIUS/blob/develop/src/unit_cell/atom_type.cpp#L412

better solution: introduce a pure function that returns l by charachter label. It can be defined for example in atom_type.hpp

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__,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

no pointers as arguments

int* const l__, const double coupling__)
{
json elem;
std::vector<int> atom_pair(atom_pair__, atom_pair__ + 2);
/* Fortran indices start from 1 */
atom_pair[0] -= 1;
atom_pair[1] -= 1;
std::vector<int> n(n__, n__ + 2);
std::vector<int> l(l__, l__ + 2);
std::vector<int> 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__)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

our convention is T const& arg

{
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<std::tuple<int, int, int, int>> 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<int, int, int, int>(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<int, int, int, int>(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<std::string> 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]);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

local variables should not have underscore

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};
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this should be r3::vector T;

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