[fix] API to pass J values for Hubbard correction (#1048)

* Convert *J__ pointer to array * use std::array; change API signature * apply format --------- Co-authored-by: Mathieu Taillefumier <mathieu.taillefumier@free.fr>
electronic-structure · Mar 9, 2025 · 1db4164 · 1db4164
1 parent dfb6bff
commit 1db4164
Show file tree

Hide file tree

Showing 14 changed files with 153 additions and 97 deletions.
diff --git a/apps/unit_tests/test_init_ctx.cpp b/apps/unit_tests/test_init_ctx.cpp
@@ -37,7 +37,7 @@ run_test(cmd_args const& args)
     }
 
     try {
-        ctx.cfg().hubbard().local(0).BE2();
+        ctx.cfg().hubbard().local(0).E2();
     } catch (nlohmann::json::exception const& e) {
     } catch (...) {
         return 4;

diff --git a/src/api/sirius.f90 b/src/api/sirius.f90
@@ -2357,7 +2357,7 @@ subroutine sirius_set_atom_type_hubbard(handler,label,l,n,occ,U,J,alpha,beta,J0,
 integer, target, intent(in) :: n
 real(8), target, intent(in) :: occ
 real(8), target, intent(in) :: U
-real(8), target, intent(in) :: J
+real(8), target, intent(in) :: J(3)
 real(8), target, intent(in) :: alpha
 real(8), target, intent(in) :: beta
 real(8), target, intent(in) :: J0

diff --git a/src/api/sirius_api.cpp b/src/api/sirius_api.cpp
@@ -2205,7 +2205,7 @@ sirius_add_atom_type_radial_function(void* const* handler__, char const* atom_ty
       doc: Hubbard U parameter.
     J:
       type: double
-      attr: in, required
+      attr: in, required, dimension(3)
       doc: Exchange J parameter for the full interaction treatment.
     alpha:
       type: double
@@ -2236,8 +2236,10 @@ sirius_set_atom_type_hubbard(void* const* handler__, char const* label__, int co
                 auto& type    = sim_ctx.unit_cell().atom_type(std::string(label__));
                 type.hubbard_correction(true);
                 if (type.file_name().empty()) {
-                    type.add_hubbard_orbital(*n__, *l__, *occ__, *U__, J__[1], J__, *alpha__, *beta__, *J0__,
-                                             std::vector<double>(), true);
+                    std::array<double, 3> hubbard_coeff({J__[0], J__[1], J__[2]});
+
+                    type.add_hubbard_orbital(*n__, *l__, *occ__, *U__, hubbard_coeff[0], hubbard_coeff, *alpha__,
+                                             *beta__, *J0__, std::vector<double>(), true);
                 } else {
                     // we use a an external file containing the potential
                     // information which means that we do not have all information

diff --git a/src/api/sirius_c_headers.h b/src/api/sirius_c_headers.h
@@ -1123,7 +1123,7 @@ sirius_add_atom_type_radial_function(void* const* handler__, char const* atom_ty
       doc: Hubbard U parameter.
     J:
       type: double
-      attr: in, required
+      attr: in, required, dimension(3)
       doc: Exchange J parameter for the full interaction treatment.
     alpha:
       type: double
@@ -3262,8 +3262,8 @@ sirius_load_state(void** gs_handler__, const char* file_name__, int* error_code_
       doc: Error code.
 */
 void
-sirius_access_density_matrix(void** gs_handler__, const char* access_type__, int const* ia__, double complex* dm__, int const* ld__,
-                          int* error_code__);
+sirius_access_density_matrix(void** gs_handler__, const char* access_type__, int const* ia__,
+                             double complex* dm__, int const* ld__, int* error_code__);
 
 /*
 sirius_access_local_occupation_matrix:
@@ -3307,8 +3307,9 @@ sirius_access_density_matrix(void** gs_handler__, const char* access_type__, int
       doc: Error code.
 */
 void
-sirius_access_local_occupation_matrix(void** handler__, const char* access_type__, int const* ia__, int const* n__, int const* l__, int const* spin__,
-                                   double complex* occ_mtrx__, int const* ld__, int* error_code__);
+sirius_access_local_occupation_matrix(void** handler__, const char* access_type__, int const* ia__, int const* n__,
+                                      int const* l__, int const* spin__, double complex* occ_mtrx__,
+                                      int const* ld__, int* error_code__);
 
 /*
 sirius_access_nonlocal_occupation_matrix:
@@ -3360,9 +3361,10 @@ sirius_access_local_occupation_matrix(void** handler__, const char* access_type_
       doc: Error code.
 */
 void
-sirius_access_nonlocal_occupation_matrix(void** handler__, const char* access_type__, int const* atom_pair__, int const* n__, int const* l__,
-                                      int const* spin__, int const* T__, double complex* occ_mtrx__,
-                                      int const* ld1__, int const* ld2__, int* error_code__);
+sirius_access_nonlocal_occupation_matrix(void** handler__, const char* access_type__, int const* atom_pair__,
+                                         int const* n__, int const* l__, int const* spin__, int const* T__,
+                                         double complex* occ_mtrx__, int const* ld1__, int const* ld2__,
+                                         int* error_code__);
 
 /*
 sirius_get_major_version:

diff --git a/src/context/config.hpp b/src/context/config.hpp
@@ -265,7 +265,7 @@ class config_t
         }
         /// Maximum orbital quantum number for which spherical coverage need to be generated.
         /**
-            This option can be used to increase shpherical coverage in muffin-tins. Impacts generation of XC potential.
+            This option can be used to increase spherical coverage in muffin-tins. Impacts generation of XC potential.
         */
         inline auto sht_lmax() const
         {
@@ -1882,9 +1882,13 @@ class config_t
             {
                 return dict_.at("J").get<double>();
             }
-            auto BE2() const
+            auto B() const
             {
-                return dict_.at("BE2").get<double>();
+                return dict_.at("B").get<double>();
+            }
+            auto E2() const
+            {
+                return dict_.at("E2").get<double>();
             }
             auto E3() const
             {

diff --git a/src/context/input_schema.json b/src/context/input_schema.json
@@ -1018,7 +1018,11 @@
                                 "type": "number",
                                 "default": 0
                             },
-                            "BE2": {
+                            "B": {
+                                "type": "number",
+                                "default": 0
+                            },
+                            "E2": {
                                 "type": "number",
                                 "default": 0
                             },

diff --git a/src/hamiltonian/initialize_subspace.hpp b/src/hamiltonian/initialize_subspace.hpp
@@ -30,19 +30,19 @@ initialize_subspace(Hamiltonian_k<T> const& Hk__, K_point<T>& kp__, int num_ao__
 
     auto& ctx = Hk__.H0().ctx();
 
-    if (ctx.cfg().control().verification() >= 2) {
-        auto eval = diag_S_davidson<T, F>(Hk__, kp__);
-        if (eval[0] <= 0) {
-            std::stringstream s;
-            s << "S-operator matrix is not positive definite\n"
-              << "  lowest eigen-value: " << eval[0];
-            RTE_WARNING(s);
-        } else {
-            std::stringstream s;
-            s << "S-matrix is OK! Minimum eigen-value : " << eval[0];
-            RTE_OUT(ctx.out(1)) << s.str() << std::endl;
-        }
-    }
+    //if (ctx.cfg().control().verification() >= 2) {
+    //    auto eval = diag_S_davidson<T, F>(Hk__, kp__);
+    //    if (eval[0] <= 0) {
+    //        std::stringstream s;
+    //        s << "S-operator matrix is not positive definite\n"
+    //          << "  lowest eigen-value: " << eval[0];
+    //        RTE_WARNING(s);
+    //    } else {
+    //        std::stringstream s;
+    //        s << "S-matrix is OK! Minimum eigen-value : " << eval[0];
+    //        RTE_OUT(ctx.out(1)) << s.str() << std::endl;
+    //    }
+    //}
 
     auto pcs = env::print_checksum();
 

diff --git a/src/hamiltonian/non_local_operator.cpp b/src/hamiltonian/non_local_operator.cpp
@@ -370,22 +370,28 @@ Q_operator<T>::initialize()
 }
 
 template <typename T>
-U_operator<T>::U_operator(Simulation_context const& ctx__, Hubbard_matrix const& um1__, std::array<double, 3> vk__)
+U_operator<T>::U_operator(Simulation_context const& ctx__, Hubbard_matrix const& um1__, r3::vector<double> vk__)
     : ctx_(ctx__)
     , vk_(vk__)
 {
     if (!ctx_.hubbard_correction()) {
         return;
     }
-    /* a pair of "total number, offests" for the Hubbard orbitals idexing */
+    /* a pair of "total number, offsets" for the Hubbard orbitals idexing */
     auto r                 = ctx_.unit_cell().num_hubbard_wf();
     this->nhwf_            = r.first;
     this->offset_          = um1__.offset();
     this->atomic_orbitals_ = um1__.atomic_orbitals();
     for (int j = 0; j < ctx_.num_mag_dims() + 1; j++) {
-        um_[j] = la::dmatrix<std::complex<T>>(r.first, r.first);
+        um_[j] = mdarray<std::complex<T>, 2>({r.first, r.first});
         um_[j].zero();
     }
+    if (env::print_checksum()) {
+        auto cs1 = um1__.local_checksum();
+        auto cs2 = um1__.nonlocal_checksum();
+        print_checksum("U_matrix_local", cs1, RTE_OUT(ctx_.out()));
+        print_checksum("U_matrix_nonlocal", cs2, RTE_OUT(ctx_.out()));
+    }
 
     /* copy local blocks */
     for (int at_lvl = 0; at_lvl < static_cast<int>(um1__.atomic_orbitals().size()); at_lvl++) {
@@ -427,12 +433,12 @@ U_operator<T>::U_operator(Simulation_context const& ctx__, Hubbard_matrix const&
         }
     }
     for (int is = 0; is < ctx_.num_spins(); is++) {
-        auto diff = check_hermitian(um_[is], r.first);
-        if (diff > 1e-10) {
-            RTE_THROW("um is not Hermitian");
-        }
+        //auto diff = check_hermitian(um_[is], r.first);
+        //if (diff > 1e-10) {
+        //    RTE_THROW("um is not Hermitian");
+        //}
         if (env::print_checksum()) {
-            print_checksum("um" + std::to_string(is), um_[is].checksum(r.first, r.first), RTE_OUT(ctx_.out()));
+            print_checksum("U_matrix_k_" + std::to_string(is), um_[is].checksum(), RTE_OUT(ctx_.out()));
         }
         if (ctx_.processing_unit() == device_t::GPU) {
             um_[is].allocate(get_memory_pool(memory_t::device)).copy_to(memory_t::device);

diff --git a/src/hamiltonian/non_local_operator.hpp b/src/hamiltonian/non_local_operator.hpp
@@ -58,14 +58,16 @@ class U_operator
 {
   private:
     Simulation_context const& ctx_;
-    std::array<la::dmatrix<std::complex<T>>, 4> um_;
+    /// Potential correction matrix for a given k-point.
+    /** Constructed from local and non-local parts of Hubbard U matrix. */
+    std::array<mdarray<std::complex<T>, 2>, 4> um_;
     std::vector<int> offset_;
     std::vector<std::pair<int, int>> atomic_orbitals_;
     int nhwf_;
     r3::vector<double> vk_;
 
   public:
-    U_operator(Simulation_context const& ctx__, Hubbard_matrix const& um1__, std::array<double, 3> vk__);
+    U_operator(Simulation_context const& ctx__, Hubbard_matrix const& um1__, r3::vector<double> vk__);
     ~U_operator() = default;
 
     inline auto
@@ -98,9 +100,9 @@ class U_operator
     }
 
     std::complex<T> const*
-    at(memory_t mem__, const int idx1, const int idx2, const int idx3) const
+    at(memory_t mem__, const int m1, const int m2, const int j) const
     {
-        return um_[idx3].at(mem__, idx1, idx2);
+        return um_[j].at(mem__, m1, m2);
     }
 
     matrix<std::complex<T>> const&

diff --git a/src/hubbard/hubbard_matrix.hpp b/src/hubbard/hubbard_matrix.hpp
@@ -23,13 +23,12 @@ class Hubbard_matrix
 {
   protected:
     Simulation_context const& ctx_;
-    /// Local part of Hubbard matrix
     int num_steps_{0};
     double constraint_error_{1.0};
-    /// table indicating if we should apply constraints on the hubbard occupation
+    /// Table indicating if we should apply constraints on the hubbard occupation
     /// to given atomic orbital group
     std::vector<bool> apply_constraints_;
-    /// occupancy matrix for each atomic level (n,l)
+    /// Local part of Hubbard matrix
     std::vector<mdarray<std::complex<double>, 3>> local_;
     /// Non-local part of Hubbard matrix.
     std::vector<mdarray<std::complex<double>, 3>> nonlocal_;
@@ -269,6 +268,26 @@ class Hubbard_matrix
 
         return -1;
     }
+
+    inline auto
+    local_checksum() const
+    {
+        std::complex<double> sum(0, 0);
+        for (auto& e : local_) {
+            sum += e.checksum();
+        }
+        return sum;
+    }
+
+    inline auto
+    nonlocal_checksum() const
+    {
+        std::complex<double> sum(0, 0);
+        for (auto& e : nonlocal_) {
+            sum += e.checksum();
+        }
+        return sum;
+    }
 };
 
 inline void

diff --git a/src/hubbard/hubbard_potential_energy.cpp b/src/hubbard/hubbard_potential_energy.cpp
@@ -57,8 +57,9 @@ generate_potential_collinear_nonlocal(Simulation_context const& ctx__, const int
 }
 
 static void
-generate_potential_collinear_local(Simulation_context const& ctx__, Atom_type const& atom_type__, const int idx_hub_wf,
-                                   mdarray<std::complex<double>, 3> const& om__, mdarray<std::complex<double>, 3>& um__)
+generate_potential_collinear_local(Simulation_context const& ctx__, Atom_type const& atom_type__,
+                                   const int idx_hub_wf__, mdarray<std::complex<double>, 3> const& om__,
+                                   mdarray<std::complex<double>, 3>& um__)
 {
     /* quick exit */
     if (!atom_type__.hubbard_correction()) {
@@ -68,7 +69,9 @@ generate_potential_collinear_local(Simulation_context const& ctx__, Atom_type co
     um__.zero();
 
     /* single orbital implementation */
-    auto& hub_wf = atom_type__.lo_descriptor_hub(idx_hub_wf);
+    auto& hub_wf = atom_type__.lo_descriptor_hub(idx_hub_wf__);
+
+    hub_wf.print_info(RTE_OUT(ctx__.out()));
 
     if (!hub_wf.use_for_calculation()) {
         return;
@@ -593,6 +596,12 @@ generate_potential(Hubbard_matrix const& om__, Hubbard_matrix& um__)
             ::sirius::generate_potential_collinear_nonlocal(ctx, i, om__.nonlocal(i), um__.nonlocal(i));
         }
     }
+    if (env::print_checksum()) {
+        print_checksum("om_local", om__.local_checksum(), RTE_OUT(ctx.out()));
+        print_checksum("om_nonlocal", om__.nonlocal_checksum(), RTE_OUT(ctx.out()));
+        print_checksum("um_local", um__.local_checksum(), RTE_OUT(ctx.out()));
+        print_checksum("um_nonlocal", um__.nonlocal_checksum(), RTE_OUT(ctx.out()));
+    }
 }
 
 double