Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2013-2017 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 "core/ActionWithMatrix.h"
      23             : #include "core/ActionRegister.h"
      24             : 
      25             : //+PLUMEDOC MCOLVAR NEIGHBORS
      26             : /*
      27             : Build a matrix with ones in for the N nearest neighbours of an atom
      28             : 
      29             : \par Examples
      30             : 
      31             : */
      32             : //+ENDPLUMEDOC
      33             : 
      34             : namespace PLMD {
      35             : namespace adjmat {
      36             : 
      37             : class Neighbors : public ActionWithMatrix {
      38             :   bool lowest;
      39             :   unsigned number;
      40             : public:
      41             :   static void registerKeywords( Keywords& keys );
      42             :   explicit Neighbors(const ActionOptions&);
      43             :   unsigned getNumberOfDerivatives() override;
      44           8 :   unsigned getNumberOfColumns() const override { return number; }
      45           0 :   bool canBeAfterInChain( ActionWithVector* av ) { return av->getLabel()!=(getPntrToArgument(0)->getPntrToAction())->getLabel(); }
      46             :   void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const ;
      47             :   void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override;
      48        2512 :   void runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const override {}
      49             :   void turnOnDerivatives() override ;
      50             : };
      51             : 
      52             : PLUMED_REGISTER_ACTION(Neighbors,"NEIGHBORS")
      53             : 
      54           7 : void Neighbors::registerKeywords( Keywords& keys ) {
      55           7 :   ActionWithMatrix::registerKeywords( keys );
      56          14 :   keys.addInputKeyword("compulsory","ARG","matrix","the label of an adjacency/distance matrix that will be used to find the nearest neighbors");
      57          14 :   keys.add("compulsory","NLOWEST","0","in each row of the output matrix set the elements that correspond to the n lowest elements in each row of the input matrix equal to one");
      58          14 :   keys.add("compulsory","NHIGHEST","0","in each row of the output matrix set the elements that correspond to the n highest elements in each row of the input matrix equal to one");
      59          14 :   keys.setValueDescription("matrix","a matrix in which the ij element is one if the ij-element of the input matrix is one of the NLOWEST/NHIGHEST elements on that row of the input matrix and zero otherwise");
      60           7 : }
      61             : 
      62           3 : Neighbors::Neighbors(const ActionOptions&ao):
      63             :   Action(ao),
      64           3 :   ActionWithMatrix(ao)
      65             : {
      66           3 :   if( getNumberOfArguments()!=1 ) error("found wrong number of arguments in input");
      67           3 :   if( getPntrToArgument(0)->getRank()!=2 ) error("input argument should be a matrix");
      68           3 :   getPntrToArgument(0)->buildDataStore();
      69             : 
      70           3 :   unsigned nlow; parse("NLOWEST",nlow);
      71           3 :   unsigned nhigh; parse("NHIGHEST",nhigh);
      72           3 :   if( nlow==0 && nhigh==0 ) error("missing NLOWEST or NHIGHEST keyword one of these two keywords must be set in input");
      73           3 :   if( nlow>0 && nhigh>0 ) error("should only be one of NLOWEST or NHIGHEST set in input");
      74           3 :   if( nlow>0 ) {
      75           3 :     number=nlow; lowest=true;
      76           3 :     log.printf("  output matrix will have non-zero values for elements that correpsond to the %d lowest elements in each row of the input matrix\n",number);
      77             :   }
      78           3 :   if( nhigh>0 ) {
      79           0 :     number=nhigh; lowest=false;
      80           0 :     log.printf("  output matrix will have non-zero values for elements that correpsond to the %d highest elements in each row of the input matrix\n",number);
      81             :   }
      82             : 
      83             :   // And get the shape
      84           3 :   std::vector<unsigned> shape( getPntrToArgument(0)->getShape() );
      85           3 :   addValue( shape ); setNotPeriodic();
      86           3 :   checkRead();
      87           3 : }
      88             : 
      89           0 : void Neighbors::turnOnDerivatives() {
      90           0 :   ActionWithValue::turnOnDerivatives();
      91           0 :   warning("think about whether your symmetry functions are continuous. If the symmetry function can be calculated from distances only then you can use NEIGHBORS. If you calculate angles between vectors or use the vectors directly then the symmetry function computed using NEIGHBORS is not continuous.  It does not make sense to use such CVs when biasing");
      92           0 : }
      93             : 
      94           0 : unsigned Neighbors::getNumberOfDerivatives() {
      95           0 :   return 0;
      96             : }
      97             : 
      98        2512 : void Neighbors::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
      99        2512 :   const Value* wval = getPntrToArgument(0); unsigned nbonds = wval->getRowLength( task_index ), ncols = wval->getShape()[1];
     100        2512 :   if( indices.size()!=1+number ) indices.resize( 1 + number );
     101        2512 :   myvals.setSplitIndex(1+number);
     102             : 
     103             :   unsigned nind=0;
     104      180890 :   for(unsigned i=0; i<nbonds; ++i) {
     105      178378 :     unsigned ipos = ncols*task_index + wval->getRowIndex( task_index, i );
     106      178378 :     double weighti = wval->get( ipos );
     107      178378 :     if( weighti<epsilon ) continue ;
     108      177866 :     nind++;
     109             :   }
     110        2512 :   if( number>nind ) plumed_merror("not enough matrix elements were stored");
     111             : 
     112             :   // Now build vectors for doing sorting
     113        2512 :   std::vector<std::pair<double,unsigned> > rows( nind ); unsigned n=0;
     114      180890 :   for(unsigned i=0; i<nbonds; ++i) {
     115             :     unsigned iind = wval->getRowIndex( task_index, i );
     116      178378 :     unsigned ipos = ncols*task_index + iind;
     117      178378 :     double weighti = wval->get( ipos );
     118      178378 :     if( weighti<epsilon ) continue ;
     119      177866 :     rows[n].first=weighti; rows[n].second=iind; n++;
     120             :   }
     121             : 
     122             :   // Now do the sort and clear all the stored values ready for recompute
     123        2512 :   std::sort( rows.begin(), rows.end() );
     124             :   // This is to make this action consistent with what in other matrix actions
     125        2512 :   unsigned start_n = getPntrToArgument(0)->getShape()[0];
     126             :   // And setup the lowest indices, which are the ones we need to calculate
     127       16560 :   for(unsigned i=0; i<number; ++i) {
     128       14048 :     indices[i+1] = start_n + rows[nind-1-i].second;
     129       14048 :     if( lowest ) indices[i+1] = start_n + rows[i].second;
     130             :   }
     131        2512 : }
     132             : 
     133       14048 : void Neighbors::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
     134       14048 :   myvals.addValue( getConstPntrToComponent(0)->getPositionInStream(), 1.0 );
     135       14048 : }
     136             : 
     137             : }
     138             : }