NeighborList.hh

Go to the documentation of this file.

#ifndef CASM_NeighborList_HH
#define CASM_NeighborList_HH

#include "casm/crystallography/UnitCellCoord.hh"

namespace CASM {

  class PrimGrid;

  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(UnitCellCoord 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);

    std::unique_ptr<PrimNeighborList> clone() const;

  private:

    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(const PrimGrid &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) {
      if(m_neighborhood.insert(it->unitcell()).second) {
        any_new = true;
      }
    }

    if(!any_new) {
      return;
    }

    // otherwise, ensure all intermediate UnitCell are included
    _expand(prev_range);
  }

}

#endif