From bd1e13e18097ed8a6f66f783f4484091b3a714f7 Mon Sep 17 00:00:00 2001 From: Anirudh Date: Mon, 1 Feb 2021 15:52:32 -0800 Subject: [PATCH] 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 | 2 + .../crystallography/BasicStructureTools.hh | 13 +- include/casm/crystallography/LatticeMap.hh | 9 +- .../crystallography/SimpleStructureTools.hh | 18 ++- .../StrucMapCalculatorInterface.hh | 14 ++ include/casm/crystallography/StrucMapping.hh | 9 +- .../crystallography/BasicStructureTools.cc | 149 +++++++++++++++++- src/casm/crystallography/LatticeMap.cc | 30 ++-- .../SimpleStrucMapCalculator.cc | 56 ++++--- .../crystallography/SimpleStructureTools.cc | 112 ++++++++++++- src/casm/crystallography/StrucMapping.cc | 1 - tests/unit/mapping/StrucMapping_test.cpp | 5 +- 12 files changed, 358 insertions(+), 60 deletions(-) diff --git a/.gitignore b/.gitignore index af14c3dda..8bbd54726 100644 --- a/.gitignore +++ b/.gitignore @@ -96,3 +96,5 @@ stamp-h1 .clang .nostyle.txt doc/doxygen_graph + +compile_commands.json \ No newline at end of file diff --git a/include/casm/crystallography/BasicStructureTools.hh b/include/casm/crystallography/BasicStructureTools.hh index 6b49b9936..61b6d7a6d 100644 --- a/include/casm/crystallography/BasicStructureTools.hh +++ b/include/casm/crystallography/BasicStructureTools.hh @@ -3,12 +3,10 @@ #include +#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 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 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 &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 diff --git a/include/casm/crystallography/LatticeMap.hh b/include/casm/crystallography/LatticeMap.hh index 53686f8ca..23a5f23dc 100644 --- a/include/casm/crystallography/LatticeMap.hh +++ b/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 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 m_parent_fsym_mats; + // parent point group in cartesian coordinates + SymOpVector m_parent_point_group; + // child point group matrices, in fractional coordinates std::vector m_child_fsym_mats; diff --git a/include/casm/crystallography/SimpleStructureTools.hh b/include/casm/crystallography/SimpleStructureTools.hh index b8d3eaef3..a545a9ff6 100644 --- a/include/casm/crystallography/SimpleStructureTools.hh +++ b/include/casm/crystallography/SimpleStructureTools.hh @@ -6,6 +6,7 @@ #include #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 &_T, SimpleStructure const &_sstruc); +/// \brief Constructs a vector containing the basis index of the ith site in the +/// supercell +std::vector superstructure_basis_idx( + Eigen::Ref 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 > atom_site_compatibility( +std::vector> 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 > mol_site_compatibility( +std::vector> 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 const &_all_dofs, SimpleStructureTools::SpeciesMode mode, - std::vector > _allowed_occupants = {}); + std::vector> _allowed_occupants = {}); + +std::vector generate_invariant_shuffle_modes( + const std::vector &factor_group, + const std::vector> &permute_group); /** @} */ } // namespace xtal diff --git a/include/casm/crystallography/StrucMapCalculatorInterface.hh b/include/casm/crystallography/StrucMapCalculatorInterface.hh index 8c53ac819..2961d8e21 100644 --- a/include/casm/crystallography/StrucMapCalculatorInterface.hh +++ b/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 &_sym_invariant_displacement_modes) { + m_sym_invariant_displacement_modes = _sym_invariant_displacement_modes; + } + + const std::vector &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 m_va_allowed; + /// \brief vector of symmetry invariant displacement modes in the parent + /// structure + std::vector m_sym_invariant_displacement_modes; + /// \brief Make an exact copy of the calculator (including any initialized /// members) virtual StrucMapCalculatorInterface *_clone() const = 0; diff --git a/include/casm/crystallography/StrucMapping.hh b/include/casm/crystallography/StrucMapping.hh index 21508e705..aba275b7d 100644 --- a/include/casm/crystallography/StrucMapping.hh +++ b/include/casm/crystallography/StrucMapping.hh @@ -5,6 +5,8 @@ #include #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> &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; } diff --git a/src/casm/crystallography/BasicStructureTools.cc b/src/casm/crystallography/BasicStructureTools.cc index 73965ac63..8f50ee431 100644 --- a/src/casm/crystallography/BasicStructureTools.cc +++ b/src/casm/crystallography/BasicStructureTools.cc @@ -2,6 +2,7 @@ #include #include +#include #include #include #include @@ -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 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 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(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 &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 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> + 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> 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 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 make_factor_group(const BasicStructure &struc, double tol) { prim_op); } } - + sort_factor_group(factor_group, struc.lattice()); return factor_group; } diff --git a/src/casm/crystallography/LatticeMap.cc b/src/casm/crystallography/LatticeMap.cc index 12ce356af..c43dd2e13 100644 --- a/src/casm/crystallography/LatticeMap.cc +++ b/src/casm/crystallography/LatticeMap.cc @@ -1,5 +1,7 @@ #include "casm/crystallography/LatticeMap.hh" +#include + #include "casm/crystallography/Lattice.hh" #include "casm/crystallography/LatticeIsEquivalent.hh" #include "casm/crystallography/Strain.hh" @@ -108,28 +110,18 @@ double StrainCostCalculator::strain_cost( } //******************************************************************************************* - double StrainCostCalculator::strain_cost( Eigen::Matrix3d const &_deformation_gradient, - Eigen::Matrix3d const &parent_lattice, - std::vector const &parent_fsym_mats) const { - // Convert the fractional sym ops back into cartesian symops - std::vector parent_sym_mats; - parent_sym_mats.reserve(parent_fsym_mats.size()); - Eigen::Matrix3d tmp_op; - for (auto const &op : parent_fsym_mats) { - tmp_op = parent_lattice * op.cast() * parent_lattice.inverse(); - parent_sym_mats.push_back(tmp_op.transpose()); - } + SymOpVector const &parent_point_group) const { // Apply the sym op to the deformation gradient and average Eigen::Matrix3d stretch_aggregate = Eigen::Matrix3d::Zero(); Eigen::Matrix3d stretch = strain::right_stretch_tensor(_deformation_gradient); - - for (auto const &op : parent_sym_mats) { - stretch_aggregate += op * stretch * op.inverse(); + for (auto const &op : parent_point_group) { + stretch_aggregate += op.matrix * stretch * op.matrix.inverse(); } - stretch_aggregate = stretch_aggregate / double(parent_sym_mats.size()); - return strain_cost(stretch - stretch_aggregate + Eigen::Matrix3d::Identity()); + stretch_aggregate = stretch_aggregate / double(parent_point_group.size()); + return strain_cost(stretch - stretch_aggregate + Eigen::Matrix3d::Identity(), + 1.0); } //******************************************************************************************* @@ -187,6 +179,9 @@ LatticeMap::LatticeMap(const Lattice &_parent, const Lattice &_child, } } + // Store the parent symmetry operations + m_parent_point_group = _parent_point_group; + // Construct fractional symops for child { IsPointGroupOp symcheck(reduced_child); @@ -301,8 +296,7 @@ const LatticeMap &LatticeMap::best_strain_mapping() const { double LatticeMap::calc_strain_cost( const Eigen::Matrix3d &deformation_gradient) const { if (symmetrize_strain_cost()) - return m_calc.strain_cost(deformation_gradient, m_parent, - m_parent_fsym_mats); + return m_calc.strain_cost(deformation_gradient, m_parent_point_group); else return m_calc.strain_cost(deformation_gradient, m_vol_factor); } diff --git a/src/casm/crystallography/SimpleStrucMapCalculator.cc b/src/casm/crystallography/SimpleStrucMapCalculator.cc index 11fd4e04d..f41a729c9 100644 --- a/src/casm/crystallography/SimpleStrucMapCalculator.cc +++ b/src/casm/crystallography/SimpleStrucMapCalculator.cc @@ -22,6 +22,7 @@ static Coordinate _make_superlattice_coordinate( UnitCell ijk = index_to_ijk_f(ijk_ix); return make_superlattice_coordinate(ijk, superlattice); } + } // namespace Local std::vector SimpleStrucMapCalculator::translations( @@ -230,37 +231,46 @@ SimpleStructure SimpleStrucMapCalculator::resolve_setting( void SimpleStrucMapCalculator::finalize( MappingNode &node, SimpleStructure const &_child_struc, bool const &symmetrize_atomic_cost) const { - // std::cout << " In SimpleStrucMapCalculator::finalize" << std::endl; - // std::cout << " symmetrize_atomic_cost : " << symmetrize_atomic_cost - // << std::endl; populate_displacement(node, _child_struc); node.is_valid = this->_assign_molecules(node, _child_struc); - if (!symmetrize_atomic_cost) + if (!symmetrize_atomic_cost || + this->sym_invariant_displacement_modes().size() == 0) node.atomic_node.cost = StrucMapping::atomic_cost(node, this->struc_info(_child_struc).size()); else { - // Construct a BasicStructure based on the atomic and lattice maps - xtal::SimpleStructure super_pstruc = make_superstructure( + // Get the parent site indices for each site in the supercell + auto basis_idx = superstructure_basis_idx( (node.lattice_node).parent.transformation_matrix_to_super().cast(), parent()); - std::vector _names(node.mol_labels.size(), "-1"); - for (Index idx = 0; idx < node.mol_labels.size(); ++idx) - _names[idx] = node.mol_labels[idx].first; - super_pstruc.atom_info.names = _names; - // std::cout << " constructed simple structure : " << std::endl; - // VaspIO::PrintPOSCAR printer(super_pstruc); - // printer.print(std::cout); - xtal::BasicStructure decorated_pstruc = - xtal::make_basic_structure(super_pstruc, std::vector({}), - SimpleStructureTools::SpeciesMode::ATOM); - auto factor_group = xtal::make_factor_group(decorated_pstruc); - - auto permute_group = - xtal::make_permutation_representation(decorated_pstruc, factor_group); - + // Expand the sym invariant basis into this supercell + std::vector supercell_sym_invariant_modes( + this->sym_invariant_displacement_modes().size(), + Eigen::MatrixXd::Zero(3, basis_idx.size())); + + Eigen::Map disp_map(node.atom_displacement.data(), + node.atom_displacement.size()); + Eigen::MatrixXd sym_removed_disp = node.atom_displacement; + for (int sym_idx = 0; + sym_idx < this->sym_invariant_displacement_modes().size(); ++sym_idx) { + for (int scel_idx = 0; scel_idx < basis_idx.size(); ++scel_idx) + supercell_sym_invariant_modes[sym_idx].col(scel_idx) = + this->sym_invariant_displacement_modes()[sym_idx].col( + basis_idx[scel_idx]); + supercell_sym_invariant_modes[sym_idx] = + supercell_sym_invariant_modes[sym_idx] / + supercell_sym_invariant_modes[sym_idx].norm(); + // Project the displacement mode on to this symmetry invariant mode + Eigen::VectorXd unrolled_sym_mode = Eigen::Map( + supercell_sym_invariant_modes[sym_idx].data(), + supercell_sym_invariant_modes[sym_idx].size()); + double projected_comp = disp_map.dot(unrolled_sym_mode); + sym_removed_disp = + sym_removed_disp - + projected_comp * supercell_sym_invariant_modes[sym_idx]; + } + node.atom_displacement = sym_removed_disp; node.atomic_node.cost = - StrucMapping::atomic_cost(node, factor_group, permute_group, - this->struc_info(_child_struc).size()); + StrucMapping::atomic_cost(node, this->struc_info(_child_struc).size()); } node.cost = node.atomic_weight * node.atomic_node.cost + node.lattice_weight * node.lattice_node.cost; diff --git a/src/casm/crystallography/SimpleStructureTools.cc b/src/casm/crystallography/SimpleStructureTools.cc index 73e755b1d..70eb86895 100644 --- a/src/casm/crystallography/SimpleStructureTools.cc +++ b/src/casm/crystallography/SimpleStructureTools.cc @@ -73,6 +73,27 @@ SimpleStructure make_superstructure(Eigen::Ref const &_T, return superstructure; } +//*************************************************************************************************** + +/// Calculates the parent basis index of each site in a supercell that is +/// generated with the make_superstructure method +/// @param _T is the transformation matrix linking `_sstruc` and the supercell +/// @param _sstruc is a SimpleStructure +/// @return std::vector with the Index in the ith entry corresponding to +/// the index of site i in _sstruc +std::vector superstructure_basis_idx( + Eigen::Ref const &_T, + SimpleStructure const &_sstruc) { + auto all_lattice_points = make_lattice_points(_T.cast()); + Index Nvol = all_lattice_points.size(); + std::vector basis_idx(_sstruc.atom_info.size() * Nvol, -1); + for (Index grid_idx = 0; grid_idx < Nvol; ++grid_idx) { + for (Index atom_idx = 0; atom_idx < _sstruc.atom_info.size(); ++atom_idx) + basis_idx[grid_idx + atom_idx * Nvol] = atom_idx; + } + return basis_idx; +} + //*************************************************************************** SimpleStructure make_simple_structure(BasicStructure const &_struc) { @@ -97,7 +118,7 @@ SimpleStructure make_simple_structure(BasicStructure const &_struc) { BasicStructure make_basic_structure( SimpleStructure const &_sstruc, std::vector const &_all_dofs, SimpleStructure::SpeciesMode mode, - std::vector > _allowed_occupants) { + std::vector> _allowed_occupants) { std::map global_dof; std::map local_dof; for (DoFKey const &dof : _all_dofs) { @@ -153,6 +174,87 @@ BasicStructure make_basic_structure( //*************************************************************************** +/** +** Calculates the invariant shuffle modes in the primitive unit cell. They +*symmetry preserving distortions are found by applying the Reynolds operator to +*the basis of displacement vectors. The average of the resulting basis vectors +*is used to form an orthonormal basis. +* +* @param factor_group use make_factor_group(struc) to obtain this group +* @param permute_group use make_permute_group(struc,factor_group) to obtain this +*group +* +* @return the vector of shuffle modes that are invariant under symmetry. Each +* element has a size `3 x basis_size`. +* +*/ +std::vector generate_invariant_shuffle_modes( + const std::vector &factor_group, + const std::vector> &permute_group) { + if (factor_group.size() != permute_group.size()) { + throw std::runtime_error( + "error, the size of the symmetry operations in " + "generate_invariant_shuffle_modes do not match"); + } + int struc_basis_size = permute_group[0].indices().size(); + // Generate a basis consisting of individual shuffles of each atom in the + // structure. + std::vector displacement_basis; + for (int basis_idx = 0; basis_idx < struc_basis_size; ++basis_idx) { + for (int dir_idx = 0; dir_idx < 3; ++dir_idx) { + Eigen::MatrixXd single_shuffle = + Eigen::MatrixXd::Zero(3, struc_basis_size); + single_shuffle(dir_idx, basis_idx) = 1.0; + displacement_basis.push_back(single_shuffle); + } + } + std::vector displacement_aggregate( + displacement_basis.size(), + Eigen::VectorXd::Zero(displacement_basis[0].cols() * + displacement_basis[0].rows())); + + for (int idx = 0; idx < factor_group.size(); ++idx) { + for (int disp_basis_idx = 0; disp_basis_idx < displacement_basis.size(); + ++disp_basis_idx) { + Eigen::MatrixXd transformed_disp = factor_group[idx].matrix * + displacement_basis[disp_basis_idx] * + permute_group[idx]; + displacement_aggregate[disp_basis_idx] += + Eigen::VectorXd(Eigen::Map( + transformed_disp.data(), + transformed_disp.cols() * transformed_disp.rows())); + } + } + Eigen::MatrixXd sym_disp_basis = Eigen::MatrixXd::Zero( + displacement_aggregate.size(), displacement_aggregate[0].size()); + for (int disp_basis_idx = 0; disp_basis_idx < displacement_basis.size(); + ++disp_basis_idx) { + displacement_aggregate[disp_basis_idx] = + displacement_aggregate[disp_basis_idx] / double(factor_group.size()); + sym_disp_basis.row(disp_basis_idx) = displacement_aggregate[disp_basis_idx]; + } + + // Perform a SVD of the resulting aggregate matrix to obtain the rank and + // space of the symmetry invariant basis vectors + sym_disp_basis.transposeInPlace(); + Eigen::JacobiSVD svd( + sym_disp_basis, Eigen::ComputeThinU | Eigen::ComputeThinV); + int matrix_rank = (svd.singularValues().array().abs() >= CASM::TOL).sum(); + Eigen::MatrixXd sym_preserving_mode_matrix = + svd.matrixV()(Eigen::all, Eigen::seq(0, matrix_rank - 1)); + std::vector sym_preserving_modes; + + for (int sym_mode_idx = 0; sym_mode_idx < matrix_rank; ++sym_mode_idx) { + Eigen::Map _tmp_mode( + sym_preserving_mode_matrix.col(sym_mode_idx).data(), 3, + struc_basis_size); + sym_preserving_modes.push_back(_tmp_mode); + } + return sym_preserving_modes; +} + +//*************************************************************************** void _atomize(SimpleStructure &_sstruc, Eigen::Ref const &_mol_occ, BasicStructure const &_reference) { @@ -228,9 +330,9 @@ void _atomize(SimpleStructure &_sstruc, //*************************************************************************** -std::vector > mol_site_compatibility( +std::vector> mol_site_compatibility( SimpleStructure const &sstruc, BasicStructure const &_prim) { - std::vector > result; + std::vector> result; result.reserve(sstruc.mol_info.names.size()); for (std::string const &sp : sstruc.mol_info.names) { result.push_back({}); @@ -243,9 +345,9 @@ std::vector > mol_site_compatibility( return result; } -std::vector > atom_site_compatibility( +std::vector> atom_site_compatibility( SimpleStructure const &sstruc, BasicStructure const &_prim) { - std::vector > result; + std::vector> result; result.reserve(sstruc.atom_info.names.size()); for (std::string const &sp : sstruc.atom_info.names) { result.push_back({}); diff --git a/src/casm/crystallography/StrucMapping.cc b/src/casm/crystallography/StrucMapping.cc index fbed0bea9..802efbcc8 100644 --- a/src/casm/crystallography/StrucMapping.cc +++ b/src/casm/crystallography/StrucMapping.cc @@ -863,7 +863,6 @@ std::set StrucMapper::_seed_k_best_from_super_lats( if (k == 0 || !valid_index(k)) { max_lattice_cost = min_lattice_cost; } - for (Lattice const &c_lat : _child_scels) { auto pg_indices = invariant_subgroup_indices(c_lat, child_factor_group); SymOpVector c_lat_factor_group; diff --git a/tests/unit/mapping/StrucMapping_test.cpp b/tests/unit/mapping/StrucMapping_test.cpp index 5f4c058a2..7003eec9d 100644 --- a/tests/unit/mapping/StrucMapping_test.cpp +++ b/tests/unit/mapping/StrucMapping_test.cpp @@ -436,6 +436,8 @@ TEST(SymInvariantMappingTest, shuffle) { // Get the factor group of the structure auto parent_fg = xtal::make_factor_group(basic_struc, 0.0001); + auto parent_permute_group = + xtal::make_permutation_representation(basic_struc, parent_fg); EXPECT_EQ(parent_fg.size(), 8) << "Expected the factor group size to be 8, instead " << parent_fg.size() << " operations were found" << std::endl; @@ -455,7 +457,8 @@ TEST(SymInvariantMappingTest, shuffle) { allowed_molecule_names(basic_struc)), lattice_weight, max_vol_change, options, cost_tol, min_va_frac, max_va_frac); - sym_struc_map.set_symmetrize_atomic_cost(true); + sym_struc_map.set_symmetrize_atomic_cost(true, parent_fg, + parent_permute_group); auto child_fg = xtal::make_factor_group(basic_struc);