bug in disjoin and adjoin when spin is not paired

gitpeterwind · Aug 24, 2023 · dd656a5 · dd656a5
1 parent 95aad36
commit dd656a5
Show file tree

Hide file tree

Showing 6 changed files with 36 additions and 17 deletions.
diff --git a/src/initial_guess/core.cpp b/src/initial_guess/core.cpp
@@ -129,8 +129,8 @@ bool initial_guess::core::setup(OrbitalVector &Phi, double prec, const Nuclei &n
     initial_guess::core::rotate_orbitals(Phi, prec, U, Psi);
     initial_guess::core::rotate_orbitals(Phi_a, prec, U, Psi);
     initial_guess::core::rotate_orbitals(Phi_b, prec, U, Psi);
-    for (auto &phi_a : Phi_a) Phi.push_back(phi_a);
-    for (auto &phi_b : Phi_b) Phi.push_back(phi_b);
+    Phi = orbital::adjoin(Phi_a, Phi_b);
+
     V.clear();
     p.clear();
 

diff --git a/src/initial_guess/cube.cpp b/src/initial_guess/cube.cpp
@@ -81,8 +81,7 @@ bool initial_guess::cube::setup(OrbitalVector &Phi, double prec, const std::stri
     success &= initial_guess::cube::project_mo(Phi_b, prec, file_b);
 
     // Collect orbitals into one vector
-    for (auto &phi_a : Phi_a) Phi.push_back(phi_a);
-    for (auto &phi_b : Phi_b) Phi.push_back(phi_b);
+    Phi = orbital::adjoin(Phi_a, Phi_b);
 
     return success;
 }

diff --git a/src/initial_guess/gto.cpp b/src/initial_guess/gto.cpp
@@ -24,9 +24,9 @@
  */
 
 #include <MRCPP/Gaussians>
+#include <MRCPP/Parallel>
 #include <MRCPP/Printer>
 #include <MRCPP/Timer>
-#include <MRCPP/Parallel>
 
 #include "gto.h"
 
@@ -94,8 +94,7 @@ bool initial_guess::gto::setup(OrbitalVector &Phi, double prec, double screen, c
     initial_guess::gto::project_mo(Phi_b, prec, bas_file, mob_file, screen);
 
     // Collect orbitals into one vector
-    for (auto &phi_a : Phi_a) Phi.push_back(phi_a);
-    for (auto &phi_b : Phi_b) Phi.push_back(phi_b);
+    Phi = orbital::adjoin(Phi_a, Phi_b);
 
     return true;
 }

diff --git a/src/initial_guess/mw.cpp b/src/initial_guess/mw.cpp
@@ -74,8 +74,7 @@ bool initial_guess::mw::setup(OrbitalVector &Phi, double prec, const std::string
     success &= initial_guess::mw::project_mo(Phi_b, prec, file_b);
 
     // Collect orbitals into one vector
-    for (auto &phi_a : Phi_a) Phi.push_back(phi_a);
-    for (auto &phi_b : Phi_b) Phi.push_back(phi_b);
+    Phi = orbital::adjoin(Phi_a, Phi_b);
 
     return success;
 }

diff --git a/src/initial_guess/sad.cpp b/src/initial_guess/sad.cpp
@@ -76,7 +76,6 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c
     print_utils::text(0, "Zeta quality", std::to_string(zeta));
     mrcpp::print::separator(0, '~', 2);
 
-
     // Make Fock operator contributions
     Timer t_tot, t_lap;
     auto P_p = std::make_shared<mrcpp::PoissonOperator>(*MRA, prec);
@@ -132,8 +131,7 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c
     initial_guess::core::rotate_orbitals(Phi, prec, U, Psi);
     initial_guess::core::rotate_orbitals(Phi_a, prec, U, Psi);
     initial_guess::core::rotate_orbitals(Phi_b, prec, U, Psi);
-    for (auto &phi_a : Phi_a) Phi.push_back(phi_a);
-    for (auto &phi_b : Phi_b) Phi.push_back(phi_b);
+    Phi = orbital::adjoin(Phi_a, Phi_b);
     p.clear();
     V.clear();
 
@@ -156,7 +154,6 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c
     print_utils::text(0, "AO basis    ", "3-21G");
     mrcpp::print::separator(0, '~', 2);
 
-
     // Make Fock operator contributions
     Timer t_tot, t_lap;
     auto P_p = std::make_shared<mrcpp::PoissonOperator>(*MRA, prec);
@@ -210,8 +207,7 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c
     initial_guess::core::rotate_orbitals(Phi, prec, U, Psi);
     initial_guess::core::rotate_orbitals(Phi_a, prec, U, Psi);
     initial_guess::core::rotate_orbitals(Phi_b, prec, U, Psi);
-    for (auto &phi_a : Phi_a) Phi.push_back(phi_a);
-    for (auto &phi_b : Phi_b) Phi.push_back(phi_b);
+    Phi = orbital::adjoin(Phi_a, Phi_b);
     p.clear();
     V.clear();
 

diff --git a/src/qmfunctions/orbital_utils.cpp b/src/qmfunctions/orbital_utils.cpp
@@ -331,8 +331,24 @@ OrbitalVector orbital::param_copy(const OrbitalVector &Phi) {
  */
 OrbitalVector orbital::adjoin(OrbitalVector &Phi_a, OrbitalVector &Phi_b) {
     OrbitalVector out;
-    for (auto &phi : Phi_a) out.push_back(phi);
-    for (auto &phi : Phi_b) out.push_back(phi);
+    for (auto &phi : Phi_a) {
+        if (phi.getRank() != out.size()) {
+            // need to send orbital from owner to new owner
+            if (mrcpp::mpi::my_orb(phi)) { mrcpp::mpi::send_function(phi, out.size() % mrcpp::mpi::wrk_size, phi.getRank(), mrcpp::mpi::comm_wrk); }
+            if (mrcpp::mpi::my_orb(out.size())) { mrcpp::mpi::recv_function(phi, phi.getRank() % mrcpp::mpi::wrk_size, phi.getRank(), mrcpp::mpi::comm_wrk); }
+        }
+        phi.setRank(out.size());
+        out.push_back(phi);
+    }
+    for (auto &phi : Phi_b) {
+        if (phi.getRank() != out.size()) {
+            // need to send orbital from owner to new owner
+            if (mrcpp::mpi::my_orb(phi)) { mrcpp::mpi::send_function(phi, out.size() % mrcpp::mpi::wrk_size, phi.getRank(), mrcpp::mpi::comm_wrk); }
+            if (mrcpp::mpi::my_orb(out.size())) { mrcpp::mpi::recv_function(phi, phi.getRank() % mrcpp::mpi::wrk_size, phi.getRank(), mrcpp::mpi::comm_wrk); }
+        }
+        phi.setRank(out.size());
+        out.push_back(phi);
+    }
     Phi_a.clear();
     Phi_b.clear();
     return out;
@@ -350,9 +366,19 @@ OrbitalVector orbital::disjoin(OrbitalVector &Phi, int spin) {
     OrbitalVector tmp;
     for (auto &i : Phi) {
         if (i.spin() == spin) {
+            if (i.getRank() != out.size()) {
+                // need to send orbital from owner to new owner
+                if (mrcpp::mpi::my_orb(i)) { mrcpp::mpi::send_function(i, out.size() % mrcpp::mpi::wrk_size, i.getRank(), mrcpp::mpi::comm_wrk); }
+                if (mrcpp::mpi::my_orb(out.size())) { mrcpp::mpi::recv_function(i, i.getRank() % mrcpp::mpi::wrk_size, i.getRank(), mrcpp::mpi::comm_wrk); }
+            }
             i.setRank(out.size());
             out.push_back(i);
         } else {
+            if (i.getRank() != tmp.size()) {
+                // need to send orbital from owner to new owner
+                if (mrcpp::mpi::my_orb(i)) { mrcpp::mpi::send_function(i, tmp.size() % mrcpp::mpi::wrk_size, i.getRank(), mrcpp::mpi::comm_wrk); }
+                if (mrcpp::mpi::my_orb(tmp.size())) { mrcpp::mpi::recv_function(i, i.getRank() % mrcpp::mpi::wrk_size, i.getRank(), mrcpp::mpi::comm_wrk); }
+            }
             i.setRank(tmp.size());
             tmp.push_back(i);
         }