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             : #ifndef __PLUMED_function_FunctionWithSingleArgument_h
      23             : #define __PLUMED_function_FunctionWithSingleArgument_h
      24             : 
      25             : #include "core/ActionWithArguments.h"
      26             : #include "core/ActionWithValue.h"
      27             : #include "core/ActionRegister.h"
      28             : #include "FunctionSetup.h"
      29             : 
      30             : namespace PLMD {
      31             : namespace function {
      32             : 
      33             : template <class T>
      34             : class FunctionWithSingleArgument :
      35             :   public ActionWithValue,
      36             :   public ActionWithArguments {
      37             : private:
      38             :   T f;
      39             :   bool ismatrix;
      40             : public:
      41             :   static void registerKeywords(Keywords& keys);
      42             :   explicit FunctionWithSingleArgument(const ActionOptions&);
      43             :   /// Get the label to write in the graph
      44           0 :   std::string writeInGraph() const override {
      45           0 :     std::size_t und = getName().find_last_of("_");
      46           0 :     return getName().substr(0,und);
      47             :   }
      48             :   unsigned getNumberOfDerivatives() override ;
      49             :   void calculate() override ;
      50             :   void apply() override ;
      51             : };
      52             : 
      53             : template <class T>
      54         123 : void FunctionWithSingleArgument<T>::registerKeywords(Keywords& keys) {
      55         123 :   Action::registerKeywords(keys);
      56         123 :   ActionWithValue::registerKeywords(keys);
      57         123 :   ActionWithArguments::registerKeywords(keys);
      58         123 :   std::string name = keys.getDisplayName();
      59         123 :   std::size_t und=name.find("_ONEARG");
      60         123 :   keys.setDisplayName( name.substr(0,und) );
      61         246 :   keys.addInputKeyword("compulsory","ARG","vector/matrix","the vector/matrix that is being used in input");
      62         123 :   T::registerKeywords( keys );
      63         123 : }
      64             : 
      65             : template <class T>
      66          53 : FunctionWithSingleArgument<T>::FunctionWithSingleArgument(const ActionOptions&ao):
      67             :   Action(ao),
      68             :   ActionWithValue(ao),
      69             :   ActionWithArguments(ao),
      70          53 :   ismatrix(false) {
      71             : 
      72          53 :   if( getNumberOfArguments()!=1 ) {
      73           0 :     error("should only be one argument in input for this action");
      74             :   }
      75          53 :   ismatrix = getPntrToArgument(0)->getRank()==2 && !getPntrToArgument(0)->hasDerivatives();
      76             : 
      77          53 :   std::vector<std::size_t> shape;
      78          53 :   FunctionData<T>::setup( f, keywords.getOutputComponents(), shape, true, this );
      79          53 : }
      80             : 
      81             : template <class T>
      82         118 : unsigned FunctionWithSingleArgument<T>::getNumberOfDerivatives() {
      83         118 :   return getPntrToArgument(0)->getNumberOfStoredValues();
      84             : }
      85             : 
      86             : template <class T>
      87        5553 : void FunctionWithSingleArgument<T>::calculate() {
      88             :   const Value* myarg = getPntrToArgument(0);
      89        5553 :   std::vector<double> args( myarg->getNumberOfStoredValues() );
      90        5553 :   if( ismatrix ) {
      91             :     unsigned nvals = 0;
      92         266 :     for(unsigned i=0; i<myarg->getShape()[0]; ++i) {
      93         247 :       nvals += myarg->getRowLength(i);
      94             :     }
      95          19 :     args.resize( nvals );
      96             :     unsigned k=0;
      97          19 :     unsigned ncols = myarg->getNumberOfColumns();
      98         266 :     for(unsigned i=0; i<myarg->getShape()[0]; ++i) {
      99         969 :       for(unsigned j=0; j<myarg->getRowLength(i); ++j) {
     100         722 :         args[k] = myarg->get(i*ncols + j, false);
     101         722 :         ++k;
     102             :       }
     103             :     }
     104             :   } else {
     105     8198545 :     for(unsigned i=0; i<args.size(); ++i) {
     106     8193011 :       args[i] = myarg->get( i );
     107             :     }
     108             :   }
     109        5553 :   std::vector<double> vals( getNumberOfComponents() ), deriv( getNumberOfComponents()*args.size() );
     110             :   auto funcout = FunctionOutput::create( getNumberOfComponents(),
     111             :                                          vals.data(),
     112             :                                          args.size(),
     113             :                                          deriv.data() );
     114        5553 :   T::calc( f,
     115        5553 :            doNotCalculateDerivatives(),
     116           0 :            View<const double>(args.data(), args.size()),
     117             :            funcout );
     118       11329 :   for(unsigned i=0; i<getNumberOfComponents(); ++i) {
     119        5776 :     Value* myval=getPntrToComponent(i);
     120        5776 :     myval->set( vals[i] );
     121             : 
     122        5776 :     if( doNotCalculateDerivatives() ) {
     123        1351 :       continue ;
     124             :     }
     125        4425 :     if( ismatrix ) {
     126          19 :       myval->resizeDerivatives( args.size() );
     127             :     }
     128      198678 :     for(unsigned j=0; j<args.size(); ++j) {
     129      194253 :       myval->setDerivative( j, funcout.derivs[i][j] );
     130             :     }
     131             :   }
     132        5553 : }
     133             : 
     134             : template <class T>
     135        5537 : void FunctionWithSingleArgument<T>::apply() {
     136             :   // Forces may have been added to any of the components (e.g. only to label.moment-3)
     137             :   bool hasforces=false;
     138        6878 :   for(unsigned i=0; i<getNumberOfComponents(); ++i) {
     139        5755 :     if( getPntrToComponent(i)->forcesWereAdded() ) {
     140             :       hasforces=true;
     141             :       break;
     142             :     }
     143             :   }
     144        5537 :   if( !hasforces ) {
     145             :     return ;
     146             :   }
     147             : 
     148             :   Value* myarg = getPntrToArgument(0);
     149        8839 :   for(unsigned i=0; i<getNumberOfComponents(); ++i) {
     150        4425 :     const Value* myval = getPntrToComponent(i);
     151        4425 :     if( !myval->forcesWereAdded() ) {
     152          11 :       continue;
     153             :     }
     154             :     double force = myval->getForce();
     155        4414 :     if( ismatrix ) {
     156             :       unsigned k=0;
     157          19 :       unsigned ncols = myarg->getNumberOfColumns();
     158         266 :       for(unsigned n=0; n<myarg->getShape()[0]; ++n) {
     159             :         const unsigned rowlen = myarg->getRowLength(n);
     160         969 :         for(unsigned j=0; j<rowlen; ++j) {
     161         722 :           myarg->addForce( n*ncols+j, force*myval->getDerivative(k), false );
     162             :           ++k;
     163             :         }
     164             :       }
     165             :     } else {
     166      197872 :       for(unsigned j=0; j<myarg->getNumberOfStoredValues(); ++j) {
     167      193477 :         myarg->addForce( j, force*myval->getDerivative(j) );
     168             :       }
     169             :     }
     170             :   }
     171             : }
     172             : 
     173             : }
     174             : }
     175             : #endif