Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "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 bandwith \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 evalulated 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 belfast-1 and \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 caclulate 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 unormalized estimate of the histogram 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 continuos 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 : GRID_WFILE=histo
190 : LABEL=hh
191 : ... HISTOGRAM
192 :
193 : DUMPGRID GRID=hh FILE=histo STRIDE=100000
194 : \endplumedfile
195 :
196 : */
197 : //+ENDPLUMEDOC
198 :
199 75 : class Histogram : public gridtools::ActionWithGrid {
200 : private:
201 : double ww;
202 : bool in_apply, mvectors;
203 : KernelFunctions* kernel;
204 : std::vector<double> forcesToApply, finalForces;
205 : std::vector<vesselbase::ActionWithVessel*> myvessels;
206 : std::vector<vesselbase::StoreDataVessel*> stashes;
207 : gridtools::HistogramOnGrid* myhist;
208 : public:
209 : static void registerKeywords( Keywords& keys );
210 : explicit Histogram(const ActionOptions&ao);
211 : unsigned getNumberOfQuantities() const ;
212 : void prepareForAveraging();
213 : void performOperations( const bool& from_update );
214 : void finishAveraging();
215 109 : bool threadSafe() const { return !in_apply; }
216 0 : bool isPeriodic() { return false; }
217 : unsigned getNumberOfDerivatives();
218 : void turnOnDerivatives();
219 : void compute( const unsigned&, MultiValue& ) const ;
220 : void apply();
221 : };
222 :
223 6477 : PLUMED_REGISTER_ACTION(Histogram,"HISTOGRAM")
224 :
225 26 : void Histogram::registerKeywords( Keywords& keys ) {
226 78 : gridtools::ActionWithGrid::registerKeywords( keys ); keys.use("ARG"); keys.remove("NORMALIZATION");
227 130 : 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");
228 104 : keys.add("optional","DATA","input data from action with vessel and compute histogram");
229 104 : keys.add("optional","VECTORS","input three dimsnional vectors for computing histogram");
230 104 : keys.add("compulsory","GRID_MIN","the lower bounds for the grid");
231 104 : keys.add("compulsory","GRID_MAX","the upper bounds for the grid");
232 104 : keys.add("optional","GRID_BIN","the number of bins for the grid");
233 104 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
234 78 : keys.use("UPDATE_FROM"); keys.use("UPDATE_UNTIL");
235 26 : }
236 :
237 25 : Histogram::Histogram(const ActionOptions&ao):
238 : Action(ao),
239 : ActionWithGrid(ao),
240 : ww(0.0),
241 : in_apply(false),
242 : mvectors(false),
243 75 : kernel(NULL)
244 : {
245 : // Read in arguments
246 50 : std::string vlab; parse("VECTORS",vlab);
247 25 : if( vlab.length()>0 ) {
248 4 : ActionWithVessel* myv = plumed.getActionSet().selectWithLabel<ActionWithVessel*>( vlab );
249 2 : if( !myv ) error("action labelled " + vlab + " does not exist or is not an ActionWithVessel");
250 4 : myvessels.push_back( myv ); stashes.push_back( myv->buildDataStashes( NULL ) );
251 2 : addDependency( myv ); mvectors=true;
252 2 : if( myv->getNumberOfQuantities()!=5 ) error("can only compute histograms for three dimensional vectors");
253 4 : log.printf(" for vector quantities calculated by %s \n", vlab.c_str() );
254 : } else {
255 69 : std::vector<std::string> mlab; parseVector("DATA",mlab);
256 23 : if( mlab.size()>0 ) {
257 34 : for(unsigned i=0; i<mlab.size(); ++i) {
258 16 : ActionWithVessel* myv = plumed.getActionSet().selectWithLabel<ActionWithVessel*>( mlab[i] );
259 8 : if( !myv ) error("action labelled " + mlab[i] + " does not exist or is not an ActionWithVessel");
260 16 : myvessels.push_back( myv ); stashes.push_back( myv->buildDataStashes( NULL ) );
261 : // log.printf(" for all base quantities calculated by %s \n",myvessel->getLabel().c_str() );
262 : // Add the dependency
263 8 : addDependency( myv );
264 : }
265 5 : unsigned nvals = myvessels[0]->getFullNumberOfTasks();
266 19 : for(unsigned i=1; i<mlab.size(); ++i) {
267 6 : if( nvals!=myvessels[i]->getFullNumberOfTasks() ) error("mismatched number of quantities calculated by actions input to histogram");
268 : }
269 15 : log.printf(" for all base quantities calculated by %s ", myvessels[0]->getLabel().c_str() );
270 22 : for(unsigned i=1; i<mlab.size(); ++i) log.printf(", %s \n", myvessels[i]->getLabel().c_str() );
271 5 : log.printf("\n");
272 : } else {
273 36 : std::vector<Value*> arg; parseArgumentList("ARG",arg);
274 18 : if(!arg.empty()) {
275 18 : log.printf(" with arguments");
276 144 : for(unsigned i=0; i<arg.size(); i++) log.printf(" %s",arg[i]->getName().c_str());
277 18 : log.printf("\n");
278 : // Retrieve the bias acting and make sure we request this also
279 18 : std::vector<Value*> bias( ActionWithArguments::getArguments() );
280 51 : for(unsigned i=0; i<bias.size(); ++i) arg.push_back( bias[i] );
281 18 : requestArguments(arg);
282 : }
283 : }
284 : }
285 :
286 : // Read stuff for grid
287 : unsigned narg = getNumberOfArguments();
288 25 : if( myvessels.size()>0 ) narg=myvessels.size();
289 : // Input of name and labels
290 25 : std::string vstring="COMPONENTS=" + getLabel();
291 25 : if( mvectors ) {
292 : vstring += " COORDINATES=x,y,z PBC=F,F,F";
293 23 : } else if( myvessels.size()>0 ) {
294 15 : vstring += " COORDINATES=" + myvessels[0]->getLabel();
295 25 : for(unsigned i=1; i<myvessels.size(); ++i) vstring +="," + myvessels[i]->getLabel();
296 : // Input for PBC
297 5 : if( myvessels[0]->isPeriodic() ) vstring+=" PBC=T";
298 : else vstring+=" PBC=F";
299 19 : for(unsigned i=1; i<myvessels.size(); ++i) {
300 3 : if( myvessels[i]->isPeriodic() ) vstring+=",T";
301 : else vstring+=",F";
302 : }
303 : } else {
304 36 : vstring += " COORDINATES=" + getPntrToArgument(0)->getName();
305 45 : for(unsigned i=1; i<getNumberOfArguments(); ++i) vstring += "," + getPntrToArgument(i)->getName();
306 : // Input for PBC
307 18 : if( getPntrToArgument(0)->isPeriodic() ) vstring+=" PBC=T";
308 : else vstring+=" PBC=F";
309 36 : for(unsigned i=1; i<getNumberOfArguments(); ++i) {
310 9 : if( getPntrToArgument(i)->isPeriodic() ) vstring+=",T";
311 : else vstring+=",F";
312 : }
313 : }
314 : // And create the grid
315 50 : createGrid( "histogram", vstring );
316 50 : if( mygrid->getType()=="flat" ) {
317 23 : if( mvectors ) error("computing histogram for three dimensional vectors but grid is not of fibonacci type - use CONCENTRATION");
318 46 : std::vector<std::string> gmin( narg ), gmax( narg );
319 69 : parseVector("GRID_MIN",gmin); parseVector("GRID_MAX",gmax);
320 46 : std::vector<unsigned> nbin; parseVector("GRID_BIN",nbin);
321 46 : std::vector<double> gspacing; parseVector("GRID_SPACING",gspacing);
322 23 : if( nbin.size()!=narg && gspacing.size()!=narg ) {
323 0 : error("GRID_BIN or GRID_SPACING must be set");
324 : }
325 23 : mygrid->setBounds( gmin, gmax, nbin, gspacing );
326 : } else {
327 4 : std::vector<unsigned> nbin; parseVector("GRID_BIN",nbin);
328 2 : if( nbin.size()!=1 ) error("should only be one index for number of bins with spherical grid");
329 6 : if( mygrid->getType()=="fibonacci" ) mygrid->setupFibonacciGrid( nbin[0] );
330 : }
331 25 : myhist = dynamic_cast<gridtools::HistogramOnGrid*>( mygrid );
332 25 : plumed_assert( myhist );
333 25 : if( myvessels.size()>0 ) {
334 : // Create a task list
335 129 : for(unsigned i=0; i<myvessels[0]->getFullNumberOfTasks(); ++i) addTaskToList(i);
336 7 : setAveragingAction( mygrid, true );
337 : } else {
338 : // Create a task list
339 9172 : for(unsigned i=0; i<mygrid->getNumberOfPoints(); ++i) addTaskToList(i);
340 18 : myhist->addOneKernelEachTimeOnly();
341 18 : setAveragingAction( mygrid, myhist->noDiscreteKernels() );
342 : }
343 25 : checkRead();
344 25 : }
345 :
346 8 : void Histogram::turnOnDerivatives() {
347 8 : ActionWithGrid::turnOnDerivatives();
348 : std::vector<AtomNumber> all_atoms, tmp_atoms;
349 52 : for(unsigned i=0; i<myvessels.size(); ++i) {
350 12 : multicolvar::MultiColvarBase* mbase=dynamic_cast<multicolvar::MultiColvarBase*>( myvessels[i] );
351 12 : if( !mbase ) error("do not know how to get histogram derivatives for actions of type " + myvessels[i]->getName() );
352 12 : tmp_atoms = mbase->getAbsoluteIndexes();
353 1434 : for(unsigned j=0; j<tmp_atoms.size(); ++j) all_atoms.push_back( tmp_atoms[j] );
354 : // Make a tempory multi value so we can avoid vector resizing
355 12 : stashes[i]->resizeTemporyMultiValues( 1 );
356 : }
357 8 : ActionAtomistic::requestAtoms( all_atoms );
358 16 : finalForces.resize( 3*all_atoms.size() + 9 );
359 24 : forcesToApply.resize( 3*all_atoms.size() + 9*myvessels.size() );
360 : // And make sure we still have the dependencies which are cleared by requestAtoms
361 52 : for(unsigned i=0; i<myvessels.size(); ++i) addDependency( myvessels[i] );
362 : // And resize the histogram so that we have a place to store the forces
363 8 : in_apply=true; mygrid->resize(); in_apply=false;
364 8 : }
365 :
366 660 : unsigned Histogram::getNumberOfDerivatives() {
367 660 : if( in_apply ) {
368 : unsigned nder=0;
369 1398 : for(unsigned i=0; i<myvessels.size(); ++i) nder += myvessels[i]->getNumberOfDerivatives();
370 : return nder;
371 : }
372 438 : return getNumberOfArguments();
373 : }
374 :
375 412 : unsigned Histogram::getNumberOfQuantities() const {
376 458 : if( mvectors ) return myvessels[0]->getNumberOfQuantities();
377 366 : else if( myvessels.size()>0 ) return myvessels.size()+2;
378 : return 2;
379 : }
380 :
381 102 : void Histogram::prepareForAveraging() {
382 102 : if( myvessels.size()>0 ) {
383 32 : deactivateAllTasks(); double norm=0;
384 576 : for(unsigned i=0; i<stashes[0]->getNumberOfStoredValues(); ++i) {
385 256 : std::vector<double> cvals( myvessels[0]->getNumberOfQuantities() );
386 256 : stashes[0]->retrieveSequentialValue( i, false, cvals );
387 512 : unsigned itask=myvessels[0]->getActiveTask(i); double tnorm = cvals[0];
388 836 : for(unsigned j=1; j<myvessels.size(); ++j) {
389 216 : if( myvessels[j]->getActiveTask(i)!=itask ) error("mismatched task identities in histogram suggests histogram is meaningless");
390 116 : if( cvals.size()!=myvessels[j]->getNumberOfQuantities() ) cvals.resize( myvessels[j]->getNumberOfQuantities() );
391 216 : stashes[j]->retrieveSequentialValue( i, false, cvals ); tnorm *= cvals[0];
392 : }
393 512 : norm += tnorm; taskFlags[i]=1;
394 : }
395 32 : lockContributors();
396 : // Sort out normalization of histogram
397 32 : if( !noNormalization() ) ww = cweight / norm;
398 5 : else ww = cweight;
399 : } else {
400 : // Now fetch the kernel and the active points
401 70 : std::vector<double> point( getNumberOfArguments() );
402 458 : for(unsigned i=0; i<point.size(); ++i) point[i]=getArgument(i);
403 70 : unsigned num_neigh; std::vector<unsigned> neighbors(1);
404 70 : kernel = myhist->getKernelAndNeighbors( point, num_neigh, neighbors );
405 :
406 70 : if( num_neigh>1 ) {
407 : // Activate relevant tasks
408 65 : deactivateAllTasks();
409 44204 : for(unsigned i=0; i<num_neigh; ++i) taskFlags[neighbors[i]]=1;
410 65 : lockContributors();
411 : } else {
412 : // This is used when we are not doing kernel density evaluation
413 10 : mygrid->addToGridElement( neighbors[0], 0, cweight );
414 : }
415 : }
416 102 : }
417 :
418 5 : void Histogram::performOperations( const bool& from_update ) { if( myvessels.size()==0 ) plumed_dbg_assert( !myhist->noDiscreteKernels() ); }
419 :
420 122 : void Histogram::finishAveraging() {
421 122 : if( myvessels.size()==0 ) delete kernel;
422 122 : }
423 :
424 15169 : void Histogram::compute( const unsigned& current, MultiValue& myvals ) const {
425 15169 : if( mvectors ) {
426 140 : std::vector<double> cvals( myvessels[0]->getNumberOfQuantities() );
427 140 : stashes[0]->retrieveSequentialValue( current, true, cvals );
428 1400 : for(unsigned i=2; i<myvessels[0]->getNumberOfQuantities(); ++i) myvals.setValue( i-1, cvals[i] );
429 140 : myvals.setValue( 0, cvals[0] ); myvals.setValue( myvessels[0]->getNumberOfQuantities() - 1, ww );
430 140 : if( in_apply ) {
431 50 : MultiValue& tmpval = stashes[0]->getTemporyMultiValue(0);
432 99 : if( tmpval.getNumberOfValues()!=myvessels[0]->getNumberOfQuantities() ||
433 49 : tmpval.getNumberOfDerivatives()!=myvessels[0]->getNumberOfDerivatives() )
434 2 : tmpval.resize( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
435 100 : stashes[0]->retrieveDerivatives( stashes[0]->getTrueIndex(current), true, tmpval );
436 650 : for(unsigned j=0; j<tmpval.getNumberActive(); ++j) {
437 600 : unsigned jder=tmpval.getActiveIndex(j); myvals.addDerivative( 0, jder, tmpval.getDerivative(0, jder) );
438 2100 : for(unsigned i=2; i<myvessels[0]->getNumberOfQuantities(); ++i) myvals.addDerivative( i-1, jder, tmpval.getDerivative(i, jder) );
439 : }
440 : myvals.updateDynamicList();
441 : }
442 15029 : } else if( myvessels.size()>0 ) {
443 316 : std::vector<double> cvals( myvessels[0]->getNumberOfQuantities() );
444 316 : stashes[0]->retrieveSequentialValue( current, false, cvals );
445 316 : unsigned derbase; double totweight=cvals[0], tnorm = cvals[0]; myvals.setValue( 1, cvals[1] );
446 : // Get the derivatives as well if we are in apply
447 316 : if( in_apply ) {
448 : // This bit gets the total weight
449 150 : double weight0 = cvals[0]; // Store the current weight
450 600 : for(unsigned j=1; j<myvessels.size(); ++j) {
451 200 : stashes[j]->retrieveSequentialValue( current, false, cvals ); totweight *= cvals[0];
452 : }
453 : // And this bit the derivatives
454 150 : MultiValue& tmpval = stashes[0]->getTemporyMultiValue(0);
455 297 : if( tmpval.getNumberOfValues()!=myvessels[0]->getNumberOfQuantities() ||
456 147 : tmpval.getNumberOfDerivatives()!=myvessels[0]->getNumberOfDerivatives() )
457 6 : tmpval.resize( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
458 300 : stashes[0]->retrieveDerivatives( stashes[0]->getTrueIndex(current), false, tmpval );
459 63030 : for(unsigned j=0; j<tmpval.getNumberActive(); ++j) {
460 31440 : unsigned jder=tmpval.getActiveIndex(j);
461 31440 : myvals.addDerivative( 1, jder, tmpval.getDerivative(1,jder) );
462 62880 : myvals.addDerivative( 0, jder, (totweight/weight0)*tmpval.getDerivative(0,jder) );
463 : }
464 150 : derbase = myvessels[0]->getNumberOfDerivatives();
465 : }
466 1048 : for(unsigned i=1; i<myvessels.size(); ++i) {
467 216 : if( cvals.size()!=myvessels[i]->getNumberOfQuantities() ) cvals.resize( myvessels[i]->getNumberOfQuantities() );
468 208 : stashes[i]->retrieveSequentialValue( current, false, cvals );
469 208 : tnorm *= cvals[0]; myvals.setValue( 1+i, cvals[1] );
470 : // Get the derivatives as well if we are in apply
471 208 : if( in_apply ) {
472 100 : MultiValue& tmpval = stashes[0]->getTemporyMultiValue(0);
473 200 : if( tmpval.getNumberOfValues()!=myvessels[0]->getNumberOfQuantities() ||
474 100 : tmpval.getNumberOfDerivatives()!=myvessels[0]->getNumberOfDerivatives() )
475 0 : tmpval.resize( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
476 200 : stashes[i]->retrieveDerivatives( stashes[i]->getTrueIndex(current), false, tmpval );
477 2080 : for(unsigned j=0; j<tmpval.getNumberActive(); ++j) {
478 : unsigned jder=tmpval.getActiveIndex(j);
479 990 : myvals.addDerivative( 1+i, derbase+jder, tmpval.getDerivative(1,jder) );
480 1980 : myvals.addDerivative( 0, derbase+jder, (totweight/cvals[0])*tmpval.getDerivative(0,jder) );
481 : }
482 100 : derbase += myvessels[i]->getNumberOfDerivatives();
483 : }
484 : }
485 316 : myvals.setValue( 0, tnorm ); myvals.setValue( 1+myvessels.size(), ww );
486 316 : if( in_apply ) myvals.updateDynamicList();
487 : } else {
488 14713 : plumed_assert( !in_apply );
489 14713 : std::vector<Value*> vv( myhist->getVectorOfValues() );
490 14713 : std::vector<double> val( getNumberOfArguments() ), der( getNumberOfArguments() );
491 : // Retrieve the location of the grid point at which we are evaluating the kernel
492 14713 : mygrid->getGridPointCoordinates( current, val );
493 14713 : if( kernel ) {
494 172269 : for(unsigned i=0; i<getNumberOfArguments(); ++i) vv[i]->set( val[i] );
495 : // Evaluate the histogram at the relevant grid point and set the values
496 14713 : double vvh = kernel->evaluate( vv, der,true); myvals.setValue( 1, vvh );
497 : } else {
498 0 : plumed_merror("normalisation of vectors does not work with arguments and spherical grids");
499 : // Evalulate dot product
500 : double dot=0; for(unsigned j=0; j<getNumberOfArguments(); ++j) { dot+=val[j]*getArgument(j); der[j]=val[j]; }
501 : // Von misses distribution for concentration parameter
502 : double newval = (myhist->von_misses_norm)*exp( (myhist->von_misses_concentration)*dot ); myvals.setValue( 1, newval );
503 : // And final derivatives
504 : for(unsigned j=0; j<getNumberOfArguments(); ++j) der[j] *= (myhist->von_misses_concentration)*newval;
505 : }
506 : // Set the derivatives and delete the vector of values
507 93491 : for(unsigned i=0; i<getNumberOfArguments(); ++i) { myvals.setDerivative( 1, i, der[i] ); delete vv[i]; }
508 : }
509 15169 : }
510 :
511 125 : void Histogram::apply() {
512 125 : if( !myhist->wasForced() ) return ;
513 20 : in_apply=true;
514 : // Run the loop to calculate the forces
515 20 : runAllTasks(); finishAveraging();
516 : // We now need to retrieve the buffer and set the forces on the atoms
517 20 : myhist->applyForce( forcesToApply );
518 : // Now make the forces make sense for the virial
519 40 : unsigned fbase=0, tbase=0, vbase = getNumberOfDerivatives() - myvessels.size()*9;
520 380 : for(unsigned i=vbase; i<vbase+9; ++i) finalForces[i]=0.0;
521 130 : for(unsigned i=0; i<myvessels.size(); ++i) {
522 7080 : for(unsigned j=0; j<myvessels[i]->getNumberOfDerivatives()-9; ++j) {
523 10575 : finalForces[fbase + j] = forcesToApply[tbase + j];
524 : }
525 : unsigned k=0;
526 600 : for(unsigned j=myvessels[i]->getNumberOfDerivatives()-9; j<myvessels[i]->getNumberOfDerivatives(); ++j) {
527 810 : finalForces[vbase + k] += forcesToApply[tbase + j]; k++;
528 : }
529 30 : fbase += myvessels[i]->getNumberOfDerivatives() - 9;
530 30 : tbase += myvessels[i]->getNumberOfDerivatives();
531 : }
532 : // And set the final forces on the atoms
533 20 : setForcesOnAtoms( finalForces );
534 : // Reset everything for next regular loop
535 20 : in_apply=false;
536 : }
537 :
538 : }
539 4839 : }
|
