Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2019 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 "vesselbase/ActionWithAveraging.h"
23 : #include "core/ActionRegister.h"
24 : #include "AverageVessel.h"
25 :
26 : //+PLUMEDOC GRIDCALC AVERAGE
27 : /*
28 : Calculate the ensemble average of a collective variable
29 :
30 : The ensemble average for a non-periodic, collective variable, \f$s\f$ is given by the following expression:
31 :
32 : \f[
33 : \langle s \rangle = \frac{ \sum_{t'=0}^t w(t') s(t') }{ \sum_{t'=0}^t w(t') }
34 : \f]
35 :
36 : Here the sum runs over a the trajectory and \f$s(t')\f$ is used to denote the value of the collective variable
37 : at time \f$t'\f$. The final quantity evalulated is a weighted
38 : average as the weights, \f$w(t')\f$, allow us to negate the effect any bias might have on the region of phase space
39 : sampled by the system. This is discussed in the section of the manual on \ref Analysis.
40 :
41 : 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:
42 :
43 : \f[
44 : \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]
45 : \f]
46 :
47 : \par Examples
48 :
49 : The following example calculates the ensemble average for the distance between atoms 1 and 2
50 : and output this to a file called COLVAR. In this example it is assumed that no bias is acting
51 : on the system and that the weights, \f$w(t')\f$ in the formulae above can thus all be set equal
52 : to one.
53 :
54 : \plumedfile
55 : d1: DISTANCE ATOMS=1,2
56 : d1a: AVERAGE ARG=d1
57 : PRINT ARG=d1a FILE=colvar STRIDE=100
58 : \endplumedfile
59 :
60 : The following example calculates the ensemble average for the torsional angle involving atoms 1, 2, 3 and 4.
61 : At variance with the previous example this quantity is periodic so the second formula in the above introduction
62 : is used to calculate the average. Furthermore, by using the CLEAR keyword we have specified that block averages
63 : are to be calculated. Consequently, after 100 steps all the information aquired thus far in the simulation is
64 : forgotten and the process of averaging is begun again. The quantities output in the colvar file are thus the
65 : block averages taken over the first 100 frames of the trajectory, the block average over the second 100 frames
66 : of trajectory and so on.
67 :
68 : \plumedfile
69 : t1: TORSION ATOMS=1,2,3,4
70 : t1a: AVERAGE ARG=t1 CLEAR=100
71 : PRINT ARG=t1a FILE=colvar STRIDE=100
72 : \endplumedfile
73 :
74 : This third example incorporates a bias. Notice that the effect the bias has on the ensemble average is removed by taking
75 : advantage of the \ref REWEIGHT_BIAS method. The final ensemble averages output to the file are thus block ensemble averages for the
76 : unbiased canononical ensemble at a temperature of 300 K.
77 :
78 : \plumedfile
79 : t1: TORSION ATOMS=1,2,3,4
80 : RESTRAINT ARG=t1 AT=pi KAPPA=100.
81 : ww: REWEIGHT_BIAS TEMP=300
82 : t1a: AVERAGE ARG=t1 LOGWEIGHTS=ww CLEAR=100
83 : PRINT ARG=t1a FILE=colvar STRIDE=100
84 : \endplumedfile
85 :
86 : */
87 : //+ENDPLUMEDOC
88 :
89 : namespace PLMD {
90 : namespace analysis {
91 :
92 4 : class Average : public vesselbase::ActionWithAveraging {
93 : private:
94 : AverageVessel* myaverage;
95 : public:
96 : static void registerKeywords( Keywords& keys );
97 : explicit Average( const ActionOptions& );
98 24 : void calculate() {}
99 24 : void apply() {}
100 : void performOperations( const bool& from_update );
101 : void finishAveraging();
102 0 : bool isPeriodic() { return false; }
103 0 : void performTask( const unsigned&, const unsigned&, MultiValue& ) const { plumed_error(); }
104 : };
105 :
106 6454 : PLUMED_REGISTER_ACTION(Average,"AVERAGE")
107 :
108 3 : void Average::registerKeywords( Keywords& keys ) {
109 6 : vesselbase::ActionWithAveraging::registerKeywords( keys ); keys.use("ARG");
110 9 : keys.remove("SERIAL"); keys.remove("LOWMEM");
111 3 : }
112 :
113 2 : Average::Average( const ActionOptions& ao ):
114 : Action(ao),
115 2 : ActionWithAveraging(ao)
116 : {
117 2 : addValue(); // Create a value so that we can output the average
118 2 : if( getNumberOfArguments()!=1 ) error("only one quantity can be averaged at a time");
119 : std::string instring;
120 2 : if( getPntrToArgument(0)->isPeriodic() ) {
121 1 : std::string min, max; getPntrToArgument(0)->getDomain(min,max);
122 4 : instring = "PERIODIC=" + min + "," + max; setPeriodic( min, max );
123 : } else {
124 1 : setNotPeriodic();
125 : }
126 : // Create a vessel to hold the average
127 8 : vesselbase::VesselOptions da("myaverage","",-1,instring,this);
128 4 : Keywords keys; AverageVessel::registerKeywords( keys );
129 4 : vesselbase::VesselOptions dar( da, keys );
130 2 : myaverage = new AverageVessel(dar); setAveragingAction( myaverage, false );
131 2 : }
132 :
133 22 : void Average::performOperations( const bool& from_update ) {
134 44 : myaverage->accumulate( cweight, getArgument(0) );
135 22 : }
136 :
137 22 : void Average::finishAveraging() {
138 22 : setValue( myaverage->getAverage() );
139 22 : }
140 :
141 : }
142 4839 : }
|
