Major changes to the mapping routine that should ensure that the

symmetrized costs are now calculated with the parent symmetry
operations. Added a few routines to SimpleStructureTools as well.

.gitignore

@@ -96,3 +96,5 @@ stamp-h1
 .clang
 .nostyle.txt
 doc/doxygen_graph
+
+compile_commands.json

include/casm/crystallography/BasicStructureTools.hh

@@ -3,12 +3,10 @@
 
 #include <vector>
 
+#include "casm/crystallography/Lattice.hh"
 #include "casm/external/Eigen/Core"
 #include "casm/global/definitions.hh"
 #include "casm/global/eigen.hh"
-// #include "submodules/eigen/Eigen/src/Core/PermutationMatrix.h"
-// #include "submodules/eigen/Eigen/src/Core/util/Constants.h"
-#include <vector>
 
 namespace CASM {
 namespace xtal {
@@ -32,6 +30,15 @@ bool is_primitive(const BasicStructure &struc, double tol = TOL);
 BasicStructure make_primitive(const BasicStructure &non_primitive_struc,
                               double tol = TOL);
 
+/// Calculates the rotation angle and axis of a symmetry operation. This
+/// function is almost exactly identical to the constructor of SymInfo::SymInfo
+std::pair<double, Eigen::Vector3d> calc_rotation_angle_and_axis(
+    const SymOp &op, const Lattice &lat);
+
+/// Sort a factor group based on a lexicographical comparison of (-det, -trace,
+/// angle, axis, tau)
+void sort_factor_group(std::vector<SymOp> &factor_group, const Lattice &lat);
+
 /// Create the factor group of the given structure. If the structure has no
 /// degrees of freedom affected by time reversal, time reversal is ignored.
 /// Otherwise symmetry operations are checked for time reversal

include/casm/crystallography/LatticeMap.hh

@@ -46,10 +46,8 @@ class StrainCostCalculator {
 
   //\brief Symmetrized strain cost; Utilizes the parent point group symmetry to
   // calculate only the symmetry breaking lattice cost
-  double strain_cost(
-      Eigen::Matrix3d const &_deformation_gradient,
-      Eigen::Matrix3d const &parent_lattice,
-      std::vector<Eigen::Matrix3i> const &parent_fsym_mats) const;
+  double strain_cost(Eigen::Matrix3d const &_deformation_gradient,
+                     SymOpVector const &parent_point_group) const;
 
  private:
   Eigen::MatrixXd m_gram_mat;
@@ -144,6 +142,9 @@ class LatticeMap {
   // parent point group matrices, in fractional coordinates
   std::vector<Eigen::Matrix3i> m_parent_fsym_mats;
 
+  // parent point group in cartesian coordinates
+  SymOpVector m_parent_point_group;
+
   // child point group matrices, in fractional coordinates
   std::vector<Eigen::Matrix3i> m_child_fsym_mats;
 

include/casm/crystallography/SimpleStructureTools.hh

@@ -6,6 +6,7 @@
 #include <vector>
 
 #include "casm/crystallography/DoFDecl.hh"
+#include "casm/crystallography/SymType.hh"
 #include "casm/external/Eigen/Dense"
 #include "casm/global/definitions.hh"
 
@@ -31,19 +32,25 @@ class BasicStructure;
 SimpleStructure make_superstructure(Eigen::Ref<const Eigen::Matrix3i> const &_T,
                                     SimpleStructure const &_sstruc);
 
+/// \brief Constructs a vector containing the basis index of the ith site in the
+/// supercell
+std::vector<Index> superstructure_basis_idx(
+    Eigen::Ref<const Eigen::Matrix3i> const &_T,
+    SimpleStructure const &_sstruc);
+
 /// \brief Construct from decorated structure
 SimpleStructure make_simple_structure(BasicStructure const &_struc);
 
 /// \brief Determine which sites of a BasicStructure can host each atom of a
 /// SimpleStructure result[i] is set of site indices in @param _prim that can
 /// host atom 'i' of @param sstruc
-std::vector<std::set<Index> > atom_site_compatibility(
+std::vector<std::set<Index>> atom_site_compatibility(
     SimpleStructure const &sstruc, BasicStructure const &_prim);
 
 /// \brief Determine which sites of a BasicStructure can host each molecule of a
 /// SimpleStructure result[i] is set of site indices in @param _prim that can
 /// host molecule 'i' of @param sstruc
-std::vector<std::set<Index> > mol_site_compatibility(
+std::vector<std::set<Index>> mol_site_compatibility(
     SimpleStructure const &sstruc, BasicStructure const &_prim);
 
 /// \brief use information in _reference to initialize atom_info from mol_info
@@ -62,7 +69,12 @@ void _atomize(SimpleStructure &_sstruc,
 BasicStructure make_basic_structure(
     SimpleStructure const &_sstruc, std::vector<DoFKey> const &_all_dofs,
     SimpleStructureTools::SpeciesMode mode,
-    std::vector<std::vector<Molecule> > _allowed_occupants = {});
+    std::vector<std::vector<Molecule>> _allowed_occupants = {});
+
+std::vector<Eigen::MatrixXd> generate_invariant_shuffle_modes(
+    const std::vector<xtal::SymOp> &factor_group,
+    const std::vector<Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic,
+                                               Index>> &permute_group);
 /** @} */
 }  // namespace xtal
 

include/casm/crystallography/StrucMapCalculatorInterface.hh

@@ -80,6 +80,16 @@ class StrucMapCalculatorInterface {
     return m_fixed_species;
   }
 
+  /// \brief Sets the sym_invariant_modes
+  void set_sym_invariant_displacement_modes(
+      const std::vector<Eigen::MatrixXd> &_sym_invariant_displacement_modes) {
+    m_sym_invariant_displacement_modes = _sym_invariant_displacement_modes;
+  }
+
+  const std::vector<Eigen::MatrixXd> &sym_invariant_displacement_modes() const {
+    return m_sym_invariant_displacement_modes;
+  };
+
   /// \brief Return maximum possible number of vacancies in underlying primitive
   /// structure
   Index max_n_va() const { return va_allowed().size(); }
@@ -227,6 +237,10 @@ class StrucMapCalculatorInterface {
 
   std::unordered_set<Index> m_va_allowed;
 
+  /// \brief vector of symmetry invariant displacement modes in the parent
+  /// structure
+  std::vector<Eigen::MatrixXd> m_sym_invariant_displacement_modes;
+
   /// \brief Make an exact copy of the calculator (including any initialized
   /// members)
   virtual StrucMapCalculatorInterface *_clone() const = 0;

include/casm/crystallography/StrucMapping.hh

@@ -5,6 +5,8 @@
 #include <vector>
 
 #include "casm/crystallography/SimpleStructure.hh"
+#include "casm/crystallography/SimpleStructureTools.hh"
+#include "casm/crystallography/StrucMapCalculatorInterface.hh"
 #include "casm/crystallography/Superlattice.hh"
 #include "casm/crystallography/SymType.hh"
 #include "casm/external/Eigen/Core"
@@ -509,8 +511,13 @@ class StrucMapper {
   /// when performing the atomic maps. This cost only accounts for that part of
   /// the displacement field that breaks the symmetry of the parent crystal
   /// structure
-  void set_symmetrize_atomic_cost(bool _sym_atomic_cost) {
+  void set_symmetrize_atomic_cost(
+      bool _sym_atomic_cost, const SymOpVector &factor_group,
+      const std::vector<Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic,
+                                                 Index>> &permutation_group) {
     m_symmetrize_atomic_cost = _sym_atomic_cost;
+    m_calc_ptr->set_sym_invariant_displacement_modes(
+        generate_invariant_shuffle_modes(factor_group, permutation_group));
   }
 
   bool symmetrize_lattice_cost() const { return m_symmetrize_lattice_cost; }

src/casm/crystallography/BasicStructureTools.cc

@@ -2,6 +2,7 @@
 
 #include <algorithm>
 #include <atomic>
+#include <tuple>
 #include <unordered_set>
 #include <utility>
 #include <vector>
@@ -12,6 +13,7 @@
 #include "casm/crystallography/Coordinate.hh"
 #include "casm/crystallography/DoFSet.hh"
 #include "casm/crystallography/IntegralCoordinateWithin.hh"
+#include "casm/crystallography/Lattice.hh"
 #include "casm/crystallography/Niggli.hh"
 #include "casm/crystallography/Site.hh"
 #include "casm/crystallography/Superlattice.hh"
@@ -315,6 +317,151 @@ BasicStructure make_primitive(const BasicStructure &non_primitive_struc,
   return primitive_struc;
 }
 
+/// \brief Calculates the rotation angle and axis of a symmetry operation. This
+/// function is almost exactly identical to the constructor of SymInfo::SymInfo
+std::pair<double, Eigen::Vector3d> calc_rotation_angle_and_axis(
+    const SymOp &op, const Lattice &lat) {
+  auto matrix = op.matrix;
+  double angle;
+  Eigen::Vector3d rotation_axis, _axis;
+
+  // Simplest case is identity: has no axis and no location
+  if (almost_equal(matrix.trace(), 3.) || almost_equal(matrix.trace(), -3.)) {
+    angle = 0;
+    _axis = Eigen::Vector3d::Zero();
+    rotation_axis = Coordinate(_axis, lat, CART).const_frac();
+    return std::make_pair(angle, rotation_axis);
+  }
+
+  // det is -1 if improper and +1 if proper
+  int det = round(matrix.determinant());
+
+  // Find eigen decomposition of proper operation (by multiplying by
+  // determinant)
+  Eigen::EigenSolver<Eigen::Matrix3d> t_eig(det * matrix);
+
+  // 'axis' is eigenvector whose eigenvalue is +1
+  for (Index i = 0; i < 3; i++) {
+    if (almost_equal(t_eig.eigenvalues()(i), std::complex<double>(1, 0))) {
+      _axis = t_eig.eigenvectors().col(i).real();
+      break;
+    }
+  }
+
+  // Sign convention for 'axis': first non-zero element is positive
+  for (Index i = 0; i < 3; i++) {
+    if (!almost_zero(_axis[i])) {
+      _axis *= float_sgn(_axis[i]);
+      break;
+    }
+  }
+
+  // get vector orthogonal to axis: ortho,
+  // apply matrix: rot
+  // and check angle between ortho and det*rot,
+  // using determinant to get the correct angle for improper
+  // (i.e. want angle before inversion for rotoinversion)
+  Eigen::Vector3d ortho = _axis.unitOrthogonal();
+  Eigen::Vector3d rot = det * (matrix * ortho);
+  angle = fmod(
+      (180. / M_PI) * atan2(_axis.dot(ortho.cross(rot)), ortho.dot(rot)) + 360.,
+      360.);
+  rotation_axis = Coordinate(_axis, lat, CART).const_frac();
+  return std::make_pair(angle, rotation_axis);
+}
+
+//*******************************************************************************************
+/// \brief Sort SymOp in the SymGroup
+///
+/// SymOp are sorted by lexicographical comparison of: (-det, -trace, angle,
+/// axis, tau)
+/// - angle is positive
+/// - axis[0] is positive
+///
+/// This function is an almost identical clone of SymGroup::sort()
+void sort_factor_group(std::vector<SymOp> &factor_group, const Lattice &lat) {
+  // floating point comparison tolerance
+  double tol = TOL;
+
+  // COORD_TYPE print_mode = CART;
+
+  // compare on vector of '-det', '-trace', 'angle', 'axis', 'tau'
+  typedef Eigen::Matrix<double, 10, 1> key_type;
+  auto make_key = [](const SymOp &op, const Lattice &lat) {
+    key_type vec;
+    int offset = 0;
+    double sym_angle;
+    Eigen::Vector3d sym_frac;
+    std::tie(sym_angle, sym_frac) = calc_rotation_angle_and_axis(op, lat);
+
+    vec[offset] = double(op.is_time_reversal_active);
+    offset++;
+
+    vec[offset] = -op.matrix.determinant();
+    offset++;
+
+    vec[offset] = -op.matrix.trace();
+    offset++;
+
+    vec[offset] = sym_angle;
+    offset++;
+
+    vec.segment<3>(offset) = sym_frac;
+    offset += 3;
+
+    vec.segment<3>(offset) = Coordinate(op.translation, lat, CART).const_frac();
+    offset += 3;
+
+    return vec;
+  };
+
+  // define symop compare function
+  auto op_compare = [tol](const key_type &A, const key_type &B) {
+    return float_lexicographical_compare(A, B, tol);
+  };
+
+  typedef std::map<key_type, SymOp,
+                   std::reference_wrapper<decltype(op_compare)>>
+      map_type;
+
+  // first put identity in position 0 in order to calculat multi_table correctly
+  for (int i = 0; i < factor_group.size(); ++i) {
+    if (factor_group[i].matrix.isIdentity() &&
+        factor_group[i].translation.norm() < TOL &&
+        !factor_group[i].is_time_reversal_active) {
+      std::swap(factor_group[0], factor_group[i]);
+      break;
+    }
+  }
+
+  // sort conjugacy class using the first symop in the sorted class
+  auto cclass_compare = [tol](const map_type &A, const map_type &B) {
+    return float_lexicographical_compare(A.begin()->first, B.begin()->first,
+                                         tol);
+  };
+
+  // sort elements in each conjugracy class (or just put all elements in the
+  // first map)
+  std::set<map_type, std::reference_wrapper<decltype(cclass_compare)>> sorter(
+      cclass_compare);
+
+  // else just sort element
+  map_type all_op(op_compare);
+  for (auto it = factor_group.begin(); it != factor_group.end(); ++it) {
+    all_op.emplace(make_key(*it, lat), *it);
+  }
+  sorter.emplace(std::move(all_op));
+
+  // copy symop back into group
+  int j = 0;
+  for (auto const &cclass : sorter) {
+    for (auto it = cclass.begin(); it != cclass.end(); ++it) {
+      factor_group[j] = it->second;
+      ++j;
+    }
+  }
+}
+
 std::vector<SymOp> make_factor_group(const BasicStructure &struc, double tol) {
   auto prim_factor_group_pair = ::make_primitive_factor_group(struc, tol);
   const BasicStructure &primitive_struc = prim_factor_group_pair.first;
@@ -348,7 +495,7 @@ std::vector<SymOp> make_factor_group(const BasicStructure &struc, double tol) {
           prim_op);
     }
   }
-
+  sort_factor_group(factor_group, struc.lattice());
   return factor_group;
 }