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         654 : LinkCells::LinkCells( Communicator& cc ) :
      29         654 :   comm(cc),
      30         654 :   cutoffwasset(false),
      31         654 :   link_cutoff(0.0),
      32         654 :   ncells(3),
      33        1308 :   nstride(3)
      34             : {
      35         654 : }
      36             : 
      37         654 : void LinkCells::setCutoff( const double& lcut ) {
      38         654 :   cutoffwasset=true; link_cutoff=lcut;
      39         654 : }
      40             : 
      41         513 : double LinkCells::getCutoff() const {
      42         513 :   plumed_assert( cutoffwasset ); return link_cutoff;
      43             : }
      44             : 
      45       11969 : void LinkCells::buildCellLists( const std::vector<Vector>& pos, const std::vector<unsigned>& indices, const Pbc& pbc ) {
      46       11969 :   plumed_assert( cutoffwasset && pos.size()==indices.size() );
      47             : 
      48             :   // Must be able to check that pbcs are not nonsensical in some way?? -- GAT
      49             : 
      50             :   // Setup the pbc object by copying it from action
      51       11969 :   mypbc.setBox( pbc.getBox() );
      52             : 
      53             :   // Setup the lists
      54       11969 :   if( pos.size()!=allcells.size() ) {
      55         460 :     allcells.resize( pos.size() ); lcell_lists.resize( pos.size() );
      56             :   }
      57             : 
      58             :   {
      59             : // This is the reciprocal lattice
      60             : // notice that reciprocal.getRow(0) is a vector that is orthogonal to b and c
      61             : // This allows to use linked cells in non orthorhomic boxes
      62       11969 :     Tensor reciprocal(transpose(mypbc.getInvBox()));
      63       11969 :     ncells[0] = std::floor( 1.0/ reciprocal.getRow(0).modulo() / link_cutoff );
      64       11969 :     if( ncells[0]==0 ) ncells[0]=1;
      65       11969 :     ncells[1] = std::floor( 1.0/ reciprocal.getRow(1).modulo() / link_cutoff );
      66       11969 :     if( ncells[1]==0 ) ncells[1]=1;
      67       11969 :     ncells[2] = std::floor( 1.0/ reciprocal.getRow(2).modulo() / link_cutoff );
      68       11969 :     if( ncells[2]==0 ) ncells[2]=1;
      69             :   }
      70             :   // Setup the strides
      71       11969 :   nstride[0]=1; nstride[1]=ncells[0]; nstride[2]=ncells[0]*ncells[1];
      72             : 
      73             :   // Setup the storage for link cells
      74       11969 :   unsigned ncellstot=ncells[0]*ncells[1]*ncells[2];
      75       11969 :   if( lcell_tots.size()!=ncellstot ) {
      76         462 :     lcell_tots.resize( ncellstot ); lcell_starts.resize( ncellstot );
      77             :   }
      78             :   // Clear nlcells
      79     1543036 :   for(unsigned i=0; i<ncellstot; ++i) lcell_tots[i]=0;
      80             :   // Clear allcells
      81       11969 :   allcells.assign( allcells.size(), 0 );
      82             : 
      83             :   // Find out what cell everyone is in
      84       11969 :   unsigned rank=comm.Get_rank(), size=comm.Get_size();
      85     1195326 :   for(unsigned i=rank; i<pos.size(); i+=size) {
      86     1183357 :     allcells[i]=findCell( pos[i] );
      87     1183357 :     lcell_tots[allcells[i]]++;
      88             :   }
      89             :   // And gather all this information on every node
      90       11969 :   comm.Sum( allcells ); comm.Sum( lcell_tots );
      91             : 
      92             :   // Now prepare the link cell lists
      93             :   unsigned tot=0;
      94     1543036 :   for(unsigned i=0; i<lcell_tots.size(); ++i) { lcell_starts[i]=tot; tot+=lcell_tots[i]; lcell_tots[i]=0; }
      95       11969 :   plumed_assert( tot==pos.size() );
      96             : 
      97             :   // And setup the link cells properly
      98     1267732 :   for(unsigned j=0; j<pos.size(); ++j) {
      99     1255763 :     unsigned myind = lcell_starts[ allcells[j] ] + lcell_tots[ allcells[j] ];
     100     1255763 :     lcell_lists[ myind ] = indices[j];
     101     1255763 :     lcell_tots[allcells[j]]++;
     102             :   }
     103       11969 : }
     104             : 
     105             : #define LINKC_MIN(n) ((n<2)? 0 : -1)
     106             : #define LINKC_MAX(n) ((n<3)? 1 : 2)
     107             : #define LINKC_PBC(n,num) ((n<0)? num-1 : n%num )
     108             : 
     109      262786 : void LinkCells::addRequiredCells( const std::array<unsigned,3>& celn, unsigned& ncells_required,
     110             :                                   std::vector<unsigned>& cells_required ) const {
     111             :   unsigned nnew_cells=0;
     112     1369765 :   for(int nx=LINKC_MIN(ncells[0]); nx<LINKC_MAX(ncells[0]); ++nx) {
     113      459037 :     int xval = celn[0] + nx;
     114      459037 :     xval=LINKC_PBC(xval,ncells[0])*nstride[0];
     115     3113069 :     for(int ny=LINKC_MIN(ncells[1]); ny<LINKC_MAX(ncells[1]); ++ny) {
     116     1042887 :       int yval = celn[1] + ny;
     117     1042887 :       yval=LINKC_PBC(yval,ncells[1])*nstride[1];
     118     8298276 :       for(int nz=LINKC_MIN(ncells[2]); nz<LINKC_MAX(ncells[2]); ++nz) {
     119     2777298 :         int zval = celn[2] + nz;
     120     2777298 :         zval=LINKC_PBC(zval,ncells[2])*nstride[2];
     121             : 
     122     2777298 :         unsigned mybox=xval+yval+zval; bool added=false;
     123     2777298 :         for(unsigned k=0; k<ncells_required; ++k) {
     124           0 :           if( mybox==cells_required[k] ) { added=true; break; }
     125             :         }
     126     2777298 :         if( !added ) { cells_required[ncells_required+nnew_cells]=mybox; nnew_cells++; }
     127             :       }
     128             :     }
     129             :   }
     130      262786 :   ncells_required += nnew_cells;
     131      262786 : }
     132             : 
     133        3375 : void LinkCells::retrieveNeighboringAtoms( const Vector& pos, std::vector<unsigned>& cell_list,
     134             :     unsigned& natomsper, std::vector<unsigned>& atoms ) const {
     135        3375 :   if( cell_list.size()!=getNumberOfCells() ) cell_list.resize( getNumberOfCells() );
     136        3375 :   unsigned ncellt=0; addRequiredCells( findMyCell( pos ), ncellt, cell_list );
     137        3375 :   retrieveAtomsInCells( ncellt, cell_list, natomsper, atoms );
     138        3375 : }
     139             : 
     140      262786 : void LinkCells::retrieveAtomsInCells( const unsigned& ncells_required,
     141             :                                       const std::vector<unsigned>& cells_required,
     142             :                                       unsigned& natomsper, std::vector<unsigned>& atoms ) const {
     143     3040084 :   for(unsigned i=0; i<ncells_required; ++i) {
     144     2777298 :     unsigned mybox=cells_required[i];
     145    26085213 :     for(unsigned k=0; k<lcell_tots[mybox]; ++k) {
     146    23307915 :       unsigned myatom = lcell_lists[lcell_starts[mybox]+k];
     147    23307915 :       if( myatom!=atoms[0] ) { // Ideally would provide an option to not do this
     148    23054122 :         atoms[natomsper]=myatom;
     149    23054122 :         natomsper++;
     150             :       }
     151             :     }
     152             :   }
     153      262786 : }
     154             : 
     155     1446143 : std::array<unsigned,3> LinkCells::findMyCell( const Vector& pos ) const {
     156     1446143 :   Vector fpos=mypbc.realToScaled( pos );
     157             :   std::array<unsigned,3> celn;
     158     5784572 :   for(unsigned j=0; j<3; ++j) {
     159     4338429 :     celn[j] = std::floor( ( Tools::pbc(fpos[j]) + 0.5 ) * ncells[j] );
     160     4338429 :     plumed_assert( celn[j]>=0 && celn[j]<ncells[j] ); // Check that atom is in box
     161             :   }
     162     1446143 :   return celn;
     163             : }
     164             : 
     165     1183357 : unsigned LinkCells::convertIndicesToIndex( const unsigned& nx, const unsigned& ny, const unsigned& nz ) const {
     166     1183357 :   return nx*nstride[0] + ny*nstride[1] + nz*nstride[2];
     167             : }
     168             : 
     169     1183357 : unsigned LinkCells::findCell( const Vector& pos ) const {
     170     1183357 :   std::array<unsigned,3> celn( findMyCell(pos ) );
     171     1183357 :   return convertIndicesToIndex( celn[0], celn[1], celn[2] );
     172             : }
     173             : 
     174       11750 : unsigned LinkCells::getMaxInCell() const {
     175       11750 :   unsigned maxn = lcell_tots[0];
     176     1527421 :   for(unsigned i=1; i<lcell_tots.size(); ++i) {
     177     1515671 :     if( lcell_tots[i]>maxn ) { maxn=lcell_tots[i]; }
     178             :   }
     179       11750 :   return maxn;
     180             : }
     181             : 
     182             : }