neighbours.hpp

Go to the documentation of this file.

//-----------------------------------------------------------------------------------
// d-SEAMS - Deferred Structural Elucidation Analysis for Molecular Simulations
//
// Copyright (c) 2018--present d-SEAMS core team
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the MIT License as published by
// the Open Source Initiative.
//
// A copy of the MIT License is included in the LICENSE file of this repository.
// You should have received a copy of the MIT License along with this program.
// If not, see <https://opensource.org/licenses/MIT>.
//-----------------------------------------------------------------------------------
 
#ifndef SEAMS_NEIGHBOURS_H_
#define SEAMS_NEIGHBOURS_H_
 
#include <algorithm>
#include <array>
#include <string>
#include <tuple>
#include <utility>
#include <vector>
 
#include <generic.hpp>
#include <mol_sys.hpp>
 
#ifdef SEAMS_HAS_LINKCELL
#include <linkcell.hpp>
#endif
#ifdef SEAMS_HAS_MINIMAGE
#include <minimage.h>
#endif



 
namespace nneigh {

inline void dumpBoundsToH(const std::vector<double> &box,
                          const std::vector<double> &boxLow, double H[3][3],
                          double origin[3]) {
  const double xspan = box.size() > 0 ? box[0] : 0.0;
  const double yspan = box.size() > 1 ? box[1] : 0.0;
  const double zspan = box.size() > 2 ? box[2] : 0.0;
  const double xlo_b = boxLow.size() > 0 ? boxLow[0] : 0.0;
  const double ylo_b = boxLow.size() > 1 ? boxLow[1] : 0.0;
  const double zlo_b = boxLow.size() > 2 ? boxLow[2] : 0.0;
  const double xy = box.size() >= 6 ? box[3] : 0.0;
  const double xz = box.size() >= 6 ? box[4] : 0.0;
  const double yz = box.size() >= 6 ? box[5] : 0.0;
  const double xmin = std::min(std::min(0.0, xy), std::min(xz, xy + xz));
  const double xmax = std::max(std::max(0.0, xy), std::max(xz, xy + xz));
  const double ymin = std::min(0.0, yz);
  const double ymax = std::max(0.0, yz);
  H[0][0] = xspan - xmax + xmin;
  H[0][1] = 0.0;
  H[0][2] = 0.0;
  H[1][0] = xy;
  H[1][1] = yspan - ymax + ymin;
  H[1][2] = 0.0;
  H[2][0] = xz;
  H[2][1] = yz;
  H[2][2] = zspan;
  origin[0] = xlo_b - xmin;
  origin[1] = ylo_b - ymin;
  origin[2] = zlo_b;
}

inline void dumpCellLengths(const std::vector<double> &box,
                            const std::vector<double> &boxLow,
                            double lengths[3]) {
  double H[3][3];
  double origin[3];
  dumpBoundsToH(box, boxLow, H, origin);
  lengths[0] = H[0][0];
  lengths[1] = H[1][1];
  lengths[2] = H[2][2];
}

inline int dumpAxialDim(const std::vector<double> &box,
                        const std::vector<double> &boxLow) {
  double lengths[3];
  dumpCellLengths(box, boxLow, lengths);
  int axial = 0;
  if (lengths[1] > lengths[axial]) {
    axial = 1;
  }
  if (lengths[2] > lengths[axial]) {
    axial = 2;
  }
  return axial;
}
 
inline int dumpAxialDim(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud) {
  return dumpAxialDim(yCloud.box, yCloud.boxLow);
}

inline void dumpToFrac(const double H[3][3], const double origin[3], double x,
                       double y, double z, double s[3]) {
  const double lx = H[0][0];
  const double ly = H[1][1];
  const double lz = H[2][2];
  const double xy = H[1][0];
  const double xz = H[2][0];
  const double yz = H[2][1];
  s[2] = (z - origin[2]) / lz;
  s[1] = (y - origin[1] - yz * s[2]) / ly;
  s[0] = (x - origin[0] - xy * s[1] - xz * s[2]) / lx;
}

inline void dumpFromFrac(const double H[3][3], const double origin[3],
                         const double s[3], double r[3]) {
  r[0] = origin[0] + H[0][0] * s[0] + H[1][0] * s[1] + H[2][0] * s[2];
  r[1] = origin[1] + H[1][1] * s[1] + H[2][1] * s[2];
  r[2] = origin[2] + H[2][2] * s[2];
}

inline double dumpVolume(const std::vector<double> &box,
                         const std::vector<double> &boxLow) {
  double H[3][3];
  double origin[3];
  dumpBoundsToH(box, boxLow, H, origin);
  const double vol = H[0][0] * H[1][1] * H[2][2];
  return vol < 0.0 ? -vol : vol;
}
 
inline double dumpVolume(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud) {
  return dumpVolume(yCloud.box, yCloud.boxLow);
}
 
#ifdef SEAMS_HAS_LINKCELL
inline lc_cell lammpsBoxToLcCell(const std::vector<double> &box,
                                 const std::vector<double> &boxLow) {
#ifdef SEAMS_HAS_MINIMAGE
  const double xspan = box.size() > 0 ? box[0] : 0.0;
  const double yspan = box.size() > 1 ? box[1] : 0.0;
  const double zspan = box.size() > 2 ? box[2] : 0.0;
  const double xy = box.size() >= 6 ? box[3] : 0.0;
  const double xz = box.size() >= 6 ? box[4] : 0.0;
  const double yz = box.size() >= 6 ? box[5] : 0.0;
  const double xlo_b = boxLow.size() > 0 ? boxLow[0] : 0.0;
  const double ylo_b = boxLow.size() > 1 ? boxLow[1] : 0.0;
  const double zlo_b = boxLow.size() > 2 ? boxLow[2] : 0.0;
  mi_cell raw;
  mi_cell_from_lammps_bounds(xspan, yspan, zspan, xy, xz, yz, xlo_b, ylo_b,
                             zlo_b, &raw);
  lc_cell c = lc_cell_ortho(raw.ax, raw.by, raw.cz);
  c.ax = raw.ax;
  c.ay = raw.ay;
  c.az = raw.az;
  c.bx = raw.bx;
  c.by = raw.by;
  c.bz = raw.bz;
  c.cx = raw.cx;
  c.cy = raw.cy;
  c.cz = raw.cz;
  c.ox = raw.ox;
  c.oy = raw.oy;
  c.oz = raw.oz;
  return c;
#else
  double H[3][3], o[3];
  dumpBoundsToH(box, boxLow, H, o);
  lc_cell c = lc_cell_ortho(H[0][0], H[1][1], H[2][2]);
  c.ax = H[0][0];
  c.ay = H[0][1];
  c.az = H[0][2];
  c.bx = H[1][0];
  c.by = H[1][1];
  c.bz = H[1][2];
  c.cx = H[2][0];
  c.cy = H[2][1];
  c.cz = H[2][2];
  c.ox = o[0];
  c.oy = o[1];
  c.oz = o[2];
  return c;
#endif
}
 
inline linkcell::Cell lammpsBoxToLinkcell(const std::vector<double> &box,
                                          const std::vector<double> &boxLow) {
  return linkcell::Cell(lammpsBoxToLcCell(box, boxLow));
}

inline void residentFrameCell(const double *box, const double *boxLow,
                              int nBox, linkcell::Cell &cell,
                              const double *&frameLens) {
  cell = linkcell::Cell::ortho(box[0], box[1], box[2]);
  frameLens = box;
  if (boxLow != nullptr || nBox >= 6) {
    std::vector<double> dump(static_cast<std::size_t>(std::max(nBox, 3)));
    for (int i = 0; i < nBox && i < static_cast<int>(dump.size()); ++i) {
      dump[static_cast<std::size_t>(i)] = box[i];
    }
    std::vector<double> lo;
    if (boxLow != nullptr) {
      lo = {boxLow[0], boxLow[1], boxLow[2]};
    }
    cell = lammpsBoxToLinkcell(dump, lo);
    frameLens = nullptr;
  }
}
#endif


std::vector<std::vector<int>> neighList(
    double rcutoff, const molSys::PointCloud<molSys::Point<double>, double> &yCloud,
    int typeI, int typeJ);

std::vector<std::vector<int>> neighListPair(
    double rcutoff, const molSys::PointCloud<molSys::Point<double>, double> &yCloud,
    int typeI, int typeJ);

bool cellListRowsThreaded(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud,
    const std::vector<int> &subset, double rcutoff,
    std::vector<std::vector<int>> &rows);

std::vector<std::vector<int>> neighListO(
    double rcutoff, const molSys::PointCloud<molSys::Point<double>, double> &yCloud,
    int typeI);

std::vector<std::vector<int>> halfNeighList(
    double rcutoff, const molSys::PointCloud<molSys::Point<double>, double> &yCloud,
    int typeI = 1);


std::vector<std::vector<int>> neighbourListByIndex(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud,
    const std::vector<std::vector<int>> &nList);

std::vector<std::vector<int>> getNewNeighbourListByIndex(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud, double cutoff);

std::vector<std::vector<int>> kNearestNeighbourList(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud, int k,
    double candidateCutoff, int typeI, bool mutual = true);

std::vector<std::vector<int>> kNearestNeighbourList(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud, int k,
    double candidateCutoff, const std::vector<int> &types, bool mutual = true);

std::pair<std::vector<std::vector<int>>, std::vector<std::vector<int>>>
kNearestNeighbourPair(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud, int k,
    double candidateCutoff, int typeI);
 
std::pair<std::vector<std::vector<int>>, std::vector<std::vector<int>>>
kNearestNeighbourPair(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud, int k,
    double candidateCutoff, const std::vector<int> &types);

std::pair<double, double> shellSeparation(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud, int k,
    int typeI);

std::vector<std::tuple<int, int, double>> nearestUnlike(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud, int typeI,
    int typeJ);

std::vector<std::pair<int, int>> mutualNearestUnlike(
    const molSys::PointCloud<molSys::Point<double>, double> &yCloud, int typeI,
    int typeJ);

[[nodiscard]] int clearNeighbourList(std::vector<std::vector<int>> &nList);

enum class BondGraph { Cutoff, KnnMutual, KnnUnion };
 
BondGraph bondGraphFromName(const std::string &name);
const char *bondGraphName(BondGraph graph);

class SkinNeighborList {
public:
  SkinNeighborList(double cutoff, double skin, int typeI,
                   BondGraph graph = BondGraph::KnnMutual, int k = 4);
 
  [[nodiscard]] BondGraph graph() const { return graph_; }

  const std::vector<std::vector<int>> &
  update(const molSys::PointCloud<molSys::Point<double>, double> &yCloud);
 
  [[nodiscard]] bool lastRebuilt() const { return rebuilt_; }
  [[nodiscard]] int lastChangedAtoms() const { return changedAtoms_; }
  [[nodiscard]] const std::vector<std::vector<int>> &bonds() const {
    return nList_;
  }
 
private:
  double cutoff_;
  double skin_;
  double cutoffSq_;
  double triggerSq_;
  int typeI_;
  int k_;
  BondGraph graph_{BondGraph::KnnMutual};
  bool mutual_{true};
  bool rebuilt_{true};
  int changedAtoms_{0};
  std::vector<double> x0_;
  std::vector<double> y0_;
  std::vector<double> z0_;
  std::vector<double> box0_;
  std::vector<std::pair<int, int>> candidates_;
  std::vector<std::pair<int, int>> bonded_;
  std::vector<std::vector<int>> nList_;
 
  bool mustRebuild(
      const molSys::PointCloud<molSys::Point<double>, double> &yCloud) const;
  void rebuildCandidates(
      const molSys::PointCloud<molSys::Point<double>, double> &yCloud);
  void refreshBonds(
      const molSys::PointCloud<molSys::Point<double>, double> &yCloud);
};
 
}  // namespace nneigh
 
#endif  // SEAMS_NEIGHBOURS_H_