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             : The following input illustrates how to use this action in tandem with [DISTANCE_MATRIX](DISTANCE_MATRIX.md) to find the six
      30             : nearest atoms to each of the first 100 atoms in the input file:
      31             : 
      32             : ```plumed
      33             : d1: DISTANCE_MATRIX GROUP=1-100
      34             : n: NEIGHBORS ARG=d1 NLOWEST=6
      35             : ```
      36             : 
      37             : Alternatively, if you would like to use a [CONTACT_MATRIX](CONTACT_MATRIX.md) to do something similar you would do the following:
      38             : 
      39             : ```plumed
      40             : c1: CONTACT_MATRIX GROUP=1-100 SWITCH={RATIONAL R_0=0.5}
      41             : n: NEIGHBORS ARG=c1 NHIGHEST=6
      42             : ```
      43             : 
      44             : This command is useful for implementing alternatives to the symmatry functions that are defined by the shortcuts
      45             : in the module symfunc.  For example, suppose that you want to calculate a variant on the [TETRAHEDRAL](TETRAHEDRAL.md) symmetry function.
      46             : In this variant on the CV the coordination sphere around each central atom is not defined using a switching function.  Instad
      47             : this coordination sphere contains only the four nearest atoms.  You can implement this CV by using the following input:
      48             : 
      49             : ```plumed
      50             : d1: DISTANCE_MATRIX GROUP=1-100 COMPONENTS
      51             : n: NEIGHBORS ARG=d1.w NLOWEST=4
      52             : f: CUSTOM ARG=n,d1.x,d1.y,d1.z,d1.w VAR=w,x,y,z,r PERIODIC=NO FUNC=w*(((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3)
      53             : ones: ONES SIZE=100
      54             : ucv: MATRIX_VECTOR_PRODUCT ARG=f,ones
      55             : cv: CUSTOM ARG=ucv PERIODIC=NO FUNC=x/4
      56             : ```
      57             : 
      58             : */
      59             : //+ENDPLUMEDOC
      60             : 
      61             : namespace PLMD {
      62             : namespace adjmat {
      63             : 
      64             : class Neighbors : public ActionWithMatrix {
      65             :   bool lowest;
      66             :   unsigned number;
      67             : public:
      68             :   static void registerKeywords( Keywords& keys );
      69             :   explicit Neighbors(const ActionOptions&);
      70             :   unsigned getNumberOfDerivatives() override;
      71           8 :   unsigned getNumberOfColumns() const override {
      72           8 :     return number;
      73             :   }
      74           0 :   bool canBeAfterInChain( ActionWithVector* av ) {
      75           0 :     return av->getLabel()!=(getPntrToArgument(0)->getPntrToAction())->getLabel();
      76             :   }
      77             :   void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const ;
      78             :   void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override;
      79        2512 :   void runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const override {}
      80             :   void turnOnDerivatives() override ;
      81             : };
      82             : 
      83             : PLUMED_REGISTER_ACTION(Neighbors,"NEIGHBORS")
      84             : 
      85           7 : void Neighbors::registerKeywords( Keywords& keys ) {
      86           7 :   ActionWithMatrix::registerKeywords( keys );
      87          14 :   keys.addInputKeyword("compulsory","ARG","matrix","the label of an adjacency/distance matrix that will be used to find the nearest neighbors");
      88           7 :   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");
      89           7 :   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");
      90          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");
      91           7 : }
      92             : 
      93           3 : Neighbors::Neighbors(const ActionOptions&ao):
      94             :   Action(ao),
      95           3 :   ActionWithMatrix(ao) {
      96           3 :   if( getNumberOfArguments()!=1 ) {
      97           0 :     error("found wrong number of arguments in input");
      98             :   }
      99           3 :   if( getPntrToArgument(0)->getRank()!=2 ) {
     100           0 :     error("input argument should be a matrix");
     101             :   }
     102           3 :   getPntrToArgument(0)->buildDataStore();
     103             : 
     104             :   unsigned nlow;
     105           3 :   parse("NLOWEST",nlow);
     106             :   unsigned nhigh;
     107           3 :   parse("NHIGHEST",nhigh);
     108           3 :   if( nlow==0 && nhigh==0 ) {
     109           0 :     error("missing NLOWEST or NHIGHEST keyword one of these two keywords must be set in input");
     110             :   }
     111           3 :   if( nlow>0 && nhigh>0 ) {
     112           0 :     error("should only be one of NLOWEST or NHIGHEST set in input");
     113             :   }
     114           3 :   if( nlow>0 ) {
     115           3 :     number=nlow;
     116           3 :     lowest=true;
     117           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);
     118             :   }
     119           3 :   if( nhigh>0 ) {
     120           0 :     number=nhigh;
     121           0 :     lowest=false;
     122           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);
     123             :   }
     124             : 
     125             :   // And get the shape
     126           3 :   std::vector<unsigned> shape( getPntrToArgument(0)->getShape() );
     127           3 :   addValue( shape );
     128           3 :   setNotPeriodic();
     129           3 :   checkRead();
     130           3 : }
     131             : 
     132           0 : void Neighbors::turnOnDerivatives() {
     133           0 :   ActionWithValue::turnOnDerivatives();
     134           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");
     135           0 : }
     136             : 
     137           0 : unsigned Neighbors::getNumberOfDerivatives() {
     138           0 :   return 0;
     139             : }
     140             : 
     141        2512 : void Neighbors::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
     142             :   const Value* wval = getPntrToArgument(0);
     143        2512 :   unsigned nbonds = wval->getRowLength( task_index ), ncols = wval->getShape()[1];
     144        2512 :   if( indices.size()!=1+number ) {
     145          26 :     indices.resize( 1 + number );
     146             :   }
     147        2512 :   myvals.setSplitIndex(1+number);
     148             : 
     149             :   unsigned nind=0;
     150      180890 :   for(unsigned i=0; i<nbonds; ++i) {
     151      178378 :     unsigned ipos = ncols*task_index + wval->getRowIndex( task_index, i );
     152      178378 :     double weighti = wval->get( ipos );
     153      178378 :     if( weighti<epsilon ) {
     154         512 :       continue ;
     155             :     }
     156      177866 :     nind++;
     157             :   }
     158        2512 :   if( number>nind ) {
     159           0 :     plumed_merror("not enough matrix elements were stored");
     160             :   }
     161             : 
     162             :   // Now build vectors for doing sorting
     163        2512 :   std::vector<std::pair<double,unsigned> > rows( nind );
     164             :   unsigned n=0;
     165      180890 :   for(unsigned i=0; i<nbonds; ++i) {
     166             :     unsigned iind = wval->getRowIndex( task_index, i );
     167      178378 :     unsigned ipos = ncols*task_index + iind;
     168      178378 :     double weighti = wval->get( ipos );
     169      178378 :     if( weighti<epsilon ) {
     170         512 :       continue ;
     171             :     }
     172      177866 :     rows[n].first=weighti;
     173      177866 :     rows[n].second=iind;
     174      177866 :     n++;
     175             :   }
     176             : 
     177             :   // Now do the sort and clear all the stored values ready for recompute
     178        2512 :   std::sort( rows.begin(), rows.end() );
     179             :   // This is to make this action consistent with what in other matrix actions
     180        2512 :   unsigned start_n = getPntrToArgument(0)->getShape()[0];
     181             :   // And setup the lowest indices, which are the ones we need to calculate
     182       16560 :   for(unsigned i=0; i<number; ++i) {
     183       14048 :     indices[i+1] = start_n + rows[nind-1-i].second;
     184       14048 :     if( lowest ) {
     185       14048 :       indices[i+1] = start_n + rows[i].second;
     186             :     }
     187             :   }
     188        2512 : }
     189             : 
     190       14048 : void Neighbors::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
     191       14048 :   myvals.addValue( getConstPntrToComponent(0)->getPositionInStream(), 1.0 );
     192       14048 : }
     193             : 
     194             : }
     195             : }