Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2014-2023 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "LinkCells.h"
      23             : #include "Communicator.h"
      24             : #include "Tools.h"
      25             : 
      26             : namespace PLMD {
      27             : 
      28         703 : LinkCells::LinkCells( Communicator& cc ) :
      29         703 :   comm(cc),
      30         703 :   cutoffwasset(false),
      31         703 :   link_cutoff(0.0),
      32         703 :   ncells(3),
      33        1406 :   nstride(3) {
      34         703 : }
      35             : 
      36         385 : void LinkCells::setCutoff( const double& lcut ) {
      37         385 :   cutoffwasset=true;
      38         385 :   link_cutoff=lcut;
      39         385 : }
      40             : 
      41           8 : double LinkCells::getCutoff() const {
      42           8 :   plumed_assert( cutoffwasset );
      43           8 :   return link_cutoff;
      44             : }
      45             : 
      46        1295 : void LinkCells::buildCellLists( const std::vector<Vector>& pos, const std::vector<unsigned>& indices, const Pbc& pbc ) {
      47        1295 :   plumed_assert( cutoffwasset && pos.size()==indices.size() );
      48             : 
      49             :   // Must be able to check that pbcs are not nonsensical in some way?? -- GAT
      50        1295 :   auto box = pbc.getBox();
      51        1295 :   if(box(0,0)==0.0 && box(0,1)==0.0 && box(0,2)==0.0 && box(1,0)==0.0 && box(1,1)==0.0 && box(1,2)==0.0 && box(2,0)==0.0 && box(2,1)==0 && box(2,2)==0) {
      52             :     // If the box is not set then we can't use link cells.  We thus set the link cell cutoff and box vectors equal to 23 (because it is the best number).
      53             :     // Setting everything this way ensures that the link cells are a 1x1x1 box.  Notice that if it is a one by one by one box then we are hard coded to return
      54             :     // 0 in findCell. GAT
      55             :     // Note: any non infinite - non zero number should work correctly here! Apparently Gareth likes numerology. I do as well, so let's keep this. GB
      56           4 :     box(0,0) = box(1,1) = box(2,2) = link_cutoff = 23;
      57             :   } else {
      58        1291 :     auto determinant = box.determinant();
      59        1291 :     plumed_assert(determinant > epsilon) <<"Cell lists cannot be built when passing a box with null volume. Volume is "<<determinant;
      60             :   }
      61             :   // Setup the pbc object by copying it from action
      62        1295 :   mypbc.setBox( box );
      63             : 
      64             :   // Setup the lists
      65        1295 :   if( pos.size()!=allcells.size() ) {
      66         313 :     allcells.resize( pos.size() );
      67         313 :     lcell_lists.resize( pos.size() );
      68             :   }
      69             : 
      70             :   {
      71             : // This is the reciprocal lattice
      72             : // notice that reciprocal.getRow(0) is a vector that is orthogonal to b and c
      73             : // This allows to use linked cells in non orthorhomic boxes
      74        1295 :     Tensor reciprocal(transpose(mypbc.getInvBox()));
      75        1295 :     ncells[0] = std::floor( 1.0/ reciprocal.getRow(0).modulo() / link_cutoff );
      76        1295 :     if( ncells[0]==0 ) {
      77         991 :       ncells[0]=1;
      78             :     }
      79        1295 :     ncells[1] = std::floor( 1.0/ reciprocal.getRow(1).modulo() / link_cutoff );
      80        1295 :     if( ncells[1]==0 ) {
      81         991 :       ncells[1]=1;
      82             :     }
      83        1295 :     ncells[2] = std::floor( 1.0/ reciprocal.getRow(2).modulo() / link_cutoff );
      84        1295 :     if( ncells[2]==0 ) {
      85         991 :       ncells[2]=1;
      86             :     }
      87             :   }
      88             :   // Setup the strides
      89        1295 :   nstride[0]=1;
      90        1295 :   nstride[1]=ncells[0];
      91        1295 :   nstride[2]=ncells[0]*ncells[1];
      92             : 
      93             :   // Setup the storage for link cells
      94        1295 :   unsigned ncellstot=ncells[0]*ncells[1]*ncells[2];
      95        1295 :   if( lcell_tots.size()!=ncellstot ) {
      96         305 :     lcell_tots.resize( ncellstot );
      97         305 :     lcell_starts.resize( ncellstot );
      98             :   }
      99             :   // Clear nlcells
     100        1295 :   lcell_tots.assign( lcell_tots.size(), 0 );
     101             :   // Clear allcells
     102        1295 :   allcells.assign( allcells.size(), 0 );
     103             : 
     104             :   // Find out what cell everyone is in
     105        1295 :   unsigned rank=comm.Get_rank(), size=comm.Get_size();
     106      411133 :   for(unsigned i=rank; i<pos.size(); i+=size) {
     107      409838 :     allcells[i]=findCell( pos[i] );
     108      409838 :     lcell_tots[allcells[i]]++;
     109             :   }
     110             :   // And gather all this information on every node
     111        1295 :   comm.Sum( allcells );
     112        1295 :   comm.Sum( lcell_tots );
     113             : 
     114             :   // Now prepare the link cell lists
     115        1295 :   unsigned tot=0;
     116      299535 :   for(unsigned i=0; i<lcell_tots.size(); ++i) {
     117      298240 :     lcell_starts[i]=tot;
     118      298240 :     tot+=lcell_tots[i];
     119      298240 :     lcell_tots[i]=0;
     120             :   }
     121        1295 :   plumed_assert( tot==pos.size() ) <<"Total number of atoms found in link cells is "<<tot<<" number of atoms is "<<pos.size();
     122             : 
     123             :   // And setup the link cells properly
     124      456412 :   for(unsigned j=0; j<pos.size(); ++j) {
     125      455117 :     unsigned myind = lcell_starts[ allcells[j] ] + lcell_tots[ allcells[j] ];
     126      455117 :     lcell_lists[ myind ] = indices[j];
     127      455117 :     lcell_tots[allcells[j]]++;
     128             :   }
     129        1295 : }
     130             : 
     131             : #define LINKC_MIN(n) ((n<2)? 0 : -1)
     132             : #define LINKC_MAX(n) ((n<3)? 1 : 2)
     133             : #define LINKC_PBC(n,num) ((n<0)? num-1 : n%num )
     134             : 
     135      228625 : void LinkCells::addRequiredCells( const std::array<unsigned,3>& celn, unsigned& ncells_required,
     136             :                                   std::vector<unsigned>& cells_required ) const {
     137             :   unsigned nnew_cells=0;
     138     1349295 :   for(int nx=LINKC_MIN(ncells[0]); nx<LINKC_MAX(ncells[0]); ++nx) {
     139      450277 :     int xval = celn[0] + nx;
     140      450277 :     xval=LINKC_PBC(xval,ncells[0])*nstride[0];
     141     3335733 :     for(int ny=LINKC_MIN(ncells[1]); ny<LINKC_MAX(ncells[1]); ++ny) {
     142     1113235 :       int yval = celn[1] + ny;
     143     1113235 :       yval=LINKC_PBC(yval,ncells[1])*nstride[1];
     144     9261962 :       for(int nz=LINKC_MIN(ncells[2]); nz<LINKC_MAX(ncells[2]); ++nz) {
     145     3090780 :         int zval = celn[2] + nz;
     146     3090780 :         zval=LINKC_PBC(zval,ncells[2])*nstride[2];
     147             : 
     148     3090780 :         unsigned mybox=xval+yval+zval;
     149             :         bool added=false;
     150     3090780 :         for(unsigned k=0; k<ncells_required; ++k) {
     151           0 :           if( mybox==cells_required[k] ) {
     152             :             added=true;
     153             :             break;
     154             :           }
     155             :         }
     156     3090780 :         if( !added ) {
     157     3090780 :           cells_required[ncells_required+nnew_cells]=mybox;
     158     3090780 :           nnew_cells++;
     159             :         }
     160             :       }
     161             :     }
     162             :   }
     163      228625 :   ncells_required += nnew_cells;
     164      228625 : }
     165             : 
     166        6576 : void LinkCells::retrieveNeighboringAtoms( const Vector& pos, std::vector<unsigned>& cell_list,
     167             :     unsigned& natomsper, std::vector<unsigned>& atoms ) const {
     168        6576 :   if( cell_list.size()!=getNumberOfCells() ) {
     169         304 :     cell_list.resize( getNumberOfCells() );
     170             :   }
     171        6576 :   unsigned ncellt=0;
     172        6576 :   addRequiredCells( findMyCell( pos ), ncellt, cell_list );
     173        6576 :   retrieveAtomsInCells( ncellt, cell_list, natomsper, atoms );
     174        6576 : }
     175             : 
     176      228625 : void LinkCells::retrieveAtomsInCells( const unsigned& ncells_required,
     177             :                                       const std::vector<unsigned>& cells_required,
     178             :                                       unsigned& natomsper, std::vector<unsigned>& atoms ) const {
     179      228625 :   plumed_assert( natomsper==1 || natomsper==2 );  // This is really a bug. If you are trying to reuse this ask GAT for help
     180     3319405 :   for(unsigned i=0; i<ncells_required; ++i) {
     181     3090780 :     unsigned mybox=cells_required[i];
     182   405180964 :     for(unsigned k=0; k<lcell_tots[mybox]; ++k) {
     183   402090184 :       unsigned myatom = lcell_lists[lcell_starts[mybox]+k];
     184   402090184 :       if( myatom!=atoms[0] ) { // Ideally would provide an option to not do this
     185   401907298 :         atoms[natomsper]=myatom;
     186   401907298 :         natomsper++;
     187             :       }
     188             :     }
     189             :   }
     190      228625 : }
     191             : 
     192      638463 : std::array<unsigned,3> LinkCells::findMyCell( const Vector& pos ) const {
     193             :   std::array<unsigned,3> celn;
     194      638463 :   if( ncells[0]*ncells[1]*ncells[2] == 1 ) {
     195      178864 :     celn[0]=celn[1]=celn[2]=0;
     196      178864 :     return celn;
     197             :   }
     198      459599 :   Vector fpos=mypbc.realToScaled( pos );
     199     1838396 :   for(unsigned j=0; j<3; ++j) {
     200     1378797 :     celn[j] = std::floor( ( Tools::pbc(fpos[j]) + 0.5 ) * ncells[j] );
     201     1378797 :     plumed_assert( celn[j]>=0 && celn[j]<ncells[j] ); // Check that atom is in box
     202             :   }
     203      459599 :   return celn;
     204             : }
     205             : 
     206      409838 : unsigned LinkCells::convertIndicesToIndex( const unsigned& nx, const unsigned& ny, const unsigned& nz ) const {
     207      409838 :   return nx*nstride[0] + ny*nstride[1] + nz*nstride[2];
     208             : }
     209             : 
     210      409838 : unsigned LinkCells::findCell( const Vector& pos ) const {
     211      409838 :   std::array<unsigned,3> celn( findMyCell(pos ) );
     212      409838 :   return convertIndicesToIndex( celn[0], celn[1], celn[2] );
     213             : }
     214             : 
     215             : 
     216             : }