Symmetrization M2: derive each 2D model's point group and report divergence from the declaration

cmake/dca_config.cmake

@@ -100,23 +100,24 @@ endif()
 # TODO: Add more point groups and lattices.
 
 # Point group
-set(DCA_POINT_GROUP "D4" CACHE STRING "Point group symmetry, options are: C6 | D4 | no_symmetry<2> | no_symmetry<3>.")
-set_property(CACHE DCA_POINT_GROUP PROPERTY STRINGS C6 D4 no_symmetry<2> no_symmetry<3>)
+set(DCA_POINT_GROUP "D4" CACHE STRING "Point group symmetry, options are: D6 | D4 | no_symmetry<2> | no_symmetry<3>.")
+set_property(CACHE DCA_POINT_GROUP PROPERTY STRINGS D6 D4 no_symmetry<2> no_symmetry<3>)
 
-if (DCA_POINT_GROUP STREQUAL "C6")
+if (DCA_POINT_GROUP STREQUAL "D6")
   set(DCA_POINT_GROUP_INCLUDE
-    "dca/phys/domains/cluster/symmetries/point_groups/2d/2d_hexagonal.hpp")
+    "dca/phys/domains/cluster/symmetries/point_groups/2d/holohedries_2d.hpp")
 
 elseif (DCA_POINT_GROUP STREQUAL "D4")
-  set(DCA_POINT_GROUP_INCLUDE "dca/phys/domains/cluster/symmetries/point_groups/2d/2d_square.hpp")
+  set(DCA_POINT_GROUP_INCLUDE
+    "dca/phys/domains/cluster/symmetries/point_groups/2d/holohedries_2d.hpp")
 
 elseif (DCA_POINT_GROUP STREQUAL "no_symmetry<2>")
   set(DCA_POINT_GROUP_INCLUDE "dca/phys/domains/cluster/symmetries/point_groups/no_symmetry.hpp")
 elseif (DCA_POINT_GROUP STREQUAL "no_symmetry<3>")
   set(DCA_POINT_GROUP_INCLUDE "dca/phys/domains/cluster/symmetries/point_groups/no_symmetry.hpp")
 
 else()
-  message(FATAL_ERROR "Please set DCA_POINT_GROUP to a valid option: C6 | D4 | no_symmetry<2> | no_symmetry<3>.")
+  message(FATAL_ERROR "Please set DCA_POINT_GROUP to a valid option: D6 | D4 | no_symmetry<2> | no_symmetry<3>.")
 endif()
 
 # Lattice type

include/dca/phys/dca_step/symmetrization/derive_point_group.hpp

@@ -0,0 +1,202 @@
+// Copyright (C) 2018 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 a 2D model's point group -- the built-in 2D holohedry pool filtered by the lattice geometry
+// and by H0 invariance -- and compares it against the group realized from the model's declared
+// DCA_point_group. Production symmetrization continues to use the DECLARED group:
+// deriveAndComparePointGroup restores the declared symmetry state before returning and only
+// reports divergences. The reports build a per-model map of models whose declaration disagrees with
+// the symmetry encoded in their own Hamiltonian, which drives later declaration-free upgrades.
+
+#ifndef DCA_PHYS_DCA_STEP_SYMMETRIZATION_DERIVE_POINT_GROUP_HPP
+#define DCA_PHYS_DCA_STEP_SYMMETRIZATION_DERIVE_POINT_GROUP_HPP
+
+#include <cmath>
+#include <complex>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+#include <type_traits>
+#include <vector>
+
+#include "dca/function/domains.hpp"
+#include "dca/function/function.hpp"
+#include "dca/phys/dca_step/symmetrization/solve_orbital_op_signs.hpp"
+#include "dca/phys/domains/cluster/cluster_symmetry.hpp"
+#include "dca/phys/domains/cluster/symmetries/point_groups/2d/holohedries_2d.hpp"
+#include "dca/phys/domains/cluster/symmetrization_algorithms/search_maximal_symmetry_group.hpp"
+#include "dca/phys/domains/cluster/symmetrization_algorithms/set_symmetry_matrices.hpp"
+#include "dca/phys/domains/quantum/electron_band_domain.hpp"
+#include "dca/phys/domains/quantum/electron_spin_domain.hpp"
+#include "dca/phys/domains/quantum/point_group_symmetry_element.hpp"
+
+namespace dca {
+namespace phys {
+// dca::phys::
+
+namespace detail {
+// dca::phys::detail::
+
+// Geometric equality of two symmetry elements (same linear part .O and fractional translation .t),
+// at the tolerance the symmetry search uses to deduplicate operations.
+inline bool isSameSymmetryElement(const domains::point_group_symmetry_element& a,
+                                  const domains::point_group_symmetry_element& b) {
+  constexpr double tol = 1e-6;
+  if (a.DIMENSION != b.DIMENSION)
+    return false;
+  double diff = 0.;
+  for (int i = 0; i < a.DIMENSION * a.DIMENSION; ++i)
+    diff = std::max(diff, std::abs(a.O[i] - b.O[i]));
+  for (int i = 0; i < a.DIMENSION; ++i)
+    diff = std::max(diff, std::abs(a.t[i] - b.t[i]));
+  return diff < tol;
+}
+
+inline bool containsSymmetryElement(const std::vector<domains::point_group_symmetry_element>& set,
+                                    const domains::point_group_symmetry_element& op) {
+  for (const auto& element : set)
+    if (isSameSymmetryElement(element, op))
+      return true;
+  return false;
+}
+
+// Detects the initializeH0 lattice interface: the legacy lattices cannot feed the
+// H0-invariance gate, so derivation skips them.
+template <typename Lattice, typename Parameters, typename H0Type, typename = void>
+struct LatticeHasInitializeH0 : std::false_type {};
+template <typename Lattice, typename Parameters, typename H0Type>
+struct LatticeHasInitializeH0<Lattice, Parameters, H0Type,
+                              std::void_t<decltype(Lattice::initializeH0(
+                                  std::declval<const Parameters&>(), std::declval<H0Type&>()))>>
+    : std::true_type {};
+
+// Makes PointGroup the live group of one cluster family, overwriting the previously realized group:
+// clear the operation lists (the search only appends), re-search under the same geometric filter,
+// resize the per-operation functions to the new count, and rebuild the symmetry tables. Assumes the
+// family was already populated once by cluster_domain_symmetry_initializer.
+template <typename KCluster, typename PointGroup>
+void installPointGroup() {
+  using CS_k = domains::cluster_symmetry<KCluster>;
+  using CS_r = domains::cluster_symmetry<typename CS_k::dual_type>;
+  using Family = typename CS_k::cluster_family_type;
+
+  CS_k::sym_unit_cell_dmn_t::parameter_type::get_elements().clear();
+  CS_k::sym_super_cell_dmn_t::parameter_type::get_elements().clear();
+  domains::search_maximal_symmetry_group<Family, PointGroup, domains::UNIT_CELL>::execute();
+  domains::search_maximal_symmetry_group<Family, PointGroup, domains::SUPER_CELL>::execute();
+  CS_k::get_symmetry_matrix().reset();
+  CS_r::get_symmetry_matrix().reset();
+  CS_k::get_mapped_point().reset();
+  CS_k::get_fold_phase().reset();
+  CS_k::get_orbital_op().reset();
+  domains::set_symmetry_matrices<Family>::execute();
+}
+
+}  // namespace detail
+// dca::phys::
+
+// Derives the point group of a 2D model for one cluster family and returns a human-readable report
+// of how it diverges from the declared group -- empty when they agree. Must be called right after
+// cluster_domain_symmetry_initializer<RClusterDmn, DCA_point_group>::execute() for the same family;
+// the declared symmetry state is restored before returning, so production behavior is unchanged.
+// 3D models (no 3D pool yet) and legacy lattices (no initializeH0) skip derivation entirely.
+template <typename RClusterDmn, typename Model, typename Parameters>
+std::string deriveAndComparePointGroup(const Parameters& parameters) {
+  using Lattice = typename Model::lattice_type;
+  using RCluster = typename RClusterDmn::parameter_type;
+  using KCluster = typename domains::cluster_symmetry<RCluster>::dual_type;
+  using BDmn = func::dmn_0<domains::electron_band_domain>;
+  using SDmn = func::dmn_0<domains::electron_spin_domain>;
+  using NuDmn = func::dmn_variadic<BDmn, SDmn>;
+  using H0Type =
+      func::function<std::complex<double>, func::dmn_variadic<NuDmn, NuDmn, func::dmn_0<KCluster>>>;
+
+  if constexpr (Lattice::DIMENSION != 2 ||
+                !detail::LatticeHasInitializeH0<Lattice, Parameters, H0Type>::value) {
+    (void)parameters;
+    return "";
+  }
+  else {
+    using KSymDmn = typename domains::cluster_symmetry<KCluster>::sym_super_cell_dmn_t;
+    using Element = domains::point_group_symmetry_element;
+
+    // The realized declared group, populated by the initializer that just ran. Only the first
+    // get_size() entries are the live group: on a re-initialization (some tests run update_domains
+    // twice per process) the search zeroes the count but appends to the never-cleared element
+    // vector, so the vector can be longer than the group. The reference tracks the singleton
+    // through the group installs below; the copy preserves the declared ops.
+    const std::vector<Element>& realized = KSymDmn::parameter_type::get_elements();
+    const int n_declared = KSymDmn::parameter_type::get_size();
+    const std::vector<Element> declared(realized.begin(), realized.begin() + n_declared);
+
+    // Derive: candidate pool -> same geometric filter -> H0-invariance gate.
+    detail::installPointGroup<KCluster, domains::holohedry_pool_2D>();
+
+    H0Type H0;
+    Model::initializeH0(parameters, H0);
+
+    std::vector<Element> derived;
+    for (const int s : verifiedSymmetryOps<KCluster>(H0))
+      derived.push_back(realized[s]);
+
+    // Restore the declared state, and verify the restoration: every consumer of the symmetry
+    // domains must see exactly the state the declared initializer produced.
+    detail::installPointGroup<KCluster, typename Lattice::DCA_point_group>();
+
+    if (KSymDmn::parameter_type::get_size() != n_declared)
+      throw std::logic_error(
+          "deriveAndComparePointGroup: failed to restore the declared symmetry state of " +
+          KCluster::get_name() + ".");
+    for (std::size_t i = 0; i < declared.size(); ++i)
+      if (!detail::isSameSymmetryElement(realized[i], declared[i]))
+        throw std::logic_error(
+            "deriveAndComparePointGroup: failed to restore the declared symmetry state of " +
+            KCluster::get_name() + ".");
+
+    // Compare as sets and report the divergence, if any, naming the ops that differ.
+    std::vector<std::string> underdeclared;
+    for (const auto& op : derived)
+      if (!detail::containsSymmetryElement(declared, op))
+        underdeclared.push_back(op.describe());
+
+    std::vector<std::string> unverified;
+    for (const auto& op : declared)
+      if (!detail::containsSymmetryElement(derived, op))
+        unverified.push_back(op.describe());
+
+    if (underdeclared.empty() && unverified.empty())
+      return "";
+
+    auto join = [](const std::vector<std::string>& ops) {
+      std::string s;
+      for (std::size_t i = 0; i < ops.size(); ++i)
+        s += (i ? ", " : "") + ops[i];
+      return s;
+    };
+
+    std::ostringstream report;
+    report << "derived-symmetry check [" << KCluster::get_name() << "]:\n"
+           << "\tdeclared group: " << declared.size() << " op(s); derived group: " << derived.size()
+           << " op(s).\n";
+    if (!underdeclared.empty())
+      report << "\tunder-declared: " << underdeclared.size()
+             << " derived H0-verified symmetry op(s) beyond the declared group "
+             << "(symmetry averaging left unexploited): " << join(underdeclared) << ".\n";
+    if (!unverified.empty())
+      report << "\tWARNING: " << unverified.size()
+             << " declared op(s) not reproduced by the derivation -- the op fails the "
+             << "H0-invariance check (or, the lattice is not in the pool's standard orientation): "
+             << join(unverified) << ".\n";
+    return report.str();
+  }
+}
+
+}  // namespace phys
+}  // namespace dca
+
+#endif  // DCA_PHYS_DCA_STEP_SYMMETRIZATION_DERIVE_POINT_GROUP_HPP

include/dca/phys/dca_step/symmetrization/solve_orbital_op_signs.hpp

@@ -29,9 +29,11 @@
 //
 // 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
+// enumerate them and keep the one that satisfies every coupling constraint. Four conditions make
 // an operation fail:
 //
+//   * a missing band image in the symmetry table (the -1 value recorded by set_symmetry_matrices
+//     when it finds no admissible image) means the op has no candidate permutation at all;
 //   * 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
@@ -78,8 +80,18 @@ void solveSignsForOp(int s, int nb, int nk, const SymFunc& sym, const FoldFunc&
   // 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<int> image(nb);
-  for (int b = 0; b < nb; ++b)
+  for (int b = 0; b < nb; ++b) {
     image[b] = sym(0, b, s).second;
+    // set_symmetry_matrices records a silent -1 sentinel when its position/flavor matching finds
+    // no admissible image (its own throw is commented out), so the table cannot be assumed
+    // complete here. Reject the op before the sentinel is used to index H0.
+    if (image[b] < 0 || image[b] >= nb)
+      throw std::out_of_range(
+          "solveOrbitalOpSignsFromH0: no band image in the symmetry table for band " +
+          std::to_string(b) + " under op " + std::to_string(s) +
+          " -- the geometric position/flavor matching found no admissible image, so the op has no "
+          "candidate permutation.");
+  }
 
   // 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