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         335 : void GridCoordinatesObject::setup( const std::string& geom, const std::vector<bool>& ipbc,
      30             :                                    const unsigned& np, const double& fib_cutoff ) {
      31         335 :   if( geom=="flat" ) { gtype=flat; dimension = ipbc.size(); }
      32          17 :   else if( geom=="fibonacci" ) { gtype=fibonacci; dimension = 3; }
      33           0 :   else plumed_merror( geom + " is invalid geometry type");
      34             : 
      35         335 :   if( gtype==flat ) {
      36         678 :     bounds_set=false; npoints=0; pbc.resize( ipbc.size() ); for(unsigned i=0; i<ipbc.size(); ++i) pbc[i]=ipbc[i];
      37          17 :   } else if( gtype==fibonacci ) {
      38          17 :     bounds_set=true; root5 = sqrt(5);
      39          17 :     npoints = np; golden = ( 1 + sqrt(5) ) / 2.0; igolden = golden - 1;
      40          17 :     fib_increment = 2*pi*igolden; log_golden2 = std::log( golden*golden );
      41          17 :     fib_offset = 2 / static_cast<double>( npoints );
      42          17 :     fib_shift = fib_offset/2 - 1;
      43             : 
      44          17 :     std::vector<double> icoord( dimension ), jcoord( dimension );
      45             :     // Find minimum distance between each pair of points
      46          17 :     std::vector<unsigned> tindices( dimension ); std::vector<double> mindists( npoints );
      47        4849 :     for(unsigned i=0; i<npoints; ++i) {
      48        4832 :       getFibonacciCoordinates( i, icoord ); mindists[i] = 0;
      49     1707040 :       for(unsigned j=0; j<npoints; ++j) {
      50     1702208 :         if( i==j ) continue ; // Points are not neighbors to themselves
      51     1697376 :         getFibonacciCoordinates( j, jcoord );
      52             :         // Calculate the dot product
      53     6789504 :         double dot=0; for(unsigned k=0; k<dimension; ++k) dot += icoord[k]*jcoord[k];
      54     1697376 :         if( dot>mindists[i] ) mindists[i]=dot;
      55             :       }
      56             :     }
      57             :     // And now take minimum of dot products
      58          17 :     double min=mindists[0];
      59        4832 :     for(unsigned i=1; i<npoints; ++i) {
      60        4815 :       if( mindists[i]<min ) min=mindists[i];
      61             :     }
      62             :     double final_cutoff;
      63          17 :     if( fib_cutoff<-1 ) final_cutoff=-1;
      64          15 :     else final_cutoff = cos( acos( fib_cutoff ) + acos( min ) );
      65             : 
      66             :     // And now construct the neighbor list
      67          17 :     fib_nlist.resize( npoints );
      68        4849 :     for(unsigned i=0; i<npoints; ++i) {
      69        4832 :       getFibonacciCoordinates( i, icoord );
      70     1707040 :       for(unsigned j=0; j<npoints; ++j) {
      71     1702208 :         if( i==j ) continue ; // Points are not neighbors to themselves
      72     1697376 :         getFibonacciCoordinates( j, jcoord );
      73             :         // Calculate the dot product
      74     6789504 :         double dot=0; for(unsigned k=0; k<dimension; ++k) dot += icoord[k]*jcoord[k];
      75     1697376 :         if( dot>final_cutoff ) { fib_nlist[i].push_back(j); }
      76             :       }
      77             :     }
      78           0 :   } else plumed_error();
      79         335 : }
      80             : 
      81         618 : void GridCoordinatesObject::setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax,
      82             :                                        const std::vector<unsigned>& binsin, std::vector<double>& spacing ) {
      83             :   plumed_dbg_assert( smin.size()==dimension && smax.size()==dimension );
      84         618 :   plumed_assert( gtype==flat && (spacing.size()==dimension || binsin.size()==dimension) );
      85         618 :   str_min.resize( dimension ); str_max.resize( dimension ); nbin.resize( dimension );
      86         618 :   min.resize( dimension ); max.resize( dimension ); dx.resize( dimension ); stride.resize( dimension );
      87             : 
      88        1164 :   npoints=1; bounds_set=(smin[0]!="auto" && smax[0]!="auto");
      89         618 :   if( bounds_set ) {
      90         596 :     for(unsigned i=1; i<dimension; ++i) {
      91         100 :       if( smin[i]=="auto" || smax[i]=="auto" ) { bounds_set=false; break; }
      92             :     }
      93             :   }
      94        1302 :   for(unsigned i=0; i<dimension; ++i) {
      95         684 :     str_min[i]=smin[i]; str_max[i]=smax[i];
      96         684 :     if( bounds_set ) {
      97         596 :       Tools::convert( str_min[i], min[i] );
      98         596 :       Tools::convert( str_max[i], max[i] );
      99             :     }
     100         684 :     if( spacing.size()==dimension && binsin.size()==dimension ) {
     101         290 :       if( spacing[i]==0 ) nbin[i] = binsin[i];
     102         273 :       else if( bounds_set ) {
     103         273 :         double range = max[i] - min[i]; nbin[i] = std::round( range / spacing[i]); dx[i]=spacing[i];
     104             :         // This check ensures that nbins is set correctly if spacing is set the same as the number of bins
     105         273 :         if( nbin[i]!=binsin[i] ) plumed_merror("mismatch between input spacing and input number of bins");
     106             :       }
     107         394 :     } else if( binsin.size()==dimension ) {
     108         394 :       nbin[i]=binsin[i]; dx[i] = ( max[i] - min[i] ) / static_cast<double>( nbin[i] );
     109           0 :     } else if( spacing.size()==dimension && bounds_set ) {
     110           0 :       nbin[i] = std::floor(( max[i] - min[i] ) / spacing[i]) + 1; dx[i]=spacing[i];
     111           0 :     } else if( bounds_set ) plumed_error();
     112         684 :     if( !pbc[i] ) { max[i] +=dx[i]; nbin[i]+=1; }
     113         684 :     stride[i]=npoints;
     114         684 :     npoints*=nbin[i];
     115             :   }
     116         618 :   if( spacing.size()!=dimension && bounds_set ) {
     117         616 :     spacing.resize(dimension); for(unsigned i=0; i<dimension; ++i) spacing[i]=dx[i];
     118             :   }
     119         618 : }
     120             : 
     121    42783311 : unsigned GridCoordinatesObject::getIndex( const std::vector<unsigned>& indices ) const {
     122             :   plumed_dbg_assert( gtype==flat && indices.size()==dimension );
     123             :   // indices are flattended using a column-major order
     124    42783311 :   unsigned index=indices[dimension-1];
     125   122880687 :   for(unsigned i=dimension-1; i>0; --i) {
     126    80097376 :     index=index*nbin[i-1]+indices[i-1];
     127             :   }
     128    42783311 :   return index;
     129             : }
     130             : 
     131      197426 : bool GridCoordinatesObject::inbounds( const std::vector<double>& point ) const {
     132      197426 :   if( gtype==fibonacci ) return true;
     133             :   plumed_dbg_assert( bounds_set && point.size()==dimension );
     134      461467 :   for(unsigned i=0; i<dimension; ++i) {
     135      296963 :     if( pbc[i] ) continue;
     136       60839 :     if( point[i]<min[i] || point[i]>(max[i]-dx[i]) ) return false;
     137             :   }
     138             :   return true;
     139             : }
     140             : 
     141     1175406 : void GridCoordinatesObject::getIndices( const std::vector<double>& point, std::vector<unsigned>& indices ) const {
     142             :   plumed_dbg_assert( gtype==flat && bounds_set && point.size()==dimension && indices.size()==dimension );
     143     3445432 :   for(unsigned i=0; i<dimension; ++i) {
     144     2270026 :     indices[i]=std::floor( (point[i] - min[i])/dx[i] );
     145     2270026 :     if( pbc[i] ) indices[i]=indices[i]%nbin[i];
     146     2179718 :     else if( indices[i]>nbin[i] ) plumed_merror("point is outside grid range");
     147             :   }
     148     1175406 : }
     149             : 
     150      615967 : unsigned GridCoordinatesObject::getIndex( const std::vector<double>& point ) const {
     151             :   plumed_dbg_assert( bounds_set && point.size()==dimension );
     152      615967 :   if( gtype==flat ) {
     153      613867 :     std::vector<unsigned> indices(dimension); getIndices( point, indices );
     154      613867 :     return getIndex( indices );
     155        2100 :   } else if( gtype==fibonacci ) {
     156        2100 :     return getFibonacciIndex( point );
     157             :   } else {
     158           0 :     plumed_error();
     159             :   }
     160             : }
     161             : 
     162       16510 : unsigned GridCoordinatesObject::getFibonacciIndex( const std::vector<double>& p ) const {
     163             :   plumed_dbg_assert( gtype==fibonacci );
     164             :   // Convert input point to coordinates on cylinder
     165       16510 :   int k=2; double phi = atan2( p[2], p[0] ), sinthet2 = 1 - p[1]*p[1];
     166             :   // Calculate power to raise golden ratio
     167       16510 :   if( sinthet2<epsilon ) { k = 2; }
     168             :   else {
     169       16510 :     k = std::floor( std::log( npoints*pi*root5*sinthet2 ) / log_golden2 );
     170             :     if( k<2 ) k = 2;
     171             :   }
     172       16510 :   double Fk = pow( golden, k ) / root5, F0 = std::round(Fk), F1 = std::round(Fk*golden);
     173       16510 :   Matrix<double> B(2,2), invB(2,2); std::vector<double> thisp(3);
     174       16510 :   B(0,0) = 2*pi*((F0+1)*igolden - std::floor((F0+1)*igolden)) - fib_increment;
     175       16510 :   B(0,1) = 2*pi*((F1+1)*igolden - std::floor((F1+1)*igolden)) - fib_increment;
     176       16510 :   B(1,0) = -2*F0/npoints; B(1,1) = -2*F1/npoints; Invert( B, invB );
     177       16510 :   std::vector<double> vv(2), rc(2); vv[0]=-phi; vv[1] = p[1] - fib_shift;
     178       16510 :   mult( invB, vv, rc ); std::vector<int> c(2); c[0]=std::floor(rc[0]); c[1]=std::floor(rc[1]);
     179             :   unsigned outind=0; double mindist = 10000000.;
     180       82550 :   for(int s=0; s<4; ++s) {
     181       66040 :     double ttt, costheta = B(1,0)*( c[0] + s%2 ) + B(1,1)*( c[1] + s/2 ) + fib_shift;
     182       66040 :     if( costheta>1 ) ttt=1; else if( costheta<-1 ) ttt=-1; else ttt=costheta;
     183       66040 :     costheta = 2*ttt - costheta;
     184       66040 :     unsigned i = std::floor( 0.5*npoints*(1+costheta) ); getFibonacciCoordinates( i, thisp );
     185      264160 :     double dist=0; for(unsigned j=0; j<3; ++j) { double tmp=thisp[j]-p[j]; dist += tmp*tmp; }
     186       66040 :     if( dist<mindist ) { outind = i; mindist = dist; }
     187             :   }
     188       16510 :   return outind;
     189             : }
     190             : 
     191    94917677 : void GridCoordinatesObject::convertIndexToIndices( const unsigned& index, const std::vector<unsigned>& nnbin, std::vector<unsigned>& indices ) const {
     192    94917677 :   plumed_dbg_assert( gtype==flat ); unsigned kk=index; indices[0]=index%nnbin[0];
     193   183492225 :   for(unsigned i=1; i<dimension-1; ++i) {
     194    88574548 :     kk=(kk-indices[i-1])/nnbin[i-1];
     195    88574548 :     indices[i]=kk%nnbin[i];
     196             :   }
     197    94917677 :   if(dimension>=2) { // I think this is wrong
     198    92728227 :     indices[dimension-1]=(kk-indices[dimension-2])/nnbin[dimension-2];
     199             :   }
     200    94917677 : }
     201             : 
     202    42760057 : void GridCoordinatesObject::getIndices( const unsigned& index, std::vector<unsigned>& indices ) const {
     203    42760057 :   plumed_dbg_assert( gtype==flat ); convertIndexToIndices( index, nbin, indices );
     204    42760057 : }
     205             : 
     206    40720264 : void GridCoordinatesObject::getGridPointCoordinates( const unsigned& ipoint, std::vector<double>& x ) const {
     207    40720264 :   std::vector<unsigned> tindices( dimension ); getGridPointCoordinates( ipoint, tindices, x );
     208    40720264 : }
     209             : 
     210    41484313 : void GridCoordinatesObject::getGridPointCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const {
     211             :   plumed_dbg_assert( bounds_set && x.size()==dimension && tindices.size()==dimension && ipoint<npoints );
     212    41484313 :   if( gtype==flat ) {
     213    39980732 :     getFlatGridCoordinates( ipoint, tindices, x );
     214     1503581 :   } else if( gtype==fibonacci ) {
     215     1503581 :     getFibonacciCoordinates( ipoint, x );
     216             :   } else {
     217           0 :     plumed_error();
     218             :   }
     219    41484313 : }
     220             : 
     221           0 : void GridCoordinatesObject::putCoordinateAtValue( const unsigned& ind, const double& val, std::vector<double>& coords ) const {
     222           0 :   std::vector<double> point( dimension ); getGridPointCoordinates( ind, point );
     223           0 :   if( gtype==flat ) {
     224           0 :     if( coords.size()!=(dimension+1) ) coords.resize( (dimension+1) );
     225           0 :     for(unsigned i=0; i<dimension; ++i) coords[i]=point[i]; coords[point.size()]=val;
     226           0 :   } else if( gtype==fibonacci ) {
     227           0 :     if( coords.size()!=3 ) coords.resize(3);
     228           0 :     for(unsigned i=0; i<3; ++i) coords[i] = val*point[i];
     229             :   } else {
     230           0 :     plumed_error();
     231             :   }
     232           0 : }
     233             : 
     234    39980732 : void GridCoordinatesObject::getFlatGridCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const {
     235    39980732 :   plumed_dbg_assert( gtype==flat ); getIndices( ipoint, tindices );
     236   156503342 :   for(unsigned i=0; i<dimension; ++i) x[i] = min[i] + dx[i]*tindices[i];
     237    39980732 : }
     238             : 
     239     4974037 : void GridCoordinatesObject::getFibonacciCoordinates( const unsigned& ipoint, std::vector<double>& x ) const {
     240             :   plumed_dbg_assert( gtype==fibonacci );
     241     4974037 :   x[1] = (ipoint*fib_offset) + fib_shift; double r = sqrt( 1 - x[1]*x[1] );
     242     4974037 :   double phi = ipoint*fib_increment; x[0] = r*cos(phi); x[2] = r*sin(phi);
     243    19896148 :   double norm=0; for(unsigned j=0; j<3; ++j) norm+=x[j]*x[j];
     244    19896148 :   norm = sqrt(norm); for(unsigned j=0; j<3; ++j) x[j] = x[j] / norm;
     245     4974037 : }
     246             : 
     247      544167 : void GridCoordinatesObject::getSplineNeighbors( const unsigned& mybox, unsigned& nneighbors, std::vector<unsigned>& mysneigh ) const {
     248      544167 :   plumed_dbg_assert( gtype==flat ); mysneigh.resize( static_cast<unsigned>(pow(2.,dimension)) );
     249             : 
     250      544167 :   unsigned inind; nneighbors = 0;
     251      544167 :   std::vector<unsigned> tmp_indices( dimension );
     252      544167 :   std::vector<unsigned> my_indices( dimension );
     253      544167 :   getIndices( mybox, my_indices );
     254     2754079 :   for(unsigned i=0; i<mysneigh.size(); ++i) {
     255             :     unsigned tmp=i; inind=0;
     256     6716896 :     for(unsigned j=0; j<dimension; ++j) {
     257     4506984 :       unsigned i0=tmp%2+my_indices[j]; tmp/=2;
     258     4506984 :       if(!pbc[j] && i0==nbin[j]) continue;
     259     4466184 :       if( pbc[j] && i0==nbin[j]) i0=0;
     260     4466184 :       tmp_indices[inind++]=i0;
     261             :     }
     262     2209912 :     if( inind==dimension ) mysneigh[nneighbors++]=getIndex( tmp_indices );
     263             :   }
     264      544167 : }
     265             : 
     266      470843 : std::vector<std::string> GridCoordinatesObject::getMin() const {
     267      470843 :   plumed_dbg_assert( gtype==flat ); return str_min;
     268             : }
     269             : 
     270      470843 : std::vector<std::string> GridCoordinatesObject::getMax() const {
     271      470843 :   plumed_dbg_assert( gtype==flat ); return str_max;
     272             : }
     273             : 
     274      483329 : std::vector<unsigned> GridCoordinatesObject::getNbin( const bool& shape ) const {
     275             :   plumed_dbg_assert( gtype==flat && nbin.size()==dimension );
     276      483329 :   std::vector<unsigned> ngrid( dimension );
     277     1591468 :   for(unsigned i=0; i<dimension; ++i) {
     278     1108139 :     if( !pbc[i] && !shape ) ngrid[i]=nbin[i] - 1;
     279      534057 :     else ngrid[i]=nbin[i];
     280             :   }
     281      483329 :   return ngrid;
     282             : }
     283             : 
     284      231602 : void GridCoordinatesObject::getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh,
     285             :     unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
     286             :   plumed_dbg_assert( bounds_set );
     287             : 
     288      231602 :   if( gtype == flat ) {
     289             :     plumed_dbg_assert( nneigh.size()==dimension );
     290      217192 :     std::vector<unsigned> indices( dimension );
     291      490941 :     for(unsigned i=0; i<dimension; ++i) indices[i] = std::floor( (pp[i]-min[i])/dx[i] );
     292      217192 :     getNeighbors( indices, nneigh, num_neighbors, neighbors );
     293       14410 :   } else if( gtype == fibonacci ) {
     294       14410 :     unsigned find = getFibonacciIndex( pp );
     295       14410 :     num_neighbors = 1 + fib_nlist[find].size();
     296       14410 :     if( neighbors.size()<num_neighbors ) neighbors.resize( num_neighbors );
     297     1500997 :     neighbors[0]=find; for(unsigned i=0; i<fib_nlist[find].size(); ++i) neighbors[1+i] = fib_nlist[find][i];
     298             :   } else {
     299           0 :     plumed_error();
     300             :   }
     301      231602 : }
     302             : 
     303      245010 : void GridCoordinatesObject::getNeighbors( const std::vector<unsigned>& indices, const std::vector<unsigned>& nneigh,
     304             :     unsigned& num_neighbors, std::vector<unsigned>& neighbors ) const {
     305             :   plumed_dbg_assert( gtype==flat && bounds_set && nneigh.size()==dimension );
     306             : 
     307      245010 :   unsigned num_neigh=1; std::vector<unsigned> small_bin( dimension );
     308      602213 :   for(unsigned i=0; i<dimension; ++i) {
     309      357203 :     small_bin[i]=(2*nneigh[i]+1);
     310      357203 :     if( pbc[i] && small_bin[i]>nbin[i] ) small_bin[i]=nbin[i];
     311      357203 :     num_neigh *=small_bin[i];
     312             :   }
     313      245010 :   if( neighbors.size()!=num_neigh ) neighbors.resize( num_neigh );
     314             : 
     315      245010 :   num_neighbors=0;
     316      245010 :   std::vector<unsigned> s_indices(dimension), t_indices(dimension);
     317    52402630 :   for(unsigned index=0; index<num_neigh; ++index) {
     318             :     bool found=true;
     319    52157620 :     convertIndexToIndices( index, small_bin, s_indices );
     320   206161849 :     for(unsigned i=0; i<dimension; ++i) {
     321   154004229 :       int i0=s_indices[i]-nneigh[i]+indices[i];
     322   154004229 :       if(!pbc[i] && i0<0)        found=false;
     323   154004229 :       if(!pbc[i] && i0>=nbin[i]) found=false;
     324   154004229 :       if( pbc[i] && i0<0)        i0=nbin[i]-(-i0)%nbin[i];
     325   154004229 :       if( pbc[i] && i0>=nbin[i]) i0%=nbin[i];
     326   154004229 :       t_indices[i]=static_cast<unsigned>(i0);
     327             :     }
     328    52157620 :     if( found ) {
     329    39399065 :       neighbors[num_neighbors]=getIndex( t_indices );
     330    39399065 :       num_neighbors++;
     331             :     }
     332             :   }
     333      245010 : }
     334             : 
     335             : }
     336             : }
     337             :