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          66 : 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           0 :   bool zeroRank() const override {
      68           7 :     return scalar_out;
      69             :   }
      70           0 :   bool doWithTasks() const override {
      71         221 :     return !scalar_out;
      72             :   }
      73             :   std::vector<std::string> getComponentsPerLabel() const override ;
      74             :   void setPeriodicityForOutputs( ActionWithValue* action ) override;
      75             :   void calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const override;
      76             : };
      77             : 
      78             : typedef FunctionShortcut<Moments> MomentsShortcut;
      79             : PLUMED_REGISTER_ACTION(MomentsShortcut,"MOMENTS")
      80             : typedef FunctionOfScalar<Moments> ScalarMoments;
      81             : PLUMED_REGISTER_ACTION(ScalarMoments,"MOMENTS_SCALAR")
      82             : typedef FunctionOfVector<Moments> VectorMoments;
      83             : PLUMED_REGISTER_ACTION(VectorMoments,"MOMENTS_VECTOR")
      84             : 
      85          33 : void Moments::registerKeywords(Keywords& keys) {
      86          66 :   keys.add("compulsory","POWERS","calculate the central moments of the distribution of collective variables. "
      87             :            "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 "
      88             :            "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 "
      89             :            "calculated values can be referenced using moment-\\f$m\\f$.");
      90          66 :   keys.addOutputComponent("moment","default","the central moments of the distribution of values. The second central moment "
      91             :                           "would be referenced elsewhere in the input file using "
      92             :                           "<em>label</em>.moment-2, the third as <em>label</em>.moment-3, etc.");
      93          33 : }
      94             : 
      95           5 : void Moments::read( ActionWithArguments* action ) {
      96           5 :   scalar_out = action->getNumberOfArguments()==1;
      97           5 :   if( scalar_out && action->getPntrToArgument(0)->getRank()==0 ) {
      98           0 :     action->error("cannot calculate moments if only given one variable");
      99             :   }
     100             : 
     101           5 :   isperiodic = (action->getPntrToArgument(0))->isPeriodic();
     102           5 :   if( isperiodic ) {
     103             :     std::string str_min, str_max;
     104           0 :     (action->getPntrToArgument(0))->getDomain( str_min, str_max );
     105           0 :     for(unsigned i=1; i<action->getNumberOfArguments(); ++i) {
     106           0 :       if( !(action->getPntrToArgument(i))->isPeriodic() ) {
     107           0 :         action->error("cannot mix periodic and non periodic variables when calculating moments");
     108             :       }
     109             :       std::string str_min2, str_max2;
     110           0 :       (action->getPntrToArgument(i))->getDomain( str_min2, str_max2);
     111           0 :       if( str_min!=str_min2 || str_max!=str_max2 ) {
     112           0 :         action->error("all input arguments should have same domain when calculating moments");
     113             :       }
     114             :     }
     115           0 :     Tools::convert(str_min,min);
     116           0 :     Tools::convert(str_max,max);
     117           0 :     pfactor = 2*pi / ( max-min );
     118             :   } else {
     119          14 :     for(unsigned i=1; i<action->getNumberOfArguments(); ++i) {
     120           9 :       if( (action->getPntrToArgument(i))->isPeriodic() ) {
     121           0 :         action->error("cannot mix periodic and non periodic variables when calculating moments");
     122             :       }
     123             :     }
     124             :   }
     125             : 
     126           5 :   parseVector(action,"POWERS",powers);
     127          13 :   for(unsigned i=0; i<powers.size(); ++i) {
     128           8 :     if( powers[i]<2 ) {
     129           0 :       action->error("first central moment is zero do you really need to calculate that");
     130             :     }
     131           8 :     action->log.printf("  computing %dth central moment of distribution of input cvs \n", powers[i]);
     132             :   }
     133           5 : }
     134             : 
     135           5 : std::vector<std::string> Moments::getComponentsPerLabel() const {
     136             :   std::vector<std::string> comp;
     137             :   std::string num;
     138          13 :   for(unsigned i=0; i<powers.size(); ++i) {
     139           8 :     Tools::convert(powers[i],num);
     140          16 :     comp.push_back( "-" + num );
     141             :   }
     142           5 :   return comp;
     143           0 : }
     144             : 
     145           5 : void Moments::setPeriodicityForOutputs( ActionWithValue* action ) {
     146          13 :   for(unsigned i=0; i<powers.size(); ++i) {
     147             :     std::string num;
     148           8 :     Tools::convert(powers[i],num);
     149          16 :     action->componentIsNotPeriodic("moment-" + num);
     150             :   }
     151           5 : }
     152             : 
     153         222 : void Moments::calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const {
     154             :   double mean=0;
     155         222 :   double inorm = 1.0 / static_cast<double>( args.size() );
     156         222 :   if( isperiodic ) {
     157             :     double sinsum=0, cossum=0, val;
     158           0 :     for(unsigned i=0; i<args.size(); ++i) {
     159           0 :       val=pfactor*( args[i] - min );
     160           0 :       sinsum+=sin(val);
     161           0 :       cossum+=cos(val);
     162             :     }
     163           0 :     mean = 0.5 + atan2( inorm*sinsum, inorm*cossum ) / (2*pi);
     164           0 :     mean = min + (max-min)*mean;
     165             :   } else {
     166        1758 :     for(unsigned i=0; i<args.size(); ++i) {
     167        1536 :       mean+=args[i];
     168             :     }
     169         222 :     mean *= inorm;
     170             :   }
     171             : 
     172             :   Value* arg0 = action->getPntrToArgument(0);
     173         656 :   for(unsigned npow=0; npow<powers.size(); ++npow) {
     174             :     double dev1=0;
     175        3406 :     for(unsigned i=0; i<args.size(); ++i) {
     176        2972 :       dev1+=pow( arg0->difference( mean, args[i] ), powers[npow] - 1 );
     177             :     }
     178         434 :     dev1*=inorm;
     179         434 :     vals[npow] = 0;
     180         434 :     double prefactor = powers[npow]*inorm;
     181        3406 :     for(unsigned i=0; i<args.size(); ++i) {
     182        2972 :       double tmp=arg0->difference( mean, args[i] );
     183        2972 :       vals[npow] += inorm*pow( tmp, powers[npow] );
     184        2972 :       derivatives(npow,i) = prefactor*(pow( tmp, powers[npow] - 1 ) - dev1);
     185             :     }
     186             :   }
     187         222 : }
     188             : 
     189             : }
     190             : }
     191             : 
     192             :