CASM  1.1.0
A Clusters Approach to Statistical Mechanics
IrrepWedge.cc
Go to the documentation of this file.
2 
5 #include "casm/misc/algorithm.hh"
8 
9 namespace CASM {
10 
11 namespace SymRepTools_v2 {
12 
13 namespace IrrepWedgeImpl {
14 
17  MatrixRep const &_rep,
18  GroupIndices const &head_group) {
19  Eigen::MatrixXd t_axes = irrep.trans_mat.transpose().real();
21  Eigen::VectorXd v = axes.col(0);
22  Eigen::VectorXd vbest;
23 
24  axes.setZero(irrep.vector_dim(), irrep.irrep_dim());
25 
26  axes.col(0) = v;
27 
28  for (Index i = 1; i < axes.cols(); ++i) {
29  double bestproj = -1;
30  for (Index element_index : head_group) {
31  v = _rep[element_index] * axes.col(0);
32  // std::cout << "v: " << v.transpose() << std::endl;
33  bool skip_op = false;
34  for (Index j = 0; j < i; ++j) {
35  if (almost_equal(v, axes.col(j))) {
36  skip_op = true;
37  break;
38  }
39  }
40  if (skip_op) continue;
41 
42  // std::cout << "bproj: " << bestproj << " proj: " <<
43  // (v.transpose()*axes).transpose() << std::endl;
44  if (bestproj < (v.transpose() * axes).sum()) {
45  bestproj = (v.transpose() * axes).sum();
46  vbest = v;
47  }
48  }
49  axes.col(i) = vbest;
50  }
51  return IrrepWedge(irrep, axes);
52 }
53 
54 } // namespace IrrepWedgeImpl
55 
57  : irrep_info(std::move(_irrep_info)), axes(std::move(_axes)) {}
58 
68  IrrepWedge irrep_wedge(make_dummy_irrep_info(axes), axes);
69  irrep_wedge.mult.reserve(axes.cols());
70  for (Index i = 0; i < axes.cols(); ++i) {
71  irrep_wedge.mult.push_back(1);
72  }
73  return irrep_wedge;
74 }
75 
76 SubWedge::SubWedge(std::vector<IrrepWedge> const &_iwedges)
77  : m_iwedges(_iwedges), m_trans_mat(_subwedge_to_trans_mat(m_iwedges)) {}
78 
80  std::vector<IrrepWedge> const &_iwedges) {
81  if (_iwedges.empty()) return Eigen::MatrixXd();
82  Eigen::MatrixXd result(_iwedges[0].axes.rows(), _iwedges[0].axes.rows());
83  Index i = 0;
84  for (Index w = 0; w < _iwedges.size(); ++w) {
85  result.block(0, i, _iwedges[w].axes.rows(), _iwedges[w].axes.cols()) =
86  _iwedges[w].axes;
87  i += _iwedges[w].axes.cols();
88  }
89  if (i < result.cols()) result.conservativeResize(Eigen::NoChange, i);
90  return result;
91 }
92 
102  return SubWedge({make_dummy_irrep_wedge(axes)});
103 }
104 
109 std::vector<IrrepWedge> make_irrep_wedges(
111  std::vector<IrrepInfo> const &irreps = irrep_decomposition.irreps;
112  MatrixRep const &fullspace_rep = irrep_decomposition.fullspace_rep;
113  GroupIndices const &head_group = irrep_decomposition.head_group;
114 
115  std::vector<IrrepWedge> wedges;
116  wedges.reserve(irreps.size());
117  double best_proj, tproj;
118 
119  // std::cout << "subspace: \n" << _subspace << std::endl;
120  // std::cout << "irrep decomposition size: " << irreps.size() << std::endl;
121 
122  for (IrrepInfo const &irrep : irreps) {
123  // 1D irreps directions can have positive and negative directions, but we
124  // only want to include one. only 1D irreps can have singly-degenerate
125  // directions and two singly-degenerate directions indicate the same vector
126  // duplicated in positive and negative direction (because they are not
127  // equivalent by symmetry) If irrep.directions[0] is singly degenerate
128  // (orbits size == 1) then irrepdim is 1 and we only need one direction to
129  // define wedge
130  wedges.emplace_back(
131  irrep, Eigen::MatrixXd::Zero(irrep.vector_dim(), irrep.irrep_dim()));
132  // std::cout << "Irrep characters: \n" << irrep.characters << std::endl;
133  // std::cout << "Irrep directions: " << irrep.directions.size() <<
134  // std::endl;
135  if (irrep.directions.empty()) {
137  irrep, fullspace_rep, head_group);
138  continue;
139  }
140 
141  // std::cout << "Irrep direction orbit" << 0 << " : " <<
142  // irrep.directions[0].size() << std::endl; std::cout << "Irrep direction: "
143  // << irrep.directions[0][0].transpose() << std::endl;
144  wedges.back().axes.col(0) = irrep.directions[0][0];
145  wedges.back().mult.push_back(irrep.directions[0].size());
146  for (Index i = 1; i < irrep.irrep_dim(); i++) {
147  // std::cout << "Irrep direction orbit" << i << " : " <<
148  // irrep.directions[i].size() << std::endl; std::cout << "Irrep direction:
149  // " << irrep.directions[i][0].transpose() << std::endl;
150  Index j_best = 0;
151  best_proj =
152  (wedges.back().axes.transpose() * irrep.directions[i][0]).sum();
153  for (Index j = 1; j < irrep.directions[i].size(); j++) {
154  tproj = (wedges.back().axes.transpose() * irrep.directions[i][j]).sum();
155  if (tproj > best_proj) {
156  best_proj = tproj;
157  j_best = j;
158  }
159  }
160 
161  wedges.back().axes.col(i) = irrep.directions[i][j_best];
162  wedges.back().mult.push_back(irrep.directions[i].size());
163  }
164  // std::cout << "New irrep wedge: \n" << wedges.back().axes.transpose() <<
165  // std::endl;
166  }
167  return wedges;
168 }
169 
172 std::vector<SubWedge> make_symrep_subwedges(
174  auto irrep_wedge_compare = [](const IrrepWedge &a,
175  const IrrepWedge &b) -> bool {
176  return Eigen::almost_equal(a.axes, b.axes);
177  };
178 
179  auto tot_wedge_compare = [irrep_wedge_compare](
180  const std::vector<IrrepWedge> &a,
181  const std::vector<IrrepWedge> &b) -> bool {
182  if (a.size() != b.size()) return false;
183  for (auto ita = a.begin(), itb = b.begin(); ita != a.end(); ++ita, ++itb)
184  if (!irrep_wedge_compare(*ita, *itb)) return false;
185  // std::cout << "Equal: \n" << SubWedge(a).trans_mat() << ",\n" <<
186  // SubWedge(b).trans_mat() << "\n\n";
187  return true;
188  };
189 
190  std::vector<IrrepWedge> init_wedges = make_irrep_wedges(irrep_decomposition);
191  MatrixRep const &fullspace_rep = irrep_decomposition.fullspace_rep;
192  GroupIndices const &head_group = irrep_decomposition.head_group;
193 
194  std::vector<SubWedge> result;
195 
196  // irrep_wedge_orbits[w] is orbit of wedges[w]
197  multivector<IrrepWedge>::X<2> irrep_wedge_orbits;
198  irrep_wedge_orbits.reserve(init_wedges.size());
199  // max_equiv[w] is irrep_wedge_orbits[w].size()-1
200  std::vector<Index> max_equiv;
201  max_equiv.reserve(init_wedges.size());
202  // std::cout << "irreducible wedges for group of order " << head_group.size()
203  // << std::endl;
204  Index imax = 0;
205  multivector<Index>::X<2> subgroups;
206  for (IrrepWedge const &wedge : init_wedges) {
207  // std::cout << "Working wedge with axes: \n" << wedge.axes.transpose() <<
208  // std::endl;
209  irrep_wedge_orbits.push_back({wedge});
210 
211  // Start getting orbit of wedges[w]
212  subgroups.push_back({});
213  for (Index element_index : head_group) {
214  IrrepWedge test_wedge{wedge};
215  test_wedge.axes = fullspace_rep[element_index] * wedge.axes;
216  Index o = 0;
217  for (; o < irrep_wedge_orbits.back().size(); ++o) {
218  if (irrep_wedge_compare(irrep_wedge_orbits.back()[o], test_wedge)) {
219  if (o == 0) {
220  subgroups.back().push_back(element_index);
221  }
222  break;
223  }
224  }
225  if (o < irrep_wedge_orbits.back().size()) continue;
226  irrep_wedge_orbits.back().push_back(test_wedge);
227  }
228  // std::cout << "wedge mult: " << irrep_wedge_orbits.back().size();
229  // std::cout << "; subgroups[" << subgroups.size() << "]: " <<
230  // subgroups.back() << std::endl; std::cout << "N equiv wedges found: " <<
231  // irrep_wedge_orbits.back().size() << std::endl;
232  max_equiv.push_back(irrep_wedge_orbits.back().size() - 1);
233  if (max_equiv.back() > max_equiv[imax]) imax = max_equiv.size() - 1;
234  }
235  max_equiv[imax] = 0;
236 
237  // Counter over combinations of equivalent wedges
238  Counter<std::vector<Index>> wcount(std::vector<Index>(init_wedges.size(), 0),
239  max_equiv,
240  std::vector<Index>(init_wedges.size(), 1));
241 
242  // std::cout << "max_equiv: " << max_equiv << std::endl;
243  // std::cout << "init wcount: " << wcount() << std::endl;
244  multivector<IrrepWedge>::X<3> tot_wedge_orbits;
245  // std::cout << "Starting slow bit!\n";
246  for (; wcount; ++wcount) {
247  std::vector<IrrepWedge> twedge = init_wedges;
248  for (Index i = 0; i < init_wedges.size(); i++)
249  twedge[i].axes = irrep_wedge_orbits[i][wcount[i]].axes;
250 
251  if (contains_if(
252  tot_wedge_orbits,
253  [tot_wedge_compare, &twedge](
254  const multivector<IrrepWedge>::X<2> &wedge_orbit) -> bool {
255  return contains(wedge_orbit, twedge, tot_wedge_compare);
256  }))
257  continue;
258 
259  tot_wedge_orbits.push_back({twedge});
260  result.emplace_back(twedge);
261  for (Index p : subgroups[imax]) {
262  for (Index i = 0; i < twedge.size(); i++)
263  twedge[i].axes =
264  fullspace_rep[p] * result.back().irrep_wedges()[i].axes;
265  if (!contains(tot_wedge_orbits.back(), twedge, tot_wedge_compare)) {
266  tot_wedge_orbits.back().push_back(twedge);
267  }
268  }
269  }
270  // std::cout << "Num subwedges: " << result.size() << std::endl;
271  return result;
272 }
273 
274 } // namespace SymRepTools_v2
275 
276 } // namespace CASM
A Counter allows looping over many incrementing variables in one loop.
Definition: Counter.hh:109
SubWedge(std::vector< IrrepWedge > const &_iwedges)
Definition: IrrepWedge.cc:76
static Eigen::MatrixXd _subwedge_to_trans_mat(std::vector< IrrepWedge > const &_iwedges)
Definition: IrrepWedge.cc:79
static IrrepWedge _wedge_from_pseudo_irrep(IrrepInfo const &irrep, MatrixRep const &_rep, GroupIndices const &head_group)
Definition: IrrepWedge.cc:16
IrrepInfo make_dummy_irrep_info(Eigen::MatrixXcd const &trans_mat)
Construct a "dummy" IrrepInfo with user specified transformtion matrix.
std::set< Index > GroupIndices
IrrepWedge make_dummy_irrep_wedge(Eigen::MatrixXd const &axes)
Construct a "dummy" IrrepWedge with user specified axes.
Definition: IrrepWedge.cc:67
std::vector< IrrepWedge > make_irrep_wedges(IrrepDecomposition const &irrep_decomposition)
Make IrrepWedges from an IrrepDecomposition.
Definition: IrrepWedge.cc:109
std::vector< SubWedge > make_symrep_subwedges(IrrepDecomposition const &irrep_decomposition)
Find full irreducible wedge of a group-represented vector space, as a vector of SubWedges,...
Definition: IrrepWedge.cc:172
SubWedge make_dummy_subwedge(Eigen::MatrixXd const &axes)
Makes a "dummy" SubWedge from a single "dummy" IrrepWedge with given axes.
Definition: IrrepWedge.cc:101
Eigen::MatrixXcd vector_space_prepare(Eigen::MatrixXcd const &vector_space, double tol)
Vector space preparation for comparison.
std::vector< Eigen::MatrixXd > MatrixRep
Main CASM namespace.
Definition: APICommand.hh:8
bool almost_equal(ClusterInvariants const &A, ClusterInvariants const &B, double tol)
Check if ClusterInvariants are equal.
Eigen::MatrixXd MatrixXd
std::vector< SymRepTools::IrrepInfo > irrep_decomposition(SymGroupRep const &_rep, SymGroup const &head_group, bool allow_complex)
Finds irreducible subspaces that comprise an underlying subspace It does not rely on the character ta...
Definition: SymRepTools.cc:761
const double TOL
Definition: definitions.hh:30
Container::value_type sum(const Container &container, typename Container::value_type init_val=0)
Definition: algorithm.hh:131
bool contains(const Container &container, const T &value)
Equivalent to container.end() != std::find(container.begin(), container.end(), value)
Definition: algorithm.hh:83
bool contains_if(const Container &container, UnaryPredicate p)
Equivalent to container.end() != std::find_if(container.begin(), container.end(), p)
Definition: algorithm.hh:98
INDEX_TYPE Index
For long integer indexing:
Definition: definitions.hh:39
Eigen::VectorXd VectorXd
bool almost_equal(const Eigen::MatrixBase< Derived1 > &val1, const Eigen::MatrixBase< Derived2 > &val2, double tol=CASM::TOL)
Definition: stream_io.hh:24
Performs irreducible subspace construction and symmetrization.
Index irrep_dim() const
Dimension of irreducible vector space (less than or equal to vector_dim())
An irreducible wedge in an irreducible vector space.
Definition: IrrepWedge.hh:14
IrrepWedge(IrrepInfo _irrep_info, Eigen::MatrixXd _axes)
Definition: IrrepWedge.cc:56
std::vector< Index > mult
Symmetric multiplicity of each column of 'axes'.
Definition: IrrepWedge.hh:27
Shortcut for multidimensional vector (std::vector< std::vector< ...)
Definition: multivector.hh:26
typename multivector_impl::multivector_tmp< T, N >::type X
Definition: multivector.hh:28