Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2018 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 "tools/KernelFunctions.h"
24 : #include "gridtools/ActionWithGrid.h"
25 : #include "vesselbase/ActionWithVessel.h"
26 : #include "vesselbase/StoreDataVessel.h"
27 : #include "core/PlumedMain.h"
28 : #include "core/ActionSet.h"
29 :
30 : namespace PLMD {
31 : namespace analysis {
32 :
33 : //+PLUMEDOC GRIDCALC HISTOGRAM
34 : /*
35 : Accumulate the average probability density along a few CVs from a trajectory.
36 :
37 : When using this method it is supposed that you have some collective variable \f$\zeta\f$ that
38 : gives a reasonable description of some physical or chemical phenomenon. As an example of what we
39 : mean by this suppose you wish to examine the following SN2 reaction:
40 :
41 : \f[
42 : \textrm{OH}^- + \textrm{CH}_3Cl \rightarrow \textrm{CH}_3OH + \textrm{Cl}^-
43 : \f]
44 :
45 : The distance between the chlorine atom and the carbon is an excellent collective variable, \f$\zeta\f$,
46 : in this case because this distance is short for the reactant, \f$\textrm{CH}_3Cl\f$, because the carbon
47 : and chlorine are chemically bonded, and because it is long for the product state when these two atoms are
48 : not chemically bonded. We thus might want to accumulate the probability density, \f$P(\zeta)\f$, as a function of this distance
49 : as this will provide us with information about the overall likelihood of the reaction. Furthermore, the
50 : free energy, \f$F(\zeta)\f$, is related to this probability density via:
51 :
52 : \f[
53 : F(\zeta) = - k_B T \ln P(\zeta)
54 : \f]
55 :
56 : Accumulating these probability densities is precisely what this Action can be used to do. Furthermore, the conversion
57 : of the histogram to the free energy can be achieved by using the method \ref CONVERT_TO_FES.
58 :
59 : We calculate histograms within PLUMED using a method known as kernel density estimation, which you can read more about here:
60 :
61 : https://en.wikipedia.org/wiki/Kernel_density_estimation
62 :
63 : In PLUMED the value of \f$\zeta\f$ at each discrete instant in time in the trajectory is accumulated. A kernel, \f$K(\zeta-\zeta(t'),\sigma)\f$,
64 : centered at the current value, \f$\zeta(t)\f$, of this quantity is generated with a bandwith \f$\sigma\f$, which
65 : is set by the user. These kernels are then used to accumulate the ensemble average for the probability density:
66 :
67 : \f[
68 : \langle P(\zeta) \rangle = \frac{ \sum_{t'=0}^t w(t') K(\zeta-\zeta(t'),\sigma) }{ \sum_{t'=0}^t w(t') }
69 : \f]
70 :
71 : Here the sums run over a portion of the trajectory specified by the user. The final quantity evalulated is a weighted
72 : average as the weights, \f$w(t')\f$, allow us to negate the effect any bias might have on the region of phase space
73 : sampled by the system. This is discussed in the section of the manual on \ref Analysis.
74 :
75 : A discrete analogue of kernel density estimation can also be used. In this analogue the kernels in the above formula
76 : are replaced by dirac delta functions. When this method is used the final function calculated is no longer a probability
77 : density - it is instead a probability mass function as each element of the function tells you the value of an integral
78 : between two points on your grid rather than the value of a (continuous) function on a grid.
79 :
80 : Additional material and examples can be also found in the tutorial \ref belfast-1.
81 :
82 : \par Examples
83 :
84 : The following input monitors two torsional angles during a simulation
85 : and outputs a continuos histogram as a function of them at the end of the simulation.
86 : \verbatim
87 : TORSION ATOMS=1,2,3,4 LABEL=r1
88 : TORSION ATOMS=2,3,4,5 LABEL=r2
89 : HISTOGRAM ...
90 : ARG=r1,r2
91 : GRID_MIN=-3.14,-3.14
92 : GRID_MAX=3.14,3.14
93 : GRID_BIN=200,200
94 : BANDWIDTH=0.05,0.05
95 : LABEL=hh
96 : ... HISTOGRAM
97 :
98 : DUMPGRID GRID=hh FILE=histo
99 : \endverbatim
100 :
101 : The following input monitors two torsional angles during a simulation
102 : and outputs a discrete histogram as a function of them at the end of the simulation.
103 : \verbatim
104 : TORSION ATOMS=1,2,3,4 LABEL=r1
105 : TORSION ATOMS=2,3,4,5 LABEL=r2
106 : HISTOGRAM ...
107 : ARG=r1,r2
108 : KERNEL=DISCRETE
109 : GRID_MIN=-3.14,-3.14
110 : GRID_MAX=3.14,3.14
111 : GRID_BIN=200,200
112 : LABEL=hh
113 : ... HISTOGRAM
114 :
115 : DUMPGRID GRID=hh FILE=histo
116 : \endverbatim
117 :
118 : The following input monitors two torsional angles during a simulation
119 : and outputs the histogram accumulated thus far every 100000 steps.
120 : \verbatim
121 : TORSION ATOMS=1,2,3,4 LABEL=r1
122 : TORSION ATOMS=2,3,4,5 LABEL=r2
123 : HISTOGRAM ...
124 : ARG=r1,r2
125 : GRID_MIN=-3.14,-3.14
126 : GRID_MAX=3.14,3.14
127 : GRID_BIN=200,200
128 : BANDWIDTH=0.05,0.05
129 : LABEL=hh
130 : ... HISTOGRAM
131 :
132 : DUMPGRID GRID=hh FILE=histo STRIDE=100000
133 : \endverbatim
134 :
135 : The following input monitors two torsional angles during a simulation
136 : and outputs a separate histogram for each 100000 steps worth of trajectory.
137 : Notice how the CLEAR keyword is used here and how it is not used in the
138 : previous example.
139 :
140 : \verbatim
141 : TORSION ATOMS=1,2,3,4 LABEL=r1
142 : TORSION ATOMS=2,3,4,5 LABEL=r2
143 : HISTOGRAM ...
144 : ARG=r1,r2 CLEAR=100000
145 : GRID_MIN=-3.14,-3.14
146 : GRID_MAX=3.14,3.14
147 : GRID_BIN=200,200
148 : BANDWIDTH=0.05,0.05
149 : GRID_WFILE=histo
150 : LABEL=hh
151 : ... HISTOGRAM
152 :
153 : DUMPGRID GRID=hh FILE=histo STRIDE=100000
154 : \endverbatim
155 :
156 : */
157 : //+ENDPLUMEDOC
158 :
159 34 : class Histogram : public gridtools::ActionWithGrid {
160 : private:
161 : double ww;
162 : KernelFunctions* kernel;
163 : vesselbase::ActionWithVessel* myvessel;
164 : vesselbase::StoreDataVessel* stash;
165 : gridtools::HistogramOnGrid* myhist;
166 : public:
167 : static void registerKeywords( Keywords& keys );
168 : explicit Histogram(const ActionOptions&ao);
169 : unsigned getNumberOfQuantities() const ;
170 : void prepareForAveraging();
171 : void performOperations( const bool& from_update );
172 : void finishAveraging();
173 0 : bool isPeriodic() { return false; }
174 : unsigned getNumberOfDerivatives();
175 : void compute( const unsigned&, MultiValue& ) const ;
176 : };
177 :
178 2540 : PLUMED_REGISTER_ACTION(Histogram,"HISTOGRAM")
179 :
180 18 : void Histogram::registerKeywords( Keywords& keys ) {
181 18 : gridtools::ActionWithGrid::registerKeywords( keys ); keys.use("ARG");
182 18 : keys.add("optional","DATA","input data from action with vessel and compute histogram");
183 18 : keys.add("compulsory","GRID_MIN","the lower bounds for the grid");
184 18 : keys.add("compulsory","GRID_MAX","the upper bounds for the grid");
185 18 : keys.add("optional","GRID_BIN","the number of bins for the grid");
186 18 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
187 18 : keys.use("UPDATE_FROM"); keys.use("UPDATE_UNTIL");
188 18 : }
189 :
190 17 : Histogram::Histogram(const ActionOptions&ao):
191 : Action(ao),
192 : ActionWithGrid(ao),
193 : ww(0.0),
194 : kernel(NULL),
195 : myvessel(NULL),
196 17 : stash(NULL)
197 : {
198 : // Read in arguments
199 17 : std::string mlab; parse("DATA",mlab);
200 17 : if( mlab.length()>0 ) {
201 1 : myvessel = plumed.getActionSet().selectWithLabel<ActionWithVessel*>(mlab);
202 1 : if(!myvessel) error("action labelled " + mlab + " does not exist or is not an ActionWithVessel");
203 1 : stash = myvessel->buildDataStashes( NULL );
204 1 : log.printf(" for all base quantities calculated by %s \n",myvessel->getLabel().c_str() );
205 : // Add the dependency
206 1 : addDependency( myvessel );
207 : } else {
208 16 : std::vector<Value*> arg; parseArgumentList("ARG",arg);
209 16 : if(!arg.empty()) {
210 16 : log.printf(" with arguments");
211 16 : for(unsigned i=0; i<arg.size(); i++) log.printf(" %s",arg[i]->getName().c_str());
212 16 : log.printf("\n");
213 : // Retrieve the bias acting and make sure we request this also
214 16 : std::vector<Value*> bias( ActionWithArguments::getArguments() );
215 16 : for(unsigned i=0; i<bias.size(); ++i) arg.push_back( bias[i] );
216 16 : requestArguments(arg);
217 16 : }
218 : }
219 :
220 : // Read stuff for grid
221 17 : unsigned narg = getNumberOfArguments();
222 17 : if( myvessel ) narg=1;
223 34 : std::vector<std::string> gmin( narg ), gmax( narg );
224 17 : parseVector("GRID_MIN",gmin); parseVector("GRID_MAX",gmax);
225 34 : std::vector<unsigned> nbin; parseVector("GRID_BIN",nbin);
226 34 : std::vector<double> gspacing; parseVector("GRID_SPACING",gspacing);
227 17 : if( nbin.size()!=narg && gspacing.size()!=narg ) {
228 0 : error("GRID_BIN or GRID_SPACING must be set");
229 : }
230 :
231 : // Input of name and labels
232 34 : std::string vstring="COMPONENTS=" + getLabel();
233 17 : if( myvessel ) {
234 1 : vstring += " COORDINATES=" + myvessel->getLabel();
235 : // Input for PBC
236 1 : if( myvessel->isPeriodic() ) vstring+=" PBC=T";
237 1 : else vstring+=" PBC=F";
238 : } else {
239 16 : vstring += " COORDINATES=" + getPntrToArgument(0)->getName();
240 16 : for(unsigned i=1; i<getNumberOfArguments(); ++i) vstring += "," + getPntrToArgument(i)->getName();
241 : // Input for PBC
242 16 : if( getPntrToArgument(0)->isPeriodic() ) vstring+=" PBC=T";
243 16 : else vstring+=" PBC=F";
244 24 : for(unsigned i=1; i<getNumberOfArguments(); ++i) {
245 8 : if( getPntrToArgument(i)->isPeriodic() ) vstring+=",T";
246 8 : else vstring+=",F";
247 : }
248 : }
249 : // And create the grid
250 17 : createGrid( "histogram", vstring );
251 17 : mygrid->setBounds( gmin, gmax, nbin, gspacing );
252 17 : myhist = dynamic_cast<gridtools::HistogramOnGrid*>( mygrid );
253 17 : plumed_assert( myhist );
254 17 : if( myvessel ) {
255 : // Create a task list
256 1 : for(unsigned i=0; i<myvessel->getFullNumberOfTasks(); ++i) addTaskToList(i);
257 1 : setAveragingAction( mygrid, true );
258 : } else {
259 : // Create a task list
260 16 : for(unsigned i=0; i<mygrid->getNumberOfPoints(); ++i) addTaskToList(i);
261 16 : myhist->addOneKernelEachTimeOnly();
262 16 : setAveragingAction( mygrid, myhist->noDiscreteKernels() );
263 : }
264 34 : checkRead();
265 17 : }
266 :
267 104 : unsigned Histogram::getNumberOfDerivatives() {
268 104 : if( myvessel) return 1;
269 96 : return getNumberOfArguments();
270 : }
271 :
272 104 : unsigned Histogram::getNumberOfQuantities() const {
273 104 : if( myvessel ) return 3;
274 96 : return 2;
275 : }
276 :
277 33 : void Histogram::prepareForAveraging() {
278 33 : if( myvessel ) {
279 4 : deactivateAllTasks(); double norm=0;
280 4 : std::vector<double> cvals( myvessel->getNumberOfQuantities() );
281 12 : for(unsigned i=0; i<stash->getNumberOfStoredValues(); ++i) {
282 8 : taskFlags[i]=1; stash->retrieveSequentialValue(i, false, cvals );
283 8 : norm += cvals[0];
284 : }
285 4 : lockContributors(); ww = cweight / norm;
286 : } else {
287 : // Now fetch the kernel and the active points
288 29 : std::vector<double> point( getNumberOfArguments() );
289 29 : for(unsigned i=0; i<point.size(); ++i) point[i]=getArgument(i);
290 58 : unsigned num_neigh; std::vector<unsigned> neighbors(1);
291 29 : kernel = myhist->getKernelAndNeighbors( point, num_neigh, neighbors );
292 :
293 29 : if( num_neigh>1 ) {
294 : // Activate relevant tasks
295 24 : deactivateAllTasks();
296 24 : for(unsigned i=0; i<num_neigh; ++i) taskFlags[neighbors[i]]=1;
297 24 : lockContributors();
298 : } else {
299 : // This is used when we are not doing kernel density evaluation
300 5 : mygrid->addToGridElement( neighbors[0], 0, cweight );
301 29 : }
302 : }
303 33 : }
304 :
305 5 : void Histogram::performOperations( const bool& from_update ) { if( !myvessel ) plumed_dbg_assert( !myhist->noDiscreteKernels() ); }
306 :
307 33 : void Histogram::finishAveraging() {
308 33 : if( !myvessel ) delete kernel;
309 33 : }
310 :
311 11052 : void Histogram::compute( const unsigned& current, MultiValue& myvals ) const {
312 11052 : if( myvessel ) {
313 8 : std::vector<double> cvals( myvessel->getNumberOfQuantities() );
314 8 : stash->retrieveSequentialValue( current, false, cvals );
315 8 : myvals.setValue( 0, cvals[0] ); myvals.setValue( 1, cvals[1] ); myvals.setValue( 2, ww );
316 : } else {
317 11044 : std::vector<Value*> vv( myhist->getVectorOfValues() );
318 22083 : std::vector<double> val( getNumberOfArguments() ), der( getNumberOfArguments() );
319 : // Retrieve the location of the grid point at which we are evaluating the kernel
320 11047 : mygrid->getGridPointCoordinates( current, val );
321 11048 : for(unsigned i=0; i<getNumberOfArguments(); ++i) vv[i]->set( val[i] );
322 : // Evaluate the histogram at the relevant grid point and set the values
323 11048 : double vvh = kernel->evaluate( vv, der,true); myvals.setValue( 1, vvh );
324 : // Set the derivatives and delete the vector of values
325 22091 : for(unsigned i=0; i<getNumberOfArguments(); ++i) { myvals.setDerivative( 1, i, der[i] ); delete vv[i]; }
326 : }
327 11055 : }
328 :
329 : }
330 2523 : }
|
