NeighborList.cc

Go to the documentation of this file.

#include "casm/clex/NeighborList.hh"
 
#include "casm/container/Counter.hh"
#include "casm/crystallography/LinearIndexConverter.hh"
#include "casm/crystallography/Superlattice.hh"
#include "casm/misc/CASM_Eigen_math.hh"
#include "casm/misc/CASM_math.hh"
#include "casm/misc/algorithm.hh"
 
namespace CASM {
 
PrimNeighborList::Matrix3Type PrimNeighborList::weight_matrix() const {
  return m_W;
}
 
void PrimNeighborList::expand(UnitCellCoord const &uccoord) {
  expand(uccoord.unitcell());
}
 
void PrimNeighborList::expand(UnitCell const &uc) {
  // save the old range
  Scalar prev_range = m_range;
 
  auto result = m_neighborhood.insert(uc);
 
  // if a new UnitCell
  if (result.second) {
    // ensure all intermediate UnitCell are included
    _expand(prev_range);
  }
}
 
bool PrimNeighborList::_expand(UnitCellCoord const &uccoord) {
  return _expand(uccoord.unitcell());
}
 
bool PrimNeighborList::_expand(UnitCell const &uc) {
  return m_neighborhood.insert(uc).second;
}
 
void PrimNeighborList::_expand(Scalar prev_range) {
  // if some are new, we need to make sure we including all intermediate
  // UnitCell in our neighborhood get the score of the last element in the set
  m_range = _dist(*m_neighborhood.rbegin());
 
  // count over possible UnitCell, adding to neighborhood if prev_range < dist
  // <= m_range we're using a std::set for m_neighborhood, so it gets sorted
 
  // We want the bounding box that contains the ellipsoid defined by
  //   m_range = (i,j,k).transpose() * W * (i,j,k)
  //           = F.t * U.t * U * F
  //           = F.t * (M.inv).t * (M.inv) * F
  //         1 = F.t * (sqrt(m_range)*M.inv).t * (sqrt(m_range)*M.inv) * F
  // has maximum bb_end
  //   bb_end(i) = M.row(i).norm() (from considerations of an affine
  //   transformation M of unit sphere) where M = sqrt(m_range)*U.inv
  // and minimum
  //   bb_begin = -bb_end
 
  Eigen::MatrixXd M = sqrt(1.0 * m_range) * m_Uinv;
 
  auto d = [&](int i) { return std::lround(std::ceil(M.row(i).norm())); };
 
  VectorXType bb_end(3);
  bb_end << d(0), d(1), d(2);
  VectorXType bb_begin = -bb_end;
  VectorXType bb_incr = VectorXType::Constant(3, 1);
 
  typedef Counter<VectorXType, Scalar, Index,
                  CASM_TMP::ParenthesesAccess<VectorXType, Scalar, Index> >
      VectorXCounter;
 
  VectorXCounter counter(bb_begin, bb_end, bb_incr);
 
  Scalar dist;
  for (; counter.valid(); ++counter) {
    dist = _dist(counter.current());
    if (prev_range < dist && dist <= m_range) {
      m_neighborhood.insert(counter.current());
    }
  }
}
 
PrimNeighborList::size_type PrimNeighborList::size() const {
  return m_neighborhood.size();
}
 
PrimNeighborList::const_iterator PrimNeighborList::begin() const {
  return m_neighborhood.cbegin();
}
 
PrimNeighborList::const_iterator PrimNeighborList::end() const {
  return m_neighborhood.cend();
}
 
PrimNeighborList::const_iterator PrimNeighborList::cbegin() const {
  return m_neighborhood.cbegin();
}
 
PrimNeighborList::const_iterator PrimNeighborList::cend() const {
  return m_neighborhood.cend();
}
 
const PrimNeighborList::SublatIndices &PrimNeighborList::sublat_indices()
    const {
  return m_sublat_indices;
}
 
PrimNeighborList::Scalar PrimNeighborList::neighbor_index(
    UnitCellCoord const &_ucc) {
  expand(_ucc);
  return _neighbor_index(_ucc);
}
 
PrimNeighborList::Scalar PrimNeighborList::_neighbor_index(
    UnitCellCoord const &_ucc) const {
  Scalar uc_ind(find_index(m_neighborhood, _ucc.unitcell()));
  Scalar sublat_dist(find_index(sublat_indices(), _ucc.sublattice()));
 
  return uc_ind * sublat_indices().size() + sublat_dist;
}
 
PrimNeighborList::Scalar PrimNeighborList::_dist(const UnitCell &A) const {
  return A.transpose() * m_W * A;
}
 
PrimNeighborList::VectorXType PrimNeighborList::_add_dist(
    const UnitCell &A) const {
  VectorXType vec(4);
  vec(0) = _dist(A);
  vec.segment(1, 3) = A;
  return vec;
}
 
bool PrimNeighborList::_compare_unitcell(const UnitCell &A,
                                         const UnitCell &B) const {
  return _compare_vec(_add_dist(A), _add_dist(B));
}
 
bool PrimNeighborList::_compare_vec(const VectorXType &A,
                                    const VectorXType &B) {
  return std::lexicographical_compare(A.data(), A.data() + A.size(), B.data(),
                                      B.data() + B.size());
}
 
PrimNeighborList::Matrix3Type PrimNeighborList::make_weight_matrix(
    const Eigen::Matrix3d lat_column_mat, Index max_element_value, double tol) {
  Eigen::Matrix3d W = lat_column_mat.transpose() * lat_column_mat;
 
  double min = std::numeric_limits<double>::max();
  for (int i = 0; i < 3; ++i) {
    for (int j = 0; j < 3; ++j) {
      if (!almost_zero(W(i, j), tol) && std::abs(W(i, j)) < std::abs(min)) {
        min = W(i, j);
      }
    }
  }
  W /= min;
 
  // make positive definite
  Eigen::LLT<Eigen::MatrixXd> llt(W);
  if (llt.info() != Eigen::Success) {
    W = -W;
  }
 
  double n = 1.0;
  while (!is_integer(n * W, tol) &&
         ((n + 1.0) * W).cwiseAbs().maxCoeff() < max_element_value) {
    n += 1.0;
  }
  return lround(n * W);
}
 
std::unique_ptr<PrimNeighborList> PrimNeighborList::clone() const {
  return std::unique_ptr<PrimNeighborList>(new PrimNeighborList(*this));
}
 
SuperNeighborList::SuperNeighborList(
    Eigen::Matrix3l const &transformation_matrix_to_super,
    const PrimNeighborList &prim_nlist) {
  xtal::UnitCellIndexConverter ijk_index_converter{
      transformation_matrix_to_super};
  m_prim_grid_size = ijk_index_converter.total_sites();
  m_site.resize(m_prim_grid_size);
  m_unitcell.resize(m_prim_grid_size);
 
  // use the PrimNeighborList to generate the UnitCell and Site indices for
  //   the neighbors of each UnitCell in the supercell
 
  // for each unit cell in the supercell
  for (Index i = 0; i < m_prim_grid_size; ++i) {
    // for each neighbor unitcell
    for (auto it = prim_nlist.begin(); it != prim_nlist.end(); ++it) {
      // get the neighbor unitcell's index
      UnitCell neighbor_unitcell = ijk_index_converter(i) + *it;
      size_type neighbor_unitcell_index =
          ijk_index_converter(neighbor_unitcell);
 
      // store the unitcell index
      m_unitcell[i].push_back(neighbor_unitcell_index);
 
      // calculate and store the site indices for all sites in the neighbor
      // unitcell that are requested
      // - Depends on Configuration sites being stored in blocks by sublattice
      // and unitcell indices
      //   determined by the UnitCellCoordIndexConverter ordering
      for (auto b_it = prim_nlist.sublat_indices().begin();
           b_it != prim_nlist.sublat_indices().end(); ++b_it) {
        m_site[i].push_back((*b_it) * m_prim_grid_size +
                            neighbor_unitcell_index);
      }
    }
  }
}
 
SuperNeighborList::SuperNeighborList(const xtal::Superlattice &superlattice,
                                     const PrimNeighborList &prim_nlist)
    : SuperNeighborList(superlattice.transformation_matrix_to_super(),
                        prim_nlist) {}
 
const std::vector<SuperNeighborList::size_type> &SuperNeighborList::sites(
    size_type unitcell_index) const {
  return m_site[unitcell_index];
}
 
const std::vector<SuperNeighborList::size_type> &SuperNeighborList::unitcells(
    size_type unitcell_index) const {
  return m_unitcell[unitcell_index];
}
 
bool SuperNeighborList::overlaps() const {
  // there is an overlap if any of the neighboring unitcell indices is repeated
  // so sort and check if any two neighboring indices are the same
 
  std::vector<size_type> nlist = m_unitcell[0];
  std::sort(nlist.begin(), nlist.end());
  return std::adjacent_find(nlist.begin(), nlist.end()) != nlist.end();
}
 
std::unique_ptr<SuperNeighborList> SuperNeighborList::clone() const {
  return std::unique_ptr<SuperNeighborList>(new SuperNeighborList(*this));
}
 
}  // namespace CASM