Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-2017 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 "FunctionShortcut.h"
      23             : #include "FunctionOfScalar.h"
      24             : #include "FunctionOfVector.h"
      25             : #include "core/ActionRegister.h"
      26             : #include "FunctionTemplateBase.h"
      27             : 
      28             : #include <cmath>
      29             : 
      30             : namespace PLMD {
      31             : namespace function {
      32             : 
      33             : //+PLUMEDOC FUNCTION MOMENTS
      34             : /*
      35             : Calculate the moments of the distribution of input quantities
      36             : 
      37             : \par Examples
      38             : 
      39             : */
      40             : //+ENDPLUMEDOC
      41             : 
      42             : //+PLUMEDOC FUNCTION MOMENTS_SCALAR
      43             : /*
      44             : Calculate the moments of the distribution of input quantities
      45             : 
      46             : \par Examples
      47             : 
      48             : */
      49             : //+ENDPLUMEDOC
      50             : 
      51             : //+PLUMEDOC FUNCTION MOMENTS_VECTOR
      52             : /*
      53             : Calculate the moments of the distribution of input vectors
      54             : 
      55             : \par Examples
      56             : 
      57             : */
      58             : //+ENDPLUMEDOC
      59             : 
      60          69 : class Moments : public FunctionTemplateBase {
      61             :   bool isperiodic, scalar_out;
      62             :   double min, max, pfactor;
      63             :   std::vector<int> powers;
      64             : public:
      65             :   void registerKeywords(Keywords& keys) override;
      66             :   void read( ActionWithArguments* action ) override;
      67           7 :   bool zeroRank() const override { return scalar_out; }
      68         221 :   bool doWithTasks() const override { return !scalar_out; }
      69             :   std::vector<std::string> getComponentsPerLabel() const override ;
      70             :   void setPeriodicityForOutputs( ActionWithValue* action ) override;
      71             :   void calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const override;
      72             : };
      73             : 
      74             : typedef FunctionShortcut<Moments> MomentsShortcut;
      75             : PLUMED_REGISTER_ACTION(MomentsShortcut,"MOMENTS")
      76             : typedef FunctionOfScalar<Moments> ScalarMoments;
      77             : PLUMED_REGISTER_ACTION(ScalarMoments,"MOMENTS_SCALAR")
      78             : typedef FunctionOfVector<Moments> VectorMoments;
      79             : PLUMED_REGISTER_ACTION(VectorMoments,"MOMENTS_VECTOR")
      80             : 
      81          27 : void Moments::registerKeywords(Keywords& keys) {
      82          54 :   keys.add("compulsory","POWERS","calculate the central moments of the distribution of collective variables. "
      83             :            "The \\f$m\\f$th central moment of a distribution is calculated using \\f$\\frac{1}{N} \\sum_{i=1}^N ( s_i - \\overline{s} )^m \\f$, where \\f$\\overline{s}\\f$ is "
      84             :            "the average for the distribution. The POWERS keyword takes a lists of integers as input or a range. Each integer is a value of \\f$m\\f$. The final "
      85             :            "calculated values can be referenced using moment-\\f$m\\f$.");
      86          54 :   keys.addOutputComponent("moment","default","scalar","the central moments of the distribution of values. The second central moment "
      87             :                           "would be referenced elsewhere in the input file using "
      88             :                           "<em>label</em>.moment-2, the third as <em>label</em>.moment-3, etc.");
      89          27 : }
      90             : 
      91           5 : void Moments::read( ActionWithArguments* action ) {
      92           5 :   scalar_out = action->getNumberOfArguments()==1;
      93           5 :   if( scalar_out && action->getPntrToArgument(0)->getRank()==0 ) action->error("cannot calculate moments if only given one variable");
      94             : 
      95           5 :   isperiodic = (action->getPntrToArgument(0))->isPeriodic();
      96           5 :   if( isperiodic ) {
      97           0 :     std::string str_min, str_max; (action->getPntrToArgument(0))->getDomain( str_min, str_max );
      98           0 :     for(unsigned i=1; i<action->getNumberOfArguments(); ++i) {
      99           0 :       if( !(action->getPntrToArgument(i))->isPeriodic() ) action->error("cannot mix periodic and non periodic variables when calculating moments");
     100           0 :       std::string str_min2, str_max2; (action->getPntrToArgument(i))->getDomain( str_min2, str_max2);
     101           0 :       if( str_min!=str_min2 || str_max!=str_max2 ) action->error("all input arguments should have same domain when calculating moments");
     102             :     }
     103           0 :     Tools::convert(str_min,min); Tools::convert(str_max,max); pfactor = 2*pi / ( max-min );
     104             :   } else {
     105          14 :     for(unsigned i=1; i<action->getNumberOfArguments(); ++i) {
     106           9 :       if( (action->getPntrToArgument(i))->isPeriodic() ) action->error("cannot mix periodic and non periodic variables when calculating moments");
     107             :     }
     108             :   }
     109             : 
     110           5 :   parseVector(action,"POWERS",powers);
     111          13 :   for(unsigned i=0; i<powers.size(); ++i) {
     112           8 :     if( powers[i]<2 ) action->error("first central moment is zero do you really need to calculate that");
     113           8 :     action->log.printf("  computing %dth central moment of distribution of input cvs \n", powers[i]);
     114             :   }
     115           5 : }
     116             : 
     117           5 : std::vector<std::string> Moments::getComponentsPerLabel() const {
     118             :   std::vector<std::string> comp; std::string num;
     119          13 :   for(unsigned i=0; i<powers.size(); ++i) {
     120          16 :     Tools::convert(powers[i],num); comp.push_back( "-" + num );
     121             :   }
     122           5 :   return comp;
     123           0 : }
     124             : 
     125           5 : void Moments::setPeriodicityForOutputs( ActionWithValue* action ) {
     126          13 :   for(unsigned i=0; i<powers.size(); ++i) { std::string num; Tools::convert(powers[i],num); action->componentIsNotPeriodic("moment-" + num); }
     127           5 : }
     128             : 
     129         222 : void Moments::calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const {
     130         222 :   double mean=0; double inorm = 1.0 / static_cast<double>( args.size() );
     131         222 :   if( isperiodic ) {
     132             :     double sinsum=0, cossum=0, val;
     133           0 :     for(unsigned i=0; i<args.size(); ++i) { val=pfactor*( args[i] - min ); sinsum+=sin(val); cossum+=cos(val); }
     134           0 :     mean = 0.5 + atan2( inorm*sinsum, inorm*cossum ) / (2*pi);
     135           0 :     mean = min + (max-min)*mean;
     136             :   } else {
     137        1758 :     for(unsigned i=0; i<args.size(); ++i) mean+=args[i];
     138         222 :     mean *= inorm;
     139             :   }
     140             : 
     141             :   Value* arg0 = action->getPntrToArgument(0);
     142         656 :   for(unsigned npow=0; npow<powers.size(); ++npow) {
     143             :     double dev1=0;
     144        3406 :     for(unsigned i=0; i<args.size(); ++i) dev1+=pow( arg0->difference( mean, args[i] ), powers[npow] - 1 );
     145         434 :     dev1*=inorm; vals[npow] = 0; double prefactor = powers[npow]*inorm;
     146        3406 :     for(unsigned i=0; i<args.size(); ++i) {
     147        2972 :       double tmp=arg0->difference( mean, args[i] ); vals[npow] += inorm*pow( tmp, powers[npow] );
     148        2972 :       derivatives(npow,i) = prefactor*(pow( tmp, powers[npow] - 1 ) - dev1);
     149             :     }
     150             :   }
     151         222 : }
     152             : 
     153             : }
     154             : }
     155             : 
     156             :