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 "ActionWithGrid.h"
      24             : #include "core/PlumedMain.h"
      25             : #include "core/ActionSet.h"
      26             : 
      27             : namespace PLMD {
      28             : namespace gridtools {
      29             : 
      30         184 : void EvaluateGridFunction::registerKeywords( Keywords& keys ) {
      31         184 :   keys.add("compulsory","INTERPOLATION_TYPE","spline","the method to use for interpolation.  Can be spline, linear, ceiling or floor.");
      32         184 :   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");
      33         368 :   keys.setValueDescription("scalar","interpolation of the input grid to get the value of the function at the input arguments");
      34         184 : }
      35             : 
      36         242 : std::vector<bool> EvaluateGridFunction::getPbc() const {
      37         242 :   std::vector<bool> ipbc( gridobject.getDimension() );
      38         484 :   for(unsigned i=0; i<ipbc.size(); ++i) {
      39         242 :     ipbc[i] = gridobject.isPeriodic(i);
      40             :   }
      41         242 :   return ipbc;
      42             : }
      43             : 
      44          88 : void EvaluateGridFunction::read( ActionWithArguments* action ) {
      45          88 :   if( action->getPntrToArgument(0)->getRank()==0 || !action->getPntrToArgument(0)->hasDerivatives() ) {
      46           0 :     action->error("should have one grid as input to this action");
      47             :   }
      48             :   // Get the input grid
      49          88 :   ActionWithGrid* ag = ActionWithGrid::getInputActionWithGrid( (action->getPntrToArgument(0))->getPntrToAction() );
      50         176 :   if( ag->getGridCoordinatesObject().getGridType()!="flat" ) {
      51           0 :     action->error("cannot interpolate on fibonacci sphere");
      52             :   }
      53          88 :   std::vector<bool> ipbc( ag->getGridCoordinatesObject().getDimension() );
      54         184 :   for(unsigned i=0; i<ipbc.size(); ++i) {
      55         192 :     ipbc[i] = ag->getGridCoordinatesObject().isPeriodic(i);
      56             :   }
      57         176 :   gridobject.setup( "flat", ipbc, 0, 0.0 );
      58             :   // Now use this information to create a gridobject
      59             :   std::vector<std::string> argn;
      60          88 :   parseFlag(action,"ZERO_OUTSIDE_GRID_RANGE",set_zero_outside_range);
      61          88 :   if( 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         176 :   parse(action,"INTERPOLATION_TYPE",itype);
      67          88 :   if( itype=="spline" ) {
      68           8 :     interpolation_type=spline;
      69           8 :     spline_interpolator=Tools::make_unique<Interpolator>( action->getPntrToArgument(0), gridobject );
      70          80 :   } else if( itype=="linear" ) {
      71          65 :     interpolation_type=linear;
      72          15 :   } else if( itype=="floor" ) {
      73           0 :     interpolation_type=floor;
      74          15 :   } else if( itype=="ceiling" ) {
      75          15 :     interpolation_type=ceiling;
      76             :   } else {
      77           0 :     action->error("type " + itype + " of interpolation is not defined");
      78             :   }
      79          88 :   action->log.printf("  generating off grid points using %s interpolation \n", itype.c_str() );
      80         176 : }
      81             : 
      82         170 : void EvaluateGridFunction::setup( ActionWithValue* action ) {
      83         170 :   FunctionTemplateBase::setup( action );
      84         170 :   ActionWithArguments* aarg = dynamic_cast<ActionWithArguments*>( action );
      85         170 :   ActionWithGrid* ag = ActionWithGrid::getInputActionWithGrid( (aarg->getPntrToArgument(0))->getPntrToAction() );
      86         170 :   const GridCoordinatesObject & ingrid = ag->getGridCoordinatesObject();
      87         170 :   std::vector<double> sp( ingrid.getGridSpacing() );
      88         170 :   gridobject.setBounds( ingrid.getMin(), ingrid.getMax(), ingrid.getNbin(false), sp );
      89         170 : }
      90             : 
      91       39145 : void EvaluateGridFunction::calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const {
      92       39145 :   if( set_zero_outside_range && !gridobject.inbounds( args ) ) {
      93           0 :     vals[0]=0.0;
      94           0 :     return;
      95             :   }
      96             :   unsigned dimension = gridobject.getDimension();
      97             :   plumed_dbg_assert( args.size()==dimension && vals.size()==1 );
      98       39145 :   if( interpolation_type==spline ) {
      99       23729 :     std::vector<double> der( dimension );
     100       23729 :     vals[0] =  spline_interpolator->splineInterpolation( args, der );
     101       74416 :     for(unsigned j=0; j<dimension; ++j) {
     102       50687 :       derivatives(0,j) = der[j];
     103             :     }
     104       15416 :   } else if( interpolation_type==linear ) {
     105             :     Value* values=action->getPntrToArgument(0);
     106        6600 :     std::vector<double> xfloor(dimension);
     107        6600 :     std::vector<unsigned> indices(dimension), nindices(dimension), ind(dimension);
     108        6600 :     gridobject.getIndices( args, indices );
     109        6600 :     unsigned nn=gridobject.getIndex(args);
     110        6600 :     gridobject.getGridPointCoordinates( nn, nindices, xfloor );
     111        6600 :     double y1 = values->get(nn);
     112        6600 :     vals[0] = y1;
     113       13200 :     for(unsigned i=0; i<args.size(); ++i) {
     114             :       int x0=1;
     115        6600 :       if(nindices[i]==indices[i]) {
     116             :         x0=0;
     117             :       }
     118        6600 :       double ddx=gridobject.getGridSpacing()[i];
     119        6600 :       double X = fabs((args[i]-xfloor[i])/ddx-(double)x0);
     120       13200 :       for(unsigned j=0; j<args.size(); ++j) {
     121        6600 :         ind[j] = indices[j];
     122             :       }
     123        6600 :       if( gridobject.isPeriodic(i) && (ind[i]+1)==gridobject.getNbin(false)[i] ) {
     124           0 :         ind[i]=0;
     125             :       } else {
     126        6600 :         ind[i] = ind[i] + 1;
     127             :       }
     128        6600 :       vals[0] += ( values->get( gridobject.getIndex(ind) ) - y1 )*X;
     129        6600 :       derivatives(0,i) = ( values->get( gridobject.getIndex(ind) ) - y1 ) / ddx;
     130             :     }
     131        8816 :   } else if( interpolation_type==floor ) {
     132             :     Value* values=action->getPntrToArgument(0);
     133           0 :     std::vector<unsigned> indices(dimension);
     134           0 :     gridobject.getIndices( args, indices );
     135           0 :     unsigned nn = gridobject.getIndex(indices);
     136             :     plumed_dbg_assert( nn<values->getNumberOfValues() );
     137           0 :     vals[0] = values->get( nn );
     138           0 :     for(unsigned j=0; j<dimension; ++j) {
     139           0 :       derivatives(0,j) = values->getGridDerivative( nn, j );
     140             :     }
     141        8816 :   } else if( interpolation_type==ceiling ) {
     142             :     Value* values=action->getPntrToArgument(0);
     143        8816 :     std::vector<unsigned> indices(dimension);
     144        8816 :     gridobject.getIndices( args, indices );
     145       17632 :     for(unsigned i=0; i<indices.size(); ++i) {
     146       17632 :       if( gridobject.isPeriodic(i) && (indices[i]+1)==gridobject.getNbin(false)[i] ) {
     147        2204 :         indices[i]=0;
     148             :       } else {
     149        6612 :         indices[i] = indices[i] + 1;
     150             :       }
     151             :     }
     152        8816 :     unsigned nn = gridobject.getIndex(indices);
     153        8816 :     vals[0] = values->get( nn );
     154       17632 :     for(unsigned j=0; j<dimension; ++j) {
     155        8816 :       derivatives(0,j) = values->getGridDerivative( nn, j );
     156             :     }
     157             :   } else {
     158           0 :     plumed_error();
     159             :   }
     160             : }
     161             : 
     162        1956 : void EvaluateGridFunction::applyForce( const ActionWithArguments* action, const std::vector<double>& args, const double& force, std::vector<double>& forcesToApply ) const {
     163             :   unsigned dimension = gridobject.getDimension();
     164        1956 :   if( interpolation_type==spline ) {
     165           0 :     action->error("can't apply forces on values interpolated using splines");
     166        1956 :   } else if( interpolation_type==linear ) {
     167             :     Value* values=action->getPntrToArgument(0);
     168         100 :     std::vector<double> xfloor(dimension);
     169         100 :     std::vector<unsigned> indices(dimension), nindices(dimension), ind(dimension);
     170         100 :     gridobject.getIndices( args, indices );
     171         100 :     unsigned nn=gridobject.getIndex(args);
     172         100 :     gridobject.getGridPointCoordinates( nn, nindices, xfloor );
     173         200 :     for(unsigned i=0; i<args.size(); ++i) {
     174             :       int x0=1;
     175         100 :       if(nindices[i]==indices[i]) {
     176             :         x0=0;
     177             :       }
     178         100 :       double ddx=gridobject.getGridSpacing()[i];
     179         100 :       double X = fabs((args[i]-xfloor[i])/ddx-(double)x0);
     180         200 :       for(unsigned j=0; j<args.size(); ++j) {
     181         100 :         ind[j] = indices[j];
     182             :       }
     183         100 :       if( gridobject.isPeriodic(i) && (ind[i]+1)==gridobject.getNbin(false)[i] ) {
     184           0 :         ind[i]=0;
     185             :       } else {
     186         100 :         ind[i] = ind[i] + 1;
     187             :       }
     188         100 :       forcesToApply[nn] += force*(1-X);
     189         100 :       forcesToApply[gridobject.getIndex(ind)] += X*force;
     190             :     }
     191        1856 :   } else if( interpolation_type==floor ) {
     192             :     Value* values=action->getPntrToArgument(0);
     193           0 :     std::vector<unsigned> indices(dimension);
     194           0 :     gridobject.getIndices( args, indices );
     195           0 :     unsigned nn = gridobject.getIndex(indices);
     196           0 :     forcesToApply[nn] += force;
     197        1856 :   } else if( interpolation_type==ceiling ) {
     198             :     Value* values=action->getPntrToArgument(0);
     199        1856 :     std::vector<unsigned> indices(dimension);
     200        1856 :     gridobject.getIndices( args, indices );
     201        3712 :     for(unsigned i=0; i<indices.size(); ++i) {
     202        3712 :       if( gridobject.isPeriodic(i) && (indices[i]+1)==gridobject.getNbin(false)[i] ) {
     203         464 :         indices[i]=0;
     204             :       } else {
     205        1392 :         indices[i] = indices[i] + 1;
     206             :       }
     207             :     }
     208        1856 :     unsigned nn = gridobject.getIndex(indices);
     209        1856 :     forcesToApply[nn] += force;
     210             :   } else {
     211           0 :     plumed_error();
     212             :   }
     213        1956 : }
     214             : 
     215             : }
     216             : }