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, $s$ is given by the following expression:
      32             : 
      33             : $$
      34             : \langle s \rangle = \frac{ \sum_{t'=0}^t w(t') s(t') }{ \sum_{t'=0}^t w(t') }
      35             : $$
      36             : 
      37             : Here the sum runs over a the trajectory and $s(t')$ is used to denote the value of the collective variable
      38             : at time $t'$.  The final quantity evaluated is a weighted
      39             : average as the weights, If the simulation is unbiassed then all the $w(t')$ values in teh expression above are
      40             : zero.  If the simulation is biased then the $w(t')$ weights are set in a way that ensures the effect any bias
      41             : has on the region of phase space sampled by the system is negated.
      42             : 
      43             : As the following example input shows you can use the AVERAGE shortcut to calculate the ensemble average of a CV using this formula:
      44             : 
      45             : ```plumed
      46             : d1: DISTANCE ATOMS=1,2
      47             : d1a: AVERAGE ARG=d1
      48             : PRINT ARG=d1a FILE=colvar STRIDE=100
      49             : ```
      50             : 
      51             : In this example no bias is acting on the system so the weights, $w(t')$ in the formulas above can thus all be set equal
      52             : to one. The shortcut illustrates how the averaging is achieved by using the [ACCUMULATE](ACCUMULATE.md) action.
      53             : 
      54             : When the variable is periodic (e.g. [TORSION](TORSION.md)) and has a value, $s$, in $a \le s \le b$ the ensemble average is evaluated using:
      55             : 
      56             : $$
      57             : \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]
      58             : $$
      59             : 
      60             : You can see how [ACCUMULATE](ACCUMULATE.md) and [CUSTOM](CUSTOM.md) can be used to implement this formula by expanding the shortcuts in the following example input:
      61             : 
      62             : ```plumed
      63             : t1: TORSION ATOMS=1,2,3,4
      64             : t1a: AVERAGE ARG=t1 CLEAR=100
      65             : PRINT ARG=t1a FILE=colvar STRIDE=100
      66             : ```
      67             : 
      68             : Notice that by using the `CLEAR` keyword we have specified that block averages
      69             : are to be calculated.  Consequently, after 100 steps all the information acquired thus far in the simulation is
      70             : forgotten and the process of averaging is begun again.  The quantities output in the colvar file are thus the
      71             : block averages taken over the first 100 frames of the trajectory, the block average over the second 100 frames
      72             : of trajectory and so on.
      73             : 
      74             : 
      75             : If a bias is acting upon the system then the $w(t')$ values in the expression above are non-zero.  You can calculate the $w(t') values
      76             : by using [REWEIGHT_BIAS](REWEIGHT_BIAS.md) or similar.  To pass these weights to the average action you would then use an input something
      77             : like the following:
      78             : 
      79             : ```plumed
      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             : ```
      86             : 
      87             : This AVERAGE action in this input is a shortcut once again so by expanding it you can obtain a better understanding of how the
      88             : formulas above are applied in this case.
      89             : 
      90             : */
      91             : //+ENDPLUMEDOC
      92             : 
      93             : namespace PLMD {
      94             : namespace generic {
      95             : 
      96             : class Average : public ActionShortcut {
      97             : public:
      98             :   static void registerKeywords( Keywords& keys );
      99             :   explicit Average( const ActionOptions& );
     100             : };
     101             : 
     102             : PLUMED_REGISTER_ACTION(Average,"AVERAGE")
     103             : 
     104           6 : void Average::registerKeywords( Keywords& keys ) {
     105           6 :   ActionShortcut::registerKeywords( keys );
     106          12 :   keys.addInputKeyword("compulsory","ARG","scalar/grid","the quantity that is being averaged");
     107           6 :   keys.add("optional","LOGWEIGHTS","the logarithm of the quantity to use as the weights when calculating averages");
     108           6 :   keys.add("compulsory","STRIDE","1","the frequency with which to store data for averaging");
     109           6 :   keys.add("compulsory","CLEAR","0","the frequency with whihc to clear the data that is being averaged");
     110           6 :   keys.add("optional","NORMALIZATION","keyword for old version of the code that is there to maintain back compatibility only. Adding this keyword does nothing");
     111          12 :   keys.setValueDescription("scalar/grid","the value of the average");
     112           6 :   keys.needsAction("COMBINE");
     113           6 :   keys.needsAction("CUSTOM");
     114           6 :   keys.needsAction("ONES");
     115           6 :   keys.needsAction("ACCUMULATE");
     116           6 : }
     117             : 
     118           4 : Average::Average( const ActionOptions& ao ):
     119             :   Action(ao),
     120           4 :   ActionShortcut(ao) {
     121             : 
     122             :   std::string lw;
     123           8 :   parse("LOGWEIGHTS",lw);
     124             :   std::string stride, clearstride;
     125           4 :   parse("STRIDE",stride);
     126           8 :   parse("CLEAR",clearstride);
     127           4 :   if( lw.length()>0 ) {
     128           2 :     readInputLine( getShortcutLabel() + "_wsum: COMBINE ARG=" + lw + " PERIODIC=NO");
     129           2 :     readInputLine( getShortcutLabel() + "_weight: CUSTOM ARG=" + getShortcutLabel() + "_wsum FUNC=exp(x) PERIODIC=NO");
     130             :   } else {
     131           6 :     readInputLine( getShortcutLabel() + "_weight: ONES SIZE=1" );
     132             :   }
     133             : 
     134             :   std::vector<std::string> arg;
     135           8 :   parseVector("ARG",arg);
     136           4 :   if( arg.size()!=1 ) {
     137           0 :     error("should only be one argument to this action");
     138             :   }
     139             :   std::vector<Value*> vals;
     140           4 :   ActionWithArguments::interpretArgumentList( arg, plumed.getActionSet(), this, vals );
     141             : 
     142           8 :   readInputLine( getShortcutLabel() + "_denom: ACCUMULATE ARG=" + getShortcutLabel() + "_weight STRIDE=" + stride + " CLEAR=" + clearstride );
     143           4 :   if( vals[0]->isPeriodic() ) {
     144             :     std::string lbound, ubound, pfactor;
     145           1 :     vals[0]->getDomain( lbound, ubound );
     146           2 :     pfactor = "((" + ubound + "-" + lbound + ")/(pi+pi))";
     147           2 :     readInputLine( getShortcutLabel() + "_sin: CUSTOM ARG=" + arg[0] + "," + getShortcutLabel() + "_weight FUNC=y*sin((x-" + lbound + ")/" + pfactor + ") PERIODIC=NO");
     148           2 :     readInputLine( getShortcutLabel() + "_cos: CUSTOM ARG=" + arg[0] + "," + getShortcutLabel() + "_weight FUNC=y*cos((x-" + lbound + ")/" + pfactor + ") PERIODIC=NO");
     149           2 :     readInputLine( getShortcutLabel() + "_sinsum: ACCUMULATE ARG=" + getShortcutLabel() + "_sin STRIDE=" + stride + " CLEAR=" + clearstride );
     150           2 :     readInputLine( getShortcutLabel() + "_cossum: ACCUMULATE ARG=" + getShortcutLabel() + "_cos STRIDE=" + stride + " CLEAR=" + clearstride );
     151           2 :     readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_sinsum," + getShortcutLabel() + "_cossum," + getShortcutLabel() + "_denom FUNC=" + lbound + "+" + pfactor + "*atan2(x/z,y/z) PERIODIC=" + lbound +"," + ubound);
     152             :   } else {
     153             :     std::string normstr;
     154           6 :     parse("NORMALIZATION",normstr);
     155           6 :     if( normstr=="true" || normstr=="false" ) {
     156           0 :       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");
     157           3 :     } else if( normstr.length()>0 ) {
     158           0 :       error("NORMALIZATION=" + normstr + " is not valid PLUMED input.  If you want an unormalised 'average' use ACCUMULATE");
     159             :     }
     160           6 :     readInputLine( getShortcutLabel() + "_prod: CUSTOM ARG=" + arg[0] + "," + getShortcutLabel() + "_weight FUNC=x*y PERIODIC=NO");
     161           3 :     if( normstr.length()==0 || normstr=="true" ) {
     162           6 :       readInputLine( getShortcutLabel() + "_numer: ACCUMULATE ARG=" + getShortcutLabel() + "_prod STRIDE=" + stride + " CLEAR=" + clearstride  );
     163           6 :       readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_numer," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
     164           0 :     } else if( normstr=="false" ) {
     165           0 :       readInputLine( getShortcutLabel() + ": ACCUMULATE ARG=" + getShortcutLabel() + "_prod STRIDE=" + stride + " CLEAR=" + clearstride  );
     166             :     } else {
     167           0 :       plumed_error();
     168             :     }
     169             :   }
     170           8 : }
     171             : 
     172             : }
     173             : }