diff --git a/src/green/impurity/inchworm_inpurity_solver.h b/src/green/impurity/inchworm_inpurity_solver.h index fec19e2..78bc395 100644 --- a/src/green/impurity/inchworm_inpurity_solver.h +++ b/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> delta_t_m(delta_t.data(), _uxl.shape()[0], ns * nio * nio); CMMatrixX> 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; } } diff --git a/test/data/transform.h5 b/test/data/transform.h5 index d380081..8f91e93 100644 Binary files a/test/data/transform.h5 and b/test/data/transform.h5 differ diff --git a/test/impurity_solver_test.cpp b/test/impurity_solver_test.cpp index 922bed1..2a500d2 100644 --- a/test/impurity_solver_test.cpp +++ b/test/impurity_solver_test.cpp @@ -26,6 +26,7 @@ #include #include +#include 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("spin_symm", "", false); p.define("bath_file", "", bath_file); - p.define("impurity_solver", "", "ED"); + p.define("impurity_solver", "", impurity_solver_type); p.define("impurity_solver_exec", "", "/bin/true"); p.define("impurity_solver_params", "", ""); p.define("dc_data_prefix", "", ""); p.define("seet_root_dir", "", TEST_PATH + ""s); p.define("seet_input", "", input_file); + if (impurity_solver_type == "INCHWORM") { + p.define("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);