[refactor] core density is computed in Density class (#1041)

Move generation of core density to Density class itself.

src/api/sirius_api.cpp

@@ -3037,7 +3037,7 @@ sirius_get_energy(void* const* gs_handler__, char const* label__, double* energy
 
                 std::map<std::string, std::function<double()>> func = {
                         {"total", [&]() { return sirius::total_energy(ctx, kset, density, potential); }},
-                        {"evalsum", [&]() { return sirius::eval_sum(unit_cell, kset); }},
+                        {"evalsum", [&]() { return sirius::eval_sum(density, kset); }},
                         {"exc", [&]() { return sirius::energy_exc(density, potential); }},
                         {"vxc", [&]() { return sirius::energy_vxc(density, potential); }},
                         {"bxc", [&]() { return sirius::energy_bxc(density, potential); }},

src/density/density.cpp

@@ -26,6 +26,128 @@
 
 namespace sirius {
 
+/// Generate all-electron core charge densiy.
+/** Solve “free” atom radial equation for core states and find total core charge density. For that, radial grid 
+ *  of the muffin-tin is extended to effective infinity and spherical potential outside muffin-tin is approximated
+ *  with alpha/r + beta tail. Bound states are found and accumulated in the core charge density. */
+static auto
+generate_core_charge_density(Atom_type const& atom_type__, relativity_t core_rel__, std::vector<double> const& vs__,
+                             std::vector<double>& rho_core__)
+{
+    std::fill(rho_core__.begin(), rho_core__.end(), 0.0);
+
+    /* nothing to do */
+    if (atom_type__.num_core_electrons() == 0.0) {
+        return std::make_pair(0.0, 0.0);
+    }
+
+    int nmtp = atom_type__.num_mt_points();
+
+    std::vector<double> free_atom_grid(nmtp);
+    for (int i = 0; i < nmtp; i++) {
+        free_atom_grid[i] = atom_type__.radial_grid(i);
+    }
+
+    int zn = atom_type__.zn();
+
+    /* extend radial grid */
+    double x  = atom_type__.radial_grid(nmtp - 1);
+    double dx = atom_type__.radial_grid().dx(nmtp - 2);
+    while (x < 10.0 + zn / 10.0) {
+        x += dx;
+        free_atom_grid.push_back(x);
+        dx *= 1.025;
+    }
+    Radial_grid_ext<double> rgrid(static_cast<int>(free_atom_grid.size()), free_atom_grid.data());
+
+    /* interpolate spherical potential inside muffin-tin */
+    Spline<double> svmt(atom_type__.radial_grid());
+    /* remove nucleus contribution from Vmt */
+    for (int ir = 0; ir < nmtp; ir++) {
+        svmt(ir) = vs__[ir] + zn * atom_type__.radial_grid().x_inv(ir);
+    }
+    svmt.interpolate();
+    /* fit tail to alpha/r + beta */
+    double alpha = -(std::pow(atom_type__.mt_radius(), 2) * svmt.deriv(1, nmtp - 1) + zn);
+    double beta  = svmt(nmtp - 1) - (zn + alpha) / atom_type__.mt_radius();
+
+    /* cook an effective potential from muffin-tin part and a tail */
+    std::vector<double> veff(rgrid.num_points());
+    for (int ir = 0; ir < nmtp; ir++) {
+        veff[ir] = vs__[ir];
+    }
+    /* simple tail alpha/r + beta */
+    for (int ir = nmtp; ir < rgrid.num_points(); ir++) {
+        veff[ir] = alpha * rgrid.x_inv(ir) + beta;
+    }
+
+    //== /* write spherical potential */
+    //== std::stringstream sstr;
+    //== sstr << "spheric_potential_" << id_ << ".dat";
+    //== FILE* fout = fopen(sstr.str().c_str(), "w");
+
+    //== for (int ir = 0; ir < rgrid.num_points(); ir++)
+    //== {
+    //==     fprintf(fout, "%18.10f %18.10f\n", rgrid[ir], veff[ir]);
+    //== }
+    //== fclose(fout);
+
+    /* atomic level energies */
+    std::vector<double> level_energy(atom_type__.num_atomic_levels());
+
+    for (int ist = 0; ist < atom_type__.num_atomic_levels(); ist++) {
+        level_energy[ist] = -1.0 * zn / 2 / std::pow(double(atom_type__.atomic_level(ist).n), 2);
+    }
+
+    mdarray<double, 2> rho_t({rgrid.num_points(), atom_type__.num_atomic_levels()});
+    rho_t.zero();
+    #pragma omp parallel for
+    for (int ist = 0; ist < atom_type__.num_atomic_levels(); ist++) {
+        if (atom_type__.atomic_level(ist).core) {
+            /* serch for the bound state */
+            Bound_state bs(core_rel__, zn, atom_type__.atomic_level(ist).n, atom_type__.atomic_level(ist).l,
+                           atom_type__.atomic_level(ist).k, rgrid, veff, level_energy[ist]);
+
+            auto& rho = bs.rho();
+            for (int i = 0; i < rgrid.num_points(); i++) {
+                rho_t(i, ist) = atom_type__.atomic_level(ist).occupancy * rho(i) / fourpi;
+            }
+
+            level_energy[ist] = bs.enu();
+        }
+    }
+
+    /* charge density */
+    Spline<double> rho(rgrid);
+
+    for (int ist = 0; ist < atom_type__.num_atomic_levels(); ist++) {
+        if (atom_type__.atomic_level(ist).core) {
+            for (int i = 0; i < rgrid.num_points(); i++) {
+                rho(i) += rho_t(i, ist);
+            }
+        }
+    }
+
+    for (int ir = 0; ir < atom_type__.num_mt_points(); ir++) {
+        rho_core__[ir] = rho(ir);
+    }
+
+    /* interpolate muffin-tin part of core density */
+    Spline<double> rho_mt(atom_type__.radial_grid(), rho_core__);
+
+    /* compute core leakage */
+    auto core_leakage = fourpi * (rho.interpolate().integrate(2) - rho_mt.integrate(2));
+
+    /* compute eigen-value sum of core states */
+    double core_eval_sum{0};
+    for (int ist = 0; ist < atom_type__.num_atomic_levels(); ist++) {
+        if (atom_type__.atomic_level(ist).core) {
+            core_eval_sum += level_energy[ist] * atom_type__.atomic_level(ist).occupancy;
+        }
+    }
+    return std::make_pair(core_leakage, core_eval_sum);
+}
+
 #if defined(SIRIUS_GPU)
 void
 update_density_rg_1_real_gpu(int size__, float const* psi_rg__, float wt__, float* density_rg__)
@@ -88,6 +210,14 @@ Density::Density(Simulation_context& ctx__)
     /* core density of the pseudopotential method */
     if (!ctx_.full_potential()) {
         rho_pseudo_core_ = std::make_unique<spf>(ctx_.spfft<double>(), ctx_.gvec_fft_sptr());
+    } else {
+        ae_core_charge_density_.resize(unit_cell_.num_atom_symmetry_classes());
+        for (int ic = 0; ic < unit_cell_.num_atom_symmetry_classes(); ic++) {
+            ae_core_charge_density_[ic] =
+                    std::vector<double>(unit_cell_.atom_symmetry_class(ic).atom_type().num_mt_points());
+        }
+        core_eval_sum_.resize(unit_cell_.num_atom_symmetry_classes());
+        core_leakage_.resize(unit_cell_.num_atom_symmetry_classes());
     }
 
     l_by_lm_ = sf::l_by_lm(ctx_.lmax_rho());
@@ -127,8 +257,22 @@ double
 Density::core_leakage() const
 {
     double sum{0};
-    for (int ic = 0; ic < unit_cell_.num_atom_symmetry_classes(); ic++) {
-        sum += core_leakage(ic) * unit_cell_.atom_symmetry_class(ic).num_atoms();
+    if (ctx_.full_potential()) {
+        for (int ic = 0; ic < unit_cell_.num_atom_symmetry_classes(); ic++) {
+            sum += this->core_leakage(ic) * unit_cell_.atom_symmetry_class(ic).num_atoms();
+        }
+    }
+    return sum;
+}
+
+double
+Density::core_eval_sum() const
+{
+    double sum{0};
+    if (ctx_.full_potential()) {
+        for (int ic = 0; ic < unit_cell_.num_atom_symmetry_classes(); ic++) {
+            sum += this->core_eval_sum(ic) * unit_cell_.atom_symmetry_class(ic).num_atoms();
+        }
     }
     return sum;
 }
@@ -1063,12 +1207,11 @@ Density::generate(K_point_set const& ks__, bool symmetrize__, bool add_core__, b
 
     if (ctx_.full_potential()) {
         if (add_core__) {
-            /* find the core states */
-            generate_core_charge_density();
             /* add core contribution */
             for (auto it : unit_cell_.spl_num_atoms()) {
                 for (int ir = 0; ir < unit_cell_.atom(it.i).num_mt_points(); ir++) {
-                    rho().mt()[it.i](0, ir) += unit_cell_.atom(it.i).symmetry_class().ae_core_charge_density(ir) / y00;
+                    int ic = unit_cell_.atom(it.i).symmetry_class().id();
+                    rho().mt()[it.i](0, ir) += ae_core_charge_density_[ic][ir] / y00;
                 }
             }
         }
@@ -1921,16 +2064,15 @@ Density::print_info(std::ostream& out__) const
 
     out__ << "Charges and magnetic moments" << std::endl << hbar(80, '-') << std::endl;
     if (ctx_.full_potential()) {
-        double total_core_leakage{0.0};
+        double total_core_leakage = this->core_leakage();
         out__ << "atom      charge    core leakage";
         if (ctx_.num_mag_dims()) {
             out__ << "                 moment                |moment|";
         }
         out__ << std::endl << hbar(80, '-') << std::endl;
 
         for (int ia = 0; ia < unit_cell_.num_atoms(); ia++) {
-            double core_leakage = unit_cell_.atom(ia).symmetry_class().core_leakage();
-            total_core_leakage += core_leakage;
+            double core_leakage = this->core_leakage(unit_cell_.atom(ia).symmetry_class().id());
             out__ << std::setw(4) << ia << std::setw(12) << std::setprecision(6) << std::fixed << mt_charge[ia]
                   << std::setw(16) << std::setprecision(6) << std::scientific << core_leakage;
             if (ctx_.num_mag_dims()) {
@@ -2055,4 +2197,47 @@ Density::load(std::string name__)
     }
 }
 
+void
+Density::generate_core_charge_density(std::vector<std::vector<double>> const& vs__)
+{
+    if (!ctx_.full_potential()) {
+        return;
+    }
+    PROFILE("sirius::Density::generate_core_charge_density");
+
+    auto& spl_idx = unit_cell_.spl_num_atom_symmetry_classes();
+
+    for (auto it : spl_idx) {
+        auto result          = ::sirius::generate_core_charge_density(unit_cell_.atom_symmetry_class(it.i).atom_type(),
+                                                                      ctx_.core_relativity(), vs__[it.i],
+                                                                      ae_core_charge_density_[it.i]);
+        core_leakage_[it.i]  = result.first;
+        core_eval_sum_[it.i] = result.second;
+    }
+
+    for (auto ic = begin_global(spl_idx); ic != end_global(spl_idx); ic++) {
+        auto rank = spl_idx.location(ic).ib;
+        ctx_.comm().bcast(ae_core_charge_density_[ic].data(),
+                          unit_cell_.atom_symmetry_class(ic).atom_type().num_mt_points(), rank);
+        ctx_.comm().bcast(&core_leakage_[ic], 1, rank);
+        ctx_.comm().bcast(&core_eval_sum_[ic], 1, rank);
+    }
+}
+
+void
+Density::generate_pseudo_core_charge_density()
+{
+    PROFILE("sirius::Density::generate_pseudo_core_charge_density");
+
+    /* get lenghts of all G shells */
+    auto q = ctx_.gvec().shells_len();
+    /* get form-factors for all G shells */
+    auto const ff = ctx_.ri().ps_core_->values(q, ctx_.comm());
+    /* make rho_core(G) */
+    auto v = make_periodic_function<true>(ctx_.unit_cell(), ctx_.gvec(), ctx_.phase_factors_t(), ff);
+
+    std::copy(v.begin(), v.end(), &rho_pseudo_core_->f_pw_local(0));
+    rho_pseudo_core_->fft_transform(1);
+}
+
 } // namespace sirius

src/density/density.hpp

@@ -244,6 +244,17 @@ class Density : public Field4D
                                  Periodic_function<double>, density_matrix_t, PAW_density<double>, Hubbard_matrix>>
             mixer_;
 
+    /// Core charge density.
+    /** All-electron core charge density of the LAPW method. It is recomputed on every SCF iteration due to
+        the change of effective potential. */
+    std::vector<std::vector<double>> ae_core_charge_density_;
+
+    /// Core eigen-value sum.
+    std::vector<double> core_eval_sum_;
+
+    /// Core leakage.
+    std::vector<double> core_leakage_;
+
     /// Generate atomic densities in the case of PAW.
     void
     generate_paw_density(paw_atom_index_t::local iapaw__);
@@ -302,39 +313,9 @@ class Density : public Field4D
     void
     generate_valence_mt();
 
-    /// Generate charge density of core states
-    void
-    generate_core_charge_density()
-    {
-        PROFILE("sirius::Density::generate_core_charge_density");
-
-        auto& spl_idx = unit_cell_.spl_num_atom_symmetry_classes();
-
-        for (auto it : spl_idx) {
-            unit_cell_.atom_symmetry_class(it.i).generate_core_charge_density(ctx_.core_relativity());
-        }
-
-        for (auto ic = begin_global(spl_idx); ic != end_global(spl_idx); ic++) {
-            auto rank = spl_idx.location(ic).ib;
-            unit_cell_.atom_symmetry_class(ic).sync_core_charge_density(ctx_.comm(), rank);
-        }
-    }
-
+    /// Generate pseudo core charge density needed for non-linear core correction.
     void
-    generate_pseudo_core_charge_density()
-    {
-        PROFILE("sirius::Density::generate_pseudo_core_charge_density");
-
-        /* get lenghts of all G shells */
-        auto q = ctx_.gvec().shells_len();
-        /* get form-factors for all G shells */
-        auto const ff = ctx_.ri().ps_core_->values(q, ctx_.comm());
-        /* make rho_core(G) */
-        auto v = make_periodic_function<true>(ctx_.unit_cell(), ctx_.gvec(), ctx_.phase_factors_t(), ff);
-
-        std::copy(v.begin(), v.end(), &rho_pseudo_core_->f_pw_local(0));
-        rho_pseudo_core_->fft_transform(1);
-    }
+    generate_pseudo_core_charge_density();
 
   public:
     /// Constructor
@@ -348,6 +329,10 @@ class Density : public Field4D
     double
     core_leakage() const;
 
+    /// Find total sum of core eigen-values.
+    double
+    core_eval_sum() const;
+
     /// Generate initial charge density and magnetization
     void
     initial_density();
@@ -365,6 +350,10 @@ class Density : public Field4D
     bool
     check_num_electrons() const;
 
+    /// Generate charge density of core states
+    void
+    generate_core_charge_density(std::vector<std::vector<double>> const& vs__);
+
     /// Generate full charge density (valence + core) and magnetization from the wave functions.
     /** This function calls generate_valence() and then in case of full-potential LAPW method adds a core density
         to get the full charge density of the system. Density is generated in spectral representation, i.e.
@@ -414,10 +403,16 @@ class Density : public Field4D
     generate_rho_aug() const;
 
     /// Return core leakage for a specific atom symmetry class
-    inline double
+    inline auto
     core_leakage(int ic) const
     {
-        return unit_cell_.atom_symmetry_class(ic).core_leakage();
+        return core_leakage_[ic];
+    }
+
+    inline auto
+    core_eval_sum(int ic) const
+    {
+        return core_eval_sum_[ic];
     }
 
     /// Return const reference to charge density (scalar functions).

src/dft/dft_ground_state.cpp

@@ -239,6 +239,9 @@ DFT_ground_state::find(double density_tol__, double energy_tol__, double iter_so
                                                          ctx_.cfg().iterative_solver().num_steps());
             /* find band occupancies */
             kset_.find_band_occupancies<double>();
+
+            auto vs = potential_.get_spherical_potential();
+            density_.generate_core_charge_density(vs);
             /* generate new density from the occupied wave-functions */
             density_.generate<double>(kset_, ctx_.use_symmetry(), true, true);
         }
@@ -432,7 +435,7 @@ void
 DFT_ground_state::print_info(std::ostream& out__) const
 {
     double evalsum1     = kset_.valence_eval_sum();
-    double evalsum2     = core_eval_sum(ctx_.unit_cell());
+    double evalsum2     = density_.core_eval_sum();
     double s_sum        = kset_.entropy_sum();
     double ekin         = energy_kin(ctx_, kset_, density_, potential_);
     double evxc         = energy_vxc(density_, potential_);