From d9e9a714ff0b89e40c99be588d537addb158ede1 Mon Sep 17 00:00:00 2001 From: toxa81 Date: Wed, 18 Jul 2018 14:16:27 +0200 Subject: [PATCH] cleanup --- src/Beta_projectors/beta_projectors.h | 2 +- .../beta_projectors_strain_deriv.h | 4 +- src/Density/density.hpp | 8 +- src/Density/initial_density.hpp | 2 +- src/Density/paw_density.hpp | 4 +- src/Density/symmetrize_density_matrix.hpp | 6 +- .../hubbard_generate_atomic_orbitals.hpp | 2 +- .../hubbard_occupancies_derivatives.hpp | 14 +- src/Hubbard/hubbard_occupancy.hpp | 20 +- ...generate_atomic_centered_wavefunctions.hpp | 95 +----- src/K_point/initialize.hpp | 2 +- src/Potential/paw_potential.hpp | 25 +- src/Potential/potential.hpp | 6 +- src/Unit_cell/atom.h | 10 +- src/Unit_cell/atom_type.h | 7 +- src/Unit_cell/basis_functions_index.hpp | 2 +- src/Unit_cell/unit_cell_symmetry.h | 1 - src/augmentation_operator.h | 10 +- src/dft_ground_state.h | 1 - src/gaunt.h | 7 +- src/k_point.h | 7 +- src/matching_coefficients.h | 2 +- src/radial_integrals.h | 45 +-- src/sht.h | 17 +- src/simulation_parameters.h | 9 +- src/spheric_function.h | 16 +- src/utils.h | 292 ----------------- src/utils/utils.hpp | 301 +++++------------- src/xc_functional.h | 1 - 29 files changed, 177 insertions(+), 741 deletions(-) delete mode 100644 src/utils.h diff --git a/src/Beta_projectors/beta_projectors.h b/src/Beta_projectors/beta_projectors.h index e2f4ed57d..e69a016ea 100644 --- a/src/Beta_projectors/beta_projectors.h +++ b/src/Beta_projectors/beta_projectors.h @@ -62,7 +62,7 @@ class Beta_projectors: public Beta_projectors_base<1> /* vs = {r, theta, phi} */ auto vs = SHT::spherical_coordinates(gkvec_.gkvec_cart(igk)); /* compute real spherical harmonics for G+k vector */ - std::vector gkvec_rlm(Utils::lmmax(ctx_.unit_cell().lmax())); + std::vector gkvec_rlm(utils::lmmax(ctx_.unit_cell().lmax())); SHT::spherical_harmonics(ctx_.unit_cell().lmax(), vs[1], vs[2], &gkvec_rlm[0]); for (int iat = 0; iat < ctx_.unit_cell().num_atom_types(); iat++) { auto& atom_type = ctx_.unit_cell().atom_type(iat); diff --git a/src/Beta_projectors/beta_projectors_strain_deriv.h b/src/Beta_projectors/beta_projectors_strain_deriv.h index 9b77a597f..caa01836e 100644 --- a/src/Beta_projectors/beta_projectors_strain_deriv.h +++ b/src/Beta_projectors/beta_projectors_strain_deriv.h @@ -21,7 +21,7 @@ class Beta_projectors_strain_deriv : public Beta_projectors_base<9> auto& beta_ri1 = ctx_.beta_ri_djl(); int lmax = ctx_.unit_cell().lmax(); - int lmmax = Utils::lmmax(lmax); + int lmmax = utils::lmmax(lmax); mdarray rlm_g(lmmax, num_gkvec_loc()); mdarray rlm_dg(lmmax, 3, num_gkvec_loc()); @@ -142,7 +142,7 @@ class Beta_projectors_strain_deriv : public Beta_projectors_base<9> // auto gvs = SHT::spherical_coordinates(gvc); // /* compute real spherical harmonics for G+k vector */ - // std::vector gkvec_rlm(Utils::lmmax(lmax_beta_ + 2)); + // std::vector gkvec_rlm(utils::lmmax(lmax_beta_ + 2)); // SHT::spherical_harmonics(lmax_beta_ + 2, gvs[1], gvs[2], &gkvec_rlm[0]); // mdarray tmp(ctx_.unit_cell().max_mt_radial_basis_size(), lmax_beta_ + 3, ctx_.unit_cell().num_atom_types()); diff --git a/src/Density/density.hpp b/src/Density/density.hpp index 2738f490c..3236e2cec 100644 --- a/src/Density/density.hpp +++ b/src/Density/density.hpp @@ -181,7 +181,7 @@ class Density : public Field4D std::unique_ptr> gaunt_coefs_{nullptr}; /// Fast mapping between composite lm index and corresponding orbital quantum number. - mdarray l_by_lm_; + std::vector l_by_lm_; /// High-frequency mixer for the pseudopotential density mixing. std::unique_ptr> hf_mixer_{nullptr}; @@ -231,9 +231,9 @@ class Density : public Field4D int l1 = atom_type__.indexr(idxrf1).l; int xi2 = atom_type__.indexb().index_by_idxrf(idxrf2); - for (int lm2 = Utils::lm_by_l_m(l2, -l2); lm2 <= Utils::lm_by_l_m(l2, l2); lm2++, xi2++) { + for (int lm2 = utils::lm(l2, -l2); lm2 <= utils::lm(l2, l2); lm2++, xi2++) { int xi1 = atom_type__.indexb().index_by_idxrf(idxrf1); - for (int lm1 = Utils::lm_by_l_m(l1, -l1); lm1 <= Utils::lm_by_l_m(l1, l1); lm1++, xi1++) { + for (int lm1 = utils::lm(l1, -l1); lm1 <= utils::lm(l1, l1); lm1++, xi1++) { for (int k = 0; k < gaunt_coeffs__.num_gaunt(lm1, lm2); k++) { int lm3 = gaunt_coeffs__.gaunt(lm1, lm2, k).lm3; auto gc = gaunt_coeffs__.gaunt(lm1, lm2, k).coef; @@ -355,7 +355,7 @@ class Density : public Field4D gaunt_coefs_ = std::unique_ptr(new gc_z(ctx_.lmax_apw(), ctx_.lmax_rho(), ctx_.lmax_apw(), SHT::gaunt_hybrid)); } - l_by_lm_ = Utils::l_by_lm(ctx_.lmax_rho()); + l_by_lm_ = utils::l_by_lm(ctx_.lmax_rho()); density_matrix_ = mdarray(unit_cell_.max_mt_basis_size(), unit_cell_.max_mt_basis_size(), ctx_.num_mag_comp(), unit_cell_.num_atoms()); diff --git a/src/Density/initial_density.hpp b/src/Density/initial_density.hpp index f55b9f701..40e57e08f 100644 --- a/src/Density/initial_density.hpp +++ b/src/Density/initial_density.hpp @@ -230,7 +230,7 @@ inline void Density::initial_density_full_pot() //} //sba.approximate(gvec_len); - auto l_by_lm = Utils::l_by_lm(lmax); + auto l_by_lm = utils::l_by_lm(lmax); std::vector zil(lmax + 1); for (int l = 0; l <= lmax; l++) { diff --git a/src/Density/paw_density.hpp b/src/Density/paw_density.hpp index b7697bf0d..3405717a4 100644 --- a/src/Density/paw_density.hpp +++ b/src/Density/paw_density.hpp @@ -20,7 +20,7 @@ inline void Density::init_paw() auto& atom_type = atom.type(); int l_max = 2 * atom_type.indexr().lmax_lo(); - int lm_max_rho = Utils::lmmax(l_max); + int lm_max_rho = utils::lmmax(l_max); paw_density_data_t pdd; @@ -88,7 +88,7 @@ inline void Density::generate_paw_atom_density(paw_density_data_t& pdd) auto& atom_type = pdd.atom_->type(); - auto l_by_lm = Utils::l_by_lm(2 * atom_type.indexr().lmax_lo()); + auto l_by_lm = utils::l_by_lm(2 * atom_type.indexr().lmax_lo()); /* get gaunt coefficients */ Gaunt_coefficients GC(atom_type.indexr().lmax_lo(), 2 * atom_type.indexr().lmax_lo(), diff --git a/src/Density/symmetrize_density_matrix.hpp b/src/Density/symmetrize_density_matrix.hpp index 415c27be8..7a3b78862 100644 --- a/src/Density/symmetrize_density_matrix.hpp +++ b/src/Density/symmetrize_density_matrix.hpp @@ -11,7 +11,7 @@ inline void Density::symmetrize_density_matrix() dm.zero(); int lmax = unit_cell_.lmax(); - int lmmax = Utils::lmmax(lmax); + int lmmax = utils::lmmax(lmax); mdarray rotm(lmmax, lmmax); @@ -42,10 +42,10 @@ inline void Density::symmetrize_density_matrix() for (int j = 0; j < ndm; j++) { for (int m3 = -l1; m3 <= l1; m3++) { - int lm3 = Utils::lm_by_l_m(l1, m3); + int lm3 = utils::lm(l1, m3); int xi3 = atom_type.indexb().index_by_lm_order(lm3, o1); for (int m4 = -l2; m4 <= l2; m4++) { - int lm4 = Utils::lm_by_l_m(l2, m4); + int lm4 = utils::lm(l2, m4); int xi4 = atom_type.indexb().index_by_lm_order(lm4, o2); dm_rot_spatial[j] += density_matrix_(xi3, xi4, j, ja) * rotm(lm1, lm3) * rotm(lm2, lm4) * alpha; } diff --git a/src/Hubbard/hubbard_generate_atomic_orbitals.hpp b/src/Hubbard/hubbard_generate_atomic_orbitals.hpp index 5ddfd3130..203da22a1 100644 --- a/src/Hubbard/hubbard_generate_atomic_orbitals.hpp +++ b/src/Hubbard/hubbard_generate_atomic_orbitals.hpp @@ -32,7 +32,7 @@ void Hubbard_potential::generate_atomic_orbitals(K_point& kp, Q_operatornumber_of_hubbard_orbitals(), num_sc); - kp.generate_atomic_centered_wavefunctions_(this->number_of_hubbard_orbitals(), sphi, this->offset, true); + kp.generate_atomic_centered_wavefunctions_aux(this->number_of_hubbard_orbitals(), sphi, this->offset, true); // check if we have a norm conserving pseudo potential only bool augment = false; diff --git a/src/Hubbard/hubbard_occupancies_derivatives.hpp b/src/Hubbard/hubbard_occupancies_derivatives.hpp index b2159b0e7..db4252117 100644 --- a/src/Hubbard/hubbard_occupancies_derivatives.hpp +++ b/src/Hubbard/hubbard_occupancies_derivatives.hpp @@ -14,10 +14,10 @@ void Hubbard_potential::compute_occupancies_derivatives(K_point& kp, // derivatives of the hubbard wave functions are needed. auto &phi = kp.hubbard_wave_functions(); - kp.generate_atomic_centered_wavefunctions_(this->number_of_hubbard_orbitals(), - phi, - this->offset, - true); + kp.generate_atomic_centered_wavefunctions_aux(this->number_of_hubbard_orbitals(), + phi, + this->offset, + true); Beta_projectors_gradient bp_grad_(ctx_, kp.gkvec(), kp.igk_loc(), kp.beta_projectors()); kp.beta_projectors().prepare(); @@ -205,7 +205,7 @@ void Hubbard_potential::compute_occupancies_stress_derivatives(K_point& kp, bool augment = false; const int lmax = ctx_.unit_cell().lmax(); - const int lmmax = Utils::lmmax(lmax); + const int lmmax = utils::lmmax(lmax); mdarray rlm_g(lmmax, kp.num_gkvec_loc()); mdarray rlm_dg(lmmax, 3, kp.num_gkvec_loc()); @@ -225,7 +225,7 @@ void Hubbard_potential::compute_occupancies_stress_derivatives(K_point& kp, bp_strain_deriv.prepare(); // compute the hubbard orbitals - kp.generate_atomic_centered_wavefunctions_(this->number_of_hubbard_orbitals(), phi, this->offset, true); + kp.generate_atomic_centered_wavefunctions_aux(this->number_of_hubbard_orbitals(), phi, this->offset, true); #ifdef __GPU if (ctx_.processing_unit() == GPU) { @@ -403,7 +403,7 @@ void Hubbard_potential::compute_gradient_strain_wavefunctions(K_point& kp__, } } else { for (int m = -l; m <= l; m++) { - int lm = Utils::lm_by_l_m(l, m); + int lm = utils::lm(l, m); auto d1 = ri_values[atom_type.id()][i] * (gvc[mu] * rlm_dg(lm, nu, igkloc) + p * rlm_g(lm, igkloc)); auto d2 = ridjl_values[atom_type.id()][i] * rlm_g(lm, igkloc) * gvc[mu] * gvc[nu] / gvs[0]; diff --git a/src/Hubbard/hubbard_occupancy.hpp b/src/Hubbard/hubbard_occupancy.hpp index 1a172ddbb..4e92de2ef 100644 --- a/src/Hubbard/hubbard_occupancy.hpp +++ b/src/Hubbard/hubbard_occupancy.hpp @@ -379,7 +379,7 @@ inline void Hubbard_potential::symmetrize_occupancy_matrix_noncolinear_case() // check if we have some symmetries if (sym.num_mag_sym()) { int lmax = unit_cell_.lmax(); - int lmmax = Utils::lmmax(lmax); + int lmmax = utils::lmmax(lmax); mdarray rotm(lmmax, lmmax); mdarray rotated_oc(lmmax, lmmax, ctx_.num_spins() * ctx_.num_spins(), unit_cell_.num_atoms()); @@ -400,9 +400,9 @@ inline void Hubbard_potential::symmetrize_occupancy_matrix_noncolinear_case() const int lmax_at = 2 * atom.type().hubbard_l() + 1; if (atom.type().hubbard_correction()) { for (int ii = 0; ii < lmax_at; ii++) { - int l1 = Utils::lm_by_l_m(atom.type().hubbard_l(), ii - atom.type().hubbard_l()); + int l1 = utils::lm(atom.type().hubbard_l(), ii - atom.type().hubbard_l()); for (int ll = 0; ll < lmax_at; ll++) { - int l2 = Utils::lm_by_l_m(atom.type().hubbard_l(), ll - atom.type().hubbard_l()); + int l2 = utils::lm(atom.type().hubbard_l(), ll - atom.type().hubbard_l()); mdarray rot_spa(ctx_.num_spins() * ctx_.num_spins()); rot_spa.zero(); for (int s1 = 0; s1 < ctx_.num_spins(); s1++) { @@ -411,9 +411,9 @@ inline void Hubbard_potential::symmetrize_occupancy_matrix_noncolinear_case() // A_ij B_jk C_kl for (int jj = 0; jj < lmax_at; jj++) { - int l3 = Utils::lm_by_l_m(atom.type().hubbard_l(), jj - atom.type().hubbard_l()); + int l3 = utils::lm(atom.type().hubbard_l(), jj - atom.type().hubbard_l()); for (int kk = 0; kk < lmax_at; kk++) { - int l4 = Utils::lm_by_l_m(atom.type().hubbard_l(), kk - atom.type().hubbard_l()); + int l4 = utils::lm(atom.type().hubbard_l(), kk - atom.type().hubbard_l()); rot_spa(2 * s1 + s2) += std::conj(rotm(l1, l3)) * occupancy_number_(jj, kk, (s1 == s2) * s1 + (s1 != s2) * (1 + 2 * s1 + s2), ia, 0) * @@ -469,7 +469,7 @@ inline void Hubbard_potential::symmetrize_occupancy_matrix() // check if we have some symmetries if (sym.num_mag_sym()) { int lmax = unit_cell_.lmax(); - int lmmax = Utils::lmmax(lmax); + int lmmax = utils::lmmax(lmax); mdarray rotm(lmmax, lmmax); mdarray rotated_oc(lmmax, lmmax, ctx_.num_spins() * ctx_.num_spins(), unit_cell_.num_atoms()); @@ -493,15 +493,15 @@ inline void Hubbard_potential::symmetrize_occupancy_matrix() rot_spa.zero(); for (int ispn = 0; ispn < ctx_.num_spins(); ispn++) { for (int ii = 0; ii < lmax_at; ii++) { - int l1 = Utils::lm_by_l_m(atom.type().hubbard_l(), ii - atom.type().hubbard_l()); + int l1 = utils::lm(atom.type().hubbard_l(), ii - atom.type().hubbard_l()); for (int ll = 0; ll < lmax_at; ll++) { - int l2 = Utils::lm_by_l_m(atom.type().hubbard_l(), ll - atom.type().hubbard_l()); + int l2 = utils::lm(atom.type().hubbard_l(), ll - atom.type().hubbard_l()); // symmetrization procedure // A_ij B_jk C_kl for (int kk = 0; kk < lmax_at; kk++) { - int l4 = Utils::lm_by_l_m(atom.type().hubbard_l(), kk - atom.type().hubbard_l()); + int l4 = utils::lm(atom.type().hubbard_l(), kk - atom.type().hubbard_l()); for (int jj = 0; jj < lmax_at; jj++) { - int l3 = Utils::lm_by_l_m(atom.type().hubbard_l(), jj - atom.type().hubbard_l()); + int l3 = utils::lm(atom.type().hubbard_l(), jj - atom.type().hubbard_l()); rot_spa(ii, kk) += std::conj(rotm(l1, l3)) * occupancy_number_(jj, kk, ispn, ia, 0) * rotm(l2, l4) * alpha; } diff --git a/src/K_point/generate_atomic_centered_wavefunctions.hpp b/src/K_point/generate_atomic_centered_wavefunctions.hpp index 1e1097b3e..d317b257c 100644 --- a/src/K_point/generate_atomic_centered_wavefunctions.hpp +++ b/src/K_point/generate_atomic_centered_wavefunctions.hpp @@ -2,7 +2,10 @@ // this list should contain: index of atom, index of wave-function and some flag to indicate if we average // wave-functions in case of spin-orbit; this should be sufficient to generate a desired sub-set of atomic wave-functions -inline void K_point::generate_atomic_centered_wavefunctions_(const int num_ao__, Wave_functions& phi, std::vector &offset, const bool hubbard) +inline void K_point::generate_atomic_centered_wavefunctions_aux(const int num_ao__, + Wave_functions& phi, + std::vector& offset, + const bool hubbard) { if (!num_ao__) { return; @@ -22,7 +25,7 @@ inline void K_point::generate_atomic_centered_wavefunctions_(const int num_ao__, /* vs = {r, theta, phi} */ auto vs = SHT::spherical_coordinates(this->gkvec().gkvec_cart(igk)); /* compute real spherical harmonics for G+k vector */ - std::vector rlm(Utils::lmmax(lmax)); + std::vector rlm(utils::lmmax(lmax)); SHT::spherical_harmonics(lmax, vs[1], vs[2], &rlm[0]); /* get values of radial integrals for a given G+k vector length */ std::vector> ri_values(unit_cell_.num_atom_types()); @@ -40,7 +43,7 @@ inline void K_point::generate_atomic_centered_wavefunctions_(const int num_ao__, auto l = std::abs(atom_type.ps_atomic_wf(i).first); auto z = std::pow(double_complex(0, -1), l) * fourpi / std::sqrt(unit_cell_.omega()); for (int m = -l; m <= l; m++) { - int lm = Utils::lm_by_l_m(l, m); + int lm = utils::lm(l, m); phi.pw_coeffs(0).prime(igk_loc, n) = z * std::conj(phase_factor) * rlm[lm] * ri_values[atom_type.id()][i]; n++; } @@ -57,7 +60,7 @@ inline void K_point::generate_atomic_centered_wavefunctions_(const int num_ao__, auto z = std::pow(double_complex(0, -1), l) * fourpi / std::sqrt(unit_cell_.omega()); if (l == atom_type.hubbard_l()) { for (int m = -l; m <= l; m++) { - int lm = Utils::lm_by_l_m(l, m); + int lm = utils::lm(l, m); if (atom_type.spin_orbit_coupling()) { phi.pw_coeffs(0).prime(igk_loc, offset[ia] + l + m) += 0.5 * z * std::conj(phase_factor) * rlm[lm] * ri_values[atom_type.id()][i]; phi.pw_coeffs(1).prime(igk_loc, offset[ia] + 3 * l + m + 1) += 0.5 * z * std::conj(phase_factor) * rlm[lm] * ri_values[atom_type.id()][i]; @@ -76,12 +79,18 @@ inline void K_point::generate_atomic_centered_wavefunctions_(const int num_ao__, } // igk_loc } -inline void K_point::generate_atomic_centered_wavefunctions(const int num_ao__, Wave_functions& phi) { - std::vector vs(1,0); - generate_atomic_centered_wavefunctions_(num_ao__, phi, vs, false); +inline void K_point::generate_atomic_centered_wavefunctions(const int num_ao__, Wave_functions& phi) +{ + std::vector vs(1, 0); + generate_atomic_centered_wavefunctions_aux(num_ao__, phi, vs, false); } -inline void K_point::compute_gradient_wavefunctions(Wave_functions &phi, const int starting_position_i, const int num_wf, Wave_functions &dphi, const int starting_position_j, const int direction) { +inline void K_point::compute_gradient_wavefunctions(Wave_functions& phi, + const int starting_position_i, + const int num_wf, + Wave_functions& dphi, + const int starting_position_j, + const int direction) { std::vector qalpha(this->num_gkvec_loc()); for (int igk_loc = 0; igk_loc < this->num_gkvec_loc(); igk_loc++) { @@ -99,73 +108,3 @@ inline void K_point::compute_gradient_wavefunctions(Wave_functions &phi, const i } } } -//inline void K_point::generate_atomic_centered_wavefunctions(const int num_ao__, Wave_functions& phi) -//{ -// int lmax{0}; -// for (int iat = 0; iat < unit_cell_.num_atom_types(); iat++) { -// auto& atom_type = unit_cell_.atom_type(iat); -// for (auto& wf : atom_type.pp_desc().atomic_pseudo_wfs_) { -// lmax = std::max(lmax, wf.first); -// } -// } -// lmax = std::max(lmax, unit_cell_.lmax()); -// -// if (num_ao__ > 0) { -// mdarray rlm_gk(this->num_gkvec_loc(), Utils::lmmax(lmax)); -// mdarray, 1> idx_gk(this->num_gkvec_loc()); -// #pragma omp parallel for schedule(static) -// for (int igk_loc = 0; igk_loc < this->num_gkvec_loc(); igk_loc++) { -// int igk = this->idxgk(igk_loc); -// /* vs = {r, theta, phi} */ -// auto vs = SHT::spherical_coordinates(this->gkvec().gkvec_cart(igk)); -// /* compute real spherical harmonics for G+k vector */ -// std::vector rlm(Utils::lmmax(lmax)); -// SHT::spherical_harmonics(lmax, vs[1], vs[2], &rlm[0]); -// for (int lm = 0; lm < Utils::lmmax(lmax); lm++) { -// rlm_gk(igk_loc, lm) = rlm[lm]; -// } -// idx_gk(igk_loc) = ctx_.atomic_wf_ri().iqdq(vs[0]); -// } -// -// /* starting index of atomic orbital block for each atom */ -// std::vector idxao; -// int n{0}; -// for (int ia = 0; ia < unit_cell_.num_atoms(); ia++) { -// auto& atom_type = unit_cell_.atom(ia).type(); -// idxao.push_back(n); -// /* increment index of atomic orbitals */ -// for (auto e : atom_type.pp_desc().atomic_pseudo_wfs_) { -// int l = e.first; -// n += (2 * l + 1); -// // } -// -// #pragma omp parallel for schedule(static) -// for (int ia = 0; ia < unit_cell_.num_atoms(); ia++) { -// double phase = twopi * geometry3d::dot(this->gkvec().vk(), unit_cell_.atom(ia).position()); -// double_complex phase_k = std::exp(double_complex(0.0, phase)); -// -// std::vector phase_gk(this->num_gkvec_loc()); -// for (int igk_loc = 0; igk_loc < this->num_gkvec_loc(); igk_loc++) { -// int igk = this->idxgk(igk_loc); -// auto G = this->gkvec().gvec(igk); -// phase_gk[igk_loc] = std::conj(ctx_.gvec_phase_factor(G, ia) * phase_k); -// } -// auto& atom_type = unit_cell_.atom(ia).type(); -// int n{0}; -// for (int i = 0; i < static_cast(atom_type.pp_desc().atomic_pseudo_wfs_.size()); i++) { -// int l = atom_type.pp_desc().atomic_pseudo_wfs_[i].first; -// double_complex z = std::pow(double_complex(0, -1), l) * fourpi / std::sqrt(unit_cell_.omega()); -// for (int m = -l; m <= l; m++) { -// int lm = Utils::lm_by_l_m(l, m); -// for (int igk_loc = 0; igk_loc < this->num_gkvec_loc(); igk_loc++) { -// phi.component(0).pw_coeffs().prime(igk_loc, idxao[ia] + n) = -// z * phase_gk[igk_loc] * rlm_gk(igk_loc, lm) * -// ctx_.atomic_wf_ri().values(i, atom_type.id())(idx_gk[igk_loc].first, -// idx_gk[igk_loc].second); -// } -// n++; -// } -// } -// } -// } -//} diff --git a/src/K_point/initialize.hpp b/src/K_point/initialize.hpp index c98ce4d3a..34580ba27 100644 --- a/src/K_point/initialize.hpp +++ b/src/K_point/initialize.hpp @@ -7,7 +7,7 @@ inline void K_point::initialize() zil_[l] = std::pow(double_complex(0, 1), l); } - l_by_lm_ = Utils::l_by_lm(ctx_.lmax_apw()); + l_by_lm_ = utils::l_by_lm(ctx_.lmax_apw()); int bs = ctx_.cyclic_block_size(); diff --git a/src/Potential/paw_potential.hpp b/src/Potential/paw_potential.hpp index 39dffb82c..7254ad444 100644 --- a/src/Potential/paw_potential.hpp +++ b/src/Potential/paw_potential.hpp @@ -19,7 +19,7 @@ inline void Potential::init_PAW() auto& atom_type = atom.type(); int l_max = 2 * atom_type.indexr().lmax_lo(); - int lm_max_rho = Utils::lmmax(l_max); + int lm_max_rho = utils::lmmax(l_max); paw_potential_data_t ppd; @@ -293,15 +293,13 @@ inline double Potential::xc_mt_PAW_noncollinear(std::vector exc_lm = transform(sht_.get(), exc_tp ); + /* transform to lm- domain */ + Spheric_function exc_lm = transform(sht_.get(), exc_tp); return inner(exc_lm, rho_u_lm + rho_d_lm); } @@ -326,7 +324,7 @@ inline double Potential::calc_PAW_hartree_potential(Atom& atom, /* calculate energy */ - auto l_by_lm = Utils::l_by_lm(utils::lmax(lmsize_rho)); + auto l_by_lm = utils::l_by_lm(utils::lmax(lmsize_rho)); /* create array for integration */ std::vector intdata(grid.num_points(),0); @@ -351,10 +349,7 @@ inline void Potential::calc_PAW_local_potential(paw_potential_data_t &ppd, std::vector> const& ae_density, std::vector> const& ps_density) { - //----------------------------------------- - //---- Calculation of Hartree potential --- - //----------------------------------------- - + /* calculation of Hartree potential */ for (int i = 0; i < ctx_.num_mag_dims() + 1; i++) { ppd.ae_potential_[i].zero(); ppd.ps_potential_[i].zero(); @@ -370,9 +365,7 @@ inline void Potential::calc_PAW_local_potential(paw_potential_data_t &ppd, ppd.hartree_energy_ = ae_hartree_energy - ps_hartree_energy; - //----------------------------------------- - //---- Calculation of XC potential --- - //----------------------------------------- + /* calculation of XC potential */ auto& ps_core = ppd.atom_->type().ps_core_charge_density(); auto& ae_core = ppd.atom_->type().paw_ae_core_charge_density(); @@ -418,9 +411,9 @@ inline void Potential::calc_PAW_local_Dij(paw_potential_data_t& pdd, mdarray GC(lmax, 2 * lmax, lmax, SHT::gaunt_rlm); diff --git a/src/Potential/potential.hpp b/src/Potential/potential.hpp index d1d64b5c2..51fa65fc5 100644 --- a/src/Potential/potential.hpp +++ b/src/Potential/potential.hpp @@ -83,7 +83,7 @@ class Potential : public Field4D std::vector zilm_; - mdarray l_by_lm_; + std::vector l_by_lm_; mdarray gvec_ylm_; @@ -332,7 +332,7 @@ class Potential : public Field4D sht_ = std::unique_ptr(new SHT(lmax_)); if (lmax_ >= 0) { - l_by_lm_ = Utils::l_by_lm(lmax_); + l_by_lm_ = utils::l_by_lm(lmax_); /* precompute i^l */ zil_.resize(lmax_ + 1); @@ -340,7 +340,7 @@ class Potential : public Field4D zil_[l] = std::pow(double_complex(0, 1), l); } - zilm_.resize(Utils::lmmax(lmax_)); + zilm_.resize(utils::lmmax(lmax_)); for (int l = 0, lm = 0; l <= lmax_; l++) { for (int m = -l; m <= l; m++, lm++) { zilm_[lm] = zil_[l]; diff --git a/src/Unit_cell/atom.h b/src/Unit_cell/atom.h index ae1a88432..571820e12 100644 --- a/src/Unit_cell/atom.h +++ b/src/Unit_cell/atom.h @@ -125,7 +125,7 @@ class Atom lmax_pot_ = type().parameters().lmax_pot(); if (type().parameters().full_potential()) { - int lmmax = Utils::lmmax(lmax_pot_); + int lmmax = utils::lmmax(lmax_pot_); int nrf = type().indexr().size(); h_radial_integrals_ = mdarray(lmmax, nrf, nrf); @@ -170,7 +170,7 @@ class Atom { PROFILE("sirius::Atom::generate_radial_integrals"); - int lmmax = Utils::lmmax(lmax_pot_); + int lmmax = utils::lmmax(lmax_pot_); int nmtp = type().num_mt_points(); int nrf = type().indexr().size(); int num_mag_dims = type().parameters().num_mag_dims(); @@ -181,7 +181,7 @@ class Atom splindex spl_lm(lmmax, comm__.size(), comm__.rank()); - auto l_by_lm = Utils::l_by_lm(lmax_pot_); + auto l_by_lm = utils::l_by_lm(lmax_pot_); h_radial_integrals_.zero(); if (num_mag_dims) { @@ -402,9 +402,9 @@ class Atom /// Set muffin-tin potential and magnetic field. inline void set_nonspherical_potential(double* veff__, double* beff__[3]) { - veff_ = mdarray(veff__, Utils::lmmax(lmax_pot_), type().num_mt_points()); + veff_ = mdarray(veff__, utils::lmmax(lmax_pot_), type().num_mt_points()); for (int j = 0; j < 3; j++) { - beff_[j] = mdarray(beff__[j], Utils::lmmax(lmax_pot_), type().num_mt_points()); + beff_[j] = mdarray(beff__[j], utils::lmmax(lmax_pot_), type().num_mt_points()); } } diff --git a/src/Unit_cell/atom_type.h b/src/Unit_cell/atom_type.h index 88f4af057..cae3bb385 100644 --- a/src/Unit_cell/atom_type.h +++ b/src/Unit_cell/atom_type.h @@ -29,7 +29,6 @@ #include "atomic_data.hpp" #include "descriptors.h" #include "geometry3d.hpp" -#include "utils.h" #include "radial_grid.h" #include "radial_solver.h" #include "xc_functional.h" @@ -1222,10 +1221,10 @@ inline void Atom_type::init(int offset_lo__) /* get number of valence electrons */ num_valence_electrons_ = zn_ - num_core_electrons_; - int lmmax_pot = Utils::lmmax(parameters_.lmax_pot()); + int lmmax_pot = utils::lmmax(parameters_.lmax_pot()); if (parameters_.full_potential()) { - auto l_by_lm = Utils::l_by_lm(parameters_.lmax_pot()); + auto l_by_lm = utils::l_by_lm(parameters_.lmax_pot()); /* index the non-zero radial integrals */ std::vector> non_zero_elements; @@ -1716,7 +1715,7 @@ inline void Atom_type::read_pseudo_paw(json const& parser) inline void Atom_type::read_input(std::string const& str__) { - json parser = Utils::read_json_from_file_or_string(str__); + json parser = utils::read_json_from_file_or_string(str__); if (parser.empty()) { return; diff --git a/src/Unit_cell/basis_functions_index.hpp b/src/Unit_cell/basis_functions_index.hpp index 545fb83b4..5010be5c4 100644 --- a/src/Unit_cell/basis_functions_index.hpp +++ b/src/Unit_cell/basis_functions_index.hpp @@ -96,7 +96,7 @@ class basis_functions_index basis_function_index_descriptor(l, m, indexr__[idxrf].j, order, idxlo, idxrf)); } } - index_by_lm_order_ = mdarray(Utils::lmmax(indexr__.lmax()), indexr__.max_num_rf()); + index_by_lm_order_ = mdarray(utils::lmmax(indexr__.lmax()), indexr__.max_num_rf()); for (int i = 0; i < (int)basis_function_index_descriptors_.size(); i++) { int lm = basis_function_index_descriptors_[i].lm; diff --git a/src/Unit_cell/unit_cell_symmetry.h b/src/Unit_cell/unit_cell_symmetry.h index 7044c3c63..deed6b129 100644 --- a/src/Unit_cell/unit_cell_symmetry.h +++ b/src/Unit_cell/unit_cell_symmetry.h @@ -35,7 +35,6 @@ extern "C" { #include "geometry3d.hpp" #include "constants.h" -#include "utils.h" #include "gvec.hpp" namespace sirius { diff --git a/src/augmentation_operator.h b/src/augmentation_operator.h index 0d65ab5b8..f76f2e2bf 100644 --- a/src/augmentation_operator.h +++ b/src/augmentation_operator.h @@ -70,7 +70,7 @@ class Augmentation_operator int lmax_beta = atom_type_.indexr().lmax(); int lmmax = utils::lmmax(2 * lmax_beta); - auto l_by_lm = Utils::l_by_lm(2 * lmax_beta); + auto l_by_lm = utils::l_by_lm(2 * lmax_beta); std::vector zilm(lmmax); for (int l = 0, lm = 0; l <= 2 * lmax_beta; l++) { @@ -87,7 +87,7 @@ class Augmentation_operator int gvec_offset = gvec_.offset(); /* array of real spherical harmonics for each G-vector */ - mdarray gvec_rlm(Utils::lmmax(2 * lmax_beta), gvec_count); + mdarray gvec_rlm(utils::lmmax(2 * lmax_beta), gvec_count); #pragma omp parallel for schedule(static) for (int igloc = 0; igloc < gvec_count; igloc++) { int ig = gvec_offset + igloc; @@ -246,7 +246,7 @@ class Augmentation_operator_gvec_deriv PROFILE("sirius::Augmentation_operator_gvec_deriv|constructor"); int lmax = lmax__; - int lmmax = Utils::lmmax(2 * lmax); + int lmmax = utils::lmmax(2 * lmax); /* Gaunt coefficients of three real spherical harmonics */ gaunt_coefs_ = std::unique_ptr>(new Gaunt_coefficients(lmax, 2 * lmax, lmax, SHT::gaunt_rlm)); @@ -293,9 +293,9 @@ class Augmentation_operator_gvec_deriv /* maximum l of beta-projectors */ int lmax_beta = atom_type__.indexr().lmax(); - int lmmax = Utils::lmmax(2 * lmax_beta); + int lmmax = utils::lmmax(2 * lmax_beta); - auto l_by_lm = Utils::l_by_lm(2 * lmax_beta); + auto l_by_lm = utils::l_by_lm(2 * lmax_beta); std::vector zilm(lmmax); for (int l = 0, lm = 0; l <= 2 * lmax_beta; l++) { diff --git a/src/dft_ground_state.h b/src/dft_ground_state.h index 156d58539..181006019 100644 --- a/src/dft_ground_state.h +++ b/src/dft_ground_state.h @@ -25,7 +25,6 @@ #ifndef __DFT_GROUND_STATE_H__ #define __DFT_GROUND_STATE_H__ -//#include "potential.h" #include "k_point_set.h" #include "utils/json.hpp" #include "Hubbard/hubbard.hpp" diff --git a/src/gaunt.h b/src/gaunt.h index d61029ba2..80adf2cc1 100644 --- a/src/gaunt.h +++ b/src/gaunt.h @@ -26,7 +26,6 @@ #define __GAUNT_H__ #include "mdarray.hpp" -#include "utils.h" namespace sirius { @@ -89,9 +88,9 @@ class Gaunt_coefficients , lmax3_(lmax3__) , lmax2_(lmax2__) { - lmmax1_ = Utils::lmmax(lmax1_); - lmmax3_ = Utils::lmmax(lmax3_); - lmmax2_ = Utils::lmmax(lmax2_); + lmmax1_ = utils::lmmax(lmax1_); + lmmax3_ = utils::lmmax(lmax3_); + lmmax2_ = utils::lmmax(lmax2_); gaunt_packed_L1_L2_ = mdarray>, 1>(lmmax3_); gaunt_L1_L2 g12; diff --git a/src/k_point.h b/src/k_point.h index cbd1ad73b..a609e16b0 100644 --- a/src/k_point.h +++ b/src/k_point.h @@ -144,7 +144,7 @@ class K_point std::vector zil_; /// Mapping between lm and l. - mdarray l_by_lm_; + std::vector l_by_lm_; /// Column rank of the processors of ScaLAPACK/ELPA diagonalization grid. int rank_col_; @@ -296,8 +296,11 @@ class K_point /// Generate two-component spinor wave functions inline void generate_spinor_wave_functions(); + inline void generate_atomic_centered_wavefunctions(const int num_ao__, Wave_functions &phi); - inline void generate_atomic_centered_wavefunctions_(const int num_ao__, Wave_functions &phi, std::vector &offset, bool hubbard); + + inline void generate_atomic_centered_wavefunctions_aux(const int num_ao__, Wave_functions &phi, std::vector &offset, bool hubbard); + void compute_gradient_wavefunctions(Wave_functions &phi, const int starting_position_i, const int num_wf, diff --git a/src/matching_coefficients.h b/src/matching_coefficients.h index a93ed65bc..7a2779ee0 100644 --- a/src/matching_coefficients.h +++ b/src/matching_coefficients.h @@ -162,7 +162,7 @@ class Matching_coefficients igk_(igk__), gkvec_(gkvec__) { - int lmmax_apw = Utils::lmmax(lmax_apw__); + int lmmax_apw = utils::lmmax(lmax_apw__); gkvec_ylm_ = mdarray(num_gkvec_, lmmax_apw); diff --git a/src/radial_integrals.h b/src/radial_integrals.h index a0a5e6414..0f065315b 100644 --- a/src/radial_integrals.h +++ b/src/radial_integrals.h @@ -25,6 +25,7 @@ #ifndef __RADIAL_INTEGRALS_H__ #define __RADIAL_INTEGRALS_H__ +#include #include "Unit_cell/unit_cell.h" #include "sbessel.h" @@ -84,7 +85,7 @@ class Radial_integrals_base /// Radial integrals of the atomic centered orbitals. /** Used in initialize_subspace and in the hubbard correction. */ template - class Radial_integrals_atomic_wf : public Radial_integrals_base<2> +class Radial_integrals_atomic_wf : public Radial_integrals_base<2> { private: @@ -482,48 +483,6 @@ class Radial_integrals_beta_jl : public Radial_integrals_base<3> } }; -/// Radial integrals for the step function of the LAPW method. -/** Radial integrals have the following expression: - * \f[ - * \Theta(\alpha, G) = \int_{0}^{R_{\alpha}} \frac{\sin(Gr)}{Gr} r^2 dr = - * \left\{ \begin{array}{ll} \displaystyle R_{\alpha}^3 / 3 & G=0 \\ - * \Big( \sin(GR_{\alpha}) - GR_{\alpha}\cos(GR_{\alpha}) \Big) / G^3 & G \ne 0 \end{array} \right. - * \f] - */ -//class Radial_integrals_theta : public Radial_integrals_base<1> -//{ -// private: -// void generate() -// { -// PROFILE("sirius::Radial_integrals|theta"); -// -// for (int iat = 0; iat < unit_cell_.num_atom_types(); iat++) { -// auto& atom_type = unit_cell_.atom_type(iat); -// auto R = atom_type.mt_radius(); -// values_(iat) = Spline(grid_q_); -// -// #pragma omp parallel for -// for (int iq = 0; iq < grid_q_.num_points(); iq++) { -// if (iq == 0) { -// values_(iat)[iq] = std::pow(R, 3) / 3.0; -// } else { -// double g = grid_q_[iq]; -// values_(iat)[iq] = (std::sin(g * R) - g * R * std::cos(g * R)) / std::pow(g, 3); -// } -// } -// values_(iat).interpolate(); -// } -// } -// -// public: -// Radial_integrals_theta(Unit_cell const& unit_cell__, double qmax__, int np__) -// : Radial_integrals_base<1>(unit_cell__, qmax__, np__) -// { -// values_ = mdarray, 1>(unit_cell_.num_atom_types()); -// generate(); -// } -//}; - template class Radial_integrals_vloc : public Radial_integrals_base<1> { diff --git a/src/sht.h b/src/sht.h index eccac768d..14a3b9e60 100644 --- a/src/sht.h +++ b/src/sht.h @@ -34,7 +34,6 @@ #include #include "typedefs.h" -#include "utils.h" #include "linalg.hpp" #include "lebedev_grids.hpp" @@ -259,7 +258,7 @@ class SHT // TODO: better name if (m == 0) { f_ylm__[lm] = f_rlm__[lm]; } else { - int lm1 = Utils::lm_by_l_m(l, -m); + int lm1 = utils::lm(l, -m); f_ylm__[lm] = ylm_dot_rlm(l, m, m) * f_rlm__[lm] + ylm_dot_rlm(l, m, -m) * f_rlm__[lm1]; } lm++; @@ -276,7 +275,7 @@ class SHT // TODO: better name if (m == 0) { f_rlm__[lm] = std::real(f_ylm__[lm]); } else { - int lm1 = Utils::lm_by_l_m(l, -m); + int lm1 = utils::lm(l, -m); f_rlm__[lm] = std::real(rlm_dot_ylm(l, m, m) * f_ylm__[lm] + rlm_dot_ylm(l, m, -m) * f_ylm__[lm1]); } lm++; @@ -356,11 +355,11 @@ class SHT // TODO: better name for (int l = 0; l <= lmax; l++) { for (int m = 0; m <= l; m++) { double_complex z = std::exp(double_complex(0.0, m * phi)); - ylm[Utils::lm_by_l_m(l, m)] = gsl_sf_legendre_sphPlm(l, m, x) * z; + ylm[utils::lm(l, m)] = gsl_sf_legendre_sphPlm(l, m, x) * z; if (m % 2) { - ylm[Utils::lm_by_l_m(l, -m)] = -std::conj(ylm[Utils::lm_by_l_m(l, m)]); + ylm[utils::lm(l, -m)] = -std::conj(ylm[utils::lm(l, m)]); } else { - ylm[Utils::lm_by_l_m(l, -m)] = std::conj(ylm[Utils::lm_by_l_m(l, m)]); + ylm[utils::lm(l, -m)] = std::conj(ylm[utils::lm(l, m)]); } } } @@ -387,13 +386,13 @@ class SHT // TODO: better name for (int l = 1; l <= lmax; l++) { for (int m = -l; m < 0; m++) { - rlm[Utils::lm_by_l_m(l, m)] = t * ylm[Utils::lm_by_l_m(l, m)].imag(); + rlm[utils::lm(l, m)] = t * ylm[utils::lm(l, m)].imag(); } - rlm[Utils::lm_by_l_m(l, 0)] = ylm[Utils::lm_by_l_m(l, 0)].real(); + rlm[utils::lm(l, 0)] = ylm[utils::lm(l, 0)].real(); for (int m = 1; m <= l; m++) { - rlm[Utils::lm_by_l_m(l, m)] = t * ylm[Utils::lm_by_l_m(l, m)].real(); + rlm[utils::lm(l, m)] = t * ylm[utils::lm(l, m)].real(); } } } diff --git a/src/simulation_parameters.h b/src/simulation_parameters.h index d829067a3..7321f0f9b 100644 --- a/src/simulation_parameters.h +++ b/src/simulation_parameters.h @@ -26,7 +26,6 @@ #define __SIMULATION_PARAMETERS_H__ #include "typedefs.h" -#include "utils.h" #include "input.h" namespace sirius { @@ -79,7 +78,7 @@ class Simulation_parameters return; } - json dict = Utils::read_json_from_file_or_string(str__); + json dict = utils::read_json_from_file_or_string(str__); /* read unit cell */ unit_cell_input_.read(dict); @@ -269,7 +268,7 @@ class Simulation_parameters inline int lmmax_apw() const { - return Utils::lmmax(parameters_input_.lmax_apw_); + return utils::lmmax(parameters_input_.lmax_apw_); } inline int lmax_rho() const @@ -279,7 +278,7 @@ class Simulation_parameters inline int lmmax_rho() const { - return Utils::lmmax(parameters_input_.lmax_rho_); + return utils::lmmax(parameters_input_.lmax_rho_); } inline int lmax_pot() const @@ -289,7 +288,7 @@ class Simulation_parameters inline int lmmax_pot() const { - return Utils::lmmax(parameters_input_.lmax_pot_); + return utils::lmmax(parameters_input_.lmax_pot_); } inline double aw_cutoff() const diff --git a/src/spheric_function.h b/src/spheric_function.h index 34cf1a94e..c58467b83 100644 --- a/src/spheric_function.h +++ b/src/spheric_function.h @@ -276,7 +276,7 @@ Spheric_function laplacian(Spheric_function const& f__ for (int l = 0; l <= lmax; l++) { int ll = l * (l + 1); for (int m = -l; m <= l; m++) { - int lm = Utils::lm_by_l_m(l, m); + int lm = utils::lm(l, m); /* get lm component */ auto s = f__.component(lm); /* compute 1st derivative */ @@ -304,7 +304,7 @@ inline Spheric_function convert(Spheric_function tpm(f__.angular_domain_size()); for (int l = 0; l <= lmax; l++) { for (int m = -l; m <= l; m++) { - int lm = Utils::lm_by_l_m(l, m); + int lm = utils::lm(l, m); tpp[lm] = SHT::rlm_dot_ylm(l, m, m); tpm[lm] = SHT::rlm_dot_ylm(l, m, -m); } @@ -319,7 +319,7 @@ inline Spheric_function convert(Spheric_function convert(Spheric_function tpm(f__.angular_domain_size()); for (int l = 0; l <= lmax; l++) { for (int m = -l; m <= l; m++) { - int lm = Utils::lm_by_l_m(l, m); + int lm = utils::lm(l, m); tpp[lm] = SHT::ylm_dot_rlm(l, m, m); tpm[lm] = SHT::ylm_dot_rlm(l, m, -m); } @@ -355,7 +355,7 @@ inline Spheric_function convert(Spheric_function gradient(Spheric_func double d2 = sqrt(double(l) / double(2 * l - 1)); for (int m = -l; m <= l; m++) { - int lm = Utils::lm_by_l_m(l, m); + int lm = utils::lm(l, m); auto s = f.component(lm); for (int mu = -1; mu <= 1; mu++) { int j = (mu + 2) % 3; // map -1,0,1 to 1,2,0 if ((l + 1) <= lmax && abs(m + mu) <= l + 1) { - int lm1 = Utils::lm_by_l_m(l + 1, m + mu); + int lm1 = utils::lm(l + 1, m + mu); double d = d1 * SHT::clebsch_gordan(l, 1, l + 1, m, mu, m + mu); for (int ir = 0; ir < f.radial_grid().num_points(); ir++) g[j](lm1, ir) += (s.deriv(1, ir) - f(lm, ir) * f.radial_grid().x_inv(ir) * double(l)) * d; } if ((l - 1) >= 0 && abs(m + mu) <= l - 1) { - int lm1 = Utils::lm_by_l_m(l - 1, m + mu); + int lm1 = utils::lm(l - 1, m + mu); double d = d2 * SHT::clebsch_gordan(l, 1, l - 1, m, mu, m + mu); for (int ir = 0; ir < f.radial_grid().num_points(); ir++) g[j](lm1, ir) -= (s.deriv(1, ir) + f(lm, ir) * f.radial_grid().x_inv(ir) * double(l + 1)) * d; diff --git a/src/utils.h b/src/utils.h deleted file mode 100644 index 4e6f0a365..000000000 --- a/src/utils.h +++ /dev/null @@ -1,292 +0,0 @@ -// Copyright (c) 2013-2018 Anton Kozhevnikov, Thomas Schulthess -// All rights reserved. -// -// Redistribution and use in source and binary forms, with or without modification, are permitted provided that -// the following conditions are met: -// -// 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the -// following disclaimer. -// 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions -// and the following disclaimer in the documentation and/or other materials provided with the distribution. -// -// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED -// WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR -// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR -// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - -/** \file utils.h - * - * \brief Contains definition and partial implementation of sirius::Utils class. - */ - -#ifndef __UTILS_H__ -#define __UTILS_H__ - -#include -#include -#include -#include -#include "typedefs.h" -#include "constants.h" -#include "sddk.hpp" - -using namespace sddk; - -/// Utility class. -class Utils // TODO: namespace utils -{ - public: - /// Maximum number of \f$ \ell, m \f$ combinations for a given \f$ \ell_{max} \f$ - static inline int lmmax(int lmax) - { - return (lmax + 1) * (lmax + 1); - } - - static inline int lm_by_l_m(int l, int m) // TODO lm_by_l_m(l, m) -> lm(l, m) - { - return (l * l + l + m); - } - - //static inline int lmax_by_lmmax(int lmmax__) // TODO: lmax_by_lmmax(lmmax) -> lmax(lmmax) - //{ - // int lmax = int(std::sqrt(double(lmmax__)) + 1e-8) - 1; - // if (lmmax(lmax) != lmmax__) { - // TERMINATE("wrong lmmax"); - // } - // return lmax; - //} - - static void write_matrix(const std::string& fname, - mdarray& matrix, - int nrow, - int ncol, - bool write_upper_only = true, - bool write_abs_only = false, - std::string fmt = "%18.12f") - { - static int icount = 0; - - if (nrow < 0 || nrow > (int)matrix.size(0) || ncol < 0 || ncol > (int)matrix.size(1)) - TERMINATE("wrong number of rows or columns"); - - icount++; - std::stringstream s; - s << icount; - std::string full_name = s.str() + "_" + fname; - - FILE* fout = fopen(full_name.c_str(), "w"); - - for (int icol = 0; icol < ncol; icol++) { - fprintf(fout, "column : %4i\n", icol); - for (int i = 0; i < 80; i++) - fprintf(fout, "-"); - fprintf(fout, "\n"); - if (write_abs_only) { - fprintf(fout, " row, absolute value\n"); - } else { - fprintf(fout, " row, real part, imaginary part, absolute value\n"); - } - for (int i = 0; i < 80; i++) - fprintf(fout, "-"); - fprintf(fout, "\n"); - - int max_row = (write_upper_only) ? std::min(icol, nrow - 1) : (nrow - 1); - for (int j = 0; j <= max_row; j++) { - if (write_abs_only) { - std::string s = "%4i " + fmt + "\n"; - fprintf(fout, s.c_str(), j, abs(matrix(j, icol))); - } else { - fprintf(fout, "%4i %18.12f %18.12f %18.12f\n", j, real(matrix(j, icol)), imag(matrix(j, icol)), - abs(matrix(j, icol))); - } - } - fprintf(fout, "\n"); - } - - fclose(fout); - } - - static void write_matrix(std::string const& fname, bool write_all, mdarray& matrix) - { - static int icount = 0; - - icount++; - std::stringstream s; - s << icount; - std::string full_name = s.str() + "_" + fname; - - FILE* fout = fopen(full_name.c_str(), "w"); - - for (int icol = 0; icol < (int)matrix.size(1); icol++) { - fprintf(fout, "column : %4i\n", icol); - for (int i = 0; i < 80; i++) - fprintf(fout, "-"); - fprintf(fout, "\n"); - fprintf(fout, " row\n"); - for (int i = 0; i < 80; i++) - fprintf(fout, "-"); - fprintf(fout, "\n"); - - int max_row = (write_all) ? ((int)matrix.size(0) - 1) : std::min(icol, (int)matrix.size(0) - 1); - for (int j = 0; j <= max_row; j++) { - fprintf(fout, "%4i %18.12f\n", j, matrix(j, icol)); - } - fprintf(fout, "\n"); - } - - fclose(fout); - } - - static void write_matrix(std::string const& fname, bool write_all, matrix const& mtrx) - { - static int icount = 0; - - icount++; - std::stringstream s; - s << icount; - std::string full_name = s.str() + "_" + fname; - - FILE* fout = fopen(full_name.c_str(), "w"); - - for (int icol = 0; icol < (int)mtrx.size(1); icol++) { - fprintf(fout, "column : %4i\n", icol); - for (int i = 0; i < 80; i++) - fprintf(fout, "-"); - fprintf(fout, "\n"); - fprintf(fout, " row\n"); - for (int i = 0; i < 80; i++) - fprintf(fout, "-"); - fprintf(fout, "\n"); - - int max_row = (write_all) ? ((int)mtrx.size(0) - 1) : std::min(icol, (int)mtrx.size(0) - 1); - for (int j = 0; j <= max_row; j++) { - fprintf(fout, "%4i %18.12f %18.12f\n", j, real(mtrx(j, icol)), imag(mtrx(j, icol))); - } - fprintf(fout, "\n"); - } - - fclose(fout); - } - - template - static void check_hermitian(const std::string& name, matrix const& mtrx, int n = -1) - { - assert(mtrx.size(0) == mtrx.size(1)); - - double maxdiff = 0.0; - int i0 = -1; - int j0 = -1; - - if (n == -1) { - n = static_cast(mtrx.size(0)); - } - - for (int i = 0; i < n; i++) { - for (int j = 0; j < n; j++) { - double diff = std::abs(mtrx(i, j) - type_wrapper::conjugate(mtrx(j, i))); - if (diff > maxdiff) { - maxdiff = diff; - i0 = i; - j0 = j; - } - } - } - - if (maxdiff > 1e-10) { - std::stringstream s; - s << name << " is not a symmetric or hermitian matrix" << std::endl - << " maximum error: i, j : " << i0 << " " << j0 << " diff : " << maxdiff; - - WARNING(s); - } - } - - static double confined_polynomial(double r, double R, int p1, int p2, int dm) - { - double t = 1.0 - std::pow(r / R, 2); - switch (dm) { - case 0: { - return (std::pow(r, p1) * std::pow(t, p2)); - } - case 2: { - return (-4 * p1 * p2 * std::pow(r, p1) * std::pow(t, p2 - 1) / std::pow(R, 2) + - p1 * (p1 - 1) * std::pow(r, p1 - 2) * std::pow(t, p2) + - std::pow(r, p1) * (4 * (p2 - 1) * p2 * std::pow(r, 2) * std::pow(t, p2 - 2) / std::pow(R, 4) - - 2 * p2 * std::pow(t, p2 - 1) / std::pow(R, 2))); - } - default: { - TERMINATE("wrong derivative order"); - return 0.0; - } - } - } - - static mdarray l_by_lm(int lmax) - { - mdarray v(lmmax(lmax)); - for (int l = 0; l <= lmax; l++) { - for (int m = -l; m <= l; m++) { - v[lm_by_l_m(l, m)] = l; - } - } - return std::move(v); - } - - static std::vector> l_m_by_lm(int lmax) - { - std::vector> v(lmmax(lmax)); - for (int l = 0; l <= lmax; l++) { - for (int m = -l; m <= l; m++) { - int lm = lm_by_l_m(l, m); - v[lm].first = l; - v[lm].second = m; - } - } - return std::move(v); - } - - /// Read json dictionary from file or string. - /** Terminate if file doesn't exist. */ - inline static json read_json_from_file_or_string(std::string const& str__) - { - json dict = {}; - if (str__.size() == 0) { - return std::move(dict); - } - - if (str__.find("{") == std::string::npos) { /* this is a file */ - if (utils::file_exists(str__)) { - try { - std::ifstream(str__) >> dict; - } catch(std::exception& e) { - std::stringstream s; - s << "wrong input json file" << std::endl - << e.what(); - TERMINATE(s); - } - } - else { - std::stringstream s; - s << "file " << str__ << " doesn't exist"; - TERMINATE(s); - } - } else { /* this is a json string */ - try { - std::istringstream(str__) >> dict; - } catch (std::exception& e) { - std::stringstream s; - s << "wrong input json string" << std::endl - << e.what(); - TERMINATE(s); - } - } - - return std::move(dict); - } - -}; - -#endif diff --git a/src/utils/utils.hpp b/src/utils/utils.hpp index a348fa9c3..24ecf54ab 100644 --- a/src/utils/utils.hpp +++ b/src/utils/utils.hpp @@ -105,6 +105,17 @@ inline int lmax(int lmmax__) return lmax; } +inline std::vector l_by_lm(int lmax) +{ + std::vector v(lmmax(lmax)); + for (int l = 0; l <= lmax; l++) { + for (int m = -l; m <= l; m++) { + v[lm(l, m)] = l; + } + } + return std::move(v); +} + inline bool file_exists(std::string file_name) { std::ifstream ifs(file_name.c_str()); @@ -252,236 +263,66 @@ inline uint32_t rand() return a; } -//== -//== static void write_matrix(const std::string& fname, -//== mdarray& matrix, -//== int nrow, -//== int ncol, -//== bool write_upper_only = true, -//== bool write_abs_only = false, -//== std::string fmt = "%18.12f") -//== { -//== static int icount = 0; -//== -//== if (nrow < 0 || nrow > (int)matrix.size(0) || ncol < 0 || ncol > (int)matrix.size(1)) -//== TERMINATE("wrong number of rows or columns"); -//== -//== icount++; -//== std::stringstream s; -//== s << icount; -//== std::string full_name = s.str() + "_" + fname; -//== -//== FILE* fout = fopen(full_name.c_str(), "w"); -//== -//== for (int icol = 0; icol < ncol; icol++) { -//== fprintf(fout, "column : %4i\n", icol); -//== for (int i = 0; i < 80; i++) -//== fprintf(fout, "-"); -//== fprintf(fout, "\n"); -//== if (write_abs_only) { -//== fprintf(fout, " row, absolute value\n"); -//== } else { -//== fprintf(fout, " row, real part, imaginary part, absolute value\n"); -//== } -//== for (int i = 0; i < 80; i++) -//== fprintf(fout, "-"); -//== fprintf(fout, "\n"); -//== -//== int max_row = (write_upper_only) ? std::min(icol, nrow - 1) : (nrow - 1); -//== for (int j = 0; j <= max_row; j++) { -//== if (write_abs_only) { -//== std::string s = "%4i " + fmt + "\n"; -//== fprintf(fout, s.c_str(), j, abs(matrix(j, icol))); -//== } else { -//== fprintf(fout, "%4i %18.12f %18.12f %18.12f\n", j, real(matrix(j, icol)), imag(matrix(j, icol)), -//== abs(matrix(j, icol))); -//== } -//== } -//== fprintf(fout, "\n"); -//== } -//== -//== fclose(fout); -//== } -//== -//== static void write_matrix(std::string const& fname, bool write_all, mdarray& matrix) -//== { -//== static int icount = 0; -//== -//== icount++; -//== std::stringstream s; -//== s << icount; -//== std::string full_name = s.str() + "_" + fname; -//== -//== FILE* fout = fopen(full_name.c_str(), "w"); -//== -//== for (int icol = 0; icol < (int)matrix.size(1); icol++) { -//== fprintf(fout, "column : %4i\n", icol); -//== for (int i = 0; i < 80; i++) -//== fprintf(fout, "-"); -//== fprintf(fout, "\n"); -//== fprintf(fout, " row\n"); -//== for (int i = 0; i < 80; i++) -//== fprintf(fout, "-"); -//== fprintf(fout, "\n"); -//== -//== int max_row = (write_all) ? ((int)matrix.size(0) - 1) : std::min(icol, (int)matrix.size(0) - 1); -//== for (int j = 0; j <= max_row; j++) { -//== fprintf(fout, "%4i %18.12f\n", j, matrix(j, icol)); -//== } -//== fprintf(fout, "\n"); -//== } -//== -//== fclose(fout); -//== } -//== -//== static void write_matrix(std::string const& fname, bool write_all, matrix const& mtrx) -//== { -//== static int icount = 0; -//== -//== icount++; -//== std::stringstream s; -//== s << icount; -//== std::string full_name = s.str() + "_" + fname; -//== -//== FILE* fout = fopen(full_name.c_str(), "w"); -//== -//== for (int icol = 0; icol < (int)mtrx.size(1); icol++) { -//== fprintf(fout, "column : %4i\n", icol); -//== for (int i = 0; i < 80; i++) -//== fprintf(fout, "-"); -//== fprintf(fout, "\n"); -//== fprintf(fout, " row\n"); -//== for (int i = 0; i < 80; i++) -//== fprintf(fout, "-"); -//== fprintf(fout, "\n"); -//== -//== int max_row = (write_all) ? ((int)mtrx.size(0) - 1) : std::min(icol, (int)mtrx.size(0) - 1); -//== for (int j = 0; j <= max_row; j++) { -//== fprintf(fout, "%4i %18.12f %18.12f\n", j, real(mtrx(j, icol)), imag(mtrx(j, icol))); -//== } -//== fprintf(fout, "\n"); -//== } -//== -//== fclose(fout); -//== } -//== -//== template -//== static void check_hermitian(const std::string& name, matrix const& mtrx, int n = -1) -//== { -//== assert(mtrx.size(0) == mtrx.size(1)); -//== -//== double maxdiff = 0.0; -//== int i0 = -1; -//== int j0 = -1; -//== -//== if (n == -1) { -//== n = static_cast(mtrx.size(0)); -//== } -//== -//== for (int i = 0; i < n; i++) { -//== for (int j = 0; j < n; j++) { -//== double diff = std::abs(mtrx(i, j) - type_wrapper::conjugate(mtrx(j, i))); -//== if (diff > maxdiff) { -//== maxdiff = diff; -//== i0 = i; -//== j0 = j; -//== } -//== } -//== } -//== -//== if (maxdiff > 1e-10) { -//== std::stringstream s; -//== s << name << " is not a symmetric or hermitian matrix" << std::endl -//== << " maximum error: i, j : " << i0 << " " << j0 << " diff : " << maxdiff; -//== -//== WARNING(s); -//== } -//== } -//== -//== static double confined_polynomial(double r, double R, int p1, int p2, int dm) -//== { -//== double t = 1.0 - std::pow(r / R, 2); -//== switch (dm) { -//== case 0: { -//== return (std::pow(r, p1) * std::pow(t, p2)); -//== } -//== case 2: { -//== return (-4 * p1 * p2 * std::pow(r, p1) * std::pow(t, p2 - 1) / std::pow(R, 2) + -//== p1 * (p1 - 1) * std::pow(r, p1 - 2) * std::pow(t, p2) + -//== std::pow(r, p1) * (4 * (p2 - 1) * p2 * std::pow(r, 2) * std::pow(t, p2 - 2) / std::pow(R, 4) - -//== 2 * p2 * std::pow(t, p2 - 1) / std::pow(R, 2))); -//== } -//== default: { -//== TERMINATE("wrong derivative order"); -//== return 0.0; -//== } -//== } -//== } -//== -//== static mdarray l_by_lm(int lmax) -//== { -//== mdarray v(lmmax(lmax)); -//== for (int l = 0; l <= lmax; l++) { -//== for (int m = -l; m <= l; m++) { -//== v[lm_by_l_m(l, m)] = l; -//== } -//== } -//== return std::move(v); -//== } -//== -//== static std::vector> l_m_by_lm(int lmax) -//== { -//== std::vector> v(lmmax(lmax)); -//== for (int l = 0; l <= lmax; l++) { -//== for (int m = -l; m <= l; m++) { -//== int lm = lm_by_l_m(l, m); -//== v[lm].first = l; -//== v[lm].second = m; -//== } -//== } -//== return std::move(v); -//== } -//== -//== /// Read json dictionary from file or string. -//== /** Terminate if file doesn't exist. */ -//== inline static json read_json_from_file_or_string(std::string const& str__) -//== { -//== json dict = {}; -//== if (str__.size() == 0) { -//== return std::move(dict); -//== } -//== -//== if (str__.find("{") == std::string::npos) { /* this is a file */ -//== if (Utils::file_exists(str__)) { -//== try { -//== std::ifstream(str__) >> dict; -//== } catch(std::exception& e) { -//== std::stringstream s; -//== s << "wrong input json file" << std::endl -//== << e.what(); -//== TERMINATE(s); -//== } -//== } -//== else { -//== std::stringstream s; -//== s << "file " << str__ << " doesn't exist"; -//== TERMINATE(s); -//== } -//== } else { /* this is a json string */ -//== try { -//== std::istringstream(str__) >> dict; -//== } catch (std::exception& e) { -//== std::stringstream s; -//== s << "wrong input json string" << std::endl -//== << e.what(); -//== TERMINATE(s); -//== } -//== } -//== -//== return std::move(dict); -//== } -//== +inline double confined_polynomial(double r, double R, int p1, int p2, int dm) +{ + double t = 1.0 - std::pow(r / R, 2); + switch (dm) { + case 0: { + return (std::pow(r, p1) * std::pow(t, p2)); + } + case 2: { + return (-4 * p1 * p2 * std::pow(r, p1) * std::pow(t, p2 - 1) / std::pow(R, 2) + + p1 * (p1 - 1) * std::pow(r, p1 - 2) * std::pow(t, p2) + + std::pow(r, p1) * (4 * (p2 - 1) * p2 * std::pow(r, 2) * std::pow(t, p2 - 2) / std::pow(R, 4) - + 2 * p2 * std::pow(t, p2 - 1) / std::pow(R, 2))); + } + default: { + TERMINATE("wrong derivative order"); + return 0.0; + } + } +} +/// Read json dictionary from file or string. +/** Terminate if file doesn't exist. */ +inline json read_json_from_file_or_string(std::string const& str__) +{ + json dict = {}; + if (str__.size() == 0) { + return std::move(dict); + } + + if (str__.find("{") == std::string::npos) { /* this is a file */ + if (file_exists(str__)) { + try { + std::ifstream(str__) >> dict; + } catch(std::exception& e) { + std::stringstream s; + s << "wrong input json file" << std::endl + << e.what(); + TERMINATE(s); + } + } + else { + std::stringstream s; + s << "file " << str__ << " doesn't exist"; + TERMINATE(s); + } + } else { /* this is a json string */ + try { + std::istringstream(str__) >> dict; + } catch (std::exception& e) { + std::stringstream s; + s << "wrong input json string" << std::endl + << e.what(); + TERMINATE(s); + } + } + + return std::move(dict); +} + +/// Get high water mark and resident space size values of a given process. inline void get_proc_status(size_t* VmHWM__, size_t* VmRSS__) { *VmHWM__ = 0; @@ -525,7 +366,7 @@ inline void get_proc_status(size_t* VmHWM__, size_t* VmRSS__) inline int get_proc_threads() { int num_threads{-1}; - + std::ifstream ifs("/proc/self/status"); if (ifs.is_open()) { std::string str; diff --git a/src/xc_functional.h b/src/xc_functional.h index 6987ee117..6928c4491 100644 --- a/src/xc_functional.h +++ b/src/xc_functional.h @@ -27,7 +27,6 @@ #include #include -#include "utils.h" namespace sirius {