add impurity prefix to output files; fix spin-symmetry in Uijkl

src/green/impurity/inchworm_inpurity_solver.h

@@ -47,10 +47,12 @@ namespace green::impurity {
       ztensor<4> sigma_new(delta_w.shape());
       size_t     nio = h_core.shape()[1];
       size_t     ns  = h_core.shape()[0];
+      bool       rhf = (ns == 1);
+      if (rhf) std::cout << "Detected ns = 1; Assuming this is RHF. Support for relativistic cases not implemented yet." << std::endl;
       // One-body term
       // Static:
       {
-        std::ofstream hooping_file("hopping.txt");
+        std::ofstream hooping_file("imp_" + std::to_string(imp_n) + "_hopping.txt");
         for (size_t i = 0; i < nio; ++i) {
           for (size_t s1 = 0; s1 < ns; ++s1) {
             for (size_t j = 0; j < nio; ++j) {
@@ -75,7 +77,7 @@ namespace green::impurity {
         MMatrixX<std::complex<double>>  delta_t_m(delta_t.data(), _uxl.shape()[0], ns * nio * nio);
         CMMatrixX<std::complex<double>> delta_w_m(delta_w.data(), delta_w.shape()[0], ns * nio * nio);
         delta_t_m = Ttc_even * Tcn * delta_w_m * std::sqrt(2.0 / ft.sd().beta());
-        std::ofstream delta_file("delta.txt");
+        std::ofstream delta_file("imp_" + std::to_string(imp_n) + "_delta.txt");
         for (size_t t = 0; t < delta_t.shape()[0]; ++t) {
           for (size_t i = 0; i < nio; ++i) {
             for (size_t s1 = 0; s1 < ns; ++s1) {
@@ -92,29 +94,37 @@ namespace green::impurity {
           }
         }
       }
-      // Two-body term
+      // Two-body term with spin symmetry
+      // NOTE: Whether or not we break spin-symmetry in G, or even in x2c relativistic cases, U will still have spin-symmetry.
       {
-        // transform interaction into physics convention
-        auto   interaction_phys = ndarray::transpose(interaction, "ijkl->ikjl");
+        size_t ndim_increment   = (rhf) ? 1 : 4; // restricted vs unrestricted
+        size_t spin_decrement   = (rhf) ? 0 : 2; // avoid double creation/annihilation in same spin-orbitals
         size_t non_zero         = 0;
-        for (size_t i = 0; i < nio * ns; ++i) {
-          for (size_t j = 0; j < nio * ns; ++j) {
-            for (size_t k = 0; k < nio * ns; ++k) {
-              for (size_t l = 0; l < nio * ns; ++l) {
-                size_t I = i / ns;
-                size_t J = j / ns;
-                size_t K = k / ns;
-                size_t L = l / ns;
-                if (std::abs(interaction_phys(I, J, K, L)) > 1e-10) {
-                  ++non_zero;
+        for (size_t I = 0; I < nio; ++I) {
+          for (size_t J = 0; J < nio; ++J) {
+            for (size_t K = 0; K < nio; ++K) {
+              for (size_t L = 0; L < nio; ++L) {
+                if (std::abs(interaction(I, J, K, L)) < 1e-10) continue; // forget all logic and skip if interaction is zero
+                if (rhf) {
+                  non_zero += ndim_increment;
+                  continue;
                 }
+                // uhf case
+                // if I == K or J == L, it warrants that the spin label of I and K will be opposite,
+                // and by construction spin of I = spin of J and spin of K = spin of L,
+                // so when (I==K) or (J==L) or both, we need to remove two combinations:
+                // (up up up up) and (dn dn dn dn)
+                non_zero += ndim_increment;
+                if (I == K || J == L) non_zero -= spin_decrement;
               }
             }
           }
         }
-        std::ofstream U_file("Uijkl.txt");
+        std::cout << "NOTE: Will write interaction tensor in the Chemist notation!" << std::endl;
+        std::cout << "i.,e., U(ijkl) cdag_i cdag_k c_l c_j" << std::endl;
+        std::ofstream U_file("imp_" + std::to_string(imp_n) + "_Uijkl.txt");
         U_file << non_zero << "\n";
-        int idx = 0;
+        size_t idx = 0;
         for (size_t i = 0; i < nio * ns; ++i) {
           for (size_t j = 0; j < nio * ns; ++j) {
             for (size_t k = 0; k < nio * ns; ++k) {
@@ -123,8 +133,20 @@ namespace green::impurity {
                 size_t J = j / ns;
                 size_t K = k / ns;
                 size_t L = l / ns;
-                if (std::abs(interaction_phys(I, J, K, L)) > 1e-10) {
-                  U_file << idx << "\t" << i << " " << j << " " << k << " " << l << " " << interaction_phys(I, J, K, L) << " " << 0.0 << "\n";
+                if (std::abs(interaction(I, J, K, L)) < 1e-10) continue;
+                if (rhf) {
+                  U_file << idx << "\t" << i << " " << j << " " << k << " " << l << " " << interaction(I, J, K, L) << " " << 0.0 << "\n";
+                  ++idx;
+                } else {
+                  // Deal with spin only for UHF
+                  size_t s1 = i % ns;
+                  size_t s2 = j % ns;
+                  size_t s3 = k % ns;
+                  size_t s4 = l % ns;
+                  if (i == k || j == l) continue;       // ignore double creation/annihilation in same spin-orbitals (different from I,J,K,L which are atomic orbitals)
+                  if (s1 != s2 || s3 != s4) continue;   // ignore spin-flip terms
+                  // what remains is a valid term in Uijkl
+                  U_file << idx << "\t" << i << " " << j << " " << k << " " << l << " " << interaction(I, J, K, L) << " " << 0.0 << "\n";
                   ++idx;
                 }
               }

test/impurity_solver_test.cpp

@@ -26,6 +26,7 @@
 #include <catch2/catch_session.hpp>
 
 #include <mpi.h>
+#include <filesystem>
 
 namespace green::impurity {
   ztensor<4> compute_local_obj(const ztensor<5>& obj, const ztensor<3>& x_k, const grids::transformer_t& ft, const bz_utils_t& bz_utils, size_t ns, size_t nso) {
@@ -57,7 +58,7 @@ namespace green::impurity {
   }
 }
 
-TEST_CASE("Impurity Solver") {
+void impurity_solver_test(std::string impurity_solver_type) {
   std::string test_file   = TEST_PATH + "/data.h5"s;
   std::string bath_file   = TEST_PATH + "/bath.txt"s;
   std::string input_file   = TEST_PATH + "/transform.h5"s;
@@ -69,12 +70,15 @@ TEST_CASE("Impurity Solver") {
   green::grids::define_parameters(p);
   p.define<bool>("spin_symm", "", false);
   p.define<std::string>("bath_file", "", bath_file);
-  p.define<std::string>("impurity_solver", "", "ED");
+  p.define<std::string>("impurity_solver", "", impurity_solver_type);
   p.define<std::string>("impurity_solver_exec", "", "/bin/true");
   p.define<std::string>("impurity_solver_params", "", "");
   p.define<std::string>("dc_data_prefix", "", "");
   p.define<std::string>("seet_root_dir", "", TEST_PATH + ""s);
   p.define<std::string>("seet_input", "", input_file);
+  if (impurity_solver_type == "INCHWORM") {
+    p.define<std::string>("itermax", "", "1");
+  }
 
   p.parse("test --BETA 100 --grid_file " + grid_file + " --input_file " + weak_input_file );
 
@@ -123,6 +127,93 @@ TEST_CASE("Impurity Solver") {
   solver.solve(mu, ovlp, h_core, sigma1, sigma, g);
 }
 
+TEST_CASE("Impurity Solver") {
+  SECTION("ED") { impurity_solver_test("ED"); }
+  SECTION("INCHWORM") {
+    impurity_solver_test("INCHWORM");
+    // Check if Hamiltonian data files were created successfully for all impurities
+    REQUIRE(std::filesystem::exists("imp_0_hopping.txt"));
+    REQUIRE(std::filesystem::exists("imp_0_delta.txt"));
+    REQUIRE(std::filesystem::exists("imp_0_Uijkl.txt"));
+    REQUIRE(std::filesystem::exists("imp_1_hopping.txt"));
+    REQUIRE(std::filesystem::exists("imp_1_delta.txt"));
+    REQUIRE(std::filesystem::exists("imp_1_Uijkl.txt"));
+    // Check hopping file data
+    {
+      std::ifstream hop_file("imp_0_hopping.txt");
+      std::string line;
+      size_t n_terms = 16;
+      size_t count = 0;
+      while (std::getline(hop_file, line)) {
+        ++count;
+      }
+      REQUIRE(n_terms == count);
+    }
+    {
+      std::ifstream hop_file("imp_1_hopping.txt");
+      std::string line;
+      size_t n_terms = 16;
+      size_t count = 0;
+      while (std::getline(hop_file, line)) {
+        ++count;
+      }
+      REQUIRE(n_terms == count);
+    }
+    // Check Delta file data
+    {
+      std::ifstream delta_file("imp_0_delta.txt");
+      std::string line;
+      size_t n_terms = 16 * 101;
+      size_t count = 0;
+      while (std::getline(delta_file, line)) {
+        ++count;
+      }
+      REQUIRE(n_terms == count);
+    }
+    {
+      std::ifstream delta_file("imp_1_delta.txt");
+      std::string line;
+      size_t n_terms = 16 * 101;
+      size_t count = 0;
+      while (std::getline(delta_file, line)) {
+        ++count;
+      }
+      REQUIRE(n_terms == count);
+    }
+    // Check Uijkl file data
+    {
+      std::ifstream U_file("imp_0_Uijkl.txt");
+      std::string line;
+      // read number of non-zero terms from the heaer
+      std::getline(U_file, line);
+      size_t n_terms = std::stol(line);
+      size_t count = 0;
+      while (std::getline(U_file, line)) {
+        ++count;
+      }
+      REQUIRE(n_terms == count);
+    }
+    {
+      std::ifstream U_file("imp_1_Uijkl.txt");
+      std::string line;
+      std::getline(U_file, line);
+      size_t n_terms = std::stol(line);
+      size_t count = 0;
+      while (std::getline(U_file, line)) {
+        ++count;
+      }
+      REQUIRE(n_terms == count);
+    }
+    // Cleanup: Remove all files (best-effort; do not fail test on cleanup)
+    std::filesystem::remove("imp_0_hopping.txt");
+    std::filesystem::remove("imp_0_delta.txt");
+    std::filesystem::remove("imp_0_Uijkl.txt");
+    std::filesystem::remove("imp_1_hopping.txt");
+    std::filesystem::remove("imp_1_delta.txt");
+    std::filesystem::remove("imp_1_Uijkl.txt");
+  }
+}
+
 int main(int argc, char** argv) {
   MPI_Init(&argc, &argv);
   int result = Catch::Session().run(argc, argv);