Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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/ActionRegister.h"
23 : #include "core/ActionShortcut.h"
24 : #include "core/ActionWithValue.h"
25 : #include "core/ActionWithArguments.h"
26 : #include "core/PlumedMain.h"
27 : #include "core/ActionSet.h"
28 :
29 : //+PLUMEDOC GRIDCALC HISTOGRAM
30 : /*
31 : Accumulate the average probability density along a few CVs from a trajectory.
32 :
33 : !!! caution ""
34 :
35 : We recommend not using this shortcut action and instead using a combination of [KDE](KDE.md) and
36 : [ACCUMULATE](ACCUMULATE.md) actions. This newer syntax will give you much more control over the
37 : calculation. HISTOGRAM works by calling these two underlying actions as you can see if you expand the shortcuts
38 : in the documentation that follows.
39 :
40 : When using this shortcut it is supposed that you have some collective variable $\zeta$ that
41 : gives a reasonable description of some physical or chemical phenomenon. As an example of what we
42 : mean by this suppose you wish to examine the following SN2 reaction:
43 :
44 : $$
45 : \textrm{OH}^- + \textrm{CH}_3Cl \rightarrow \textrm{CH}_3OH + \textrm{Cl}^-
46 : $$
47 :
48 : The distance between the chlorine atom and the carbon is an excellent collective variable, $\zeta$,
49 : in this case because this distance is short for the reactant, $\textrm{CH}_3Cl$, because the carbon
50 : and chlorine are chemically bonded, and because it is long for the product state when these two atoms are
51 : not chemically bonded. We thus might want to accumulate the probability density, $P(\zeta)$, as a function of this distance
52 : as this will provide us with information about the overall likelihood of the reaction. Furthermore, the
53 : free energy, $F(\zeta)$, is related to this probability density via:
54 :
55 : $$
56 : F(\zeta) = - k_B T \ln P(\zeta)
57 : $$
58 :
59 : Accumulating these probability densities is precisely what this shortcut can be used to do. Furthermore, the conversion
60 : of the histogram to the free energy can be achieved by using the method [CONVERT_TO_FES](CONVERT_TO_FES.md).
61 :
62 : We calculate histograms within PLUMED using a method known as [kernel density estimation](https://en.wikipedia.org/wiki/Kernel_density_estimation).
63 : This shortcut action thus uses the [KDE](KDE.md) and [ACCUMULATE](ACCUMULATE.md) actions to build up the time average of the histogram.
64 : __You will have much more control over what PLUMED is computing if you take some time to learn to use the new syntax with [KDE](KDE.md) and
65 : [ACCUMULATE](ACCUMULATE.md) and if you stop using the HISTOGRAM shortcut action.__
66 :
67 :
68 : In PLUMED the value of $\zeta$ at each discrete instant in time in the trajectory is accumulated. A kernel, $K(\zeta-\zeta(t'),\sigma)$,
69 : centered at the current value, $\zeta(t)$, of this quantity is generated with a bandwidth $\sigma$, which
70 : is set by the user. These kernels are then used to accumulate the ensemble average for the probability density:
71 :
72 : $$
73 : \langle P(\zeta) \rangle = \frac{ \sum_{t'=0}^t w(t') K(\zeta-\zeta(t'),\sigma) }{ \sum_{t'=0}^t w(t') }
74 : $$
75 :
76 : Here the sums run over a portion of the trajectory specified by the user. The final quantity evaluated is a weighted
77 : average as the weights, $w(t')$, allow us to negate the effect any bias might have on the region of phase space
78 : sampled by the system.
79 :
80 : A discrete analogue of kernel density estimation can also be used. In this analogue the kernels in the above formula
81 : are replaced by Dirac delta functions. When this method is used the final function calculated is no longer a probability
82 : density - it is instead a probability mass function as each element of the function tells you the value of an integral
83 : between two points on your grid rather than the value of a (continuous) function on a grid.
84 :
85 : Additional material and examples can be also found in [this tutorial](https://www.plumed-tutorials.org/lessons/21/002/data/NAVIGATION.html),
86 : which also introduces the technique known as block averaging.
87 :
88 : ## Examples
89 :
90 : The following input monitors two torsional angles during a simulation
91 : and outputs a continuous histogram as a function of them at the end of the simulation.
92 :
93 : ```plumed
94 : r1: TORSION ATOMS=1,2,3,4
95 : r2: TORSION ATOMS=2,3,4,5
96 : hh: HISTOGRAM ...
97 : ARG=r1,r2
98 : GRID_MIN=-pi,-pi
99 : GRID_MAX=pi,pi
100 : GRID_BIN=200,200
101 : BANDWIDTH=0.05,0.05
102 : ...
103 :
104 : DUMPGRID ARG=hh FILE=histo
105 : ```
106 :
107 : Instead of the number of bins to use when constructing your histogram you can specify the grid spacing to be used as shown below:
108 :
109 : ```plumed
110 : r1: TORSION ATOMS=1,2,3,4
111 : r2: TORSION ATOMS=2,3,4,5
112 : hh: HISTOGRAM ...
113 : ARG=r1,r2
114 : GRID_MIN=-pi,-pi
115 : GRID_MAX=pi,pi
116 : GRID_SPACING=0.01,0.01
117 : BANDWIDTH=0.05,0.05
118 : ...
119 :
120 : DUMPGRID ARG=hh FILE=histo
121 : ```
122 :
123 : The following input monitors two torsional angles during a simulation
124 : and outputs a discrete histogram as a function of them at the end of the simulation.
125 :
126 : ```plumed
127 : r1: TORSION ATOMS=1,2,3,4
128 : r2: TORSION ATOMS=2,3,4,5
129 : hh: HISTOGRAM ...
130 : ARG=r1,r2
131 : KERNEL=DISCRETE
132 : GRID_MIN=-pi,-pi
133 : GRID_MAX=pi,pi
134 : GRID_BIN=200,200
135 : ...
136 :
137 : DUMPGRID ARG=hh FILE=histo
138 : ```
139 :
140 : The following input monitors two torsional angles during a simulation
141 : and outputs the histogram accumulated thus far every 100000 steps.
142 :
143 : ```plumed
144 : r1: TORSION ATOMS=1,2,3,4
145 : r2: TORSION ATOMS=2,3,4,5
146 : hh: HISTOGRAM ...
147 : ARG=r1,r2
148 : GRID_MIN=-pi,-pi
149 : GRID_MAX=pi,pi
150 : GRID_BIN=200,200
151 : BANDWIDTH=0.05,0.05
152 : ...
153 :
154 : DUMPGRID ARG=hh FILE=histo STRIDE=100000
155 : ```
156 :
157 : The following input monitors two torsional angles during a simulation
158 : and outputs a separate histogram for each 100000 steps worth of trajectory.
159 : Notice how the CLEAR keyword is used here and how it is not used in the
160 : previous example.
161 :
162 : ```plumed
163 : r1: TORSION ATOMS=1,2,3,4
164 : r2: TORSION ATOMS=2,3,4,5
165 : hh: HISTOGRAM ...
166 : ARG=r1,r2 CLEAR=100000
167 : GRID_MIN=-pi,-pi
168 : GRID_MAX=pi,pi
169 : GRID_BIN=200,200
170 : BANDWIDTH=0.05,0.05
171 : ...
172 :
173 : DUMPGRID ARG=hh FILE=histo STRIDE=100000
174 : ```
175 :
176 : The following input accumulates a histogram using a subset of the data in this trajectory.
177 : When you use this input the first 10 ps of simulation time is discarded. Data is then collected
178 : in the 190 ps after that first 10 ps. The remainder of the trajectory is not used to update the histogram.
179 :
180 : ```plumed
181 : r1: TORSION ATOMS=1,2,3,4
182 : r2: TORSION ATOMS=2,3,4,5
183 : hh: HISTOGRAM ...
184 : ARG=r1,r2
185 : GRID_MIN=-pi,-pi
186 : GRID_MAX=pi,pi
187 : GRID_BIN=200,200
188 : BANDWIDTH=0.05,0.05
189 : UPDATE_FROM=10
190 : UPDATE_UNTIL=200
191 : ...
192 :
193 : DUMPGRID ARG=hh FILE=histo STRIDE=0
194 : ```
195 :
196 : In this final example there is a fixed restraint on the distance between atoms 1 and 2. Clearly, this
197 : restraint will have an effect on the region of phase space that will be sampled when an MD simulation is
198 : run using this variable. Consequently, when the histogram as a function of the distance, $x$, is accumulated,
199 : we use reweighting into order to discount the effect of the bias from our final histogram.
200 :
201 : ```plumed
202 : x: DISTANCE ATOMS=1,2
203 : RESTRAINT ARG=x SLOPE=1.0 AT=0.0
204 : bias: REWEIGHT_BIAS TEMP=300
205 :
206 : hB: HISTOGRAM ...
207 : ARG=x
208 : GRID_MIN=0.0
209 : GRID_MAX=3.0
210 : GRID_BIN=100
211 : BANDWIDTH=0.1
212 : LOGWEIGHTS=bias
213 : NORMALIZATION=true
214 : ...
215 :
216 : DUMPGRID ARG=hB FILE=histoB STRIDE=1
217 : ```
218 :
219 : */
220 : //+ENDPLUMEDOC
221 :
222 : namespace PLMD {
223 : namespace gridtools {
224 :
225 : class Histogram : public ActionShortcut {
226 : public:
227 : static void registerKeywords( Keywords& keys );
228 : explicit Histogram( const ActionOptions& );
229 : };
230 :
231 : PLUMED_REGISTER_ACTION(Histogram,"HISTOGRAM")
232 :
233 42 : void Histogram::registerKeywords( Keywords& keys ) {
234 42 : ActionShortcut::registerKeywords( keys );
235 42 : keys.use("UPDATE_FROM");
236 42 : keys.use("UPDATE_UNTIL");
237 42 : keys.add("compulsory","NORMALIZATION","ndata","This controls how the data is normalized it can be set equal to true, false or ndata. See above for an explanation");
238 84 : keys.addInputKeyword("optional","ARG","scalar/vector/matrix","the quantities that are being used to construct the histogram");
239 42 : keys.addDeprecatedKeyword("DATA","ARG");
240 42 : keys.add("compulsory","GRID_MIN","auto","the lower bounds for the grid");
241 42 : keys.add("compulsory","GRID_MAX","auto","the upper bounds for the grid");
242 42 : keys.add("optional","BANDWIDTH","the bandwidths for kernel density esimtation");
243 : std::string strcut;
244 42 : Tools::convert( sqrt(2*dp2cutoff), strcut );
245 42 : keys.add("compulsory","KERNEL","GAUSSIAN CUTOFF=" + strcut,"the switching function which is used to transform the distance from the kernel center");
246 84 : keys.linkActionInDocs("KERNEL","LESS_THAN");
247 42 : keys.add("optional","GRID_BIN","the number of bins for the grid");
248 42 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
249 42 : keys.add("optional","LOGWEIGHTS","the logarithm of the quantity to use as the weights when calculating averages");
250 42 : keys.add("compulsory","STRIDE","1","the frequency with which to store data for averaging");
251 42 : keys.add("compulsory","CLEAR","0","the frequency with which to clear the data that is being averaged");
252 84 : keys.setValueDescription("grid","the estimate of the histogram as a function of the argument that was obtained");
253 42 : keys.needsAction("COMBINE");
254 42 : keys.needsAction("CUSTOM");
255 42 : keys.needsAction("ONES");
256 42 : keys.needsAction("KDE");
257 84 : keys.needsAction("ACCUMULATE");
258 84 : }
259 :
260 33 : Histogram::Histogram( const ActionOptions& ao ):
261 : Action(ao),
262 33 : ActionShortcut(ao) {
263 : std::string normflag;
264 68 : parse("NORMALIZATION",normflag);
265 : std::string lw;
266 68 : parse("LOGWEIGHTS",lw);
267 : std::string stride, clearstride;
268 33 : parse("STRIDE",stride);
269 66 : parse("CLEAR",clearstride);
270 33 : if( lw.length()>0 && normflag!="ndata" ) {
271 12 : readInputLine( getShortcutLabel() + "_wsum: COMBINE ARG=" + lw + " PERIODIC=NO");
272 18 : readInputLine( getShortcutLabel() + "_weight: CUSTOM ARG=" + getShortcutLabel() + "_wsum FUNC=exp(x) PERIODIC=NO");
273 : } else {
274 27 : if( lw.length()>0 ) {
275 0 : error("set NORMALIZATION=true/false when using LOGWEIGHTS as otherwise the weights are ignored. Alternatively, learn to use the new syntax for histograms with KDE/ACCUMULATE to have more control over what PLUMED is calculating");
276 : }
277 56 : readInputLine( getShortcutLabel() + "_weight: ONES SIZE=1" );
278 : }
279 :
280 33 : std::vector<std::string> arglist;
281 66 : parseVector("ARG",arglist);
282 33 : if( arglist.size()==0 ) {
283 0 : error("arguments have not been specified use ARG");
284 : }
285 33 : std::vector<Value*> theargs;
286 33 : ActionWithArguments::interpretArgumentList( arglist, plumed.getActionSet(), this, theargs );
287 33 : plumed_assert( theargs.size()>0 );
288 33 : std::string argstr=theargs[0]->getName();
289 44 : for(unsigned i=1; i<theargs.size(); ++i) {
290 24 : argstr += "," + theargs[i]->getName();
291 : }
292 : std::string strnum;
293 33 : Tools::convert( theargs[0]->getNumberOfValues(), strnum );
294 33 : if( theargs[0]->getNumberOfValues()==1 ) {
295 : // Create the KDE object
296 90 : readInputLine( getShortcutLabel() + "_kde: KDE ARG=" + argstr + " " + convertInputLineToString() );
297 : } else {
298 : // Create the KDE object
299 9 : readInputLine( getShortcutLabel() + "_kde_u: KDE ARG=" + argstr + " " + convertInputLineToString() );
300 : // Normalise the KDE object
301 9 : readInputLine( getShortcutLabel() + "_kde: CUSTOM ARG=" + getShortcutLabel() + "_kde_u PERIODIC=NO FUNC=x/" + strnum );
302 : }
303 : // Now get the quantity to accumulate
304 93 : readInputLine( getShortcutLabel() + "_kdep: CUSTOM ARG=" + getShortcutLabel() + "_kde," + getShortcutLabel() + "_weight FUNC=x*y PERIODIC=NO");
305 : // And accumulate the average
306 31 : if( normflag=="false" ) {
307 35 : readInputLine( getShortcutLabel() + ": ACCUMULATE ARG=" + getShortcutLabel() + "_kdep STRIDE=" + stride + " CLEAR=" + clearstride + " " + getUpdateLimits() );
308 : } else {
309 120 : readInputLine( getShortcutLabel() + "_u: ACCUMULATE ARG=" + getShortcutLabel() + "_kdep STRIDE=" + stride + " CLEAR=" + clearstride + " " + getUpdateLimits() );
310 120 : readInputLine( getShortcutLabel() + "_nsum: ACCUMULATE ARG=" + getShortcutLabel() + "_weight STRIDE=" + stride + " CLEAR=" + clearstride + " " + getUpdateLimits() );
311 : // And divide by the total weight
312 98 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_u," + getShortcutLabel() + "_nsum FUNC=x/y PERIODIC=NO");
313 : }
314 68 : }
315 :
316 : }
317 : }
|
