Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-2020 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_gridtools_FunctionOfGrid_h
      23             : #define __PLUMED_gridtools_FunctionOfGrid_h
      24             : 
      25             : #include "ActionWithGrid.h"
      26             : #include "function/Custom.h"
      27             : #include "tools/Matrix.h"
      28             : 
      29             : namespace PLMD {
      30             : namespace gridtools {
      31             : 
      32             : template <class T>
      33             : class FunctionOfGrid : public ActionWithGrid {
      34             : private:
      35             : /// The function that is being computed
      36             :   T myfunc;
      37             : public:
      38             :   static void registerKeywords(Keywords&);
      39             :   explicit FunctionOfGrid(const ActionOptions&);
      40             : /// This does setup required on first step
      41             :   void setupOnFirstStep( const bool incalc ) override ;
      42             : /// Get the number of derivatives for this action
      43             :   unsigned getNumberOfDerivatives() override ;
      44             : /// Get the label to write in the graph
      45           0 :   std::string writeInGraph() const override { return myfunc.getGraphInfo( getName() ); }
      46             : /// Get the underlying names
      47             :   std::vector<std::string> getGridCoordinateNames() const override ;
      48             : /// Get the underlying grid coordinates object
      49             :   const GridCoordinatesObject& getGridCoordinatesObject() const override ;
      50             : /// Calculate the function
      51             :   void performTask( const unsigned& current, MultiValue& myvals ) const override ;
      52             : ///
      53             :   void gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
      54             :                           const unsigned& bufstart, std::vector<double>& buffer ) const override ;
      55             : /// Add the forces
      56             :   void apply() override;
      57             : };
      58             : 
      59             : template <class T>
      60        1047 : void FunctionOfGrid<T>::registerKeywords(Keywords& keys ) {
      61        1047 :   ActionWithGrid::registerKeywords(keys); std::string name = keys.getDisplayName();
      62        1047 :   std::size_t und=name.find("_GRID"); keys.setDisplayName( name.substr(0,und) );
      63        2094 :   keys.reserve("compulsory","PERIODIC","if the output of your function is periodic then you should specify the periodicity of the function.  If the output is not periodic you must state this using PERIODIC=NO");
      64        1047 :   T tfunc; tfunc.registerKeywords( keys ); if( typeid(tfunc)==typeid(function::Custom()) ) keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
      65        2094 :   if( keys.getDisplayName()=="INTEGRATE") {
      66         294 :     keys.setValueDescription("scalar","the numerical integral of the input function over its whole domain");
      67         294 :     keys.addInputKeyword("compulsory","ARG","grid","the label of the function on a grid that is being integrated");
      68        1800 :   } else if( keys.getDisplayName()=="SUM") {
      69          58 :     keys.setValueDescription("scalar","the sum of the value of the function over all the grid points where it has been evaluated");
      70          58 :     keys.addInputKeyword("compulsory","ARG","grid","the label of the function on a grid from which we are computing a sum");
      71        1742 :   } else if( keys.outputComponentExists(".#!value") ) {
      72        1742 :     keys.setValueDescription("grid","the grid obtained by doing an element-wise application of " + keys.getOutputComponentDescription(".#!value") + " to the input grid");
      73        1742 :     keys.addInputKeyword("compulsory","ARG","scalar/grid","the labels of the scalars and functions on a grid that we are using to compute the required function");
      74             :   }
      75        1918 : }
      76             : 
      77             : template <class T>
      78         523 : FunctionOfGrid<T>::FunctionOfGrid(const ActionOptions&ao):
      79             :   Action(ao),
      80         523 :   ActionWithGrid(ao)
      81             : {
      82         523 :   if( getNumberOfArguments()==0 ) error("found no arguments");
      83             :   // This will require a fix
      84         523 :   if( getPntrToArgument(0)->getRank()==0 || !getPntrToArgument(0)->hasDerivatives() ) error("first input to this action must be a grid");
      85             :   // Get the shape of the input grid
      86         523 :   std::vector<unsigned> shape( getPntrToArgument(0)->getShape() );
      87         937 :   for(unsigned i=1; i<getNumberOfArguments(); ++i ) {
      88         414 :     if( getPntrToArgument(i)->getRank()==0 ) continue;
      89         292 :     std::vector<unsigned> s( getPntrToArgument(i)->getShape() );
      90         292 :     if( s.size()!=shape.size() ) error("mismatch between dimensionalities of input grids");
      91             :   }
      92             :   // Read the input and do some checks
      93         523 :   myfunc.read( this );
      94             :   // Check we are not calculating an integral
      95         523 :   if( myfunc.zeroRank() ) { shape.resize(0); }
      96             :   // Check that derivatives are available
      97          90 :   if( !myfunc.derivativesImplemented() ) error("derivatives have not been implemended for " + getName() );
      98             :   // Get the names of the components
      99         523 :   std::vector<std::string> components( keywords.getOutputComponents() );
     100             :   // Create the values to hold the output
     101        1046 :   if( components.size()!=1 || components[0]!=".#!value" ) error("functions of grid should only output one grid");
     102         523 :   addValueWithDerivatives( shape );
     103             :   // Set the periodicities of the output components
     104         523 :   myfunc.setPeriodicityForOutputs( this );
     105             :   // Check if we can turn off the derivatives when they are zero
     106         433 :   if( myfunc.getDerivativeZeroIfValueIsZero() )  {
     107         492 :     for(int i=0; i<getNumberOfComponents(); ++i) getPntrToComponent(i)->setDerivativeIsZeroWhenValueIsZero();
     108             :   }
     109         523 :   setupOnFirstStep( false );
     110        1046 : }
     111             : 
     112             : template <class T>
     113        1045 : void FunctionOfGrid<T>::setupOnFirstStep( const bool incalc ) {
     114             :   double volume = 1.0;
     115        1045 :   const GridCoordinatesObject& mygrid = getGridCoordinatesObject();
     116        1045 :   unsigned npoints = getPntrToArgument(0)->getNumberOfValues();
     117        2090 :   if( mygrid.getGridType()=="flat" ) {
     118         999 :     std::vector<unsigned> shape( getGridCoordinatesObject().getNbin(true) );
     119        1794 :     for(unsigned i=1; i<getNumberOfArguments(); ++i ) {
     120         795 :       if( getPntrToArgument(i)->getRank()==0 ) continue;
     121         573 :       std::vector<unsigned> s( getPntrToArgument(i)->getShape() );
     122        1160 :       for(unsigned j=0; j<shape.size(); ++j) {
     123         587 :         if( shape[j]!=s[j] ) error("mismatch between sizes of input grids");
     124             :       }
     125             :     }
     126        1998 :     for(int i=0; i<getNumberOfComponents(); ++i) {
     127         999 :       if( getPntrToComponent(i)->getRank()>0 ) getPntrToComponent(i)->setShape(shape);
     128             :     }
     129         999 :     std::vector<double> vv( getGridCoordinatesObject().getGridSpacing() );
     130        1081 :     volume=vv[0]; for(unsigned i=1; i<vv.size(); ++i) volume *=vv[i];
     131          14 :   } else volume=4*pi / static_cast<double>( npoints );
     132             :   // This resizes the scalars
     133        2090 :   for(int i=0; i<getNumberOfComponents(); ++i) {
     134        1045 :     if( getPntrToComponent(i)->getRank()==0 ) getPntrToComponent(i)->resizeDerivatives( npoints );
     135             :   }
     136        1045 :   if( getName()=="SUM_GRID" ) volume = 1.0;
     137             :   // This sets the prefactor to the volume which converts integrals to sums
     138        1045 :   myfunc.setup( this ); myfunc.setPrefactor( this, volume );
     139        1045 : }
     140             : 
     141             : template <class T>
     142       15460 : const GridCoordinatesObject& FunctionOfGrid<T>::getGridCoordinatesObject() const {
     143       15460 :   ActionWithGrid* ag=ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
     144       15460 :   plumed_assert( ag ); return ag->getGridCoordinatesObject();
     145             : }
     146             : 
     147             : template <class T>
     148         198 : std::vector<std::string> FunctionOfGrid<T>::getGridCoordinateNames() const {
     149         198 :   ActionWithGrid* ag=ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
     150         198 :   plumed_assert( ag ); return ag->getGridCoordinateNames();
     151             : }
     152             : 
     153             : template <class T>
     154        1814 : unsigned FunctionOfGrid<T>::getNumberOfDerivatives() {
     155         514 :   if( myfunc.zeroRank() ) return getPntrToArgument(0)->getNumberOfValues();
     156        1300 :   unsigned nder = getGridCoordinatesObject().getDimension();
     157        1300 :   return getGridCoordinatesObject().getDimension() + getNumberOfArguments() - myfunc.getArgStart();
     158             : }
     159             : 
     160             : template <class T>
     161      974814 : void FunctionOfGrid<T>::performTask( const unsigned& current, MultiValue& myvals ) const {
     162      974814 :   unsigned argstart=myfunc.getArgStart(); std::vector<double> args( getNumberOfArguments() - argstart );
     163     2629182 :   for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
     164     1654368 :     if( getPntrToArgument(i)->getRank()==0 ) args[i-argstart]=getPntrToArgument(i)->get();
     165     1110136 :     else args[i-argstart] = getPntrToArgument(i)->get(current);
     166             :   }
     167             :   // Calculate the function and its derivatives
     168      974814 :   std::vector<double> vals(1); Matrix<double> derivatives( 1, getNumberOfArguments()-argstart );
     169      974814 :   myfunc.calc( this, args, vals, derivatives ); unsigned np = myvals.getTaskIndex();
     170             :   // And set the values and derivatives
     171      974814 :   unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream();
     172      974814 :   myvals.addValue( ostrn, vals[0] );
     173       23665 :   if( !myfunc.zeroRank() ) {
     174             :     // Add the derivatives for a grid
     175     2581852 :     for(unsigned j=argstart; j<getNumberOfArguments(); ++j) {
     176             :       // We store all the derivatives of all the input values - i.e. the grid points these are used in apply
     177     1630703 :       myvals.addDerivative( ostrn, getConstPntrToComponent(0)->getRank()+j-argstart, derivatives(0,j-argstart) );
     178             :       // And now we calculate the derivatives of the value that is stored on the grid correctly so that we can interpolate functions
     179     1630703 :       if( getPntrToArgument(j)->getRank()!=0 ) {
     180     3413677 :         for(unsigned k=0; k<getPntrToArgument(j)->getRank(); ++k) myvals.addDerivative( ostrn, k, derivatives(0,j-argstart)*getPntrToArgument(j)->getGridDerivative( np, k ) );
     181             :       }
     182             :     }
     183      951149 :     unsigned nderivatives = getConstPntrToComponent(0)->getNumberOfGridDerivatives();
     184     4575808 :     for(unsigned j=0; j<nderivatives; ++j) myvals.updateIndex( ostrn, j );
     185       23665 :   } else if( !doNotCalculateDerivatives() ) {
     186             :     // These are the derivatives of the integral
     187        8161 :     myvals.addDerivative( ostrn, current, derivatives(0,0) ); myvals.updateIndex( ostrn, current );
     188             :   }
     189      974814 : }
     190             : 
     191             : template <class T>
     192      951149 : void FunctionOfGrid<T>::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
     193             :     const unsigned& bufstart, std::vector<double>& buffer ) const {
     194      951149 :   if( getConstPntrToComponent(0)->getRank()>0 && getConstPntrToComponent(0)->hasDerivatives() ) {
     195             :     plumed_dbg_assert( getNumberOfComponents()==1 && valindex==0 );
     196      951149 :     unsigned nder = getConstPntrToComponent(0)->getNumberOfGridDerivatives();
     197      951149 :     unsigned ostr = getConstPntrToComponent(0)->getPositionInStream();
     198      951149 :     unsigned kp = bufstart + code*(1+nder); buffer[kp] += myvals.get( ostr );
     199     4575808 :     for(unsigned i=0; i<nder; ++i) buffer[kp + 1 + i] += myvals.getDerivative( ostr, i );
     200           0 :   } else ActionWithVector::gatherStoredValue( valindex, code, myvals, bufstart, buffer );
     201      951149 : }
     202             : 
     203             : template <class T>
     204        6067 : void FunctionOfGrid<T>::apply() {
     205        6561 :   if( doNotCalculateDerivatives() || !getPntrToComponent(0)->forcesWereAdded() ) return;
     206             : 
     207             :   // This applies forces for the integral
     208         494 :   if( myfunc.zeroRank() ) { ActionWithVector::apply(); return; }
     209             : 
     210             :   // Work out how to deal with arguments
     211             :   unsigned nscalars=0, argstart=myfunc.getArgStart();
     212        1870 :   for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
     213        1245 :     if( getPntrToArgument(i)->getRank()==0 ) { nscalars++; }
     214             :   }
     215             : 
     216         625 :   std::vector<double> totv(nscalars,0); Value* outval=getPntrToComponent(0);
     217        7575 :   for(unsigned i=0; i<outval->getNumberOfValues(); ++i) {
     218             :     nscalars=0;
     219       19845 :     for(unsigned j=argstart; j<getNumberOfArguments(); ++j) {
     220       12895 :       double fforce = outval->getForce(i);
     221       12895 :       if( getPntrToArgument(j)->getRank()==0 ) {
     222        3205 :         totv[nscalars] += fforce*outval->getGridDerivative( i, outval->getRank()+j ); nscalars++;
     223             :       } else {
     224        9690 :         double vval = outval->getGridDerivative( i, outval->getRank()+j  );
     225        9690 :         getPntrToArgument(j)->addForce( i, fforce*vval );
     226             :       }
     227             :     }
     228             :   }
     229             :   nscalars=0;
     230        1870 :   for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
     231        1245 :     if( getPntrToArgument(i)->getRank()==0 ) { getPntrToArgument(i)->addForce( 0, totv[nscalars] ); nscalars++; }
     232             :   }
     233             : 
     234             : }
     235             : 
     236             : }
     237             : }
     238             : #endif