Merge pull request #567 from electronic-structure/fix/ldapu

Bugfix in LDA+U occupancy matrix generation

apps/dft_loop/sirius.scf.cpp

@@ -337,7 +337,7 @@ void run_tasks(cmd_args const& args)
             band.initialize_subspace(ks, H0);
             if (ctx->hubbard_correction()) {
                 TERMINATE("fix me");
-                potential.U().hubbard_compute_occupation_numbers(ks); // TODO: this is wrong; U matrix should come form the saved file
+                potential.U().compute_occupation_matrix(ks); // TODO: this is wrong; U matrix should come form the saved file
                 potential.U().calculate_hubbard_potential_and_energy();
             }
         }

src/SDDK/dmatrix.cpp

@@ -75,6 +75,19 @@ dmatrix<T>::dmatrix(int num_rows__, int num_cols__, memory_t mem_type__)
 {
 }
 
+template <typename T>
+dmatrix<T>::dmatrix(int num_rows__, int num_cols__, memory_pool& mp__, std::string const& label__)
+    : matrix<T>(num_rows__, num_cols__, mp__, label__)
+    , num_rows_(num_rows__)
+    , num_cols_(num_cols__)
+    , bs_row_(1)
+    , bs_col_(1)
+    , spl_row_(num_rows_, 1, 0, bs_row_)
+    , spl_col_(num_cols_, 1, 0, bs_col_)
+{
+}
+
+
 template <typename T>
 dmatrix<T>::dmatrix(T* ptr__, int num_rows__, int num_cols__)
     : matrix<T>(ptr__, num_rows__, num_cols__)

src/SDDK/dmatrix.hpp

@@ -88,6 +88,8 @@ class dmatrix : public matrix<T>
 
     dmatrix(int num_rows__, int num_cols__, memory_t mem_type__ = memory_t::host);
 
+    dmatrix(int num_rows__, int num_cols__, memory_pool& mp__, std::string const& label__ = "");
+
     dmatrix(T* ptr__, int num_rows__, int num_cols__, BLACS_grid const& blacs_grid__, int bs_row__, int bs_col__);
 
     dmatrix(int num_rows__, int num_cols__, BLACS_grid const& blacs_grid__, int bs_row__, int bs_col__,

src/api/sirius_api.cpp

@@ -3418,7 +3418,7 @@ void sirius_get_wave_functions(void*          const* ks_handler__,
 void sirius_calculate_hubbard_occupancies(void* const* handler__)
 {
     auto& gs = get_gs(handler__);
-    gs.potential().U().hubbard_compute_occupation_numbers(gs.k_point_set());
+    gs.potential().U().compute_occupation_matrix(gs.k_point_set());
 }
 
 
@@ -3441,12 +3441,10 @@ void sirius_calculate_hubbard_occupancies(void* const* handler__)
       doc: Leading dimensions of the occupation matrix.
 @api end
 */
-void sirius_set_hubbard_occupancies(void* const* handler__,
-                                    std::complex<double>* occ__,
-                                    int   const *ld__)
+void sirius_set_hubbard_occupancies(void* const* handler__, std::complex<double>* occ__, int const *ld__)
 {
     auto& gs = get_gs(handler__);
-    gs.potential().U().access_hubbard_occupancies("set", occ__, ld__);
+    gs.potential().U().access_hubbard_occupancies("set", occ__, *ld__);
 }
 
 /*
@@ -3468,12 +3466,10 @@ void sirius_set_hubbard_occupancies(void* const* handler__,
       doc: Leading dimensions of the occupation matrix.
 @api end
 */
-void sirius_get_hubbard_occupancies(void* const* handler__,
-                                    std::complex<double>*      occ__,
-                                    int   const *ld__)
+void sirius_get_hubbard_occupancies(void* const* handler__, std::complex<double>* occ__, int const *ld__)
 {
     auto& gs = get_gs(handler__);
-    gs.potential().U().access_hubbard_occupancies("get", occ__, ld__);
+    gs.potential().U().access_hubbard_occupancies("get", occ__, *ld__);
 }
 
 /*
@@ -3495,12 +3491,10 @@ void sirius_get_hubbard_occupancies(void* const* handler__,
       doc: Leading dimensions of the matrix.
 @api end
 */
-void sirius_set_hubbard_potential(void* const* handler__,
-                                  std::complex<double>*      pot__,
-                                  int   const *ld__)
+void sirius_set_hubbard_potential(void* const* handler__, std::complex<double>* pot__, int const *ld__)
 {
     auto& gs = get_gs(handler__);
-    gs.potential().U().access_hubbard_potential("set", pot__, ld__);
+    gs.potential().U().access_hubbard_potential("set", pot__, *ld__);
 }
 
 
@@ -3523,12 +3517,10 @@ void sirius_set_hubbard_potential(void* const* handler__,
       doc: Leading dimensions of the matrix.
 @api end
 */
-void sirius_get_hubbard_potential(void* const* handler__,
-                                  std::complex<double>*      pot__,
-                                  int   const *ld__)
+void sirius_get_hubbard_potential(void* const* handler__, std::complex<double>* pot__, int const *ld__)
 {
     auto& gs = get_gs(handler__);
-    gs.potential().U().access_hubbard_potential("get", pot__, ld__);
+    gs.potential().U().access_hubbard_potential("get", pot__, *ld__);
 }
 
 /*

src/dft/dft_ground_state.cpp

@@ -195,7 +195,7 @@ json DFT_ground_state::find(double rms_tol, double energy_tol, double initial_to
     std::vector<double> etot_hist;
 
     if (ctx_.hubbard_correction()) { // TODO: move to inititialization functions
-        potential_.U().hubbard_compute_occupation_numbers(kset_);
+        potential_.U().compute_occupation_matrix(kset_);
         potential_.U().calculate_hubbard_potential_and_energy();
     }
 
@@ -285,7 +285,7 @@ json DFT_ground_state::find(double rms_tol, double energy_tol, double initial_to
 
         /* Compute the hubbard correction */
         if (ctx_.hubbard_correction()) {
-            potential_.U().hubbard_compute_occupation_numbers(kset_);
+            potential_.U().compute_occupation_matrix(kset_);
             potential_.U().calculate_hubbard_potential_and_energy();
         }
 

src/hubbard/hubbard.cpp

@@ -41,10 +41,12 @@ Hubbard::Hubbard(Simulation_context& ctx__)
         approximation_ = false;
     }
 
-    occupancy_number_  = mdarray<double_complex, 4>(indexb_max, indexb_max, 4, ctx_.unit_cell().num_atoms());
-    hubbard_potential_ = mdarray<double_complex, 4>(indexb_max, indexb_max, 4, ctx_.unit_cell().num_atoms());
+    occupation_matrix_ = sddk::mdarray<double_complex, 4>(indexb_max, indexb_max, 4, ctx_.unit_cell().num_atoms(),
+            memory_t::host, "occupation_matrix_");
+    hubbard_potential_ = sddk::mdarray<double_complex, 4>(indexb_max, indexb_max, 4, ctx_.unit_cell().num_atoms(),
+            memory_t::host, "hubbard_potential_");
 
-    auto r  = ctx_.unit_cell().num_wf_with_U();
+    auto r = ctx_.unit_cell().num_wf_with_U();
 
     number_of_hubbard_orbitals_ = r.first;
     offset_ = r.second;

src/hubbard/hubbard.hpp

@@ -54,7 +54,7 @@ class Hubbard
 
     int number_of_hubbard_orbitals_{0};
 
-    mdarray<double_complex, 4> occupancy_number_;
+    sddk::mdarray<double_complex, 4> occupation_matrix_;
 
     double hubbard_energy_{0.0};
     double hubbard_energy_u_{0.0};
@@ -166,14 +166,9 @@ class Hubbard
                                  Wave_functions& phi,
                                  Wave_functions& ophi);
 
-    /// Generate the atomic orbitals.
-    //void generate_atomic_orbitals(K_point& kp, Q_operator& q_op);
+    void compute_occupation_matrix(K_point_set& kset_);
 
-    void hubbard_compute_occupation_numbers(K_point_set& kset_);
-
-    void compute_occupancies_derivatives(K_point&                    kp,
-                                         Q_operator& q_op,
-                                         mdarray<double_complex, 6>& dn);
+    void compute_occupancies_derivatives(K_point& kp, Q_operator& q_op, mdarray<double_complex, 6>& dn);
 
     /// Compute derivatives of the occupancy matrix w.r.t.atomic displacement.
     /** \param [in]  kp   K-point.
@@ -216,22 +211,17 @@ class Hubbard
 
     mdarray<double_complex, 4>& occupation_matrix()
     {
-        return occupancy_number_;
+        return occupation_matrix_;
     }
 
     mdarray<double_complex, 4>& potential_matrix()
     {
         return hubbard_potential_;
     }
 
-    void access_hubbard_potential(char const*     what,
-                                  double_complex* occ,
-                                  int const*      ld);
-
-    void access_hubbard_occupancies(char const*     what,
-                                    double_complex* occ,
-                                    int const*      ld);
+    void access_hubbard_potential(std::string const& what, double_complex* occ, int ld);
 
+    void access_hubbard_occupancies(std::string const& what, double_complex* occ, int ld);
 };
 
 } // namespace sirius