StrucMapping.hh

Go to the documentation of this file.

#ifndef CASM_StrucMapping
#define CASM_StrucMapping
 
#include <unordered_set>
#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"
#include "casm/global/definitions.hh"
#include "casm/misc/CASM_Eigen_math.hh"
#include "casm/misc/CASM_math.hh"
#include "casm/misc/cloneable_ptr.hh"
 
namespace CASM {
namespace xtal {
class Lattice;
class LatticeMap;
class SimpleStructure;
class StrucMapCalculatorInterface;
 
// In this file:
struct LatticeNode;
struct AssignmentNode;
struct MappingNode;
class StrucMapper;
 
namespace StrucMapping {
 
typedef std::vector<std::vector<Index>> PermuteOpVector;
inline double big_inf() { return 10E20; }
 
inline double small_inf() { return 10E10; }
 
inline bool is_inf(double _val) { return _val > small_inf() / 2.; }
 
double atomic_cost_child(const MappingNode &mapped_result, Index Nsites);
 
double atomic_cost_parent(const MappingNode &mapped_result, Index Nsites);
 
double atomic_cost(const MappingNode &mapped_config, Index Nsites);
 
double atomic_cost(
    const MappingNode &basic_mapping_node, SymOpVector &factor_group,
    const std::vector<Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic,
                                               Index>> &permutation_group,
    Index Nsites);
 
}  // namespace StrucMapping
 
struct LatticeNode {
  Eigen::Matrix3d stretch;
 
  Eigen::Matrix3d isometry;
 
  Superlattice parent;
 
  Superlattice child;
 
  double cost;
 
  LatticeNode(Lattice const &parent_prim, Lattice const &parent_scel,
              Lattice const &child_prim, Lattice const &child_scel,
              Index child_N_atom, double _cost = StrucMapping::big_inf());
 
  LatticeNode(LatticeMap const &_lat_map, Lattice const &parent_prim,
              Lattice const &child_prim);
};
 
bool less(LatticeNode const &A, LatticeNode const &B, double cost_tol);
 
bool identical(LatticeNode const &A, LatticeNode const &B, double cost_tol);
 
struct AssignmentNode {
  AssignmentNode(double _cost_tol = 1e-6)
      : time_reversal(false), cost(0), m_cost_tol(_cost_tol) {}
 
  Eigen::Vector3d translation;
 
  bool time_reversal;
 
  std::set<std::pair<Index, Index>> forced_on;
 
  std::vector<Index> irow;
 
  std::vector<Index> icol;
 
  std::vector<Index> assignment;
 
  Eigen::MatrixXd cost_mat;
 
  double cost;
 
  double cost_tol() const { return m_cost_tol; }
 
  bool empty() const { return cost_mat.size() == 0 && assignment.empty(); }
 
  bool operator<(AssignmentNode const &other) const;
 
  std::vector<Index> permutation() const {
    std::vector<Index> result(assignment.size() + forced_on.size(), 0);
    for (auto const &pair : forced_on) {
      result[pair.first] = pair.second;
    }
    for (Index i = 0; i < assignment.size(); ++i) {
      result[irow[i]] = icol[assignment[i]];
    }
    return result;
  }
 
 private:
  double m_cost_tol;
};
 
bool identical(AssignmentNode const &A, AssignmentNode const &B);
 
struct MappingNode {
  // typedefs to provide flexibility if we eventually change to a
  // Eigen::Matrix<3,Eigen::Dynamic>
  typedef Eigen::MatrixXd DisplacementMatrix;
 
  // Can treat as a Eigen::VectorXd
  using Displacement = DisplacementMatrix::ColXpr;
  using ConstDisplacement = DisplacementMatrix::ConstColXpr;
 
  // Label molecules as name and occupant index (optional, default 0)
  using MoleculeLabel = std::pair<std::string, Index>;
 
  using AtomIndexSet = std::set<Index>;
 
  using MoleculeMap = std::vector<AtomIndexSet>;
 
  LatticeNode lattice_node;
  AssignmentNode atomic_node;
  double lattice_weight;
  double atomic_weight;
 
  bool is_viable;
 
  mutable bool is_valid;
 
  mutable bool is_partitioned;
 
  double cost;
 
  Eigen::MatrixXd atom_displacement;
 
  std::vector<Index> atom_permutation;
 
  MoleculeMap mol_map;
 
  std::vector<MoleculeLabel> mol_labels;
 
  static MappingNode invalid();
 
  MappingNode(LatticeNode _lattice_node, double _lattice_weight)
      : lattice_node(std::move(_lattice_node)),
        is_viable(true),
        is_valid(false),
        is_partitioned(false) {
    set_lattice_weight(_lattice_weight);
    cost = lattice_weight * lattice_node.cost;
  }
 
  double cost_tol() const { return atomic_node.cost_tol(); }
 
  void set_lattice_weight(double _lw) {
    lattice_weight = max(min(_lw, 1.0), 1e-9);
    atomic_weight = 1. - lattice_weight;
  }
 
  std::pair<Index, Index> vol_pair() const {
    return std::pair<Index, Index>(lattice_node.parent.size(),
                                   lattice_node.child.size());
  }
 
  void calc();
 
  void clear() {
    throw std::runtime_error("MappingNode::clear() not implemented");
  }
 
  Eigen::Matrix3d const &isometry() const { return lattice_node.isometry; }
 
  Eigen::Matrix3d const &stretch() const { return lattice_node.stretch; }
 
  Eigen::Vector3d const &translation() const { return atomic_node.translation; }
 
  bool time_reversal() const { return atomic_node.time_reversal; }
 
  Displacement disp(Index i) { return atom_displacement.col(i); }
 
  ConstDisplacement disp(Index i) const { return atom_displacement.col(i); }
 
  bool operator<(MappingNode const &other) const;
};
 
inline Eigen::Matrix3d const &get_matrix(MappingNode const &_node) {
  return _node.isometry();
}
 
inline Eigen::Vector3d const &get_translation(MappingNode const &_node) {
  return _node.translation();
}
 
inline bool get_time_reversal(MappingNode const &_node) {
  return _node.time_reversal();
}
 
class StrucMapper {
 public:
  using LatMapType = std::map<Index, std::vector<Lattice>>;
 
  enum Options {
    none = 0,
    strict = (1u << 0),
    robust = (1u << 1),
    sym_strain = (1u << 2),
    sym_basis = (1u << 3),
    // soft_va_limit ensures that if no supercell volume satisfies vacancy
    // constraints, the smallest possible volume is used. Default behavior
    // results in no valid mapping
    soft_va_limit = (1u << 4)
  };
 
  StrucMapper(StrucMapCalculatorInterface const &_calculator,
              double _lattice_weight = 0.5, double _max_volume_change = 0.5,
              int _options = 0,  // this should actually be a bitwise-OR of
                                 // StrucMapper::Options
              double _cost_tol = TOL, double _min_va_frac = 0.,
              double _max_va_frac = 1.);
 
  double cost_tol() const { return m_cost_tol; }
 
  double xtal_tol() const { return m_xtal_tol; }
 
  double lattice_weight() const { return m_lattice_weight; }
 
  double atomic_weight() const { return 1. - m_lattice_weight; }
 
  void set_lattice_weight(double _lw) {
    m_lattice_weight = max(min(_lw, 1.0), 1e-9);
  }
 
  Index lattice_transformation_range() const {
    return m_lattice_transformation_range;
  }
 
  void set_lattice_transformation_range(Index _new_range) {
    m_lattice_transformation_range = _new_range;
  }
 
  void set_symmetrize_lattice_cost(bool _sym_lat_cost) {
    m_symmetrize_lattice_cost = _sym_lat_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; }
  bool symmetrize_atomic_cost() const { return m_symmetrize_atomic_cost; }
 
  double min_va_frac() const { return m_min_va_frac; }
 
  void set_min_va_frac(double _min_va) { m_min_va_frac = max(_min_va, 0.); }
 
  double max_va_frac() const { return m_max_va_frac; }
 
  void set_max_va_frac(double _max_va) { m_max_va_frac = min(_max_va, 0.99); }
 
  int options() const { return m_options; }
 
  SimpleStructure const &parent() const;
 
  void add_allowed_lattice(Lattice const &_lat) {
    m_allowed_superlat_map
        [std::abs(round(volume(_lat) / parent().lat_column_mat.determinant()))]
            .push_back(_lat);
  }
 
  void clear_allowed_lattices() const { m_allowed_superlat_map.clear(); }
 
  bool lattices_constrained() const { return m_allowed_superlat_map.size(); }
 
  void set_filter(
      std::function<bool(Lattice const &, Lattice const &)> _filter_f) {
    m_filtered = true;
    m_filter_f = _filter_f;
    m_superlat_map.clear();
  }
 
  void unset_filter() {
    m_filtered = false;
    m_superlat_map.clear();
  }
 
  std::set<MappingNode> map_ideal_struc(
      const SimpleStructure &child_struc, Index k = 1,
      double max_cost = StrucMapping::big_inf(), double min_cost = -TOL,
      bool keep_invalid = false) const;
 
  std::set<MappingNode> map_deformed_struc(
      const SimpleStructure &child_struc, Index k = 1,
      double max_cost = StrucMapping::big_inf(), double min_cost = -TOL,
      bool keep_invalid = false,
      SymOpVector const &child_factor_group = {SymOp::identity()}) const;
 
  std::set<MappingNode> map_deformed_struc_impose_lattice_vols(
      const SimpleStructure &child_struc, Index min_vol, Index max_vol,
      Index k = 1, double max_cost = StrucMapping::big_inf(),
      double min_cost = -TOL, bool keep_invalid = false,
      SymOpVector const &child_factor_group = {SymOp::identity()}) const;
 
  std::set<MappingNode> map_deformed_struc_impose_lattice(
      const SimpleStructure &child_struc, const Lattice &imposed_lat,
      Index k = 1, double max_cost = StrucMapping::big_inf(),
      double min_cost = -TOL, bool keep_invalid = false,
      SymOpVector const &child_factor_group = {SymOp::identity()}) const;
 
 
  std::set<MappingNode> map_deformed_struc_impose_lattice_node(
      const SimpleStructure &child_struc, const LatticeNode &imposed_node,
      Index k = 1, double max_cost = StrucMapping::big_inf(),
      double min_cost = -TOL, bool keep_invalid = false) const;
 
 
  Index k_best_maps_better_than(SimpleStructure const &child_struc,
                                std::set<MappingNode> &queue, Index k = 1,
                                double max_cost = StrucMapping::big_inf(),
                                double min_cost = -TOL,
                                bool keep_invalid = false,
                                bool keep_tail = false,
                                bool no_partiton = false) const;
 
  StrucMapCalculatorInterface const &calculator() const { return *m_calc_ptr; }
 
 private:
  std::set<MappingNode> _seed_k_best_from_super_lats(
      SimpleStructure const &child_struc,
      std::vector<Lattice> const &_parent_scels,
      std::vector<Lattice> const &_child_scels, Index k, double max_strain_cost,
      double min_strain_cost,
      SymOpVector const &child_factor_group = {SymOp::identity()}) const;
 
  std::set<MappingNode> _seed_from_vol_range(
      SimpleStructure const &child_struc, Index k, Index min_vol, Index max_vol,
      double max_strain_cost, double min_strain_cost,
      SymOpVector const &child_factor_group = {SymOp::identity()}) const;
 
  Index _n_species(SimpleStructure const &sstruc) const;
 
  // Implement filter as std::function<bool(Lattice const&,Lattice const&)>
  bool _filter_lat(Lattice const &_parent_lat,
                   Lattice const &_child_lat) const {
    return m_filter_f(_parent_lat, _child_lat);
  }
 
  std::pair<Index, Index> _vol_range(const SimpleStructure &child_struc) const;
 
  notstd::cloneable_ptr<StrucMapCalculatorInterface> m_calc_ptr;
 
  Eigen::MatrixXd m_strain_gram_mat;
 
  double m_lattice_weight;
  double m_max_volume_change;
  int m_options;
  double m_cost_tol;
  double m_xtal_tol;
  double m_min_va_frac;
  double m_max_va_frac;
 
  Index m_lattice_transformation_range;
 
  bool m_symmetrize_lattice_cost;
  bool m_symmetrize_atomic_cost;
 
  bool m_filtered;
  std::function<bool(Lattice const &, Lattice const &)> m_filter_f;
 
  mutable LatMapType m_superlat_map;
  mutable LatMapType m_allowed_superlat_map;
 
  std::vector<Lattice> _lattices_of_vol(Index prim_vol) const;
};
 
}  // namespace xtal
}  // namespace CASM
#endif