diff --git a/cmake/sources.cmake b/cmake/sources.cmake index d3103d21..2aa1ff70 100644 --- a/cmake/sources.cmake +++ b/cmake/sources.cmake @@ -81,6 +81,7 @@ set(XDIAG_SOURCES symmetries/group_action/group_action_sublattice.cpp symmetries/group_action/sublattice_stability.cpp + operators/coupling.cpp operators/bond.cpp operators/bondlist.cpp operators/compiler.cpp diff --git a/misc/data/kitaev_gamma/lattice-files/honeycomb.8.HeisenbergKitaevGamma.fsl.lat b/misc/data/kitaev_gamma/lattice-files/honeycomb.8.HeisenbergKitaevGamma.fsl.lat index cf18fa0f..970670eb 100644 --- a/misc/data/kitaev_gamma/lattice-files/honeycomb.8.HeisenbergKitaevGamma.fsl.lat +++ b/misc/data/kitaev_gamma/lattice-files/honeycomb.8.HeisenbergKitaevGamma.fsl.lat @@ -37,30 +37,30 @@ HB J 2 5 HB J 6 1 HB J 3 4 HB J 7 0 -GENERICKITAEVX K 2 0 -GENERICKITAEVX K 6 4 -GENERICKITAEVX K 3 1 -GENERICKITAEVX K 7 5 -GENERICKITAEVY K 0 7 -GENERICKITAEVY K 4 3 -GENERICKITAEVY K 1 6 -GENERICKITAEVY K 5 2 -GENERICKITAEVZ K 0 3 -GENERICKITAEVZ K 4 7 -GENERICKITAEVZ K 1 2 -GENERICKITAEVZ K 5 6 -GENERICGAMMAX G 2 0 -GENERICGAMMAX G 6 4 -GENERICGAMMAX G 3 1 -GENERICGAMMAX G 7 5 -GENERICGAMMAY G 0 7 -GENERICGAMMAY G 4 3 -GENERICGAMMAY G 1 6 -GENERICGAMMAY G 5 2 -GENERICGAMMAZ G 0 3 -GENERICGAMMAZ G 4 7 -GENERICGAMMAZ G 1 2 -GENERICGAMMAZ G 5 6 +GENERICKITAEVX KX 2 0 +GENERICKITAEVX KX 6 4 +GENERICKITAEVX KX 3 1 +GENERICKITAEVX KX 7 5 +GENERICKITAEVY KY 0 7 +GENERICKITAEVY KY 4 3 +GENERICKITAEVY KY 1 6 +GENERICKITAEVY KY 5 2 +GENERICKITAEVZ KZ 0 3 +GENERICKITAEVZ KZ 4 7 +GENERICKITAEVZ KZ 1 2 +GENERICKITAEVZ KZ 5 6 +GENERICGAMMAX GX 2 0 +GENERICGAMMAX GX 6 4 +GENERICGAMMAX GX 3 1 +GENERICGAMMAX GX 7 5 +GENERICGAMMAY GY 0 7 +GENERICGAMMAY GY 4 3 +GENERICGAMMAY GY 1 6 +GENERICGAMMAY GY 5 2 +GENERICGAMMAZ GZ 0 3 +GENERICGAMMAZ GZ 4 7 +GENERICGAMMAZ GZ 1 2 +GENERICGAMMAZ GZ 5 6 [SymmetryOps]=8 [S0] 0 1 2 3 4 5 6 7 [S1] 1 0 3 2 5 4 7 6 diff --git a/tests/algorithms/arnoldi/test_arnoldi.cpp b/tests/algorithms/arnoldi/test_arnoldi.cpp index f8f1faad..69fa7655 100644 --- a/tests/algorithms/arnoldi/test_arnoldi.cpp +++ b/tests/algorithms/arnoldi/test_arnoldi.cpp @@ -55,8 +55,8 @@ TEST_CASE("ritz_vecs_arnoldi", "[arnoldi]") { BondList bonds; for (int i = 0; i < n_sites; ++i) { - bonds << Bond("HB", "J1", {i, (i + 1) % n_sites}); - bonds << Bond("HB", "J2", {i, (i + 2) % n_sites}); + bonds += Bond("HB", "J1", {i, (i + 1) % n_sites}); + bonds += Bond("HB", "J2", {i, (i + 2) % n_sites}); } bonds["J1"] = 1; bonds["J2"] = .5; diff --git a/tests/algorithms/lanczos/test_lanczos_pro.cpp b/tests/algorithms/lanczos/test_lanczos_pro.cpp index 9846c185..edb4ebcf 100644 --- a/tests/algorithms/lanczos/test_lanczos_pro.cpp +++ b/tests/algorithms/lanczos/test_lanczos_pro.cpp @@ -139,8 +139,8 @@ TEST_CASE("lanczos_pro", "[lanczos]") { std::string lfilename = XDIAG_DIRECTORY "/misc/data/shastry.16.HB.J.Jd.fsl.toml"; auto lfile = FileToml(lfilename, 'r'); - auto bonds = BondList(lfile["Interactions"]); - int n_sites = bonds.n_sites(); + auto bonds = lfile["Interactions"].as(); + int n_sites = 16; bonds["J"] = 0.63; bonds["Jd"] = 1.00; auto group = PermutationGroup(lfile["Symmetries"]); @@ -195,8 +195,8 @@ TEST_CASE("lanczos_pro", "[lanczos]") { std::string lfilename = XDIAG_DIRECTORY "/misc/data/shastry.20.HB.J.Jd.fsl.toml"; auto lfile = FileToml(lfilename, 'r'); - auto bonds = BondList(lfile["Interactions"]); - int n_sites = bonds.n_sites(); + auto bonds = lfile["Interactions"].as(); + int n_sites = 20; bonds["J"] = 0.63; bonds["Jd"] = 1.00; auto group = PermutationGroup(lfile["Symmetries"]); diff --git a/tests/algorithms/time_evolution/test_time_evolution.cpp b/tests/algorithms/time_evolution/test_time_evolution.cpp index 89eef3d7..a23e4113 100644 --- a/tests/algorithms/time_evolution/test_time_evolution.cpp +++ b/tests/algorithms/time_evolution/test_time_evolution.cpp @@ -49,7 +49,7 @@ TEST_CASE("analytic_case_free_particle_1D", "[time_evolution]") try { auto block = Electron(n_sites, nup, ndn); BondList bonds; for (int i = 0; i < n_sites; i++) { - bonds << Bond("HOP", "t", {i, (i + 1) % n_sites}); + bonds += Bond("HOP", "t", {i, (i + 1) % n_sites}); } bonds["t"] = t; bonds["U"] = 1; @@ -141,14 +141,14 @@ TEST_CASE("analytic_case_free_particle_2D", "[time_evolution]") try { // nearest neighbour hops int s_up_hop = L * ((i + 1) % L) + j; int s_dn_hop = L * i + (j + 1) % L; - bonds << Bond("HOP", "t", {s, s_up_hop}); - bonds << Bond("HOP", "t", {s, s_dn_hop}); + bonds += Bond("HOP", "t", {s, s_up_hop}); + bonds += Bond("HOP", "t", {s, s_dn_hop}); // next nearest neighbour hops int s_across_up = L * ((i + 1) % L) + (j + 1) % L; int s_across_dn = L * ((L + (i - 1) % L) % L) + (j + 1) % L; - bonds << Bond("HOP", "t1", {s, s_across_up}); - bonds << Bond("HOP", "t1", {s, s_across_dn}); + bonds += Bond("HOP", "t1", {s, s_across_up}); + bonds += Bond("HOP", "t1", {s, s_across_dn}); } } bonds["t"] = t; @@ -206,10 +206,10 @@ TEST_CASE("tj_complex_timeevo", "[time_evolution]") try { int site = y * L + x; int right = y * L + nx; int top = ny * L + x; - bonds << Bond("HOP", "T", {site, right}); - bonds << Bond("TJISING", "J", {site, right}); - bonds << Bond("HOP", "T", {site, top}); - bonds << Bond("TJISING", "J", {site, top}); + bonds += Bond("HOP", "T", {site, right}); + bonds += Bond("TJISING", "J", {site, right}); + bonds += Bond("HOP", "T", {site, top}); + bonds += Bond("TJISING", "J", {site, top}); } } bonds["T"] = (std::complex)(1.0 + 0.2i); diff --git a/tests/blocks/electron/test_electron_matrix.cpp b/tests/blocks/electron/test_electron_matrix.cpp index 1f0727a6..f69d6eff 100644 --- a/tests/blocks/electron/test_electron_matrix.cpp +++ b/tests/blocks/electron/test_electron_matrix.cpp @@ -53,7 +53,7 @@ TEST_CASE("electron_matrix", "[electron]") { auto block = Electron(n_sites, n_up, n_dn); for (int i = 0; i < n_sites; ++i) - bondlist << Bond("HOP", "T", {i, (i + 1) % n_sites}); + bondlist += Bond("HOP", "T", {i, (i + 1) % n_sites}); bondlist["T"] = 1.0; bondlist["U"] = 5.0; auto H1 = matrix(bondlist, block, block); @@ -121,22 +121,24 @@ TEST_CASE("electron_matrix", "[electron]") { ////////////////////////////////// // Test two site exact solution - bondlist.clear(); - bondlist << Bond("HOP", "T", {0, 1}); - auto block2 = Electron(2, 1, 1); - for (int i = 0; i < 20; ++i) { - double U = 1.234 * i; - Log("electron_matrix: two-site exact solution test, U={}", U); - bondlist["T"] = 1.0; - bondlist["U"] = U; - double e0_exact = 0.5 * (U - sqrt(U * U + 16)); - auto H = matrix(bondlist, block2, block2); - REQUIRE(H.is_hermitian(1e-8)); - arma::Col eigs; - arma::eig_sym(eigs, H); - double e0 = eigs(0); - // printf("e0: %f, e0_exact: %f\n", e0, e0_exact); - REQUIRE(close(e0_exact, e0)); + { + BondList bondlist; + bondlist += Bond("HOP", "T", {0, 1}); + auto block2 = Electron(2, 1, 1); + for (int i = 0; i < 20; ++i) { + double U = 1.234 * i; + Log("electron_matrix: two-site exact solution test, U={}", U); + bondlist["T"] = 1.0; + bondlist["U"] = U; + double e0_exact = 0.5 * (U - sqrt(U * U + 16)); + auto H = matrix(bondlist, block2, block2); + REQUIRE(H.is_hermitian(1e-8)); + arma::Col eigs; + arma::eig_sym(eigs, H); + double e0 = eigs(0); + // printf("e0: %f, e0_exact: %f\n", e0, e0_exact); + REQUIRE(close(e0_exact, e0)); + } } /////////////////////////////////////////////////// @@ -150,10 +152,10 @@ TEST_CASE("electron_matrix", "[electron]") { // Create single particle matrix arma::Mat Hs(n_sites, n_sites, arma::fill::zeros); for (auto bond : bondlist) { - assert(bond.size() == 2); - int s1 = bond.site(0); - int s2 = bond.site(1); - auto name = bond.coupling_name(); + REQUIRE(bond.size() == 2); + int s1 = bond[0]; + int s2 = bond[1]; + auto name = bond.coupling().as(); Hs(s1, s2) = -bondlist[name].as(); Hs(s2, s1) = -bondlist[name].as(); } @@ -196,11 +198,11 @@ TEST_CASE("electron_matrix", "[electron]") { // Create single particle matrix for upspins arma::cx_mat Hs_up(n_sites, n_sites, arma::fill::zeros); - for (auto bond : bondlist.bonds_of_type("HOPUP")) { + for (auto bond : bonds_of_type("HOPUP", bondlist)) { assert(bond.size() == 2); - int s1 = bond.site(0); - int s2 = bond.site(1); - auto name = bond.coupling_name(); + int s1 = bond[0]; + int s2 = bond[1]; + auto name = bond.coupling().as(); Hs_up(s1, s2) = -bondlist[name].as(); Hs_up(s2, s1) = -conj(bondlist[name].as()); } @@ -209,11 +211,11 @@ TEST_CASE("electron_matrix", "[electron]") { // Create single particle matrix for dnspins arma::cx_mat Hs_dn(n_sites, n_sites, arma::fill::zeros); - for (auto bond : bondlist.bonds_of_type("HOPDN")) { + for (auto bond : bonds_of_type("HOPDN", bondlist)) { assert(bond.size() == 2); - int s1 = bond.site(0); - int s2 = bond.site(1); - auto name = bond.coupling_name(); + int s1 = bond[0]; + int s2 = bond[1]; + auto name = bond.coupling().as(); Hs_dn(s1, s2) = -bondlist[name].as(); Hs_dn(s2, s1) = -conj(bondlist[name].as()); } diff --git a/tests/blocks/electron/test_electron_raiselower.cpp b/tests/blocks/electron/test_electron_raiselower.cpp index cd2d740a..ebad2802 100644 --- a/tests/blocks/electron/test_electron_raiselower.cpp +++ b/tests/blocks/electron/test_electron_raiselower.cpp @@ -26,8 +26,8 @@ TEST_CASE("electron_raise_lower", "[electron]") { for (auto op_i_str : op_strs) { for (auto op_j_str : op_strs) { - auto op_i = Bond(op_i_str, i); - auto op_j = Bond(op_j_str, j); + auto op_i = Bond(op_i_str, 1.0, i); + auto op_j = Bond(op_j_str, 1.0, j); auto op_i_m = matrixC(op_i, block); auto op_j_m = matrixC(op_j, block); cx_mat anti_comm = op_i_m * op_j_m + op_j_m * op_i_m; @@ -61,20 +61,20 @@ TEST_CASE("electron_raise_lower", "[electron]") { // Check whether NUMBERUP agrees { - auto bond = Bond("NUMBERUP", i); + auto bond = Bond("NUMBERUP", 1.0, i); auto m = matrix(bond, block); - auto mi1 = matrix(Bond("CDAGUP", i), block); - auto mi2 = matrix(Bond("CUP", i), block); + auto mi1 = matrix(Bond("CDAGUP", 1.0, i), block); + auto mi2 = matrix(Bond("CUP", 1.0, i), block); auto mm = mi1 * mi2; REQUIRE(norm(mm - m) < 1e-12); } // Check whether NUMBERDN agrees { - auto bond = Bond("NUMBERDN", i); + auto bond = Bond("NUMBERDN", 1.0, i); auto m = matrix(bond, block); - auto mi1 = matrix(Bond("CDAGDN", i), block); - auto mi2 = matrix(Bond("CDN", i), block); + auto mi1 = matrix(Bond("CDAGDN", 1.0, i), block); + auto mi2 = matrix(Bond("CDN", 1.0, i), block); auto mm = mi1 * mi2; REQUIRE(norm(mm - m) < 1e-12); } @@ -84,43 +84,43 @@ TEST_CASE("electron_raise_lower", "[electron]") { // Check whether HOPUP agrees { - auto bond = Bond("HOPUP", {i, j}); + auto bond = Bond("HOPUP", 1.0, {i, j}); auto m = matrix(bond, block); - auto m1 = matrix(Bond("CDAGUP", i), block); - auto m2 = matrix(Bond("CUP", j), block); - auto m3 = matrix(Bond("CDAGUP", j), block); - auto m4 = matrix(Bond("CUP", i), block); + auto m1 = matrix(Bond("CDAGUP", 1.0, i), block); + auto m2 = matrix(Bond("CUP", 1.0, j), block); + auto m3 = matrix(Bond("CDAGUP", 1.0, j), block); + auto m4 = matrix(Bond("CUP", 1.0, i), block); auto mm = -m1 * m2 - m3 * m4; REQUIRE(norm(mm - m) < 1e-12); } // Check whether HOPDN agrees { - auto bond = Bond("HOPDN", {i, j}); + auto bond = Bond("HOPDN", 1.0, {i, j}); auto m = matrix(bond, block); - auto m1 = matrix(Bond("CDAGDN", i), block); - auto m2 = matrix(Bond("CDN", j), block); - auto m3 = matrix(Bond("CDAGDN", j), block); - auto m4 = matrix(Bond("CDN", i), block); + auto m1 = matrix(Bond("CDAGDN", 1.0, i), block); + auto m2 = matrix(Bond("CDN", 1.0, j), block); + auto m3 = matrix(Bond("CDAGDN", 1.0, j), block); + auto m4 = matrix(Bond("CDN", 1.0, i), block); auto mm = -m1 * m2 - m3 * m4; REQUIRE(norm(mm - m) < 1e-12); } // Check whether ISING agrees { - auto bond = Bond("ISING", {i, j}); + auto bond = Bond("ISING", 1.0, {i, j}); auto m = matrix(bond, block); - auto mi1 = matrix(Bond("CDAGUP", i), block); - auto mi2 = matrix(Bond("CUP", i), block); - auto mi3 = matrix(Bond("CDAGDN", i), block); - auto mi4 = matrix(Bond("CDN", i), block); + auto mi1 = matrix(Bond("CDAGUP", 1.0, i), block); + auto mi2 = matrix(Bond("CUP", 1.0, i), block); + auto mi3 = matrix(Bond("CDAGDN", 1.0, i), block); + auto mi4 = matrix(Bond("CDN", 1.0, i), block); auto szi = 0.5 * mi1 * mi2 - 0.5 * mi3 * mi4; - auto mj1 = matrix(Bond("CDAGUP", j), block); - auto mj2 = matrix(Bond("CUP", j), block); - auto mj3 = matrix(Bond("CDAGDN", j), block); - auto mj4 = matrix(Bond("CDN", j), block); + auto mj1 = matrix(Bond("CDAGUP", 1.0, j), block); + auto mj2 = matrix(Bond("CUP", 1.0, j), block); + auto mj3 = matrix(Bond("CDAGDN", 1.0, j), block); + auto mj4 = matrix(Bond("CDN", 1.0, j), block); auto szj = 0.5 * mj1 * mj2 - 0.5 * mj3 * mj4; auto mm = szi * szj; @@ -129,20 +129,20 @@ TEST_CASE("electron_raise_lower", "[electron]") { // Check whether EXCHANGE agrees { - auto bond = Bond("EXCHANGE", {i, j}); + auto bond = Bond("EXCHANGE", 1.0, {i, j}); auto m = matrix(bond, block); - auto mi1 = matrix(Bond("CDAGUP", i), block); - auto mi2 = matrix(Bond("CUP", i), block); - auto mi3 = matrix(Bond("CDAGDN", i), block); - auto mi4 = matrix(Bond("CDN", i), block); + auto mi1 = matrix(Bond("CDAGUP", 1.0, i), block); + auto mi2 = matrix(Bond("CUP", 1.0, i), block); + auto mi3 = matrix(Bond("CDAGDN", 1.0, i), block); + auto mi4 = matrix(Bond("CDN", 1.0, i), block); auto spi = mi1 * mi4; auto smi = mi3 * mi2; - auto mj1 = matrix(Bond("CDAGUP", j), block); - auto mj2 = matrix(Bond("CUP", j), block); - auto mj3 = matrix(Bond("CDAGDN", j), block); - auto mj4 = matrix(Bond("CDN", j), block); + auto mj1 = matrix(Bond("CDAGUP", 1.0, j), block); + auto mj2 = matrix(Bond("CUP", 1.0, j), block); + auto mj3 = matrix(Bond("CDAGDN", 1.0, j), block); + auto mj4 = matrix(Bond("CDN", 1.0, j), block); auto spj = mj1 * mj4; auto smj = mj3 * mj2; diff --git a/tests/blocks/electron/testcases_electron.cpp b/tests/blocks/electron/testcases_electron.cpp index cde0566d..3c4fbefd 100644 --- a/tests/blocks/electron/testcases_electron.cpp +++ b/tests/blocks/electron/testcases_electron.cpp @@ -6,7 +6,7 @@ BondList get_linear_chain(int64_t n_sites, double t, double U) { // Create model BondList bondlist; for (int64_t s = 0; s < n_sites; ++s) - bondlist << Bond("HOP", "T", {s, (s + 1) % n_sites}); + bondlist += Bond("HOP", "T", {s, (s + 1) % n_sites}); bondlist["T"] = t; bondlist["U"] = U; return bondlist; @@ -16,9 +16,9 @@ BondList get_linear_chain_hb(int64_t n_sites, double J) { // Create model BondList bondlist; // for (int64_t s = 0; s < n_sites; ++s) - // bondlist << Bond("HOP", "T", {s, (s + 1) % n_sites}); + // bondlist += Bond("HOP", "T", {s, (s + 1) % n_sites}); for (int64_t s = 0; s < n_sites; ++s) - bondlist << Bond("EXCHANGE", "J", {s, (s + 1) % n_sites}); + bondlist += Bond("EXCHANGE", "J", {s, (s + 1) % n_sites}); bondlist["J"] = J; return bondlist; } @@ -83,9 +83,9 @@ get_cyclic_group_irreps_mult(int64_t n_sites) { BondList heisenberg_triangle() { BondList bondlist; - bondlist << Bond("HEISENBERG", "J", {0, 1}); - bondlist << Bond("HEISENBERG", "J", {1, 2}); - bondlist << Bond("HEISENBERG", "J", {2, 0}); + bondlist += Bond("HB", "J", {0, 1}); + bondlist += Bond("HB", "J", {1, 2}); + bondlist += Bond("HB", "J", {2, 0}); bondlist["J"] = 1.0; return bondlist; } @@ -101,7 +101,7 @@ BondList heisenberg_alltoall(int64_t n_sites) { ss << "J" << s1 << "_" << s2; std::string name = ss.str(); double value = distribution(generator); - bondlist << Bond("HEISENBERG", name, {s1, s2}); + bondlist += Bond("HB", name, {s1, s2}); bondlist[name] = value; } return bondlist; @@ -110,99 +110,99 @@ BondList heisenberg_alltoall(int64_t n_sites) { BondList heisenberg_kagome15() { BondList bondlist; - bondlist << Bond("HEISENBERG", "J", {0, 1}); - bondlist << Bond("HEISENBERG", "J", {0, 5}); - bondlist << Bond("HEISENBERG", "J", {1, 2}); - bondlist << Bond("HEISENBERG", "J", {1, 6}); - bondlist << Bond("HEISENBERG", "J", {2, 3}); - bondlist << Bond("HEISENBERG", "J", {2, 6}); - bondlist << Bond("HEISENBERG", "J", {2, 10}); - bondlist << Bond("HEISENBERG", "J", {3, 4}); - bondlist << Bond("HEISENBERG", "J", {3, 10}); - bondlist << Bond("HEISENBERG", "J", {3, 14}); - bondlist << Bond("HEISENBERG", "J", {4, 5}); - bondlist << Bond("HEISENBERG", "J", {4, 14}); - bondlist << Bond("HEISENBERG", "J", {6, 7}); - bondlist << Bond("HEISENBERG", "J", {7, 8}); - bondlist << Bond("HEISENBERG", "J", {8, 9}); - bondlist << Bond("HEISENBERG", "J", {9, 10}); - bondlist << Bond("HEISENBERG", "J", {9, 11}); - bondlist << Bond("HEISENBERG", "J", {10, 11}); - bondlist << Bond("HEISENBERG", "J", {11, 12}); - bondlist << Bond("HEISENBERG", "J", {12, 13}); - bondlist << Bond("HEISENBERG", "J", {13, 14}); + bondlist += Bond("HB", "J", {0, 1}); + bondlist += Bond("HB", "J", {0, 5}); + bondlist += Bond("HB", "J", {1, 2}); + bondlist += Bond("HB", "J", {1, 6}); + bondlist += Bond("HB", "J", {2, 3}); + bondlist += Bond("HB", "J", {2, 6}); + bondlist += Bond("HB", "J", {2, 10}); + bondlist += Bond("HB", "J", {3, 4}); + bondlist += Bond("HB", "J", {3, 10}); + bondlist += Bond("HB", "J", {3, 14}); + bondlist += Bond("HB", "J", {4, 5}); + bondlist += Bond("HB", "J", {4, 14}); + bondlist += Bond("HB", "J", {6, 7}); + bondlist += Bond("HB", "J", {7, 8}); + bondlist += Bond("HB", "J", {8, 9}); + bondlist += Bond("HB", "J", {9, 10}); + bondlist += Bond("HB", "J", {9, 11}); + bondlist += Bond("HB", "J", {10, 11}); + bondlist += Bond("HB", "J", {11, 12}); + bondlist += Bond("HB", "J", {12, 13}); + bondlist += Bond("HB", "J", {13, 14}); bondlist["J"] = 1.0; return bondlist; } BondList heisenberg_kagome39() { BondList bondlist; - bondlist << Bond("HEISENBERG", "J", {0, 1}); - bondlist << Bond("HEISENBERG", "J", {0, 5}); - bondlist << Bond("HEISENBERG", "J", {0, 27}); - bondlist << Bond("HEISENBERG", "J", {0, 28}); - bondlist << Bond("HEISENBERG", "J", {1, 2}); - bondlist << Bond("HEISENBERG", "J", {1, 6}); - bondlist << Bond("HEISENBERG", "J", {1, 28}); - bondlist << Bond("HEISENBERG", "J", {2, 3}); - bondlist << Bond("HEISENBERG", "J", {2, 6}); - bondlist << Bond("HEISENBERG", "J", {2, 10}); - bondlist << Bond("HEISENBERG", "J", {3, 4}); - bondlist << Bond("HEISENBERG", "J", {3, 10}); - bondlist << Bond("HEISENBERG", "J", {3, 14}); - bondlist << Bond("HEISENBERG", "J", {4, 5}); - bondlist << Bond("HEISENBERG", "J", {4, 14}); - bondlist << Bond("HEISENBERG", "J", {4, 38}); - bondlist << Bond("HEISENBERG", "J", {5, 27}); - bondlist << Bond("HEISENBERG", "J", {5, 38}); - bondlist << Bond("HEISENBERG", "J", {6, 7}); - bondlist << Bond("HEISENBERG", "J", {6, 29}); - bondlist << Bond("HEISENBERG", "J", {7, 8}); - bondlist << Bond("HEISENBERG", "J", {7, 19}); - bondlist << Bond("HEISENBERG", "J", {7, 29}); - bondlist << Bond("HEISENBERG", "J", {8, 9}); - bondlist << Bond("HEISENBERG", "J", {8, 15}); - bondlist << Bond("HEISENBERG", "J", {8, 19}); - bondlist << Bond("HEISENBERG", "J", {9, 10}); - bondlist << Bond("HEISENBERG", "J", {9, 11}); - bondlist << Bond("HEISENBERG", "J", {9, 15}); - bondlist << Bond("HEISENBERG", "J", {10, 11}); - bondlist << Bond("HEISENBERG", "J", {11, 12}); - bondlist << Bond("HEISENBERG", "J", {11, 18}); - bondlist << Bond("HEISENBERG", "J", {12, 13}); - bondlist << Bond("HEISENBERG", "J", {12, 18}); - bondlist << Bond("HEISENBERG", "J", {12, 26}); - bondlist << Bond("HEISENBERG", "J", {13, 14}); - bondlist << Bond("HEISENBERG", "J", {13, 26}); - bondlist << Bond("HEISENBERG", "J", {13, 37}); - bondlist << Bond("HEISENBERG", "J", {14, 37}); - bondlist << Bond("HEISENBERG", "J", {15, 16}); - bondlist << Bond("HEISENBERG", "J", {15, 22}); - bondlist << Bond("HEISENBERG", "J", {16, 17}); - bondlist << Bond("HEISENBERG", "J", {16, 22}); - bondlist << Bond("HEISENBERG", "J", {16, 33}); - bondlist << Bond("HEISENBERG", "J", {17, 18}); - bondlist << Bond("HEISENBERG", "J", {17, 23}); - bondlist << Bond("HEISENBERG", "J", {17, 33}); - bondlist << Bond("HEISENBERG", "J", {18, 23}); - bondlist << Bond("HEISENBERG", "J", {19, 20}); - bondlist << Bond("HEISENBERG", "J", {19, 30}); - bondlist << Bond("HEISENBERG", "J", {20, 21}); - bondlist << Bond("HEISENBERG", "J", {20, 30}); - bondlist << Bond("HEISENBERG", "J", {20, 31}); - bondlist << Bond("HEISENBERG", "J", {21, 22}); - bondlist << Bond("HEISENBERG", "J", {21, 31}); - bondlist << Bond("HEISENBERG", "J", {21, 32}); - bondlist << Bond("HEISENBERG", "J", {22, 32}); - bondlist << Bond("HEISENBERG", "J", {23, 24}); - bondlist << Bond("HEISENBERG", "J", {23, 34}); - bondlist << Bond("HEISENBERG", "J", {24, 25}); - bondlist << Bond("HEISENBERG", "J", {24, 34}); - bondlist << Bond("HEISENBERG", "J", {24, 35}); - bondlist << Bond("HEISENBERG", "J", {25, 26}); - bondlist << Bond("HEISENBERG", "J", {25, 35}); - bondlist << Bond("HEISENBERG", "J", {25, 36}); - bondlist << Bond("HEISENBERG", "J", {26, 36}); + bondlist += Bond("HB", "J", {0, 1}); + bondlist += Bond("HB", "J", {0, 5}); + bondlist += Bond("HB", "J", {0, 27}); + bondlist += Bond("HB", "J", {0, 28}); + bondlist += Bond("HB", "J", {1, 2}); + bondlist += Bond("HB", "J", {1, 6}); + bondlist += Bond("HB", "J", {1, 28}); + bondlist += Bond("HB", "J", {2, 3}); + bondlist += Bond("HB", "J", {2, 6}); + bondlist += Bond("HB", "J", {2, 10}); + bondlist += Bond("HB", "J", {3, 4}); + bondlist += Bond("HB", "J", {3, 10}); + bondlist += Bond("HB", "J", {3, 14}); + bondlist += Bond("HB", "J", {4, 5}); + bondlist += Bond("HB", "J", {4, 14}); + bondlist += Bond("HB", "J", {4, 38}); + bondlist += Bond("HB", "J", {5, 27}); + bondlist += Bond("HB", "J", {5, 38}); + bondlist += Bond("HB", "J", {6, 7}); + bondlist += Bond("HB", "J", {6, 29}); + bondlist += Bond("HB", "J", {7, 8}); + bondlist += Bond("HB", "J", {7, 19}); + bondlist += Bond("HB", "J", {7, 29}); + bondlist += Bond("HB", "J", {8, 9}); + bondlist += Bond("HB", "J", {8, 15}); + bondlist += Bond("HB", "J", {8, 19}); + bondlist += Bond("HB", "J", {9, 10}); + bondlist += Bond("HB", "J", {9, 11}); + bondlist += Bond("HB", "J", {9, 15}); + bondlist += Bond("HB", "J", {10, 11}); + bondlist += Bond("HB", "J", {11, 12}); + bondlist += Bond("HB", "J", {11, 18}); + bondlist += Bond("HB", "J", {12, 13}); + bondlist += Bond("HB", "J", {12, 18}); + bondlist += Bond("HB", "J", {12, 26}); + bondlist += Bond("HB", "J", {13, 14}); + bondlist += Bond("HB", "J", {13, 26}); + bondlist += Bond("HB", "J", {13, 37}); + bondlist += Bond("HB", "J", {14, 37}); + bondlist += Bond("HB", "J", {15, 16}); + bondlist += Bond("HB", "J", {15, 22}); + bondlist += Bond("HB", "J", {16, 17}); + bondlist += Bond("HB", "J", {16, 22}); + bondlist += Bond("HB", "J", {16, 33}); + bondlist += Bond("HB", "J", {17, 18}); + bondlist += Bond("HB", "J", {17, 23}); + bondlist += Bond("HB", "J", {17, 33}); + bondlist += Bond("HB", "J", {18, 23}); + bondlist += Bond("HB", "J", {19, 20}); + bondlist += Bond("HB", "J", {19, 30}); + bondlist += Bond("HB", "J", {20, 21}); + bondlist += Bond("HB", "J", {20, 30}); + bondlist += Bond("HB", "J", {20, 31}); + bondlist += Bond("HB", "J", {21, 22}); + bondlist += Bond("HB", "J", {21, 31}); + bondlist += Bond("HB", "J", {21, 32}); + bondlist += Bond("HB", "J", {22, 32}); + bondlist += Bond("HB", "J", {23, 24}); + bondlist += Bond("HB", "J", {23, 34}); + bondlist += Bond("HB", "J", {24, 25}); + bondlist += Bond("HB", "J", {24, 34}); + bondlist += Bond("HB", "J", {24, 35}); + bondlist += Bond("HB", "J", {25, 26}); + bondlist += Bond("HB", "J", {25, 35}); + bondlist += Bond("HB", "J", {25, 36}); + bondlist += Bond("HB", "J", {26, 36}); bondlist["J"] = 1.0; return bondlist; } @@ -218,7 +218,7 @@ BondList freefermion_alltoall(int64_t n_sites) { ss << "T" << s1 << "_" << s2; std::string name = ss.str(); double value = distribution(generator); - bondlist << Bond("HOP", name, {s1, s2}); + bondlist += Bond("HOP", name, {s1, s2}); bondlist[name] = value; } return bondlist; @@ -238,7 +238,7 @@ BondList freefermion_alltoall_complex_updn(int64_t n_sites) { std::string name_up = ss_up.str(); complex value_up = complex(distribution(generator), distribution(generator)); - bondlist << Bond("HOPUP", name_up, {s1, s2}); + bondlist += Bond("HOPUP", name_up, {s1, s2}); bondlist[name_up] = value_up; // Hopping on dnspins @@ -247,7 +247,7 @@ BondList freefermion_alltoall_complex_updn(int64_t n_sites) { std::string name_dn = ss_dn.str(); complex value_dn = complex(distribution(generator), distribution(generator)); - bondlist << Bond("HOPDN", name_dn, {s1, s2}); + bondlist += Bond("HOPDN", name_dn, {s1, s2}); bondlist[name_dn] = value_dn; } return bondlist; @@ -268,18 +268,18 @@ std::tuple> randomAlltoAll4NoU() { bondlist["T23"] = -2; bondlist["J23"] = 1; - bondlist << Bond("HOP", "T01", {0, 1}); - bondlist << Bond("HOP", "T02", {0, 2}); - bondlist << Bond("HOP", "T03", {0, 3}); - bondlist << Bond("HOP", "T12", {1, 2}); - bondlist << Bond("HOP", "T13", {1, 3}); - bondlist << Bond("HOP", "T23", {2, 3}); - bondlist << Bond("HEISENBERG", "J01", {0, 1}); - bondlist << Bond("HEISENBERG", "J02", {0, 2}); - bondlist << Bond("HEISENBERG", "J03", {0, 3}); - bondlist << Bond("HEISENBERG", "J12", {1, 2}); - bondlist << Bond("HEISENBERG", "J13", {1, 3}); - bondlist << Bond("HEISENBERG", "J23", {2, 3}); + bondlist += Bond("HOP", "T01", {0, 1}); + bondlist += Bond("HOP", "T02", {0, 2}); + bondlist += Bond("HOP", "T03", {0, 3}); + bondlist += Bond("HOP", "T12", {1, 2}); + bondlist += Bond("HOP", "T13", {1, 3}); + bondlist += Bond("HOP", "T23", {2, 3}); + bondlist += Bond("HB", "J01", {0, 1}); + bondlist += Bond("HB", "J02", {0, 2}); + bondlist += Bond("HB", "J03", {0, 3}); + bondlist += Bond("HB", "J12", {1, 2}); + bondlist += Bond("HB", "J13", {1, 3}); + bondlist += Bond("HB", "J23", {2, 3}); arma::Col eigs = {-17.035603173216636, -16.054529653295518, @@ -558,18 +558,18 @@ std::tuple> randomAlltoAll4() { bondlist["T23"] = 0; bondlist["J23"] = -4; - bondlist << Bond("HOP", "T01", {0, 1}); - bondlist << Bond("HOP", "T02", {0, 2}); - bondlist << Bond("HOP", "T03", {0, 3}); - bondlist << Bond("HOP", "T12", {1, 2}); - bondlist << Bond("HOP", "T13", {1, 3}); - bondlist << Bond("HOP", "T23", {2, 3}); - bondlist << Bond("HEISENBERG", "J01", {0, 1}); - bondlist << Bond("HEISENBERG", "J02", {0, 2}); - bondlist << Bond("HEISENBERG", "J03", {0, 3}); - bondlist << Bond("HEISENBERG", "J12", {1, 2}); - bondlist << Bond("HEISENBERG", "J13", {1, 3}); - bondlist << Bond("HEISENBERG", "J23", {2, 3}); + bondlist += Bond("HOP", "T01", {0, 1}); + bondlist += Bond("HOP", "T02", {0, 2}); + bondlist += Bond("HOP", "T03", {0, 3}); + bondlist += Bond("HOP", "T12", {1, 2}); + bondlist += Bond("HOP", "T13", {1, 3}); + bondlist += Bond("HOP", "T23", {2, 3}); + bondlist += Bond("HB", "J01", {0, 1}); + bondlist += Bond("HB", "J02", {0, 2}); + bondlist += Bond("HB", "J03", {0, 3}); + bondlist += Bond("HB", "J12", {1, 2}); + bondlist += Bond("HB", "J13", {1, 3}); + bondlist += Bond("HB", "J23", {2, 3}); arma::Col eigs = { -12.270231830055396, -12.270231830055389, -10.733666336755952, @@ -670,12 +670,12 @@ BondList randomAlltoAll3() { bondlist["J02"] = -1; bondlist["T12"] = -5; bondlist["J12"] = -3; - bondlist << Bond("HUBBARDHOP", "T01", {0, 1}); - bondlist << Bond("HUBBARDHOP", "T02", {0, 2}); - bondlist << Bond("HUBBARDHOP", "T12", {1, 2}); - bondlist << Bond("HEISENBERG", "J01", {0, 1}); - bondlist << Bond("HEISENBERG", "J02", {0, 2}); - bondlist << Bond("HEISENBERG", "J12", {1, 2}); + bondlist += Bond("HOP", "T01", {0, 1}); + bondlist += Bond("HOP", "T02", {0, 2}); + bondlist += Bond("HOP", "T12", {1, 2}); + bondlist += Bond("HB", "J01", {0, 1}); + bondlist += Bond("HB", "J02", {0, 2}); + bondlist += Bond("HB", "J12", {1, 2}); return bondlist; } @@ -683,22 +683,22 @@ BondList square2x2(double t, double J) { BondList bondlist; bondlist["T"] = t; bondlist["J"] = J; - bondlist << Bond("HUBBARDHOP", "T", {0, 1}); - bondlist << Bond("HUBBARDHOP", "T", {1, 0}); - bondlist << Bond("HUBBARDHOP", "T", {2, 3}); - bondlist << Bond("HUBBARDHOP", "T", {3, 2}); - bondlist << Bond("HUBBARDHOP", "T", {0, 2}); - bondlist << Bond("HUBBARDHOP", "T", {2, 0}); - bondlist << Bond("HUBBARDHOP", "T", {1, 3}); - bondlist << Bond("HUBBARDHOP", "T", {3, 1}); - bondlist << Bond("HEISENBERG", "J", {0, 1}); - bondlist << Bond("HEISENBERG", "J", {1, 0}); - bondlist << Bond("HEISENBERG", "J", {2, 3}); - bondlist << Bond("HEISENBERG", "J", {3, 2}); - bondlist << Bond("HEISENBERG", "J", {0, 2}); - bondlist << Bond("HEISENBERG", "J", {2, 0}); - bondlist << Bond("HEISENBERG", "J", {1, 3}); - bondlist << Bond("HEISENBERG", "J", {3, 1}); + bondlist += Bond("HOP", "T", {0, 1}); + bondlist += Bond("HOP", "T", {1, 0}); + bondlist += Bond("HOP", "T", {2, 3}); + bondlist += Bond("HOP", "T", {3, 2}); + bondlist += Bond("HOP", "T", {0, 2}); + bondlist += Bond("HOP", "T", {2, 0}); + bondlist += Bond("HOP", "T", {1, 3}); + bondlist += Bond("HOP", "T", {3, 1}); + bondlist += Bond("HB", "J", {0, 1}); + bondlist += Bond("HB", "J", {1, 0}); + bondlist += Bond("HB", "J", {2, 3}); + bondlist += Bond("HB", "J", {3, 2}); + bondlist += Bond("HB", "J", {0, 2}); + bondlist += Bond("HB", "J", {2, 0}); + bondlist += Bond("HB", "J", {1, 3}); + bondlist += Bond("HB", "J", {3, 1}); return bondlist; } @@ -706,42 +706,42 @@ BondList square3x3(double t, double J) { BondList bondlist; bondlist["T"] = t; bondlist["J"] = J; - bondlist << Bond("HOP", "T", {0, 1}); - bondlist << Bond("HOP", "T", {1, 2}); - bondlist << Bond("HOP", "T", {2, 0}); - bondlist << Bond("HOP", "T", {3, 4}); - bondlist << Bond("HOP", "T", {4, 5}); - bondlist << Bond("HOP", "T", {5, 3}); - bondlist << Bond("HOP", "T", {6, 7}); - bondlist << Bond("HOP", "T", {7, 8}); - bondlist << Bond("HOP", "T", {8, 6}); - bondlist << Bond("HOP", "T", {0, 3}); - bondlist << Bond("HOP", "T", {3, 6}); - bondlist << Bond("HOP", "T", {6, 0}); - bondlist << Bond("HOP", "T", {1, 4}); - bondlist << Bond("HOP", "T", {4, 7}); - bondlist << Bond("HOP", "T", {7, 1}); - bondlist << Bond("HOP", "T", {2, 5}); - bondlist << Bond("HOP", "T", {5, 8}); - bondlist << Bond("HOP", "T", {8, 2}); - bondlist << Bond("HEISENBERG", "J", {0, 1}); - bondlist << Bond("HEISENBERG", "J", {1, 2}); - bondlist << Bond("HEISENBERG", "J", {2, 0}); - bondlist << Bond("HEISENBERG", "J", {3, 4}); - bondlist << Bond("HEISENBERG", "J", {4, 5}); - bondlist << Bond("HEISENBERG", "J", {5, 3}); - bondlist << Bond("HEISENBERG", "J", {6, 7}); - bondlist << Bond("HEISENBERG", "J", {7, 8}); - bondlist << Bond("HEISENBERG", "J", {8, 6}); - bondlist << Bond("HEISENBERG", "J", {0, 3}); - bondlist << Bond("HEISENBERG", "J", {3, 6}); - bondlist << Bond("HEISENBERG", "J", {6, 0}); - bondlist << Bond("HEISENBERG", "J", {1, 4}); - bondlist << Bond("HEISENBERG", "J", {4, 7}); - bondlist << Bond("HEISENBERG", "J", {7, 1}); - bondlist << Bond("HEISENBERG", "J", {2, 5}); - bondlist << Bond("HEISENBERG", "J", {5, 8}); - bondlist << Bond("HEISENBERG", "J", {8, 2}); + bondlist += Bond("HOP", "T", {0, 1}); + bondlist += Bond("HOP", "T", {1, 2}); + bondlist += Bond("HOP", "T", {2, 0}); + bondlist += Bond("HOP", "T", {3, 4}); + bondlist += Bond("HOP", "T", {4, 5}); + bondlist += Bond("HOP", "T", {5, 3}); + bondlist += Bond("HOP", "T", {6, 7}); + bondlist += Bond("HOP", "T", {7, 8}); + bondlist += Bond("HOP", "T", {8, 6}); + bondlist += Bond("HOP", "T", {0, 3}); + bondlist += Bond("HOP", "T", {3, 6}); + bondlist += Bond("HOP", "T", {6, 0}); + bondlist += Bond("HOP", "T", {1, 4}); + bondlist += Bond("HOP", "T", {4, 7}); + bondlist += Bond("HOP", "T", {7, 1}); + bondlist += Bond("HOP", "T", {2, 5}); + bondlist += Bond("HOP", "T", {5, 8}); + bondlist += Bond("HOP", "T", {8, 2}); + bondlist += Bond("HB", "J", {0, 1}); + bondlist += Bond("HB", "J", {1, 2}); + bondlist += Bond("HB", "J", {2, 0}); + bondlist += Bond("HB", "J", {3, 4}); + bondlist += Bond("HB", "J", {4, 5}); + bondlist += Bond("HB", "J", {5, 3}); + bondlist += Bond("HB", "J", {6, 7}); + bondlist += Bond("HB", "J", {7, 8}); + bondlist += Bond("HB", "J", {8, 6}); + bondlist += Bond("HB", "J", {0, 3}); + bondlist += Bond("HB", "J", {3, 6}); + bondlist += Bond("HB", "J", {6, 0}); + bondlist += Bond("HB", "J", {1, 4}); + bondlist += Bond("HB", "J", {4, 7}); + bondlist += Bond("HB", "J", {7, 1}); + bondlist += Bond("HB", "J", {2, 5}); + bondlist += Bond("HB", "J", {5, 8}); + bondlist += Bond("HB", "J", {8, 2}); return bondlist; } diff --git a/tests/blocks/electron_symmetric/test_electron_symmetric_apply.cpp b/tests/blocks/electron_symmetric/test_electron_symmetric_apply.cpp index d28dc052..a5becd66 100644 --- a/tests/blocks/electron_symmetric/test_electron_symmetric_apply.cpp +++ b/tests/blocks/electron_symmetric/test_electron_symmetric_apply.cpp @@ -133,7 +133,7 @@ TEST_CASE("electron_symmetric_apply", "[electron]") { "electron_symmetric_apply: Hubbard 3x3 triangular (+ Heisenberg terms)"); auto bondlist_hb = bondlist; for (auto bond : bondlist) { - bondlist_hb << Bond("HB", "J", {bond[0], bond[1]}); + bondlist_hb += Bond("HB", "J", {bond[0], bond[1]}); } bondlist_hb["J"] = 0.4; test_electron_symmetric_apply(bondlist_hb, space_group, irreps); diff --git a/tests/blocks/electron_symmetric/test_electron_symmetric_matrix.cpp b/tests/blocks/electron_symmetric/test_electron_symmetric_matrix.cpp index 8a86692c..780ec59b 100644 --- a/tests/blocks/electron_symmetric/test_electron_symmetric_matrix.cpp +++ b/tests/blocks/electron_symmetric/test_electron_symmetric_matrix.cpp @@ -256,7 +256,7 @@ TEST_CASE("electron_symmetric_matrix", "[electron]") { Log("electron_symmetric_matrix: Hubbard 3x3 triangular(+ Heisenberg terms)"); auto bondlist_hb = bondlist; for (auto bond : bondlist) { - bondlist_hb << Bond("HB", "J", {bond[0], bond[1]}); + bondlist_hb += Bond("HB", "J", {bond[0], bond[1]}); } bondlist_hb["J"] = 0.4; test_electron_symmetric_spectra(bondlist_hb, space_group, irreps, diff --git a/tests/blocks/spinhalf/test_spinhalf_apply.cpp b/tests/blocks/spinhalf/test_spinhalf_apply.cpp index 4e641db9..e0d94ea5 100644 --- a/tests/blocks/spinhalf/test_spinhalf_apply.cpp +++ b/tests/blocks/spinhalf/test_spinhalf_apply.cpp @@ -12,8 +12,7 @@ using namespace xdiag; -void test_apply(BondList bonds) { - int N = bonds.n_sites(); +void test_apply(int N, BondList bonds) { for (int nup = 1; nup <= N; ++nup) { auto block = Spinhalf(N, nup); auto H = matrix(bonds, block, block); @@ -46,13 +45,13 @@ TEST_CASE("spinhalf_apply", "[spinhalf]") { Log.out("spinhalf_apply: Heisenberg chain apply test, J=1.0, N=2,..,6"); for (int N = 2; N <= 6; ++N) { auto bonds = HBchain(N, 1.0); - test_apply(bonds); + test_apply(N, bonds); } Log.out("spinhalf_apply: Heisenberg alltoall apply test, N=2,..,6"); for (int N = 2; N <= 6; ++N) { auto bonds = HB_alltoall(N); - test_apply(bonds); + test_apply(N, bonds); } Log.out("spinhalf_apply: Heisenberg all-to-all Sz <-> NoSz comparison"); diff --git a/tests/blocks/spinhalf/test_spinhalf_matrix.cpp b/tests/blocks/spinhalf/test_spinhalf_matrix.cpp index aa1f19b9..b0057d05 100644 --- a/tests/blocks/spinhalf/test_spinhalf_matrix.cpp +++ b/tests/blocks/spinhalf/test_spinhalf_matrix.cpp @@ -140,21 +140,21 @@ TEST_CASE("spinhalf_matrix", "[spinhalf]") { for (int j = 0; j < n_sites; ++j) { BondList sp_i_m; - sp_i_m << Bond("S+", i); + sp_i_m += Bond("S+", 1.0, i); auto sp_i_m_mat = matrix(sp_i_m, blockm, block); auto sp_i_mat = matrix(sp_i_m, block_raw, block_raw); BondList sm_j_m; - sm_j_m << Bond("S-", j); + sm_j_m += Bond("S-", 1.0, j); auto sm_j_m_mat = matrix(sm_j_m, block, blockm); auto sm_j_mat = matrix(sm_j_m, block_raw, block_raw); BondList sp_i_p; - sp_i_p << Bond("S+", i); + sp_i_p += Bond("S+", 1.0, i); auto sp_i_p_mat = matrix(sp_i_p, block, blockp); BondList sm_j_p; - sm_j_p << Bond("S-", j); + sm_j_p += Bond("S-", 1.0, j); auto sm_j_p_mat = matrix(sm_j_p, blockp, block); auto C1 = sp_i_m_mat * sm_j_m_mat; @@ -166,7 +166,7 @@ TEST_CASE("spinhalf_matrix", "[spinhalf]") { if (i == j) { BondList sz; - sz << Bond("SZ", i); + sz += Bond("SZ", 1.0, i); auto sz_mat = matrix(sz, block, block); auto sz_matr = matrix(sz, block_raw, block_raw); REQUIRE(close(comm, arma::mat(2.0 * sz_mat))); diff --git a/tests/blocks/spinhalf/testcases_spinhalf.cpp b/tests/blocks/spinhalf/testcases_spinhalf.cpp index a3cb4591..d08c1246 100644 --- a/tests/blocks/spinhalf/testcases_spinhalf.cpp +++ b/tests/blocks/spinhalf/testcases_spinhalf.cpp @@ -11,10 +11,10 @@ BondList HBchain(int64_t n_sites, double J1, double J2) { bondlist["J2"] = J2; for (int64_t s = 0; s < n_sites; ++s) { - bondlist << Bond("HB", "J1", {s, (s + 1) % n_sites}); + bondlist += Bond("HB", "J1", {s, (s + 1) % n_sites}); } for (int64_t s = 0; s < n_sites; ++s) { - bondlist << Bond("HB", "J2", {s, (s + 2) % n_sites}); + bondlist += Bond("HB", "J2", {s, (s + 2) % n_sites}); } return bondlist; } @@ -107,7 +107,7 @@ BondList HB_alltoall(int64_t n_sites) { ss << "J" << s1 << "_" << s2; std::string name = ss.str(); double value = distribution(generator); - bondlist << Bond("HB", name, {s1, s2}); + bondlist += Bond("HB", name, {s1, s2}); bondlist[name] = value; } return bondlist; @@ -116,78 +116,78 @@ BondList HB_alltoall(int64_t n_sites) { std::tuple triangular_12_complex(int64_t nup, double eta) { BondList bonds; - bonds << Bond("ISING", "Jz", {0, 5}); - bonds << Bond("ISING", "Jz", {8, 2}); - bonds << Bond("ISING", "Jz", {2, 7}); - bonds << Bond("ISING", "Jz", {1, 4}); - bonds << Bond("ISING", "Jz", {4, 8}); - bonds << Bond("ISING", "Jz", {6, 10}); - bonds << Bond("ISING", "Jz", {10, 0}); - bonds << Bond("ISING", "Jz", {5, 9}); - bonds << Bond("ISING", "Jz", {9, 3}); - bonds << Bond("ISING", "Jz", {3, 6}); - bonds << Bond("ISING", "Jz", {7, 11}); - bonds << Bond("ISING", "Jz", {11, 1}); - bonds << Bond("ISING", "Jz", {0, 8}); - bonds << Bond("ISING", "Jz", {8, 6}); - bonds << Bond("ISING", "Jz", {2, 10}); - bonds << Bond("ISING", "Jz", {1, 9}); - bonds << Bond("ISING", "Jz", {4, 3}); - bonds << Bond("ISING", "Jz", {6, 1}); - bonds << Bond("ISING", "Jz", {10, 4}); - bonds << Bond("ISING", "Jz", {5, 2}); - bonds << Bond("ISING", "Jz", {9, 7}); - bonds << Bond("ISING", "Jz", {3, 11}); - bonds << Bond("ISING", "Jz", {7, 0}); - bonds << Bond("ISING", "Jz", {11, 5}); - bonds << Bond("ISING", "Jz", {0, 4}); - bonds << Bond("ISING", "Jz", {8, 3}); - bonds << Bond("ISING", "Jz", {2, 6}); - bonds << Bond("ISING", "Jz", {1, 5}); - bonds << Bond("ISING", "Jz", {4, 9}); - bonds << Bond("ISING", "Jz", {6, 11}); - bonds << Bond("ISING", "Jz", {10, 1}); - bonds << Bond("ISING", "Jz", {5, 8}); - bonds << Bond("ISING", "Jz", {9, 2}); - bonds << Bond("ISING", "Jz", {3, 7}); - bonds << Bond("ISING", "Jz", {7, 10}); - bonds << Bond("ISING", "Jz", {11, 0}); - bonds << Bond("EXCHANGE", "Jx", {0, 8}); - bonds << Bond("EXCHANGE", "Jx", {8, 6}); - bonds << Bond("EXCHANGE", "Jx", {2, 10}); - bonds << Bond("EXCHANGE", "Jx", {1, 9}); - bonds << Bond("EXCHANGE", "Jx", {4, 3}); - bonds << Bond("EXCHANGE", "Jx", {6, 1}); - bonds << Bond("EXCHANGE", "Jx", {10, 4}); - bonds << Bond("EXCHANGE", "Jx", {5, 2}); - bonds << Bond("EXCHANGE", "Jx", {9, 7}); - bonds << Bond("EXCHANGE", "Jx", {3, 11}); - bonds << Bond("EXCHANGE", "Jx", {7, 0}); - bonds << Bond("EXCHANGE", "Jx", {11, 5}); - bonds << Bond("EXCHANGE", "Jx", {0, 10}); - bonds << Bond("EXCHANGE", "Jx", {8, 4}); - bonds << Bond("EXCHANGE", "Jx", {2, 8}); - bonds << Bond("EXCHANGE", "Jx", {1, 11}); - bonds << Bond("EXCHANGE", "Jx", {4, 1}); - bonds << Bond("EXCHANGE", "Jx", {6, 3}); - bonds << Bond("EXCHANGE", "Jx", {10, 6}); - bonds << Bond("EXCHANGE", "Jx", {5, 0}); - bonds << Bond("EXCHANGE", "Jx", {9, 5}); - bonds << Bond("EXCHANGE", "Jx", {3, 9}); - bonds << Bond("EXCHANGE", "Jx", {7, 2}); - bonds << Bond("EXCHANGE", "Jx", {11, 7}); - bonds << Bond("EXCHANGE", "Jx", {0, 11}); - bonds << Bond("EXCHANGE", "Jx", {8, 5}); - bonds << Bond("EXCHANGE", "Jx", {2, 9}); - bonds << Bond("EXCHANGE", "Jx", {1, 10}); - bonds << Bond("EXCHANGE", "Jx", {4, 0}); - bonds << Bond("EXCHANGE", "Jx", {6, 2}); - bonds << Bond("EXCHANGE", "Jx", {10, 7}); - bonds << Bond("EXCHANGE", "Jx", {5, 1}); - bonds << Bond("EXCHANGE", "Jx", {9, 4}); - bonds << Bond("EXCHANGE", "Jx", {3, 8}); - bonds << Bond("EXCHANGE", "Jx", {7, 3}); - bonds << Bond("EXCHANGE", "Jx", {11, 6}); + bonds += Bond("ISING", "Jz", {0, 5}); + bonds += Bond("ISING", "Jz", {8, 2}); + bonds += Bond("ISING", "Jz", {2, 7}); + bonds += Bond("ISING", "Jz", {1, 4}); + bonds += Bond("ISING", "Jz", {4, 8}); + bonds += Bond("ISING", "Jz", {6, 10}); + bonds += Bond("ISING", "Jz", {10, 0}); + bonds += Bond("ISING", "Jz", {5, 9}); + bonds += Bond("ISING", "Jz", {9, 3}); + bonds += Bond("ISING", "Jz", {3, 6}); + bonds += Bond("ISING", "Jz", {7, 11}); + bonds += Bond("ISING", "Jz", {11, 1}); + bonds += Bond("ISING", "Jz", {0, 8}); + bonds += Bond("ISING", "Jz", {8, 6}); + bonds += Bond("ISING", "Jz", {2, 10}); + bonds += Bond("ISING", "Jz", {1, 9}); + bonds += Bond("ISING", "Jz", {4, 3}); + bonds += Bond("ISING", "Jz", {6, 1}); + bonds += Bond("ISING", "Jz", {10, 4}); + bonds += Bond("ISING", "Jz", {5, 2}); + bonds += Bond("ISING", "Jz", {9, 7}); + bonds += Bond("ISING", "Jz", {3, 11}); + bonds += Bond("ISING", "Jz", {7, 0}); + bonds += Bond("ISING", "Jz", {11, 5}); + bonds += Bond("ISING", "Jz", {0, 4}); + bonds += Bond("ISING", "Jz", {8, 3}); + bonds += Bond("ISING", "Jz", {2, 6}); + bonds += Bond("ISING", "Jz", {1, 5}); + bonds += Bond("ISING", "Jz", {4, 9}); + bonds += Bond("ISING", "Jz", {6, 11}); + bonds += Bond("ISING", "Jz", {10, 1}); + bonds += Bond("ISING", "Jz", {5, 8}); + bonds += Bond("ISING", "Jz", {9, 2}); + bonds += Bond("ISING", "Jz", {3, 7}); + bonds += Bond("ISING", "Jz", {7, 10}); + bonds += Bond("ISING", "Jz", {11, 0}); + bonds += Bond("EXCHANGE", "Jx", {0, 8}); + bonds += Bond("EXCHANGE", "Jx", {8, 6}); + bonds += Bond("EXCHANGE", "Jx", {2, 10}); + bonds += Bond("EXCHANGE", "Jx", {1, 9}); + bonds += Bond("EXCHANGE", "Jx", {4, 3}); + bonds += Bond("EXCHANGE", "Jx", {6, 1}); + bonds += Bond("EXCHANGE", "Jx", {10, 4}); + bonds += Bond("EXCHANGE", "Jx", {5, 2}); + bonds += Bond("EXCHANGE", "Jx", {9, 7}); + bonds += Bond("EXCHANGE", "Jx", {3, 11}); + bonds += Bond("EXCHANGE", "Jx", {7, 0}); + bonds += Bond("EXCHANGE", "Jx", {11, 5}); + bonds += Bond("EXCHANGE", "Jx", {0, 10}); + bonds += Bond("EXCHANGE", "Jx", {8, 4}); + bonds += Bond("EXCHANGE", "Jx", {2, 8}); + bonds += Bond("EXCHANGE", "Jx", {1, 11}); + bonds += Bond("EXCHANGE", "Jx", {4, 1}); + bonds += Bond("EXCHANGE", "Jx", {6, 3}); + bonds += Bond("EXCHANGE", "Jx", {10, 6}); + bonds += Bond("EXCHANGE", "Jx", {5, 0}); + bonds += Bond("EXCHANGE", "Jx", {9, 5}); + bonds += Bond("EXCHANGE", "Jx", {3, 9}); + bonds += Bond("EXCHANGE", "Jx", {7, 2}); + bonds += Bond("EXCHANGE", "Jx", {11, 7}); + bonds += Bond("EXCHANGE", "Jx", {0, 11}); + bonds += Bond("EXCHANGE", "Jx", {8, 5}); + bonds += Bond("EXCHANGE", "Jx", {2, 9}); + bonds += Bond("EXCHANGE", "Jx", {1, 10}); + bonds += Bond("EXCHANGE", "Jx", {4, 0}); + bonds += Bond("EXCHANGE", "Jx", {6, 2}); + bonds += Bond("EXCHANGE", "Jx", {10, 7}); + bonds += Bond("EXCHANGE", "Jx", {5, 1}); + bonds += Bond("EXCHANGE", "Jx", {9, 4}); + bonds += Bond("EXCHANGE", "Jx", {3, 8}); + bonds += Bond("EXCHANGE", "Jx", {7, 3}); + bonds += Bond("EXCHANGE", "Jx", {11, 6}); bonds["Jz"] = 1.0; bonds["Jx"] = complex(cos(2 * M_PI * eta), sin(2 * M_PI * eta)); diff --git a/tests/blocks/spinhalf_symmetric/test_kitaev_gamma.cpp b/tests/blocks/spinhalf_symmetric/test_kitaev_gamma.cpp index 6f72932e..3177ad5b 100644 --- a/tests/blocks/spinhalf_symmetric/test_kitaev_gamma.cpp +++ b/tests/blocks/spinhalf_symmetric/test_kitaev_gamma.cpp @@ -1,9 +1,9 @@ #include "../../catch.hpp" -#include #include #include #include +#include #include void run_kitaev_gamma_test( @@ -22,23 +22,21 @@ void run_kitaev_gamma_test( cx_mat sy(mat({{0., 0.}, {0., 0.}}), mat({{0., -0.5}, {0.5, 0.}})); cx_mat sz(mat({{0.5, 0.0}, {0.0, -0.5}}), mat({{0., 0.}, {0., 0.0}})); - cx_mat kx = kron(sx, sx); - cx_mat ky = kron(sy, sy); - cx_mat kz = kron(sz, sz); - cx_mat gx = kron(sy, sz) + kron(sz, sy); - cx_mat gy = kron(sx, sz) + kron(sz, sx); - cx_mat gz = kron(sx, sy) + kron(sy, sx); - - bonds.set_matrix("GENERICKITAEVX", kx); - bonds.set_matrix("GENERICKITAEVY", ky); - bonds.set_matrix("GENERICKITAEVZ", kz); - bonds.set_matrix("GENERICGAMMAX", gx); - bonds.set_matrix("GENERICGAMMAY", gy); - bonds.set_matrix("GENERICGAMMAZ", gz); + cx_mat kx = K * kron(sx, sx); + cx_mat ky = K * kron(sy, sy); + cx_mat kz = K * kron(sz, sz); + cx_mat gx = G * (kron(sy, sz) + kron(sz, sy)); + cx_mat gy = G * (kron(sx, sz) + kron(sz, sx)); + cx_mat gz = G * (kron(sx, sy) + kron(sy, sx)); bonds["J"] = 0.; - bonds["K"] = K; - bonds["G"] = G; + bonds["KX"] = kx; + bonds["KY"] = ky; + bonds["KZ"] = kz; + bonds["GX"] = gx; + bonds["GY"] = gy; + bonds["GZ"] = gz; + for (auto [name, e0_reference] : irrep_names_e0) { auto irrep = read_representation(lfile, name); auto block = Spinhalf(8, group, irrep); diff --git a/tests/blocks/tj/test_tj_matrix.cpp b/tests/blocks/tj/test_tj_matrix.cpp index 1c953346..af7499da 100644 --- a/tests/blocks/tj/test_tj_matrix.cpp +++ b/tests/blocks/tj/test_tj_matrix.cpp @@ -8,11 +8,12 @@ #include #include #include +#include using namespace xdiag; -void test_tjmodel_e0_real(BondList bonds, int nup, int ndn, double e0) { - int n_sites = bonds.n_sites(); +void test_tjmodel_e0_real(BondList bonds, int n_sites, int nup, int ndn, + double e0) { auto block = tJ(n_sites, nup, ndn); auto H = matrix(bonds, block, block); arma::vec eigs; @@ -22,8 +23,8 @@ void test_tjmodel_e0_real(BondList bonds, int nup, int ndn, double e0) { REQUIRE(std::abs(e0 - eigs(0)) < 1e-6); } -void test_tjmodel_fulleigs(BondList bonds, arma::Col exact_eigs) { - int n_sites = bonds.n_sites(); +void test_tjmodel_fulleigs(BondList bonds, int n_sites, + arma::Col exact_eigs) { std::vector all_eigs; for (int ndn = 0; ndn <= n_sites; ++ndn) { @@ -48,7 +49,7 @@ void test_tjmodel_fulleigs(BondList bonds, arma::Col exact_eigs) { REQUIRE(close(arma::vec(all_eigs), exact_eigs)); } -TEST_CASE("tj_matrix", "[tj]") { +TEST_CASE("tj_matrix", "[tj]") try { using namespace xdiag::testcases::tj; { @@ -96,7 +97,7 @@ TEST_CASE("tj_matrix", "[tj]") { {4, 2, -2.11803398}, {5, 0, -0.99999999}, {5, 1, -0.49999999}, {6, 0, 1.500000000}}; for (auto [nup, ndn, e0] : nup_ndn_e0) - test_tjmodel_e0_real(bonds, nup, ndn, e0); + test_tjmodel_e0_real(bonds, 6, nup, ndn, e0); } { @@ -114,26 +115,33 @@ TEST_CASE("tj_matrix", "[tj]") { {4, 2, 0.000000000}, {5, 0, -2.00000000}, {5, 1, 0.000000000}, {6, 0, 0.000000000}}; for (auto [nup, ndn, e0] : nup_ndn_e0) - test_tjmodel_e0_real(bonds, nup, ndn, e0); + test_tjmodel_e0_real(bonds, 6, nup, ndn, e0); } for (int L = 3; L <= 6; ++L) { Log.out("tj_matrix: ALPS full spectrum test, chain N={}", L); auto [bonds, eigs] = tJchain_fullspectrum_alps(L); - test_tjmodel_fulleigs(bonds, eigs); + test_tjmodel_fulleigs(bonds, L, eigs); } { Log.out("tj_matrix: ALPS full spectrum test, square 2x2"); auto [bonds, eigs] = tj_square2x2_fullspectrum_alps(); - test_tjmodel_fulleigs(bonds, eigs); + test_tjmodel_fulleigs(bonds, 4, eigs); } for (int N = 3; N <= 6; ++N) { Log.out("tj_matrix: random all-to-all complex exchange test, N={}", N); auto bonds = tj_alltoall_complex(N); - auto bonds_hb = bonds.bonds_of_type("HB"); + BondList bonds_hb; + for (auto bond : bonds) { + if (bond.type() == "HB") { + bonds_hb += bond; + std::string name = bond.coupling().as(); + bonds_hb[name] = bonds[name]; + } + } for (int nup = 0; nup <= N; ++nup) for (int ndn = 0; ndn <= N - nup; ++ndn) { @@ -147,6 +155,7 @@ TEST_CASE("tj_matrix", "[tj]") { int ndn = N - nup; auto block1 = tJ(N, nup, ndn); auto block2 = Spinhalf(N, nup); + auto H1 = matrixC(bonds_hb, block1, block1); auto H2 = matrixC(bonds_hb, block2, block2); arma::vec eigs1; @@ -161,12 +170,14 @@ TEST_CASE("tj_matrix", "[tj]") { { Log.out("tj_matrix: Henry's Matlab test, random 3"); auto [bonds, eigs] = randomAlltoAll3(); - test_tjmodel_fulleigs(bonds, eigs); + test_tjmodel_fulleigs(bonds, 3, eigs); } { Log.out("tj_matrix: Henry's Matlab test, random 4"); auto [bonds, eigs] = randomAlltoAll4(); - test_tjmodel_fulleigs(bonds, eigs); + test_tjmodel_fulleigs(bonds, 4, eigs); } +} catch (Error e) { + xdiag::error_trace(e); } diff --git a/tests/blocks/tj/test_tj_raiselower.cpp b/tests/blocks/tj/test_tj_raiselower.cpp index 07f2d7c1..f223478f 100644 --- a/tests/blocks/tj/test_tj_raiselower.cpp +++ b/tests/blocks/tj/test_tj_raiselower.cpp @@ -112,8 +112,8 @@ TEST_CASE("tj_raise_lower", "[tj]") { auto block_j = tJ(n_sites, nup_j, ndn_j); auto block_ij = tJ(n_sites, nup_ij, ndn_ij); - auto op_i = Bond(op_i_str, i); - auto op_j = Bond(op_j_str, j); + auto op_i = Bond(op_i_str, 1.0, i); + auto op_j = Bond(op_j_str, 1.0, j); auto op_i_m = matrixC(op_i, block, block_i); auto op_ij_m = matrixC(op_j, block_i, block_ij); @@ -141,12 +141,12 @@ TEST_CASE("tj_raise_lower", "[tj]") { REQUIRE(norm(anti_comm) < 1e-12); } else if (((op_i_str == "CDAGUP") && (op_j_str == "CUP")) || ((op_i_str == "CUP") && (op_j_str == "CDAGUP"))) { - auto ndn_op = Bond("NUMBERDN", i); + auto ndn_op = Bond("NUMBERDN", 1.0, i); auto ndn_op_m = matrix(ndn_op, block); REQUIRE(norm(anti_comm + ndn_op_m - id) < 1e-12); } else if (((op_i_str == "CDAGDN") && (op_j_str == "CDN")) || ((op_i_str == "CDN") && (op_j_str == "CDAGDN"))) { - auto nup_op = Bond("NUMBERUP", i); + auto nup_op = Bond("NUMBERUP", 1.0, i); auto nup_op_m = matrix(nup_op, block); REQUIRE(norm(anti_comm + nup_op_m - id) < 1e-12); } else if ((op_i_str == "CDAGUP") && (op_j_str == "CDN")) { diff --git a/tests/blocks/tj/testcases_tj.cpp b/tests/blocks/tj/testcases_tj.cpp index 4e2e2997..75b3e1b0 100644 --- a/tests/blocks/tj/testcases_tj.cpp +++ b/tests/blocks/tj/testcases_tj.cpp @@ -1,4 +1,5 @@ #include "testcases_tj.hpp" +#include namespace xdiag::testcases::tj { @@ -8,8 +9,8 @@ BondList tJchain(int n_sites, double t, double J) { bondlist["T"] = t; bondlist["J"] = J; for (int s = 0; s < n_sites; ++s) { - bondlist << Bond("HOP", "T", {s, (s + 1) % n_sites}); - bondlist << Bond("HB", "J", {s, (s + 1) % n_sites}); + bondlist += Bond("HOP", "T", {s, (s + 1) % n_sites}); + bondlist += Bond("HB", "J", {s, (s + 1) % n_sites}); } return bondlist; } @@ -20,8 +21,8 @@ std::tuple> tJchain_fullspectrum_alps(int L) { bonds["T"] = 1.0; bonds["J"] = 1.0; for (int s = 0; s < L; ++s) { - bonds << Bond("HOP", "T", {s, (s + 1) % L}); - bonds << Bond("TJHB", "J", {s, (s + 1) % L}); + bonds += Bond("HOP", "T", {s, (s + 1) % L}); + bonds += Bond("TJHB", "J", {s, (s + 1) % L}); } arma::Col eigs; if (L == 3) { @@ -581,22 +582,22 @@ tj_square2x2_fullspectrum_alps() { BondList bondlist; bondlist["T"] = 1.0; bondlist["J"] = 1.0; - bondlist << Bond("HOP", "T", {0, 1}); - bondlist << Bond("HOP", "T", {1, 0}); - bondlist << Bond("HOP", "T", {2, 3}); - bondlist << Bond("HOP", "T", {3, 2}); - bondlist << Bond("HOP", "T", {0, 2}); - bondlist << Bond("HOP", "T", {2, 0}); - bondlist << Bond("HOP", "T", {1, 3}); - bondlist << Bond("HOP", "T", {3, 1}); - bondlist << Bond("TJHB", "J", {0, 1}); - bondlist << Bond("TJHB", "J", {1, 0}); - bondlist << Bond("TJHB", "J", {2, 3}); - bondlist << Bond("TJHB", "J", {3, 2}); - bondlist << Bond("TJHB", "J", {0, 2}); - bondlist << Bond("TJHB", "J", {2, 0}); - bondlist << Bond("TJHB", "J", {1, 3}); - bondlist << Bond("TJHB", "J", {3, 1}); + bondlist += Bond("HOP", "T", {0, 1}); + bondlist += Bond("HOP", "T", {1, 0}); + bondlist += Bond("HOP", "T", {2, 3}); + bondlist += Bond("HOP", "T", {3, 2}); + bondlist += Bond("HOP", "T", {0, 2}); + bondlist += Bond("HOP", "T", {2, 0}); + bondlist += Bond("HOP", "T", {1, 3}); + bondlist += Bond("HOP", "T", {3, 1}); + bondlist += Bond("TJHB", "J", {0, 1}); + bondlist += Bond("TJHB", "J", {1, 0}); + bondlist += Bond("TJHB", "J", {2, 3}); + bondlist += Bond("TJHB", "J", {3, 2}); + bondlist += Bond("TJHB", "J", {0, 2}); + bondlist += Bond("TJHB", "J", {2, 0}); + bondlist += Bond("TJHB", "J", {1, 3}); + bondlist += Bond("TJHB", "J", {3, 1}); auto eigs = arma::Col({-6.744562646538028616e+00, -6.000000000000000000e+00, -5.605551275463989569e+00, -5.605551275463988681e+00, @@ -653,7 +654,7 @@ BondList tj_alltoall(int n_sites) { ss << "T" << s1 << "_" << s2; std::string name = ss.str(); double value = distribution(generator); - bondlist << Bond("HOP", name, {s1, s2}); + bondlist += Bond("HOP", name, {s1, s2}); bondlist[name] = value; } for (int s1 = 0; s1 < n_sites; ++s1) @@ -662,7 +663,7 @@ BondList tj_alltoall(int n_sites) { ss << "J" << s1 << "_" << s2; std::string name = ss.str(); double value = distribution(generator); - bondlist << Bond("HB", name, {s1, s2}); + bondlist += Bond("HB", name, {s1, s2}); bondlist[name] = value; } return bondlist; @@ -679,7 +680,7 @@ BondList tj_alltoall_complex(int n_sites) { ss << "T" << s1 << "_" << s2; std::string name = ss.str(); complex value = complex(distribution(generator), distribution(generator)); - bondlist << Bond("HOP", name, {s1, s2}); + bondlist += Bond("HOP", name, {s1, s2}); bondlist[name] = value; } for (int s1 = 0; s1 < n_sites; ++s1) @@ -688,7 +689,7 @@ BondList tj_alltoall_complex(int n_sites) { ss << "J" << s1 << "_" << s2; std::string name = ss.str(); double value = distribution(generator); - bondlist << Bond("HB", name, {s1, s2}); + bondlist += Bond("HB", name, {s1, s2}); bondlist[name] = value; } return bondlist; @@ -708,18 +709,18 @@ std::tuple> randomAlltoAll4() { bondlist["J13"] = 1; bondlist["T23"] = -4; bondlist["J23"] = 4; - bondlist << Bond("HOP", "T01", {0, 1}); - bondlist << Bond("HOP", "T02", {0, 2}); - bondlist << Bond("HOP", "T03", {0, 3}); - bondlist << Bond("HOP", "T12", {1, 2}); - bondlist << Bond("HOP", "T13", {1, 3}); - bondlist << Bond("HOP", "T23", {2, 3}); - bondlist << Bond("HB", "J01", {0, 1}); - bondlist << Bond("HB", "J02", {0, 2}); - bondlist << Bond("HB", "J03", {0, 3}); - bondlist << Bond("HB", "J12", {1, 2}); - bondlist << Bond("HB", "J13", {1, 3}); - bondlist << Bond("HB", "J23", {2, 3}); + bondlist += Bond("HOP", "T01", {0, 1}); + bondlist += Bond("HOP", "T02", {0, 2}); + bondlist += Bond("HOP", "T03", {0, 3}); + bondlist += Bond("HOP", "T12", {1, 2}); + bondlist += Bond("HOP", "T13", {1, 3}); + bondlist += Bond("HOP", "T23", {2, 3}); + bondlist += Bond("HB", "J01", {0, 1}); + bondlist += Bond("HB", "J02", {0, 2}); + bondlist += Bond("HB", "J03", {0, 3}); + bondlist += Bond("HB", "J12", {1, 2}); + bondlist += Bond("HB", "J13", {1, 3}); + bondlist += Bond("HB", "J23", {2, 3}); arma::Col eigs = { 11.248037068163532, 11.248037068163532, 11.248037068163525, @@ -762,12 +763,12 @@ std::tuple> randomAlltoAll3() { bondlist["J02"] = -1; bondlist["T12"] = -5; bondlist["J12"] = -3; - bondlist << Bond("HOP", "T01", {0, 1}); - bondlist << Bond("HOP", "T02", {0, 2}); - bondlist << Bond("HOP", "T12", {1, 2}); - bondlist << Bond("HB", "J01", {0, 1}); - bondlist << Bond("HB", "J02", {0, 2}); - bondlist << Bond("HB", "J12", {1, 2}); + bondlist += Bond("HOP", "T01", {0, 1}); + bondlist += Bond("HOP", "T02", {0, 2}); + bondlist += Bond("HOP", "T12", {1, 2}); + bondlist += Bond("HB", "J01", {0, 1}); + bondlist += Bond("HB", "J02", {0, 2}); + bondlist += Bond("HB", "J12", {1, 2}); arma::Col eigs = { 6.256066270684332, 5.099019513592784, 5.099019513592784, diff --git a/tests/blocks/tj_symmetric/test_tj_symmetric_matrix.cpp b/tests/blocks/tj_symmetric/test_tj_symmetric_matrix.cpp index 2f84747f..0cf9b4aa 100644 --- a/tests/blocks/tj_symmetric/test_tj_symmetric_matrix.cpp +++ b/tests/blocks/tj_symmetric/test_tj_symmetric_matrix.cpp @@ -100,7 +100,7 @@ TEST_CASE("tj_symmetric_matrix", "[tj]") { n_sites); BondList bonds; for (int64_t s = 0; s < n_sites; ++s) { - bonds << Bond("TJHB", {s, (s + 1) % n_sites}); + bonds += Bond("TJHB", 1.0, {s, (s + 1) % n_sites}); } auto [space_group, irreps, multiplicities] = get_cyclic_group_irreps_mult(n_sites); diff --git a/tests/io/test_file_toml.cpp b/tests/io/test_file_toml.cpp index 083a2778..773e9b88 100644 --- a/tests/io/test_file_toml.cpp +++ b/tests/io/test_file_toml.cpp @@ -1,20 +1,21 @@ #include "../catch.hpp" #include -#include #include #include #include +#include #include #include #include -#include +#include +#include -template void test_write_read(T val) { +template void test_write_read(T val) try { using namespace xdiag; // Test writing - std::string filename = XDIAG_DIRECTORY"/misc/data/toml/write.toml"; + std::string filename = XDIAG_DIRECTORY "/misc/data/toml/write.toml"; auto fl = FileToml(filename, 'w'); std::string key = "val"; @@ -27,14 +28,16 @@ template void test_write_read(T val) { // XDIAG_SHOW(val); // XDIAG_SHOW(val_r); REQUIRE(val == val_r); + } catch (xdiag::Error const &e) { + XDIAG_RETHROW(e); } -TEST_CASE("file_toml", "[io]") { +TEST_CASE("file_toml", "[io]") try { using namespace xdiag; using namespace arma; // Just try to parse everything in the example toml file - std::string filename = XDIAG_DIRECTORY"/misc/data/toml/read.toml"; + std::string filename = XDIAG_DIRECTORY "/misc/data/toml/read.toml"; auto fl = FileToml(filename, 'r'); // Parse String @@ -157,7 +160,7 @@ TEST_CASE("file_toml", "[io]") { test_write_read(f); test_write_read(i); - filename = XDIAG_DIRECTORY"/misc/data/toml/write.toml"; + filename = XDIAG_DIRECTORY "/misc/data/toml/write.toml"; fl = FileToml(filename, 'w'); fl["a"] = a; fl["b"] = b; @@ -206,7 +209,8 @@ TEST_CASE("file_toml", "[io]") { fl["perm"] = p; test_write_read(p); - std::string lfile = XDIAG_DIRECTORY"/misc/data/triangular.j1j2jch/" + std::string lfile = + XDIAG_DIRECTORY "/misc/data/triangular.j1j2jch/" "triangular.12.j1j2jch.sublattices.fsl.lat"; auto group = PermutationGroup(read_permutations(lfile)); test_write_read(group); @@ -234,42 +238,18 @@ TEST_CASE("file_toml", "[io]") { auto matr = arma::mat(3, 3, arma::fill::randu); - bond = Bond(matr, "H", 1); - test_write_read(bond); - - bond = Bond(matr, "J1", {1, 2}); + bond = Bond("H", matr, 1); test_write_read(bond); - bond = Bond(matr, 1.0, 1); - test_write_read(bond); - - bond = Bond(matr, 1.2, {1, 2}); - test_write_read(bond); - - bond = Bond(matr, 2.1 + 1.2i, 1); - test_write_read(bond); - - bond = Bond(matr, 2.1 + 1.2i, {1, 2}); + bond = Bond("J1", matr, {1, 2}); test_write_read(bond); auto matc = arma::cx_mat(3, 3, arma::fill::randu); - bond = Bond(matc, "H", 1); - test_write_read(bond); - - bond = Bond(matc, "J1", {1, 2}); - test_write_read(bond); - - bond = Bond(matc, 1.0, 1); + bond = Bond("H", matc, 1); test_write_read(bond); - bond = Bond(matc, 1.2, {1, 2}); - test_write_read(bond); - - bond = Bond(matc, 2.1 + 1.2i, 1); - test_write_read(bond); - - bond = Bond(matc, 2.1 + 1.2i, {1, 2}); + bond = Bond("J1", matc, {1, 2}); test_write_read(bond); auto bonds = read_bondlist(lfile); @@ -282,5 +262,7 @@ TEST_CASE("file_toml", "[io]") { bonds["M1"] = arma::mat(3, 4, arma::fill::randu); bonds["M2"] = arma::cx_mat(3, 4, arma::fill::randu); test_write_read(bonds); - + + } catch (xdiag::Error const &e) { + XDIAG_RETHROW(e); } diff --git a/tests/operators/test_non_branching_bonds.cpp b/tests/operators/test_non_branching_bonds.cpp index a8c28a23..f028a007 100644 --- a/tests/operators/test_non_branching_bonds.cpp +++ b/tests/operators/test_non_branching_bonds.cpp @@ -22,7 +22,7 @@ TEST_CASE("non_branching_bonds", "[operators]") try { cx_mat ones(mat({{1.0, 1.0}, {1.0, 1.0}}), mat({{1.0, 1.0}, {1.0, 1.0}})); for (auto ss : {ones}) { - auto bond = Bond(ss, 0); + auto bond = Bond("MAT", ss, 0); auto block = Spinhalf(1); auto h = matrixC(bond, block); REQUIRE(close(h, ss)); @@ -35,10 +35,10 @@ TEST_CASE("non_branching_bonds", "[operators]") try { BondList bonds; for (int64_t i = 0; i < N - 1; ++i) { - bonds << Bond("ISING", J, {i, (i + 1) % N}); + bonds += Bond("ISING", J, {i, (i + 1) % N}); } for (int64_t i = 0; i < N; ++i) { - bonds << Bond(sx, H, i); + bonds += Bond("SX", sx, i); } auto block = Spinhalf(N); @@ -57,12 +57,12 @@ TEST_CASE("non_branching_bonds", "[operators]") try { std::vector sites(k); std::iota(sites.begin(), sites.end(), 0); - auto bondr = Bond(mr, sites); + auto bondr = Bond("MR", mr, sites); auto hr = matrix(bondr, block); REQUIRE(close(hr, mr)); auto mc = cx_mat(p2, p2, fill::randn); - auto bondc = Bond(mc, sites); + auto bondc = Bond("MC", mc, sites); auto hc = matrixC(bondc, block); REQUIRE(close(hc, mc)); } @@ -73,10 +73,10 @@ TEST_CASE("non_branching_bonds", "[operators]") try { auto block6 = Spinhalf(6); BondList bonds1; - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {0, 1, 2}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {1, 2, 3}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {2, 3, 4}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {3, 4, 5}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {0, 1, 2}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {1, 2, 3}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {2, 3, 4}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {3, 4, 5}); bonds1["Jchi"] = 1.0; auto H1 = matrixC(bonds1, block6); @@ -85,11 +85,10 @@ TEST_CASE("non_branching_bonds", "[operators]") try { kron(sz, kron(sx, sy) - kron(sy, sx)); BondList bonds2; - bonds2 << Bond(jchimat, "Jchi", {0, 1, 2}); - bonds2 << Bond(jchimat, "Jchi", {1, 2, 3}); - bonds2 << Bond(jchimat, "Jchi", {2, 3, 4}); - bonds2 << Bond(jchimat, "Jchi", {3, 4, 5}); - bonds2["Jchi"] = 1.0; + bonds2 += Bond("Jchi", jchimat, {0, 1, 2}); + bonds2 += Bond("Jchi", jchimat, {1, 2, 3}); + bonds2 += Bond("Jchi", jchimat, {2, 3, 4}); + bonds2 += Bond("Jchi", jchimat, {3, 4, 5}); auto H2 = matrixC(bonds2, block6); REQUIRE(close(H1, H2)); @@ -99,30 +98,30 @@ TEST_CASE("non_branching_bonds", "[operators]") try { auto block12 = Spinhalf(12); BondList bonds1; - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {0, 4, 6}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {3, 1, 9}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {9, 7, 4}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {4, 2, 10}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {10, 8, 5}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {6, 10, 1}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {1, 5, 7}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {7, 11, 2}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {2, 3, 8}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {8, 9, 0}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {5, 0, 11}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {11, 6, 3}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {4, 10, 6}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {1, 7, 9}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {7, 2, 4}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {2, 8, 10}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {8, 0, 5}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {10, 5, 1}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {5, 11, 7}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {11, 3, 2}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {3, 9, 8}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {9, 4, 0}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {0, 6, 11}); - bonds1 << Bond("SCALARCHIRALITY", "Jchi", {6, 1, 3}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {0, 4, 6}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {3, 1, 9}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {9, 7, 4}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {4, 2, 10}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {10, 8, 5}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {6, 10, 1}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {1, 5, 7}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {7, 11, 2}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {2, 3, 8}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {8, 9, 0}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {5, 0, 11}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {11, 6, 3}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {4, 10, 6}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {1, 7, 9}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {7, 2, 4}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {2, 8, 10}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {8, 0, 5}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {10, 5, 1}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {5, 11, 7}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {11, 3, 2}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {3, 9, 8}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {9, 4, 0}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {0, 6, 11}); + bonds1 += Bond("SCALARCHIRALITY", "Jchi", {6, 1, 3}); bonds1["Jchi"] = 1.0; auto H1 = matrixC(bonds1, block12); @@ -131,31 +130,30 @@ TEST_CASE("non_branching_bonds", "[operators]") try { kron(kron(sy, sx), sz) - kron(kron(sz, sy), sx); BondList bonds2; - bonds2 << Bond(jchimat, "Jchi", {0, 4, 6}); - bonds2 << Bond(jchimat, "Jchi", {3, 1, 9}); - bonds2 << Bond(jchimat, "Jchi", {9, 7, 4}); - bonds2 << Bond(jchimat, "Jchi", {4, 2, 10}); - bonds2 << Bond(jchimat, "Jchi", {10, 8, 5}); - bonds2 << Bond(jchimat, "Jchi", {6, 10, 1}); - bonds2 << Bond(jchimat, "Jchi", {1, 5, 7}); - bonds2 << Bond(jchimat, "Jchi", {7, 11, 2}); - bonds2 << Bond(jchimat, "Jchi", {2, 3, 8}); - bonds2 << Bond(jchimat, "Jchi", {8, 9, 0}); - bonds2 << Bond(jchimat, "Jchi", {5, 0, 11}); - bonds2 << Bond(jchimat, "Jchi", {11, 6, 3}); - bonds2 << Bond(jchimat, "Jchi", {4, 10, 6}); - bonds2 << Bond(jchimat, "Jchi", {1, 7, 9}); - bonds2 << Bond(jchimat, "Jchi", {7, 2, 4}); - bonds2 << Bond(jchimat, "Jchi", {2, 8, 10}); - bonds2 << Bond(jchimat, "Jchi", {8, 0, 5}); - bonds2 << Bond(jchimat, "Jchi", {10, 5, 1}); - bonds2 << Bond(jchimat, "Jchi", {5, 11, 7}); - bonds2 << Bond(jchimat, "Jchi", {11, 3, 2}); - bonds2 << Bond(jchimat, "Jchi", {3, 9, 8}); - bonds2 << Bond(jchimat, "Jchi", {9, 4, 0}); - bonds2 << Bond(jchimat, "Jchi", {0, 6, 11}); - bonds2 << Bond(jchimat, "Jchi", {6, 1, 3}); - bonds2["Jchi"] = 1.0; + bonds2 += Bond("Jchi", jchimat, {0, 4, 6}); + bonds2 += Bond("Jchi", jchimat, {3, 1, 9}); + bonds2 += Bond("Jchi", jchimat, {9, 7, 4}); + bonds2 += Bond("Jchi", jchimat, {4, 2, 10}); + bonds2 += Bond("Jchi", jchimat, {10, 8, 5}); + bonds2 += Bond("Jchi", jchimat, {6, 10, 1}); + bonds2 += Bond("Jchi", jchimat, {1, 5, 7}); + bonds2 += Bond("Jchi", jchimat, {7, 11, 2}); + bonds2 += Bond("Jchi", jchimat, {2, 3, 8}); + bonds2 += Bond("Jchi", jchimat, {8, 9, 0}); + bonds2 += Bond("Jchi", jchimat, {5, 0, 11}); + bonds2 += Bond("Jchi", jchimat, {11, 6, 3}); + bonds2 += Bond("Jchi", jchimat, {4, 10, 6}); + bonds2 += Bond("Jchi", jchimat, {1, 7, 9}); + bonds2 += Bond("Jchi", jchimat, {7, 2, 4}); + bonds2 += Bond("Jchi", jchimat, {2, 8, 10}); + bonds2 += Bond("Jchi", jchimat, {8, 0, 5}); + bonds2 += Bond("Jchi", jchimat, {10, 5, 1}); + bonds2 += Bond("Jchi", jchimat, {5, 11, 7}); + bonds2 += Bond("Jchi", jchimat, {11, 3, 2}); + bonds2 += Bond("Jchi", jchimat, {3, 9, 8}); + bonds2 += Bond("Jchi", jchimat, {9, 4, 0}); + bonds2 += Bond("Jchi", jchimat, {0, 6, 11}); + bonds2 += Bond("Jchi", jchimat, {6, 1, 3}); auto H2 = matrixC(bonds2, block12); // XDIAG_SHOW(norm(H1)); diff --git a/tests/operators/test_symmetrized_operator.cpp b/tests/operators/test_symmetrized_operator.cpp index 99e9b41d..4d18586e 100644 --- a/tests/operators/test_symmetrized_operator.cpp +++ b/tests/operators/test_symmetrized_operator.cpp @@ -22,9 +22,10 @@ TEST_CASE("symmetrized_operator", "[symmetries]") try { // int n_sites = 6; auto bondlist = testcases::electron::get_linear_chain(n_sites, 1.0, 5.0); for (int i = 0; i < n_sites; ++i) { - bondlist << Bond("HB", "J2", {i, (i + 2) % n_sites}); + bondlist += Bond("HB", "J2", {i, (i + 2) % n_sites}); } bondlist["J2"] = 0.321; + bondlist["T"] = 0; auto [space_group, irreps] = testcases::electron::get_cyclic_group_irreps(n_sites); @@ -35,6 +36,7 @@ TEST_CASE("symmetrized_operator", "[symmetries]") try { if (block_nosym.size() == 0) { continue; } + // XDIAG_SHOW(bondlist); auto [e0_nosym, v0_nosym] = eig0(bondlist, block_nosym); auto res = eigs_lanczos(bondlist, block_nosym); @@ -86,7 +88,9 @@ TEST_CASE("symmetrized_operator", "[symmetries]") try { { auto &v = v0_sym; auto Hv = v; + // Log("UUU"); apply(bondlist, v, Hv); + // Log("VVV"); auto e = dot(v, Hv); // XDIAG_SHOW(e); // XDIAG_SHOW(e0_nosym); @@ -96,19 +100,23 @@ TEST_CASE("symmetrized_operator", "[symmetries]") try { // Measure correlators for (int i = 1; i < n_sites; ++i) { BondList corr_nosym; - corr_nosym << Bond("HB", "J", {0, i}); + corr_nosym += Bond("HB", "J", {0, i}); corr_nosym["J"] = 1.0; auto corr_sym = symmetrized_operator(corr_nosym, space_group); + // Log("XXX"); auto val_nosym = inner(corr_nosym, v0_nosym); - auto val_sym = inner(corr_sym, v0_sym); + // Log("YYY {} {}", corr_sym.isreal(), v0_sym.isreal()); + // XDIAG_SHOW(corr_sym); + auto val_sym = innerC(corr_sym, v0_sym); + // Log("ZZZ"); // Log.out("(0,{}) nosym: {}, sym: {}", i, real(val_nosym), // real(val_sym)); REQUIRE(close(val_nosym, val_sym, 1e-6, 1e-6)); } // for (int j=0; j #include +#include #include using namespace xdiag; @@ -23,15 +23,15 @@ TEST_CASE("product_state", "[states]") { auto block = Electron(n_sites); auto psi = State(block, true); fill(psi, pstate); - + int nup = 0; int ndn = 0; for (int i = 0; i < n_sites; ++i) { std::string p = pstate[i]; - auto sz = Bond("SZ", i); + auto sz = Bond("SZ", 1.0, i); auto szc = inner(sz, psi); - auto n = Bond("NUMBER", i); + auto n = Bond("NUMBER", 1.0, i); auto nc = inner(n, psi); if (p == "Emp") { @@ -56,13 +56,13 @@ TEST_CASE("product_state", "[states]") { auto block2 = Electron(n_sites, nup, ndn); auto psi2 = State(block2); fill(psi2, pstate); - + for (int i = 0; i < n_sites; ++i) { std::string p = pstate[i]; - auto sz = Bond("SZ", i); + auto sz = Bond("SZ", 1.0, i); auto szc = inner(sz, psi2); - auto n = Bond("NUMBER", i); + auto n = Bond("NUMBER", 1.0, i); auto nc = inner(n, psi2); if (p == "Emp") { @@ -112,10 +112,10 @@ TEST_CASE("product_state", "[states]") { fill(psi2, pstate); for (int i = 0; i < n_sites; ++i) { std::string p = pstate[i]; - auto sz = Bond("SZ", i); + auto sz = Bond("SZ", 1.0, i); auto szc = inner(sz, psi2); - auto n = Bond("NUMBER", i); + auto n = Bond("NUMBER", 1.0, i); auto nc = inner(n, psi2); if (p == "Emp") { @@ -151,14 +151,14 @@ TEST_CASE("product_state", "[states]") { for (int i = 0; i < n_sites; ++i) { std::string p = pstate[i]; - auto sz = Bond("SZ", i); + auto sz = Bond("SZ", 1.0, i); auto szc = inner(sz, psi); - + if (p == "Up") { REQUIRE(close(szc, 0.5)); ++nup; } else if (p == "Dn") { - REQUIRE(close(szc, -0.5)); + REQUIRE(close(szc, -0.5)); } } @@ -167,7 +167,7 @@ TEST_CASE("product_state", "[states]") { fill(psi2, pstate); for (int i = 0; i < n_sites; ++i) { std::string p = pstate[i]; - auto sz = Bond("SZ", i); + auto sz = Bond("SZ", 1.0, i); auto szc = inner(sz, psi2); if (p == "Up") { diff --git a/xdiag/algebra/algebra.cpp b/xdiag/algebra/algebra.cpp index 6a48c121..d991e83e 100644 --- a/xdiag/algebra/algebra.cpp +++ b/xdiag/algebra/algebra.cpp @@ -34,7 +34,6 @@ double dot(State const &v, State const &w) try { if ((v.n_cols() > 1) || (w.n_cols() > 1)) { XDIAG_THROW( "Cannot compute dot product of state with more than one column"); - return 0; } if ((v.isreal()) && (w.isreal())) { @@ -42,11 +41,9 @@ double dot(State const &v, State const &w) try { } else { XDIAG_THROW("Unable to compute real dot product of a complex " "state, consider using dotC instead."); - return 0; } } catch (Error const &error) { XDIAG_RETHROW(error); - return 0; } complex dotC(State const &v, State const &w) try { @@ -60,7 +57,7 @@ complex dotC(State const &v, State const &w) try { } if ((v.isreal()) && (w.isreal())) { return dot(v.block(), v.vector(0, false), w.vector(0, false)); - } else if ((v.isreal()) && (w.iscomplex())) { + } else if ((v.isreal()) && (!w.isreal())) { State v2; try { v2 = v; @@ -69,7 +66,7 @@ complex dotC(State const &v, State const &w) try { XDIAG_THROW("Unable to create intermediate complex State"); } return dot(v.block(), v2.vectorC(0, false), w.vectorC(0, false)); - } else if ((w.isreal()) && (v.iscomplex())) { + } else if ((w.isreal()) && (!v.isreal())) { State w2; try { w2 = w; @@ -77,78 +74,87 @@ complex dotC(State const &v, State const &w) try { } catch (...) { XDIAG_THROW("Unable to create intermediate complex State"); } - return dot(v.block(), v.vectorC(0, false), w.vectorC(0, false)); + return dot(v.block(), v.vectorC(0, false), w2.vectorC(0, false)); } else { return dot(v.block(), v.vectorC(0, false), w.vectorC(0, false)); } } catch (Error const &error) { XDIAG_RETHROW(error); - return 0; } double inner(BondList const &bonds, State const &v) try { - auto w = v; - apply(bonds, v, w); - return dot(w, v); + if (v.isreal() && bonds.isreal()) { + auto w = v; + apply(bonds, v, w); + return dot(w, v); + } else { + XDIAG_THROW("\"inner\" function computing product can only " + "be called if both the state and the bonds are real. Maybe use " + "innerC(...)"); + } } catch (Error const &error) { XDIAG_RETHROW(error); - return 0; } -complex innerC(BondList const &bonds, State const &v) try { - auto w = v; - apply(bonds, v, w); - return dotC(w, v); +double inner(Bond const &bond, State const &v) try { + return inner(BondList({bond}), v); } catch (Error const &error) { XDIAG_RETHROW(error); - return 0; } -double inner(Bond const &bond, State const &v) try { - return inner(BondList({bond}), v); +complex innerC(BondList const &bonds, State const &v) try { + if (v.isreal() && bonds.isreal()) { + auto w = v; + apply(bonds, v, w); + return (complex)dot(w, v); + } else if (v.isreal() && !bonds.isreal()) { + auto v2 = v; + auto w = v2; + v2.make_complex(); + apply(bonds, v2, w); + return dotC(w, v); + } else { + auto w = v; + apply(bonds, v, w); + return dotC(w, v); + } } catch (Error const &error) { XDIAG_RETHROW(error); - return 0; } complex innerC(Bond const &bond, State const &v) try { return innerC(BondList({bond}), v); } catch (Error const &error) { XDIAG_RETHROW(error); - return 0; } -double inner(State const &v, BondList const &bonds, State const &w) try { - auto bw = zeros_like(w); - apply(bonds, w, bw); - return dot(v, bw); -} catch (Error const &error) { - XDIAG_RETHROW(error); - return 0; -} +// double inner(State const &v, BondList const &bonds, State const &w) try { +// auto bw = zeros_like(w); +// apply(bonds, w, bw); +// return dot(v, bw); +// } catch (Error const &error) { +// XDIAG_RETHROW(error); +// } -complex innerC(State const &v, BondList const &bonds, State const &w) try { - auto bw = zeros_like(w); - apply(bonds, w, bw); - return dotC(v, bw); -} catch (Error const &error) { - XDIAG_RETHROW(error); - return 0; -} +// complex innerC(State const &v, BondList const &bonds, State const &w) try { +// auto bw = zeros_like(w); +// apply(bonds, w, bw); +// return dotC(v, bw); +// } catch (Error const &error) { +// XDIAG_RETHROW(error); +// } -double inner(State const &v, Bond const &bond, State const &w) try { - return inner(v, BondList({bond}), w); -} catch (Error const &error) { - XDIAG_RETHROW(error); - return 0; -} +// double inner(State const &v, Bond const &bond, State const &w) try { +// return inner(v, BondList({bond}), w); +// } catch (Error const &error) { +// XDIAG_RETHROW(error); +// } -complex innerC(State const &v, Bond const &bond, State const &w) try { - return innerC(v, BondList({bond}), w); -} catch (Error const &error) { - XDIAG_RETHROW(error); - return 0; -} +// complex innerC(State const &v, Bond const &bond, State const &w) try { +// return innerC(v, BondList({bond}), w); +// } catch (Error const &error) { +// XDIAG_RETHROW(error); +// } State &operator*=(State &X, complex alpha) try { if (X.isreal()) { @@ -160,7 +166,6 @@ State &operator*=(State &X, complex alpha) try { return X; } catch (Error const &error) { XDIAG_RETHROW(error); - return X; } State &operator*=(State &X, double alpha) try { @@ -172,7 +177,6 @@ State &operator*=(State &X, double alpha) try { return X; } catch (Error const &error) { XDIAG_RETHROW(error); - return X; } State &operator/=(State &X, complex alpha) try { @@ -185,7 +189,6 @@ State &operator/=(State &X, complex alpha) try { return X; } catch (Error const &error) { XDIAG_RETHROW(error); - return X; } State &operator/=(State &X, double alpha) try { @@ -197,7 +200,6 @@ State &operator/=(State &X, double alpha) try { return X; } catch (Error const &error) { XDIAG_RETHROW(error); - return X; } // template @@ -375,7 +377,6 @@ double dot(block_variant_t const &block, arma::vec const &v, #endif } catch (Error const &error) { XDIAG_RETHROW(error); - return 0; } complex dot(block_variant_t const &block, arma::cx_vec const &v, @@ -393,7 +394,6 @@ complex dot(block_variant_t const &block, arma::cx_vec const &v, #endif } catch (Error const &error) { XDIAG_RETHROW(error); - return 0; } template @@ -401,7 +401,6 @@ double norm(block_variant_t const &block, arma::Col const &v) try { return std::sqrt(xdiag::real(dot(block, v, v))); } catch (Error const &error) { XDIAG_RETHROW(error); - return 0; } template double norm(block_variant_t const &, arma::Col const &); @@ -418,7 +417,6 @@ double norm1(block_variant_t const &block, arma::Col const &v) try { return nrm; } catch (Error const &error) { XDIAG_RETHROW(error); - return 0; } template double norm1(block_variant_t const &, arma::Col const &); @@ -435,7 +433,6 @@ double norminf(block_variant_t const &block, arma::Col const &v) try { return nrm; } catch (Error const &error) { XDIAG_RETHROW(error); - return 0; } template double norminf(block_variant_t const &, arma::Col const &); diff --git a/xdiag/algebra/algebra.hpp b/xdiag/algebra/algebra.hpp index 691f7591..6643e9f2 100644 --- a/xdiag/algebra/algebra.hpp +++ b/xdiag/algebra/algebra.hpp @@ -17,32 +17,17 @@ double dot(State const &v, State const &w); complex dotC(State const &v, State const &w); double inner(BondList const &bonds, State const &v); -complex innerC(BondList const &bonds, State const &v); -double inner(Bond const &bonds, State const &v); -complex innerC(Bond const &bonds, State const &v); +double inner(Bond const &bond, State const &v); -double inner(State const &v, BondList const &bonds, State const &w); -complex innerC(State const &v, BondList const &bonds, State const &w); -double inner(State const &v, Bond const &bonds, State const &w); -complex innerC(State const &v, Bond const &bonds, State const &w); +complex innerC(BondList const &bonds, State const &v); +complex innerC(Bond const &bond, State const &v); +// arithmetic operators State &operator*=(State &X, complex alpha); State &operator*=(State &X, double alpha); State &operator/=(State &X, complex alpha); State &operator/=(State &X, double alpha); -// // inner -// double inner(Bond const &bond, State const &v); -// double inner(State const &w, BondList const &bonds, State const &v); -// double inner(State const &w, Bond const &bond, State const &v); -// complex innerC(BondList const &bonds, State const &v); -// complex innerC(Bond const &bond, State const &v); -// complex innerC(State const &w, BondList const &bonds, State const &v); -// complex innerC(State const &w, Bond const &bond, State const &v); - -// State &operator*=(State &X, complex alpha); -// State &operator*=(State &X, double alpha); - double dot(block_variant_t const &block, arma::vec const &v, arma::vec const &w); complex dot(block_variant_t const &block, arma::cx_vec const &v, diff --git a/xdiag/algebra/apply.cpp b/xdiag/algebra/apply.cpp index 7875ace2..0c9accf2 100644 --- a/xdiag/algebra/apply.cpp +++ b/xdiag/algebra/apply.cpp @@ -26,19 +26,41 @@ template void apply(BondList const &, block_variant_t const &, void apply(BondList const &bonds, State const &v, State &w) try { if ((v.n_cols() == 1) && (w.n_cols() == 1)) { - if (bonds.isreal() && v.isreal() && w.isreal()) { - arma::vec vvec = v.vector(0, false); - arma::vec wvec = w.vector(0, false); - apply(bonds, v.block(), vvec, w.block(), wvec); - } else if (v.iscomplex() && w.iscomplex()) { - arma::cx_vec vvec = v.vectorC(0, false); - arma::cx_vec wvec = w.vectorC(0, false); - apply(bonds, v.block(), vvec, w.block(), wvec); + if (bonds.isreal()) { + if (v.isreal() && w.isreal()) { + arma::vec vvec = v.vector(0, false); + arma::vec wvec = w.vector(0, false); + apply(bonds, v.block(), vvec, w.block(), wvec); + } else if (v.isreal() && !w.isreal()) { + auto w2 = State(w.block(), true); + arma::vec vvec = v.vector(0, false); + arma::vec wvec = w2.vector(0, false); + apply(bonds, v.block(), vvec, w.block(), wvec); + w = w2; + } else if (!v.isreal() && w.isreal()) { + w.make_complex(); + arma::cx_vec vvec = v.vectorC(0, false); + arma::cx_vec wvec = w.vectorC(0, false); + apply(bonds, v.block(), vvec, w.block(), wvec); + } else if (!v.isreal() && !w.isreal()) { + arma::cx_vec vvec = v.vectorC(0, false); + arma::cx_vec wvec = w.vectorC(0, false); + apply(bonds, v.block(), vvec, w.block(), wvec); + } } else { - XDIAG_THROW( - "Apply operator only works if both states are complex or both states " - "are real and the operator is real. Consider Making the vectors " - "complex first by using method .make_complex()."); + if (v.isreal()) { + auto v2 = v; + v2.make_complex(); + w.make_complex(); + arma::cx_vec vvec = v2.vectorC(0, false); + arma::cx_vec wvec = w.vectorC(0, false); + apply(bonds, v.block(), vvec, w.block(), wvec); + } else { + w.make_complex(); + arma::cx_vec vvec = v.vectorC(0, false); + arma::cx_vec wvec = w.vectorC(0, false); + apply(bonds, v.block(), vvec, w.block(), wvec); + } } } else { XDIAG_THROW("Applying a BondList to a state with multiple " diff --git a/xdiag/algorithms/lanczos/eigs_lanczos.cpp b/xdiag/algorithms/lanczos/eigs_lanczos.cpp index c0e9ddee..4fdf0fbe 100644 --- a/xdiag/algorithms/lanczos/eigs_lanczos.cpp +++ b/xdiag/algorithms/lanczos/eigs_lanczos.cpp @@ -20,12 +20,12 @@ eigs_lanczos_result_t eigs_lanczos(BondList const &bonds, if (neigvals < 1) { XDIAG_THROW("Argument \"neigvals\" needs to be >= 1"); } - if (!bonds.ishermitian()) { - XDIAG_THROW("Input BondList is not hermitian"); - } + // if (!bonds.ishermitian()) { + // XDIAG_THROW("Input BondList is not hermitian"); + // } - bool cplx = bonds.iscomplex() || iscomplex(block) || force_complex || - state0.iscomplex(); + bool cplx = + !bonds.isreal() || iscomplex(block) || force_complex || !state0.isreal(); if (cplx) { state0.make_complex(); } @@ -118,11 +118,11 @@ eigs_lanczos_result_t eigs_lanczos(BondList const &bonds, if (neigvals < 1) { XDIAG_THROW("Argument \"neigvals\" needs to be >= 1"); } - if (!bonds.ishermitian()) { - XDIAG_THROW("Input BondList is not hermitian"); - } + // if (!bonds.ishermitian()) { + // XDIAG_THROW("Input BondList is not hermitian"); + // } - bool cplx = bonds.iscomplex() || iscomplex(block) || force_complex; + bool cplx = (!bonds.isreal()) || iscomplex(block) || force_complex; State state0(block, !cplx); fill(state0, RandomState(random_seed)); diff --git a/xdiag/algorithms/lanczos/eigvals_lanczos.cpp b/xdiag/algorithms/lanczos/eigvals_lanczos.cpp index 0f72425b..e1bcb9c3 100644 --- a/xdiag/algorithms/lanczos/eigvals_lanczos.cpp +++ b/xdiag/algorithms/lanczos/eigvals_lanczos.cpp @@ -23,13 +23,14 @@ eigvals_lanczos(BondList const &bonds, block_variant_t const &block, if (neigvals < 1) { XDIAG_THROW("Argument \"neigvals\" needs to be >= 1"); } - if (!bonds.ishermitian()) { - XDIAG_THROW("Input BondList is not hermitian"); - } - bool cplx = bonds.iscomplex() || iscomplex(block) || force_complex || - state0.iscomplex(); - if (state0.iscomplex() && !iscomplex(block)) { - Log(1, "warning: starting REAL block diagonalization with COMPLEX startstate"); + // if (!bonds.ishermitian()) { + // XDIAG_THROW("Input BondList is not hermitian"); + // } + bool cplx = + !bonds.isreal() || iscomplex(block) || force_complex || !state0.isreal(); + if (!state0.isreal() && !iscomplex(block)) { + Log(1, + "warning: starting REAL block diagonalization with COMPLEX startstate"); } auto converged = [neigvals, precision](Tmatrix const &tmat) -> bool { @@ -92,11 +93,11 @@ eigvals_lanczos(BondList const &bonds, block_variant_t const &block, if (neigvals < 1) { XDIAG_THROW("Argument \"neigvals\" needs to be >= 1"); } - if (!bonds.ishermitian()) { - XDIAG_THROW("Input BondList is not hermitian"); - } + // if (!bonds.ishermitian()) { + // XDIAG_THROW("Input BondList is not hermitian"); + // } - bool cplx = bonds.iscomplex() || iscomplex(block) || force_complex; + bool cplx = (!bonds.isreal()) || iscomplex(block) || force_complex; State state0(block, !cplx); fill(state0, RandomState(random_seed)); diff --git a/xdiag/blocks/electron/compile.cpp b/xdiag/blocks/electron/compile.cpp index 6f83fc29..92913088 100644 --- a/xdiag/blocks/electron/compile.cpp +++ b/xdiag/blocks/electron/compile.cpp @@ -4,11 +4,13 @@ namespace xdiag::electron { -BondList compile(BondList const &bonds, int64_t n_sites) try { +BondList compile(BondList const &bonds, int64_t n_sites, double precision) try { using namespace operators; BondList bonds_explicit = make_explicit(bonds); + BondList bonds_clean = clean_zeros(bonds_explicit, precision); + BondList bonds_compiled; - for (auto bond : bonds_explicit) { + for (auto bond : bonds_clean) { std::string type = bond.type(); // Exchange and Ising terms @@ -49,16 +51,16 @@ BondList compile(BondList const &bonds, int64_t n_sites) try { bonds_compiled += Bond("NUMBERDN", bond.coupling(), bond.sites()); } else if (type == "SZ") { check_bond(bond, n_sites, 1, false, "number"); - if (bond.coupling_is()) { - bonds_compiled += - Bond("NUMBERUP", 0.5 * bond.coupling(), bond.sites()); + if (bond.coupling().is()) { bonds_compiled += - Bond("NUMBERDN", -0.5 * bond.coupling(), bond.sites()); - } else if (bond.coupling_is()) { + Bond("NUMBERUP", 0.5 * bond.coupling().as(), bond.sites()); bonds_compiled += - Bond("NUMBERUP", 0.5 * bond.coupling(), bond.sites()); + Bond("NUMBERDN", -0.5 * bond.coupling().as(), bond.sites()); + } else if (bond.coupling().is()) { bonds_compiled += - Bond("NUMBERDN", -0.5 * bond.coupling(), bond.sites()); + Bond("NUMBERUP", 0.5 * bond.coupling().as(), bond.sites()); + bonds_compiled += Bond("NUMBERDN", -0.5 * bond.coupling().as(), + bond.sites()); } } @@ -80,7 +82,7 @@ BondList compile(BondList const &bonds, int64_t n_sites) try { } } // Set Hubbbbard U term again - if (bonds.coupling_defined("U")) { + if (bonds.defined("U")) { bonds_compiled["U"] = bonds["U"]; } diff --git a/xdiag/blocks/electron/compile.hpp b/xdiag/blocks/electron/compile.hpp index 0c7d983e..c3226333 100644 --- a/xdiag/blocks/electron/compile.hpp +++ b/xdiag/blocks/electron/compile.hpp @@ -5,10 +5,6 @@ namespace xdiag::electron { -const std::vector special_bond_types = { - "HB", "EXCHANGE", "ISING", "HOP", "HOPUP", "HOPDN", "NUMBER", - "NUMBERUP", "NUMBERDN", "SZ", "CDAGUP", "CDAGDN", "CUP", "CDN"}; - -BondList compile(BondList const &bonds, double precision = 1e-12); +BondList compile(BondList const &bonds, int64_t n_sites, double precision); } // namespace xdiag::electron diff --git a/xdiag/blocks/electron/dispatch.hpp b/xdiag/blocks/electron/dispatch.hpp index 81039d89..d7bedd28 100644 --- a/xdiag/blocks/electron/dispatch.hpp +++ b/xdiag/blocks/electron/dispatch.hpp @@ -9,7 +9,8 @@ namespace xdiag::electron { template inline void dispatch(BondList const &bonds, Electron const &block_in, - Electron const &block_out, Fill &&fill) try { + Electron const &block_out, Fill &&fill, + double zero_precision) try { using namespace basis::electron; auto const &basis_in = block_in.basis(); @@ -20,55 +21,55 @@ inline void dispatch(BondList const &bonds, Electron const &block_in, // uint16_t [&](BasisNp const &idx_in, BasisNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); + apply_terms(bonds, idx_in, idx_out, fill, zero_precision); }, [&](BasisNoNp const &idx_in, BasisNoNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); + apply_terms(bonds, idx_in, idx_out, fill, zero_precision); }, [&](BasisSymmetricNp const &idx_in, BasisSymmetricNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); + apply_terms(bonds, idx_in, idx_out, fill, zero_precision); }, [&](BasisSymmetricNoNp const &idx_in, BasisSymmetricNoNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); + apply_terms(bonds, idx_in, idx_out, fill, zero_precision); }, // uint32_t [&](BasisNp const &idx_in, BasisNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); + apply_terms(bonds, idx_in, idx_out, fill, zero_precision); }, [&](BasisNoNp const &idx_in, BasisNoNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); + apply_terms(bonds, idx_in, idx_out, fill, zero_precision); }, [&](BasisSymmetricNp const &idx_in, BasisSymmetricNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); + apply_terms(bonds, idx_in, idx_out, fill, zero_precision); }, [&](BasisSymmetricNoNp const &idx_in, BasisSymmetricNoNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); + apply_terms(bonds, idx_in, idx_out, fill, zero_precision); }, // uint64_t [&](BasisNp const &idx_in, BasisNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); + apply_terms(bonds, idx_in, idx_out, fill, zero_precision); }, [&](BasisNoNp const &idx_in, BasisNoNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); + apply_terms(bonds, idx_in, idx_out, fill, zero_precision); }, [&](BasisSymmetricNp const &idx_in, BasisSymmetricNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); + apply_terms(bonds, idx_in, idx_out, fill, zero_precision); }, [&](BasisSymmetricNoNp const &idx_in, BasisSymmetricNoNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); + apply_terms(bonds, idx_in, idx_out, fill, zero_precision); }, [&](auto const &idx_in, auto const &idx_out) { XDIAG_THROW("Invalid basis or combination of bases"); diff --git a/xdiag/blocks/electron/electron_apply.cpp b/xdiag/blocks/electron/electron_apply.cpp index 9fda7b14..abf83259 100644 --- a/xdiag/blocks/electron/electron_apply.cpp +++ b/xdiag/blocks/electron/electron_apply.cpp @@ -9,24 +9,22 @@ namespace xdiag { template void apply(BondList const &bonds, Electron const &block_in, arma::Col const &vec_in, Electron const &block_out, - arma::Col &vec_out) try { + arma::Col &vec_out, double zero_precision) try { vec_out.zeros(); - - BondList bondsc = electron::compile(bonds, 1e-12); auto fill = [&](int64_t idx_in, int64_t idx_out, coeff_t val) { return fill_apply(vec_in, vec_out, idx_in, idx_out, val); }; - electron::dispatch(bondsc, block_in, block_out, fill); + electron::dispatch(bonds, block_in, block_out, fill, zero_precision); } catch (Error const &e) { XDIAG_RETHROW(e); } template void apply(BondList const &, Electron const &, arma::Col const &, Electron const &, - arma::Col &); + arma::Col &, double); template void apply(BondList const &, Electron const &, arma::Col const &, Electron const &, - arma::Col &); + arma::Col &, double); } // namespace xdiag diff --git a/xdiag/blocks/electron/electron_apply.hpp b/xdiag/blocks/electron/electron_apply.hpp index e62725e5..975d4d5d 100644 --- a/xdiag/blocks/electron/electron_apply.hpp +++ b/xdiag/blocks/electron/electron_apply.hpp @@ -10,6 +10,6 @@ namespace xdiag { template void apply(BondList const &bonds, Electron const &block_in, arma::Col const &vec_in, Electron const &block_out, - arma::Col &vec_out); + arma::Col &vec_out, double zero_precision = 1e-12); } // namespace xdiag diff --git a/xdiag/blocks/electron/electron_matrix.cpp b/xdiag/blocks/electron/electron_matrix.cpp index 68bbe2d3..5b4ff937 100644 --- a/xdiag/blocks/electron/electron_matrix.cpp +++ b/xdiag/blocks/electron/electron_matrix.cpp @@ -9,40 +9,35 @@ namespace xdiag { template arma::Mat matrix_gen(BondList const &bonds, Electron const &block_in, - Electron const &block_out) try { - BondList bondsc = electron::compile(bonds, 1e-12); - + Electron const &block_out, + double zero_precision) try { int64_t m = block_out.size(); int64_t n = block_in.dim(); arma::Mat mat(m, n, arma::fill::zeros); - matrix_gen(mat.memptr(), bonds, block_in, block_out); + matrix_gen(mat.memptr(), bonds, block_in, block_out, zero_precision); return mat; } catch (Error const &error) { XDIAG_RETHROW(error); - return arma::Mat(); } arma::mat matrix(BondList const &bonds, Electron const &block_in, - Electron const &block_out) try { - return matrix_gen(bonds, block_in, block_out); + Electron const &block_out, double zero_precision) try { + return matrix_gen(bonds, block_in, block_out, zero_precision); } catch (Error const &error) { XDIAG_RETHROW(error); - return arma::mat(); } arma::cx_mat matrixC(BondList const &bonds, Electron const &block_in, - Electron const &block_out) try { - return matrix_gen(bonds, block_in, block_out); + Electron const &block_out, double zero_precision) try { + return matrix_gen(bonds, block_in, block_out, zero_precision); } catch (Error const &error) { XDIAG_RETHROW(error); - return arma::cx_mat(); } template void matrix_gen(coeff_t *mat, BondList const &bonds, Electron const &block_in, - Electron const &block_out) try { - BondList bondsc = electron::compile(bonds, 1e-12); + Electron const &block_out, double zero_precision) try { int64_t m = block_out.size(); int64_t n = block_in.size(); std::fill(mat, mat + m * n, 0); @@ -51,21 +46,21 @@ void matrix_gen(coeff_t *mat, BondList const &bonds, Electron const &block_in, return fill_matrix(mat, idx_in, idx_out, m, val); }; - electron::dispatch(bondsc, block_in, block_out, fill); + electron::dispatch(bonds, block_in, block_out, fill, zero_precision); } catch (Error const &error) { XDIAG_RETHROW(error); } void matrix(double *mat, BondList const &bonds, Electron const &block_in, - Electron const &block_out) try { - return matrix_gen(mat, bonds, block_in, block_out); + Electron const &block_out, double zero_precision) try { + return matrix_gen(mat, bonds, block_in, block_out, zero_precision); } catch (Error const &error) { XDIAG_RETHROW(error); } void matrixC(complex *mat, BondList const &bonds, Electron const &block_in, - Electron const &block_out) try { - return matrix_gen(mat, bonds, block_in, block_out); + Electron const &block_out, double zero_precision) try { + return matrix_gen(mat, bonds, block_in, block_out, zero_precision); } catch (Error const &error) { XDIAG_RETHROW(error); } diff --git a/xdiag/blocks/electron/electron_matrix.hpp b/xdiag/blocks/electron/electron_matrix.hpp index af43eeeb..50d96742 100644 --- a/xdiag/blocks/electron/electron_matrix.hpp +++ b/xdiag/blocks/electron/electron_matrix.hpp @@ -8,16 +8,16 @@ namespace xdiag { arma::mat matrix(BondList const &bonds, Electron const &block_in, - Electron const &block_out); + Electron const &block_out, double zero_precision=1e-12); arma::cx_mat matrixC(BondList const &bonds, Electron const &block_in, - Electron const &block_out); + Electron const &block_out, double zero_precision=1e-12); // Developer functions void matrix(double *mat, BondList const &bonds, Electron const &block_in, - Electron const &block_out); + Electron const &block_out, double zero_precision=1e-12); void matrixC(complex *mat, BondList const &bonds, Electron const &block_in, - Electron const &block_out); + Electron const &block_out, double zero_precision=1e-12); } // namespace xdiag diff --git a/xdiag/blocks/electron/terms/apply_exchange.hpp b/xdiag/blocks/electron/terms/apply_exchange.hpp index 9620588f..69a0e600 100644 --- a/xdiag/blocks/electron/terms/apply_exchange.hpp +++ b/xdiag/blocks/electron/terms/apply_exchange.hpp @@ -39,7 +39,9 @@ void apply_exchange(Bond const &bond, Basis &&basis, Filler &&fill) { assert(bond.size() == 2); assert(bond.sites_disjoint()); - coeff_t J = bond.coupling(); + Coupling cpl = bond.coupling(); + assert(cpl.isexplicit() && !cpl.ismatrix()); + coeff_t J = cpl.as(); int64_t s1 = bond[0]; int64_t s2 = bond[1]; diff --git a/xdiag/blocks/electron/terms/apply_hopping.hpp b/xdiag/blocks/electron/terms/apply_hopping.hpp index 1c7a8629..2ee94e41 100644 --- a/xdiag/blocks/electron/terms/apply_hopping.hpp +++ b/xdiag/blocks/electron/terms/apply_hopping.hpp @@ -26,7 +26,10 @@ void apply_hopping(Bond const &bond, Basis &&basis, Fill &&fill) { int64_t l = std::min(s1, s2); int64_t u = std::max(s1, s2); bit_t fermimask = (((bit_t)1 << (u - l - 1)) - 1) << (l + 1); - coeff_t t = bond.coupling(); + + Coupling cpl = bond.coupling(); + assert(cpl.isexplicit() && !cpl.ismatrix()); + coeff_t t = cpl.as(); auto non_zero_term = [&flipmask](bit_t const &spins) -> bool { return bits::popcnt(spins & flipmask) & 1; diff --git a/xdiag/blocks/electron/terms/apply_ising.hpp b/xdiag/blocks/electron/terms/apply_ising.hpp index 0fadfecd..b8290262 100644 --- a/xdiag/blocks/electron/terms/apply_ising.hpp +++ b/xdiag/blocks/electron/terms/apply_ising.hpp @@ -13,7 +13,9 @@ void apply_ising(Bond const &bond, Basis &&basis, Fill &&fill) { assert(bond.size() == 2); assert(bond.sites_disjoint()); - coeff_t J = bond.coupling(); + Coupling cpl = bond.coupling(); + assert(cpl.isexplicit() && !cpl.ismatrix()); + coeff_t J = cpl.as(); int64_t s1 = bond[0]; int64_t s2 = bond[1]; diff --git a/xdiag/blocks/electron/terms/apply_number.hpp b/xdiag/blocks/electron/terms/apply_number.hpp index f9d84dcf..56bd886b 100644 --- a/xdiag/blocks/electron/terms/apply_number.hpp +++ b/xdiag/blocks/electron/terms/apply_number.hpp @@ -12,7 +12,9 @@ void apply_number(Bond const &bond, Basis &&basis, Fill fill) { std::string type = bond.type(); assert((type == "NUMBERUP") || (type == "NUMBERDN")); - coeff_t mu = bond.coupling(); + Coupling cpl = bond.coupling(); + assert(cpl.isexplicit() && !cpl.ismatrix()); + coeff_t mu = cpl.as(); int64_t s = bond[0]; bit_t mask = (bit_t)1 << s; diff --git a/xdiag/blocks/electron/terms/apply_raise_lower.hpp b/xdiag/blocks/electron/terms/apply_raise_lower.hpp index 09cff29b..00e34110 100644 --- a/xdiag/blocks/electron/terms/apply_raise_lower.hpp +++ b/xdiag/blocks/electron/terms/apply_raise_lower.hpp @@ -23,7 +23,10 @@ void apply_raise_lower(Bond const &bond, BasisIn &&basis_in, int64_t s = bond[0]; bit_t site_mask = (bit_t)1 << s; bit_t fermi_mask = site_mask - 1; - coeff_t c = bond.coupling(); + + Coupling cpl = bond.coupling(); + assert(cpl.isexplicit() && !cpl.ismatrix()); + coeff_t c = cpl.as(); // Raising operators if ((type == "CDAGUP") || (type == "CDAGDN")) { diff --git a/xdiag/blocks/electron/terms/apply_terms.hpp b/xdiag/blocks/electron/terms/apply_terms.hpp index ee09946a..495fbc7f 100644 --- a/xdiag/blocks/electron/terms/apply_terms.hpp +++ b/xdiag/blocks/electron/terms/apply_terms.hpp @@ -15,9 +15,11 @@ namespace xdiag::electron { template void apply_terms(BondList const &bonds, BasisIn const &basis_in, - BasisOut const &basis_out, Fill &fill) try { + BasisOut const &basis_out, Fill &fill, + double zero_precision) try { - BondList bonds_compiled = electron::compile(bonds); + BondList bonds_compiled = + electron::compile(bonds, basis_in.n_sites(), zero_precision); // Separate bond acting on ups, dns, diagonally or fully mixing BondList bonds_ups; @@ -50,11 +52,12 @@ void apply_terms(BondList const &bonds, BasisIn const &basis_in, } } - if (bonds.coupling_defined("U")) { - coupling_t U = bonds["U"]; - if (std::holds_alternative(U)){ - electron::apply_u(std::get(U), - basis_in, fill); + if (bonds.defined("U")) { + Coupling cpl = bonds["U"]; + + if (cpl.is()) { + double U = cpl.as(); + electron::apply_u(U, basis_in, fill); } else { XDIAG_THROW("Coupling U must either be a real number"); } diff --git a/xdiag/blocks/spinhalf/compile.cpp b/xdiag/blocks/spinhalf/compile.cpp index b7bd1669..2bce5156 100644 --- a/xdiag/blocks/spinhalf/compile.cpp +++ b/xdiag/blocks/spinhalf/compile.cpp @@ -10,10 +10,11 @@ namespace xdiag::spinhalf { BondList compile(BondList const &bonds, int64_t n_sites, double precision) try { using namespace operators; - BondList bonds_explicit = make_explicit(bonds); + BondList bonds_clean = clean_zeros(bonds_explicit, precision); + BondList bonds_compiled; - for (auto bond : bonds_explicit) { + for (auto bond : bonds_clean) { if (bond.ismatrix()) { BondList bonds_nb = non_branching_bonds(bond, precision); @@ -51,7 +52,6 @@ BondList compile(BondList const &bonds, int64_t n_sites, double precision) try { return bonds_compiled; } catch (Error const &e) { XDIAG_RETHROW(e); - return BondList(); } } // namespace xdiag::spinhalf diff --git a/xdiag/blocks/spinhalf/compile.hpp b/xdiag/blocks/spinhalf/compile.hpp index 2337c113..627b2f9b 100644 --- a/xdiag/blocks/spinhalf/compile.hpp +++ b/xdiag/blocks/spinhalf/compile.hpp @@ -5,7 +5,6 @@ namespace xdiag::spinhalf { -BondList compile(BondList const &bonds, int64_t n_sites, - double precision = 1e-12); +BondList compile(BondList const &bonds, int64_t n_sites, double precision); } // namespace xdiag::spinhalf diff --git a/xdiag/blocks/spinhalf/dispatch.hpp b/xdiag/blocks/spinhalf/dispatch.hpp index 098b1249..802fad12 100644 --- a/xdiag/blocks/spinhalf/dispatch.hpp +++ b/xdiag/blocks/spinhalf/dispatch.hpp @@ -11,75 +11,88 @@ namespace xdiag::spinhalf { template inline void dispatch(BondList const &bonds, Spinhalf const &block_in, - Spinhalf const &block_out, fill_f fill) try { + Spinhalf const &block_out, fill_f fill, + double zero_precision) try { using namespace basis::spinhalf; auto const &basis_in = block_in.basis(); auto const &basis_out = block_out.basis(); std::visit( - overload{ - // uint32_t - [&](BasisSz const &idx_in, - BasisSz const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - [&](BasisNoSz const &idx_in, - BasisNoSz const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - [&](BasisSymmetricSz const &idx_in, - BasisSymmetricSz const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - [&](BasisSymmetricNoSz const &idx_in, - BasisSymmetricNoSz const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, + overload{// uint32_t + [&](BasisSz const &idx_in, + BasisSz const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + [&](BasisNoSz const &idx_in, + BasisNoSz const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + [&](BasisSymmetricSz const &idx_in, + BasisSymmetricSz const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + [&](BasisSymmetricNoSz const &idx_in, + BasisSymmetricNoSz const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, - // uint64_t - [&](BasisSz const &idx_in, - BasisSz const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - [&](BasisNoSz const &idx_in, - BasisNoSz const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - [&](BasisSymmetricSz const &idx_in, - BasisSymmetricSz const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - [&](BasisSymmetricNoSz const &idx_in, - BasisSymmetricNoSz const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - [&](BasisSublattice const &idx_in, - BasisSublattice const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - [&](BasisSublattice const &idx_in, - BasisSublattice const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - [&](BasisSublattice const &idx_in, - BasisSublattice const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - [&](BasisSublattice const &idx_in, - BasisSublattice const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - [&](BasisSublattice const &idx_in, - BasisSublattice const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, + // uint64_t + [&](BasisSz const &idx_in, + BasisSz const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + [&](BasisNoSz const &idx_in, + BasisNoSz const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + [&](BasisSymmetricSz const &idx_in, + BasisSymmetricSz const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + [&](BasisSymmetricNoSz const &idx_in, + BasisSymmetricNoSz const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + [&](BasisSublattice const &idx_in, + BasisSublattice const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + [&](BasisSublattice const &idx_in, + BasisSublattice const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + [&](BasisSublattice const &idx_in, + BasisSublattice const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + [&](BasisSublattice const &idx_in, + BasisSublattice const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + [&](BasisSublattice const &idx_in, + BasisSublattice const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, - [&](auto const &idx_in, auto const &idx_out) { - XDIAG_THROW("Invalid basis or combination of bases"); - (void)idx_in; - (void)idx_out; - }}, + [&](auto const &idx_in, auto const &idx_out) { + XDIAG_THROW("Invalid basis or combination of bases"); + (void)idx_in; + (void)idx_out; + }}, basis_in, basis_out); } catch (Error const &e) { XDIAG_RETHROW(e); diff --git a/xdiag/blocks/spinhalf/spinhalf_apply.cpp b/xdiag/blocks/spinhalf/spinhalf_apply.cpp index 8bd7faf5..20b23c69 100644 --- a/xdiag/blocks/spinhalf/spinhalf_apply.cpp +++ b/xdiag/blocks/spinhalf/spinhalf_apply.cpp @@ -9,24 +9,24 @@ namespace xdiag { template void apply(BondList const &bonds, Spinhalf const &block_in, arma::Col const &vec_in, Spinhalf const &block_out, - arma::Col &vec_out) try { + arma::Col &vec_out, double zero_precision) try { vec_out.zeros(); - BondList bondsc = spinhalf::compile(bonds, 1e-12); auto fill = [&](int64_t idx_in, int64_t idx_out, coeff_t val) { return fill_apply(vec_in, vec_out, idx_in, idx_out, val); }; - spinhalf::dispatch(bondsc, block_in, block_out, fill); + + spinhalf::dispatch(bonds, block_in, block_out, fill, zero_precision); } catch (Error const &e) { XDIAG_RETHROW(e); } template void apply(BondList const &, Spinhalf const &, arma::Col const &, Spinhalf const &, - arma::Col &); + arma::Col &, double); template void apply(BondList const &, Spinhalf const &, arma::Col const &, Spinhalf const &, - arma::Col &); + arma::Col &, double); } // namespace xdiag diff --git a/xdiag/blocks/spinhalf/spinhalf_apply.hpp b/xdiag/blocks/spinhalf/spinhalf_apply.hpp index c6466ff9..7686c0de 100644 --- a/xdiag/blocks/spinhalf/spinhalf_apply.hpp +++ b/xdiag/blocks/spinhalf/spinhalf_apply.hpp @@ -11,6 +11,6 @@ namespace xdiag { template void apply(BondList const &bonds, Spinhalf const &block_in, arma::Col const &vec_in, Spinhalf const &block_out, - arma::Col &vec_out); + arma::Col &vec_out, double zero_precision = 1e-12); } // namespace xdiag diff --git a/xdiag/blocks/spinhalf/spinhalf_matrix.cpp b/xdiag/blocks/spinhalf/spinhalf_matrix.cpp index 06ced342..2ce85999 100644 --- a/xdiag/blocks/spinhalf/spinhalf_matrix.cpp +++ b/xdiag/blocks/spinhalf/spinhalf_matrix.cpp @@ -8,39 +8,34 @@ namespace xdiag { template arma::Mat matrix_gen(BondList const &bonds, Spinhalf const &block_in, - Spinhalf const &block_out) try { - BondList bondsc = spinhalf::compile(bonds, 1e-12); - + Spinhalf const &block_out, + double zero_precision) try { int64_t m = block_out.size(); int64_t n = block_in.dim(); arma::Mat mat(m, n, arma::fill::zeros); - matrix_gen(mat.memptr(), bonds, block_in, block_out); + matrix_gen(mat.memptr(), bonds, block_in, block_out, zero_precision); return mat; } catch (Error const &e) { XDIAG_RETHROW(e); - return arma::Mat(); } arma::mat matrix(BondList const &bonds, Spinhalf const &block_in, - Spinhalf const &block_out) try { - return matrix_gen(bonds, block_in, block_out); + Spinhalf const &block_out, double zero_precision) try { + return matrix_gen(bonds, block_in, block_out, zero_precision); } catch (Error const &e) { XDIAG_RETHROW(e); - return arma::mat(); } arma::cx_mat matrixC(BondList const &bonds, Spinhalf const &block_in, - Spinhalf const &block_out) try { - return matrix_gen(bonds, block_in, block_out); + Spinhalf const &block_out, double zero_precision) try { + return matrix_gen(bonds, block_in, block_out, zero_precision); } catch (Error const &e) { XDIAG_RETHROW(e); - return arma::cx_mat(); } template void matrix_gen(coeff_t *mat, BondList const &bonds, Spinhalf const &block_in, - Spinhalf const &block_out) try { - BondList bondsc = spinhalf::compile(bonds, 1e-12); + Spinhalf const &block_out, double zero_precision) try { int64_t m = block_out.size(); int64_t n = block_in.size(); std::fill(mat, mat + m * n, 0); @@ -48,21 +43,21 @@ void matrix_gen(coeff_t *mat, BondList const &bonds, Spinhalf const &block_in, auto fill = [&](int64_t idx_in, int64_t idx_out, coeff_t val) { return fill_matrix(mat, idx_in, idx_out, m, val); }; - spinhalf::dispatch(bondsc, block_in, block_out, fill); + spinhalf::dispatch(bonds, block_in, block_out, fill, zero_precision); } catch (Error const &e) { XDIAG_RETHROW(e); } void matrix(double *mat, BondList const &bonds, Spinhalf const &block_in, - Spinhalf const &block_out) try { - return matrix_gen(mat, bonds, block_in, block_out); + Spinhalf const &block_out, double zero_precision) try { + return matrix_gen(mat, bonds, block_in, block_out, zero_precision); } catch (Error const &e) { XDIAG_RETHROW(e); } void matrixC(complex *mat, BondList const &bonds, Spinhalf const &block_in, - Spinhalf const &block_out) try { - return matrix_gen(mat, bonds, block_in, block_out); + Spinhalf const &block_out, double zero_precision) try { + return matrix_gen(mat, bonds, block_in, block_out, zero_precision); } catch (Error const &e) { XDIAG_RETHROW(e); } diff --git a/xdiag/blocks/spinhalf/spinhalf_matrix.hpp b/xdiag/blocks/spinhalf/spinhalf_matrix.hpp index 48355fdc..beab5a23 100644 --- a/xdiag/blocks/spinhalf/spinhalf_matrix.hpp +++ b/xdiag/blocks/spinhalf/spinhalf_matrix.hpp @@ -2,23 +2,23 @@ #include +#include #include #include -#include namespace xdiag { arma::mat matrix(BondList const &bonds, Spinhalf const &block_in, - Spinhalf const &block_out); + Spinhalf const &block_out, double zero_precision = 1e-12); arma::cx_mat matrixC(BondList const &bonds, Spinhalf const &block_in, - Spinhalf const &block_out); + Spinhalf const &block_out, double zero_precision = 1e-12); // Developer functions void matrix(double *mat, BondList const &bonds, Spinhalf const &block_in, - Spinhalf const &block_out); + Spinhalf const &block_out, double zero_precision = 1e-12); void matrixC(complex *mat, BondList const &bonds, Spinhalf const &block_in, - Spinhalf const &block_out); + Spinhalf const &block_out, double zero_precision = 1e-12); } // namespace xdiag diff --git a/xdiag/blocks/spinhalf/terms/apply_exchange.hpp b/xdiag/blocks/spinhalf/terms/apply_exchange.hpp index d1e23e4f..9aae5b32 100644 --- a/xdiag/blocks/spinhalf/terms/apply_exchange.hpp +++ b/xdiag/blocks/spinhalf/terms/apply_exchange.hpp @@ -21,7 +21,10 @@ void apply_exchange(Bond const &bond, BasisIn &&basis_in, BasisOut &&basis_out, assert(bond.size() == 2); assert(bond.sites_disjoint()); - coeff_t J = bond.coupling(); + Coupling cpl = bond.coupling(); + assert(cpl.isexplicit() && !cpl.ismatrix()); + + coeff_t J = cpl.as(); int64_t s1 = bond[0]; int64_t s2 = bond[1]; bit_t flipmask = ((bit_t)1 << s1) | ((bit_t)1 << s2); diff --git a/xdiag/blocks/spinhalf/terms/apply_ising.hpp b/xdiag/blocks/spinhalf/terms/apply_ising.hpp index b74e5a73..78946028 100644 --- a/xdiag/blocks/spinhalf/terms/apply_ising.hpp +++ b/xdiag/blocks/spinhalf/terms/apply_ising.hpp @@ -17,12 +17,14 @@ void apply_ising(Bond const &bond, BasisIn &&basis_in, BasisOut &&basis_out, assert(bond.size() == 2); assert(bond.sites_disjoint()); + Coupling cpl = bond.coupling(); + assert(cpl.isexplicit() && !cpl.ismatrix()); + + coeff_t J = cpl.as(); int64_t s1 = bond[0]; int64_t s2 = bond[1]; bit_t mask = ((bit_t)1 << s1) | ((bit_t)1 << s2); - coeff_t J = bond.coupling(); - coeff_t val_same = J / 4.0; coeff_t val_diff = -J / 4.0; diff --git a/xdiag/blocks/spinhalf/terms/apply_scalar_chirality.hpp b/xdiag/blocks/spinhalf/terms/apply_scalar_chirality.hpp index eab576e6..61e7f5ed 100644 --- a/xdiag/blocks/spinhalf/terms/apply_scalar_chirality.hpp +++ b/xdiag/blocks/spinhalf/terms/apply_scalar_chirality.hpp @@ -21,7 +21,10 @@ void apply_scalar_chirality(Bond const &bond, BasisIn &&basis_in, assert(bond.sites_disjoint()); assert(iscomplex()); - complex J = bond.coupling(); + Coupling cpl = bond.coupling(); + assert(cpl.isexplicit() && !cpl.ismatrix()); + + complex J = cpl.as(); coeff_t Jquarter = 0.; coeff_t Jquarter_conj = 0.; if constexpr (isreal()) { diff --git a/xdiag/blocks/spinhalf/terms/apply_spsm.hpp b/xdiag/blocks/spinhalf/terms/apply_spsm.hpp index 88b1f7e0..efaed94c 100644 --- a/xdiag/blocks/spinhalf/terms/apply_spsm.hpp +++ b/xdiag/blocks/spinhalf/terms/apply_spsm.hpp @@ -21,7 +21,10 @@ void apply_spsm(Bond const &bond, BasisIn &&basis_in, BasisOut &&basis_out, assert((bond.type() == "S+") || (bond.type() == "S-")); assert(bond.size() == 1); - coeff_t J = bond.coupling(); + Coupling cpl = bond.coupling(); + assert(cpl.isexplicit() && !cpl.ismatrix()); + + coeff_t J = cpl.as(); int64_t s = bond[0]; bit_t mask = ((bit_t)1 << s); diff --git a/xdiag/blocks/spinhalf/terms/apply_sz.hpp b/xdiag/blocks/spinhalf/terms/apply_sz.hpp index e100c8c3..91f6f6e7 100644 --- a/xdiag/blocks/spinhalf/terms/apply_sz.hpp +++ b/xdiag/blocks/spinhalf/terms/apply_sz.hpp @@ -16,10 +16,13 @@ void apply_sz(Bond const &bond, BasisIn &&basis_in, BasisOut &&basis_out, assert(bond.type_defined() && (bond.type() == "SZ")); assert(bond.size() == 1); + Coupling cpl = bond.coupling(); + assert(cpl.isexplicit() && !cpl.ismatrix()); + coeff_t H = cpl.as(); + int64_t s = bond[0]; bit_t mask = ((bit_t)1 << s); - coeff_t H = bond.coupling(); coeff_t val_up = H / 2.; coeff_t val_dn = -H / 2.; diff --git a/xdiag/blocks/spinhalf/terms/apply_terms.hpp b/xdiag/blocks/spinhalf/terms/apply_terms.hpp index 23f1f9c6..149f14b2 100644 --- a/xdiag/blocks/spinhalf/terms/apply_terms.hpp +++ b/xdiag/blocks/spinhalf/terms/apply_terms.hpp @@ -15,8 +15,12 @@ namespace xdiag::spinhalf { template void apply_terms(BondList const &bonds, BasisIn const &basis_in, - BasisOut const &basis_out, Fill &fill) try { - for (auto bond : bonds) { + BasisOut const &basis_out, Fill &fill, + double zero_precision) try { + BondList bonds_compiled = + spinhalf::compile(bonds, basis_in.n_sites(), zero_precision); + + for (auto bond : bonds_compiled) { if (bond.type() == "EXCHANGE") { spinhalf::apply_exchange(bond, basis_in, diff --git a/xdiag/blocks/spinhalf_distributed/spinhalf_distributed_apply.cpp b/xdiag/blocks/spinhalf_distributed/spinhalf_distributed_apply.cpp index 0ff9e7ad..87363413 100644 --- a/xdiag/blocks/spinhalf_distributed/spinhalf_distributed_apply.cpp +++ b/xdiag/blocks/spinhalf_distributed/spinhalf_distributed_apply.cpp @@ -8,10 +8,10 @@ namespace xdiag { template void apply(BondList const &bonds, tJ const &block_in, arma::Col const &vec_in, tJ const &block_out, - arma::Col &vec_out) try { + arma::Col &vec_out, double zero_precision) try { vec_out.zeros(); - BondList bondsc = tj::compile(bonds, 1e-12); + BondList bondsc = spinhalf::compile(bonds, block_in, zero_precision); auto fill = [&](int64_t idx_in, int64_t idx_out, coeff_t val) { return fill_apply(vec_in, vec_out, idx_in, idx_out, val); }; @@ -22,10 +22,10 @@ void apply(BondList const &bonds, tJ const &block_in, template void apply(BondList const &, tJ const &, arma::Col const &, tJ const &, - arma::Col &); + arma::Col &, double); template void apply(BondList const &, tJ const &, arma::Col const &, tJ const &, - arma::Col &); + arma::Col &, double); } // namespace xdiag diff --git a/xdiag/blocks/tj/compile.cpp b/xdiag/blocks/tj/compile.cpp index 84b35c1e..899801e7 100644 --- a/xdiag/blocks/tj/compile.cpp +++ b/xdiag/blocks/tj/compile.cpp @@ -6,11 +6,13 @@ namespace xdiag::tj { -BondList compile(BondList const &bonds, int64_t n_sites) try { +BondList compile(BondList const &bonds, int64_t n_sites, double precision) try { using namespace operators; BondList bonds_explicit = make_explicit(bonds); + BondList bonds_clean = clean_zeros(bonds_explicit, precision); + BondList bonds_compiled; - for (auto bond : bonds_explicit) { + for (auto bond : bonds_clean) { std::string type = bond.type(); // Exchange and Ising terms @@ -59,16 +61,16 @@ BondList compile(BondList const &bonds, int64_t n_sites) try { bonds_compiled += Bond("NUMBERDN", bond.coupling(), bond.sites()); } else if (type == "SZ") { check_bond(bond, n_sites, 1, false, "number"); - if (bond.coupling_is()) { - bonds_compiled += - Bond("NUMBERUP", 0.5 * bond.coupling(), bond.sites()); + if (bond.coupling().is()) { bonds_compiled += - Bond("NUMBERDN", -0.5 * bond.coupling(), bond.sites()); - } else if (bond.coupling_is()) { + Bond("NUMBERUP", 0.5 * bond.coupling().as(), bond.sites()); bonds_compiled += - Bond("NUMBERUP", 0.5 * bond.coupling(), bond.sites()); + Bond("NUMBERDN", -0.5 * bond.coupling().as(), bond.sites()); + } else if (bond.coupling().is()) { bonds_compiled += - Bond("NUMBERDN", -0.5 * bond.coupling(), bond.sites()); + Bond("NUMBERUP", 0.5 * bond.coupling().as(), bond.sites()); + bonds_compiled += Bond("NUMBERDN", -0.5 * bond.coupling().as(), + bond.sites()); } } diff --git a/xdiag/blocks/tj/compile.hpp b/xdiag/blocks/tj/compile.hpp index d352efd8..165709b8 100644 --- a/xdiag/blocks/tj/compile.hpp +++ b/xdiag/blocks/tj/compile.hpp @@ -6,6 +6,6 @@ namespace xdiag::tj { -BondList compile(BondList const &bonds, double precision = 1e-12); +BondList compile(BondList const &bonds, int64_t n_sites, double precision); } // namespace xdiag::tj diff --git a/xdiag/blocks/tj/dispatch.hpp b/xdiag/blocks/tj/dispatch.hpp index d1f7c9f8..a0499ae5 100644 --- a/xdiag/blocks/tj/dispatch.hpp +++ b/xdiag/blocks/tj/dispatch.hpp @@ -9,47 +9,53 @@ namespace xdiag::tj { template inline void dispatch(BondList const &bonds, tJ const &block_in, - tJ const &block_out, fill_f &&fill) try { + tJ const &block_out, fill_f &&fill, + double zero_precision) try { using namespace basis::tj; auto const &basis_in = block_in.basis(); auto const &basis_out = block_out.basis(); std::visit( - overload{ - // uint16_t - [&](BasisNp const &idx_in, - BasisNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - [&](BasisSymmetricNp const &idx_in, - BasisSymmetricNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - // uint32_t - [&](BasisNp const &idx_in, - BasisNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - [&](BasisSymmetricNp const &idx_in, - BasisSymmetricNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - // uint64_t - [&](BasisNp const &idx_in, - BasisNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - [&](BasisSymmetricNp const &idx_in, - BasisSymmetricNp const &idx_out) { - apply_terms(bonds, idx_in, idx_out, fill); - }, - [&](auto const &idx_in, auto const &idx_out) { - XDIAG_THROW("Invalid basis or combination of bases"); - (void)idx_in; - (void)idx_out; - }}, + overload{// uint16_t + [&](BasisNp const &idx_in, + BasisNp const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + [&](BasisSymmetricNp const &idx_in, + BasisSymmetricNp const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + // uint32_t + [&](BasisNp const &idx_in, + BasisNp const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + [&](BasisSymmetricNp const &idx_in, + BasisSymmetricNp const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + // uint64_t + [&](BasisNp const &idx_in, + BasisNp const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + [&](BasisSymmetricNp const &idx_in, + BasisSymmetricNp const &idx_out) { + apply_terms(bonds, idx_in, idx_out, + fill, zero_precision); + }, + [&](auto const &idx_in, auto const &idx_out) { + XDIAG_THROW("Invalid basis or combination of bases"); + (void)idx_in; + (void)idx_out; + }}, basis_in, basis_out); -} catch (Error const& e) { +} catch (Error const &e) { XDIAG_RETHROW(e); } diff --git a/xdiag/blocks/tj/terms/apply_exchange.hpp b/xdiag/blocks/tj/terms/apply_exchange.hpp index 0fe68ae9..c77cdfc0 100644 --- a/xdiag/blocks/tj/terms/apply_exchange.hpp +++ b/xdiag/blocks/tj/terms/apply_exchange.hpp @@ -17,12 +17,15 @@ void apply_exchange(Bond const &bond, Basis &&basis, Filler &&fill) { std::string type = bond.type(); assert(type == "EXCHANGE"); - int64_t s1 = bond[0]; - int64_t s2 = bond[1]; - coeff_t J = bond.coupling(); + Coupling cpl = bond.coupling(); + assert(cpl.isexplicit() && !cpl.ismatrix()); + coeff_t J = cpl.as(); coeff_t Jhalf = J / 2.; coeff_t Jhalf_conj = conj(Jhalf); + int64_t s1 = bond[0]; + int64_t s2 = bond[1]; + // Prepare bitmasks bit_t flipmask = ((bit_t)1 << s1) | ((bit_t)1 << s2); int64_t l = std::min(s1, s2); diff --git a/xdiag/blocks/tj/terms/apply_hopping.hpp b/xdiag/blocks/tj/terms/apply_hopping.hpp index 1d8b2579..36872bb4 100644 --- a/xdiag/blocks/tj/terms/apply_hopping.hpp +++ b/xdiag/blocks/tj/terms/apply_hopping.hpp @@ -1,10 +1,10 @@ #pragma once #include -#include -#include #include #include +#include +#include namespace xdiag::tj { @@ -19,13 +19,16 @@ void apply_hopping(Bond const &bond, Basis &&basis, Fill &&fill) { std::string type = bond.type(); assert((type == "HOPUP") || (type == "HOPDN")); + Coupling cpl = bond.coupling(); + assert(cpl.isexplicit() && !cpl.ismatrix()); + coeff_t t = cpl.as(); + int64_t s1 = bond[0]; int64_t s2 = bond[1]; bit_t flipmask = ((bit_t)1 << s1) | ((bit_t)1 << s2); int64_t l = std::min(s1, s2); int64_t u = std::max(s1, s2); bit_t fermimask = (((bit_t)1 << (u - l - 1)) - 1) << (l + 1); - coeff_t t = bond.coupling(); auto term_action = [&](bit_t spins) -> std::pair { bool fermi = bits::popcnt(spins & fermimask) & 1; @@ -42,8 +45,8 @@ void apply_hopping(Bond const &bond, Basis &&basis, Fill &&fill) { auto non_zero_term = [&flipmask](bit_t const &ups) -> bool { return bits::popcnt(ups & flipmask) & 1; }; - tj::generic_term_ups( - basis, basis, non_zero_term, term_action, fill); + tj::generic_term_ups(basis, basis, non_zero_term, + term_action, fill); } else if (type == "HOPDN") { auto non_zero_term_ups = [&flipmask](bit_t const &ups) -> bool { return (ups & flipmask) == 0; @@ -52,8 +55,7 @@ void apply_hopping(Bond const &bond, Basis &&basis, Fill &&fill) { return bits::popcnt(dns & flipmask) & 1; }; tj::generic_term_dns( - basis, basis, non_zero_term_ups, non_zero_term_dns, term_action, - fill); + basis, basis, non_zero_term_ups, non_zero_term_dns, term_action, fill); } } diff --git a/xdiag/blocks/tj/terms/apply_ising.hpp b/xdiag/blocks/tj/terms/apply_ising.hpp index d063e348..6dd72b81 100644 --- a/xdiag/blocks/tj/terms/apply_ising.hpp +++ b/xdiag/blocks/tj/terms/apply_ising.hpp @@ -17,7 +17,10 @@ void apply_ising(Bond const &bond, Basis &&basis, Fill &&fill) { std::string type = bond.type(); assert((type == "ISING") || (type == "TJISING")); - coeff_t J = bond.coupling(); + Coupling cpl = bond.coupling(); + assert(cpl.isexplicit() && !cpl.ismatrix()); + coeff_t J = cpl.as(); + int64_t s1 = bond[0]; int64_t s2 = bond[1]; bit_t s1_mask = (bit_t)1 << s1; diff --git a/xdiag/blocks/tj/terms/apply_number.hpp b/xdiag/blocks/tj/terms/apply_number.hpp index db9cf39b..2f0f8d12 100644 --- a/xdiag/blocks/tj/terms/apply_number.hpp +++ b/xdiag/blocks/tj/terms/apply_number.hpp @@ -16,7 +16,10 @@ void apply_number(Bond const &bond, Basis &&basis, Fill &&fill) { std::string type = bond.type(); assert((type == "NUMBERUP") || (type == "NUMBERDN")); - coeff_t mu = bond.coupling(); + Coupling cpl = bond.coupling(); + assert(cpl.isexplicit() && !cpl.ismatrix()); + coeff_t mu = cpl.as(); + int64_t s = bond[0]; bit_t mask = (bit_t)1 << s; diff --git a/xdiag/blocks/tj/terms/apply_raise_lower.hpp b/xdiag/blocks/tj/terms/apply_raise_lower.hpp index 5d6c4fbd..c9899172 100644 --- a/xdiag/blocks/tj/terms/apply_raise_lower.hpp +++ b/xdiag/blocks/tj/terms/apply_raise_lower.hpp @@ -23,7 +23,10 @@ void apply_raise_lower(Bond const &bond, BasisIn &&basis_in, int64_t s = bond[0]; bit_t site_mask = (bit_t)1 << s; bit_t fermi_mask = site_mask - 1; - coeff_t c = bond.coupling(); + + Coupling cpl = bond.coupling(); + assert(cpl.isexplicit() && !cpl.ismatrix()); + coeff_t c = cpl.as(); // Raising operators if ((type == "CDAGUP") || (type == "CDAGDN")) { diff --git a/xdiag/blocks/tj/terms/apply_terms.hpp b/xdiag/blocks/tj/terms/apply_terms.hpp index cf31fea8..79f91e02 100644 --- a/xdiag/blocks/tj/terms/apply_terms.hpp +++ b/xdiag/blocks/tj/terms/apply_terms.hpp @@ -13,9 +13,10 @@ namespace xdiag::tj { template void apply_terms(BondList const &bonds, BasisIn const &basis_in, - BasisOut const &basis_out, Fill &fill) { + BasisOut const &basis_out, Fill &fill, double zero_precision) { - BondList bonds_compiled = tj::compile(bonds); + BondList bonds_compiled = + tj::compile(bonds, basis_in.n_sites(), zero_precision); for (auto bond : bonds_compiled) { std::string type = bond.type(); if ((type == "ISING") || (type == "TJISING")) { diff --git a/xdiag/blocks/tj/tj_apply.cpp b/xdiag/blocks/tj/tj_apply.cpp index 8ddd0e09..e338a7bb 100644 --- a/xdiag/blocks/tj/tj_apply.cpp +++ b/xdiag/blocks/tj/tj_apply.cpp @@ -8,24 +8,24 @@ namespace xdiag { template void apply(BondList const &bonds, tJ const &block_in, arma::Col const &vec_in, tJ const &block_out, - arma::Col &vec_out) try { + arma::Col &vec_out, double zero_precision) try { vec_out.zeros(); - BondList bondsc = tj::compile(bonds, 1e-12); + BondList bondsc = tj::compile(bonds, block_in.n_sites(), zero_precision); auto fill = [&](int64_t idx_in, int64_t idx_out, coeff_t val) { return fill_apply(vec_in, vec_out, idx_in, idx_out, val); }; - tj::dispatch(bondsc, block_in, block_out, fill); + tj::dispatch(bondsc, block_in, block_out, fill, zero_precision); } catch (Error const &e) { XDIAG_RETHROW(e); } template void apply(BondList const &, tJ const &, arma::Col const &, tJ const &, - arma::Col &); + arma::Col &, double); template void apply(BondList const &, tJ const &, arma::Col const &, tJ const &, - arma::Col &); + arma::Col &, double); } // namespace xdiag diff --git a/xdiag/blocks/tj/tj_apply.hpp b/xdiag/blocks/tj/tj_apply.hpp index 2aff3a9d..9671fbe9 100644 --- a/xdiag/blocks/tj/tj_apply.hpp +++ b/xdiag/blocks/tj/tj_apply.hpp @@ -10,6 +10,6 @@ namespace xdiag { template void apply(BondList const &bonds, tJ const &block_in, arma::Col const &vec_in, tJ const &block_out, - arma::Col &vec_out); + arma::Col &vec_out, double zero_precision = 1e-12); } // namespace xdiag diff --git a/xdiag/blocks/tj/tj_matrix.cpp b/xdiag/blocks/tj/tj_matrix.cpp index 3d398fa7..dded36e1 100644 --- a/xdiag/blocks/tj/tj_matrix.cpp +++ b/xdiag/blocks/tj/tj_matrix.cpp @@ -8,39 +8,33 @@ namespace xdiag { template arma::Mat matrix_gen(BondList const &bonds, tJ const &block_in, - tJ const &block_out) try { - BondList bondsc = tj::compile(bonds, 1e-12); - + tJ const &block_out, double zero_precision) try { int64_t m = block_out.size(); int64_t n = block_in.dim(); arma::Mat mat(m, n, arma::fill::zeros); - matrix_gen(mat.memptr(), bonds, block_in, block_out); + matrix_gen(mat.memptr(), bonds, block_in, block_out, zero_precision); return mat; } catch (Error const &e) { XDIAG_RETHROW(e); - return arma::Mat(); } -arma::mat matrix(BondList const &bonds, tJ const &block_in, - tJ const &block_out) try { - return matrix_gen(bonds, block_in, block_out); +arma::mat matrix(BondList const &bonds, tJ const &block_in, tJ const &block_out, + double zero_precision) try { + return matrix_gen(bonds, block_in, block_out, zero_precision); } catch (Error const &e) { XDIAG_RETHROW(e); - return arma::mat(); } arma::cx_mat matrixC(BondList const &bonds, tJ const &block_in, - tJ const &block_out) try { - return matrix_gen(bonds, block_in, block_out); + tJ const &block_out, double zero_precision) try { + return matrix_gen(bonds, block_in, block_out, zero_precision); } catch (Error const &e) { XDIAG_RETHROW(e); - return arma::cx_mat(); } template void matrix_gen(coeff_t *mat, BondList const &bonds, tJ const &block_in, - tJ const &block_out) try { - BondList bondsc = tj::compile(bonds, 1e-12); + tJ const &block_out, double zero_precision) try { int64_t m = block_out.size(); int64_t n = block_in.size(); @@ -50,22 +44,22 @@ void matrix_gen(coeff_t *mat, BondList const &bonds, tJ const &block_in, return fill_matrix(mat, idx_in, idx_out, m, val); }; - tj::dispatch(bondsc, block_in, block_out, fill); + tj::dispatch(bonds, block_in, block_out, fill, zero_precision); } catch (Error const &e) { XDIAG_RETHROW(e); } void matrix(double *mat, BondList const &bonds, tJ const &block_in, - tJ const &block_out) try { - return matrix_gen(mat, bonds, block_in, block_out); + tJ const &block_out, double zero_precision) try { + return matrix_gen(mat, bonds, block_in, block_out, zero_precision); } catch (Error const &e) { XDIAG_RETHROW(e); } void matrixC(complex *mat, BondList const &bonds, tJ const &block_in, - tJ const &block_out) try { - return matrix_gen(mat, bonds, block_in, block_out); + tJ const &block_out, double zero_precision) try { + return matrix_gen(mat, bonds, block_in, block_out, zero_precision); } catch (Error const &e) { XDIAG_RETHROW(e); } diff --git a/xdiag/blocks/tj/tj_matrix.hpp b/xdiag/blocks/tj/tj_matrix.hpp index a84d6335..39b78c50 100644 --- a/xdiag/blocks/tj/tj_matrix.hpp +++ b/xdiag/blocks/tj/tj_matrix.hpp @@ -7,17 +7,17 @@ namespace xdiag { -arma::mat matrix(BondList const &bonds, tJ const &block_in, - tJ const &block_out); +arma::mat matrix(BondList const &bonds, tJ const &block_in, tJ const &block_out, + double zero_precision = 1e-12); arma::cx_mat matrixC(BondList const &bonds, tJ const &block_in, - tJ const &block_out); + tJ const &block_out, double zero_precision = 1e-12); // Developer functions void matrix(double *mat, BondList const &bonds, tJ const &block_in, - tJ const &block_out); + tJ const &block_out, double zero_precision = 1e-12); void matrixC(complex *mat, BondList const &bonds, tJ const &block_in, - tJ const &block_out); + tJ const &block_out, double zero_precision = 1e-12); } // namespace xdiag diff --git a/xdiag/blocks/tj_distributed/tj_distributed_apply.cpp b/xdiag/blocks/tj_distributed/tj_distributed_apply.cpp index 52d7e31d..5622a2f1 100644 --- a/xdiag/blocks/tj_distributed/tj_distributed_apply.cpp +++ b/xdiag/blocks/tj_distributed/tj_distributed_apply.cpp @@ -9,9 +9,9 @@ namespace xdiag { template void apply(BondList const &bonds, tJDistributed const &block_in, arma::Col const &vec_in, tJDistributed const &block_out, - arma::Col &vec_out) try { + arma::Col &vec_out, double zero_precision) try { vec_out.zeros(); - BondList bondsc = tj::compile(bonds, 1e-12); + BondList bondsc = tj::compile(bonds, block_in.n_sites(), zero_precision); tj_distributed::dispatch(bondsc, block_in, vec_in, block_out, vec_out); } catch (Error const &e) { @@ -20,10 +20,10 @@ void apply(BondList const &bonds, tJDistributed const &block_in, template void apply(BondList const &, tJDistributed const &, arma::Col const &, tJDistributed const &, - arma::Col &); + arma::Col &, double); template void apply(BondList const &, tJDistributed const &, arma::Col const &, tJDistributed const &, - arma::Col &); + arma::Col &, double); } // namespace xdiag diff --git a/xdiag/io/toml/toml_conversion.cpp b/xdiag/io/toml/toml_conversion.cpp index 4be16544..af42ab5b 100644 --- a/xdiag/io/toml/toml_conversion.cpp +++ b/xdiag/io/toml/toml_conversion.cpp @@ -16,7 +16,6 @@ template T get_toml_value(toml_const_node_view const &node) try { } } catch (Error const &e) { XDIAG_RETHROW(e); - return T(); } template <> @@ -34,7 +33,6 @@ complex get_toml_value(toml_const_node_view const &node) try { } } catch (Error const &e) { XDIAG_RETHROW(e); - return 0.; } template bool get_toml_value(toml_const_node_view const &); @@ -58,7 +56,6 @@ template T get_toml_value(toml_node_view const &node) try { } } catch (Error const &e) { XDIAG_RETHROW(e); - return T(); } template <> complex get_toml_value(toml_node_view const &node) try { @@ -71,11 +68,9 @@ template <> complex get_toml_value(toml_node_view const &node) try { return complex{real_imag[0], real_imag[1]}; } else { XDIAG_THROW("Error parsing toml: failed trying to parse complex number!"); - return complex(); } } catch (Error const &e) { XDIAG_RETHROW(e); - return 0.; } template bool get_toml_value(toml_node_view const &); @@ -99,7 +94,6 @@ template T get_toml_value(toml::node const &node) try { } } catch (Error const &e) { XDIAG_RETHROW(e); - return T(); } template <> complex get_toml_value(toml::node const &node) try { @@ -115,7 +109,6 @@ template <> complex get_toml_value(toml::node const &node) try { } } catch (Error const &e) { XDIAG_RETHROW(e); - return 0.; } template bool get_toml_value(toml::node const &); template int8_t get_toml_value(toml::node const &); @@ -138,7 +131,6 @@ toml::array get_toml_array(toml_const_node_view const &node) try { } } catch (Error const &e) { XDIAG_RETHROW(e); - return toml::array(); } toml::array get_toml_array(toml_node_view const &node) try { @@ -147,11 +139,9 @@ toml::array get_toml_array(toml_node_view const &node) try { return *arr; } else { XDIAG_THROW("Error parsing toml: node does not contain an array!"); - return toml::array(); } } catch (Error const &e) { XDIAG_RETHROW(e); - return toml::array(); } toml::array get_toml_array(toml::node const &node) try { @@ -163,7 +153,6 @@ toml::array get_toml_array(toml::node const &node) try { } } catch (Error const &e) { XDIAG_RETHROW(e); - return toml::array(); } template @@ -182,7 +171,6 @@ std::array toml_array_to_std_array(toml::array const &toml_array) try { return array; } catch (Error const &e) { XDIAG_RETHROW(e); - return std::array(); } template std::array @@ -198,7 +186,6 @@ std::vector toml_array_to_std_vector(toml::array const &toml_array) try { return vector; } catch (Error const &e) { XDIAG_RETHROW(e); - return std::vector(); } template std::vector @@ -234,7 +221,6 @@ arma::Col toml_array_to_arma_vector(toml::array const &toml_array) try { return vector; } catch (Error const &e) { XDIAG_RETHROW(e); - return arma::Col(); } template arma::Col @@ -283,7 +269,6 @@ arma::Mat toml_array_to_arma_matrix(toml::array const &toml_array) try { } } catch (Error const &e) { XDIAG_RETHROW(e); - return arma::Mat(); } template arma::Mat @@ -304,7 +289,6 @@ toml::array std_vector_to_toml_array(std::vector const &value) try { return arr; } catch (Error const &e) { XDIAG_RETHROW(e); - return toml::array(); } template <> @@ -316,7 +300,6 @@ toml::array std_vector_to_toml_array(std::vector const &value) try { return arr; } catch (Error const &e) { XDIAG_RETHROW(e); - return toml::array(); } template <> @@ -328,7 +311,6 @@ toml::array std_vector_to_toml_array(std::vector const &value) try { return arr; } catch (Error const &e) { XDIAG_RETHROW(e); - return toml::array(); } template toml::array std_vector_to_toml_array(std::vector const &); @@ -350,7 +332,6 @@ toml::array arma_vector_to_toml_array(arma::Col const &value) try { return arr; } catch (Error const &e) { XDIAG_RETHROW(e); - return toml::array(); } template <> @@ -362,7 +343,6 @@ toml::array arma_vector_to_toml_array(arma::Col const &value) try { return arr; } catch (Error const &e) { XDIAG_RETHROW(e); - return toml::array(); } template <> @@ -374,7 +354,6 @@ toml::array arma_vector_to_toml_array(arma::Col const &value) try { return arr; } catch (Error const &e) { XDIAG_RETHROW(e); - return toml::array(); } template toml::array arma_vector_to_toml_array(arma::Col const &); @@ -393,7 +372,6 @@ toml::array arma_matrix_to_toml_array(arma::Mat const &mat) try { return arr; } catch (Error const &e) { XDIAG_RETHROW(e); - return toml::array(); } template <> @@ -409,7 +387,6 @@ toml::array arma_matrix_to_toml_array(arma::Mat const &mat) try { return arr; } catch (Error const &e) { XDIAG_RETHROW(e); - return toml::array(); } template <> @@ -438,36 +415,35 @@ toml::array bond_to_toml_array(Bond const &bond) try { array.push_back(bond.type()); // Coupling - if (bond.coupling_is()) { - array.push_back(bond.coupling()); - } else if (bond.coupling_is()) { - array.push_back(bond.coupling()); - } else if (bond.coupling_is()) { - complex cpl = bond.coupling(); + if (bond.coupling().is()) { + array.push_back(bond.coupling().as()); + } else if (bond.coupling().is()) { + array.push_back(bond.coupling().as()); + } else if (bond.coupling().is()) { + complex cpl = bond.coupling().as(); array.push_back(toml::array{cpl.real(), cpl.imag()}); - } else if (bond.coupling_is()) { - arma::mat mat = bond.coupling(); + } else if (bond.coupling().is()) { + arma::mat mat = bond.coupling().as(); array.push_back(arma_matrix_to_toml_array(mat)); - } else if (bond.coupling_is()) { - arma::cx_mat mat = bond.coupling(); + } else if (bond.coupling().is()) { + arma::cx_mat mat = bond.coupling().as(); array.push_back(arma_matrix_to_toml_array(mat)); } else { XDIAG_THROW("Unable to convert Bond to toml array "); } // Sites - if (bond.size() == 1) { - array.push_back(bond[0]); - } else { - array.push_back(std_vector_to_toml_array(bond.sites())); + for (int i = 0; i < bond.size(); ++i) { + array.push_back(bond[i]); } + return array; } catch (Error const &e) { XDIAG_RETHROW(e); } Bond toml_array_to_bond(toml::array const &array) try { - if (array.size() != 3) { + if (array.size() < 3) { XDIAG_THROW( "Error parsing toml to xdiag::Bond: toml array must have exactly three " "entries.\n" @@ -490,17 +466,16 @@ Bond toml_array_to_bond(toml::array const &array) try { // then get the sites std::vector sites; - auto site_node = array[2].value(); - auto sites_node = array[2].as_array(); - if (site_node) { - sites = {*site_node}; - } else if (sites_node) { - sites = toml_array_to_std_vector(*sites_node); - } else { - XDIAG_THROW("Error parsing toml to xdiag::Bond: third entry must either " - "be a single integer or a 1D array of integers"); - } + for (int i = 2; i < array.size(); ++i) { + auto site_node = array[i].value(); + if (site_node) { + sites.push_back(*site_node); + } else { + XDIAG_THROW("Error parsing toml to xdiag::Bond: third and onward entries " + "must be integers"); + } + } // then get the coupling auto coupling_node_string = array[1].value(); auto coupling_node_double = array[1].value(); @@ -534,14 +509,27 @@ Bond toml_array_to_bond(toml::array const &array) try { "the second entry defining the coupling. However, it is " "found to be empty"); } - auto real_array = mat_array.value(); - auto cplx_array = mat_array.as_array(); - if (real_array) { - arma::mat coupling = toml_array_to_arma_matrix(mat_array); - return Bond(type, coupling, sites); - } else if (cplx_array) { - arma::cx_mat coupling = toml_array_to_arma_matrix(mat_array); - return Bond(type, coupling, sites); + auto mat2_array = mat_array[0].as_array(); + if (mat2_array) { + if (mat2_array->size() == 0) { + XDIAG_THROW( + "Error parsing toml to xdiag::Bond: an array was handed to " + "the second entry defining the coupling. However, it is " + "found to be empty"); + } + auto real_array = (*mat2_array)[0].value(); + auto cplx_array = (*mat2_array)[0].as_array(); + + if (real_array) { + arma::mat coupling = toml_array_to_arma_matrix(mat_array); + return Bond(type, coupling, sites); + } else if (cplx_array) { + arma::cx_mat coupling = toml_array_to_arma_matrix(mat_array); + return Bond(type, coupling, sites); + } else { + XDIAG_THROW("Error parsing toml to xdiag::Bond: Invalid form of the " + "coupling."); + } } else { XDIAG_THROW("Error parsing toml to xdiag::Bond: Invalid form of the " "coupling."); @@ -550,13 +538,15 @@ Bond toml_array_to_bond(toml::array const &array) try { XDIAG_THROW("Error parsing toml to xdiag::Bond: Invalid form of the " "coupling."); } + } else { + XDIAG_THROW("Error parsing toml to xdiag::Bond: Invalid form of the " + "coupling."); } } catch (Error const &e) { XDIAG_RETHROW(e); - return Bond(); } -BondList toml_array_to_bond_list(toml::array const &array) { +BondList toml_array_to_bond_list(toml::array const &array) try { BondList bonds; for (std::size_t i = 0; i < array.size(); ++i) { auto bond_array = array[i].as_array(); @@ -570,9 +560,11 @@ BondList toml_array_to_bond_list(toml::array const &array) { } } return bonds; +} catch (Error const &e) { + XDIAG_RETHROW(e); } -BondList toml_table_to_bond_list(toml::table const &table) { +BondList toml_table_to_bond_list(toml::table const &table) try { BondList bonds; auto interactions_opt = table["Interactions"].as_array(); if (interactions_opt) { @@ -584,20 +576,58 @@ BondList toml_table_to_bond_list(toml::table const &table) { toml::table couplings = *couplings_opt; for (auto [key_node, value_node] : couplings) { std::string key(key_node.str()); - auto value = get_toml_value(value_node); - bonds[key] = value; + + auto val = value_node.value(); + if (val) { + auto value = get_toml_value(value_node); + bonds[key] = value; + } else { + if (value_node.as_array()->size() == 0) { + XDIAG_THROW("Invalid format of coupling in BondList") + } + + auto val = *value_node.as_array(); + auto val_real = val[0].value(); + auto val_array = val[0].as_array(); + if (val_real) { + auto value = get_toml_value(value_node); + bonds[key] = value; + } + if (val_array) { + auto val = *(val_array->as_array()); + + auto val2_real = val[0].value(); + auto val2_arr = val[0].as_array(); + if (val2_real) { + auto value = + toml_array_to_arma_matrix(*value_node.as_array()); + bonds[key] = value; + } else if (val2_arr) { + auto value = + toml_array_to_arma_matrix(*value_node.as_array()); + bonds[key] = value; + + } else { + XDIAG_THROW("Invalid format of coupling in BondList") + } + } + } } } return bonds; +} catch (Error const &e) { + XDIAG_RETHROW(e); } -toml::array bond_list_to_toml_array(BondList const &bonds) { +toml::array bond_list_to_toml_array(BondList const &bonds) try { toml::array array; for (auto &&bond : bonds) { array.push_back(bond_to_toml_array(bond)); } return array; +} catch (Error const &e) { + XDIAG_RETHROW(e); } toml::table bond_list_to_toml_table(BondList const &bonds) try { @@ -607,20 +637,20 @@ toml::table bond_list_to_toml_table(BondList const &bonds) try { if (bonds.couplings().size() > 0) { toml::table couplings; for (auto name : bonds.couplings()) { - coupling_t cpl = bonds[name]; - if (std::holds_alternative(cpl)) { - couplings.insert_or_assign(name, std::get(cpl)); - } else if (std::holds_alternative(cpl)) { - couplings.insert_or_assign(name, std::get(cpl)); - } else if (std::holds_alternative(cpl)) { - complex c = std::get(cpl); + Coupling cpl = bonds[name]; + if (cpl.is()) { + couplings.insert_or_assign(name, cpl.as()); + } else if (cpl.is()) { + couplings.insert_or_assign(name, cpl.as()); + } else if (cpl.is()) { + complex c = cpl.as(); couplings.insert_or_assign(name, toml::array{c.real(), c.imag()}); - } else if (std::holds_alternative(cpl)) { + } else if (cpl.is()) { couplings.insert_or_assign( - name, arma_matrix_to_toml_array(std::get(cpl))); - } else if (std::holds_alternative(cpl)) { + name, arma_matrix_to_toml_array(cpl.as())); + } else if (cpl.is()) { couplings.insert_or_assign( - name, arma_matrix_to_toml_array(std::get(cpl))); + name, arma_matrix_to_toml_array(cpl.as())); } else { XDIAG_THROW("Unable to convert BondList to toml"); } @@ -630,7 +660,6 @@ toml::table bond_list_to_toml_table(BondList const &bonds) try { return table; } catch (Error const &e) { XDIAG_RETHROW(e); - return toml::table(); } } // namespace xdiag::io diff --git a/xdiag/operators/bond.cpp b/xdiag/operators/bond.cpp index ef9370be..9281c0ae 100644 --- a/xdiag/operators/bond.cpp +++ b/xdiag/operators/bond.cpp @@ -8,120 +8,130 @@ namespace xdiag { -Bond::Bond(std::string type, coupling_t coupling, +Bond::Bond(std::string type, Coupling coupling, std::vector const &sites) try : type_(type), coupling_(coupling), sites_(sites) { } catch (Error const &e) { XDIAG_RETHROW(e); } -Bond::Bond(std::string type, coupling_t coupling, int64_t site) try +Bond::Bond(std::string type, Coupling coupling, int64_t site) try : type_(type), coupling_(coupling), sites_({site}) { } catch (Error const &e) { XDIAG_RETHROW(e); } -std::string Bond::type() const { return type_; } +Bond::Bond(std::string type, const char *coupling, + std::vector const &sites) try + : Bond(type, Coupling(std::string(coupling)), sites) { +} catch (Error const &e) { + XDIAG_RETHROW(e); +} +Bond::Bond(std::string type, const char *coupling, int64_t site) try + : Bond(type, coupling, std::vector({site})) { +} catch (Error const &e) { + XDIAG_RETHROW(e); +} -template bool Bond::coupling_is() const { - return std::holds_alternative(coupling_); -}; -template bool Bond::coupling_is() const; -template bool Bond::coupling_is() const; -template bool Bond::coupling_is() const; -template bool Bond::coupling_is() const; -template bool Bond::coupling_is() const; - -template T Bond::coupling() const try { - return std::get(coupling_); -} catch (std::bad_variant_access const &ex) { - XDIAG_THROW("Invalid type given when retrieving coupling"); - return T(); -} -template std::string Bond::coupling() const; -template double Bond::coupling() const; -template arma::mat Bond::coupling() const; - -template <> complex Bond::coupling() const try { - if (coupling_is()) { - return complex(coupling(), 0.0); - } else { - return std::get(coupling_); - } -} catch (std::bad_variant_access const &ex) { - XDIAG_THROW("Invalid type given when retrieving coupling"); - return complex(); +Bond::Bond(std::string type, std::string coupling, + std::vector const &sites) try + : Bond(type, Coupling(coupling), sites) { +} catch (Error const &e) { + XDIAG_RETHROW(e); +} +Bond::Bond(std::string type, std::string coupling, int64_t site) try + : Bond(type, coupling, std::vector({site})) { +} catch (Error const &e) { + XDIAG_RETHROW(e); +} +Bond::Bond(std::string type, double coupling, + std::vector const &sites) try + : Bond(type, Coupling(coupling), sites) { +} catch (Error const &e) { + XDIAG_RETHROW(e); +} +Bond::Bond(std::string type, double coupling, int64_t site) try + : Bond(type, coupling, std::vector({site})) { +} catch (Error const &e) { + XDIAG_RETHROW(e); +} +Bond::Bond(std::string type, complex coupling, + std::vector const &sites) try + : Bond(type, Coupling(coupling), sites) { +} catch (Error const &e) { + XDIAG_RETHROW(e); +} +Bond::Bond(std::string type, complex coupling, int64_t site) try + : Bond(type, coupling, std::vector({site})) { +} catch (Error const &e) { + XDIAG_RETHROW(e); } -template <> arma::cx_mat Bond::coupling() const try { - if (coupling_is()) { - arma::mat mat_r = coupling(); - arma::mat mat_i = arma::mat(mat_r.n_rows, mat_r.n_cols, arma::fill::zeros); - return arma::cx_mat(mat_r, mat_i); - } else { - std::get(coupling_); - } -} catch (std::bad_variant_access const &ex) { - XDIAG_THROW("Invalid type given when retrieving coupling"); - return arma::cx_mat(); -} -coupling_t Bond::coupling() const { return coupling_; } -std::string Bond::coupling_type_string() const try { - if (std::holds_alternative(coupling_)) { - return "std::string"; - } else if (std::holds_alternative(coupling_)) { - return "double"; - } else if (std::holds_alternative(coupling_)) { - return "complex"; - } else if (std::holds_alternative(coupling_)) { - return "arma::mat"; - } else if (std::holds_alternative(coupling_)) { - return "arma::cx_mat"; - } else { - XDIAG_THROW("Invalid coupling type found"); - } +Bond::Bond(std::string type, arma::mat const &coupling, + std::vector const &sites) try + : Bond(type, Coupling(coupling), sites) { +} catch (Error const &e) { + XDIAG_RETHROW(e); +} +Bond::Bond(std::string type, arma::mat const &coupling, int64_t site) try + : Bond(type, coupling, std::vector({site})) { } catch (Error const &e) { XDIAG_RETHROW(e); - return ""; } +Bond::Bond(std::string type, arma::cx_mat const &coupling, + std::vector const &sites) try + : Bond(type, Coupling(coupling), sites) { +} catch (Error const &e) { + XDIAG_RETHROW(e); +} +Bond::Bond(std::string type, arma::cx_mat const &coupling, int64_t site) try + : Bond(type, coupling, std::vector({site})) { +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +std::string Bond::type() const { return type_; } + +Coupling const &Bond::coupling() const { return coupling_; } + int64_t Bond::size() const { return sites_.size(); } int64_t Bond::operator[](int64_t idx) const try { return sites_.at(idx); } catch (std::out_of_range const &exc) { XDIAG_THROW("Site index out of range for Bond"); - return 0; } std::vector const &Bond::sites() const { return sites_; } bool Bond::isexplicit() const try { - return !coupling_is(); + return coupling_.isexplicit(); } catch (Error const &e) { XDIAG_RETHROW(e); - return false; } bool Bond::isreal() const try { - if (coupling_is()) { + if (coupling_.is()) { XDIAG_THROW("Cannot determine whether coupling is real, since coupling is " "set to a std::string value"); + } else if (type_ == "SCALARCHIRALITY") { + return false; } else { - return coupling_is() || coupling_is(); + return coupling_.isreal(); } } catch (Error const &e) { XDIAG_RETHROW(e); - return false; } -bool Bond::ismatrix() const { - return coupling_is() || coupling_is(); -} +bool Bond::ismatrix() const { return coupling_.ismatrix(); } -bool Bond::operator==(const Bond &rhs) const { +bool Bond::operator==(Bond const &rhs) const try { return (type_ == rhs.type_) && (coupling_ == rhs.coupling_) && (sites_ == rhs.sites_); +} catch (Error const &e) { + XDIAG_RETHROW(e); } -bool Bond::operator!=(const Bond &rhs) const { return !operator==(rhs); } + +bool Bond::operator!=(Bond const &rhs) const { return !operator==(rhs); } bool sites_disjoint(Bond const &bond) { auto const &sites = bond.sites(); @@ -142,28 +152,14 @@ std::vector common_sites(Bond const &b1, Bond const &b2) { std::ostream &operator<<(std::ostream &out, const Bond &bond) try { - if (bond.coupling_is()) { - out << bond.type() << " " << bond.coupling() << " "; - } else if (bond.coupling_is()) { - out << bond.type() << " " << bond.coupling() << " "; - } else if (bond.coupling_is()) { - out << bond.type() << " " << bond.coupling() << " "; - } else if (bond.coupling_is()) { - out << bond.type() << " " << bond.coupling() << " "; - } else if (bond.coupling_is()) { - out << bond.type() << " " << bond.coupling() << " "; - } else { - XDIAG_THROW("Unable to print Bond"); - } - + out << bond.type() << " " << bond.coupling() << " "; for (int64_t site : bond.sites()) { out << site << " "; } + return out; } catch (Error const &e) { XDIAG_RETHROW(e); - return out; } - } // namespace xdiag diff --git a/xdiag/operators/bond.hpp b/xdiag/operators/bond.hpp index 87562be9..68d284f0 100644 --- a/xdiag/operators/bond.hpp +++ b/xdiag/operators/bond.hpp @@ -1,45 +1,53 @@ #pragma once -#include -#include #include - #include #include +#include namespace xdiag { -using coupling_t = - std::variant; - class Bond { public: // Constructors with type name Bond() = default; - Bond(std::string type, coupling_t coupling, + Bond(std::string type, Coupling coupling, std::vector const &sites); + Bond(std::string type, Coupling coupling, int64_t site); + + Bond(std::string type, const char *coupling, + std::vector const &sites); + Bond(std::string type, const char *coupling, int64_t site); + Bond(std::string type, std::string coupling, + std::vector const &sites); + Bond(std::string type, std::string coupling, int64_t site); + Bond(std::string type, double coupling, std::vector const &sites); + Bond(std::string type, double coupling, int64_t site); + Bond(std::string type, complex coupling, std::vector const &sites); + Bond(std::string type, complex coupling, int64_t site); + Bond(std::string type, arma::mat const &coupling, std::vector const &sites); - Bond(std::string type, coupling_t coupling, int64_t site); + Bond(std::string type, arma::mat const &coupling, int64_t site); + Bond(std::string type, arma::cx_mat const &coupling, + std::vector const &sites); + Bond(std::string type, arma::cx_mat const &coupling, int64_t site); std::string type() const; - template bool coupling_is() const; - template T coupling() const; - coupling_t coupling() const; - std::string coupling_type_string() const; + Coupling const &coupling() const; int64_t size() const; int64_t operator[](int64_t idx) const; std::vector const &sites() const; - bool isexplicit() const; bool isreal() const; bool ismatrix() const; + bool isexplicit() const; bool operator==(const Bond &rhs) const; bool operator!=(const Bond &rhs) const; private: std::string type_; - coupling_t coupling_; + Coupling coupling_; std::vector sites_; }; @@ -47,5 +55,4 @@ bool sites_disjoint(Bond const &bond); std::vector common_sites(Bond const &b1, Bond const &b2); std::ostream &operator<<(std::ostream &out, Bond const &bond); - } // namespace xdiag diff --git a/xdiag/operators/bondlist.cpp b/xdiag/operators/bondlist.cpp index 8a2b15fb..443e44a5 100644 --- a/xdiag/operators/bondlist.cpp +++ b/xdiag/operators/bondlist.cpp @@ -12,15 +12,16 @@ namespace xdiag { BondList::BondList(std::vector const &bonds) : bonds_(bonds) {} int64_t BondList::size() const { return bonds_.size(); } -coupling_t &BondList::operator[](std::string name) { return coupling_[name]; } -coupling_t const &BondList::operator[](std::string name) const { - return coupling_.at(name); +Coupling &BondList::operator[](std::string name) { return couplings_[name]; } +Coupling const &BondList::operator[](std::string name) const { + return couplings_.at(name); } -bool BondList::coupling_defined(std::string name) { - return couplings_.find(name) == couplings_.end(); + +bool BondList::defined(std::string name) const { + return couplings_.find(name) != couplings_.end(); } -std::vector const &couplings() const { +std::vector BondList::couplings() const { std::vector names; for (auto const &cpl : couplings_) { names.push_back(cpl.first); @@ -28,14 +29,19 @@ std::vector const &couplings() const { return names; } -bool BondList::isreal() const { - std::all_of(bonds_.begin(), bonds_.end(), - [](Bond const &b) { return b.isreal(); }); +bool BondList::isreal() const try { + BondList bonds_explicit = make_explicit(*this); + return std::all_of(bonds_explicit.begin(), bonds_explicit.end(), + [](Bond const &b) { return b.isreal(); }); +} catch (Error const &e) { + XDIAG_RETHROW(e); } -bool BondList::isexplicit() const { - std::all_of(bonds_.begin(), bonds_.end(), - [](Bond const &b) { return b.isexplicit(); }); +bool BondList::isexplicit() const try { + return std::all_of(bonds_.begin(), bonds_.end(), + [](Bond const &b) { return b.isexplicit(); }); +} catch (Error const &e) { + XDIAG_RETHROW(e); } void BondList::operator+=(Bond const &bond) { bonds_.push_back(bond); } @@ -48,7 +54,7 @@ void BondList::operator+=(BondList const &bonds) try { // Add possible couplings for (auto it = bonds.couplings_.begin(); it != bonds.couplings_.end(); it++) { std::string name = it->first; - coupling_t cpl = it->second; + Coupling cpl = it->second; // if name does not exist yet, add coupling if (couplings_.find(name) == couplings_.end()) { @@ -67,21 +73,35 @@ void BondList::operator+=(BondList const &bonds) try { XDIAG_RETHROW(e); } -BondList operator+(Bond const &bond) const { +BondList BondList::operator+(Bond const &bond) const { BondList new_bonds = *this; new_bonds += bond; return new_bonds; } -BondList operator+(BondList const &bonds) const try { +BondList BondList::operator+(BondList const &bonds) const try { BondList new_bonds = *this; new_bonds += bonds; + return new_bonds; } catch (Error const &e) { XDIAG_RETHROW(e); } bool BondList::operator==(BondList const &other) const { - return (bonds_ == other.bonds_) && (couplings_ == other.couplings_); + if (bonds_ != other.bonds_) { + return false; + } + auto c1s = couplings(); + auto c2s = other.couplings(); + if (c1s != c2s) { + return false; + } + for (std::string c : c1s) { + if (couplings_.at(c) != other.couplings_.at(c)) { + return false; + } + } + return true; } bool BondList::operator!=(BondList const &other) const { return !operator==(other); @@ -98,18 +118,18 @@ BondList make_explicit(BondList const &bonds) try { BondList explicit_bonds; for (auto const &bond : bonds) { if (bond.isexplicit()) { - explicit_bonds.push_back(bond); + explicit_bonds += bond; } else { - std::string name = bond.coupling(); - if (bonds.coupling_defined(name)) { + std::string name = bond.coupling().as(); + if (bonds.defined(name)) { std::string type = bond.type(); - coupling_t cpl = bonds[name]; - if (std::holds_alternative(cpl)) { + Coupling cpl = bonds[name]; + if (!cpl.isexplicit()) { XDIAG_THROW(fmt::format("Unable to make BondList explicit: coupling " "\"{}\" is yet a string", name)); } - explicit_bonds.push_back(type, cpl, bond.sites()); + explicit_bonds += Bond(type, cpl, bond.sites()); } else { XDIAG_THROW(fmt::format( "Unable to make BondList explicit: coupling \"{}\" is not defined", @@ -117,10 +137,21 @@ BondList make_explicit(BondList const &bonds) try { } } } + return explicit_bonds; } catch (Error const &e) { XDIAG_RETHROW(e); } +BondList bonds_of_type(std::string type, BondList const &bonds) { + BondList new_bonds; + for (auto const &bond : bonds) { + if (bond.type() == type) { + new_bonds += bond; + } + } + return new_bonds; +} + BondList read_bondlist(std::string filename) { std::vector bonds; diff --git a/xdiag/operators/bondlist.hpp b/xdiag/operators/bondlist.hpp index f40c9ab2..138bf310 100644 --- a/xdiag/operators/bondlist.hpp +++ b/xdiag/operators/bondlist.hpp @@ -12,10 +12,10 @@ class BondList { explicit BondList(std::vector const &bonds); int64_t size() const; - coupling_t &operator[](std::string name); - coupling_t const &operator[](std::string name) const; - std::vector const &couplings() const; - bool coupling_defined(std::string name) const; + bool defined(std::string name) const; + Coupling &operator[](std::string name); + Coupling const &operator[](std::string name) const; + std::vector couplings() const; bool isreal() const; bool isexplicit() const; @@ -40,10 +40,13 @@ class BondList { private: std::vector bonds_; - std::map couplings_; + std::map couplings_; }; +BondList bonds_of_type(std::string type, BondList const &bonds); BondList make_explicit(BondList const &bonds); + +// legacy function, can be removed later BondList read_bondlist(std::string filename); } // namespace xdiag diff --git a/xdiag/operators/compiler.cpp b/xdiag/operators/compiler.cpp index e2980d0a..f8cbe2b6 100644 --- a/xdiag/operators/compiler.cpp +++ b/xdiag/operators/compiler.cpp @@ -3,6 +3,41 @@ #include namespace xdiag::operators { + + + +BondList clean_zeros(BondList const &bonds, double precision) try { + BondList clean_bonds; + for (auto const &bond : bonds) { + + // if bond is only defined by string, cannot be cleaned + if (!bond.isexplicit()) { + XDIAG_THROW("Cannot clean zero bond since is is not explicit. Maybe call " + "\"make_explicit\" first"); + } + + Coupling cpl = bond.coupling(); + if (cpl.is()) { + double cval = cpl.as(); + if (std::abs(cval) > precision) { + clean_bonds += bond; + } + } else if (cpl.is()) { + complex cval = cpl.as(); + if (std::abs(cval) > precision) { + clean_bonds += bond; + } + } + // coupling is matrix + else { + clean_bonds += bond; + } + } + return clean_bonds; +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + void check_bond(Bond const &bond, int64_t n_sites_total, int64_t n_sites_bond, bool disjoint, std::string type) { check_bond_in_range(bond, n_sites_total); @@ -54,7 +89,7 @@ void check_bond_has_correct_number_of_sites(Bond const &bond, int64_t ns) try { } void check_bond_has_disjoint_sites(Bond const &bond) try { - if (!bond.sites_disjoint()) { + if (!sites_disjoint(bond)) { XDIAG_THROW(std::string("bond of type \"") + bond.type() + std::string("\": sites are not disjoint as expected.")); } @@ -63,7 +98,7 @@ void check_bond_has_disjoint_sites(Bond const &bond) try { } void check_bond_coupling_has_type(Bond const &bond, std::string type) try { - std::string bond_type = bond.coupling_type_string(); + std::string bond_type = bond.coupling().type(); if (!(bond_type == type)) { XDIAG_THROW(fmt::format("Coupling of bond is expected to be of type \"{}\" " "but received a type \"{}\"", @@ -75,7 +110,7 @@ void check_bond_coupling_has_type(Bond const &bond, std::string type) try { void check_bond_coupling_has_type(Bond const &bond, std::string type1, std::string type2) try { - std::string bond_type = bond.coupling_type_string(); + std::string bond_type = bond.coupling().type(); if (!(bond_type == type1) && !(bond_type == type2)) { XDIAG_THROW(fmt::format("Coupling of bond is expected to be of type either " "\"{}\" or \"{}\" but received a type \"{}\"", @@ -85,174 +120,4 @@ void check_bond_coupling_has_type(Bond const &bond, std::string type1, XDIAG_RETHROW(e); } -bool coupling_defined(Bond const &bond, BondList const &bonds) { - if (bond.coupling_defined()) { - return true; - } else { - std::string name = bond.coupling_name(); - if (bonds.coupling_defined(name)) { - return true; - } else { - return false; - } - } -} - -bool matrix_defined(Bond const &bond, BondList const &bonds) { - if (bond.matrix_defined()) { - return true; - } else { - std::string type = bond.type(); - if (bonds.matrix_defined(type)) { - return true; - } else { - return false; - } - } -} - -complex coupling(Bond const &bond, BondList const &bonds) try { - if (!coupling_defined(bond, bonds)) { - XDIAG_THROW( - fmt::format("coupling \"{}\" is neither defined by Bond nor BondList", - bond.coupling_name())); - } - if (bond.coupling_defined()) { - return bond.coupling(); - } else { - std::string name = bond.coupling_name(); - return bonds.coupling(name); - } -} catch (Error const &e) { - XDIAG_RETHROW(e); - return 0.; -} - -arma::cx_mat matrix(Bond const &bond, BondList const &bonds) try { - if (!matrix_defined(bond, bonds)) { - XDIAG_THROW( - std::string( - "Error: the following matrix type is neither defined by Bond nor " - "BondList") + - bond.type()); - } - if (bond.matrix_defined()) { - return bond.matrix(); - } else { - std::string name = bond.type(); - return bonds.matrix(name); - } -} catch (Error const &e) { - XDIAG_RETHROW(e); - return arma::cx_mat(); -} - -BondList compile_explicit_couplings(BondList const &bonds, double precision, - std::string undefined_behavior) try { - BondList bonds_compiled; - - for (auto it : bonds.matrices()) { - bonds_compiled.set_matrix(it.first, it.second); - } - for (auto it : bonds.couplings()) { - bonds_compiled.set_coupling(it.first, it.second); - } - - for (Bond bond : bonds) { - - if (coupling_defined(bond, bonds)) { - complex cpl = coupling(bond, bonds); - - if (std::abs(cpl) < precision) { - continue; - } - - if (bond.type_defined()) { - bonds_compiled << Bond(bond.type(), cpl, bond.sites()); - } else { - bonds_compiled << Bond(bond.matrix(), cpl, bond.sites()); - } - } else { - if (undefined_behavior == "error") { - XDIAG_THROW(std::string("undefined coupling in bond: ") + - bond.coupling_name()); - } else if (undefined_behavior == "keep") { - bonds_compiled << bond; - } else if (undefined_behavior == "ignore") { - } else { - XDIAG_THROW("undefined_behaviour must be one of error/keep/ignore"); - } - } - } - return bonds_compiled; -} catch (Error const &e) { - XDIAG_RETHROW(e); - return BondList(); -} - -BondList compile_explicit_matrices(BondList const &bonds, double precision, - std::string undefined_behavior) try { - BondList bonds_compiled; - for (auto it : bonds.matrices()) { - bonds_compiled.set_matrix(it.first, it.second); - } - for (auto it : bonds.couplings()) { - bonds_compiled.set_coupling(it.first, it.second); - } - - for (Bond bond : bonds) { - - if (matrix_defined(bond, bonds)) { - arma::cx_mat mat = matrix(bond, bonds); - - // Go through matrix and set small elements to zero - for (arma::uword j = 0; j < mat.n_cols; ++j) { - for (arma::uword i = 0; i < mat.n_rows; ++i) { - if (std::abs(mat(i, j)) < precision) { - mat(i, j) = 0.; - } - } - } - - if (arma::norm(mat) < precision) { - continue; - } - - if (bond.coupling_defined()) { - bonds_compiled << Bond(mat, bond.coupling(), bond.sites()); - } else { - bonds_compiled << Bond(mat, bond.coupling_name(), bond.sites()); - } - } else { - if (undefined_behavior == "error") { - XDIAG_THROW(std::string("undefined matrix in bond: ") + - bond.coupling_name()); - } else if (undefined_behavior == "keep") { - bonds_compiled << bond; - } else if (undefined_behavior == "ignore") { - } else { - XDIAG_THROW( - "Error: undefined_behaviour must be one of error/keep/ignore"); - } - } - } - return bonds_compiled; -} catch (Error const &e) { - XDIAG_RETHROW(e); - return BondList(); -} - -BondList compile_explicit(BondList const &bonds, double precision, - std::string undefined_behavior) try { - BondList bonds_compiled; - bonds_compiled = - compile_explicit_couplings(bonds, precision, undefined_behavior); - bonds_compiled = - compile_explicit_matrices(bonds_compiled, precision, undefined_behavior); - return bonds_compiled; -} catch (Error const &e) { - XDIAG_RETHROW(e); - return BondList(); -} - } // namespace xdiag::operators diff --git a/xdiag/operators/compiler.hpp b/xdiag/operators/compiler.hpp index b9de18fc..76b8cc65 100644 --- a/xdiag/operators/compiler.hpp +++ b/xdiag/operators/compiler.hpp @@ -5,6 +5,8 @@ namespace xdiag::operators { +BondList clean_zeros(BondList const &bonds, double precision = 1e-12); + // Checks for the correct format of bonds void check_bond(Bond const &bond, int64_t n_sites_total, int64_t n_sites_bond, bool disjoint, std::string type); @@ -16,17 +18,4 @@ void check_bond_coupling_has_type(Bond const &bond, std::string type); void check_bond_coupling_has_type(Bond const &bond, std::string type1, std::string type2); -// bool coupling_defined(Bond const &bond, BondList const &bonds); -// bool matrix_defined(Bond const &bond, BondList const &bonds); - -// complex coupling(Bond const &bond, BondList const &bonds); -// arma::cx_mat matrix(Bond const &bond, BondList const &bonds); - -// BondList compile_explicit_couplings(BondList const &bonds, double precision, -// std::string undefined_behavior); -// BondList compile_explicit_matrices(BondList const &bonds, double precision, -// std::string undefined_behavior); -// BondList compile_explicit(BondList const &bonds, double precision, -// std::string undefined_behavior); - } // namespace xdiag::operators diff --git a/xdiag/operators/coupling.cpp b/xdiag/operators/coupling.cpp new file mode 100644 index 00000000..c041d8c5 --- /dev/null +++ b/xdiag/operators/coupling.cpp @@ -0,0 +1,355 @@ +#include "coupling.hpp" + +namespace xdiag { + +Coupling::Coupling(std::string value) : value_(value) {} +Coupling::Coupling(double value) : value_(value) {} +Coupling::Coupling(complex value) : value_(value) {} +Coupling::Coupling(arma::mat const &value) : value_(value) {} +Coupling::Coupling(arma::cx_mat const &value) : value_(value) {} +Coupling::Coupling(variant_t const &value) : value_(value) {} + +bool Coupling::isreal() const try { + return std::visit(overload{ + [](std::string const &) { + XDIAG_THROW("Cannot determine whether coupling is " + "real, since it is of std::string type"); + return false; + }, + [](double const &) { return true; }, + [](complex const &) { return false; }, + [](arma::mat const &) { return true; }, + [](arma::cx_mat const &) { return false; }, + }, + value_); +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +bool Coupling::ismatrix() const { + return std::visit(overload{ + [](std::string const &) { return false; }, + [](double const &) { return false; }, + [](complex const &) { return false; }, + [](arma::mat const &) { return true; }, + [](arma::cx_mat const &) { return true; }, + }, + value_); +} + +bool Coupling::isexplicit() const { + return std::visit(overload{ + [](std::string const &) { return false; }, + [](double const &) { return true; }, + [](complex const &) { return true; }, + [](arma::mat const &) { return true; }, + [](arma::cx_mat const &) { return true; }, + }, + value_); +} + +template bool Coupling::is() const { + return std::holds_alternative(value_); +} +template bool Coupling::is() const; +template bool Coupling::is() const; +template bool Coupling::is() const; +template bool Coupling::is() const; +template bool Coupling::is() const; + +template coeff_t Coupling::as() const try { + if (const coeff_t *pval = std::get_if(&value_)) { + return *pval; + } else { + XDIAG_THROW(fmt::format( + "Cannot retrieve desired value, since Coupling is of type \"{}\"", + type())); + } +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +template std::string Coupling::as() const; +template double Coupling::as() const; +template arma::mat Coupling::as() const; + +template <> complex Coupling::as() const try { + if (const double *pval = std::get_if(&value_)) { + return (complex)*pval; + } else if (const complex *pval = std::get_if(&value_)) { + return *pval; + } else { + XDIAG_THROW(fmt::format( + "Cannot retrieve desired value, since Coupling is of type \"{}\"", + type())); + } +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +template <> arma::cx_mat Coupling::as() const try { + if (const arma::mat *pval = std::get_if(&value_)) { + return arma::cx_mat( + *pval, arma::mat(pval->n_rows, pval->n_cols, arma::fill::zeros)); + } else if (const arma::cx_mat *pval = std::get_if(&value_)) { + return *pval; + } else { + XDIAG_THROW(fmt::format( + "Cannot retrieve desired value, since Coupling is of type \"{}\"", + type())); + } +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +std::string Coupling::type() const { + return std::visit(overload{ + [](std::string const &) { return "std::string"; }, + [](double const &) { return "double"; }, + [](complex const &) { return "complex"; }, + [](arma::mat const &) { return "arma::mat"; }, + [](arma::cx_mat const &) { return "arma::cx_mat"; }, + }, + value_); +} + +bool Coupling::operator==(Coupling const &rhs) const { + return std::visit( + overload{ + [](std::string const &a, std::string const &b) { return a == b; }, + [](double const &a, double const &b) { return a == b; }, + [](complex const &a, complex const &b) { return a == b; }, + [](arma::mat const &a, arma::mat const &b) { + return arma::all(vectorise(a) == vectorise(b)); + }, + [](arma::cx_mat const &a, arma::cx_mat const &b) { + return arma::all(vectorise(a) == vectorise(b)); + }, + [](auto const &, auto const &) { return false; }}, + value_, rhs.value_); +} + +bool Coupling::operator!=(Coupling const &rhs) const { + return !operator==(rhs); +} + +Coupling &Coupling::operator=(std::string rhs) { + value_ = rhs; + return *this; +} +Coupling &Coupling::operator=(double rhs) { + value_ = rhs; + return *this; +} +Coupling &Coupling::operator=(complex rhs) { + value_ = rhs; + return *this; +} +Coupling &Coupling::operator=(arma::mat const &rhs) { + value_ = rhs; + return *this; +} +Coupling &Coupling::operator=(arma::cx_mat const &rhs) { + value_ = rhs; + return *this; +} + +Coupling &Coupling::operator+=(Coupling const &rhs) try { + std::visit( + overload{[](double &a, double const &b) { a += b; }, + [&](double &a, complex const &b) { value_ = (complex)a + b; }, + [](complex &a, double const &b) { a += b; }, + [](complex &a, complex const &b) { a += b; }, + [](arma::mat &a, arma::mat const &b) { a += b; }, + [&](arma::mat &a, arma::cx_mat const &b) { + auto acplx = arma::cx_mat( + a, arma::mat(a.n_rows, a.n_cols, arma::fill::zeros)); + value_ = arma::cx_mat(acplx + b); + }, + [](arma::cx_mat &a, arma::mat const &b) { + auto bcplx = arma::cx_mat( + b, arma::mat(a.n_rows, a.n_cols, arma::fill::zeros)); + a += bcplx; + }, + [](arma::cx_mat &a, arma::cx_mat const &b) { a += b; }, + [&](auto &&, auto &&) { + XDIAG_THROW(fmt::format( + "Cannot add couplings of type \"{}\" and \"{}\"", type(), + rhs.type())); + }}, + value_, rhs.value_); + return *this; +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +Coupling &Coupling::operator-=(Coupling const &rhs) try { + std::visit( + overload{[](double &a, double const &b) { a -= b; }, + [&](double &a, complex const &b) { value_ = (complex)a - b; }, + [](complex &a, double const &b) { a -= b; }, + [](complex &a, complex const &b) { a -= b; }, + [](arma::mat &a, arma::mat const &b) { a -= b; }, + [&](arma::mat &a, arma::cx_mat const &b) { + auto acplx = arma::cx_mat( + a, arma::mat(a.n_rows, a.n_cols, arma::fill::zeros)); + value_ = arma::cx_mat(acplx - b); + }, + [](arma::cx_mat &a, arma::mat const &b) { + auto bcplx = arma::cx_mat( + b, arma::mat(a.n_rows, a.n_cols, arma::fill::zeros)); + a -= bcplx; + }, + [](arma::cx_mat &a, arma::cx_mat const &b) { a -= b; }, + [&](auto &&, auto &&) { + XDIAG_THROW(fmt::format( + "Cannot subtract couplings of type \"{}\" and \"{}\"", + type(), rhs.type())); + }}, + value_, rhs.value_); + return *this; +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +Coupling &Coupling::operator*=(double rhs) try { + std::visit( + overload{ + [](std::string const &) { + XDIAG_THROW( + "Cannot scalar multiply a coupling of type \"std::string\""); + }, + [&](double &a) { a *= rhs; }, + [&](complex &a) { a *= rhs; }, + [&](arma::mat &a) { a *= rhs; }, + [&](arma::cx_mat &a) { a *= rhs; }, + }, + value_); + return *this; +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +Coupling &Coupling::operator*=(complex rhs) try { + std::visit( + overload{ + [](std::string const &) { + XDIAG_THROW( + "Cannot scalar multiply a coupling of type \"std::string\""); + }, + [&](double &a) { value_ = (complex)a * rhs; }, + [&](complex &a) { a *= rhs; }, + [&](arma::mat &a) { + value_ = + arma::cx_mat(arma::cx_mat(a, arma::mat(a.n_rows, a.n_cols, + arma::fill::zeros)) * + rhs); + }, + [&](arma::cx_mat &a) { a *= rhs; }, + }, + value_); + return *this; +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +Coupling::variant_t Coupling::value() const { return value_; } + +Coupling &operator/=(Coupling &a, double b) try { + return a *= (1 / b); +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +Coupling &operator/=(Coupling &a, complex b) try { + return a *= (1 / b); +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +Coupling operator+(Coupling const &a, Coupling const &b) try { + Coupling ret = a; + ret += b; + return ret; +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +Coupling operator-(Coupling const &a, Coupling const &b) try { + Coupling ret = a; + ret -= b; + return ret; +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +Coupling operator*(double scalar, Coupling const &x) try { + Coupling ret = x; + ret *= scalar; + return ret; +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +Coupling operator*(complex scalar, Coupling const &x) try { + Coupling ret = x; + ret *= scalar; + return ret; +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +Coupling operator*(Coupling const &x, double scalar) try { + Coupling ret = x; + ret *= scalar; + return ret; +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +Coupling operator*(Coupling const &x, complex scalar) try { + Coupling ret = x; + ret *= scalar; + return ret; +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +Coupling operator/(double scalar, Coupling const &x) try { + Coupling ret = x; + ret /= scalar; + return ret; +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +Coupling operator/(complex scalar, Coupling const &x) try { + Coupling ret = x; + ret /= scalar; + return ret; +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +Coupling operator/(Coupling const &x, double scalar) try { + Coupling ret = x; + ret /= scalar; + return ret; +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +Coupling operator/(Coupling const &x, complex scalar) try { + Coupling ret = x; + ret /= scalar; + return ret; +} catch (Error const &e) { + XDIAG_RETHROW(e); +} + +std::ostream &operator<<(std::ostream &out, Coupling const &cpl) { + std::visit([&out](auto &&val) { out << val; }, cpl.value()); + return out; +} + +} // namespace xdiag diff --git a/xdiag/operators/coupling.hpp b/xdiag/operators/coupling.hpp new file mode 100644 index 00000000..ef159c72 --- /dev/null +++ b/xdiag/operators/coupling.hpp @@ -0,0 +1,71 @@ +#pragma once + +#include +#include +#include +#include + +namespace xdiag { + +class Coupling { +public: + using variant_t = + std::variant; + Coupling() = default; + explicit Coupling(std::string value); + explicit Coupling(double value); + explicit Coupling(complex value); + explicit Coupling(arma::mat const &value); + explicit Coupling(arma::cx_mat const &value); + explicit Coupling(variant_t const &value); + + bool isreal() const; + bool ismatrix() const; + bool isexplicit() const; + + template bool is() const; + template coeff_t as() const; + + std::string type() const; + + bool operator==(Coupling const &rhs) const; + bool operator!=(Coupling const &rhs) const; + + Coupling &operator=(std::string rhs); + Coupling &operator=(double rhs); + Coupling &operator=(complex rhs); + Coupling &operator=(arma::mat const &rhs); + Coupling &operator=(arma::cx_mat const &rhs); + + // vector space operations + Coupling &operator+=(Coupling const &rhs); + Coupling &operator-=(Coupling const &rhs); + + Coupling &operator*=(double rhs); + Coupling &operator*=(complex rhs); + + variant_t value() const; + +private: + variant_t value_; +}; + +Coupling &operator/=(Coupling &a, double b); +Coupling &operator/=(Coupling &a, complex b); + +Coupling operator+(Coupling const &a, Coupling const &b); +Coupling operator-(Coupling const &a, Coupling const &b); + +Coupling operator*(double scalar, Coupling const &x); +Coupling operator*(complex scalar, Coupling const &x); +Coupling operator*(Coupling const &x, double scalar); +Coupling operator*(Coupling const &x, complex scalar); + +Coupling operator/(double scalar, Coupling const &x); +Coupling operator/(complex scalar, Coupling const &x); +Coupling operator/(Coupling const &x, double scalar); +Coupling operator/(Coupling const &x, complex scalar); + +std::ostream &operator<<(std::ostream &out, Coupling const &cpl); + +} // namespace xdiag diff --git a/xdiag/operators/non_branching_bonds.cpp b/xdiag/operators/non_branching_bonds.cpp index 2b048680..133368d8 100644 --- a/xdiag/operators/non_branching_bonds.cpp +++ b/xdiag/operators/non_branching_bonds.cpp @@ -72,33 +72,34 @@ decompose_matrix_to_nonbranching(arma::Mat const &mat, } mats_nb.push_back(mat_nb); } - - // Create the bonds from the nonbranchin matrices - BondList bonds_nb; - for (auto mat_nb : mats_nb) { - bonds_nb << Bond("NONBRANCHINGBOND", mat_nb, bond.sites()); - } - return bonds_nb; - + return mats_nb; } catch (Error const &e) { XDIAG_RETHROW(e); - return BondListstd::vector>(); } BondList non_branching_bonds(Bond const &bond, double precision) try { - if (bond.coupling_is()) { - arma::mat mat = bond.coupling(); - return decompose_matrix_to_nonbranching(mat, precision); - } else if (bond.coupling_is()) { - arma::cx_mat mat = bond.coupling(); - return decompose_matrix_to_nonbranching(mat, precision); + if (bond.coupling().is()) { + arma::mat mat = bond.coupling().as(); + auto mats_nb = decompose_matrix_to_nonbranching(mat, precision); + BondList bonds_nb; + for (auto mat_nb : mats_nb) { + bonds_nb += Bond("NONBRANCHINGBOND", mat_nb, bond.sites()); + } + return bonds_nb; + } else if (bond.coupling().is()) { + arma::cx_mat mat = bond.coupling().as(); + auto mats_nb = decompose_matrix_to_nonbranching(mat, precision); + BondList bonds_nb; + for (auto mat_nb : mats_nb) { + bonds_nb += Bond("NONBRANCHINGBOND", mat_nb, bond.sites()); + } + return bonds_nb; } else { XDIAG_THROW( "Cannot convert bond to nonbranching bonds. Coupling is not a matrix"); } } catch (Error const &e) { XDIAG_RETHROW(e); - return BondList(); } BondList non_branching_bonds(BondList const &bonds, double precision) { @@ -110,7 +111,7 @@ BondList non_branching_bonds(BondList const &bonds, double precision) { } template -bool is_non_branching_matrix(arma::Mat const &mat) { +bool is_non_branching_matrix(arma::Mat const &mat, double precision) { for (arma::uword i = 0; i < mat.n_rows; ++i) { int64_t non_zero_in_row = 0; for (arma::uword j = 0; j < mat.n_cols; ++j) { @@ -122,15 +123,16 @@ bool is_non_branching_matrix(arma::Mat const &mat) { return false; } } + return true; } bool is_non_branching_bond(Bond const &bond, double precision) { - if (bond.coupling_is()) { - arma::mat mat = bond.coupling(); - return is_non_branchin_matrix(mat, precision); - } else if (bond.coupling_is()) { - arma::cx_mat mat = bond.coupling(); - return is_non_branchin_matrix(mat, precision); + if (bond.coupling().is()) { + arma::mat mat = bond.coupling().as(); + return is_non_branching_matrix(mat, precision); + } else if (bond.coupling().is()) { + arma::cx_mat mat = bond.coupling().as(); + return is_non_branching_matrix(mat, precision); } else { return false; } @@ -140,18 +142,15 @@ template NonBranchingBond::NonBranchingBond(Bond const &bond, double precision) try : sites_(bond.sites()), dim_(1 << sites_.size()), mask_(0) { - if (!is_non_branching_bond(bond, precision)) { - XDIAG_THROW( - "Error: trying to create a NonBranchingBond from a Bond which is " - "branching"); + if (!bond.ismatrix()) { + XDIAG_THROW("Error constructing NonBranchingBond: the coupling of the bond " + "must be an explicit matrix"); } - if (!bond.coupling_defined()) { + if (!is_non_branching_bond(bond, precision)) { XDIAG_THROW( - "Error: cannot create Nonbranching bond from Bond without having " - "its coupling defined."); + "trying to create a NonBranchingBond from a Bond which is branching"); } - coeff_t cpl = bond.coupling(); for (auto s : bond.sites()) { mask_ |= ((bit_t)1 << s); @@ -162,7 +161,7 @@ NonBranchingBond::NonBranchingBond(Bond const &bond, mask_ = ~mask_; // Log("Y: {}", BSTR(mask_)); - arma::cx_mat matrix_ = bond.matrix(); + arma::cx_mat matrix_ = bond.coupling().as(); // Matrix dimension is 2**(no. sites of bond) if ((matrix_.n_cols != dim_) || (matrix_.n_rows != dim_)) { @@ -186,9 +185,9 @@ NonBranchingBond::NonBranchingBond(Bond const &bond, XDIAG_THROW("Error: trying to create a real NonBranchingBond, but " "found a truly complex matrix entry"); } - coeff_[in] = cpl * real(matrix_(out, in)); + coeff_[in] = real(matrix_(out, in)); } else { - coeff_[in] = cpl * matrix_(out, in); + coeff_[in] = matrix_(out, in); } // ++non_zero_in_row; } diff --git a/xdiag/operators/non_branching_bonds.hpp b/xdiag/operators/non_branching_bonds.hpp index cd9abe92..72e6048b 100644 --- a/xdiag/operators/non_branching_bonds.hpp +++ b/xdiag/operators/non_branching_bonds.hpp @@ -7,7 +7,6 @@ namespace xdiag::operators { BondList non_branching_bonds(Bond const &bond, double precision = 1e-12); BondList non_branching_bonds(BondList const &bonds, double precision = 1e-12); - bool is_non_branching_bond(Bond const &bond, double precision = 1e-12); template class NonBranchingBond { diff --git a/xdiag/operators/symmetrized_operator.cpp b/xdiag/operators/symmetrized_operator.cpp index ab16ce9e..b3ff8854 100644 --- a/xdiag/operators/symmetrized_operator.cpp +++ b/xdiag/operators/symmetrized_operator.cpp @@ -1,5 +1,7 @@ #include "symmetrized_operator.hpp" +#include + namespace xdiag { BondList symmetrized_operator(Bond const &bond, PermutationGroup const &group) { @@ -21,63 +23,31 @@ BondList symmetrized_operator(BondList const &bonds, BondList symmetrized_operator(BondList const &bonds, PermutationGroup const &group, - Representation const &irrep) { + Representation const &irrep) try { BondList bonds_sym; int64_t N_group = group.size(); + BondList bonds_explicit = make_explicit(bonds); + for (auto bond : bonds_explicit) { - for (auto bond : bonds) { - - complex coupling = 0.; - - if (bond.coupling_named()) { - std::string name = bond.coupling_name(); - if (bonds.coupling_defined(name)) { - coupling = bonds.coupling(name); - } else { - Log.err( - "Error in symmetrized_operator: coupling with name {} undefined " - "in BondList", - name); - } - } else { - coupling = bond.coupling(); - } - - if (bond.type_defined()) { - std::string type = bond.type(); - - // Create all symmetrized bonds - for (int64_t i = 0; i < N_group; ++i) { - Permutation perm = group[i]; - complex bloch = irrep.character(i); - - std::vector sites_sym(bond.size(), 0); - for (int64_t site_idx = 0; site_idx < bond.size(); ++site_idx) { - sites_sym[site_idx] = perm[bond[site_idx]]; - } - complex cpl_sym = bloch * coupling / (complex)N_group; - bonds_sym << Bond(type, cpl_sym, sites_sym); - } - } else { - arma::cx_mat mat = bond.matrix(); - - // Create all symmetrized bonds - for (int64_t i = 0; i < N_group; ++i) { - Permutation perm = group[i]; - complex bloch = irrep.character(i); + std::string type = bond.type(); + Coupling coupling = bond.coupling(); - std::vector sites_sym(bond.size(), 0); - for (int64_t site_idx = 0; site_idx < bond.size(); ++site_idx) { - sites_sym[site_idx] = perm[bond[site_idx]]; - } + // Create all symmetrized bonds + for (int64_t i = 0; i < N_group; ++i) { + Permutation perm = group[i]; + complex bloch = irrep.character(i); - complex cpl_sym = bloch * coupling / (complex)N_group; - bonds_sym << Bond(mat, cpl_sym, sites_sym); + std::vector sites_sym(bond.size(), 0); + for (int64_t site_idx = 0; site_idx < bond.size(); ++site_idx) { + sites_sym[site_idx] = perm[bond[site_idx]]; } + Coupling cpl_sym = bloch * coupling / (complex)N_group; + bonds_sym += Bond(type, cpl_sym, sites_sym); } } - return bonds_sym; +} catch (Error const &e) { + XDIAG_RETHROW(e); } } // namespace xdiag diff --git a/xdiag/states/state.cpp b/xdiag/states/state.cpp index 6b3c8ac7..354a2b58 100644 --- a/xdiag/states/state.cpp +++ b/xdiag/states/state.cpp @@ -128,7 +128,6 @@ State::State(block_t const &block, arma::Mat const &matrix) int64_t State::n_sites() const { return xdiag::n_sites(block_); } bool State::isreal() const { return real_; } -bool State::iscomplex() const { return !isreal(); } State State::real() const try { if (isreal()) { diff --git a/xdiag/states/state.hpp b/xdiag/states/state.hpp index a81f1260..f297ebb6 100644 --- a/xdiag/states/state.hpp +++ b/xdiag/states/state.hpp @@ -8,7 +8,6 @@ namespace xdiag { class State { public: - // Interface State() = default; State(block_variant_t const &block, bool real = true, int64_t n_cols = 1); @@ -24,7 +23,6 @@ class State { int64_t n_sites() const; bool isreal() const; - bool iscomplex() const; State real() const; State imag() const; diff --git a/xdiag/utils/error.hpp b/xdiag/utils/error.hpp index 76051b2b..3844e396 100644 --- a/xdiag/utils/error.hpp +++ b/xdiag/utils/error.hpp @@ -37,8 +37,12 @@ void check_dimension_works_with_blas_int_size(int64_t dim); } // namespace xdiag +// comment: throw 0; necessary to let compiler know there need not be a return +// value #define XDIAG_THROW(message) \ - xdiag::throw_error(message, __FILE__, __func__, __LINE__); + xdiag::throw_error(message, __FILE__, __func__, __LINE__); \ + throw 0; #define XDIAG_RETHROW(error) \ - xdiag::rethrow_error(error, __FILE__, __func__, __LINE__); + xdiag::rethrow_error(error, __FILE__, __func__, __LINE__); \ + throw 0; diff --git a/xdiag/utils/print.cpp b/xdiag/utils/print.cpp index dafa3918..705c6076 100644 --- a/xdiag/utils/print.cpp +++ b/xdiag/utils/print.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include @@ -73,8 +74,8 @@ void print_pretty(const char *identifier, BondList const &bondlist) { if (bondlist.couplings().size() > 0) { printf("Couplings:\n"); for (auto name : bondlist.couplings()) { - auto cpl = bondlist[name]; - std::visit([](auto &&arg) { std::cout << arg << "\n"; }, cpl); + Coupling cpl = bondlist[name]; + std::cout << name << ": " << cpl << "\n"; } } }