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 ); keys.use("ARG");
      56          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");
      57          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");
      58           7 :   keys.setValueDescription("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");
      59           7 : }
      60             : 
      61           3 : Neighbors::Neighbors(const ActionOptions&ao):
      62             :   Action(ao),
      63           3 :   ActionWithMatrix(ao)
      64             : {
      65           3 :   if( getNumberOfArguments()!=1 ) error("found wrong number of arguments in input");
      66           3 :   if( getPntrToArgument(0)->getRank()!=2 ) error("input argument should be a matrix");
      67           3 :   getPntrToArgument(0)->buildDataStore();
      68             : 
      69           3 :   unsigned nlow; parse("NLOWEST",nlow);
      70           3 :   unsigned nhigh; parse("NHIGHEST",nhigh);
      71           3 :   if( nlow==0 && nhigh==0 ) error("missing NLOWEST or NHIGHEST keyword one of these two keywords must be set in input");
      72           3 :   if( nlow>0 && nhigh>0 ) error("should only be one of NLOWEST or NHIGHEST set in input");
      73           3 :   if( nlow>0 ) {
      74           3 :     number=nlow; lowest=true;
      75           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);
      76             :   }
      77           3 :   if( nhigh>0 ) {
      78           0 :     number=nhigh; lowest=false;
      79           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);
      80             :   }
      81             : 
      82             :   // And get the shape
      83           3 :   std::vector<unsigned> shape( getPntrToArgument(0)->getShape() );
      84           3 :   addValue( shape ); setNotPeriodic();
      85           3 :   checkRead();
      86           3 : }
      87             : 
      88           0 : void Neighbors::turnOnDerivatives() {
      89           0 :   ActionWithValue::turnOnDerivatives();
      90           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");
      91           0 : }
      92             : 
      93           0 : unsigned Neighbors::getNumberOfDerivatives() {
      94           0 :   return 0;
      95             : }
      96             : 
      97        2512 : void Neighbors::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
      98        2512 :   const Value* wval = getPntrToArgument(0); unsigned nbonds = wval->getRowLength( task_index ), ncols = wval->getShape()[1];
      99        2512 :   if( indices.size()!=1+number ) indices.resize( 1 + number );
     100        2512 :   myvals.setSplitIndex(1+number);
     101             : 
     102             :   unsigned nind=0;
     103      180890 :   for(unsigned i=0; i<nbonds; ++i) {
     104      178378 :     unsigned ipos = ncols*task_index + wval->getRowIndex( task_index, i );
     105      178378 :     double weighti = wval->get( ipos );
     106      178378 :     if( weighti<epsilon ) continue ;
     107      177866 :     nind++;
     108             :   }
     109        2512 :   if( number>nind ) plumed_merror("not enough matrix elements were stored");
     110             : 
     111             :   // Now build vectors for doing sorting
     112        2512 :   std::vector<std::pair<double,unsigned> > rows( nind ); unsigned n=0;
     113      180890 :   for(unsigned i=0; i<nbonds; ++i) {
     114             :     unsigned iind = wval->getRowIndex( task_index, i );
     115      178378 :     unsigned ipos = ncols*task_index + iind;
     116      178378 :     double weighti = wval->get( ipos );
     117      178378 :     if( weighti<epsilon ) continue ;
     118      177866 :     rows[n].first=weighti; rows[n].second=iind; n++;
     119             :   }
     120             : 
     121             :   // Now do the sort and clear all the stored values ready for recompute
     122        2512 :   std::sort( rows.begin(), rows.end() );
     123             :   // This is to make this action consistent with what in other matrix actions
     124        2512 :   unsigned start_n = getPntrToArgument(0)->getShape()[0];
     125             :   // And setup the lowest indices, which are the ones we need to calculate
     126       16560 :   for(unsigned i=0; i<number; ++i) {
     127       14048 :     indices[i+1] = start_n + rows[nind-1-i].second;
     128       14048 :     if( lowest ) indices[i+1] = start_n + rows[i].second;
     129             :   }
     130        2512 : }
     131             : 
     132       14048 : void Neighbors::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
     133       14048 :   myvals.addValue( getConstPntrToComponent(0)->getPositionInStream(), 1.0 );
     134       14048 : }
     135             : 
     136             : }
     137             : }