Program Listing for File cell_lists.hpp#
↰ Return to documentation for file (source/pele/cell_lists.hpp)
#ifndef _PELE_CELL_LISTS_H_
#define _PELE_CELL_LISTS_H_
#include <algorithm>
#include <array>
#include <cassert>
#include <cstddef>
#include <exception>
#include <iostream>
#include <limits>
#include <memory>
#include <omp.h>
#include <ostream>
#include <sstream>
#include <vector>
#include "array.hpp"
#include "base_potential.hpp"
#include "distance.hpp"
#include "queue.hpp"
#include "vecn.hpp"
#ifndef NDEBUG
#define ASSERT_EX(condition, statement) \
do { \
if (!(condition)) { \
statement; \
assert(condition); \
} \
} while (false)
#else
#define ASSERT_EX(condition, statement) ((void)0)
#endif
namespace pele {
// Type of each cell inside the cell list
using cell_t = std::vector<size_t>;
template <size_t ndim>
class CellListsContainer {
protected:
void setup_cells(const size_t nsubdoms,
const std::vector<size_t> subdom_ncells) {
#ifdef _OPENMP
#pragma omp parallel
{
size_t isubdom = omp_get_thread_num();
m_cell_atoms[isubdom] = std::vector<cell_t>(subdom_ncells[isubdom]);
}
#else
for (size_t isubdom = 0; isubdom < nsubdoms; ++isubdom) {
m_cell_atoms[isubdom] = std::vector<cell_t>(subdom_ncells[isubdom]);
}
#endif
size_t total_cells = 0;
for (const size_t subdom_ncell : subdom_ncells) {
total_cells += subdom_ncell;
}
m_cell_neighbors = std::vector<std::vector<cell_t *>>(total_cells);
for (std::vector<cell_t *> &neighbors : m_cell_neighbors) {
neighbors = std::vector<cell_t *>();
}
}
public:
std::vector<std::vector<cell_t>> m_cell_atoms;
std::vector<std::vector<std::array<cell_t *, 2>>> m_cell_neighbor_pairs_inner;
std::vector<std::vector<std::array<cell_t *, 2>>>
m_cell_neighbor_pairs_boundary;
std::vector<std::vector<cell_t *>> m_cell_neighbors;
CellListsContainer(const std::vector<size_t> subdom_ncells)
: m_cell_atoms(subdom_ncells.size()),
m_cell_neighbor_pairs_inner(subdom_ncells.size()),
m_cell_neighbor_pairs_boundary(subdom_ncells.size()) {
setup_cells(subdom_ncells.size(), subdom_ncells);
}
void reset() {
#ifdef _OPENMP
#pragma omp parallel
{
size_t isubdom = omp_get_thread_num();
for (auto &v : m_cell_atoms[isubdom]) {
v.clear();
}
}
#else
size_t nsubdoms = m_cell_atoms.size();
for (size_t isubdom = 0; isubdom < nsubdoms; ++isubdom) {
for (auto &v : m_cell_atoms[isubdom]) {
v.clear();
}
}
#endif
}
inline void add_atom_to_cell(const size_t iatom, const size_t icell,
const size_t isubdom) {
m_cell_atoms[isubdom][icell].push_back(iatom);
}
inline void remove_atom_from_cell_ind(const size_t index, const size_t icell,
const size_t isubdom) {
std::swap(m_cell_atoms[isubdom][icell][index],
m_cell_atoms[isubdom][icell].back());
m_cell_atoms[isubdom][icell].pop_back();
}
inline void remove_atom_from_cell(const size_t iatom, const size_t icell,
const size_t isubdom) {
auto atom_ptr = std::find(m_cell_atoms[isubdom][icell].begin(),
m_cell_atoms[isubdom][icell].end(), iatom);
std::swap(*atom_ptr, m_cell_atoms[isubdom][icell].back());
m_cell_atoms[isubdom][icell].pop_back();
}
};
inline double get_ncellx_scale(pele::Array<double> radii,
pele::Array<double> boxv, size_t omp_threads) {
double ndim = boxv.size();
double ncellsx_max =
std::max<size_t>(omp_threads, pow(radii.size(), 1 / ndim));
double rcut = radii.get_max() * 2;
size_t ncellsx = boxv[0] / rcut;
int ncellsx_scale;
if (ncellsx <= ncellsx_max) {
if (ncellsx >= omp_threads) {
ncellsx_scale = 1;
} else {
ncellsx_scale = ceil(omp_threads / ncellsx);
}
} else {
ncellsx_scale = ncellsx_max / ncellsx;
}
return ncellsx_scale;
}
} // namespace pele
namespace pele {
template <typename distance_policy>
class LatticeNeighbors {
protected:
void calc_subdom_stats() {
size_t remainder = m_ncells_vec[1] % m_nsubdoms;
m_subdoms_balanced = remainder == 0;
size_t subdom_len;
m_subdom_avg_len = 0;
m_subdom_ncells = std::vector<size_t>(m_nsubdoms);
m_subdom_limits = std::vector<size_t>(m_nsubdoms + 1);
m_subdom_limits[0] = 0;
std::stringstream ncells_stream;
ncells_stream << "Length per subdomain: ";
for (size_t isubdom = 0; isubdom < m_nsubdoms; ++isubdom) {
if (remainder > isubdom) {
subdom_len = m_ncells_vec[1] / m_nsubdoms + 1;
} else {
subdom_len = m_ncells_vec[1] / m_nsubdoms;
}
m_subdom_avg_len += subdom_len;
m_subdom_limits[isubdom + 1] = m_subdom_limits[isubdom] + subdom_len;
m_subdom_ncells[isubdom] = subdom_len * m_ncells / m_ncells_vec[1];
ncells_stream << subdom_len;
if (isubdom < m_nsubdoms - 1) {
ncells_stream << ", ";
}
}
m_subdom_avg_len /= m_nsubdoms;
std::cout << ncells_stream.str() << std::endl;
#ifdef _OPENMP
if ((m_boxvec[1] / m_rcut) / m_nsubdoms < 2 && omp_get_max_threads() > 1) {
throw std::runtime_error(
"Some subdomains are thinner than 2*rcut. The "
"parallelization would not be stable. Reduce the number of threads "
"(environment variable OMP_NUM_THREADS)!");
}
#else
if (m_ncells_vec[1] / m_nsubdoms < 1) {
throw std::runtime_error("Some subdomains are thinner than 1 row!");
}
#endif
}
public:
static const size_t ndim = distance_policy::_ndim;
const std::shared_ptr<distance_policy> m_dist;
typedef VecN<ndim, size_t> cell_vec_t;
pele::VecN<ndim> m_boxvec;
pele::VecN<ndim> m_inv_boxvec;
double m_rcut;
cell_vec_t m_ncells_vec;
pele::VecN<ndim> m_rcell_vec;
size_t m_ncells;
size_t m_nsubdoms;
std::vector<size_t> m_subdom_limits;
std::vector<size_t> m_subdom_ncells;
bool m_subdoms_balanced;
size_t m_subdom_avg_len;
LatticeNeighbors(std::shared_ptr<distance_policy> const &dist,
pele::Array<double> const &boxvec, const double rcut,
pele::Array<size_t> const &ncells_vec)
: m_dist(dist),
m_boxvec(boxvec),
m_rcut(rcut),
m_ncells_vec(ncells_vec),
m_inv_boxvec(ndim) {
#ifdef _OPENMP
m_nsubdoms = omp_get_max_threads();
#else
m_nsubdoms = 1;
#endif
m_ncells = m_ncells_vec.prod();
for (size_t idim = 0; idim < ndim; ++idim) {
m_inv_boxvec[idim] = 1 / m_boxvec[idim];
m_rcell_vec[idim] = m_boxvec[idim] / m_ncells_vec[idim];
}
calc_subdom_stats();
}
size_t cell_vec_to_global_ind(cell_vec_t const &v) const {
// note: the speed of this function is important because it is called
// once per atom each time get_energy is called.
size_t cum = 1;
size_t index = 0;
for (size_t idim = 0; idim < ndim; ++idim) {
index += v[idim] * cum;
cum *= m_ncells_vec[idim];
}
return index;
}
size_t cell_vec_to_local_ind(cell_vec_t const &v,
const size_t isubdom) const {
// note: the speed of this function is important because it is called
// once per atom each time get_energy is called.
size_t cum = 1;
size_t index = 0;
for (size_t idim = 0; idim < ndim; ++idim) {
if (idim == 1) {
index += (v[1] - m_subdom_limits[isubdom]) * cum;
cum *= m_subdom_limits[isubdom + 1] - m_subdom_limits[isubdom];
} else {
index += v[idim] * cum;
cum *= m_ncells_vec[idim];
}
}
return index;
}
cell_vec_t global_ind_to_cell_vec(const size_t icell) const {
cell_vec_t v;
size_t remaining_icell = icell;
size_t fraction;
for (size_t idim = 0; idim < ndim - 1; ++idim) {
fraction = (size_t)remaining_icell / m_ncells_vec[idim];
// should be optimized by compiler to use remainder of previous division
v[idim] = remaining_icell % m_ncells_vec[idim];
remaining_icell = fraction;
}
v[ndim - 1] = remaining_icell;
return v;
}
cell_vec_t local_ind_to_cell_vec(const size_t icell,
const size_t isubdom) const {
cell_vec_t v;
size_t remaining_icell = icell;
size_t fraction;
for (size_t idim = 0; idim < ndim - 1; ++idim) {
if (idim == 1) {
size_t split_len =
m_subdom_limits[isubdom + 1] - m_subdom_limits[isubdom];
fraction = (size_t)remaining_icell / split_len;
// should be optimized by compiler to use remainder of previous division
v[idim] = (remaining_icell % split_len) + m_subdom_limits[isubdom];
} else {
fraction = (size_t)remaining_icell / m_ncells_vec[idim];
// should be optimized by compiler to use remainder of previous division
v[idim] = remaining_icell % m_ncells_vec[idim];
}
remaining_icell = fraction;
}
v[ndim - 1] = remaining_icell;
if (ndim - 1 == 1) {
v[ndim - 1] += m_subdom_limits[isubdom];
}
return v;
}
VecN<ndim> cell_vec_to_position(cell_vec_t const &v) const {
VecN<ndim> x;
for (size_t idim = 0; idim < ndim; ++idim) {
x[idim] = m_rcell_vec[idim] * v[idim];
}
return x;
}
cell_vec_t position_to_cell_vec(const double *const x) const {
double x_in_box[ndim];
m_dist->put_atom_in_box(x_in_box, x);
cell_vec_t cell_vec;
for (size_t idim = 0; idim < ndim; ++idim) {
cell_vec[idim] =
m_ncells_vec[idim] * (x_in_box[idim] * m_inv_boxvec[idim] + 0.5 -
std::numeric_limits<double>::epsilon());
ASSERT_EX(
cell_vec[idim] < m_ncells_vec[idim],
std::cout << "This error occurs when the coordinates are so "
"unreasonably large that the floating point precision "
"is not enough to figure out the cell index."
<< std::endl;
std::cout << "x_in_box[idim] = " << x_in_box[idim] << std::endl;
std::cout << "m_inv_boxvec[idim] = " << m_inv_boxvec[idim]
<< std::endl;
std::cout << "cell_vec[idim] = " << cell_vec[idim] << std::endl;
std::cout << "m_ncells_vec[idim] = " << m_ncells_vec[idim]
<< std::endl;);
}
return cell_vec;
}
void position_to_local_ind(const double *const x, size_t &icell,
size_t &isubdom) const {
// note: the speed of this function is important because it is called
// once per atom each time update_container is called.
cell_vec_t cell_vec = position_to_cell_vec(x);
isubdom = get_subdomain(cell_vec);
icell = cell_vec_to_local_ind(cell_vec, isubdom);
}
void position_to_global_ind(const double *const x, size_t &icell) const {
// note: the speed of this function is important because it is called
// once per atom each time update_container is called.
cell_vec_t cell_vec = position_to_cell_vec(x);
icell = cell_vec_to_global_ind(cell_vec);
}
void global_ind_to_local_ind(const size_t icell_global, size_t &icell_local,
size_t &isubdom) const {
cell_vec_t cell_vec = global_ind_to_cell_vec(icell_global);
isubdom = get_subdomain(cell_vec);
icell_local = cell_vec_to_local_ind(cell_vec, isubdom);
}
size_t local_ind_to_global_ind(const size_t icell_local,
const size_t isubdom) const {
cell_vec_t cell_vec = local_ind_to_cell_vec(icell_local, isubdom);
return cell_vec_to_global_ind(cell_vec);
}
size_t get_subdomain(cell_vec_t const &v) const {
if (m_subdoms_balanced) {
return v[1] / m_subdom_avg_len;
} else {
return std::upper_bound(m_subdom_limits.begin(), m_subdom_limits.end(),
v[1]) -
m_subdom_limits.begin() - 1;
}
}
size_t get_subdomain(const size_t icell_global) const {
cell_vec_t cell_vec = global_ind_to_cell_vec(icell_global);
return get_subdomain(cell_vec);
}
double minimum_distance(cell_vec_t const &v1, cell_vec_t const &v2) const {
// copy them so we don't accidentally change them
auto lower_left1 = cell_vec_to_position(v1);
auto lower_left2 = cell_vec_to_position(v2);
pele::VecN<ndim> ll1, ll2, dr;
pele::VecN<ndim>
minimum_dist; // the minimum possible distance in each direction
for (size_t i = 0; i < ndim; ++i) {
double min_dist = std::numeric_limits<double>::max();
double dri;
double rcell = m_rcell_vec[i];
// find the minimum distance in the i'th direction.
ll1 = lower_left1;
ll2 = lower_left2;
m_dist->get_rij(dr.data(), ll1.data(), ll2.data());
dri = std::abs(dr[i]);
if (dri < min_dist) {
min_dist = dri;
}
ll1 = lower_left1;
ll2 = lower_left2;
ll1[i] += rcell;
m_dist->get_rij(dr.data(), ll1.data(), ll2.data());
dri = std::abs(dr[i]);
if (dri < min_dist) {
min_dist = dri;
}
ll1 = lower_left1;
ll2 = lower_left2;
ll2[i] += rcell;
m_dist->get_rij(dr.data(), ll1.data(), ll2.data());
dri = std::abs(dr[i]);
if (dri < min_dist) {
min_dist = dri;
}
ll1 = lower_left1;
ll2 = lower_left2;
ll1[i] += rcell;
ll2[i] += rcell;
m_dist->get_rij(dr.data(), ll1.data(), ll2.data());
dri = std::abs(dr[i]);
if (dri < min_dist) {
min_dist = dri;
}
minimum_dist[i] = min_dist;
}
double r2_min = dot<ndim>(minimum_dist, minimum_dist);
return std::sqrt(r2_min);
}
void find_global_neighbor_inds(const size_t idim, cell_vec_t const &v0,
std::vector<size_t> &neighbors,
cell_vec_t const &vorigin) const {
if (idim == ndim) {
double rmin = minimum_distance(v0, vorigin);
if (rmin <= m_rcut) {
neighbors.push_back(cell_vec_to_global_ind(v0));
}
} else {
auto v = v0;
for (v[idim] = 0; v[idim] < m_ncells_vec[idim]; ++v[idim]) {
find_global_neighbor_inds(idim + 1, v, neighbors, vorigin);
}
}
}
void find_global_neighbor_inds_nearest(const size_t idim,
cell_vec_t const &possible_neighbor,
std::vector<size_t> &neighbors,
cell_vec_t const &cell) const {
if (idim >= ndim) {
double rmin = minimum_distance(possible_neighbor, cell);
if (rmin <= m_rcut) {
size_t global_ind = cell_vec_to_global_ind(possible_neighbor);
auto found = std::find(neighbors.begin(), neighbors.end(), global_ind);
if (found == neighbors.end()) {
neighbors.push_back(global_ind);
}
}
} else {
auto v = possible_neighbor;
int search_range = std::ceil(m_rcut / m_rcell_vec[idim]);
// A little extra distance to make sure we cover corner to corner cases
for (int delta = -2 * search_range; delta <= 2 * search_range; ++delta) {
// v is size_t, so we need to be careful about negative values
int m_n_cells = m_ncells_vec[idim];
int expected = cell[idim] + delta;
expected = (expected % m_n_cells + m_n_cells) % m_n_cells;
v[idim] = expected;
find_global_neighbor_inds_nearest(idim + 1, v, neighbors, cell);
}
}
}
std::vector<size_t> find_all_global_neighbor_inds(const size_t icell) const {
cell_vec_t vcell = global_ind_to_cell_vec(icell);
std::vector<size_t> neighbors;
if (typeid(*m_dist) == typeid(periodic_distance<2>) or
typeid(*m_dist) == typeid(periodic_distance<3>)) {
// For a regular lattice does not iterate over all the cells
// instead only the nearest ones
find_global_neighbor_inds_nearest(0, vcell, neighbors, vcell);
return neighbors;
}
find_global_neighbor_inds(0, vcell, neighbors, vcell);
return neighbors;
}
void find_neighbor_pairs(
std::vector<std::vector<std::array<cell_t *, 2>>> &cell_neighbors_inner,
std::vector<std::vector<std::array<cell_t *, 2>>>
&cell_neighbors_boundary,
std::vector<std::vector<cell_t *>> &cell_neighbors,
std::vector<std::vector<cell_t>> &cells) const {
#ifdef _OPENMP
#pragma omp parallel
{
size_t isubdom = omp_get_thread_num();
find_neighbor_pairs_subdom(cell_neighbors_inner[isubdom],
cell_neighbors_boundary[isubdom],
cell_neighbors, cells, isubdom);
}
#else
for (size_t isubdom = 0; isubdom < m_nsubdoms; ++isubdom) {
find_neighbor_pairs_subdom(cell_neighbors_inner[isubdom],
cell_neighbors_boundary[isubdom],
cell_neighbors, cells, isubdom);
}
#endif
}
void find_neighbor_pairs_subdom(
std::vector<std::array<cell_t *, 2>> &neighbor_pairs_inner,
std::vector<std::array<cell_t *, 2>> &neighbor_pairs_boundary,
std::vector<std::vector<cell_t *>> &cell_neighbors,
std::vector<std::vector<cell_t>> &cells, const size_t isubdom) const {
neighbor_pairs_inner = std::vector<std::array<cell_t *, 2>>();
neighbor_pairs_boundary = std::vector<std::array<cell_t *, 2>>();
// Reserve memory for the cell neighbors (only for performance)
// This calculates the exact number of cells if ncellx_scale <= 1
if (m_nsubdoms == 1) {
neighbor_pairs_inner.reserve(m_subdom_ncells[isubdom] *
((std::pow(3, ndim) + 1) * 0.5));
} else {
const size_t subdom_length =
m_subdom_limits[isubdom + 1] - m_subdom_limits[isubdom];
const size_t surface_ncells = m_subdom_ncells[isubdom] / subdom_length;
if (isubdom == 0) {
// For Lees-Edwards Boundary Conditions
neighbor_pairs_boundary.reserve(surface_ncells * 4 *
std::pow(3, ndim - 2));
} else {
neighbor_pairs_boundary.reserve(surface_ncells * std::pow(3, ndim - 1));
}
neighbor_pairs_inner.reserve(m_subdom_ncells[isubdom] *
((std::pow(3, ndim) + 1) * 0.5) -
surface_ncells * std::pow(3, ndim - 1));
}
for (size_t local_icell = 0; local_icell < m_subdom_ncells[isubdom];
++local_icell) {
size_t global_icell = local_ind_to_global_ind(local_icell, isubdom);
auto global_neighbors = find_all_global_neighbor_inds(global_icell);
for (const size_t global_jcell : global_neighbors) {
size_t local_jcell, jsubdom;
global_ind_to_local_ind(global_jcell, local_jcell, jsubdom);
cell_neighbors[global_icell].push_back(&cells[jsubdom][local_jcell]);
if (isubdom == jsubdom) {
if (local_jcell >= local_icell) { // avoid duplicates
std::array<cell_t *, 2> neighbors = {&cells[isubdom][local_icell],
&cells[jsubdom][local_jcell]};
neighbor_pairs_inner.push_back(neighbors);
}
} else {
if (pos_direction_y(
global_icell,
global_jcell)) { // avoid duplicates, balance load
std::array<cell_t *, 2> neighbors = {&cells[isubdom][local_icell],
&cells[jsubdom][local_jcell]};
neighbor_pairs_boundary.push_back(neighbors);
}
}
}
cell_neighbors[global_icell].shrink_to_fit();
}
neighbor_pairs_inner.shrink_to_fit();
neighbor_pairs_boundary.shrink_to_fit();
}
bool pos_direction_y(const size_t icell, const size_t jcell) const {
auto ipos = cell_vec_to_position(global_ind_to_cell_vec(icell));
auto jpos = cell_vec_to_position(global_ind_to_cell_vec(jcell));
pele::VecN<ndim> dr;
m_dist->get_rij(dr.data(), ipos.data(), jpos.data());
return dr[1] >= 0;
}
};
template <class visitor_t, typename distance_policy = periodic_distance<3>>
class CellListsLoop {
protected:
static const size_t m_ndim = distance_policy::_ndim;
typedef VecN<m_ndim, size_t> cell_vec_t;
visitor_t &m_visitor;
CellListsContainer<m_ndim> const &m_container;
LatticeNeighbors<distance_policy> m_lattice_tool;
virtual void loop_cell_pairs(
std::vector<std::array<cell_t *, 2>> const &neighbor_pairs,
const size_t isubdom) {
for (auto const &ijpair : neighbor_pairs) {
cell_t *icell = ijpair[0];
cell_t *jcell = ijpair[1];
// do double loop through atoms, avoiding duplicate pairs
for (auto iatom = icell->begin(); iatom != icell->end(); ++iatom) {
// if icell==jcell we need to avoid duplicate atom pairs
auto jend = (icell == jcell) ? iatom : jcell->end();
for (auto jatom = jcell->begin(); jatom != jend; ++jatom) {
m_visitor.insert_atom_pair(*iatom, *jatom, isubdom);
}
}
}
}
virtual void loop_cell_pairs_specific(std::vector<size_t> const &icells,
std::vector<size_t> const &iatoms) {
for (size_t i = 0; i < icells.size(); ++i) {
size_t isubdom = m_lattice_tool.get_subdomain(icells[i]);
for (cell_t *jcell : m_container.m_cell_neighbors[icells[i]]) {
for (auto jatom = jcell->begin(); jatom != jcell->end(); ++jatom) {
if (iatoms[i] != *jatom) {
m_visitor.insert_atom_pair(iatoms[i], *jatom, isubdom);
}
}
}
}
}
virtual void loop_through_cell_neighbors_for_unique_pairs(
size_t cell_containing_iatom, size_t iatom,
std::vector<size_t> const &iatoms) {
size_t isubdom = m_lattice_tool.get_subdomain(cell_containing_iatom);
bool atom_in_iatoms = false;
for (cell_t *jcell : m_container.m_cell_neighbors[cell_containing_iatom]) {
for (auto jatom = jcell->begin(); jatom != jcell->end(); ++jatom) {
if (iatom != *jatom) {
atom_in_iatoms =
std::find(iatoms.begin(), iatoms.end(), *jatom) != iatoms.end();
// Conditional to avoid double counting
if ((atom_in_iatoms && iatom < *jatom) || (!atom_in_iatoms)) {
m_visitor.insert_atom_pair(iatom, *jatom, isubdom);
}
}
}
}
}
public:
CellListsLoop(visitor_t &visitor, CellListsContainer<m_ndim> const &container,
pele::LatticeNeighbors<distance_policy> lattice_tool)
: m_visitor(visitor),
m_container(container),
m_lattice_tool(lattice_tool) {}
void loop_through_atom_pairs() {
#ifdef _OPENMP
#pragma omp parallel
{
size_t isubdom = omp_get_thread_num();
loop_cell_pairs(m_container.m_cell_neighbor_pairs_inner[isubdom],
isubdom);
#pragma omp barrier
loop_cell_pairs(m_container.m_cell_neighbor_pairs_boundary[isubdom],
isubdom);
}
#else
size_t nsubdoms = m_container.m_cell_atoms.size();
for (size_t isubdom = 0; isubdom < nsubdoms; ++isubdom) {
loop_cell_pairs(m_container.m_cell_neighbor_pairs_inner[isubdom],
isubdom);
loop_cell_pairs(m_container.m_cell_neighbor_pairs_boundary[isubdom],
isubdom);
}
#endif
}
void loop_through_atom_pairs_specific(pele::Array<double> const &coords,
std::vector<size_t> const &iatoms) {
std::vector<size_t> icells = std::vector<size_t>(iatoms.size());
for (size_t i = 0; i < iatoms.size(); ++i) {
m_lattice_tool.position_to_global_ind(coords.data() + m_ndim * iatoms[i],
icells[i]);
}
loop_cell_pairs_specific(icells, iatoms);
}
void loop_through_unique_pairs_containing_iatoms(
pele::Array<double> const &coords, std::vector<size_t> const &iatoms) {
size_t cell_containing_iatom;
for (size_t i = 0; i < iatoms.size(); ++i) {
m_lattice_tool.position_to_global_ind(coords.data() + m_ndim * iatoms[i],
cell_containing_iatom);
loop_through_cell_neighbors_for_unique_pairs(cell_containing_iatom,
iatoms[i], iatoms);
}
}
};
template <typename distance_policy = periodic_distance<3>>
class CellLists {
public:
static const size_t m_ndim = distance_policy::_ndim;
protected:
bool m_initialized; // flag for whether the cell lists have been
// initialized with coordinates
pele::LatticeNeighbors<distance_policy> m_lattice_tool;
std::vector<SafePushQueue<std::array<size_t, 2>>> add_atom_queue;
CellListsContainer<m_ndim> m_container;
public:
~CellLists() {}
CellLists(std::shared_ptr<distance_policy> const &dist,
pele::Array<double> const &boxv, const double rcut,
const double ncellx_scale = 1.0, const bool balance_omp = true);
template <class callback_class>
inline CellListsLoop<callback_class, distance_policy> get_atom_pair_looper(
callback_class &callback) const {
return CellListsLoop<callback_class, distance_policy>(callback, m_container,
m_lattice_tool);
}
size_t get_nr_cells() const { return m_lattice_tool.m_ncells; }
size_t get_nr_cellsx() const { return m_lattice_tool.m_ncells_vec[0]; }
pele::Array<size_t> get_order(const size_t natoms) const {
auto order = pele::Array<size_t>(natoms);
size_t ind = 0;
for (auto const &cell_atoms : m_container.m_cell_atoms) {
for (auto const &atoms : cell_atoms) {
for (auto const &iatom : atoms) {
order[ind] = iatom;
ind++;
}
}
}
return order;
}
void update(pele::Array<double> const &coords);
void update_specific(pele::Array<double> const &coords,
std::vector<size_t> const &iatoms,
std::vector<double> const &old_coords);
protected:
void print_warnings(const size_t natoms);
void build_cell_neighbors_list();
void reset_container(pele::Array<double> const &coords);
void update_container(pele::Array<double> const &coords);
void update_container_specific(pele::Array<double> const &coords,
std::vector<size_t> const &iatoms,
std::vector<double> const &old_coords);
private:
static Array<size_t> get_ncells_vec(Array<double> const &boxv,
const double rcut,
const double ncellx_scale,
const bool balance_omp);
};
template <typename distance_policy>
CellLists<distance_policy>::CellLists(
std::shared_ptr<distance_policy> const &dist,
pele::Array<double> const &boxv, const double rcut,
const double ncellx_scale, const bool balance_omp)
: m_initialized(false),
m_lattice_tool(dist, boxv, rcut,
get_ncells_vec(boxv, rcut, ncellx_scale, balance_omp)),
m_container(m_lattice_tool.m_subdom_ncells),
add_atom_queue(m_lattice_tool.m_nsubdoms) {
build_cell_neighbors_list();
if (m_ndim < 2) {
throw std::runtime_error(
"CellLists::CellLists: Cell lists only work with more than 1 "
"dimension "
"(due to the split of subdomains in y-dimension)");
}
if (boxv.size() != m_ndim) {
throw std::runtime_error(
"CellLists::CellLists: distance policy boxv and "
"cell list boxv differ in size");
}
if (*std::min_element(boxv.begin(), boxv.end()) < rcut) {
throw std::runtime_error("CellLists::CellLists: illegal rcut");
}
for (size_t i = 1; i < boxv.size(); ++i) {
if (boxv[i] < 0) {
throw std::runtime_error(
"CellLists::CellLists: illegal input: boxvector");
}
}
if (ncellx_scale < 0) {
throw std::runtime_error(
"CellLists::CellLists: illegal input: ncellx_scale");
}
#ifdef _OPENMP
if (std::floor(ncellx_scale) != ncellx_scale && ncellx_scale > 1 &&
omp_get_max_threads() > 1) {
throw std::runtime_error(
"CellLists::CellLists: Non-integer values > 1 of "
"ncellx_scale can break the parallelization!");
}
#endif
size_t ncell_min = *std::min_element(m_lattice_tool.m_ncells_vec.begin(),
m_lattice_tool.m_ncells_vec.end());
if ((m_ndim == 2 && ncell_min <= 4) || (m_ndim == 3 && ncell_min <= 3)) {
// If there are only a few cells in any direction then it doesn't make
// sense to use cell lists because so many cells will be neighbors with
// each other. It would be better to use simple loops over atom pairs.
std::cout << "CellLists: efficiency warning: there are not many cells ("
<< ncell_min << ") in at least one direction.\n";
}
if (m_lattice_tool.m_rcut >
0.5 * *std::min_element(m_lattice_tool.m_boxvec.begin(),
m_lattice_tool.m_boxvec.end())) {
// an atom can interact with more than just the nearest image of it's
// neighbor
std::cerr << "CellLists: warning: rcut > half the box length. This might "
"cause errors with periodic boundaries.\n";
}
}
template <typename distance_policy>
Array<size_t> CellLists<distance_policy>::get_ncells_vec(
Array<double> const &boxv, const double rcut, const double ncellx_scale,
const bool balance_omp) {
pele::Array<size_t> res(boxv.size());
for (size_t i = 0; i < res.size(); ++i) {
res[i] = std::max<size_t>(1, ncellx_scale * boxv[i] / rcut);
}
#ifdef _OPENMP
if (omp_get_max_threads() > res[1]) {
throw std::runtime_error(
"More threads than cells in y-direction. "
"Reduce the number of threads "
"(environment variable OMP_NUM_THREADS)!");
}
if (balance_omp) {
res[1] = (res[1] / omp_get_max_threads()) * omp_get_max_threads();
}
#endif
return res;
}
template <typename distance_policy>
void CellLists<distance_policy>::print_warnings(const size_t natoms) {
if (m_lattice_tool.m_ncells > natoms) {
// It would be more efficient (I think) to reduce the number of cells.
std::cout << "CellLists: efficiency warning: the number of cells ("
<< m_lattice_tool.m_ncells << ")"
<< " is greater than the number of atoms (" << natoms << ").\n";
}
}
template <typename distance_policy>
void CellLists<distance_policy>::update(pele::Array<double> const &coords) {
if (m_initialized) {
if (m_lattice_tool.m_nsubdoms == 1) {
reset_container(coords);
} else {
update_container(coords);
}
} else {
m_initialized = true;
const size_t natoms = coords.size() / m_ndim;
print_warnings(natoms);
reset_container(coords);
}
}
template <typename distance_policy>
void CellLists<distance_policy>::update_specific(
pele::Array<double> const &coords, std::vector<size_t> const &iatoms,
std::vector<double> const &old_coords) {
if (m_initialized) {
update_container_specific(coords, iatoms, old_coords);
} else {
m_initialized = true;
size_t natoms = coords.size() / m_ndim;
print_warnings(natoms);
reset_container(coords);
}
}
template <typename distance_policy>
void CellLists<distance_policy>::build_cell_neighbors_list() {
m_lattice_tool.find_neighbor_pairs(m_container.m_cell_neighbor_pairs_inner,
m_container.m_cell_neighbor_pairs_boundary,
m_container.m_cell_neighbors,
m_container.m_cell_atoms);
}
template <typename distance_policy>
void CellLists<distance_policy>::reset_container(
pele::Array<double> const &coords) {
m_container.reset();
const size_t natoms = coords.size() / m_ndim;
// std::cout << "CellLists: resetting container with " << natoms << "
// atoms\n";
for (size_t iatom = 0; iatom < natoms; ++iatom) {
const double *const x = coords.data() + m_ndim * iatom;
size_t icell, isubdom;
m_lattice_tool.position_to_local_ind(x, icell, isubdom);
m_container.add_atom_to_cell(iatom, icell, isubdom);
}
}
template <typename distance_policy>
void CellLists<distance_policy>::update_container(
pele::Array<double> const &coords) {
#ifdef _OPENMP
#pragma omp parallel
{
size_t isubdom = omp_get_thread_num();
for (size_t icell = 0; icell < m_lattice_tool.m_subdom_ncells[isubdom];
++icell) {
size_t atom_nr = 0;
while (atom_nr < m_container.m_cell_atoms[isubdom][icell].size()) {
size_t iatom = m_container.m_cell_atoms[isubdom][icell][atom_nr];
const double *const new_x = coords.data() + m_ndim * iatom;
size_t new_cell, new_subdom;
m_lattice_tool.position_to_local_ind(new_x, new_cell, new_subdom);
if (new_cell != icell || new_subdom != isubdom) {
m_container.remove_atom_from_cell_ind(atom_nr, icell, isubdom);
if (isubdom == new_subdom) {
m_container.add_atom_to_cell(iatom, new_cell, isubdom);
} else {
std::array<size_t, 2> add_info = {iatom, new_cell};
add_atom_queue[new_subdom].push(add_info);
}
} else {
++atom_nr;
}
}
}
#pragma omp barrier
while (!add_atom_queue[isubdom].empty()) {
auto add_info = add_atom_queue[isubdom].front();
add_atom_queue[isubdom].pop();
m_container.add_atom_to_cell(add_info[0], add_info[1], isubdom);
}
}
#else
for (size_t isubdom = 0; isubdom < m_lattice_tool.m_nsubdoms; ++isubdom) {
for (size_t icell = 0; icell < m_lattice_tool.m_subdom_ncells[isubdom];
++icell) {
size_t atom_nr = 0;
while (atom_nr < m_container.m_cell_atoms[isubdom][icell].size()) {
size_t iatom = m_container.m_cell_atoms[isubdom][icell][atom_nr];
const double *const new_x = coords.data() + m_ndim * iatom;
size_t new_cell, new_subdom;
m_lattice_tool.position_to_local_ind(new_x, new_cell, new_subdom);
if (new_cell != icell || new_subdom != isubdom) {
m_container.remove_atom_from_cell_ind(atom_nr, icell, isubdom);
m_container.add_atom_to_cell(iatom, new_cell, new_subdom);
} else {
++atom_nr;
}
}
}
}
#endif
}
template <typename distance_policy>
void CellLists<distance_policy>::update_container_specific(
pele::Array<double> const &coords, std::vector<size_t> const &iatoms,
std::vector<double> const &old_coords) {
for (size_t iatom = 0; iatom < iatoms.size(); iatom++) {
size_t old_cell, old_subdom;
m_lattice_tool.position_to_local_ind(old_coords.data() + m_ndim * iatom,
old_cell, old_subdom);
size_t new_cell, new_subdom;
m_lattice_tool.position_to_local_ind(coords.data() + m_ndim * iatoms[iatom],
new_cell, new_subdom);
if (new_cell != old_cell || new_subdom != old_subdom) {
m_container.remove_atom_from_cell(iatoms[iatom], old_cell, old_subdom);
m_container.add_atom_to_cell(iatoms[iatom], new_cell, new_subdom);
}
}
}
} // namespace pele
#endif // #ifndef _PELE_CELL_LISTS_H_