Skip to content
New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

[refactor] valence D-matrix in Potential class #1037

Open
wants to merge 10 commits into
base: develop
Choose a base branch
from
Open

[refactor] valence D-matrix in Potential class #1037

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
a3e77f6
refactoring
toxa81 Dec 10, 2024
6b38b94
refactoring
toxa81 Dec 10, 2024
bf3699a
order of const keyword
toxa81 Dec 11, 2024
3252248
move valence D-matrix to Potential class
toxa81 Dec 11, 2024
16c9b3d
delete commented code
toxa81 Dec 11, 2024
63c2f50
Merge branch 'develop' of github.com:electronic-structure/SIRIUS into…
toxa81 Dec 11, 2024
84e5349
Merge branch 'develop' into refactor/d_mtrx_in_potential
toxa81 Dec 11, 2024
8aefe14
apply format
toxa81 Dec 19, 2024
99d9312
rename file; move some code to cpp
toxa81 Dec 21, 2024
b9e84c1
fix
toxa81 Dec 21, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ set(_SOURCES
"hubbard/hubbard_occupancies_derivatives.cpp"
"hubbard/hubbard_potential_energy.cpp"
"hubbard/hubbard_matrix.cpp"
"potential/generate_d_operator_matrix.cpp"
"potential/generate_d_mtrx.cpp"
"potential/generate_pw_coeffs.cpp"
"potential/paw_potential.cpp"
"potential/poisson.cpp"
Expand Down
2 changes: 1 addition & 1 deletion src/api/sirius_api.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6421,7 +6421,7 @@ sirius_generate_d_operator_matrix(void* const* gs_handler__, int* error_code__)
call_sirius(
[&]() {
auto& gs = get_gs(gs_handler__);
gs.potential().generate_D_operator_matrix();
gs.potential().generate_d_mtrx();
},
error_code__);
}
Expand Down
11 changes: 3 additions & 8 deletions src/context/simulation_context.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -253,7 +253,8 @@ class Simulation_context : public Simulation_parameters
memory_t host_memory_t_{memory_t::none};

/// SPLA library context.
std::shared_ptr<::spla::Context> spla_ctx_{new ::spla::Context{SPLA_PU_HOST}};
/** Context is mutable to allow SPLA library to change it. */
mutable std::shared_ptr<::spla::Context> spla_ctx_{new ::spla::Context{SPLA_PU_HOST}};

std::ostream* output_stream_{nullptr};
std::ofstream output_file_stream_;
Expand Down Expand Up @@ -706,14 +707,8 @@ class Simulation_context : public Simulation_parameters
return fft_coarse_grid_;
}

auto const&
spla_context() const
{
return *spla_ctx_;
}

auto&
spla_context()
spla_context() const
{
return *spla_ctx_;
}
Expand Down
4 changes: 2 additions & 2 deletions src/density/density.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ inline auto
get_rho_up_dn(int num_mag_dims__, double rho__, r3::vector<double> mag__)
{
if (rho__ < 0.0) {
return std::make_pair<double, double>(0.0, 0.0);
return std::make_pair(0.0, 0.0);
}

double mag{0};
Expand All @@ -85,7 +85,7 @@ get_rho_up_dn(int num_mag_dims__, double rho__, r3::vector<double> mag__)
mag = std::min(mag__.length(), rho__);
}

return std::make_pair<double, double>(0.5 * (rho__ + mag), 0.5 * (rho__ - mag));
return std::make_pair(0.5 * (rho__ + mag), 0.5 * (rho__ - mag));
}

/// PAW density storage.
Expand Down
22 changes: 11 additions & 11 deletions src/geometry/force.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,7 +95,7 @@ Force::add_k_point_contribution(K_point<T>& kp__, mdarray<double, 2>& forces__)
mdarray<real_type<F>, 2> f({3, ctx_.unit_cell().num_atoms()});
f.zero();

add_k_point_contribution_nonlocal<T, F>(ctx_, bp_grad, kp__, f);
add_k_point_contribution_nonlocal<T, F>(potential_, bp_grad, kp__, f);

for (int ia = 0; ia < ctx_.unit_cell().num_atoms(); ia++) {
for (int x : {0, 1, 2}) {
Expand Down Expand Up @@ -275,7 +275,7 @@ Force::calc_forces_hubbard()

if (ctx_.hubbard_correction()) {
/* recompute the hubbard potential */
::sirius::generate_potential(density_.occupation_matrix(), potential_.hubbard_potential());
//::sirius::generate_potential(density_.occupation_matrix(), potential_.hubbard_potential());

Q_operator<double> q_op(ctx_);

Expand Down Expand Up @@ -385,9 +385,9 @@ Force::calc_forces_us()
forces_us_ = mdarray<double, 2>({3, ctx_.unit_cell().num_atoms()});
forces_us_.zero();

potential_.fft_transform(-1);
//potential_.fft_transform(-1);

Unit_cell& unit_cell = ctx_.unit_cell();
auto const& unit_cell = ctx_.unit_cell();

double reduce_g_fact = ctx_.gvec().reduced() ? 2.0 : 1.0;

Expand Down Expand Up @@ -484,11 +484,11 @@ Force::calc_forces_scf_corr()
forces_scf_corr_.zero();

/* get main arrays */
auto& dveff = potential_.dveff();
auto const& dveff = potential_.dveff();

Unit_cell& unit_cell = ctx_.unit_cell();
auto const& unit_cell = ctx_.unit_cell();

fft::Gvec const& gvec = ctx_.gvec();
auto const& gvec = ctx_.gvec();

int gvec_count = gvec.count();
int gvec_offset = gvec.offset();
Expand All @@ -499,7 +499,7 @@ Force::calc_forces_scf_corr()

#pragma omp parallel for
for (int ia = 0; ia < unit_cell.num_atoms(); ia++) {
Atom& atom = unit_cell.atom(ia);
auto const& atom = unit_cell.atom(ia);

int iat = atom.type_id();

Expand Down Expand Up @@ -545,9 +545,9 @@ Force::calc_forces_core()
/* transform from real space to reciprocal */
xc_pot.rg().fft_transform(-1);

Unit_cell& unit_cell = ctx_.unit_cell();
auto const& unit_cell = ctx_.unit_cell();

fft::Gvec const& gvecs = ctx_.gvec();
auto const& gvecs = ctx_.gvec();

int gvec_count = gvecs.count();
int gvec_offset = gvecs.offset();
Expand All @@ -557,7 +557,7 @@ Force::calc_forces_core()
/* here the calculations are in lattice vectors space */
#pragma omp parallel for
for (int ia = 0; ia < unit_cell.num_atoms(); ia++) {
Atom& atom = unit_cell.atom(ia);
auto const& atom = unit_cell.atom(ia);
if (atom.type().ps_core_charge_density().empty()) {
continue;
}
Expand Down
2 changes: 1 addition & 1 deletion src/geometry/force.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ class Force
private:
Simulation_context& ctx_;

const Density& density_;
Density const& density_;

Potential& potential_;

Expand Down
42 changes: 22 additions & 20 deletions src/geometry/non_local_functor.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,19 +17,22 @@
#include "function3d/periodic_function.hpp"
#include "k_point/k_point.hpp"
#include "density/augmentation_operator.hpp"
#include "potential/potential.hpp"
#include "beta_projectors/beta_projectors_base.hpp"

namespace sirius {

/** \tparam T Precision type of the wave-functions */
template <typename T, typename F>
void
add_k_point_contribution_nonlocal(Simulation_context& ctx__, Beta_projectors_base<T>& bp_base__, K_point<T>& kp__,
add_k_point_contribution_nonlocal(Potential& potential__, Beta_projectors_base<T>& bp_base__, K_point<T>& kp__,
mdarray<real_type<F>, 2>& collect_res__)
{
PROFILE("sirius::add_k_point_contribution_nonlocal");

auto& uc = ctx__.unit_cell();
auto& ctx = potential__.ctx();

auto& uc = ctx.unit_cell();

if (uc.max_mt_basis_size() == 0) {
return;
Expand All @@ -39,7 +42,7 @@ add_k_point_contribution_nonlocal(Simulation_context& ctx__, Beta_projectors_bas

double main_two_factor{-2};

auto mt = ctx__.processing_unit_memory_t();
auto mt = ctx.processing_unit_memory_t();
auto bp_gen = bp.make_generator(mt);
auto beta_coeffs = bp_gen.prepare();

Expand All @@ -54,25 +57,25 @@ add_k_point_contribution_nonlocal(Simulation_context& ctx__, Beta_projectors_bas
/* store <beta|psi> for spin up and down */
matrix<F> beta_phi_chunks[2];

for (int ispn = 0; ispn < ctx__.num_spins(); ispn++) {
int nbnd = kp__.num_occupied_bands(ispn);
beta_phi_chunks[ispn] = inner_prod_beta<F>(ctx__.spla_context(), mt, ctx__.host_memory_t(),
is_device_memory(mt), beta_coeffs, kp__.spinor_wave_functions(),
wf::spin_index(ispn), wf::band_range(0, nbnd));
for (int ispn = 0; ispn < ctx.num_spins(); ispn++) {
int nbnd = kp__.num_occupied_bands(ispn);
beta_phi_chunks[ispn] =
inner_prod_beta<F>(ctx.spla_context(), mt, ctx.host_memory_t(), is_device_memory(mt), beta_coeffs,
kp__.spinor_wave_functions(), wf::spin_index(ispn), wf::band_range(0, nbnd));
}

for (int x = 0; x < bp_base__.num_comp(); x++) {
/* generate chunk for inner product of beta gradient */
bp_base_gen.generate(beta_coeffs_base, icnk, x);

for (int ispn = 0; ispn < ctx__.num_spins(); ispn++) {
for (int ispn = 0; ispn < ctx.num_spins(); ispn++) {
int spin_factor = (ispn == 0 ? 1 : -1);

int nbnd = kp__.num_occupied_bands(ispn);

/* inner product of beta gradient and WF */
auto bp_base_phi_chunk = inner_prod_beta<F>(
ctx__.spla_context(), mt, ctx__.host_memory_t(), is_device_memory(mt), beta_coeffs_base,
ctx.spla_context(), mt, ctx.host_memory_t(), is_device_memory(mt), beta_coeffs_base,
kp__.spinor_wave_functions(), wf::spin_index(ispn), wf::band_range(0, nbnd));

splindex_block<> spl_nbnd(nbnd, n_blocks(kp__.comm().size()), block_id(kp__.comm().rank()));
Expand All @@ -86,6 +89,8 @@ add_k_point_contribution_nonlocal(Simulation_context& ctx__, Beta_projectors_bas
int nbf = beta_coeffs_base.beta_chunk_->desc_(beta_desc_idx::nbf, ia_chunk);
int iat = uc.atom(ia).type_id();

auto& d_mtrx = potential__.d_mtrx(ia);

if (uc.atom(ia).type().spin_orbit_coupling()) {
RTE_THROW("stress and forces with SO coupling are not upported");
}
Expand Down Expand Up @@ -120,14 +125,14 @@ add_k_point_contribution_nonlocal(Simulation_context& ctx__, Beta_projectors_bas
/* Qij exists only in the case of ultrasoft/PAW */
real_type<F> qij{0};
if (uc.atom(ia).type().augment()) {
qij = ctx__.augmentation_op(iat).q_mtrx(ibf, jbf);
qij = ctx.augmentation_op(iat).q_mtrx(ibf, jbf);
}
std::complex<real_type<F>> dij{0};

/* get non-magnetic or collinear spin parts of dij*/
switch (ctx__.num_spins()) {
switch (ctx.num_spins()) {
case 1: {
dij = uc.atom(ia).d_mtrx(ibf, jbf, 0);
dij = d_mtrx(ibf, jbf, 0);
if (lm1 == lm2) {
dij += uc.atom(ia).type().d_mtrx_ion()(idxrf1, idxrf2);
}
Expand All @@ -136,8 +141,7 @@ add_k_point_contribution_nonlocal(Simulation_context& ctx__, Beta_projectors_bas

case 2: {
/* Dij(00) = dij + dij_Z ; Dij(11) = dij - dij_Z*/
dij = (uc.atom(ia).d_mtrx(ibf, jbf, 0) +
spin_factor * uc.atom(ia).d_mtrx(ibf, jbf, 1));
dij = (d_mtrx(ibf, jbf, 0) + spin_factor * d_mtrx(ibf, jbf, 1));
if (lm1 == lm2) {
dij += uc.atom(ia).type().d_mtrx_ion()(idxrf1, idxrf2);
}
Expand All @@ -154,10 +158,10 @@ add_k_point_contribution_nonlocal(Simulation_context& ctx__, Beta_projectors_bas
for_bnd(ibf, jbf, dij, qij, beta_phi_chunks[ispn]);

/* for non-collinear case*/
if (ctx__.num_mag_dims() == 3) {
if (ctx.num_mag_dims() == 3) {
/* Dij(10) = dij_X + i dij_Y ; Dij(01) = dij_X - i dij_Y */
dij = std::complex<real_type<F>>(uc.atom(ia).d_mtrx(ibf, jbf, 2),
spin_factor * uc.atom(ia).d_mtrx(ibf, jbf, 3));
dij = std::complex<real_type<F>>(d_mtrx(ibf, jbf, 2),
spin_factor * d_mtrx(ibf, jbf, 3));
/* add non-diagonal spin components*/
for_bnd(ibf, jbf, dij, 0.0, beta_phi_chunks[ispn + spin_factor]);
}
Expand All @@ -167,8 +171,6 @@ add_k_point_contribution_nonlocal(Simulation_context& ctx__, Beta_projectors_bas
} // ispn
} // x
}

// bp_base__.dismiss();
}

} // namespace sirius
Expand Down
2 changes: 1 addition & 1 deletion src/geometry/stress.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ Stress::calc_stress_nonloc_aux()
auto mg = kp->spinor_wave_functions().memory_guard(mem, wf::copy_to::device);
Beta_projectors_strain_deriv<T> bp_strain_deriv(ctx_, kp->gkvec());

add_k_point_contribution_nonlocal<T, F>(ctx_, bp_strain_deriv, *kp, collect_result);
add_k_point_contribution_nonlocal<T, F>(potential_, bp_strain_deriv, *kp, collect_result);
}

#pragma omp parallel
Expand Down
2 changes: 0 additions & 2 deletions src/hamiltonian/diagonalize_pp.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -146,8 +146,6 @@ diagonalize_pp_exact(int ispn__, Hamiltonian_k<T> const& Hk__, K_point<T>& kp__)
}
} // i (atoms in chunk)
}
// kp.beta_projectors_row().dismiss();
// kp.beta_projectors_col().dismiss();

if (ctx.cfg().control().verification() >= 1) {
double max_diff = check_hermitian(ovlp, kp__.num_gkvec());
Expand Down
2 changes: 1 addition & 1 deletion src/hamiltonian/hamiltonian.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ Hamiltonian0<T>::Hamiltonian0(Potential& potential__, bool precompute_lapw__, bo
new Local_operator<T>(ctx_, ctx_.spfft_coarse<T>(), ctx_.gvec_coarse_fft_sptr(), &potential__));

if (!ctx_.full_potential()) {
d_op_ = std::unique_ptr<D_operator<T>>(new D_operator<T>(ctx_));
d_op_ = std::unique_ptr<D_operator<T>>(new D_operator<T>(potential__));
q_op_ = std::unique_ptr<Q_operator<T>>(new Q_operator<T>(ctx_));
}
if (ctx_.full_potential()) {
Expand Down
27 changes: 14 additions & 13 deletions src/hamiltonian/non_local_operator.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include "non_local_operator.hpp"
#include "beta_projectors/beta_projectors_base.hpp"
#include "hubbard/hubbard_matrix.hpp"
#include "potential/potential.hpp"

namespace sirius {

Expand Down Expand Up @@ -93,21 +94,21 @@ Non_local_operator<T>::get_matrix(int ispn, memory_t mem) const
}

template <typename T>
D_operator<T>::D_operator(Simulation_context const& ctx_)
: Non_local_operator<T>(ctx_)
D_operator<T>::D_operator(Potential& potential__)
: Non_local_operator<T>(potential__.ctx())
{
if (ctx_.gamma_point()) {
this->op_ = mdarray<T, 3>({1, this->packed_mtrx_size_, ctx_.num_mag_dims() + 1});
if (potential__.ctx().gamma_point()) {
this->op_ = mdarray<T, 3>({1, this->packed_mtrx_size_, this->ctx_.num_mag_dims() + 1});
} else {
this->op_ = mdarray<T, 3>({2, this->packed_mtrx_size_, ctx_.num_mag_dims() + 1});
this->op_ = mdarray<T, 3>({2, this->packed_mtrx_size_, this->ctx_.num_mag_dims() + 1});
}
this->op_.zero();
initialize();
initialize(potential__);
}

template <typename T>
void
D_operator<T>::initialize()
D_operator<T>::initialize(Potential& potential__)
{
PROFILE("sirius::D_operator::initialize");

Expand Down Expand Up @@ -147,7 +148,7 @@ D_operator<T>::initialize()
for (int alpha = 0; alpha < 4; alpha++) {
for (int sigma1 = 0; sigma1 < 2; sigma1++) {
for (int sigma2 = 0; sigma2 < 2; sigma2++) {
result += atom.d_mtrx(xi1p, xi2p, alpha) *
result += potential__.d_mtrx(ia)(xi1p, xi2p, alpha) *
pauli_matrix[alpha][sigma1][sigma2] *
atom.type().f_coefficients(xi1, xi1p, sigma, sigma1) *
atom.type().f_coefficients(xi2p, xi2, sigma2, sigmap);
Expand Down Expand Up @@ -215,8 +216,8 @@ D_operator<T>::initialize()
int idx = xi2 * nbf + xi1;
switch (this->ctx_.num_mag_dims()) {
case 3: {
T bx = uc.atom(ia).d_mtrx(xi1, xi2, 2);
T by = uc.atom(ia).d_mtrx(xi1, xi2, 3);
T bx = potential__.d_mtrx(ia)(xi1, xi2, 2);
T by = potential__.d_mtrx(ia)(xi1, xi2, 3);

this->op_(0, this->packed_mtrx_offset_(ia) + idx, 2) = bx;
this->op_(1, this->packed_mtrx_offset_(ia) + idx, 2) = -by;
Expand All @@ -225,8 +226,8 @@ D_operator<T>::initialize()
this->op_(1, this->packed_mtrx_offset_(ia) + idx, 3) = by;
}
case 1: {
T v = uc.atom(ia).d_mtrx(xi1, xi2, 0);
T bz = uc.atom(ia).d_mtrx(xi1, xi2, 1);
T v = potential__.d_mtrx(ia)(xi1, xi2, 0);
T bz = potential__.d_mtrx(ia)(xi1, xi2, 1);

/* add ionic part */
if (lm1 == lm2) {
Expand All @@ -238,7 +239,7 @@ D_operator<T>::initialize()
break;
}
case 0: {
this->op_(0, this->packed_mtrx_offset_(ia) + idx, 0) = uc.atom(ia).d_mtrx(xi1, xi2, 0);
this->op_(0, this->packed_mtrx_offset_(ia) + idx, 0) = potential__.d_mtrx(ia)(xi1, xi2, 0);
/* add ionic part */
if (lm1 == lm2) {
this->op_(0, this->packed_mtrx_offset_(ia) + idx, 0) += dion(idxrf1, idxrf2);
Expand Down
Loading
Loading