Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-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 "EvaluateGridFunction.h"
      23             : #include "core/PlumedMain.h"
      24             : #include "core/ActionSet.h"
      25             : 
      26             : namespace PLMD {
      27             : namespace gridtools {
      28             : 
      29          91 : void EvaluateGridFunction::registerKeywords( Keywords& keys ) {
      30          91 :   keys.add("compulsory","INTERPOLATION_TYPE","spline","the method to use for interpolation.  Can be spline, linear, ceiling or floor.");
      31          91 :   keys.addFlag("ZERO_OUTSIDE_GRID_RANGE",false,"if we are asked to evaluate the function for a number that is outside the range of the grid set it to zero");
      32         182 :   keys.setValueDescription("scalar","interpolation of the input grid to get the value of the function at the input arguments");
      33          91 : }
      34             : 
      35          23 : std::vector<bool> EvaluateGridFunction::getPbc() const {
      36             :   const GridCoordinatesObject & gridobject = getGridObject();
      37          23 :   std::vector<bool> ipbc( gridobject.getDimension() );
      38          47 :   for(unsigned i=0; i<ipbc.size(); ++i) {
      39          24 :     ipbc[i] = gridobject.isPeriodic(i);
      40             :   }
      41          23 :   return ipbc;
      42             : }
      43             : 
      44          33 : void EvaluateGridFunction::read( EvaluateGridFunction& func, ActionWithArguments* action, function::FunctionOptions& options ) {
      45          33 :   func.function = action->getPntrToArgument(0);
      46          33 :   if( func.function->getRank()==0 || !func.function->hasDerivatives() ) {
      47           0 :     action->error("should have one grid as input to this action");
      48             :   }
      49             :   // Get the input grid
      50          33 :   func.gridact = ActionWithGrid::getInputActionWithGrid( func.function->getPntrToAction() );
      51          33 :   plumed_assert( func.gridact );
      52          33 :   if( func.gridact->getGridCoordinatesObject().getGridType()!="flat" ) {
      53           0 :     action->error("cannot interpolate on fibonacci sphere");
      54             :   }
      55          33 :   std::vector<std::string> argn;
      56          33 :   if( action->keywords.exists("ZERO_OUTSIDE_GRID_RANGE") ) {
      57          12 :     action->parseFlag("ZERO_OUTSIDE_GRID_RANGE",func.set_zero_outside_range);
      58             :   } else {
      59          27 :     func.set_zero_outside_range=false;
      60             :   }
      61          33 :   if( func.set_zero_outside_range ) {
      62           0 :     action->log.printf("  function is zero outside grid range \n");
      63             :   }
      64             :   // Get the type of interpolation that we are doing
      65             :   std::string itype;
      66          33 :   action->parse("INTERPOLATION_TYPE",itype);
      67          33 :   if( itype=="spline" ) {
      68          12 :     func.interpolation_type=spline;
      69          24 :     func.spline_interpolator=Tools::make_unique<Interpolator>( func.function, func.gridact->getGridCoordinatesObject() );
      70          21 :   } else if( itype=="linear" ) {
      71           5 :     func.interpolation_type=linear;
      72          16 :   } else if( itype=="floor" ) {
      73           1 :     func.interpolation_type=floor;
      74          15 :   } else if( itype=="ceiling" ) {
      75          15 :     func.interpolation_type=ceiling;
      76             :   } else {
      77           0 :     action->error("type " + itype + " of interpolation is not defined");
      78             :   }
      79          33 :   action->log.printf("  generating off grid points using %s interpolation \n", itype.c_str() );
      80          33 : }
      81             : 
      82       47872 : void EvaluateGridFunction::calc( const EvaluateGridFunction& func, bool noderiv, const View<const double> args, function::FunctionOutput& funcout ) {
      83             :   const GridCoordinatesObject & gridobject = func.getGridObject();
      84       47872 :   if( func.set_zero_outside_range && !gridobject.inbounds( args ) ) {
      85           0 :     funcout.values[0]=0.0;
      86           0 :     return;
      87             :   }
      88             :   unsigned dimension = gridobject.getDimension();
      89             :   plumed_dbg_assert( args.size()==dimension && funcout.values.size()==1 );
      90       47872 :   if( func.interpolation_type==spline ) {
      91       32155 :     std::vector<double> der( dimension );
      92       32155 :     funcout.values[0] =  func.spline_interpolator->splineInterpolation( args, der );
      93       32155 :     if( noderiv ) {
      94             :       return ;
      95             :     }
      96       74485 :     for(unsigned j=0; j<dimension; ++j) {
      97       50735 :       funcout.derivs[0][j] = der[j];
      98             :     }
      99       15717 :   } else if( func.interpolation_type==linear ) {
     100        6840 :     std::vector<double> xfloor(dimension);
     101        6840 :     std::vector<unsigned> indices(dimension), nindices(dimension), ind(dimension);
     102        6840 :     gridobject.getIndices( args, indices );
     103        6840 :     unsigned nn=gridobject.getIndex(args);
     104        6840 :     gridobject.getGridPointCoordinates( nn, nindices, xfloor );
     105        6840 :     double y1 = func.function->get(nn);
     106        6840 :     funcout.values[0] = y1;
     107       20080 :     for(unsigned i=0; i<args.size(); ++i) {
     108             :       int x0=1;
     109       13240 :       if(nindices[i]==indices[i]) {
     110             :         x0=0;
     111             :       }
     112       13240 :       double ddx=gridobject.getGridSpacing()[i];
     113       13240 :       double X = fabs((args[i]-xfloor[i])/ddx-(double)x0);
     114       39280 :       for(unsigned j=0; j<args.size(); ++j) {
     115       26040 :         ind[j] = indices[j];
     116             :       }
     117       13240 :       if( gridobject.isPeriodic(i) && (ind[i]+1)==gridobject.getNbin(false)[i] ) {
     118           0 :         ind[i]=0;
     119             :       } else {
     120       13240 :         ind[i] = ind[i] + 1;
     121             :       }
     122       13240 :       if( X<epsilon ) {
     123        6400 :         continue ;
     124             :       }
     125        6840 :       funcout.values[0] += ( func.function->get( gridobject.getIndex(ind) ) - y1 )*X;
     126        6840 :       if( noderiv ) {
     127           0 :         continue ;
     128             :       }
     129        6840 :       funcout.derivs[0][i] = ( func.function->get( gridobject.getIndex(ind) ) - y1 ) / ddx;
     130             :     }
     131        8877 :   } else if( func.interpolation_type==floor ) {
     132          61 :     std::vector<unsigned> indices(dimension);
     133          61 :     gridobject.getIndices( args, indices );
     134          61 :     unsigned nn = gridobject.getIndex(indices);
     135             :     plumed_dbg_assert( nn<func.function->getNumberOfValues() );
     136          61 :     funcout.values[0] = func.function->get( nn );
     137          61 :     if( noderiv ) {
     138             :       return;
     139             :     }
     140         122 :     for(unsigned j=0; j<dimension; ++j) {
     141          61 :       funcout.derivs[0][j] = func.function->getGridDerivative( nn, j );
     142             :     }
     143        8816 :   } else if( func.interpolation_type==ceiling ) {
     144        8816 :     std::vector<unsigned> indices(dimension);
     145        8816 :     gridobject.getIndices( args, indices );
     146       17632 :     for(unsigned i=0; i<indices.size(); ++i) {
     147       11020 :       if( gridobject.isPeriodic(i) && (indices[i]+1)==gridobject.getNbin(false)[i] ) {
     148        2204 :         indices[i]=0;
     149             :       } else {
     150        6612 :         indices[i] = indices[i] + 1;
     151             :       }
     152             :     }
     153        8816 :     unsigned nn = gridobject.getIndex(indices);
     154        8816 :     funcout.values[0] = func.function->get( nn );
     155        8816 :     if( noderiv ) {
     156             :       return;
     157             :     }
     158       17632 :     for(unsigned j=0; j<dimension; ++j) {
     159        8816 :       funcout.derivs[0][j] = func.function->getGridDerivative( nn, j );
     160             :     }
     161             :   } else {
     162           0 :     plumed_error();
     163             :   }
     164             : }
     165             : 
     166             : }
     167             : }