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 "tools/KernelFunctions.h"
24 : #include "gridtools/ActionWithGrid.h"
25 : #include "vesselbase/ActionWithVessel.h"
26 : #include "vesselbase/StoreDataVessel.h"
27 : #include "multicolvar/MultiColvarBase.h"
28 : #include "core/PlumedMain.h"
29 : #include "core/ActionSet.h"
30 :
31 : namespace PLMD {
32 : namespace analysis {
33 :
34 : //+PLUMEDOC GRIDCALC HISTOGRAM
35 : /*
36 : Accumulate the average probability density along a few CVs from a trajectory.
37 :
38 : When using this method it is supposed that you have some collective variable \f$\zeta\f$ that
39 : gives a reasonable description of some physical or chemical phenomenon. As an example of what we
40 : mean by this suppose you wish to examine the following SN2 reaction:
41 :
42 : \f[
43 : \textrm{OH}^- + \textrm{CH}_3Cl \rightarrow \textrm{CH}_3OH + \textrm{Cl}^-
44 : \f]
45 :
46 : The distance between the chlorine atom and the carbon is an excellent collective variable, \f$\zeta\f$,
47 : in this case because this distance is short for the reactant, \f$\textrm{CH}_3Cl\f$, because the carbon
48 : and chlorine are chemically bonded, and because it is long for the product state when these two atoms are
49 : not chemically bonded. We thus might want to accumulate the probability density, \f$P(\zeta)\f$, as a function of this distance
50 : as this will provide us with information about the overall likelihood of the reaction. Furthermore, the
51 : free energy, \f$F(\zeta)\f$, is related to this probability density via:
52 :
53 : \f[
54 : F(\zeta) = - k_B T \ln P(\zeta)
55 : \f]
56 :
57 : Accumulating these probability densities is precisely what this Action can be used to do. Furthermore, the conversion
58 : of the histogram to the free energy can be achieved by using the method \ref CONVERT_TO_FES.
59 :
60 : We calculate histograms within PLUMED using a method known as kernel density estimation, which you can read more about here:
61 :
62 : https://en.wikipedia.org/wiki/Kernel_density_estimation
63 :
64 : 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$,
65 : centered at the current value, \f$\zeta(t)\f$, of this quantity is generated with a bandwidth \f$\sigma\f$, which
66 : is set by the user. These kernels are then used to accumulate the ensemble average for the probability density:
67 :
68 : \f[
69 : \langle P(\zeta) \rangle = \frac{ \sum_{t'=0}^t w(t') K(\zeta-\zeta(t'),\sigma) }{ \sum_{t'=0}^t w(t') }
70 : \f]
71 :
72 : Here the sums run over a portion of the trajectory specified by the user. The final quantity evaluated is a weighted
73 : average as the weights, \f$w(t')\f$, allow us to negate the effect any bias might have on the region of phase space
74 : sampled by the system. This is discussed in the section of the manual on \ref Analysis.
75 :
76 : A discrete analogue of kernel density estimation can also be used. In this analogue the kernels in the above formula
77 : are replaced by Dirac delta functions. When this method is used the final function calculated is no longer a probability
78 : density - it is instead a probability mass function as each element of the function tells you the value of an integral
79 : between two points on your grid rather than the value of a (continuous) function on a grid.
80 :
81 : Additional material and examples can be also found in the tutorials \ref lugano-1.
82 :
83 : \par A note on block averaging and errors
84 :
85 : Some particularly important
86 : issues related to the convergence of histograms and the estimation of error bars around the ensemble averages you calculate are covered in \ref trieste-2.
87 : The technique for estimating error bars that is known as block averaging is introduced in this tutorial. The essence of this technique is that
88 : the trajectory is split into a set of blocks and separate ensemble averages are calculated from each separate block of data. If \f$\{A_i\}\f$ is
89 : the set of \f$N\f$ block averages that are obtained from this technique then the final error bar is calculated as:
90 :
91 : \f[
92 : \textrm{error} = \sqrt{ \frac{1}{N} \frac{1}{N-1} \sum_{i=1}^N (A_i^2 - \langle A \rangle )^2 } \qquad \textrm{where} \qquad \langle A \rangle = \frac{1}{N} \sum_{i=1}^N A_i
93 : \f]
94 :
95 : If the simulation is biased and reweighting is performed then life is a little more complex as each of the block averages should be calculated as a
96 : weighted average. Furthermore, the weights should be taken into account when the final ensemble and error bars are calculated. As such the error should be:
97 :
98 : \f[
99 : \textrm{error} = \sqrt{ \frac{1}{N} \frac{\sum_{i=1}^N W_i }{\sum_{i=1}^N W_i - \sum_{i=1}^N W_i^2 / \sum_{i=1}^N W_i} \sum_{i=1}^N W_i (A_i^2 - \langle A \rangle )^2 }
100 : \f]
101 :
102 : where \f$W_i\f$ is the sum of all the weights for the \f$i\f$th block of data.
103 :
104 : If we wish to calculate a normalized histogram we must calculate ensemble averages from our biased simulation using:
105 : \f[
106 : \langle H(x) \rangle = \frac{\sum_{t=1}^M w_t K( x - x_t,\sigma) }{\sum_{t=1}^M w_t}
107 : \f]
108 : where the sums runs over the trajectory, \f$w_t\f$ is the weight of the \f$t\f$th trajectory frame, \f$x_t\f$ is the value of the CV for the \f$t\f$th
109 : trajectory frame and \f$K\f$ is a kernel function centered on \f$x_t\f$ with bandwidth \f$\sigma\f$. The quantity that is evaluated is the value of the
110 : normalized histogram at point \f$x\f$. The following ensemble average will be calculated if you use the NORMALIZATION=true option in HISTOGRAM.
111 : If the ensemble average is calculated in this way we must calculate the associated error bars from our block averages using the second of the expressions
112 : above.
113 :
114 : A number of works have shown that when biased simulations are performed it is often better to calculate an estimate of the histogram that is not normalized using:
115 : \f[
116 : \langle H(x) \rangle = \frac{1}{M} \sum_{t=1}^M w_t K( x - x_t,\sigma)
117 : \f]
118 : instead of the expression above. As such this is what is done by default in HISTOGRAM or if the NORMALIZATION=ndata option is used.
119 : When the histogram is calculated in this second way the first of the two formula above can be used when calculating error bars from
120 : block averages.
121 :
122 : \par Examples
123 :
124 : The following input monitors two torsional angles during a simulation
125 : and outputs a continuous histogram as a function of them at the end of the simulation.
126 : \plumedfile
127 : TORSION ATOMS=1,2,3,4 LABEL=r1
128 : TORSION ATOMS=2,3,4,5 LABEL=r2
129 : HISTOGRAM ...
130 : ARG=r1,r2
131 : GRID_MIN=-3.14,-3.14
132 : GRID_MAX=3.14,3.14
133 : GRID_BIN=200,200
134 : BANDWIDTH=0.05,0.05
135 : LABEL=hh
136 : ... HISTOGRAM
137 :
138 : DUMPGRID GRID=hh FILE=histo
139 : \endplumedfile
140 :
141 : The following input monitors two torsional angles during a simulation
142 : and outputs a discrete histogram as a function of them at the end of the simulation.
143 : \plumedfile
144 : TORSION ATOMS=1,2,3,4 LABEL=r1
145 : TORSION ATOMS=2,3,4,5 LABEL=r2
146 : HISTOGRAM ...
147 : ARG=r1,r2
148 : KERNEL=DISCRETE
149 : GRID_MIN=-3.14,-3.14
150 : GRID_MAX=3.14,3.14
151 : GRID_BIN=200,200
152 : LABEL=hh
153 : ... HISTOGRAM
154 :
155 : DUMPGRID GRID=hh FILE=histo
156 : \endplumedfile
157 :
158 : The following input monitors two torsional angles during a simulation
159 : and outputs the histogram accumulated thus far every 100000 steps.
160 : \plumedfile
161 : TORSION ATOMS=1,2,3,4 LABEL=r1
162 : TORSION ATOMS=2,3,4,5 LABEL=r2
163 : HISTOGRAM ...
164 : ARG=r1,r2
165 : GRID_MIN=-3.14,-3.14
166 : GRID_MAX=3.14,3.14
167 : GRID_BIN=200,200
168 : BANDWIDTH=0.05,0.05
169 : LABEL=hh
170 : ... HISTOGRAM
171 :
172 : DUMPGRID GRID=hh FILE=histo STRIDE=100000
173 : \endplumedfile
174 :
175 : The following input monitors two torsional angles during a simulation
176 : and outputs a separate histogram for each 100000 steps worth of trajectory.
177 : Notice how the CLEAR keyword is used here and how it is not used in the
178 : previous example.
179 :
180 : \plumedfile
181 : TORSION ATOMS=1,2,3,4 LABEL=r1
182 : TORSION ATOMS=2,3,4,5 LABEL=r2
183 : HISTOGRAM ...
184 : ARG=r1,r2 CLEAR=100000
185 : GRID_MIN=-3.14,-3.14
186 : GRID_MAX=3.14,3.14
187 : GRID_BIN=200,200
188 : BANDWIDTH=0.05,0.05
189 : LABEL=hh
190 : ... HISTOGRAM
191 :
192 : DUMPGRID GRID=hh FILE=histo STRIDE=100000
193 : \endplumedfile
194 :
195 : */
196 : //+ENDPLUMEDOC
197 :
198 : class Histogram : public gridtools::ActionWithGrid {
199 : private:
200 : double ww;
201 : bool in_apply, mvectors;
202 : std::unique_ptr<KernelFunctions> kernel;
203 : std::vector<double> forcesToApply, finalForces;
204 : std::vector<vesselbase::ActionWithVessel*> myvessels;
205 : std::vector<vesselbase::StoreDataVessel*> stashes;
206 : // this is a plain pointer since this object is now owned
207 : gridtools::HistogramOnGrid* myhist;
208 : public:
209 : static void registerKeywords( Keywords& keys );
210 : explicit Histogram(const ActionOptions&ao);
211 : unsigned getNumberOfQuantities() const override;
212 : void prepareForAveraging() override;
213 : void performOperations( const bool& from_update ) override;
214 : void finishAveraging() override;
215 121 : bool threadSafe() const override {
216 121 : return !in_apply;
217 : }
218 0 : bool isPeriodic() override {
219 0 : return false;
220 : }
221 : unsigned getNumberOfDerivatives() override;
222 : void turnOnDerivatives() override;
223 : void compute( const unsigned&, MultiValue& ) const override;
224 : void apply() override;
225 : };
226 :
227 13853 : PLUMED_REGISTER_ACTION(Histogram,"HISTOGRAM")
228 :
229 38 : void Histogram::registerKeywords( Keywords& keys ) {
230 38 : gridtools::ActionWithGrid::registerKeywords( keys );
231 38 : keys.use("ARG");
232 38 : keys.remove("NORMALIZATION");
233 76 : 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");
234 76 : keys.add("optional","DATA","input data from action with vessel and compute histogram");
235 76 : keys.add("optional","VECTORS","input three dimensional vectors for computing histogram");
236 76 : keys.add("compulsory","GRID_MIN","the lower bounds for the grid");
237 76 : keys.add("compulsory","GRID_MAX","the upper bounds for the grid");
238 76 : keys.add("optional","GRID_BIN","the number of bins for the grid");
239 76 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
240 38 : keys.use("UPDATE_FROM");
241 38 : keys.use("UPDATE_UNTIL");
242 38 : }
243 :
244 34 : Histogram::Histogram(const ActionOptions&ao):
245 : Action(ao),
246 : ActionWithGrid(ao),
247 34 : ww(0.0),
248 34 : in_apply(false),
249 34 : mvectors(false) {
250 : // Read in arguments
251 34 : if( getNumberOfArguments()==0 ) {
252 : std::string vlab;
253 16 : parse("VECTORS",vlab);
254 8 : if( vlab.length()>0 ) {
255 2 : ActionWithVessel* myv = plumed.getActionSet().selectWithLabel<ActionWithVessel*>( vlab );
256 2 : if( !myv ) {
257 0 : error("action labelled " + vlab + " does not exist or is not an ActionWithVessel");
258 : }
259 2 : myvessels.push_back( myv );
260 2 : stashes.push_back( myv->buildDataStashes( NULL ) );
261 2 : addDependency( myv );
262 2 : mvectors=true;
263 2 : if( myv->getNumberOfQuantities()!=5 ) {
264 0 : error("can only compute histograms for three dimensional vectors");
265 : }
266 2 : log.printf(" for vector quantities calculated by %s \n", vlab.c_str() );
267 : } else {
268 : std::vector<std::string> mlab;
269 12 : parseVector("DATA",mlab);
270 6 : if( mlab.size()>0 ) {
271 16 : for(unsigned i=0; i<mlab.size(); ++i) {
272 10 : ActionWithVessel* myv = plumed.getActionSet().selectWithLabel<ActionWithVessel*>( mlab[i] );
273 10 : if( !myv ) {
274 0 : error("action labelled " + mlab[i] + " does not exist or is not an ActionWithVessel");
275 : }
276 10 : myvessels.push_back( myv );
277 10 : stashes.push_back( myv->buildDataStashes( NULL ) );
278 : // log.printf(" for all base quantities calculated by %s \n",myvessel->getLabel().c_str() );
279 : // Add the dependency
280 10 : addDependency( myv );
281 : }
282 6 : unsigned nvals = myvessels[0]->getFullNumberOfTasks();
283 10 : for(unsigned i=1; i<mlab.size(); ++i) {
284 4 : if( nvals!=myvessels[i]->getFullNumberOfTasks() ) {
285 0 : error("mismatched number of quantities calculated by actions input to histogram");
286 : }
287 : }
288 6 : log.printf(" for all base quantities calculated by %s ", myvessels[0]->getLabel().c_str() );
289 10 : for(unsigned i=1; i<mlab.size(); ++i) {
290 4 : log.printf(", %s \n", myvessels[i]->getLabel().c_str() );
291 : }
292 6 : log.printf("\n");
293 : } else {
294 0 : error("input data is missing use ARG/VECTORS/DATA");
295 : }
296 6 : }
297 : }
298 :
299 : // Read stuff for grid
300 : unsigned narg = getNumberOfArguments();
301 34 : if( myvessels.size()>0 ) {
302 8 : narg=myvessels.size();
303 : }
304 : // Input of name and labels
305 34 : std::string vstring="COMPONENTS=" + getLabel();
306 34 : if( mvectors ) {
307 : vstring += " COORDINATES=x,y,z PBC=F,F,F";
308 32 : } else if( myvessels.size()>0 ) {
309 6 : vstring += " COORDINATES=" + myvessels[0]->getLabel();
310 10 : for(unsigned i=1; i<myvessels.size(); ++i) {
311 8 : vstring +="," + myvessels[i]->getLabel();
312 : }
313 : // Input for PBC
314 6 : if( myvessels[0]->isPeriodic() ) {
315 : vstring+=" PBC=T";
316 : } else {
317 : vstring+=" PBC=F";
318 : }
319 10 : for(unsigned i=1; i<myvessels.size(); ++i) {
320 4 : if( myvessels[i]->isPeriodic() ) {
321 : vstring+=",T";
322 : } else {
323 : vstring+=",F";
324 : }
325 : }
326 : } else {
327 26 : vstring += " COORDINATES=" + getPntrToArgument(0)->getName();
328 35 : for(unsigned i=1; i<getNumberOfArguments(); ++i) {
329 18 : vstring += "," + getPntrToArgument(i)->getName();
330 : }
331 : // Input for PBC
332 26 : if( getPntrToArgument(0)->isPeriodic() ) {
333 : vstring+=" PBC=T";
334 : } else {
335 : vstring+=" PBC=F";
336 : }
337 35 : for(unsigned i=1; i<getNumberOfArguments(); ++i) {
338 9 : if( getPntrToArgument(i)->isPeriodic() ) {
339 : vstring+=",T";
340 : } else {
341 : vstring+=",F";
342 : }
343 : }
344 : }
345 : // And create the grid
346 34 : auto grid=createGrid( "histogram", vstring );
347 : // notice that createGrid also sets mygrid=grid.get()
348 68 : if( mygrid->getType()=="flat" ) {
349 32 : if( mvectors ) {
350 0 : error("computing histogram for three dimensional vectors but grid is not of fibonacci type - use CONCENTRATION");
351 : }
352 32 : std::vector<std::string> gmin( narg ), gmax( narg );
353 32 : parseVector("GRID_MIN",gmin);
354 64 : parseVector("GRID_MAX",gmax);
355 : std::vector<unsigned> nbin;
356 64 : parseVector("GRID_BIN",nbin);
357 : std::vector<double> gspacing;
358 64 : parseVector("GRID_SPACING",gspacing);
359 32 : if( nbin.size()!=narg && gspacing.size()!=narg ) {
360 0 : error("GRID_BIN or GRID_SPACING must be set");
361 : }
362 32 : mygrid->setBounds( gmin, gmax, nbin, gspacing );
363 32 : } else {
364 : std::vector<unsigned> nbin;
365 4 : parseVector("GRID_BIN",nbin);
366 2 : if( nbin.size()!=1 ) {
367 0 : error("should only be one index for number of bins with spherical grid");
368 : }
369 4 : if( mygrid->getType()=="fibonacci" ) {
370 2 : mygrid->setupFibonacciGrid( nbin[0] );
371 : }
372 : }
373 34 : myhist = dynamic_cast<gridtools::HistogramOnGrid*>( mygrid );
374 34 : plumed_assert( myhist );
375 34 : if( myvessels.size()>0 ) {
376 : // Create a task list
377 72 : for(unsigned i=0; i<myvessels[0]->getFullNumberOfTasks(); ++i) {
378 64 : addTaskToList(i);
379 : }
380 8 : setAveragingAction( std::move(grid), true );
381 : } else if( storeThenAverage() ) {
382 7 : setAveragingAction( std::move(grid), true );
383 : } else {
384 : // Create a task list
385 9256 : for(unsigned i=0; i<mygrid->getNumberOfPoints(); ++i) {
386 9237 : addTaskToList(i);
387 : }
388 19 : myhist->addOneKernelEachTimeOnly();
389 19 : setAveragingAction( std::move(grid), myhist->noDiscreteKernels() );
390 : }
391 34 : checkRead();
392 68 : }
393 :
394 8 : void Histogram::turnOnDerivatives() {
395 8 : ActionWithGrid::turnOnDerivatives();
396 : std::vector<AtomNumber> all_atoms, tmp_atoms;
397 20 : for(unsigned i=0; i<myvessels.size(); ++i) {
398 12 : multicolvar::MultiColvarBase* mbase=dynamic_cast<multicolvar::MultiColvarBase*>( myvessels[i] );
399 12 : if( !mbase ) {
400 0 : error("do not know how to get histogram derivatives for actions of type " + myvessels[i]->getName() );
401 : }
402 12 : tmp_atoms = mbase->getAbsoluteIndexes();
403 482 : for(unsigned j=0; j<tmp_atoms.size(); ++j) {
404 470 : all_atoms.push_back( tmp_atoms[j] );
405 : }
406 : // Make a tempory multi value so we can avoid vector resizing
407 12 : stashes[i]->resizeTemporyMultiValues( 1 );
408 : }
409 8 : ActionAtomistic::requestAtoms( all_atoms );
410 8 : finalForces.resize( 3*all_atoms.size() + 9 );
411 8 : forcesToApply.resize( 3*all_atoms.size() + 9*myvessels.size() );
412 : // And make sure we still have the dependencies which are cleared by requestAtoms
413 20 : for(unsigned i=0; i<myvessels.size(); ++i) {
414 12 : addDependency( myvessels[i] );
415 : }
416 : // And resize the histogram so that we have a place to store the forces
417 8 : in_apply=true;
418 8 : mygrid->resize();
419 8 : in_apply=false;
420 8 : }
421 :
422 728 : unsigned Histogram::getNumberOfDerivatives() {
423 728 : if( in_apply ) {
424 : unsigned nder=0;
425 540 : for(unsigned i=0; i<myvessels.size(); ++i) {
426 318 : nder += myvessels[i]->getNumberOfDerivatives();
427 : }
428 : return nder;
429 : }
430 506 : return getNumberOfArguments();
431 : }
432 :
433 462 : unsigned Histogram::getNumberOfQuantities() const {
434 462 : if( mvectors ) {
435 46 : return myvessels[0]->getNumberOfQuantities();
436 416 : } else if( myvessels.size()>0 ) {
437 122 : return myvessels.size()+2;
438 : }
439 294 : return ActionWithAveraging::getNumberOfQuantities();
440 : }
441 :
442 115 : void Histogram::prepareForAveraging() {
443 115 : if( myvessels.size()>0 ) {
444 36 : deactivateAllTasks();
445 : double norm=0;
446 332 : for(unsigned i=0; i<stashes[0]->getNumberOfStoredValues(); ++i) {
447 296 : std::vector<double> cvals( myvessels[0]->getNumberOfQuantities() );
448 296 : stashes[0]->retrieveSequentialValue( i, false, cvals );
449 296 : unsigned itask=myvessels[0]->getActiveTask(i);
450 296 : double tnorm = cvals[0];
451 444 : for(unsigned j=1; j<myvessels.size(); ++j) {
452 148 : if( myvessels[j]->getActiveTask(i)!=itask ) {
453 0 : error("mismatched task identities in histogram suggests histogram is meaningless");
454 : }
455 148 : if( cvals.size()!=myvessels[j]->getNumberOfQuantities() ) {
456 8 : cvals.resize( myvessels[j]->getNumberOfQuantities() );
457 : }
458 148 : stashes[j]->retrieveSequentialValue( i, false, cvals );
459 148 : tnorm *= cvals[0];
460 : }
461 296 : norm += tnorm;
462 296 : taskFlags[i]=1;
463 : }
464 36 : lockContributors();
465 : // Sort out normalization of histogram
466 36 : if( !noNormalization() ) {
467 31 : ww = cweight / norm;
468 : } else {
469 5 : ww = cweight;
470 : }
471 : } else if( !storeThenAverage() ) {
472 : // Now fetch the kernel and the active points
473 72 : std::vector<double> point( getNumberOfArguments() );
474 180 : for(unsigned i=0; i<point.size(); ++i) {
475 108 : point[i]=getArgument(i);
476 : }
477 : unsigned num_neigh;
478 72 : std::vector<unsigned> neighbors(1);
479 144 : kernel=myhist->getKernelAndNeighbors( point, num_neigh, neighbors );
480 :
481 72 : if( num_neigh>1 ) {
482 : // Activate relevant tasks
483 67 : deactivateAllTasks();
484 14975 : for(unsigned i=0; i<num_neigh; ++i) {
485 14908 : taskFlags[neighbors[i]]=1;
486 : }
487 67 : lockContributors();
488 : } else {
489 : // This is used when we are not doing kernel density evaluation
490 5 : mygrid->addToGridElement( neighbors[0], 0, cweight );
491 : }
492 : }
493 115 : }
494 :
495 5 : void Histogram::performOperations( const bool& from_update ) {
496 : if( myvessels.size()==0 ) {
497 : plumed_dbg_assert( !myhist->noDiscreteKernels() );
498 : }
499 5 : }
500 :
501 135 : void Histogram::finishAveraging() {
502 135 : if( myvessels.size()==0 ) {
503 : kernel.reset();
504 : }
505 135 : }
506 :
507 15404 : void Histogram::compute( const unsigned& current, MultiValue& myvals ) const {
508 15404 : if( mvectors ) {
509 140 : std::vector<double> cvals( myvessels[0]->getNumberOfQuantities() );
510 140 : stashes[0]->retrieveSequentialValue( current, true, cvals );
511 560 : for(unsigned i=2; i<myvessels[0]->getNumberOfQuantities(); ++i) {
512 420 : myvals.setValue( i-1, cvals[i] );
513 : }
514 : myvals.setValue( 0, cvals[0] );
515 140 : myvals.setValue( myvessels[0]->getNumberOfQuantities() - 1, ww );
516 140 : if( in_apply ) {
517 50 : MultiValue& tmpval = stashes[0]->getTemporyMultiValue(0);
518 99 : if( tmpval.getNumberOfValues()!=myvessels[0]->getNumberOfQuantities() ||
519 49 : tmpval.getNumberOfDerivatives()!=myvessels[0]->getNumberOfDerivatives() ) {
520 1 : tmpval.resize( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
521 : }
522 50 : stashes[0]->retrieveDerivatives( stashes[0]->getTrueIndex(current), true, tmpval );
523 350 : for(unsigned j=0; j<tmpval.getNumberActive(); ++j) {
524 300 : unsigned jder=tmpval.getActiveIndex(j);
525 300 : myvals.addDerivative( 0, jder, tmpval.getDerivative(0, jder) );
526 1200 : for(unsigned i=2; i<myvessels[0]->getNumberOfQuantities(); ++i) {
527 900 : myvals.addDerivative( i-1, jder, tmpval.getDerivative(i, jder) );
528 : }
529 : }
530 : myvals.updateDynamicList();
531 : }
532 15264 : } else if( myvessels.size()>0 ) {
533 356 : std::vector<double> cvals( myvessels[0]->getNumberOfQuantities() );
534 356 : stashes[0]->retrieveSequentialValue( current, false, cvals );
535 : unsigned derbase=0;
536 356 : double totweight=cvals[0], tnorm = cvals[0];
537 : myvals.setValue( 1, cvals[1] );
538 : // Get the derivatives as well if we are in apply
539 356 : if( in_apply ) {
540 : // This bit gets the total weight
541 150 : double weight0 = cvals[0]; // Store the current weight
542 250 : for(unsigned j=1; j<myvessels.size(); ++j) {
543 100 : stashes[j]->retrieveSequentialValue( current, false, cvals );
544 100 : totweight *= cvals[0];
545 : }
546 : // And this bit the derivatives
547 150 : MultiValue& tmpval = stashes[0]->getTemporyMultiValue(0);
548 297 : if( tmpval.getNumberOfValues()!=myvessels[0]->getNumberOfQuantities() ||
549 147 : tmpval.getNumberOfDerivatives()!=myvessels[0]->getNumberOfDerivatives() ) {
550 3 : tmpval.resize( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
551 : }
552 150 : stashes[0]->retrieveDerivatives( stashes[0]->getTrueIndex(current), false, tmpval );
553 31590 : for(unsigned j=0; j<tmpval.getNumberActive(); ++j) {
554 31440 : unsigned jder=tmpval.getActiveIndex(j);
555 31440 : myvals.addDerivative( 1, jder, tmpval.getDerivative(1,jder) );
556 31440 : myvals.addDerivative( 0, jder, (totweight/weight0)*tmpval.getDerivative(0,jder) );
557 : }
558 150 : derbase = myvessels[0]->getNumberOfDerivatives();
559 : }
560 604 : for(unsigned i=1; i<myvessels.size(); ++i) {
561 248 : if( cvals.size()!=myvessels[i]->getNumberOfQuantities() ) {
562 8 : cvals.resize( myvessels[i]->getNumberOfQuantities() );
563 : }
564 248 : stashes[i]->retrieveSequentialValue( current, false, cvals );
565 248 : tnorm *= cvals[0];
566 248 : myvals.setValue( 1+i, cvals[1] );
567 : // Get the derivatives as well if we are in apply
568 248 : if( in_apply ) {
569 100 : MultiValue& tmpval = stashes[0]->getTemporyMultiValue(0);
570 200 : if( tmpval.getNumberOfValues()!=myvessels[0]->getNumberOfQuantities() ||
571 100 : tmpval.getNumberOfDerivatives()!=myvessels[0]->getNumberOfDerivatives() ) {
572 0 : tmpval.resize( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
573 : }
574 100 : stashes[i]->retrieveDerivatives( stashes[i]->getTrueIndex(current), false, tmpval );
575 1090 : for(unsigned j=0; j<tmpval.getNumberActive(); ++j) {
576 : unsigned jder=tmpval.getActiveIndex(j);
577 990 : myvals.addDerivative( 1+i, derbase+jder, tmpval.getDerivative(1,jder) );
578 990 : myvals.addDerivative( 0, derbase+jder, (totweight/cvals[0])*tmpval.getDerivative(0,jder) );
579 : }
580 100 : derbase += myvessels[i]->getNumberOfDerivatives();
581 : }
582 : }
583 : myvals.setValue( 0, tnorm );
584 356 : myvals.setValue( 1+myvessels.size(), ww );
585 356 : if( in_apply ) {
586 : myvals.updateDynamicList();
587 : }
588 : } else {
589 14908 : plumed_assert( !in_apply );
590 14908 : std::vector<std::unique_ptr<Value>> vv( myhist->getVectorOfValues() );
591 14908 : std::vector<double> val( getNumberOfArguments() ), der( getNumberOfArguments() );
592 : // Retrieve the location of the grid point at which we are evaluating the kernel
593 14908 : mygrid->getGridPointCoordinates( current, val );
594 14908 : if( kernel ) {
595 54492 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
596 39584 : vv[i]->set( val[i] );
597 : }
598 : // Evaluate the histogram at the relevant grid point and set the values
599 29816 : double vvh = kernel->evaluate( Tools::unique2raw(vv), der,true);
600 : myvals.setValue( 1, vvh );
601 : } else {
602 0 : plumed_merror("normalisation of vectors does not work with arguments and spherical grids");
603 : // Evalulate dot product
604 : double dot=0;
605 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
606 : dot+=val[j]*getArgument(j);
607 : der[j]=val[j];
608 : }
609 : // Von misses distribution for concentration parameter
610 : double newval = (myhist->von_misses_norm)*std::exp( (myhist->von_misses_concentration)*dot );
611 : myvals.setValue( 1, newval );
612 : // And final derivatives
613 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
614 : der[j] *= (myhist->von_misses_concentration)*newval;
615 : }
616 : }
617 : // Set the derivatives and delete the vector of values
618 54492 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
619 39584 : myvals.setDerivative( 1, i, der[i] );
620 : }
621 14908 : }
622 15404 : }
623 :
624 133 : void Histogram::apply() {
625 133 : if( !myhist->wasForced() ) {
626 : return ;
627 : }
628 20 : in_apply=true;
629 : // Run the loop to calculate the forces
630 20 : runAllTasks();
631 20 : finishAveraging();
632 : // We now need to retrieve the buffer and set the forces on the atoms
633 20 : myhist->applyForce( forcesToApply );
634 : // Now make the forces make sense for the virial
635 20 : unsigned fbase=0, tbase=0, vbase = getNumberOfDerivatives() - myvessels.size()*9;
636 200 : for(unsigned i=vbase; i<vbase+9; ++i) {
637 180 : finalForces[i]=0.0;
638 : }
639 50 : for(unsigned i=0; i<myvessels.size(); ++i) {
640 3555 : for(unsigned j=0; j<myvessels[i]->getNumberOfDerivatives()-9; ++j) {
641 3525 : finalForces[fbase + j] = forcesToApply[tbase + j];
642 : }
643 : unsigned k=0;
644 300 : for(unsigned j=myvessels[i]->getNumberOfDerivatives()-9; j<myvessels[i]->getNumberOfDerivatives(); ++j) {
645 270 : finalForces[vbase + k] += forcesToApply[tbase + j];
646 270 : k++;
647 : }
648 30 : fbase += myvessels[i]->getNumberOfDerivatives() - 9;
649 30 : tbase += myvessels[i]->getNumberOfDerivatives();
650 : }
651 : // And set the final forces on the atoms
652 20 : setForcesOnAtoms( finalForces );
653 : // Reset everything for next regular loop
654 20 : in_apply=false;
655 : }
656 :
657 : }
658 : }
|
