libhoomd/computes/NeighborListBinned.cc

   1 /*
   2 Highly Optimized Object-oriented Many-particle Dynamics -- Blue Edition
   3 (HOOMD-blue) Open Source Software License Copyright 2009-2014 The Regents of
   4 the University of Michigan All rights reserved.
   5
   6 HOOMD-blue may contain modifications ("Contributions") provided, and to which
   7 copyright is held, by various Contributors who have granted The Regents of the
   8 University of Michigan the right to modify and/or distribute such Contributions.
   9
  10 You may redistribute, use, and create derivate works of HOOMD-blue, in source
  11 and binary forms, provided you abide by the following conditions:
  12
  13 * Redistributions of source code must retain the above copyright notice, this
  14 list of conditions, and the following disclaimer both in the code and
  15 prominently in any materials provided with the distribution.
  16
  17 * Redistributions in binary form must reproduce the above copyright notice, this
  18 list of conditions, and the following disclaimer in the documentation and/or
  19 other materials provided with the distribution.
  20
  21 * All publications and presentations based on HOOMD-blue, including any reports
  22 or published results obtained, in whole or in part, with HOOMD-blue, will
  23 acknowledge its use according to the terms posted at the time of submission on:
  24 http://codeblue.umich.edu/hoomd-blue/citations.html
  25
  26 * Any electronic documents citing HOOMD-Blue will link to the HOOMD-Blue website:
  27 http://codeblue.umich.edu/hoomd-blue/
  28
  29 * Apart from the above required attributions, neither the name of the copyright
  30 holder nor the names of HOOMD-blue's contributors may be used to endorse or
  31 promote products derived from this software without specific prior written
  32 permission.
  33
  34 Disclaimer
  35
  36 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS ``AS IS'' AND
  37 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  38 WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, AND/OR ANY
  39 WARRANTIES THAT THIS SOFTWARE IS FREE OF INFRINGEMENT ARE DISCLAIMED.
  40
  41 IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
  42 INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  43 BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  44 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  45 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
  46 OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
  47 ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  48 */
  49
  50 // Maintainer: joaander
  51
  52 /*! \file NeighborListBinned.cc
  53     \brief Defines NeighborListBinned
  54 */
  55
  56 #include "NeighborListBinned.h"
  57
  58 #include <boost/python.hpp>
  59 using namespace boost::python;
  60
  61 #ifdef ENABLE_MPI
  62 #include "Communicator.h"
  63 #endif
  64
  65 NeighborListBinned::NeighborListBinned(boost::shared_ptr<SystemDefinition> sysdef,
  66                                        Scalar r_cut,
  67                                        Scalar r_buff,
  68                                        boost::shared_ptr<CellList> cl)
  69     : NeighborList(sysdef, r_cut, r_buff), m_cl(cl)
  70     {
  71     m_exec_conf->msg->notice(5) << "Constructing NeighborListBinned" << endl;
  72
  73     // create a default cell list if one was not specified
  74     if (!m_cl)
  75         m_cl = boost::shared_ptr<CellList>(new CellList(sysdef));
  76
  77     m_cl->setNominalWidth(r_cut + r_buff + m_d_max - Scalar(1.0));
  78     m_cl->setRadius(1);
  79     m_cl->setComputeTDB(false);
  80     m_cl->setFlagIndex();
  81     }
  82
  83 NeighborListBinned::~NeighborListBinned()
  84     {
  85     m_exec_conf->msg->notice(5) << "Destroying NeighborListBinned" << endl;
  86
  87     }
  88
  89 void NeighborListBinned::setRCut(Scalar r_cut, Scalar r_buff)
  90     {
  91     NeighborList::setRCut(r_cut, r_buff);
  92
  93     m_cl->setNominalWidth(r_cut + r_buff + m_d_max - Scalar(1.0));
  94     }
  95
  96 void NeighborListBinned::setMaximumDiameter(Scalar d_max)
  97     {
  98     NeighborList::setMaximumDiameter(d_max);
  99
 100     // need to update the cell list settings appropriately
 101     m_cl->setNominalWidth(m_r_cut + m_r_buff + m_d_max - Scalar(1.0));
 102     }
 103
 104 void NeighborListBinned::buildNlist(unsigned int timestep)
 105     {
 106     m_cl->compute(timestep);
 107
 108     uint3 dim = m_cl->getDim();
 109     Scalar3 ghost_width = m_cl->getGhostWidth();
 110
 111     if (m_prof)
 112         m_prof->push(exec_conf, "compute");
 113
 114
 115     // acquire the particle data and box dimension
 116     ArrayHandle<Scalar4> h_pos(m_pdata->getPositions(), access_location::host, access_mode::read);
 117     ArrayHandle<unsigned int> h_body(m_pdata->getBodies(), access_location::host, access_mode::read);
 118     ArrayHandle<Scalar> h_diameter(m_pdata->getDiameters(), access_location::host, access_mode::read);
 119
 120     const BoxDim& box = m_pdata->getBox();
 121     Scalar3 nearest_plane_distance = box.getNearestPlaneDistance();
 122
 123     // start by creating a temporary copy of r_cut sqaured
 124     Scalar rmax = m_r_cut + m_r_buff;
 125     // add d_max - 1.0, if diameter filtering is not already taking care of it
 126     if (!m_filter_diameter)
 127         rmax += m_d_max - Scalar(1.0);
 128     Scalar rmaxsq = rmax*rmax;
 129
 130     if ((box.getPeriodic().x && nearest_plane_distance.x <= rmax * 2.0) ||
 131         (box.getPeriodic().y && nearest_plane_distance.y <= rmax * 2.0) ||
 132         (this->m_sysdef->getNDimensions() == 3 && box.getPeriodic().z && nearest_plane_distance.z <= rmax * 2.0))
 133         {
 134         m_exec_conf->msg->error() << "nlist: Simulation box is too small! Particles would be interacting with themselves." << endl;
 135         throw runtime_error("Error updating neighborlist bins");
 136         }
 137
 138     // access the cell list data arrays
 139     ArrayHandle<unsigned int> h_cell_size(m_cl->getCellSizeArray(), access_location::host, access_mode::read);
 140     ArrayHandle<Scalar4> h_cell_xyzf(m_cl->getXYZFArray(), access_location::host, access_mode::read);
 141     ArrayHandle<unsigned int> h_cell_adj(m_cl->getCellAdjArray(), access_location::host, access_mode::read);
 142
 143     ArrayHandle<unsigned int> h_nlist(m_nlist, access_location::host, access_mode::overwrite);
 144     ArrayHandle<unsigned int> h_n_neigh(m_n_neigh, access_location::host, access_mode::overwrite);
 145
 146     unsigned int conditions = 0;
 147
 148     // access indexers
 149     Index3D ci = m_cl->getCellIndexer();
 150     Index2D cli = m_cl->getCellListIndexer();
 151     Index2D cadji = m_cl->getCellAdjIndexer();
 152
 153     // get periodic flags
 154     uchar3 periodic = box.getPeriodic();
 155
 156     // for each local particle
 157     unsigned int nparticles = m_pdata->getN();
 158
 159     for (int i = 0; i < (int)nparticles; i++)
 160         {
 161         unsigned int cur_n_neigh = 0;
 162
 163         Scalar3 my_pos = make_scalar3(h_pos.data[i].x, h_pos.data[i].y, h_pos.data[i].z);
 164         unsigned int bodyi = h_body.data[i];
 165         Scalar di = h_diameter.data[i];
 166
 167         // find the bin each particle belongs in
 168         Scalar3 f = box.makeFraction(my_pos,ghost_width);
 169         int ib = (unsigned int)(f.x * dim.x);
 170         int jb = (unsigned int)(f.y * dim.y);
 171         int kb = (unsigned int)(f.z * dim.z);
 172
 173         // need to handle the case where the particle is exactly at the box hi
 174         if (ib == (int)dim.x && periodic.x)
 175             ib = 0;
 176         if (jb == (int)dim.y && periodic.y)
 177             jb = 0;
 178         if (kb == (int)dim.z && periodic.z)
 179             kb = 0;
 180
 181         // identify the bin
 182         unsigned int my_cell = ci(ib,jb,kb);
 183
 184         // loop through all neighboring bins
 185         for (unsigned int cur_adj = 0; cur_adj < cadji.getW(); cur_adj++)
 186             {
 187             unsigned int neigh_cell = h_cell_adj.data[cadji(cur_adj, my_cell)];
 188
 189             // check against all the particles in that neighboring bin to see if it is a neighbor
 190             unsigned int size = h_cell_size.data[neigh_cell];
 191             for (unsigned int cur_offset = 0; cur_offset < size; cur_offset++)
 192                 {
 193                 Scalar4& cur_xyzf = h_cell_xyzf.data[cli(cur_offset, neigh_cell)];
 194                 unsigned int cur_neigh = __scalar_as_int(cur_xyzf.w);
 195
 196                 Scalar3 neigh_pos = make_scalar3(cur_xyzf.x, cur_xyzf.y, cur_xyzf.z);
 197
 198                 Scalar3 dx = my_pos - neigh_pos;
 199
 200                 dx = box.minImage(dx);
 201
 202                 bool excluded = (i == (int)cur_neigh);
 203
 204                 if (m_filter_body && bodyi != NO_BODY)
 205                     excluded = excluded | (bodyi == h_body.data[cur_neigh]);
 206
 207                 Scalar sqshift = Scalar(0.0);
 208                 if (m_filter_diameter)
 209                     {
 210                     // compute the shift in radius to accept neighbors based on their diameters
 211                     Scalar delta = (di + h_diameter.data[cur_neigh]) * Scalar(0.5) - Scalar(1.0);
 212                     // r^2 < (r_max + delta)^2
 213                     // r^2 < r_maxsq + delta^2 + 2*r_max*delta
 214                     sqshift = (delta + Scalar(2.0) * rmax) * delta;
 215                     }
 216
 217                 Scalar dr_sq = dot(dx,dx);
 218
 219                 if (dr_sq <= (rmaxsq + sqshift) && !excluded)
 220                     {
 221                     if (m_storage_mode == full || i < (int)cur_neigh)
 222                         {
 223                         // local neighbor
 224                         if (cur_n_neigh < m_nlist_indexer.getH())
 225                             h_nlist.data[m_nlist_indexer(i, cur_n_neigh)] = cur_neigh;
 226                         else
 227                             conditions = max(conditions, cur_n_neigh+1);
 228
 229                         cur_n_neigh++;
 230                         }
 231                     }
 232                 }
 233             }
 234
 235         h_n_neigh.data[i] = cur_n_neigh;
 236         }
 237
 238     // write out conditions
 239     m_conditions.resetFlags(conditions);
 240
 241     if (m_prof)
 242         m_prof->pop(exec_conf);
 243     }
 244
 245 void export_NeighborListBinned()
 246     {
 247     class_<NeighborListBinned, boost::shared_ptr<NeighborListBinned>, bases<NeighborList>, boost::noncopyable >
 248                      ("NeighborListBinned", init< boost::shared_ptr<SystemDefinition>, Scalar, Scalar, boost::shared_ptr<CellList> >())
 249                      ;
 250     }