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