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/ActionShortcut.h"
      23             : #include "core/ActionWithArguments.h"
      24             : #include "core/PlumedMain.h"
      25             : #include "core/ActionRegister.h"
      26             : 
      27             : //+PLUMEDOC GRIDCALC AVERAGE
      28             : /*
      29             : Calculate the ensemble average of a collective variable
      30             : 
      31             : The ensemble average for a non-periodic, collective variable, \f$s\f$ is given by the following expression:
      32             : 
      33             : \f[
      34             : \langle s \rangle = \frac{ \sum_{t'=0}^t w(t') s(t') }{ \sum_{t'=0}^t w(t') }
      35             : \f]
      36             : 
      37             : Here the sum runs over a the trajectory and \f$s(t')\f$ is used to denote the value of the collective variable
      38             : at time \f$t'\f$.  The final quantity evaluated is a weighted
      39             : average as the weights, \f$w(t')\f$, allow us to negate the effect any bias might have on the region of phase space
      40             : sampled by the system.  This is discussed in the section of the manual on \ref Analysis.
      41             : 
      42             : When the variable is periodic (e.g. \ref TORSION) and has a value, \f$s\f$, in \f$a \le s \le b\f$ the ensemble average is evaluated using:
      43             : 
      44             : \f[
      45             : \langle s \rangle = a + \frac{b - a}{2\pi} \arctan \left[ \frac{ \sum_{t'=0}^t w(t') \sin\left( \frac{2\pi [s(t')-a]}{b - a} \right) }{ \sum_{t'=0}^t w(t') \cos\left( \frac{2\pi [s(t')-a]}{b - a} \right) } \right]
      46             : \f]
      47             : 
      48             : \par Examples
      49             : 
      50             : The following example calculates the ensemble average for the distance between atoms 1 and 2
      51             : and output this to a file called COLVAR.  In this example it is assumed that no bias is acting
      52             : on the system and that the weights, \f$w(t')\f$ in the formulas above can thus all be set equal
      53             : to one.
      54             : 
      55             : \plumedfile
      56             : d1: DISTANCE ATOMS=1,2
      57             : d1a: AVERAGE ARG=d1
      58             : PRINT ARG=d1a FILE=colvar STRIDE=100
      59             : \endplumedfile
      60             : 
      61             : The following example calculates the ensemble average for the torsional angle involving atoms 1, 2, 3 and 4.
      62             : At variance with the previous example this quantity is periodic so the second formula in the above introduction
      63             : is used to calculate the average.  Furthermore, by using the CLEAR keyword we have specified that block averages
      64             : are to be calculated.  Consequently, after 100 steps all the information acquired thus far in the simulation is
      65             : forgotten and the process of averaging is begun again.  The quantities output in the colvar file are thus the
      66             : block averages taken over the first 100 frames of the trajectory, the block average over the second 100 frames
      67             : of trajectory and so on.
      68             : 
      69             : \plumedfile
      70             : t1: TORSION ATOMS=1,2,3,4
      71             : t1a: AVERAGE ARG=t1 CLEAR=100
      72             : PRINT ARG=t1a FILE=colvar STRIDE=100
      73             : \endplumedfile
      74             : 
      75             : This third example incorporates a bias.  Notice that the effect the bias has on the ensemble average is removed by taking
      76             : advantage of the \ref REWEIGHT_BIAS method.  The final ensemble averages output to the file are thus block ensemble averages for the
      77             : unbiased canonical ensemble at a temperature of 300 K.
      78             : 
      79             : \plumedfile
      80             : t1: TORSION ATOMS=1,2,3,4
      81             : RESTRAINT ARG=t1 AT=pi KAPPA=100.
      82             : ww: REWEIGHT_BIAS TEMP=300
      83             : t1a: AVERAGE ARG=t1 LOGWEIGHTS=ww CLEAR=100
      84             : PRINT ARG=t1a FILE=colvar STRIDE=100
      85             : \endplumedfile
      86             : 
      87             : */
      88             : //+ENDPLUMEDOC
      89             : 
      90             : namespace PLMD {
      91             : namespace generic {
      92             : 
      93             : class Average : public ActionShortcut {
      94             : public:
      95             :   static void registerKeywords( Keywords& keys );
      96             :   explicit Average( const ActionOptions& );
      97             : };
      98             : 
      99             : PLUMED_REGISTER_ACTION(Average,"AVERAGE")
     100             : 
     101           6 : void Average::registerKeywords( Keywords& keys ) {
     102           6 :   ActionShortcut::registerKeywords( keys );
     103          12 :   keys.addInputKeyword("compulsory","ARG","scalar/grid","the quantity that is being averaged");
     104          12 :   keys.add("optional","LOGWEIGHTS","the logarithm of the quantity to use as the weights when calculating averages");
     105          12 :   keys.add("compulsory","STRIDE","1","the frequency with which to store data for averaging");
     106          12 :   keys.add("compulsory","CLEAR","0","the frequency with whihc to clear the data that is being averaged");
     107          12 :   keys.add("optional","NORMALIZATION","keyword for old version of the code that is there to maintain back compatibility only. Adding this keyword does nothing");
     108          12 :   keys.setValueDescription("scalar/grid","the value of the average");
     109          24 :   keys.needsAction("COMBINE"); keys.needsAction("CUSTOM"); keys.needsAction("ONES"); keys.needsAction("ACCUMULATE");
     110           6 : }
     111             : 
     112           4 : Average::Average( const ActionOptions& ao ):
     113             :   Action(ao),
     114           4 :   ActionShortcut(ao)
     115             : {
     116             : 
     117          16 :   std::string lw; parse("LOGWEIGHTS",lw); std::string stride, clearstride; parse("STRIDE",stride); parse("CLEAR",clearstride);
     118           4 :   if( lw.length()>0 ) {
     119           2 :     readInputLine( getShortcutLabel() + "_wsum: COMBINE ARG=" + lw + " PERIODIC=NO");
     120           2 :     readInputLine( getShortcutLabel() + "_weight: CUSTOM ARG=" + getShortcutLabel() + "_wsum FUNC=exp(x) PERIODIC=NO");
     121           6 :   } else readInputLine( getShortcutLabel() + "_weight: ONES SIZE=1" );
     122             : 
     123           8 :   std::vector<std::string> arg; parseVector("ARG",arg);
     124           4 :   if( arg.size()!=1 ) error("should only be one argument to this action");
     125           4 :   std::vector<Value*> vals; ActionWithArguments::interpretArgumentList( arg, plumed.getActionSet(), this, vals );
     126             : 
     127           8 :   readInputLine( getShortcutLabel() + "_denom: ACCUMULATE ARG=" + getShortcutLabel() + "_weight STRIDE=" + stride + " CLEAR=" + clearstride );
     128           4 :   if( vals[0]->isPeriodic() ) {
     129           2 :     std::string lbound, ubound, pfactor; vals[0]->getDomain( lbound, ubound ); pfactor = "((" + ubound + "-" + lbound + ")/(pi+pi))";
     130           2 :     readInputLine( getShortcutLabel() + "_sin: CUSTOM ARG=" + arg[0] + "," + getShortcutLabel() + "_weight FUNC=y*sin((x-" + lbound + ")/" + pfactor + ") PERIODIC=NO");
     131           2 :     readInputLine( getShortcutLabel() + "_cos: CUSTOM ARG=" + arg[0] + "," + getShortcutLabel() + "_weight FUNC=y*cos((x-" + lbound + ")/" + pfactor + ") PERIODIC=NO");
     132           2 :     readInputLine( getShortcutLabel() + "_sinsum: ACCUMULATE ARG=" + getShortcutLabel() + "_sin STRIDE=" + stride + " CLEAR=" + clearstride );
     133           2 :     readInputLine( getShortcutLabel() + "_cossum: ACCUMULATE ARG=" + getShortcutLabel() + "_cos STRIDE=" + stride + " CLEAR=" + clearstride );
     134           2 :     readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_sinsum," + getShortcutLabel() + "_cossum," + getShortcutLabel() + "_denom FUNC=" + lbound + "+" + pfactor + "*atan2(x/z,y/z) PERIODIC=" + lbound +"," + ubound);
     135             :   } else {
     136           6 :     std::string normstr; parse("NORMALIZATION",normstr);
     137           6 :     if( normstr=="true" || normstr=="false" ) warning("NORMALIZATION is deprecated. You are advised to take this out of input files in future and use the new syntax with ACCUMULATE for unormalized data rather than the shortcut AVERAGE");
     138           3 :     else if( normstr.length()>0 ) error("NORMALIZATION=" + normstr + " is not valid PLUMED input.  If you want an unormalised 'average' use ACCUMULATE");
     139           6 :     readInputLine( getShortcutLabel() + "_prod: CUSTOM ARG=" + arg[0] + "," + getShortcutLabel() + "_weight FUNC=x*y PERIODIC=NO");
     140           3 :     if( normstr.length()==0 || normstr=="true" ) {
     141           6 :       readInputLine( getShortcutLabel() + "_numer: ACCUMULATE ARG=" + getShortcutLabel() + "_prod STRIDE=" + stride + " CLEAR=" + clearstride  );
     142           6 :       readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_numer," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
     143           0 :     } else if( normstr=="false" ) readInputLine( getShortcutLabel() + ": ACCUMULATE ARG=" + getShortcutLabel() + "_prod STRIDE=" + stride + " CLEAR=" + clearstride  );
     144           0 :     else plumed_error();
     145             :   }
     146           8 : }
     147             : 
     148             : }
     149             : }