diff --git a/include/dca/phys/dca_step/symmetrization/solve_orbital_op_signs.hpp b/include/dca/phys/dca_step/symmetrization/solve_orbital_op_signs.hpp new file mode 100644 index 000000000..ae2d5768f --- /dev/null +++ b/include/dca/phys/dca_step/symmetrization/solve_orbital_op_signs.hpp @@ -0,0 +1,214 @@ +// Copyright (C) 2026 UT-Battelle, LLC +// All rights reserved. +// +// See LICENSE for terms of usage. +// See CITATION.md for citation guidelines, if DCA++ is used for scientific publications. +// +// Author: Tyler Sax (tylersax@gmail.com) +// +// Derives the orbital-operation matrix U_S of each point-group operation from H0(k) and populates +// cluster_symmetry::get_orbital_op(), without reading hand-coded Lattice::transformationSignOf*. +// +// Scope: for every analytic model, U_S is a real signed permutation U_S = D P, where the +// permutation P is already encoded in the symmetry table (the band image .second) and the only +// unknown is the diagonal of +/-1 signs D = diag(sigma). +// +// The symmetry table stores not the exact momentum image S k but its folded representative +// k_new = S k - G, and with intra-cell orbital offsets a_b the fold is not free: H0 picks up the +// diagonal folding phase V(G) = diag(phi_b), phi_b = e^{i G.a_b} (cluster_symmetry::get_fold_phase(), +// a real +/-1 for every shipped model). Substituting U_S = D P into the symmetry relation +// H0(S k) = U_S H0(k) U_S^dagger and folding S k back to k_new gives, entry by entry, +// +// sigma(b0) sigma(b1) H0(k)(b0,b1) = phi(image(b0)) phi(image(b1)) H0(k_new)(image(b0), image(b1)). +// +// So each nonzero coupling fixes the product sigma(b0) sigma(b1) to the +/-1 fold-corrected ratio +// of the two H0 entries. The phi factors matter: without them the solve absorbs the folding phase +// into the signs and returns a representative-convention-dependent dressing of U_S (e.g. masking +// the odd parity of p_x under a mirror at a zone-boundary momentum) instead of the intrinsic +// orbital transform. +// +// The analytic models shipped with DCA++ have at most n_b = 3 orbitals, so once the gauge +// sigma(0) = +1 is fixed there are at most 2^(n_b-1) <= 4 candidate sign vectors. We simply +// enumerate them and keep the one that satisfies every coupling constraint. Three conditions make +// an operation fail: +// +// * a magnitude mismatch (an entry vanishing on only one side) means the permutation is not a +// symmetry of H0; +// * a non-real or non-unit ratio means the op needs a non-signed-permutation U_S (genuine orbital +// mixing, not yet supported) or is not a symmetry; +// * no candidate sign vector satisfying all couplings means no signed-permutation U_S exists. + +#ifndef DCA_PHYS_DCA_STEP_SYMMETRIZATION_SOLVE_ORBITAL_OP_SIGNS_HPP +#define DCA_PHYS_DCA_STEP_SYMMETRIZATION_SOLVE_ORBITAL_OP_SIGNS_HPP + +#include +#include +#include +#include +#include + +#include "dca/function/domains.hpp" +#include "dca/function/function.hpp" +#include "dca/phys/domains/cluster/cluster_symmetry.hpp" +#include "dca/phys/domains/quantum/electron_band_domain.hpp" + +namespace dca { +namespace phys { +// dca::phys:: + +namespace detail { +// dca::phys::detail:: + +// Absolute threshold below which an H0 entry counts as a zero (no coupling). +constexpr double kCouplingZeroTol = 1e-10; +// Tolerance on the H0 ratio being a real +/-1: the imaginary part and the deviation of the +// magnitude from 1 must both be below this. +constexpr double kSignRealTol = 1e-9; + +// Solves the +/-1 signs of one operation s and writes that op's U_S block into u_s. fold_phase is +// the diagonal folding phase phi (callable as fold_phase(k, band, s)) that undresses the folded +// H0 samples, so the solved signs are the intrinsic orbital transform. Throws if op s is not a +// signed-permutation symmetry of H0. Shared by the whole-group populator -- which lets the throw +// propagate, since a model must be a genuine symmetry group -- and by the group verification, +// which catches it to classify the op. +template +void solveSignsForOp(int s, int nb, int nk, const SymFunc& sym, const FoldFunc& fold_phase, + const H0Function& H0, UFunc& u_s) { + // Band permutation P: row b's single entry lands in column image[b]. This is the geometry-derived + // half of U_S (which orbital maps to which); it is k-independent for a point-group op, so read it + // at k = 0. Only the signs that decorate it are unknown. + std::vector image(nb); + for (int b = 0; b < nb; ++b) + image[b] = sym(0, b, s).second; + + // Gather one sign-product constraint per nonzero coupling, over all cluster momenta. Each nonzero + // coupling forces sigma(b0) sigma(b1) to the +/-1 fold-corrected ratio of the two H0 entries. The + // ratio must be a real +/-1 for a signed permutation to exist at all; the two checks below + // (magnitude match, real +/-1) are independent of the signs themselves, so they are tested here + // and reported as the first two rejection branches. + struct Constraint { + int b0, b1, product; + }; + std::vector constraints; + for (int k = 0; k < nk; ++k) { + const int k_new = sym(k, 0, s).first; + for (int b0 = 0; b0 < nb; ++b0) + for (int b1 = 0; b1 < nb; ++b1) { + const std::complex lhs = H0(b0, 0, b1, 0, k); + const std::complex rhs = H0(image[b0], 0, image[b1], 0, k_new); + + const bool zero_lhs = std::abs(lhs) < kCouplingZeroTol; + const bool zero_rhs = std::abs(rhs) < kCouplingZeroTol; + if (zero_lhs && zero_rhs) + continue; // coupling absent on both sides: consistent, but no constraint. + if (zero_lhs != zero_rhs) + throw std::logic_error( + "solveOrbitalOpSignsFromH0: H0 magnitude mismatch for op " + std::to_string(s) + + " -- the band permutation is not a symmetry of H0 (a coupling vanishes on only one " + "side)."); + // Undress the fold before taking the ratio. phi lives at the image bands -- the rhs entry + // sits at (image(b0), image(b1)) -- with G determined by (k, s); the stored table carries + // the band as a free index, so the image-band lookup is direct. + const double fold = fold_phase(k, image[b0], s) * fold_phase(k, image[b1], s); + const std::complex ratio = lhs / (fold * rhs); + if (std::abs(std::imag(ratio)) > kSignRealTol || + std::abs(std::abs(ratio) - 1.) > kSignRealTol) + throw std::domain_error( + "solveOrbitalOpSignsFromH0: H0 ratio is not a real +/-1 for op " + std::to_string(s) + + " -- the operation requires a non-signed-permutation U_S (genuine orbital mixing) or " + "is not a symmetry of H0."); + + constraints.push_back({b0, b1, std::real(ratio) > 0. ? 1 : -1}); + } + } + + // Find the signs by trying every assignment with the gauge sigma(0) = +1 -- 2^(n_b-1) <= 4 + // candidates for the supported models. A candidate is U_S iff it satisfies every gathered product + // constraint; for a connected coupling graph (every analytic model) it is unique. The undetermined + // global sign of any orbital that never couples is left at +1 by starting the search from the + // all-plus assignment, matching the lowest-index-+1 gauge. + std::vector sigma(nb, 1); + bool found = false; + for (int mask = 0; mask < (1 << (nb - 1)) && !found; ++mask) { + for (int b = 1; b < nb; ++b) + sigma[b] = (mask & (1 << (b - 1))) ? -1 : 1; // sigma[0] stays +1, the gauge. + + found = true; + for (const Constraint& c : constraints) + if (sigma[c.b0] * sigma[c.b1] != c.product) { + found = false; + break; + } + } + // No sign vector satisfied every constraint, so no signed permutation reproduces H0 under this op + if (!found) + throw std::logic_error( + "solveOrbitalOpSignsFromH0: no consistent +/-1 signs for op " + std::to_string(s) + + " -- the H0 couplings admit no signed-permutation U_S."); + + // Materialize U_S = D P one entry at a time: sign sigma(b) at (row b, column image[b]). + for (int b = 0; b < nb; ++b) + u_s(b, image[b], s) = static_cast(sigma[b]); +} + +} // namespace detail +// dca::phys:: + +// Populates cluster_symmetry::get_orbital_op() by solving each op's signs from H0. +// Throws if any op is not a signed-permutation symmetry of H0 (a model handed to this must be a +// genuine symmetry group). +template +void solveOrbitalOpSignsFromH0(const H0Function& H0) { + using BDmn = func::dmn_0; + using SymDmn = typename domains::cluster_symmetry::sym_super_cell_dmn_t; + using CDmn = typename domains::cluster_symmetry::c_dmn_t; + + const auto& sym = domains::cluster_symmetry::get_symmetry_matrix(); + const auto& fold_phase = domains::cluster_symmetry::get_fold_phase(); + auto& u_s = domains::cluster_symmetry::get_orbital_op(); + u_s = 0.; + + const int nb = BDmn::dmn_size(); + const int nk = CDmn::dmn_size(); + const int n_ops = SymDmn::dmn_size(); + + for (int s = 0; s < n_ops; ++s) + detail::solveSignsForOp(s, nb, nk, sym, fold_phase, H0, u_s); +} + +// Returns the sorted op indices that pass the sign-consistency (H0-invariance) check -- the +// H0-verified symmetry group. Non-throwing: an op that fails any branch is simply omitted. The +// candidate ops come from the symmetry table, which today is the declared point group filtered by +// lattice geometry. +template +std::vector verifiedSymmetryOps(const H0Function& H0) { + using BDmn = func::dmn_0; + using SymDmn = typename domains::cluster_symmetry::sym_super_cell_dmn_t; + using CDmn = typename domains::cluster_symmetry::c_dmn_t; + + const auto& sym = domains::cluster_symmetry::get_symmetry_matrix(); + const auto& fold_phase = domains::cluster_symmetry::get_fold_phase(); + auto& u_s = domains::cluster_symmetry::get_orbital_op(); + u_s = 0.; + + const int nb = BDmn::dmn_size(); + const int nk = CDmn::dmn_size(); + const int n_ops = SymDmn::dmn_size(); + + std::vector verified; + for (int s = 0; s < n_ops; ++s) { + try { + detail::solveSignsForOp(s, nb, nk, sym, fold_phase, H0, u_s); + verified.push_back(s); + } + catch (const std::exception&) { + // op s is not a signed-permutation symmetry of H0 -> not part of the verified group. + } + } + return verified; +} + +} // namespace phys +} // namespace dca + +#endif // DCA_PHYS_DCA_STEP_SYMMETRIZATION_SOLVE_ORBITAL_OP_SIGNS_HPP diff --git a/include/dca/phys/domains/cluster/cluster_symmetry.hpp b/include/dca/phys/domains/cluster/cluster_symmetry.hpp index acbff272c..14697b9f9 100644 --- a/include/dca/phys/domains/cluster/cluster_symmetry.hpp +++ b/include/dca/phys/domains/cluster/cluster_symmetry.hpp @@ -54,11 +54,48 @@ class cluster_symmetry> { typedef func::dmn_variadic, sym_super_cell_dmn_t> symmetry_matrix_dmn_t; + // (band_row, band_col) x super-cell op: the orbital operation matrix U_S, one per op. + typedef func::dmn_variadic, sym_super_cell_dmn_t> orbital_op_dmn_t; + + // (cluster point) x super-cell op: the band-independent image of each point under each op -- today's + // get_symmetry_matrix().first promoted out of the pair. (Momentum images do not depend on the band.) + typedef func::dmn_variadic mapped_point_dmn_t; + static func::function, symmetry_matrix_dmn_t>& get_symmetry_matrix() { static func::function, symmetry_matrix_dmn_t> symmetry_matrix( "symmetry_matrix_super_cell"); return symmetry_matrix; } + + // Image of each cluster point under each super-cell op: the ".first" of get_symmetry_matrix(), + // split into its own band-independent accessor. Populated for the momentum cluster in + // set_symmetry_matrices, where the k-image (from linear_transform of k) does not depend on the + // band. The real-space point image IS band-dependent -- intra-cell orbital offsets move with the + // site, so it cannot be separated this way -- and is not stored here. M1 shadow accessor; the + // uniform symmetry relation reads G(k) = U_S G(mapped_point) U_S^dagger. + static func::function& get_mapped_point() { + static func::function mapped_point("mapped_point_super_cell"); + return mapped_point; + } + + // Diagonal folding phase V(G) = diag(e^{i G.a_alpha}) per (cluster point, band, super-cell op). + // Populated alongside get_symmetry_matrix() in set_symmetry_matrices. Runtime guard rejects (throws) + // any phase that is not +/-1, which holds for every shipped (symmorphic, point-symmetry-only) model. + // \todo Lift the storage to a complex phase to support a genuine non-+/-1 fold. + static func::function& get_fold_phase() { + static func::function fold_phase("fold_phase_super_cell"); + return fold_phase; + } + + // Orbital operation U_S per super-cell op: the band x band matrix s.t. the symmetry relation is + // G(k) = U_S G(S^-1 k) U_S^dagger. M1 shadow accessor; populated by a model-aware step (it needs + // the Hamiltonian / the model's sign data, neither visible in this cluster-domain layer). v1 stores + // it real -- the values are a signed permutation (+/-1) for every shipped symmorphic model; complex + // U_S to support genuine orbital mixing is a future generalization. + static func::function& get_orbital_op() { + static func::function orbital_op("orbital_op_super_cell"); + return orbital_op; + } }; } // domains diff --git a/include/dca/phys/domains/cluster/symmetrization_algorithms/set_symmetry_matrices.hpp b/include/dca/phys/domains/cluster/symmetrization_algorithms/set_symmetry_matrices.hpp index 7a231ad70..b4b2d14c8 100644 --- a/include/dca/phys/domains/cluster/symmetrization_algorithms/set_symmetry_matrices.hpp +++ b/include/dca/phys/domains/cluster/symmetrization_algorithms/set_symmetry_matrices.hpp @@ -12,6 +12,7 @@ #ifndef DCA_PHYS_DOMAINS_CLUSTER_SYMMETRIZATION_ALGORITHMS_SET_SYMMETRY_MATRICES_HPP #define DCA_PHYS_DOMAINS_CLUSTER_SYMMETRIZATION_ALGORITHMS_SET_SYMMETRY_MATRICES_HPP +#include #include #include @@ -151,6 +152,13 @@ void set_symmetry_matrices::set_k_symmetry_matrix() { sym_super_cell_dmn_t>>& k_symmetry_matrix = cluster_symmetry::get_symmetry_matrix(); // k_cluster_type::get_symmetry_matrix(); + func::function, sym_super_cell_dmn_t>>& + fold_phase = cluster_symmetry::get_fold_phase(); + + func::function>& mapped_point = + cluster_symmetry::get_mapped_point(); + // r_symmetry_matrix.print_fingerprint(); // k_symmetry_matrix.print_fingerprint(); @@ -159,12 +167,34 @@ void set_symmetry_matrices::set_k_symmetry_matrix() { for (int l = 0; l < sym_super_cell_dmn_t::dmn_size(); ++l) { std::vector k = k_dmn_t::get_elements()[i]; std::vector trafo_k(DIMENSION, 0); - auto sym_elem = sym_super_cell_dmn_t::get_elements()[l]; - auto kelems = k_dmn_t::get_elements(); sym_super_cell_dmn_t::get_elements()[l].linear_transform(&k[0], &trafo_k[0]); - k_symmetry_matrix(i, j, l).first = find_k_index(trafo_k); + const int k_image = find_k_index(trafo_k); + k_symmetry_matrix(i, j, l).first = k_image; k_symmetry_matrix(i, j, l).second = r_symmetry_matrix(i, j, l).second; + + // mapped_point is the same k-image, promoted to a band-independent accessor. trafo_k (hence + // k_image) is computed from the momentum alone, so it does not depend on the band j; store it + // once (j == 0) to make that independence explicit. + if (j == 0) + mapped_point(i, l) = k_image; + + // Folding phase V(G) = e^{i G.a_band}. find_k_index folds (O k) back into the BZ by + // subtracting the reciprocal vector G. Here, we recover G = (O k) - k[image] + // (it is otherwise discarded) and dot it with the band's intra-cell offset a_vec. + const std::vector& a_vec = b_dmn_t::get_elements()[j].a_vec; + double G_dot_a = 0.; + for (int d = 0; d < DIMENSION; ++d) + G_dot_a += (trafo_k[d] - k_dmn_t::get_elements()[k_image][d]) * a_vec[d]; + + // Fold phase is assumed to be +/-1 for every model currently supported. Store it real; + // reject a genuine non-+/-1 fold. + // \todo Lift the storage to a complex phase to support a genuine non-+/-1 fold. + if (std::abs(std::sin(G_dot_a)) > 1.e-6) + throw std::domain_error( + "set_symmetry_matrices: encountered a non-+/-1 folding phase, which is out of scope " + "for the current symmorphic / point-symmetry-only symmetrization (real fold phase)."); + fold_phase(i, j, l) = std::cos(G_dot_a); } } } diff --git a/test/integration/phys/symmetrization/symmetrize_characterization_test.cpp b/test/integration/phys/symmetrization/symmetrize_characterization_test.cpp index b80d1b4d8..558f10bf1 100644 --- a/test/integration/phys/symmetrization/symmetrize_characterization_test.cpp +++ b/test/integration/phys/symmetrization/symmetrize_characterization_test.cpp @@ -40,9 +40,12 @@ #include "dca/platform/dca_gpu.h" #include "dca/phys/dca_step/symmetrization/symmetrize.hpp" +#include "dca/phys/dca_step/symmetrization/solve_orbital_op_signs.hpp" #include #include +#include +#include #include #include #include @@ -103,6 +106,7 @@ struct SquareD4 { static std::vector expectedFailingReps() { return {}; } static std::vector expectedFailingKOps() { return {}; } static std::vector expectedFailingROps() { return {}; } + static std::vector expectedUnverifiedKOps() { return {}; } }; // CASE 2 -- three-band Hubbard (Emery/CuO2) model on D4. This case fails the @@ -121,6 +125,7 @@ struct ThreebandD4 { static std::vector expectedFailingReps() { return {"k_iw", "r_iw", "r_tau"}; } static std::vector expectedFailingKOps() { return {}; } static std::vector expectedFailingROps() { return {}; } + static std::vector expectedUnverifiedKOps() { return {}; } }; // CASE 3 -- 1D chain. Since this model has identity-only symmetry, this test exercises @@ -135,6 +140,7 @@ struct SinglebandChain { static std::vector expectedFailingReps() { return {}; } static std::vector expectedFailingKOps() { return {}; } static std::vector expectedFailingROps() { return {}; } + static std::vector expectedUnverifiedKOps() { return {}; } }; // Templated test fixture @@ -204,7 +210,43 @@ double maxAbsDiff(const F& a, const F& b) { return m; } -// TEST 1 of 4: DomainActivation -- sanity gate. Reads the op counts off the symmetry +// Runs fn and returns the what() string of any std::exception it throws, or "" if it ran +// to completion. The rejection fixtures use this to assert not just that the sign solver +// throws, but that it throws for the intended reason (the three branches catch three +// distinct corruptions and share two exception types, so the message is what distinguishes +// them). +template +std::string captureThrowMessage(F&& fn) { + try { + fn(); + } + catch (const std::exception& e) { + return e.what(); + } + return ""; +} + +// Location of an off-diagonal H0 coupling: the band pair (b0, b1) at momentum k, and its magnitude. +struct OffDiagTarget { + int b0, b1, k; + double mag; +}; + +// Finds the strongest off-diagonal coupling H0(b0, b1, k) (b0 != b1) over all momenta. Shared by the +// magnitude-mismatch and non-unit-ratio rejection fixtures, which perturb that one entry (and its +// Hermitian partner) in different ways. +template +OffDiagTarget findStrongestOffDiagonal(const F& H0, int nb, int nk) { + OffDiagTarget t{-1, -1, -1, 0.}; + for (int k = 0; k < nk; ++k) + for (int b0 = 0; b0 < nb; ++b0) + for (int b1 = 0; b1 < nb; ++b1) + if (b0 != b1 && std::abs(H0(b0, 0, b1, 0, k)) > t.mag) + t = {b0, b1, k, std::abs(H0(b0, 0, b1, 0, k))}; + return t; +} + +// TEST 1: DomainActivation -- sanity gate. Reads the op counts off the symmetry // domains and asserts they match the trait's expected_num_symmetries. Without this, a // misconfigured case could run green simply because it had only the identity op to check. TYPED_TEST(SymmetrizeCharacterizationTest, DomainActivation) { @@ -220,7 +262,7 @@ TYPED_TEST(SymmetrizeCharacterizationTest, DomainActivation) { EXPECT_EQ(n_r_ops, TypeParam::expected_num_symmetries); } -// TEST 2 of 4: CoarseNoOp -- the full production pipeline, asserted at 1e-12. Covers the +// TEST 2: CoarseNoOp -- the full production pipeline, asserted at 1e-12. Covers the // four representations the existing test exercises (k/r x w/t). For each, run // Symmetrize::execute on a copy of G0 and diff against the original. The set of // representations that are not left invariant is compared to the expected failing set @@ -287,7 +329,7 @@ TYPED_TEST(SymmetrizeCharacterizationTest, CoarseNoOp) { << "Coarse no-op red/green map changed vs the expected failing set."; } -// TEST 3 of 4: PerOpMapMomentum -- per-operation map (momentum channel). For each +// TEST 3: PerOpMapMomentum -- per-operation map (momentum channel). For each // operation s, mirror the momentum band-aware executeCluster relation on G0(k, iw): // G0(b0, b1, k) == sign_s * G0(b0_new, b1_new, k_new) when sign_s != 0, // where (k_new, b*_new) come from get_symmetry_matrix() and sign_s from @@ -354,7 +396,149 @@ TYPED_TEST(SymmetrizeCharacterizationTest, PerOpMapMomentum) { << "Momentum per-op red/green map changed vs the expected failing set."; } -// TEST 4 of 4: PerOpMapRealSpace -- mirror of TEST 3 but in real space / tau, and with one +// TEST 4: PerOpMapMomentumUS -- the shadow orbital-op record, derived from H0. Solve each U_S as +// a signed permutation directly from H0(k) by enumerating the sign assignments consistent with the +// fold-corrected H0 couplings (no transformationSignOf*), verify it is a well-formed unitary signed +// permutation, and re-run the per-op momentum no-op as the matrix conjugation +// G0(k) == (U_S V) G0(k_new) (U_S V)^dagger, where V = diag(fold_phase(k, ., s)) redresses the +// intrinsic U_S for the folded representatives the table stores -- rather than the element-wise +// scalar relation of TEST 3. It must reproduce TEST 3's red/green map; its purpose is to prove the +// consolidated, H0-derived matrix record is faithful and correct, not to find new bugs. +TYPED_TEST(SymmetrizeCharacterizationTest, PerOpMapMomentumUS) { + using Fixture = SymmetrizeCharacterizationTest; + using KClusterDmn = typename Fixture::KClusterDmn; + using WDmn = typename Fixture::WDmn; + using KCluster = typename Fixture::KCluster; + + // Populate the shadow orbital-op table from H0 alone; throws if the op is not a signed-permutation + // symmetry of H0. + dca::phys::solveOrbitalOpSignsFromH0(this->H0_); + const auto& u_s = dca::phys::domains::cluster_symmetry::get_orbital_op(); + const auto& sym = dca::phys::domains::cluster_symmetry::get_symmetry_matrix(); + const auto& fold_phase = dca::phys::domains::cluster_symmetry::get_fold_phase(); + + const int n_ops = Fixture::KSymDmn::dmn_size(); + const int nb = BDmn::dmn_size(); + const int ns = SDmn::dmn_size(); + const int nk = KClusterDmn::dmn_size(); + const int nw = WDmn::dmn_size(); + + // ---- Check 1 (structural): each U_S is a signed permutation -- exactly one +/-1 per row and column. + // (A real matrix with one unit-magnitude entry per row and per column is orthogonal by + // construction, so there is no separate U_S U_S^T = I check to run.) + for (int s = 0; s < n_ops; ++s) + for (int r = 0; r < nb; ++r) { + int row_nnz = 0, col_nnz = 0; + for (int c = 0; c < nb; ++c) { + if (u_s(r, c, s) != 0.) { + ++row_nnz; + EXPECT_DOUBLE_EQ(std::abs(u_s(r, c, s)), 1.0); + } + if (u_s(c, r, s) != 0.) + ++col_nnz; + } + EXPECT_EQ(row_nnz, 1) << "U_S op " << s << " row " << r << " is not a permutation row"; + EXPECT_EQ(col_nnz, 1) << "U_S op " << s << " col " << r << " is not a permutation column"; + } + + // ---- Check 2 (action): per-op no-op via conjugation G0(k) == (U_S V) G0(k_new) (U_S V)^dagger + // (all real). U_S is the intrinsic orbital transform; V = diag(fold_phase(k, ., s)) restores the + // folding phase of the representative k_new = S k - G the table stores, at the image (contracted) + // band indices. + function, dmn_variadic> G0; + dca::phys::compute_G0_k_w(this->H0_, this->parameters_.get_chemical_potential(), 1, G0); + + std::vector failing_ops; + std::cout << "[per-op momentum map, U_S conjugation] " << TypeParam::Input << "\n"; + for (int s = 0; s < n_ops; ++s) { + double max_viol = 0.; + for (int k = 0; k < nk; ++k) { + const int k_new = sym(k, 0, s).first; // k-image is band-independent + for (int sp = 0; sp < ns; ++sp) + for (int w = 0; w < nw; ++w) + for (int b0 = 0; b0 < nb; ++b0) + for (int b1 = 0; b1 < nb; ++b1) { + std::complex conj_val = 0.; + for (int a = 0; a < nb; ++a) + for (int b = 0; b < nb; ++b) + conj_val += u_s(b0, a, s) * fold_phase(k, a, s) * G0(a, sp, b, sp, k_new, w) * + fold_phase(k, b, s) * u_s(b1, b, s); + max_viol = std::max(max_viol, std::abs(G0(b0, sp, b1, sp, k, w) - conj_val)); + } + } + const bool failed = max_viol >= kTolerance; + std::cout << " op " << s << ": " << (failed ? "FAIL" : "PASS") + << " max|G0 - U_S.G0.U_S^dag| = " << max_viol << "\n"; + if (failed) + failing_ops.push_back(s); + } + + std::sort(failing_ops.begin(), failing_ops.end()); + auto expected_ops = TypeParam::expectedFailingKOps(); + std::sort(expected_ops.begin(), expected_ops.end()); + EXPECT_EQ(failing_ops, expected_ops) + << "U_S conjugation per-op momentum map changed vs the expected failing set."; +} + +// TEST 5: GroupShadowCrossCheck -- the headline U_S verification statement. The candidate ops +// in the symmetry table come from the declared point group (DCA_point_group), filtered by lattice +// geometry; DomainActivation already pins their count to the declared order. This test closes the +// loop at the Hamiltonian level: the set of ops the sign-consistency solve accepts (the H0-verified +// group) must equal the full declared/realized group. Equivalently, no declared op may be merely a +// lattice symmetry that H0 does not respect. The set of rejected ops is compared to the per-case +// expectation (empty on every model included in this test). When candidate ops are derived from H0, +// this same predicate becomes the actual filter that carves the real group out of a candidate table; +// here it is validated against the declared group as a known-correct oracle. +TYPED_TEST(SymmetrizeCharacterizationTest, GroupShadowCrossCheck) { + using Fixture = SymmetrizeCharacterizationTest; + using KCluster = typename Fixture::KCluster; + const int n_ops = Fixture::KSymDmn::dmn_size(); + + // The H0-verified group: ops the sign-consistency solve accepts (non-throwing). + const std::vector verified = dca::phys::verifiedSymmetryOps(this->H0_); + + // Anything the table declares but H0 does not verify. + std::vector rejected; + for (int s = 0; s < n_ops; ++s) + if (std::find(verified.begin(), verified.end(), s) == verified.end()) + rejected.push_back(s); + + std::cout << "[group shadow] " << TypeParam::Input << ": " << verified.size() << "/" << n_ops + << " declared ops H0-verified\n"; + + std::sort(rejected.begin(), rejected.end()); + auto expected_rejected = TypeParam::expectedUnverifiedKOps(); + std::sort(expected_rejected.begin(), expected_rejected.end()); + EXPECT_EQ(rejected, expected_rejected) + << "H0-verified group differs from the declared group: the set of declared ops the " + "sign-consistency check rejects changed vs the expectation."; +} + +// TEST 6: MappedPointShadow -- the mapped_point accessor. The band-independent ".first" of the +// legacy symmetry matrix (the image of each momentum under each op) is promoted into its own +// (cluster-point x op) accessor. Verify it reproduces .first for every band, which confirms both +// that the shadow accessor is faithful and that the momentum image really is band-independent. +TYPED_TEST(SymmetrizeCharacterizationTest, MappedPointShadow) { + using Fixture = SymmetrizeCharacterizationTest; + using KCluster = typename Fixture::KCluster; + using KClusterDmn = typename Fixture::KClusterDmn; + + const auto& sym = dca::phys::domains::cluster_symmetry::get_symmetry_matrix(); + const auto& mapped_point = dca::phys::domains::cluster_symmetry::get_mapped_point(); + + const int n_ops = Fixture::KSymDmn::dmn_size(); + const int nb = BDmn::dmn_size(); + const int nk = KClusterDmn::dmn_size(); + + for (int s = 0; s < n_ops; ++s) + for (int k = 0; k < nk; ++k) + for (int b = 0; b < nb; ++b) + EXPECT_EQ(mapped_point(k, s), sym(k, b, s).first) + << "mapped_point disagrees with legacy .first at (k=" << k << ", b=" << b << ", op=" << s + << ")"; +} + +// TEST 7: PerOpMapRealSpace -- mirror of TEST 3 but in real space / tau, and with one // key difference: it skips off-diagonal band pairs (b0!=b1) because production's real-space // executeCluster does the same. That skip is itself one of the suspected bugs, so this test // characterizes production's behavior, off-diag blind spot included. @@ -422,4 +606,163 @@ TYPED_TEST(SymmetrizeCharacterizationTest, PerOpMapRealSpace) { << "Real-space per-op red/green map changed vs the expected failing set."; } +// TEST 8 (a/b/c): rejection fixtures for the sign-consistency solver. +// +// solveOrbitalOpSignsFromH0 has three throw branches. No shipped model exercises any of the +// three on a correct H0 -- the spine cases all pass -- so we perturb a copy of the (exactly +// symmetric) H0 to break one operation and assert the solver throws for the intended reason. +// +// These need a multi-orbital model. A single band has no off-diagonal coupling to corrupt, +// and on the spine's square cluster every nonzero single-band entry sits at a k that is a +// fixed point of every op (the BZ center and corner), so no single-entry perturbation can +// break a symmetry there. Threeband (d, px, py) is the driver; square and singleband skip. + +// 8a/8b -- the two ways a single broken off-diagonal coupling is rejected. Both perturb the +// strongest off-diagonal H0 entry (and its Hermitian partner) so its symmetry image no longer +// matches, differing only in how: +// * vanish it (-> magnitude mismatch: a coupling present on one side of the relation but gone on +// the other, so the band permutation is no longer a symmetry of H0); +// * scale it by 1.5 (-> non-unit ratio: the H0 ratio that would fix the sign product is no longer +// a real +/-1, the signature of an op needing a non-signed-permutation U_S -- genuine orbital +// mixing -- or simply not a symmetry). +// They share everything but the perturbation, so they live in one test with two assertion blocks. +TYPED_TEST(SymmetrizeCharacterizationTest, RejectsBrokenOffDiagonalCoupling) { + using Fixture = SymmetrizeCharacterizationTest; + using KCluster = typename Fixture::KCluster; + const int nb = BDmn::dmn_size(); + const int nk = Fixture::KClusterDmn::dmn_size(); + if (nb < 2) { + std::cout << "[skip] " << TypeParam::Input + << ": needs a multi-orbital model with off-diagonal couplings.\n"; + return; + } + + const auto t = findStrongestOffDiagonal(this->H0_, nb, nk); + ASSERT_GT(t.mag, kTolerance) << "no off-diagonal coupling to perturb."; + + // Vanish: the op connecting this entry to its (still nonzero) image sees it gone on one side only. + { + auto H0p = this->H0_; + H0p(t.b0, 0, t.b1, 0, t.k) = 0.; + H0p(t.b1, 0, t.b0, 0, t.k) = 0.; + const std::string msg = + captureThrowMessage([&] { dca::phys::solveOrbitalOpSignsFromH0(H0p); }); + ASSERT_FALSE(msg.empty()) << "sign solver did not reject a vanished coupling."; + EXPECT_NE(msg.find("magnitude mismatch"), std::string::npos) << "actual: " << msg; + } + + // Scale by 1.5: the entry stays nonzero but its ratio to the image is not +/-1. + { + auto H0p = this->H0_; + H0p(t.b0, 0, t.b1, 0, t.k) *= 1.5; + H0p(t.b1, 0, t.b0, 0, t.k) *= 1.5; + const std::string msg = + captureThrowMessage([&] { dca::phys::solveOrbitalOpSignsFromH0(H0p); }); + ASSERT_FALSE(msg.empty()) << "sign solver did not reject a non-unit coupling ratio."; + EXPECT_NE(msg.find("not a real"), std::string::npos) << "actual: " << msg; + } +} + +// 8c -- inconsistent signs: pairwise sign products that are each a valid +/-1 but cannot be +// jointly satisfied by any assignment. At the zone corner all three bands couple (d-px, d-py, +// px-py), so the d-px and d-py signs force the px-py sign product, and a px<->py swap op +// reads the px-py sign off H0 directly. Flipping one direction of the px-py coupling +// makes the direct reading contradict the forced product: no sign vector works. +TYPED_TEST(SymmetrizeCharacterizationTest, RejectsInconsistentSigns) { + using Fixture = SymmetrizeCharacterizationTest; + using KCluster = typename Fixture::KCluster; + const int nb = BDmn::dmn_size(); + const int nk = Fixture::KClusterDmn::dmn_size(); + if (nb < 3) { + std::cout << "[skip] " << TypeParam::Input << ": needs >=3 coupled bands to form a sign cycle.\n"; + return; + } + + auto H0p = this->H0_; + // The zone corner is the only momentum where the px-py coupling H0(1,2) is nonzero (both + // sin(k/2) factors nonzero); there all three bands couple. + int kc = -1; + for (int k = 0; k < nk; ++k) + if (std::abs(H0p(1, 0, 2, 0, k)) > kTolerance) { + kc = k; + break; + } + ASSERT_GE(kc, 0) << "no zone-corner px-py coupling found."; + + // An op that swaps bands 1<->2 and fixes the corner reads the px-py sign as the raw H0 + // ratio; without one the contradiction never arises. + const auto& sym = dca::phys::domains::cluster_symmetry::get_symmetry_matrix(); + const int n_ops = Fixture::KSymDmn::dmn_size(); + bool have_swap = false; + for (int s = 0; s < n_ops; ++s) + if (sym(0, 1, s).second == 2 && sym(0, 2, s).second == 1 && sym(kc, 0, s).first == kc) { + have_swap = true; + break; + } + if (!have_swap) { + std::cout << "[skip] " << TypeParam::Input << ": no px<->py swap op fixing the zone corner.\n"; + return; + } + + H0p(1, 0, 2, 0, kc) *= -1.; // one direction only, on purpose. + + const std::string msg = + captureThrowMessage([&] { dca::phys::solveOrbitalOpSignsFromH0(H0p); }); + ASSERT_FALSE(msg.empty()) << "sign solver did not reject an inconsistent sign system."; + EXPECT_NE(msg.find("no consistent"), std::string::npos) << "actual: " << msg; +} + +// TEST 9: the folding phase must not be absorbed into U_S (PR #368 review). At a zone-boundary +// momentum the folded representative of the mirror image coincides with the original point, so the +// raw H0 ratio is +1 there; only the fold correction phi = e^{i G.a_b} exposes the intrinsic odd +// parity of p_x. The scenario is the threeband d-p_x block at k = (pi, 0) under the k_x -> -k_x +// mirror: the image (-pi, 0) folds back to (pi, 0) with G = (2 pi, 0), and a_px = (1/2, 0) gives +// phi_px = cos(pi) = -1 while phi_d = +1. Without the correction the solve returns sigma_px = +1 +// -- the representative-convention-dressed sign -- instead of the intrinsic -1. Mock-based and +// case-independent (unlike the fixtures above): the dressing is uniform across every momentum of +// the shipped clusters, so no perturbation of a real H0 can expose it; the single boundary point +// with an explicit fold functor is the minimal reproducer. +TEST(SolveOrbitalOpSigns, BoundaryFoldDoesNotMaskPxParity) { + // One momentum, one op: the mirror image folds back onto the same representative, no band + // permutation. + struct Sym { + std::pair operator()(int /*k*/, int b, int /*s*/) const { + return {0, b}; + } + }; + // The d-p_x block of the threeband H0 at k = (pi, 0) with t_pd = 1: 2 i t_pd sin(k_x / 2). + struct H0 { + std::complex operator()(int b0, int /*s0*/, int b1, int /*s1*/, int /*k*/) const { + const std::complex i(0., 1.); + if (b0 == 0 && b1 == 1) + return 2. * i; + if (b0 == 1 && b1 == 0) + return -2. * i; // Hermitian partner. + return 0.; + } + }; + // phi_b = cos(G . a_b) with G = (2 pi, 0): -1 for p_x (band 1, a = (1/2, 0)), +1 for d (at the + // origin) and p_y (a = (0, 1/2), orthogonal to G). + struct Fold { + double operator()(int /*k*/, int b, int /*s*/) const { + return b == 1 ? -1. : 1.; + } + }; + // Output stand-in for cluster_symmetry::get_orbital_op(): the solver writes the solved U_S + // here (band row, band column; the single op collapses the op index). + struct USTable { + double data[3][3] = {}; + double& operator()(int r, int c, int /*s*/) { + return data[r][c]; + } + }; + + USTable u_s; + dca::phys::detail::solveSignsForOp(0, 3, 1, Sym{}, Fold{}, H0{}, u_s); + + EXPECT_DOUBLE_EQ(u_s(0, 0, 0), 1.); + EXPECT_DOUBLE_EQ(u_s(1, 1, 0), -1.); // The intrinsic p_x parity, not the fold-dressed +1. + EXPECT_DOUBLE_EQ(u_s(2, 2, 0), 1.); +} + } // namespace