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