fix bug in reallocating HContainter when not all the processors have …

…elements of (0,0) atom pair
maki49 · Aug 30, 2024 · 69e3487 · 69e3487
1 parent a582874
commit 69e3487
Showing 1 changed file with 31 additions and 32 deletions.
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp
@@ -24,18 +24,20 @@ namespace hamilt
             for (auto& Htmp2 : Htmp1.second)
             {
                 const int& iat1 = Htmp2.first.first;
-                const Abfs::Vector3_Order<int>& R = RI_Util::array3_to_Vector3(Htmp2.first.second);
-                BaseMatrix<TR>* HlocR = hR->find_matrix(iat0, iat1, R.x, R.y, R.z);
-                if (HlocR == nullptr)
-                { // add R to HContainer
-                    need_allocate = true;
-                    AtomPair<TR> tmp(iat0, iat1, R.x, R.y, R.z, hR->find_pair(0, 0)->get_paraV());
-                    hR->insert_pair(tmp);
+                if (hR->find_pair(iat0, iat1))
+                {
+                    const Abfs::Vector3_Order<int>& R = RI_Util::array3_to_Vector3(Htmp2.first.second);
+                    BaseMatrix<TR>* HlocR = hR->find_matrix(iat0, iat1, R.x, R.y, R.z);
+                    if (HlocR == nullptr)
+                    { // add R to HContainer
+                        need_allocate = true;
+                        AtomPair<TR> tmp(iat0, iat1, R.x, R.y, R.z, hR->get_paraV());
+                        hR->insert_pair(tmp);
+                    }
                 }
             }
         }
-        if (need_allocate) { hR->allocate(nullptr, true);
-}
+        if (need_allocate) { hR->allocate(nullptr, true); }
     }
     /// allocate according to BvK cells, used in scf
     template <typename TR>
@@ -45,28 +47,27 @@ namespace hamilt
     {
         auto Rs = RI_Util::get_Born_von_Karmen_cells(Rs_period);
         bool need_allocate = false;
-        for (int iat0 = 0;iat0 < GlobalC::ucell.nat;++iat0)
+        for (int iap = 0;iap < hR->size_atom_pairs();++iap)
         {
-            for (int iat1 = 0;iat1 < GlobalC::ucell.nat;++iat1)
+            hamilt::AtomPair<TR>& ap = hR->get_atom_pair(iap);
+            const int iat0 = ap.get_atom_i();
+            const int iat1 = ap.get_atom_j();
+            for (auto& cell : Rs)
             {
-                for (auto& cell : Rs)
-                {
-                    const Abfs::Vector3_Order<int>& R = RI_Util::array3_to_Vector3(
-                        (cell_nearest ?
-                            cell_nearest->get_cell_nearest_discrete(iat0, iat1, cell)
-                            : cell));
-                    BaseMatrix<TR>* HlocR = hR->find_matrix(iat0, iat1, R.x, R.y, R.z);
-                    if (HlocR == nullptr)
-                    { // add R to HContainer
-                        need_allocate = true;
-                        AtomPair<TR> tmp(iat0, iat1, R.x, R.y, R.z, hR->find_pair(0, 0)->get_paraV());
-                        hR->insert_pair(tmp);
-                    }
+                const Abfs::Vector3_Order<int>& R = RI_Util::array3_to_Vector3(
+                    (cell_nearest ?
+                        cell_nearest->get_cell_nearest_discrete(iat0, iat1, cell)
+                        : cell));
+                BaseMatrix<TR>* HlocR = hR->find_matrix(iat0, iat1, R.x, R.y, R.z);
+                if (HlocR == nullptr)
+                { // add R to HContainer
+                    need_allocate = true;
+                    AtomPair<TR> tmp(iat0, iat1, R.x, R.y, R.z, hR->get_paraV());
+                    hR->insert_pair(tmp);
                 }
             }
         }
-        if (need_allocate) { hR->allocate(nullptr, true);
-}
+        if (need_allocate) { hR->allocate(nullptr, true); }
     }
 
 template <typename TK, typename TR>
@@ -150,8 +151,7 @@ OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
                         this->restart = GlobalC::restart.load_disk(
                             "Hexx", ik,
                             pv->get_local_size(), this->Hexxd_k_load[ik].data(), false);
-                        if (!this->restart) { break;
-}
+                        if (!this->restart) { break; }
                     }
                 }
                 else
@@ -163,8 +163,7 @@ OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
                         this->restart = GlobalC::restart.load_disk(
                             "Hexx", ik,
                             pv->get_local_size(), this->Hexxc_k_load[ik].data(), false);
-                        if (!this->restart) { break;
-}
+                        if (!this->restart) { break; }
                     }
                 }
             }
@@ -182,8 +181,8 @@ OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
 
             if (!this->restart) {
                 std::cout << "WARNING: Hexx not found, restart from the non-exx loop." << std::endl
-                << "If the loaded charge density is EXX-solved, this may lead to poor convergence." << std::endl;
-}
+                    << "If the loaded charge density is EXX-solved, this may lead to poor convergence." << std::endl;
+            }
             GlobalC::restart.info_load.load_H_finish = this->restart;
         }
     }