From e4c7408a58fea0e5ffc491e7611825ca8ee3ce8f Mon Sep 17 00:00:00 2001 From: Anton Kozhevnikov Date: Mon, 2 Dec 2024 18:50:25 +0100 Subject: [PATCH] [feature] proper I/O from/to hdf5 (#1028) * save config to hdf5 * restart from hdf5 * serialize species to json * simplify hdf5 write of density matrix * do not save potential as it can be generated from the density; use less strict convergence check --- apps/mini_app/sirius.scf.cpp | 38 ++++- apps/utils/unit_cell_tools.cpp | 4 +- examples/pp-pw/Si7Ge/sirius.json | 5 +- python_module/py_sirius.cpp | 2 - src/api/sirius_api.cpp | 4 +- src/context/simulation_context.cpp | 25 ++- src/context/simulation_context.hpp | 48 ++++-- src/context/simulation_parameters.cpp | 13 +- src/core/hdf5_tree.hpp | 67 +++++++- src/core/json.hpp | 11 ++ src/density/density.cpp | 62 +++++++ src/density/density.hpp | 31 +--- src/dft/dft_ground_state.cpp | 5 +- src/dft/dft_ground_state.hpp | 2 +- src/hubbard/hubbard_matrix.cpp | 18 +- src/hubbard/hubbard_matrix.hpp | 5 +- src/potential/potential.cpp | 50 +----- src/potential/potential.hpp | 6 - src/unit_cell/atom_type.cpp | 237 +++++++++++++++++++++++++- src/unit_cell/atom_type.hpp | 188 ++++---------------- src/unit_cell/unit_cell.cpp | 5 +- src/unit_cell/unit_cell.hpp | 2 +- 22 files changed, 528 insertions(+), 300 deletions(-) diff --git a/apps/mini_app/sirius.scf.cpp b/apps/mini_app/sirius.scf.cpp index b5b5caeb8..7777db5ca 100644 --- a/apps/mini_app/sirius.scf.cpp +++ b/apps/mini_app/sirius.scf.cpp @@ -82,19 +82,24 @@ preprocess_json_input(std::string fname__) std::unique_ptr create_sim_ctx(std::string fname__, cmd_args const& args__) { - auto json = preprocess_json_input(fname__); + std::string config_string; + if (isHDF5(fname__)) { + config_string = fname__; + } else { + auto json = preprocess_json_input(fname__); + config_string = json.dump(); + } - auto ctx_ptr = std::make_unique(json.dump(), mpi::Communicator::world()); - Simulation_context& ctx = *ctx_ptr; + auto ctx = std::make_unique(config_string); - auto& inp = ctx.cfg().parameters(); + auto& inp = ctx->cfg().parameters(); if (inp.gamma_point() && !(inp.ngridk()[0] * inp.ngridk()[1] * inp.ngridk()[2] == 1)) { RTE_THROW("this is not a Gamma-point calculation") } - ctx.import(args__); + ctx->import(args__); - return ctx_ptr; + return ctx; } auto @@ -110,6 +115,12 @@ ground_state(Simulation_context& ctx, int task_id, cmd_args const& args, int wri << "+----------------------+" << std::endl; break; } + case task_t::ground_state_restart: { + ctx.out() << "+--------------------------+" << std::endl + << "| restart SCF ground state |" << std::endl + << "+--------------------------+" << std::endl; + break; + } case task_t::ground_state_new_relax: { ctx.out() << "+---------------------------------------------+" << std::endl << "| new SCF ground state with atomic relaxation |" << std::endl @@ -145,11 +156,18 @@ ground_state(Simulation_context& ctx, int task_id, cmd_args const& args, int wri auto& density = dft.density(); if (task_id == task_t::ground_state_restart) { - if (!file_exists(storage_file_name)) { + auto fname = args.value("input", storage_file_name); + if (!isHDF5(fname)) { + fname = storage_file_name; + } + if (!file_exists(fname)) { RTE_THROW("storage file is not found"); } - density.load(storage_file_name); - potential.load(storage_file_name); + density.load(fname); + density.generate_paw_density(); + potential.generate(density, ctx.use_symmetry(), true); + Hamiltonian0 H0(potential, true); + initialize_subspace(kset, H0); } else { dft.initial_state(); } @@ -579,6 +597,8 @@ main(int argn, char** argv) {"iterative_solver.orthogonalize=", ""}, {"iterative_solver.early_restart=", "{double} value between 0 and 1 to control the early restart ratio in Davidson"}, + {"iterative_solver.energy_tolerance=", "{double} starting tolerance of iterative solver"}, + {"iterative_solver.num_steps=", "{int} number of steps in iterative solver"}, {"mixer.type=", "{string} mixer name (anderson, anderson_stable, broyden2, linear)"}, {"mixer.beta=", "{double} mixing parameter"}, {"volume_scale0=", "{double} starting volume scale for EOS calculation"}, diff --git a/apps/utils/unit_cell_tools.cpp b/apps/utils/unit_cell_tools.cpp index 224fc8624..7afec764d 100644 --- a/apps/utils/unit_cell_tools.cpp +++ b/apps/utils/unit_cell_tools.cpp @@ -34,7 +34,7 @@ create_supercell(cmd_args const& args__) std::cout << std::endl; } - Simulation_context ctx("sirius.json", mpi::Communicator::self()); + Simulation_context ctx(std::string("sirius.json"), mpi::Communicator::self()); auto scell_lattice_vectors = dot(ctx.unit_cell().lattice_vectors(), r3::matrix(scell)); @@ -113,7 +113,7 @@ create_supercell(cmd_args const& args__) void find_primitive() { - Simulation_context ctx("sirius.json", mpi::Communicator::self()); + Simulation_context ctx(std::string("sirius.json"), mpi::Communicator::self()); double lattice[3][3]; for (int i : {0, 1, 2}) { diff --git a/examples/pp-pw/Si7Ge/sirius.json b/examples/pp-pw/Si7Ge/sirius.json index bb1ee2f70..e7287bf77 100644 --- a/examples/pp-pw/Si7Ge/sirius.json +++ b/examples/pp-pw/Si7Ge/sirius.json @@ -16,7 +16,7 @@ "gk_cutoff" : 5.0, "pw_cutoff" : 20.00, - "energy_tol" : 1e-8, + "energy_tol" : 1e-7, "density_tol" : 1e-7, "num_dft_iter" : 100, @@ -33,7 +33,8 @@ "mixer" : { "beta" : 0.8, "type" : "anderson", - "max_history" : 8 + "max_history" : 8, + "use_hartree" : true }, "unit_cell": { diff --git a/python_module/py_sirius.cpp b/python_module/py_sirius.cpp index 889714268..195f97f6b 100644 --- a/python_module/py_sirius.cpp +++ b/python_module/py_sirius.cpp @@ -301,8 +301,6 @@ PYBIND11_MODULE(py_sirius, m) .def(py::init(), py::keep_alive<1, 2>(), "ctx"_a) .def("generate", &Potential::generate, "density"_a, "use_sym"_a, "transform_to_rg"_a) .def("fft_transform", &Potential::fft_transform) - .def("save", &Potential::save) - .def("load", &Potential::load) .def_property("vxc", py::overload_cast<>(&Potential::xc_potential), py::overload_cast<>(&Potential::xc_potential), py::return_value_policy::reference_internal) .def_property("exc", py::overload_cast<>(&Potential::xc_energy_density), diff --git a/src/api/sirius_api.cpp b/src/api/sirius_api.cpp index f6bf3c624..1e30e761b 100644 --- a/src/api/sirius_api.cpp +++ b/src/api/sirius_api.cpp @@ -2819,7 +2819,7 @@ sirius_generate_effective_potential(void* const* gs_handler__, int* error_code__ call_sirius( [&]() { auto& gs = get_gs(gs_handler__); - gs.potential().generate(gs.density(), gs.ctx().use_symmetry(), false); + gs.potential().generate(gs.density(), gs.ctx().use_symmetry(), true); }, error_code__); } @@ -6454,7 +6454,6 @@ sirius_save_state(void** gs_handler__, const char* file_name__, int* error_code_ auto& gs = get_gs(gs_handler__); std::string file_name(file_name__); gs.ctx().create_storage_file(file_name); - gs.potential().save(file_name); gs.density().save(file_name); }, error_code__); @@ -6486,7 +6485,6 @@ sirius_load_state(void** gs_handler__, const char* file_name__, int* error_code_ [&]() { auto& gs = get_gs(gs_handler__); std::string file_name(file_name__); - gs.potential().load(file_name); gs.density().load(file_name); }, error_code__); diff --git a/src/context/simulation_context.cpp b/src/context/simulation_context.cpp index e8d3e4586..0f59a030a 100644 --- a/src/context/simulation_context.cpp +++ b/src/context/simulation_context.cpp @@ -1166,18 +1166,18 @@ Simulation_context::update() void Simulation_context::create_storage_file(std::string name__) const { + if (!initialized_) { + RTE_THROW("Simulation_context must be initialized first"); + } if (comm_.rank() == 0) { /* create new hdf5 file */ HDF5_tree fout(name__, hdf5_access_t::truncate); fout.create_node("parameters"); - fout.create_node("effective_potential"); - fout.create_node("effective_magnetic_field"); fout.create_node("density"); fout.create_node("magnetization"); for (int j = 0; j < num_mag_dims(); j++) { fout["magnetization"].create_node(j); - fout["effective_magnetic_field"].create_node(j); } fout["parameters"].write("num_spins", num_spins()); @@ -1196,10 +1196,23 @@ Simulation_context::create_storage_file(std::string name__) const fout.create_node("unit_cell"); fout["unit_cell"].create_node("atoms"); - for (int j = 0; j < unit_cell().num_atoms(); j++) { - fout["unit_cell"]["atoms"].create_node(j); - fout["unit_cell"]["atoms"][j].write("mt_basis_size", unit_cell().atom(j).mt_basis_size()); + fout["unit_cell"].create_node("atom_types"); + for (int ia = 0; ia < unit_cell().num_atoms(); ia++) { + fout["unit_cell"]["atoms"].create_node(ia); + fout["unit_cell"]["atoms"][ia].write("mt_basis_size", unit_cell().atom(ia).mt_basis_size()); } + for (int iat = 0; iat < unit_cell().num_atom_types(); iat++) { + fout["unit_cell"]["atom_types"].create_node(iat); + fout["unit_cell"]["atom_types"][iat].write("config", unit_cell().atom_type(iat).serialize().dump()); + } + auto config = this->serialize()["config"]; + config.erase("locked"); + config["control"].erase("mpi_grid_dims"); + config["control"].erase("fft_mode"); + config["control"].erase("gen_evp_solver_name"); + config["control"].erase("std_evp_solver_name"); + config["settings"].erase("fft_grid_size"); + fout.write("config", config.dump()); } comm_.barrier(); } diff --git a/src/context/simulation_context.hpp b/src/context/simulation_context.hpp index e83335037..6c5bcdb83 100644 --- a/src/context/simulation_context.hpp +++ b/src/context/simulation_context.hpp @@ -318,6 +318,8 @@ class Simulation_context : public Simulation_parameters init_common(); } + /// Create an empty simulation context with an explicit communicators for k-point and + /// band parallelisation. Simulation_context(mpi::Communicator const& comm__, mpi::Communicator const& comm_k__, mpi::Communicator const& comm_band__) : comm_(comm__) @@ -327,29 +329,41 @@ class Simulation_context : public Simulation_parameters init_common(); } - /// Create a simulation context with world communicator and load parameters from JSON string or JSON file. - Simulation_context(std::string const& str__) - : comm_(mpi::Communicator::world()) - { - init_common(); - import(str__); - unit_cell_->import(cfg().unit_cell()); - } - - explicit Simulation_context(nlohmann::json const& dict__) - : comm_(mpi::Communicator::world()) + /// Create a simulation context with world communicator and load parameters from JSON string or a file. + explicit Simulation_context(std::string const& str__, mpi::Communicator const& comm__ = mpi::Communicator::world()) + : comm_(comm__) { init_common(); - import(dict__); - unit_cell_->import(cfg().unit_cell()); + if (!is_json_string(str__) && isHDF5(str__)) { + HDF5_tree fin(str__, hdf5_access_t::read_only); + std::string json_string; + fin.read("config", json_string); + auto dict = read_json_from_file_or_string(json_string); + for (auto& e : dict["unit_cell"]["atom_types"]) { + auto label = e.get(); + dict["unit_cell"]["atom_files"][label] = ""; + } + import(dict); + unit_cell_->import(cfg().unit_cell()); + + /* need to set type of calculation before parsing species */ + this->electronic_structure_method(cfg().parameters().electronic_structure_method()); + for (int iat = 0; iat < unit_cell_->num_atom_types(); iat++) { + fin["unit_cell"]["atom_types"][iat].read("config", json_string); + unit_cell_->atom_type(iat).read_input(json_string); + } + } else { + import(str__); + unit_cell_->import(cfg().unit_cell()); + } } - // /// Create a simulation context with world communicator and load parameters from JSON string or JSON file. - Simulation_context(std::string const& str__, mpi::Communicator const& comm__) + explicit Simulation_context(nlohmann::json const& dict__, + mpi::Communicator const& comm__ = mpi::Communicator::world()) : comm_(comm__) { init_common(); - import(str__); + import(dict__); unit_cell_->import(cfg().unit_cell()); } @@ -764,7 +778,7 @@ class Simulation_context : public Simulation_parameters /// Export parameters of simulation context as a JSON dictionary. nlohmann::json - serialize() + serialize() const { nlohmann::json dict; dict["config"] = cfg().dict(); diff --git a/src/context/simulation_parameters.cpp b/src/context/simulation_parameters.cpp index bff4d4127..d31700e3e 100644 --- a/src/context/simulation_parameters.cpp +++ b/src/context/simulation_parameters.cpp @@ -142,16 +142,16 @@ get_section_options(std::string const& section__) } void -Simulation_parameters::import(std::string const& str__) +Simulation_parameters::import(nlohmann::json const& dict__) { - auto json = read_json_from_file_or_string(str__); - import(json); + cfg_.import(dict__); } void -Simulation_parameters::import(nlohmann::json const& dict__) +Simulation_parameters::import(std::string const& str__) { - cfg_.import(dict__); + auto dict = read_json_from_file_or_string(str__); + this->import(dict); } void @@ -174,6 +174,9 @@ Simulation_parameters::import(cmd_args const& args__) cfg_.iterative_solver().early_restart( args__.value("iterative_solver.early_restart", cfg_.iterative_solver().early_restart())); + cfg_.iterative_solver().energy_tolerance( + args__.value("iterative_solver.energy_tolerance", cfg_.iterative_solver().energy_tolerance())); + cfg_.iterative_solver().num_steps(args__.value("iterative_solver.num_steps", cfg_.iterative_solver().num_steps())); cfg_.mixer().beta(args__.value("mixer.beta", cfg_.mixer().beta())); cfg_.mixer().type(args__.value("mixer.type", cfg_.mixer().type())); } diff --git a/src/core/hdf5_tree.hpp b/src/core/hdf5_tree.hpp index cf9565095..c5e9e262a 100644 --- a/src/core/hdf5_tree.hpp +++ b/src/core/hdf5_tree.hpp @@ -70,6 +70,12 @@ struct hdf5_type_wrapper }; }; +inline bool +isHDF5(std::string const& filename__) +{ + return H5Fis_hdf5(filename__.c_str()) > 0; +} + /// Interface to the HDF5 library. class HDF5_tree { @@ -366,6 +372,38 @@ class HDF5_tree } } + /// Get dimensions of the dataset. + std::vector + dims(std::string const& name__) const + { + HDF5_group group(file_id_, path_); + + HDF5_dataset dataset(group.id(), name__); + + // Get the dataspace of the dataset + hid_t dataspace_id = H5Dget_space(dataset.id()); + if (dataspace_id < 0) { + RTE_THROW("Error getting dataspace"); + } + + // Get the number of dimensions (rank) of the dataset + int ndims = H5Sget_simple_extent_ndims(dataspace_id); + if (ndims < 0) { + RTE_THROW("Error getting number of dimensions"); + } + + // Get the dimensions of the dataset + hsize_t dims[ndims]; + if (H5Sget_simple_extent_dims(dataspace_id, dims, NULL) < 0) { + RTE_THROW("Error getting dimensions"); + } + std::vector result(ndims); + for (int i = 0; i < ndims; i++) { + result[i] = dims[i]; + } + return result; + } + public: /// Constructor to create the HDF5 tree. HDF5_tree(std::string const& file_name__, hdf5_access_t access__) @@ -504,6 +542,16 @@ class HDF5_tree write(name, &vec[0], (int)vec.size()); } + void + write(std::string const& name, std::string const& str) + { + mdarray s_char({str.size()}); + for (size_t i = 0; i < str.size(); i++) { + s_char[i] = str[i]; + } + this->write(name, s_char); + } + /// Write attribute template void @@ -640,14 +688,29 @@ class HDF5_tree read(name, &vec[0], (int)vec.size()); } - HDF5_tree + inline void + read(std::string const& name, std::string& str) + { + auto d = this->dims(name); + if (d.size() != 1) { + RTE_THROW("wrong size of config dataset"); + } + std::vector s_char(d[0]); + this->read(name, s_char); + str = std::string(d[0], ' '); + for (int i = 0; i < d[0]; i++) { + str[i] = s_char[i]; + } + } + + inline HDF5_tree operator[](std::string const& path__) { auto new_path = path_ + path__ + "/"; return HDF5_tree(file_id_, new_path); } - HDF5_tree + inline HDF5_tree operator[](int idx__) { auto new_path = path_ + std::to_string(idx__) + "/"; diff --git a/src/core/json.hpp b/src/core/json.hpp index 963d11c69..79f4ec6a1 100644 --- a/src/core/json.hpp +++ b/src/core/json.hpp @@ -39,6 +39,17 @@ try_parse(std::istream& is) return dict; } +inline bool +is_json_string(std::string const& str__) +{ + try { + auto json = nlohmann::json::parse(str__); + return true; + } catch (nlohmann::json::parse_error const&) { + } + return false; +} + inline nlohmann::json read_json_from_file(std::string const& filename) { diff --git a/src/density/density.cpp b/src/density/density.cpp index 9e9136faf..258130964 100644 --- a/src/density/density.cpp +++ b/src/density/density.cpp @@ -2040,4 +2040,66 @@ Density::print_info(std::ostream& out__) const // } } +void +Density::save(std::string name__) const +{ + rho().hdf5_write(name__, "density"); + for (int j = 0; j < ctx_.num_mag_dims(); j++) { + mag(j).hdf5_write(name__, "magnetization/" + std::to_string(j)); + } + ctx_.comm().barrier(); + if (ctx_.comm().rank() == 0) { + HDF5_tree fout(name__, hdf5_access_t::read_write); + fout.create_node("density_matrix"); + for (int ia = 0; ia < unit_cell_.num_atoms(); ia++) { + fout["density_matrix"].create_node(ia).write("data", this->density_matrix(ia)); + } + if (ctx_.hubbard_correction()) { + fout.create_node("occupation_matrix"); + fout["occupation_matrix"].create_node("local"); + fout["occupation_matrix"].create_node("nonlocal"); + for (size_t i = 0; i < this->occupation_matrix().local().size(); i++) { + fout["occupation_matrix"]["local"].create_node(i).write("data", this->occupation_matrix().local(i)); + } + for (size_t i = 0; i < this->occupation_matrix().nonlocal().size(); i++) { + fout["occupation_matrix"]["nonlocal"].create_node(i).write("data", + this->occupation_matrix().nonlocal(i)); + } + } + } +} + +void +Density::load(std::string name__) +{ + HDF5_tree fin(name__, hdf5_access_t::read_only); + + int ngv; + fin.read("/parameters/num_gvec", &ngv, 1); + if (ngv != ctx_.gvec().num_gvec()) { + RTE_THROW("wrong number of G-vectors"); + } + mdarray gv({3, ngv}); + fin.read("/parameters/gvec", gv); + + rho().hdf5_read(name__, "density", gv); + rho().rg().fft_transform(1); + for (int j = 0; j < ctx_.num_mag_dims(); j++) { + mag(j).hdf5_read(name__, "magnetization/" + std::to_string(j), gv); + mag(j).rg().fft_transform(1); + } + + for (int ia = 0; ia < unit_cell_.num_atoms(); ia++) { + fin["density_matrix"][ia].read("data", this->density_matrix(ia)); + } + if (ctx_.hubbard_correction()) { + for (size_t i = 0; i < this->occupation_matrix().local().size(); i++) { + fin["occupation_matrix"]["local"][i].read("data", this->occupation_matrix().local(i)); + } + for (size_t i = 0; i < this->occupation_matrix().nonlocal().size(); i++) { + fin["occupation_matrix"]["nonlocal"][i].read("data", this->occupation_matrix().nonlocal(i)); + } + } +} + } // namespace sirius diff --git a/src/density/density.hpp b/src/density/density.hpp index 63e6d2b6f..ca9d40014 100644 --- a/src/density/density.hpp +++ b/src/density/density.hpp @@ -600,37 +600,10 @@ class Density : public Field4D } void - save(std::string name__) const - { - rho().hdf5_write(name__, "density"); - for (int j = 0; j < ctx_.num_mag_dims(); j++) { - std::stringstream s; - s << "magnetization/" << j; - mag(j).hdf5_write(name__, s.str()); - } - ctx_.comm().barrier(); - } + save(std::string name__) const; void - load(std::string name__) - { - HDF5_tree fin(name__, hdf5_access_t::read_only); - - int ngv; - fin.read("/parameters/num_gvec", &ngv, 1); - if (ngv != ctx_.gvec().num_gvec()) { - RTE_THROW("wrong number of G-vectors"); - } - mdarray gv({3, ngv}); - fin.read("/parameters/gvec", gv); - - rho().hdf5_read(name__, "density", gv); - rho().rg().fft_transform(1); - for (int j = 0; j < ctx_.num_mag_dims(); j++) { - mag(j).hdf5_read(name__, "magnetization/" + std::to_string(j), gv); - mag(j).rg().fft_transform(1); - } - } + load(std::string name__); void save_to_xsf() diff --git a/src/dft/dft_ground_state.cpp b/src/dft/dft_ground_state.cpp index 68118b415..975cf2757 100644 --- a/src/dft/dft_ground_state.cpp +++ b/src/dft/dft_ground_state.cpp @@ -348,7 +348,8 @@ DFT_ground_state::find(double density_tol__, double energy_tol__, double iter_so ctx_.message(2, __func__, out); /* check if the calculation has converged */ bool converged{true}; - converged = (std::abs(eold - etot) < energy_tol__) && result.converged && iter_solver_converged; + // converged = (std::abs(eold - etot) < energy_tol__) && result.converged && iter_solver_converged; + converged = (std::abs(eold - etot) < energy_tol__) && iter_solver_converged; if (ctx_.cfg().mixer().use_hartree()) { converged = converged && (eha_res < density_tol__); } else { @@ -379,7 +380,7 @@ DFT_ground_state::find(double density_tol__, double energy_tol__, double iter_so density_.mag(j).rg().fft_transform(-1); } } - potential_.save(storage_file_name); + // potential_.save(storage_file_name); density_.save(storage_file_name); // kset_.save(storage_file_name); } diff --git a/src/dft/dft_ground_state.hpp b/src/dft/dft_ground_state.hpp index 511bafa39..9ee79b969 100644 --- a/src/dft/dft_ground_state.hpp +++ b/src/dft/dft_ground_state.hpp @@ -89,7 +89,7 @@ class DFT_ground_state n = ctx_.num_itsol_steps(); kset_.comm().allreduce(&n, 1); if (ctx_.verbosity() >= 2) { - RTE_OUT(ctx_.out()) << "numbef of iterative solver steps: " << n << std::endl; + RTE_OUT(ctx_.out()) << "number of iterative solver steps: " << n << std::endl; } } diff --git a/src/hubbard/hubbard_matrix.cpp b/src/hubbard/hubbard_matrix.cpp index f811dc2cf..bf3a9ead1 100644 --- a/src/hubbard/hubbard_matrix.cpp +++ b/src/hubbard/hubbard_matrix.cpp @@ -323,11 +323,11 @@ Hubbard_matrix::print_nonlocal(int idx__, std::ostream& out__) const void Hubbard_matrix::zero() { - for (int ia = 0; ia < static_cast(local_.size()); ia++) { - local_[ia].zero(); + for (int i = 0; i < static_cast(local_.size()); i++) { + local_[i].zero(); } - for (int i = 0; i < static_cast(ctx_.cfg().hubbard().nonlocal().size()); i++) { + for (int i = 0; i < static_cast(nonlocal_.size()); i++) { nonlocal_[i].zero(); } @@ -338,4 +338,16 @@ Hubbard_matrix::zero() } } +void +Hubbard_matrix::print(std::ostream& out__) const +{ + for (int i = 0; i < static_cast(local_.size()); i++) { + this->print_local(i, out__); + } + + for (int i = 0; i < static_cast(nonlocal_.size()); i++) { + this->print_nonlocal(i, out__); + } +} + } // namespace sirius diff --git a/src/hubbard/hubbard_matrix.hpp b/src/hubbard/hubbard_matrix.hpp index 7c6715eae..636f2a684 100644 --- a/src/hubbard/hubbard_matrix.hpp +++ b/src/hubbard/hubbard_matrix.hpp @@ -63,11 +63,14 @@ class Hubbard_matrix access(std::string const& what__, std::complex* ptr__, int ld__); void - print_local(int ia__, std::ostream& out__) const; + print_local(int idx__, std::ostream& out__) const; void print_nonlocal(int idx__, std::ostream& out__) const; + void + print(std::ostream& out__) const; + void zero(); diff --git a/src/potential/potential.cpp b/src/potential/potential.cpp index cbda1287c..328af402b 100644 --- a/src/potential/potential.cpp +++ b/src/potential/potential.cpp @@ -296,7 +296,7 @@ Potential::generate(Density const& density__, bool use_symmetry__, bool transfor * 1) compute D-matrix * 2) establish a mapping between fine and coarse FFT grid for the Hloc operator * 3) symmetrize effective potential */ - fft_transform(-1); + this->fft_transform(-1); } if (use_symmetry__) { @@ -372,54 +372,6 @@ Potential::generate(Density const& density__, bool use_symmetry__, bool transfor } } -void -Potential::save(std::string name__) -{ - effective_potential().hdf5_write(name__, "effective_potential"); - for (int j = 0; j < ctx_.num_mag_dims(); j++) { - effective_magnetic_field(j).hdf5_write(name__, "effective_magnetic_field/" + std::to_string(j)); - } - if (ctx_.comm().rank() == 0 && !ctx_.full_potential()) { - HDF5_tree fout(name__, hdf5_access_t::read_write); - for (int j = 0; j < ctx_.unit_cell().num_atoms(); j++) { - if (ctx_.unit_cell().atom(j).mt_basis_size() != 0) { - fout["unit_cell"]["atoms"][j].write("D_operator", ctx_.unit_cell().atom(j).d_mtrx()); - } - } - } - comm_.barrier(); -} - -void -Potential::load(std::string name__) -{ - HDF5_tree fin(name__, hdf5_access_t::read_only); - - int ngv; - fin.read("/parameters/num_gvec", &ngv, 1); - if (ngv != ctx_.gvec().num_gvec()) { - RTE_THROW("wrong number of G-vectors"); - } - mdarray gv({3, ngv}); - fin.read("/parameters/gvec", gv); - - effective_potential().hdf5_read(name__, "effective_potential", gv); - - for (int j = 0; j < ctx_.num_mag_dims(); j++) { - effective_magnetic_field(j).hdf5_read(name__, "effective_magnetic_field/" + std::to_string(j), gv); - } - - if (ctx_.full_potential()) { - update_atomic_potential(); - } - - if (!ctx_.full_potential()) { - for (int j = 0; j < ctx_.unit_cell().num_atoms(); j++) { - fin["unit_cell"]["atoms"][j].read("D_operator", ctx_.unit_cell().atom(j).d_mtrx()); - } - } -} - void Potential::update_atomic_potential() { diff --git a/src/potential/potential.hpp b/src/potential/potential.hpp index 364834a1a..a3b83bd5d 100644 --- a/src/potential/potential.hpp +++ b/src/potential/potential.hpp @@ -608,12 +608,6 @@ class Potential : public Field4D void generate(Density const& density__, bool use_sym__, bool transform_to_rg__); - void - save(std::string name__); - - void - load(std::string name__); - void update_atomic_potential(); diff --git a/src/unit_cell/atom_type.cpp b/src/unit_cell/atom_type.cpp index 6c5ae50a8..7afb82237 100644 --- a/src/unit_cell/atom_type.cpp +++ b/src/unit_cell/atom_type.cpp @@ -621,6 +621,9 @@ Atom_type::read_pseudo_uspp(nlohmann::json const& parser) std::istringstream iss(std::string(1, c1[0])); iss >> n; } + if (dict[k].count("principal_quantum_number")) { + n = dict[k]["principal_quantum_number"].get(); + } if (spin_orbit_coupling() && dict[k].count("total_angular_momentum") && l != 0) { @@ -643,7 +646,7 @@ Atom_type::read_pseudo_uspp(nlohmann::json const& parser) void Atom_type::read_pseudo_paw(nlohmann::json const& parser) { - is_paw_ = true; + this->is_paw_ = true; auto& header = parser["pseudo_potential"]["header"]; /* read core energy */ @@ -654,7 +657,7 @@ Atom_type::read_pseudo_paw(nlohmann::json const& parser) } /* cutoff index */ - int cutoff_radius_index = parser["pseudo_potential"]["header"]["cutoff_radius_index"].get(); + int cutoff_radius_index = parser["pseudo_potential"]["header"].value("cutoff_radius_index", -1); /* read core density and potential */ paw_ae_core_charge_density( @@ -679,7 +682,9 @@ Atom_type::read_pseudo_paw(nlohmann::json const& parser) RTE_THROW(s); } - add_ae_paw_wf(std::vector(wfc.begin(), wfc.begin() + cutoff_radius_index)); + int last = (cutoff_radius_index == -1) ? wfc.size() : cutoff_radius_index; + + add_ae_paw_wf({wfc.begin(), wfc.begin() + last}); wfc = parser["pseudo_potential"]["paw_data"]["ps_wfc"][i]["radial_function"].get>(); @@ -691,14 +696,15 @@ Atom_type::read_pseudo_paw(nlohmann::json const& parser) RTE_THROW(s); } - add_ps_paw_wf(std::vector(wfc.begin(), wfc.begin() + cutoff_radius_index)); + last = (cutoff_radius_index == -1) ? wfc.size() : cutoff_radius_index; + + add_ps_paw_wf({wfc.begin(), wfc.begin() + last}); } } void Atom_type::read_input(nlohmann::json const& parser) { - if (!parameters_.full_potential()) { read_pseudo_uspp(parser); @@ -1256,4 +1262,225 @@ Atom_type::add_hubbard_orbital(int n__, int l__, double occ__, double U, double std::move(s.interpolate()), use_for_calculations__, idx_rf); } +nlohmann::json +Atom_type::serialize() const +{ + nlohmann::json dict = {{"header", nlohmann::json::object()}}; + + dict["header"]["z_valence"] = zn_; + dict["header"]["mesh_size"] = this->num_mt_points(); + if (is_paw()) { + dict["header"]["pseudo_type"] = "PAW"; + } else if (augment()) { + dict["header"]["pseudo_type"] = "US"; + } else { + dict["header"]["pseudo_type"] = "NC"; + } + dict["header"]["number_of_proj"] = num_beta_radial_functions(); + dict["header"]["element"] = symbol_; + dict["radial_grid"] = radial_grid().values(); + dict["local_potential"] = local_potential(); + dict["core_charge_density"] = ps_core_charge_density(); + dict["total_charge_density"] = ps_total_charge_density(); + dict["atomic_wave_functions"] = nlohmann::json::array(); + for (auto& e : ps_atomic_wfs_) { + auto o = nlohmann::json::object(); + o["angular_momentum"] = e.am.l(); + o["radial_function"] = e.f.values(); + o["principal_quantum_number"] = e.n; + dict["atomic_wave_functions"].push_back(o); + } + dict["beta_projectors"] = nlohmann::json::array(); + for (auto& e : beta_radial_functions_) { + auto o = nlohmann::json::object(); + o["angular_momentum"] = e.first.l(); + o["radial_function"] = e.second.values(); + dict["beta_projectors"].push_back(o); + } + int nbf = num_beta_radial_functions(); + std::vector v(nbf * nbf); + for (int i = 0; i < nbf; i++) { + for (int j = 0; j < nbf; j++) { + v[j * nbf + i] = this->d_mtrx_ion()(i, j); + } + } + dict["D_ion"] = v; + if (augment_) { + dict["augmentation"] = nlohmann::json::array(); + for (int i = 0; i < nbf; i++) { + for (int j = i; j < nbf; j++) { + for (int l = 0; l <= 2 * lmax_beta(); l++) { + auto o = nlohmann::json::object(); + o["i"] = i; + o["j"] = j; + o["angular_momentum"] = l; + o["radial_function"] = q_radial_function(i, j, l).values(); + dict["augmentation"].push_back(o); + } + } + } + } + if (is_paw()) { + auto paw = nlohmann::json::object(); + paw["ae_core_charge_density"] = paw_ae_core_charge_density(); + paw["ae_wfc"] = nlohmann::json::array(); + for (auto& e : ae_paw_wfs_) { + auto o = nlohmann::json::object(); + o["radial_function"] = e; + paw["ae_wfc"].push_back(o); + } + paw["pw_wfc"] = nlohmann::json::array(); + for (auto& e : ps_paw_wfs_) { + auto o = nlohmann::json::object(); + o["radial_function"] = e; + paw["ps_wfc"].push_back(o); + } + paw["occupations"] = paw_wf_occ_; + dict["paw_data"] = paw; + } + + nlohmann::json result = {{"pseudo_potential", dict}}; + + return result; +} + +void +Atom_type::init_aw_descriptors() +{ + RTE_ASSERT(this->lmax_apw() >= -1); + + if (this->lmax_apw() >= 0 && aw_default_l_.size() == 0) { + RTE_THROW("default AW descriptor is empty"); + } + + aw_descriptors_.clear(); + for (int l = 0; l <= this->lmax_apw(); l++) { + aw_descriptors_.push_back(aw_default_l_); + for (size_t ord = 0; ord < aw_descriptors_[l].size(); ord++) { + aw_descriptors_[l][ord].n = l + 1; + aw_descriptors_[l][ord].l = l; + } + } + + for (size_t i = 0; i < aw_specific_l_.size(); i++) { + int l = aw_specific_l_[i][0].l; + if (l < this->lmax_apw()) { + aw_descriptors_[l] = aw_specific_l_[i]; + } + } +} + +void +Atom_type::add_aw_descriptor(int n, int l, double enu, int dme, int auto_enu) +{ + if (static_cast(aw_descriptors_.size()) < (l + 1)) { + aw_descriptors_.resize(l + 1, radial_solution_descriptor_set()); + } + + radial_solution_descriptor rsd; + + rsd.n = n; + if (n == -1) { + /* default principal quantum number value for any l */ + rsd.n = l + 1; + for (int ist = 0; ist < num_atomic_levels(); ist++) { + /* take next level after the core */ + if (atomic_level(ist).core && atomic_level(ist).l == l) { + rsd.n = atomic_level(ist).n + 1; + } + } + } + + rsd.l = l; + rsd.dme = dme; + rsd.enu = enu; + rsd.auto_enu = auto_enu; + aw_descriptors_[l].push_back(rsd); +} + +void +Atom_type::add_lo_descriptor(int ilo, int n, int l, double enu, int dme, int auto_enu) +{ + if ((int)lo_descriptors_.size() == ilo) { + angular_momentum am(l); + lo_descriptors_.push_back(local_orbital_descriptor(am)); + } else { + if (l != lo_descriptors_[ilo].am.l()) { + std::stringstream s; + s << "wrong angular quantum number" << std::endl + << "atom type id: " << id() << " (" << symbol_ << ")" << std::endl + << "idxlo: " << ilo << std::endl + << "n: " << l << std::endl + << "l: " << n << std::endl + << "expected l: " << lo_descriptors_[ilo].am.l() << std::endl; + RTE_THROW(s); + } + } + + radial_solution_descriptor rsd; + + rsd.n = n; + if (n == -1) { + /* default value for any l */ + rsd.n = l + 1; + for (int ist = 0; ist < num_atomic_levels(); ist++) { + if (atomic_level(ist).core && atomic_level(ist).l == l) { + /* take next level after the core */ + rsd.n = atomic_level(ist).n + 1; + } + } + } + + rsd.l = l; + rsd.dme = dme; + rsd.enu = enu; + rsd.auto_enu = auto_enu; + lo_descriptors_[ilo].rsd_set.push_back(rsd); +} + +void +Atom_type::add_ps_atomic_wf(int n__, angular_momentum am__, std::vector f__, double occ__) +{ + Spline rwf(radial_grid_, f__); + auto d = std::sqrt(inner(rwf, rwf, 0, radial_grid_.num_points())); + if (d < 1e-4) { + std::stringstream s; + s << "small norm (" << d << ") of radial atomic pseudo wave-function for n=" << n__ << " and j=" << am__.j(); + RTE_THROW(s); + } + + ps_atomic_wfs_.emplace_back(n__, am__, occ__, std::move(rwf)); +} + +void +Atom_type::add_q_radial_function(int idxrf1__, int idxrf2__, int l__, std::vector qrf__) +{ + /* sanity check */ + if (l__ > 2 * lmax_beta()) { + std::stringstream s; + s << "wrong l for Q radial functions of atom type " << label_ << std::endl + << "current l: " << l__ << std::endl + << "lmax_beta: " << lmax_beta() << std::endl + << "maximum allowed l: " << 2 * lmax_beta(); + + RTE_THROW(s); + } + + if (!augment_) { + /* once we add a Q-radial function, we need to augment the charge */ + augment_ = true; + /* number of radial beta-functions */ + int nbrf = num_beta_radial_functions(); + q_radial_functions_l_ = mdarray, 2>({nbrf * (nbrf + 1) / 2, 2 * lmax_beta() + 1}); + + for (int l = 0; l <= 2 * lmax_beta(); l++) { + for (int idx = 0; idx < nbrf * (nbrf + 1) / 2; idx++) { + q_radial_functions_l_(idx, l) = Spline(radial_grid_); + } + } + } + + q_radial_functions_l_(packed_index(idxrf1__, idxrf2__), l__) = Spline(radial_grid_, qrf__); +} + } // namespace sirius diff --git a/src/unit_cell/atom_type.hpp b/src/unit_cell/atom_type.hpp index d6eaaa91b..121eedcc3 100644 --- a/src/unit_cell/atom_type.hpp +++ b/src/unit_cell/atom_type.hpp @@ -281,13 +281,6 @@ class Atom_type inline void read_pseudo_paw(nlohmann::json const& parser); - /// Read atomic parameters from json file. - inline void - read_input(std::string const& str__); - - inline void - read_input(nlohmann::json const& parser); - /// Read atomic parameters directly from UPF v2 files #ifdef SIRIUS_USE_PUGIXML inline void @@ -301,31 +294,8 @@ class Atom_type #endif /// Initialize descriptors of the augmented-wave radial functions. - inline void - init_aw_descriptors() - { - RTE_ASSERT(this->lmax_apw() >= -1); - - if (this->lmax_apw() >= 0 && aw_default_l_.size() == 0) { - RTE_THROW("default AW descriptor is empty"); - } - - aw_descriptors_.clear(); - for (int l = 0; l <= this->lmax_apw(); l++) { - aw_descriptors_.push_back(aw_default_l_); - for (size_t ord = 0; ord < aw_descriptors_[l].size(); ord++) { - aw_descriptors_[l][ord].n = l + 1; - aw_descriptors_[l][ord].l = l; - } - } - - for (size_t i = 0; i < aw_specific_l_.size(); i++) { - int l = aw_specific_l_[i][0].l; - if (l < this->lmax_apw()) { - aw_descriptors_[l] = aw_specific_l_[i]; - } - } - } + void + init_aw_descriptors(); /* forbid copy constructor */ Atom_type(Atom_type const& src) = delete; @@ -376,6 +346,13 @@ class Atom_type /// Move constructor. Atom_type(Atom_type&& src) = default; + /// Read atomic parameters from json file. + void + read_input(std::string const& str__); + + void + read_input(nlohmann::json const& parser); + /// Initialize the atom type. /** Once the unit cell is populated with all atom types and atoms, each atom type can be initialized. */ void @@ -399,6 +376,28 @@ class Atom_type void print_info(std::ostream& out__) const; + /// Serialize atom type information to json + nlohmann::json + serialize() const; + + /// Add augmented-wave descriptor. + void + add_aw_descriptor(int n, int l, double enu, int dme, int auto_enu); + + /// Add local orbital descriptor + void + add_lo_descriptor(int ilo, int n, int l, double enu, int dme, int auto_enu); + + /// Add atomic radial function to the list. + void + add_ps_atomic_wf(int n__, angular_momentum am__, std::vector f__, double occ__ = 0.0); + + /// Add radial function of the augmentation charge. + /** Radial functions of beta projectors must be added already. Their total number will be used to + deterimine the storage size for the radial functions of the augmented charge. */ + inline void + add_q_radial_function(int idxrf1__, int idxrf2__, int l__, std::vector qrf__); + /// Set the radial grid of the given type. inline void set_radial_grid(radial_grid_t grid_type__, int num_points__, double rmin__, double rmax__, double p__) @@ -442,76 +441,6 @@ class Atom_type return atomic_levels_[idx]; } - /// Add augmented-wave descriptor. - inline void - add_aw_descriptor(int n, int l, double enu, int dme, int auto_enu) - { - if (static_cast(aw_descriptors_.size()) < (l + 1)) { - aw_descriptors_.resize(l + 1, radial_solution_descriptor_set()); - } - - radial_solution_descriptor rsd; - - rsd.n = n; - if (n == -1) { - /* default principal quantum number value for any l */ - rsd.n = l + 1; - for (int ist = 0; ist < num_atomic_levels(); ist++) { - /* take next level after the core */ - if (atomic_level(ist).core && atomic_level(ist).l == l) { - rsd.n = atomic_level(ist).n + 1; - } - } - } - - rsd.l = l; - rsd.dme = dme; - rsd.enu = enu; - rsd.auto_enu = auto_enu; - aw_descriptors_[l].push_back(rsd); - } - - /// Add local orbital descriptor - inline void - add_lo_descriptor(int ilo, int n, int l, double enu, int dme, int auto_enu) - { - if ((int)lo_descriptors_.size() == ilo) { - angular_momentum am(l); - lo_descriptors_.push_back(local_orbital_descriptor(am)); - } else { - if (l != lo_descriptors_[ilo].am.l()) { - std::stringstream s; - s << "wrong angular quantum number" << std::endl - << "atom type id: " << id() << " (" << symbol_ << ")" << std::endl - << "idxlo: " << ilo << std::endl - << "n: " << l << std::endl - << "l: " << n << std::endl - << "expected l: " << lo_descriptors_[ilo].am.l() << std::endl; - RTE_THROW(s); - } - } - - radial_solution_descriptor rsd; - - rsd.n = n; - if (n == -1) { - /* default value for any l */ - rsd.n = l + 1; - for (int ist = 0; ist < num_atomic_levels(); ist++) { - if (atomic_level(ist).core && atomic_level(ist).l == l) { - /* take next level after the core */ - rsd.n = atomic_level(ist).n + 1; - } - } - } - - rsd.l = l; - rsd.dme = dme; - rsd.enu = enu; - rsd.auto_enu = auto_enu; - lo_descriptors_[ilo].rsd_set.push_back(rsd); - } - /// Add the entire local orbital descriptor. inline void add_lo_descriptor(local_orbital_descriptor const& lod__) @@ -519,22 +448,6 @@ class Atom_type lo_descriptors_.push_back(lod__); } - /// Add atomic radial function to the list. - inline void - add_ps_atomic_wf(int n__, angular_momentum am__, std::vector f__, double occ__ = 0.0) - { - Spline rwf(radial_grid_, f__); - auto d = std::sqrt(inner(rwf, rwf, 0, radial_grid_.num_points())); - if (d < 1e-4) { - std::stringstream s; - s << "small norm (" << d << ") of radial atomic pseudo wave-function for n=" << n__ - << " and j=" << am__.j(); - RTE_THROW(s); - } - - ps_atomic_wfs_.emplace_back(n__, am__, occ__, std::move(rwf)); - } - /// Return a tuple describing a given atomic radial function auto const& ps_atomic_wf(int idx__) const @@ -548,7 +461,10 @@ class Atom_type add_beta_radial_function(angular_momentum am__, std::vector beta__) { if (augment_) { - RTE_THROW("can't add more beta projectors"); + std::stringstream s; + s << "augmentation charge has already been added" << std::endl + << "this fixes the number of beta projectors to " << this->num_beta_radial_functions() << std::endl + << "can't add more beta projectors" << std::endl; } Spline s(radial_grid_, beta__); beta_radial_functions_.push_back(std::make_pair(am__, std::move(s))); @@ -568,40 +484,6 @@ class Atom_type return beta_radial_functions_[idxrf__]; } - /// Add radial function of the augmentation charge. - /** Radial functions of beta projectors must be added already. Their total number will be used to - deterimine the storage size for the radial functions of the augmented charge. */ - inline void - add_q_radial_function(int idxrf1__, int idxrf2__, int l__, std::vector qrf__) - { - /* sanity check */ - if (l__ > 2 * lmax_beta()) { - std::stringstream s; - s << "wrong l for Q radial functions of atom type " << label_ << std::endl - << "current l: " << l__ << std::endl - << "lmax_beta: " << lmax_beta() << std::endl - << "maximum allowed l: " << 2 * lmax_beta(); - - RTE_THROW(s); - } - - if (!augment_) { - /* once we add a Q-radial function, we need to augment the charge */ - augment_ = true; - /* number of radial beta-functions */ - int nbrf = num_beta_radial_functions(); - q_radial_functions_l_ = mdarray, 2>({nbrf * (nbrf + 1) / 2, 2 * lmax_beta() + 1}); - - for (int l = 0; l <= 2 * lmax_beta(); l++) { - for (int idx = 0; idx < nbrf * (nbrf + 1) / 2; idx++) { - q_radial_functions_l_(idx, l) = Spline(radial_grid_); - } - } - } - - q_radial_functions_l_(packed_index(idxrf1__, idxrf2__), l__) = Spline(radial_grid_, qrf__); - } - /// Return true if this atom type has an augementation charge. inline bool augment() const diff --git a/src/unit_cell/unit_cell.cpp b/src/unit_cell/unit_cell.cpp index afbe3e398..e9cc91624 100644 --- a/src/unit_cell/unit_cell.cpp +++ b/src/unit_cell/unit_cell.cpp @@ -379,7 +379,8 @@ Unit_cell::serialize(bool cart_pos__) const if (cart_pos__) { v = dot(lattice_vectors_, v); } - dict["atoms"][atom_type(iat).label()].push_back({v[0], v[1], v[2]}); + auto f = atom(ia).vector_field(); + dict["atoms"][atom_type(iat).label()].push_back({v[0], v[1], v[2], f[0], f[1], f[2]}); } } return dict; @@ -596,7 +597,7 @@ Unit_cell::generate_radial_integrals() } std::string -Unit_cell::chemical_formula() +Unit_cell::chemical_formula() const { std::string name; for (int iat = 0; iat < num_atom_types(); iat++) { diff --git a/src/unit_cell/unit_cell.hpp b/src/unit_cell/unit_cell.hpp index 881533cc3..b56d64829 100644 --- a/src/unit_cell/unit_cell.hpp +++ b/src/unit_cell/unit_cell.hpp @@ -242,7 +242,7 @@ class Unit_cell /// Get a simple simple chemical formula bases on the total unit cell. /** Atoms of each type are counted and packed in a string. For example, O2Ni2 or La2O4Cu */ std::string - chemical_formula(); + chemical_formula() const; /// Update the parameters that depend on atomic positions or lattice vectors. void