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 {
      46           0 :     return myfunc.getGraphInfo( getName() );
      47             :   }
      48             : /// Get the underlying names
      49             :   std::vector<std::string> getGridCoordinateNames() const override ;
      50             : /// Get the underlying grid coordinates object
      51             :   const GridCoordinatesObject& getGridCoordinatesObject() const override ;
      52             : /// Calculate the function
      53             :   void performTask( const unsigned& current, MultiValue& myvals ) const override ;
      54             : ///
      55             :   void gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
      56             :                           const unsigned& bufstart, std::vector<double>& buffer ) const override ;
      57             : /// Add the forces
      58             :   void apply() override;
      59             : };
      60             : 
      61             : template <class T>
      62        1047 : void FunctionOfGrid<T>::registerKeywords(Keywords& keys ) {
      63        1047 :   ActionWithGrid::registerKeywords(keys);
      64        1047 :   keys.use("ARG");
      65        1047 :   std::string name = keys.getDisplayName();
      66        1047 :   std::size_t und=name.find("_GRID");
      67        1047 :   keys.setDisplayName( name.substr(0,und) );
      68        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");
      69         871 :   T tfunc;
      70        1047 :   tfunc.registerKeywords( keys );
      71        1047 :   if( typeid(tfunc)==typeid(function::Custom()) ) {
      72           0 :     keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
      73             :   }
      74        2094 :   if( keys.getDisplayName()=="INTEGRATE") {
      75         294 :     keys.setValueDescription("the numerical integral of the input function over its whole domain");
      76        1800 :   } else if( keys.outputComponentExists(".#!value") ) {
      77        1800 :     keys.setValueDescription("the grid obtained by doing an element-wise application of " + keys.getOutputComponentDescription(".#!value") + " to the input grid");
      78             :   }
      79        1918 : }
      80             : 
      81             : template <class T>
      82         523 : FunctionOfGrid<T>::FunctionOfGrid(const ActionOptions&ao):
      83             :   Action(ao),
      84         523 :   ActionWithGrid(ao) {
      85         523 :   if( getNumberOfArguments()==0 ) {
      86           0 :     error("found no arguments");
      87             :   }
      88             :   // This will require a fix
      89         523 :   if( getPntrToArgument(0)->getRank()==0 || !getPntrToArgument(0)->hasDerivatives() ) {
      90           0 :     error("first input to this action must be a grid");
      91             :   }
      92             :   // Get the shape of the input grid
      93         523 :   std::vector<unsigned> shape( getPntrToArgument(0)->getShape() );
      94         937 :   for(unsigned i=1; i<getNumberOfArguments(); ++i ) {
      95         414 :     if( getPntrToArgument(i)->getRank()==0 ) {
      96         122 :       continue;
      97             :     }
      98         292 :     std::vector<unsigned> s( getPntrToArgument(i)->getShape() );
      99         292 :     if( s.size()!=shape.size() ) {
     100           0 :       error("mismatch between dimensionalities of input grids");
     101             :     }
     102             :   }
     103             :   // Read the input and do some checks
     104         523 :   myfunc.read( this );
     105             :   // Check we are not calculating an integral
     106         523 :   if( myfunc.zeroRank() ) {
     107          90 :     shape.resize(0);
     108             :   }
     109             :   // Check that derivatives are available
     110          90 :   if( !myfunc.derivativesImplemented() ) {
     111             :     error("derivatives have not been implemended for " + getName() );
     112             :   }
     113             :   // Get the names of the components
     114         523 :   std::vector<std::string> components( keywords.getOutputComponents() );
     115             :   // Create the values to hold the output
     116        1046 :   if( components.size()!=1 || components[0]!=".#!value" ) {
     117           0 :     error("functions of grid should only output one grid");
     118             :   }
     119         523 :   addValueWithDerivatives( shape );
     120             :   // Set the periodicities of the output components
     121         523 :   myfunc.setPeriodicityForOutputs( this );
     122             :   // Check if we can turn off the derivatives when they are zero
     123         433 :   if( myfunc.getDerivativeZeroIfValueIsZero() )  {
     124         492 :     for(int i=0; i<getNumberOfComponents(); ++i) {
     125         246 :       getPntrToComponent(i)->setDerivativeIsZeroWhenValueIsZero();
     126             :     }
     127             :   }
     128         523 :   setupOnFirstStep( false );
     129        1046 : }
     130             : 
     131             : template <class T>
     132        1045 : void FunctionOfGrid<T>::setupOnFirstStep( const bool incalc ) {
     133             :   double volume = 1.0;
     134        1045 :   const GridCoordinatesObject& mygrid = getGridCoordinatesObject();
     135        1045 :   unsigned npoints = getPntrToArgument(0)->getNumberOfValues();
     136        2090 :   if( mygrid.getGridType()=="flat" ) {
     137         999 :     std::vector<unsigned> shape( getGridCoordinatesObject().getNbin(true) );
     138        1794 :     for(unsigned i=1; i<getNumberOfArguments(); ++i ) {
     139         795 :       if( getPntrToArgument(i)->getRank()==0 ) {
     140         222 :         continue;
     141             :       }
     142         573 :       std::vector<unsigned> s( getPntrToArgument(i)->getShape() );
     143        1160 :       for(unsigned j=0; j<shape.size(); ++j) {
     144         587 :         if( shape[j]!=s[j] ) {
     145           0 :           error("mismatch between sizes of input grids");
     146             :         }
     147             :       }
     148             :     }
     149        1998 :     for(int i=0; i<getNumberOfComponents(); ++i) {
     150         999 :       if( getPntrToComponent(i)->getRank()>0 ) {
     151         833 :         getPntrToComponent(i)->setShape(shape);
     152             :       }
     153             :     }
     154         999 :     std::vector<double> vv( getGridCoordinatesObject().getGridSpacing() );
     155         166 :     volume=vv[0];
     156        1081 :     for(unsigned i=1; i<vv.size(); ++i) {
     157           4 :       volume *=vv[i];
     158             :     }
     159             :   } else {
     160          14 :     volume=4*pi / static_cast<double>( npoints );
     161             :   }
     162             :   // This resizes the scalars
     163        2090 :   for(int i=0; i<getNumberOfComponents(); ++i) {
     164        1045 :     if( getPntrToComponent(i)->getRank()==0 ) {
     165         180 :       getPntrToComponent(i)->resizeDerivatives( npoints );
     166             :     }
     167             :   }
     168        1045 :   if( getName()=="SUM_GRID" ) {
     169             :     volume = 1.0;
     170             :   }
     171             :   // This sets the prefactor to the volume which converts integrals to sums
     172        1045 :   myfunc.setup( this );
     173         180 :   myfunc.setPrefactor( this, volume );
     174        1045 : }
     175             : 
     176             : template <class T>
     177       15460 : const GridCoordinatesObject& FunctionOfGrid<T>::getGridCoordinatesObject() const {
     178       15460 :   ActionWithGrid* ag=ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
     179       15460 :   plumed_assert( ag );
     180       15460 :   return ag->getGridCoordinatesObject();
     181             : }
     182             : 
     183             : template <class T>
     184         198 : std::vector<std::string> FunctionOfGrid<T>::getGridCoordinateNames() const {
     185         198 :   ActionWithGrid* ag=ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
     186         198 :   plumed_assert( ag );
     187         198 :   return ag->getGridCoordinateNames();
     188             : }
     189             : 
     190             : template <class T>
     191        1814 : unsigned FunctionOfGrid<T>::getNumberOfDerivatives() {
     192         514 :   if( myfunc.zeroRank() ) {
     193         514 :     return getPntrToArgument(0)->getNumberOfValues();
     194             :   }
     195        1300 :   unsigned nder = getGridCoordinatesObject().getDimension();
     196        1300 :   return getGridCoordinatesObject().getDimension() + getNumberOfArguments() - myfunc.getArgStart();
     197             : }
     198             : 
     199             : template <class T>
     200      974814 : void FunctionOfGrid<T>::performTask( const unsigned& current, MultiValue& myvals ) const {
     201             :   unsigned argstart=myfunc.getArgStart();
     202      974814 :   std::vector<double> args( getNumberOfArguments() - argstart );
     203     2629182 :   for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
     204     1654368 :     if( getPntrToArgument(i)->getRank()==0 ) {
     205      544232 :       args[i-argstart]=getPntrToArgument(i)->get();
     206             :     } else {
     207     1110136 :       args[i-argstart] = getPntrToArgument(i)->get(current);
     208             :     }
     209             :   }
     210             :   // Calculate the function and its derivatives
     211      974814 :   std::vector<double> vals(1);
     212      974814 :   Matrix<double> derivatives( 1, getNumberOfArguments()-argstart );
     213      974814 :   myfunc.calc( this, args, vals, derivatives );
     214      974814 :   unsigned np = myvals.getTaskIndex();
     215             :   // And set the values and derivatives
     216      974814 :   unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream();
     217      974814 :   myvals.addValue( ostrn, vals[0] );
     218       23665 :   if( !myfunc.zeroRank() ) {
     219             :     // Add the derivatives for a grid
     220     2581852 :     for(unsigned j=argstart; j<getNumberOfArguments(); ++j) {
     221             :       // We store all the derivatives of all the input values - i.e. the grid points these are used in apply
     222     1630703 :       myvals.addDerivative( ostrn, getConstPntrToComponent(0)->getRank()+j-argstart, derivatives(0,j-argstart) );
     223             :       // And now we calculate the derivatives of the value that is stored on the grid correctly so that we can interpolate functions
     224     1630703 :       if( getPntrToArgument(j)->getRank()!=0 ) {
     225     3413677 :         for(unsigned k=0; k<getPntrToArgument(j)->getRank(); ++k) {
     226     2327206 :           myvals.addDerivative( ostrn, k, derivatives(0,j-argstart)*getPntrToArgument(j)->getGridDerivative( np, k ) );
     227             :         }
     228             :       }
     229             :     }
     230      951149 :     unsigned nderivatives = getConstPntrToComponent(0)->getNumberOfGridDerivatives();
     231     4575808 :     for(unsigned j=0; j<nderivatives; ++j) {
     232     3624659 :       myvals.updateIndex( ostrn, j );
     233             :     }
     234       23665 :   } else if( !doNotCalculateDerivatives() ) {
     235             :     // These are the derivatives of the integral
     236        8161 :     myvals.addDerivative( ostrn, current, derivatives(0,0) );
     237        8161 :     myvals.updateIndex( ostrn, current );
     238             :   }
     239      974814 : }
     240             : 
     241             : template <class T>
     242      951149 : void FunctionOfGrid<T>::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
     243             :     const unsigned& bufstart, std::vector<double>& buffer ) const {
     244      951149 :   if( getConstPntrToComponent(0)->getRank()>0 && getConstPntrToComponent(0)->hasDerivatives() ) {
     245             :     plumed_dbg_assert( getNumberOfComponents()==1 && valindex==0 );
     246      951149 :     unsigned nder = getConstPntrToComponent(0)->getNumberOfGridDerivatives();
     247      951149 :     unsigned ostr = getConstPntrToComponent(0)->getPositionInStream();
     248      951149 :     unsigned kp = bufstart + code*(1+nder);
     249      951149 :     buffer[kp] += myvals.get( ostr );
     250     4575808 :     for(unsigned i=0; i<nder; ++i) {
     251     3624659 :       buffer[kp + 1 + i] += myvals.getDerivative( ostr, i );
     252             :     }
     253             :   } else {
     254           0 :     ActionWithVector::gatherStoredValue( valindex, code, myvals, bufstart, buffer );
     255             :   }
     256      951149 : }
     257             : 
     258             : template <class T>
     259        6067 : void FunctionOfGrid<T>::apply() {
     260        6067 :   if( doNotCalculateDerivatives() || !getPntrToComponent(0)->forcesWereAdded() ) {
     261        5442 :     return;
     262             :   }
     263             : 
     264             :   // This applies forces for the integral
     265         494 :   if( myfunc.zeroRank() ) {
     266         494 :     ActionWithVector::apply();
     267         494 :     return;
     268             :   }
     269             : 
     270             :   // Work out how to deal with arguments
     271             :   unsigned nscalars=0, argstart=myfunc.getArgStart();
     272        1870 :   for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
     273        1245 :     if( getPntrToArgument(i)->getRank()==0 ) {
     274          20 :       nscalars++;
     275             :     }
     276             :   }
     277             : 
     278         625 :   std::vector<double> totv(nscalars,0);
     279         625 :   Value* outval=getPntrToComponent(0);
     280        7575 :   for(unsigned i=0; i<outval->getNumberOfValues(); ++i) {
     281             :     nscalars=0;
     282       19845 :     for(unsigned j=argstart; j<getNumberOfArguments(); ++j) {
     283       12895 :       double fforce = outval->getForce(i);
     284       12895 :       if( getPntrToArgument(j)->getRank()==0 ) {
     285        3205 :         totv[nscalars] += fforce*outval->getGridDerivative( i, outval->getRank()+j );
     286        3205 :         nscalars++;
     287             :       } else {
     288        9690 :         double vval = outval->getGridDerivative( i, outval->getRank()+j  );
     289        9690 :         getPntrToArgument(j)->addForce( i, fforce*vval );
     290             :       }
     291             :     }
     292             :   }
     293             :   nscalars=0;
     294        1870 :   for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
     295        1245 :     if( getPntrToArgument(i)->getRank()==0 ) {
     296          20 :       getPntrToArgument(i)->addForce( 0, totv[nscalars] );
     297          20 :       nscalars++;
     298             :     }
     299             :   }
     300             : 
     301             : }
     302             : 
     303             : }
     304             : }
     305             : #endif