NeighborList.hh

Go to the documentation of this file.

#ifndef CASM_NeighborList_HH
#define CASM_NeighborList_HH
 
#include <memory>
#include <set>
#include <vector>
 
#include "casm/crystallography/UnitCellCoord.hh"
 
namespace CASM {
 
namespace xtal {
class UnitCell;
class UnitCellCoord;
class Superlattice;
}  // namespace xtal
using xtal::UnitCell;
using xtal::UnitCellCoord;
 
class PrimNeighborList {
 public:
  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
 
  typedef Index Scalar;
  typedef Eigen::Matrix3l Matrix3Type;
  typedef Eigen::Matrix<long, Eigen::Dynamic, 1> VectorXType;
  typedef std::set<UnitCell, std::function<bool(UnitCell, UnitCell)> >
      NeighborSet;
  typedef NeighborSet::size_type size_type;
  typedef NeighborSet::const_iterator const_iterator;
  typedef std::set<int> SublatIndices;
 
  template <typename SublatIterator>
  PrimNeighborList(const Matrix3Type &W, SublatIterator begin,
                   SublatIterator end);
 
  Matrix3Type weight_matrix() const;
 
  void expand(UnitCell const &uc);
 
  void expand(UnitCellCoord const &uccoord);
 
  template <typename UnitCellCoordIterator>
  void expand(UnitCellCoordIterator begin, UnitCellCoordIterator end);
 
  size_type size() const;
 
  const_iterator begin() const;
 
  const_iterator end() const;
 
  const_iterator cbegin() const;
 
  const_iterator cend() const;
 
  const SublatIndices &sublat_indices() const;
 
  static Matrix3Type make_weight_matrix(const Eigen::Matrix3d lat_column_mat,
                                        Index max_element_value, double tol);
 
  Scalar neighbor_index(UnitCellCoord const &_ucc);
 
  template <typename UnitCellCoordIterator>
  std::vector<Scalar> neighbor_indices(UnitCellCoordIterator _begin,
                                       UnitCellCoordIterator _end);
 
  std::unique_ptr<PrimNeighborList> clone() const;
 
 private:
  Scalar _neighbor_index(UnitCellCoord const &_ucc) const;
 
  bool _expand(UnitCell const &uc);
 
  bool _expand(UnitCellCoord const &uccoord);
 
  void _expand(Scalar prev_range);
 
  Scalar _dist(const UnitCell &A) const;
 
  VectorXType _add_dist(const UnitCell &A) const;
 
  bool _compare_unitcell(const UnitCell &A, const UnitCell &B) const;
 
  static bool _compare_vec(const VectorXType &A, const VectorXType &B);
 
  Matrix3Type m_W;
 
  Eigen::MatrixXd m_Uinv;
 
  NeighborSet m_neighborhood;
 
  Scalar m_range;
 
  SublatIndices m_sublat_indices;
};
 
class SuperNeighborList {
 public:
  typedef Index size_type;
 
  SuperNeighborList(Eigen::Matrix3l const &transformation_matrix_to_super,
                    PrimNeighborList const &prim_nlist);
 
  SuperNeighborList(const xtal::Superlattice &prim_grid,
                    const PrimNeighborList &prim_nlist);
 
  size_type unitcell_index(size_type site_index) const {
    return site_index % m_prim_grid_size;
  }
 
  size_type sublat_index(size_type site_index) const {
    return site_index / m_prim_grid_size;
  }
 
  const std::vector<size_type> &sites(size_type unitcell_index) const;
 
  const std::vector<size_type> &unitcells(size_type unitcell_index) const;
 
  bool overlaps() const;
 
  std::unique_ptr<SuperNeighborList> clone() const;
 
 private:
  size_type m_prim_grid_size;
 
  std::vector<std::vector<size_type> > m_site;
 
  std::vector<std::vector<size_type> > m_unitcell;
};
 
template <typename SublatIterator>
PrimNeighborList::PrimNeighborList(const Matrix3Type &W, SublatIterator begin,
                                   SublatIterator end)
    : m_W(W),
      m_neighborhood(std::bind(&PrimNeighborList::_compare_unitcell, this,
                               std::placeholders::_1, std::placeholders::_2)),
      m_range(0),
      m_sublat_indices(begin, end) {
  // throw if W is not positive definite
  Eigen::LLT<Eigen::MatrixXd> llt(W.cast<double>());
  if (llt.info() != Eigen::Success) {
    std::cerr << "Error constructing PrimNeighborList: weight matrix is not "
                 "positive definite."
              << std::endl;
    std::cerr << "weight matrix: \n" << W << std::endl;
    throw std::runtime_error(
        "Error constructing PrimNeighborList: weight matrix is not positive "
        "definite.");
  }
 
  // store U.inverse()
  Eigen::MatrixXd U = llt.matrixU();
  m_Uinv = U.inverse();
 
  // add origin unit cell to neighborhood
  m_neighborhood.insert(UnitCell(0, 0, 0));
}
 
template <typename UnitCellCoordIterator>
void PrimNeighborList::expand(UnitCellCoordIterator begin,
                              UnitCellCoordIterator end) {
  // save the old range
  Scalar prev_range = m_range;
 
  bool any_new = false;
  for (auto it = begin; it != end; ++it) {
    any_new = _expand(*it) || any_new;
  }
 
  if (!any_new) {
    return;
  }
 
  // otherwise, ensure all intermediate UnitCell are included
  _expand(prev_range);
}
 
template <typename UnitCellCoordIterator>
std::vector<PrimNeighborList::Scalar> PrimNeighborList::neighbor_indices(
    UnitCellCoordIterator _begin, UnitCellCoordIterator _end) {
  expand(_begin, _end);
  std::vector<Scalar> result;
  std::transform(_begin, _end, std::back_inserter(result),
                 [this](UnitCellCoord const &A) -> Scalar {
                   return this->_neighbor_index(A);
                 });
  return result;
}
 
}  // namespace CASM
 
#endif