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 "MatrixOperationBase.h" 23 : 24 : namespace PLMD { 25 : namespace matrixtools { 26 : 27 363 : void MatrixOperationBase::registerKeywords( Keywords& keys ) { 28 363 : Action::registerKeywords( keys ); ActionWithArguments::registerKeywords( keys ); ActionWithValue::registerKeywords( keys ); 29 726 : keys.use("ARG"); keys.remove("NUMERICAL_DERIVATIVES"); 30 726 : keys.add("optional","MATRIX","the input matrix (can use ARG instead)"); 31 363 : } 32 : 33 211 : MatrixOperationBase::MatrixOperationBase(const ActionOptions&ao): 34 : Action(ao), 35 : ActionWithArguments(ao), 36 211 : ActionWithValue(ao) 37 : { 38 211 : if( getNumberOfArguments()==0 ) { 39 0 : std::vector<Value*> args; parseArgumentList("MATRIX",args); requestArguments(args); 40 : } 41 211 : if( getNumberOfArguments()!=1 ) error("should only be one argument to this action"); 42 211 : if( getPntrToArgument(0)->getRank()!=2 || getPntrToArgument(0)->hasDerivatives() ) { 43 17 : if( getName()=="TRANSPOSE" ) { 44 17 : if (getPntrToArgument(0)->getRank()!=1 || getPntrToArgument(0)->hasDerivatives() ) error("input to this argument should be a matrix or vector"); 45 0 : } else error("input to this argument should be a matrix"); 46 : } 47 211 : getPntrToArgument(0)->buildDataStore(); 48 211 : } 49 : 50 131 : void MatrixOperationBase::retrieveFullMatrix( Matrix<double>& mymatrix ) { 51 131 : if( mymatrix.nrows()!=getPntrToArgument(0)->getShape()[0] || mymatrix.nrows()!=getPntrToArgument(0)->getShape()[1] ) { 52 0 : mymatrix.resize( getPntrToArgument(0)->getShape()[0], getPntrToArgument(0)->getShape()[1] ); 53 : } 54 131 : unsigned nedge; getPntrToArgument(0)->retrieveEdgeList( nedge, pairs, vals ); mymatrix=0; 55 : bool symmetric = getPntrToArgument(0)->isSymmetric(); 56 373270 : for(unsigned i=0; i<nedge; ++i ) { 57 373139 : mymatrix( pairs[i].first, pairs[i].second ) = vals[i]; 58 373139 : if( symmetric ) mymatrix( pairs[i].second, pairs[i].first ) = vals[i]; 59 : } 60 131 : } 61 : 62 : 63 1448 : void MatrixOperationBase::apply() { 64 1448 : if( doNotCalculateDerivatives() ) return; 65 : 66 : bool forces=false; 67 1546 : for(int i=0; i<getNumberOfComponents(); ++i) { 68 1546 : if( getPntrToComponent(i)->forcesWereAdded() ) { forces=true; break; } 69 : } 70 1444 : if( !forces ) return; 71 : 72 : Value* mat = getPntrToArgument(0); 73 : unsigned ncols=mat->getNumberOfColumns(); 74 147286 : for(unsigned i=0; i<mat->getShape()[0]; ++i) { 75 : unsigned ncol = mat->getRowLength(i); 76 4282909 : for(unsigned j=0; j<ncol; ++j) mat->addForce( i*ncols+j, getForceOnMatrixElement( i, mat->getRowIndex(i,j) ), false ); 77 : } 78 : } 79 : 80 : 81 : } 82 : }