LCOV - code coverage report
Current view: top level - gridtools - GridCoordinatesObject.cpp (source / functions) Hit Total Coverage
Test: plumed test coverage Lines: 266 314 84.7 %
Date: 2026-11-05 10:19:04 Functions: 19 24 79.2 %

          Line data    Source code
       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-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 "GridCoordinatesObject.h"
      23             : #include "tools/Random.h"
      24             : #include "tools/Matrix.h"
      25             : 
      26             : namespace PLMD {
      27             : namespace gridtools {
      28             : 
      29         175 : void GridCoordinatesObject::setup( const std::string& geom, const std::vector<bool>& ipbc,
      30             :                                    const unsigned& np, const double& fib_cutoff ) {
      31         175 :   if( geom=="flat" ) {
      32         143 :     gtype=flat;
      33         143 :     dimension = ipbc.size();
      34          32 :   } else if( geom=="fibonacci" ) {
      35          32 :     gtype=fibonacci;
      36          32 :     dimension = 3;
      37             :   } else {
      38           0 :     plumed_merror( geom + " is invalid geometry type");
      39             :   }
      40             : 
      41         175 :   if( gtype==flat ) {
      42         143 :     bounds_set=false;
      43         143 :     npoints=0;
      44         143 :     pbc.resize( ipbc.size() );
      45         328 :     for(unsigned i=0; i<ipbc.size(); ++i) {
      46         185 :       pbc[i]=ipbc[i];
      47             :     }
      48          32 :   } else if( gtype==fibonacci ) {
      49          32 :     if( np==0 ) {
      50           4 :       plumed_merror("the number of points on the sphere should be at least 1");
      51             :     }
      52          30 :     bounds_set=true;
      53          30 :     root5 = sqrt(5);
      54          30 :     npoints = np;
      55          30 :     golden = ( 1 + sqrt(5) ) / 2.0;
      56          30 :     igolden = golden - 1;
      57          30 :     fib_increment = 2*pi*igolden;
      58          30 :     log_golden2 = std::log( golden*golden );
      59          30 :     fib_offset = 2 / static_cast<double>( npoints );
      60          30 :     fib_shift = fib_offset/2 - 1;
      61             : 
      62          30 :     std::vector<double> icoord( dimension ), jcoord( dimension );
      63             :     // Find minimum distance between each pair of points
      64          30 :     std::vector<unsigned> tindices( dimension );
      65          30 :     std::vector<double> mindists( npoints );
      66        7792 :     for(unsigned i=0; i<npoints; ++i) {
      67        7762 :       getFibonacciCoordinates( i, icoord );
      68        7762 :       mindists[i] = 0;
      69     2633750 :       for(unsigned j=0; j<npoints; ++j) {
      70     2625988 :         if( i==j ) {
      71        7762 :           continue ;  // Points are not neighbors to themselves
      72             :         }
      73     2618226 :         getFibonacciCoordinates( j, jcoord );
      74             :         // Calculate the dot product
      75             :         double dot=0;
      76    10472904 :         for(unsigned k=0; k<dimension; ++k) {
      77     7854678 :           dot += icoord[k]*jcoord[k];
      78             :         }
      79     2618226 :         if( dot>mindists[i] ) {
      80      133331 :           mindists[i]=dot;
      81             :         }
      82             :       }
      83             :     }
      84             :     // And now take minimum of dot products
      85          30 :     const double minDist=*std::min_element(mindists.begin(),mindists.end());
      86             :     double final_cutoff;
      87          30 :     if( fib_cutoff<-1 ) {
      88             :       final_cutoff=-1;
      89             :     } else {
      90          27 :       final_cutoff = cos( acos( fib_cutoff ) + acos( minDist ) );
      91             :     }
      92             : 
      93             :     // And now construct the neighbor list
      94          30 :     fib_nlist.resize( npoints );
      95        7792 :     for(unsigned i=0; i<npoints; ++i) {
      96        7762 :       fib_nlist[i].resize(0);
      97        7762 :       getFibonacciCoordinates( i, icoord );
      98     2633750 :       for(unsigned j=0; j<npoints; ++j) {
      99     2625988 :         if( i==j ) {
     100        7762 :           continue ;  // Points are not neighbors to themselves
     101             :         }
     102     2618226 :         getFibonacciCoordinates( j, jcoord );
     103             :         // Calculate the dot product
     104             :         double dot=0;
     105    10472904 :         for(unsigned k=0; k<dimension; ++k) {
     106     7854678 :           dot += icoord[k]*jcoord[k];
     107             :         }
     108     2618226 :         if( dot>final_cutoff ) {
     109     1347586 :           fib_nlist[i].push_back(j);
     110             :         }
     111             :       }
     112             :     }
     113             :   } else {
     114           0 :     plumed_error();
     115             :   }
     116         173 : }
     117             : 
     118         139 : void GridCoordinatesObject::setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax,
     119             :                                        const std::vector<std::size_t>& binsin, std::vector<double>& spacing ) {
     120             :   plumed_dbg_assert( smin.size()==dimension && smax.size()==dimension );
     121         139 :   plumed_assert( gtype==flat && (spacing.size()==dimension || binsin.size()==dimension) );
     122         139 :   str_min.resize( dimension );
     123         139 :   str_max.resize( dimension );
     124         139 :   nbin.resize( dimension );
     125         139 :   min.resize( dimension );
     126         139 :   max.resize( dimension );
     127         139 :   dx.resize( dimension );
     128         139 :   stride.resize( dimension );
     129             : 
     130         139 :   npoints=1;
     131         139 :   bounds_set=(smin[0]!="auto" && smax[0]!="auto");
     132         139 :   if( bounds_set ) {
     133         177 :     for(unsigned i=1; i<dimension; ++i) {
     134          38 :       if( smin[i]=="auto" || smax[i]=="auto" ) {
     135           0 :         bounds_set=false;
     136           0 :         break;
     137             :       }
     138             :     }
     139             :   }
     140         312 :   for(unsigned i=0; i<dimension; ++i) {
     141         177 :     str_min[i]=smin[i];
     142             :     str_max[i]=smax[i];
     143         177 :     if( bounds_set ) {
     144         177 :       Tools::convert( str_min[i], min[i] );
     145         177 :       Tools::convert( str_max[i], max[i] );
     146         177 :       if( max[i]<=min[i] ) {
     147           2 :         plumed_merror("GRID_MAX should be larger than GRID_MIN");
     148             :       }
     149             :     }
     150         176 :     if( spacing.size()==dimension && (spacing[i]<0 || (spacing[i]==0 && binsin.size()!=dimension)) ) {
     151           4 :       plumed_merror("GRID_SPACING should be positive");
     152             :     }
     153         174 :     if( spacing.size()==dimension && binsin.size()==dimension ) {
     154           0 :       if( spacing[i]==0 ) {
     155           0 :         nbin[i] = binsin[i];
     156           0 :       } else if( bounds_set ) {
     157           0 :         double range = max[i] - min[i];
     158           0 :         nbin[i] = std::round( range / spacing[i]);
     159           0 :         dx[i]=spacing[i];
     160             :         // This check ensures that nbins is set correctly if spacing is set the same as the number of bins
     161           0 :         if( nbin[i]!=binsin[i] ) {
     162           0 :           plumed_merror("mismatch between input spacing and input number of bins");
     163             :         }
     164             :       }
     165         174 :     } else if( binsin.size()==dimension ) {
     166         171 :       nbin[i]=binsin[i];
     167         171 :       dx[i] = ( max[i] - min[i] ) / static_cast<double>( nbin[i] );
     168           3 :     } else if( spacing.size()==dimension && bounds_set ) {
     169             :       // as in METAD: the grid ends at GRID_MAX, with a spacing not larger than GRID_SPACING
     170           3 :       nbin[i] = std::ceil(( max[i] - min[i] ) / spacing[i]);
     171           3 :       dx[i] = ( max[i] - min[i] ) / static_cast<double>( nbin[i] );
     172           0 :     } else if( bounds_set ) {
     173           0 :       plumed_error();
     174             :     }
     175         174 :     if( bounds_set && nbin[i]==0 ) {
     176           2 :       plumed_merror("GRID_BIN should be at least 1");
     177             :     }
     178         173 :     if( !pbc[i] ) {
     179          93 :       max[i] +=dx[i];
     180          93 :       nbin[i]+=1;
     181             :     }
     182         173 :     stride[i]=npoints;
     183         173 :     npoints*=nbin[i];
     184             :   }
     185         135 :   if( spacing.size()!=dimension && bounds_set ) {
     186         132 :     spacing.resize(dimension);
     187         302 :     for(unsigned i=0; i<dimension; ++i) {
     188         170 :       spacing[i]=dx[i];
     189             :     }
     190             :   }
     191         135 : }
     192             : 
     193    42387699 : unsigned GridCoordinatesObject::getIndex( const std::vector<unsigned>& indices ) const {
     194             :   plumed_dbg_assert( gtype==flat && indices.size()==dimension );
     195             :   // indices are flattended using a column-major order
     196    42387699 :   unsigned index=indices[dimension-1];
     197   122020144 :   for(unsigned i=dimension-1; i>0; --i) {
     198    79632445 :     index=index*nbin[i-1]+indices[i-1];
     199             :   }
     200    42387699 :   return index;
     201             : }
     202             : 
     203           0 : bool GridCoordinatesObject::inbounds( const std::vector<double>& point ) const {
     204           0 :   return inbounds( View<const double>( point.data(), point.size() ) );
     205             : }
     206             : 
     207           0 : bool GridCoordinatesObject::inbounds( const View<const double> point ) const {
     208           0 :   if( gtype==fibonacci ) {
     209             :     return true;
     210             :   }
     211             :   plumed_dbg_assert( bounds_set && point.size()==dimension );
     212           0 :   for(unsigned i=0; i<dimension; ++i) {
     213           0 :     if( pbc[i] ) {
     214           0 :       continue;
     215             :     }
     216           0 :     if( point[i]<min[i] || point[i]>(max[i]-dx[i]) ) {
     217             :       return false;
     218             :     }
     219             :   }
     220             :   return true;
     221             : }
     222             : 
     223           0 : void GridCoordinatesObject::getIndices( const std::vector<double>& point, std::vector<unsigned>& indices ) const {
     224           0 :   getIndices( View<const double>( point.data(), point.size() ), indices );
     225           0 : }
     226             : 
     227     1169232 : void GridCoordinatesObject::getIndices( const View<const double> point, std::vector<unsigned>& indices ) const {
     228             :   plumed_dbg_assert( gtype==flat && bounds_set && point.size()==dimension && indices.size()==dimension );
     229     3453599 :   for(unsigned i=0; i<dimension; ++i) {
     230     2284368 :     indices[i]=std::floor( (point[i] - min[i])/dx[i] );
     231     2284368 :     if( pbc[i] ) {
     232       88672 :       indices[i]=indices[i]%nbin[i];
     233     2195696 :     } else if( indices[i]>=nbin[i] ) {
     234           2 :       plumed_merror("point is outside grid range");
     235             :     }
     236             :   }
     237     1169231 : }
     238             : 
     239        2520 : unsigned GridCoordinatesObject::getIndex( const std::vector<double>& point ) const {
     240        2520 :   return getIndex( View<const double>(point.data(),point.size()) );
     241             : }
     242             : 
     243      594217 : unsigned GridCoordinatesObject::getIndex( const View<const double> point ) const {
     244             :   plumed_dbg_assert( bounds_set && point.size()==dimension );
     245      594217 :   if( gtype==flat ) {
     246      592117 :     std::vector<unsigned> indices(dimension);
     247      592117 :     getIndices( point, indices );
     248      592116 :     return getIndex( indices );
     249        2100 :   } else if( gtype==fibonacci ) {
     250        2100 :     return getFibonacciIndex( point );
     251             :   } else {
     252           0 :     plumed_error();
     253             :   }
     254             : }
     255             : 
     256      134280 : unsigned GridCoordinatesObject::getFibonacciIndex( const View<const double> p ) const {
     257             :   plumed_dbg_assert( gtype==fibonacci );
     258             :   // Convert input point to coordinates on cylinder
     259             :   int k=2;
     260      134280 :   double phi = atan2( p[2], p[0] ), sinthet2 = 1 - p[1]*p[1];
     261             :   // Calculate power to raise golden ratio
     262      134280 :   if( sinthet2<epsilon ) {
     263             :     k = 2;
     264             :   } else {
     265      134280 :     k = std::floor( std::log( npoints*pi*root5*sinthet2 ) / log_golden2 );
     266      134280 :     if( k<2 ) {
     267             :       k = 2;
     268             :     }
     269             :   }
     270      134280 :   double Fk = pow( golden, k ) / root5, F0 = std::round(Fk), F1 = std::round(Fk*golden);
     271      134280 :   Matrix<double> B(2,2), invB(2,2);
     272      134280 :   std::vector<double> thisp(3);
     273      134280 :   B(0,0) = 2*pi*((F0+1)*igolden - std::floor((F0+1)*igolden)) - fib_increment;
     274      134280 :   B(0,1) = 2*pi*((F1+1)*igolden - std::floor((F1+1)*igolden)) - fib_increment;
     275      134280 :   B(1,0) = -2*F0/npoints;
     276      134280 :   B(1,1) = -2*F1/npoints;
     277      134280 :   Invert( B, invB );
     278      134280 :   std::vector<double> vv(2), rc(2);
     279      134280 :   vv[0]=-phi;
     280      134280 :   vv[1] = p[1] - fib_shift;
     281      134280 :   mult( invB, vv, rc );
     282      134280 :   std::vector<int> c(2);
     283      134280 :   c[0]=std::floor(rc[0]);
     284      134280 :   c[1]=std::floor(rc[1]);
     285             :   unsigned outind=0;
     286             :   double mindist = 10000000.;
     287      671400 :   for(int s=0; s<4; ++s) {
     288      537120 :     double ttt, costheta = B(1,0)*( c[0] + s%2 ) + B(1,1)*( c[1] + s/2 ) + fib_shift;
     289      537120 :     if( costheta>1 ) {
     290             :       ttt=1;
     291      530862 :     } else if( costheta<-1 ) {
     292             :       ttt=-1;
     293             :     } else {
     294             :       ttt=costheta;
     295             :     }
     296      537120 :     costheta = 2*ttt - costheta;
     297      537120 :     unsigned i = std::floor( 0.5*npoints*(1+costheta) );
     298      537120 :     getFibonacciCoordinates( i, thisp );
     299             :     double dist=0;
     300     2148480 :     for(unsigned j=0; j<3; ++j) {
     301     1611360 :       double tmp=thisp[j]-p[j];
     302     1611360 :       dist += tmp*tmp;
     303             :     }
     304      537120 :     if( dist<mindist ) {
     305      294295 :       outind = i;
     306             :       mindist = dist;
     307             :     }
     308             :   }
     309      134280 :   return outind;
     310      134280 : }
     311             : 
     312    56848157 : void GridCoordinatesObject::convertIndexToIndices( const unsigned& index, const std::vector<unsigned>& nnbin, std::vector<unsigned>& indices ) const {
     313             :   plumed_dbg_assert( gtype==flat );
     314    56848157 :   unsigned kk=index;
     315    56848157 :   indices[0]=index%nnbin[0];
     316   107861389 :   for(unsigned i=1; i<dimension-1; ++i) {
     317    51013232 :     kk=(kk-indices[i-1])/nnbin[i-1];
     318    51013232 :     indices[i]=kk%nnbin[i];
     319             :   }
     320    56848157 :   if(dimension>=2) { // I think this is wrong
     321    55883264 :     indices[dimension-1]=(kk-indices[dimension-2])/nnbin[dimension-2];
     322             :   }
     323    56848157 : }
     324             : 
     325     5177699 : void GridCoordinatesObject::getIndices( const unsigned& index, std::vector<unsigned>& indices ) const {
     326             :   plumed_dbg_assert( gtype==flat );
     327     5177699 :   convertIndexToIndices( index, nbin, indices );
     328     5177699 : }
     329             : 
     330     1175350 : void GridCoordinatesObject::getGridPointCoordinates( const unsigned& ipoint, std::vector<double>& x ) const {
     331     1175350 :   std::vector<unsigned> tindices( dimension );
     332     1175350 :   getGridPointCoordinates( ipoint, tindices, x );
     333     1175350 : }
     334             : 
     335     1985967 : void GridCoordinatesObject::getGridPointCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const {
     336             :   plumed_dbg_assert( bounds_set && x.size()==dimension && tindices.size()==dimension && ipoint<npoints );
     337     1985967 :   if( gtype==flat ) {
     338     1975071 :     getFlatGridCoordinates( ipoint, tindices, x );
     339       10896 :   } else if( gtype==fibonacci ) {
     340       10896 :     getFibonacciCoordinates( ipoint, x );
     341             :   } else {
     342           0 :     plumed_error();
     343             :   }
     344     1985967 : }
     345             : 
     346           0 : void GridCoordinatesObject::putCoordinateAtValue( const unsigned& ind, const double& val, std::vector<double>& coords ) const {
     347           0 :   std::vector<double> point( dimension );
     348           0 :   getGridPointCoordinates( ind, point );
     349           0 :   if( gtype==flat ) {
     350           0 :     if( coords.size()!=(dimension+1) ) {
     351           0 :       coords.resize( (dimension+1) );
     352             :     }
     353           0 :     for(unsigned i=0; i<dimension; ++i) {
     354           0 :       coords[i]=point[i];
     355             :     }
     356           0 :     coords[point.size()]=val;
     357           0 :   } else if( gtype==fibonacci ) {
     358           0 :     if( coords.size()!=3 ) {
     359           0 :       coords.resize(3);
     360             :     }
     361           0 :     for(unsigned i=0; i<3; ++i) {
     362           0 :       coords[i] = val*point[i];
     363             :     }
     364             :   } else {
     365           0 :     plumed_error();
     366             :   }
     367           0 : }
     368             : 
     369     1975071 : void GridCoordinatesObject::getFlatGridCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const {
     370             :   plumed_dbg_assert( gtype==flat );
     371     1975071 :   getIndices( ipoint, tindices );
     372     6193793 :   for(unsigned i=0; i<dimension; ++i) {
     373     4218722 :     x[i] = min[i] + dx[i]*tindices[i];
     374             :   }
     375     1975071 : }
     376             : 
     377     5799992 : void GridCoordinatesObject::getFibonacciCoordinates( const unsigned& ipoint, std::vector<double>& x ) const {
     378             :   plumed_dbg_assert( gtype==fibonacci );
     379     5799992 :   x[1] = (ipoint*fib_offset) + fib_shift;
     380     5799992 :   double r = sqrt( 1 - x[1]*x[1] );
     381     5799992 :   double phi = ipoint*fib_increment;
     382     5799992 :   x[0] = r*cos(phi);
     383     5799992 :   x[2] = r*sin(phi);
     384             :   double norm=0;
     385    23199968 :   for(unsigned j=0; j<3; ++j) {
     386    17399976 :     norm+=x[j]*x[j];
     387             :   }
     388     5799992 :   norm = sqrt(norm);
     389    23199968 :   for(unsigned j=0; j<3; ++j) {
     390    17399976 :     x[j] = x[j] / norm;
     391             :   }
     392     5799992 : }
     393             : 
     394      552562 : void GridCoordinatesObject::getSplineNeighbors( const unsigned& mybox, unsigned& nneighbors, std::vector<unsigned>& mysneigh ) const {
     395             :   plumed_dbg_assert( gtype==flat );
     396      552562 :   mysneigh.resize( static_cast<unsigned>(pow(2.,dimension)) );
     397             : 
     398             :   unsigned inind;
     399      552562 :   nneighbors = 0;
     400      552562 :   std::vector<unsigned> tmp_indices( dimension );
     401      552562 :   std::vector<unsigned> my_indices( dimension );
     402      552562 :   getIndices( mybox, my_indices );
     403     2779388 :   for(unsigned i=0; i<mysneigh.size(); ++i) {
     404             :     unsigned tmp=i;
     405             :     inind=0;
     406     6750996 :     for(unsigned j=0; j<dimension; ++j) {
     407     4524170 :       unsigned i0=tmp%2+my_indices[j];
     408     4524170 :       tmp/=2;
     409     4524170 :       if(!pbc[j] && i0==nbin[j]) {
     410       40800 :         continue;
     411             :       }
     412     4483370 :       if( pbc[j] && i0==nbin[j]) {
     413             :         i0=0;
     414             :       }
     415     4483370 :       tmp_indices[inind++]=i0;
     416             :     }
     417     2226826 :     if( inind==dimension ) {
     418     2186226 :       mysneigh[nneighbors++]=getIndex( tmp_indices );
     419             :     }
     420             :   }
     421      552562 : }
     422             : 
     423      561694 : std::vector<std::string> GridCoordinatesObject::getMin() const {
     424             :   plumed_dbg_assert( gtype==flat );
     425      561694 :   return str_min;
     426             : }
     427             : 
     428      561694 : std::vector<std::string> GridCoordinatesObject::getMax() const {
     429             :   plumed_dbg_assert( gtype==flat );
     430      561694 :   return str_max;
     431             : }
     432             : 
     433      573146 : std::vector<std::size_t> GridCoordinatesObject::getNbin( const bool& shape ) const {
     434             :   plumed_dbg_assert( gtype==flat && nbin.size()==dimension );
     435      573146 :   std::vector<std::size_t> ngrid( dimension );
     436     1951806 :   for(unsigned i=0; i<dimension; ++i) {
     437     1378660 :     if( !pbc[i] && !shape ) {
     438      665970 :       ngrid[i]=nbin[i] - 1;
     439             :     } else {
     440      712690 :       ngrid[i]=nbin[i];
     441             :     }
     442             :   }
     443      573146 :   return ngrid;
     444             : }
     445             : 
     446           0 : void GridCoordinatesObject::getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh,
     447             :     unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
     448           0 :   getNeighbors( View( pp.data(), pp.size() ), nneigh, num_neighbors, neighbors );
     449           0 : }
     450             : 
     451      343692 : void GridCoordinatesObject::getNeighbors( const View<const double> pp, const std::vector<unsigned>& nneigh,
     452             :     unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
     453             :   plumed_dbg_assert( bounds_set );
     454             : 
     455      343692 :   if( gtype == flat ) {
     456             :     plumed_dbg_assert( nneigh.size()==dimension );
     457      211512 :     std::vector<unsigned> indices( dimension );
     458      490397 :     for(unsigned i=0; i<dimension; ++i) {
     459      278885 :       indices[i] = std::floor( (pp[i]-min[i])/dx[i] );
     460             :     }
     461      211512 :     getNeighbors( indices, nneigh, num_neighbors, neighbors );
     462      132180 :   } else if( gtype == fibonacci ) {
     463      132180 :     unsigned find = getFibonacciIndex( View<const double>(pp.data(), pp.size()) );
     464      132180 :     num_neighbors = 1 + fib_nlist[find].size();
     465      132180 :     if( neighbors.size()<num_neighbors ) {
     466      132180 :       neighbors.resize( num_neighbors );
     467             :     }
     468      132180 :     neighbors[0]=find;
     469    18175217 :     for(unsigned i=0; i<fib_nlist[find].size(); ++i) {
     470    18043037 :       neighbors[1+i] = fib_nlist[find][i];
     471             :     }
     472             :   } else {
     473           0 :     plumed_error();
     474             :   }
     475      343692 : }
     476             : 
     477      227976 : void GridCoordinatesObject::getNeighbors( const std::vector<unsigned>& indices, const std::vector<unsigned>& nneigh,
     478             :     unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
     479             :   plumed_dbg_assert( gtype==flat && bounds_set && nneigh.size()==dimension );
     480             : 
     481             :   unsigned num_neigh=1;
     482      227976 :   std::vector<unsigned> small_bin( dimension );
     483      556253 :   for(unsigned i=0; i<dimension; ++i) {
     484      328277 :     small_bin[i]=(2*nneigh[i]+1);
     485      328277 :     if( pbc[i] && small_bin[i]>nbin[i] ) {
     486      104400 :       small_bin[i]=nbin[i];
     487             :     }
     488      328277 :     num_neigh *=small_bin[i];
     489             :   }
     490      227976 :   if( neighbors.size()!=num_neigh ) {
     491      211515 :     neighbors.resize( num_neigh );
     492             :   }
     493             : 
     494      227976 :   num_neighbors=0;
     495      227976 :   std::vector<unsigned> s_indices(dimension);
     496      227976 :   std::vector<unsigned> t_indices(dimension);
     497    51898434 :   for(unsigned index=0; index<num_neigh; ++index) {
     498             :     bool found=true;
     499    51670458 :     convertIndexToIndices( index, small_bin, s_indices );
     500   204692505 :     for(unsigned i=0; i<dimension; ++i) {
     501   153022047 :       int i0=s_indices[i]-nneigh[i]+indices[i];
     502   153022047 :       if (i0<0) {
     503    14398084 :         if (!pbc[i]) {
     504             :           found = false;
     505             :         } else {
     506     8410766 :           i0=nbin[i]-(-i0)%nbin[i];
     507             :         }
     508   138623963 :       } else if(static_cast<unsigned>(i0)>=nbin[i]) { // i0 is >=0
     509    15504616 :         if (!pbc[i]) {
     510             :           found = false;
     511             :         } else {
     512     8653206 :           i0%=nbin[i];
     513             :         }
     514             :       }
     515   153022047 :       if(!pbc[i]) {
     516    35502337 :         if( i0<0 || static_cast<unsigned>(i0)>=nbin[i] ) { // i0 is >=0
     517             :           found=false;
     518             :         }
     519             :       } else {
     520   117519710 :         if(i0<0) {
     521           0 :           i0=nbin[i]-(-i0)%nbin[i];
     522             :         } else {
     523   117519710 :           i0%=nbin[i];
     524             :         }
     525             :       }
     526   153022047 :       t_indices[i]=static_cast<unsigned>(i0);
     527             :     }
     528    51670458 :     if( found ) {
     529    38975870 :       neighbors[num_neighbors]=getIndex( t_indices );
     530    38975870 :       num_neighbors++;
     531             :     }
     532             :   }
     533      227976 : }
     534             : 
     535             : }
     536             : }
     537             : 

Generated by: LCOV version 1.16