Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 : Copyright (c) 2011-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 "core/ActionWithMatrix.h" 23 : #include "core/ActionRegister.h" 24 : 25 : //+PLUMEDOC MCOLVAR VORONOI 26 : /* 27 : Do a voronoi analysis 28 : 29 : \par Examples 30 : 31 : */ 32 : //+ENDPLUMEDOC 33 : 34 : namespace PLMD { 35 : namespace matrixtools { 36 : 37 : class Voronoi : public ActionWithMatrix { 38 : public: 39 : static void registerKeywords( Keywords& keys ); 40 : explicit Voronoi(const ActionOptions&); 41 : void prepare() override ; 42 0 : unsigned getNumberOfDerivatives() override { return 0; } 43 12 : unsigned getNumberOfColumns() const override { return getConstPntrToComponent(0)->getShape()[1]; } 44 1036 : void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const override {} 45 0 : void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override {} 46 1036 : void runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const override {} 47 : void gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals, 48 : const unsigned& bufstart, std::vector<double>& buffer ) const override ; 49 : }; 50 : 51 : PLUMED_REGISTER_ACTION(Voronoi,"VORONOI") 52 : 53 11 : void Voronoi::registerKeywords( Keywords& keys ) { 54 11 : ActionWithMatrix::registerKeywords(keys); keys.use("ARG"); 55 11 : keys.setValueDescription("a matrix in which element ij is equal to one if the ij component of the input matrix is lower than all the ik elements of the matrix where k is not j and zero otherwise"); 56 11 : } 57 : 58 6 : Voronoi::Voronoi(const ActionOptions&ao): 59 : Action(ao), 60 6 : ActionWithMatrix(ao) 61 : { 62 6 : if( getNumberOfArguments()!=1 ) error("should be one arguments to this action, a matrix"); 63 6 : if( getPntrToArgument(0)->getRank()!=2 || getPntrToArgument(0)->hasDerivatives() ) error("argument to this action should be a matrix"); 64 6 : if( getPntrToArgument(0)->getShape()[1]>getPntrToArgument(0)->getShape()[0] ) warning("would expect number of columns in matrix to exceed number of rows"); 65 6 : getPntrToArgument(0)->buildDataStore(); std::vector<unsigned> shape( getPntrToArgument(0)->getShape() ); 66 6 : addValue( shape ); setNotPeriodic(); getPntrToComponent(0)->buildDataStore(); 67 6 : } 68 : 69 6 : void Voronoi::prepare() { 70 6 : Value* myval = getPntrToComponent(0); 71 6 : if( myval->getShape()[0]==getPntrToArgument(0)->getShape()[0] && myval->getShape()[1]==getPntrToArgument(0)->getShape()[1] ) return; 72 6 : std::vector<unsigned> shape( getPntrToArgument(0)->getShape() ); myval->setShape(shape); 73 : } 74 : 75 1036 : void Voronoi::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals, 76 : const unsigned& bufstart, std::vector<double>& buffer ) const { 77 1036 : Value* arg0 = getPntrToArgument(0); unsigned nv = 0; std::size_t cc=code; double minmax = arg0->get( cc*arg0->getShape()[1] ); 78 251175 : for(unsigned i=0; i<arg0->getShape()[1]; ++i) { 79 250139 : double value = arg0->get( code*arg0->getShape()[1] + i ); 80 250139 : if( value<minmax ) { minmax = value; nv = i; } 81 : } 82 1036 : buffer[bufstart + code*arg0->getShape()[1] + nv] = 1; 83 1036 : } 84 : 85 : } 86 : }