Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-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 "core/ActionRegister.h"
      23             : #include "ActionWithGrid.h"
      24             : 
      25             : //+PLUMEDOC GRIDANALYSIS GET_VOLUME_ELEMENT
      26             : /*
      27             : Get the volume element from the input grid
      28             : 
      29             : This action is used in the shortcut action that is used within PLUMED to calculate integrals as you can see below
      30             : 
      31             : ```plumed
      32             : c1: COORDINATIONNUMBER SPECIESA=1-10 SPECIESB=1-200 SWITCH={RATIONAL R_0=1.0}
      33             : hu: KDE ARG=c1 GRID_BIN=200 GRID_MIN=0 GRID_MAX=10 BANDWIDTH=0.1
      34             : h: CUSTOM ARG=hu FUNC=x/10 PERIODIC=NO
      35             : iv: INTEGRATE_GRID ARG=h PERIODIC=NO
      36             : ```
      37             : 
      38             : As discused in the documentation for [INTEGRATE_GRID](INTEGRATE_GRID.md) we compute (approximate) integrals in PLUMED using:
      39             : 
      40             : $$
      41             : y = \frac{b-a}{n} \sum_{i=0}^{n} f\left( a + \frac{i(b-a)}{n} \right)
      42             : $$
      43             : 
      44             : This action gets the value of $\frac{b-a}{n}$ so that we can multiply by this value by the sum in the expression
      45             : above, which is computed by a [SUM](SUM.md) action.
      46             : 
      47             : */
      48             : //+ENDPLUMEDOC
      49             : 
      50             : namespace PLMD {
      51             : namespace gridtools {
      52             : 
      53             : class GetGridVolumeElement :
      54             :   public ActionWithValue,
      55             :   public ActionWithArguments {
      56             : public:
      57             :   static void registerKeywords( Keywords& keys );
      58             :   explicit GetGridVolumeElement(const ActionOptions&ao);
      59             :   unsigned getNumberOfDerivatives() override ;
      60             :   void calculate() override ;
      61          60 :   void apply() override {}
      62             : };
      63             : 
      64             : PLUMED_REGISTER_ACTION(GetGridVolumeElement,"GET_VOLUME_ELEMENT")
      65             : 
      66          36 : void GetGridVolumeElement::registerKeywords( Keywords& keys ) {
      67          36 :   Action::registerKeywords( keys );
      68          36 :   ActionWithValue::registerKeywords( keys );
      69          36 :   ActionWithArguments::registerKeywords( keys );
      70          72 :   keys.addInputKeyword("compulsory","ARG","grid","the label for function on the grid that you need the volume element for");
      71          72 :   keys.setValueDescription("scalar","the volume element ");
      72          36 :   keys.remove("NUMERICAL_DERIVATIVES");
      73          36 : }
      74             : 
      75          16 : GetGridVolumeElement::GetGridVolumeElement(const ActionOptions&ao):
      76             :   Action(ao),
      77             :   ActionWithValue(ao),
      78          16 :   ActionWithArguments(ao) {
      79          16 :   if( getNumberOfArguments()!=1 ) {
      80           0 :     error("should only be one argument to this action");
      81             :   }
      82          16 :   if( getPntrToArgument(0)->getRank()==0 || !getPntrToArgument(0)->hasDerivatives() ) {
      83           0 :     error("input to this action should be a grid");
      84             :   }
      85          16 :   ActionWithGrid* ag = dynamic_cast<ActionWithGrid*>( getPntrToArgument(0)->getPntrToAction() );
      86          16 :   if( !ag ) {
      87           0 :     error("could not find grid details from input value");
      88             :   }
      89          16 :   addValue();
      90          16 :   setNotPeriodic();
      91          16 : }
      92             : 
      93          18 : unsigned GetGridVolumeElement::getNumberOfDerivatives() {
      94          18 :   return 0;
      95             : }
      96             : 
      97          60 : void GetGridVolumeElement::calculate() {
      98             :   double volume;
      99          60 :   ActionWithGrid* ag = dynamic_cast<ActionWithGrid*>( getPntrToArgument(0)->getPntrToAction() );
     100          60 :   plumed_assert( ag );
     101          60 :   const GridCoordinatesObject& mygrid = ag->getGridCoordinatesObject();
     102          60 :   unsigned npoints = getPntrToArgument(0)->getNumberOfValues();
     103          60 :   if( mygrid.getGridType()=="flat" ) {
     104          40 :     std::vector<double> vv( mygrid.getGridSpacing() );
     105          40 :     volume=vv[0];
     106          60 :     for(unsigned i=1; i<vv.size(); ++i) {
     107          20 :       volume *=vv[i];
     108             :     }
     109             :   } else {
     110          20 :     volume = 4*pi / static_cast<double>( npoints );
     111             :   }
     112          60 :   getPntrToComponent(0)->set( volume );
     113          60 : }
     114             : 
     115             : }
     116             : }