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 { return !in_apply; }
216 0 : bool isPeriodic() override { return false; }
217 : unsigned getNumberOfDerivatives() override;
218 : void turnOnDerivatives() override;
219 : void compute( const unsigned&, MultiValue& ) const override;
220 : void apply() override;
221 : };
222 :
223 10487 : PLUMED_REGISTER_ACTION(Histogram,"HISTOGRAM")
224 :
225 35 : void Histogram::registerKeywords( Keywords& keys ) {
226 70 : gridtools::ActionWithGrid::registerKeywords( keys ); keys.use("ARG"); keys.remove("NORMALIZATION");
227 70 : 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 70 : keys.add("optional","DATA","input data from action with vessel and compute histogram");
229 70 : keys.add("optional","VECTORS","input three dimensional vectors for computing histogram");
230 70 : keys.add("compulsory","GRID_MIN","the lower bounds for the grid");
231 70 : keys.add("compulsory","GRID_MAX","the upper bounds for the grid");
232 70 : keys.add("optional","GRID_BIN","the number of bins for the grid");
233 70 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
234 70 : 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 34 : ww(0.0),
241 34 : in_apply(false),
242 34 : 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 2 : ActionWithVessel* myv = plumed.getActionSet().selectWithLabel<ActionWithVessel*>( vlab );
249 2 : if( !myv ) error("action labelled " + vlab + " does not exist or is not an ActionWithVessel");
250 2 : 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 2 : log.printf(" for vector quantities calculated by %s \n", vlab.c_str() );
254 : } else {
255 12 : std::vector<std::string> mlab; parseVector("DATA",mlab);
256 6 : if( mlab.size()>0 ) {
257 16 : for(unsigned i=0; i<mlab.size(); ++i) {
258 10 : 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 10 : 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 10 : for(unsigned i=1; i<mlab.size(); ++i) {
267 4 : if( nvals!=myvessels[i]->getFullNumberOfTasks() ) error("mismatched number of quantities calculated by actions input to histogram");
268 : }
269 6 : log.printf(" for all base quantities calculated by %s ", myvessels[0]->getLabel().c_str() );
270 10 : 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 6 : }
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 6 : vstring += " COORDINATES=" + myvessels[0]->getLabel();
287 10 : 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 10 : for(unsigned i=1; i<myvessels.size(); ++i) {
292 4 : if( myvessels[i]->isPeriodic() ) vstring+=",T";
293 : else vstring+=",F";
294 : }
295 : } else {
296 26 : vstring += " COORDINATES=" + getPntrToArgument(0)->getName();
297 35 : 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 35 : 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 34 : 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 32 : 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 32 : } 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 4 : 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 72 : for(unsigned i=0; i<myvessels[0]->getFullNumberOfTasks(); ++i) addTaskToList(i);
329 8 : setAveragingAction( std::move(grid), true );
330 : } else if( storeThenAverage() ) {
331 7 : setAveragingAction( std::move(grid), true );
332 : } else {
333 : // Create a task list
334 9256 : for(unsigned i=0; i<mygrid->getNumberOfPoints(); ++i) addTaskToList(i);
335 19 : myhist->addOneKernelEachTimeOnly();
336 19 : setAveragingAction( std::move(grid), myhist->noDiscreteKernels() );
337 : }
338 34 : checkRead();
339 68 : }
340 :
341 8 : void Histogram::turnOnDerivatives() {
342 8 : ActionWithGrid::turnOnDerivatives();
343 : std::vector<AtomNumber> all_atoms, tmp_atoms;
344 20 : 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 : // cppcheck-suppress nullPointerRedundantCheck
348 12 : tmp_atoms = mbase->getAbsoluteIndexes();
349 482 : for(unsigned j=0; j<tmp_atoms.size(); ++j) all_atoms.push_back( tmp_atoms[j] );
350 : // Make a tempory multi value so we can avoid vector resizing
351 12 : stashes[i]->resizeTemporyMultiValues( 1 );
352 : }
353 8 : ActionAtomistic::requestAtoms( all_atoms );
354 8 : finalForces.resize( 3*all_atoms.size() + 9 );
355 8 : forcesToApply.resize( 3*all_atoms.size() + 9*myvessels.size() );
356 : // And make sure we still have the dependencies which are cleared by requestAtoms
357 20 : for(unsigned i=0; i<myvessels.size(); ++i) addDependency( myvessels[i] );
358 : // And resize the histogram so that we have a place to store the forces
359 8 : in_apply=true; mygrid->resize(); in_apply=false;
360 8 : }
361 :
362 728 : unsigned Histogram::getNumberOfDerivatives() {
363 728 : if( in_apply ) {
364 : unsigned nder=0;
365 540 : for(unsigned i=0; i<myvessels.size(); ++i) nder += myvessels[i]->getNumberOfDerivatives();
366 : return nder;
367 : }
368 506 : return getNumberOfArguments();
369 : }
370 :
371 462 : unsigned Histogram::getNumberOfQuantities() const {
372 462 : if( mvectors ) return myvessels[0]->getNumberOfQuantities();
373 416 : else if( myvessels.size()>0 ) return myvessels.size()+2;
374 294 : return ActionWithAveraging::getNumberOfQuantities();
375 : }
376 :
377 115 : void Histogram::prepareForAveraging() {
378 115 : if( myvessels.size()>0 ) {
379 36 : deactivateAllTasks(); double norm=0;
380 332 : for(unsigned i=0; i<stashes[0]->getNumberOfStoredValues(); ++i) {
381 296 : std::vector<double> cvals( myvessels[0]->getNumberOfQuantities() );
382 296 : stashes[0]->retrieveSequentialValue( i, false, cvals );
383 296 : unsigned itask=myvessels[0]->getActiveTask(i); double tnorm = cvals[0];
384 444 : for(unsigned j=1; j<myvessels.size(); ++j) {
385 148 : if( myvessels[j]->getActiveTask(i)!=itask ) error("mismatched task identities in histogram suggests histogram is meaningless");
386 148 : if( cvals.size()!=myvessels[j]->getNumberOfQuantities() ) cvals.resize( myvessels[j]->getNumberOfQuantities() );
387 148 : stashes[j]->retrieveSequentialValue( i, false, cvals ); tnorm *= cvals[0];
388 : }
389 296 : norm += tnorm; taskFlags[i]=1;
390 : }
391 36 : lockContributors();
392 : // Sort out normalization of histogram
393 36 : if( !noNormalization() ) ww = cweight / norm;
394 5 : else ww = cweight;
395 : } else if( !storeThenAverage() ) {
396 : // Now fetch the kernel and the active points
397 72 : std::vector<double> point( getNumberOfArguments() );
398 180 : for(unsigned i=0; i<point.size(); ++i) point[i]=getArgument(i);
399 72 : unsigned num_neigh; std::vector<unsigned> neighbors(1);
400 144 : kernel=myhist->getKernelAndNeighbors( point, num_neigh, neighbors );
401 :
402 72 : if( num_neigh>1 ) {
403 : // Activate relevant tasks
404 67 : deactivateAllTasks();
405 14975 : for(unsigned i=0; i<num_neigh; ++i) taskFlags[neighbors[i]]=1;
406 67 : lockContributors();
407 : } else {
408 : // This is used when we are not doing kernel density evaluation
409 5 : mygrid->addToGridElement( neighbors[0], 0, cweight );
410 : }
411 : }
412 115 : }
413 :
414 5 : void Histogram::performOperations( const bool& from_update ) { if( myvessels.size()==0 ) plumed_dbg_assert( !myhist->noDiscreteKernels() ); }
415 :
416 135 : void Histogram::finishAveraging() {
417 135 : if( myvessels.size()==0 ) kernel.reset();
418 135 : }
419 :
420 15404 : void Histogram::compute( const unsigned& current, MultiValue& myvals ) const {
421 15404 : if( mvectors ) {
422 140 : std::vector<double> cvals( myvessels[0]->getNumberOfQuantities() );
423 140 : stashes[0]->retrieveSequentialValue( current, true, cvals );
424 560 : for(unsigned i=2; i<myvessels[0]->getNumberOfQuantities(); ++i) myvals.setValue( i-1, cvals[i] );
425 140 : myvals.setValue( 0, cvals[0] ); myvals.setValue( myvessels[0]->getNumberOfQuantities() - 1, ww );
426 140 : if( in_apply ) {
427 50 : MultiValue& tmpval = stashes[0]->getTemporyMultiValue(0);
428 99 : if( tmpval.getNumberOfValues()!=myvessels[0]->getNumberOfQuantities() ||
429 49 : tmpval.getNumberOfDerivatives()!=myvessels[0]->getNumberOfDerivatives() )
430 1 : tmpval.resize( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
431 50 : stashes[0]->retrieveDerivatives( stashes[0]->getTrueIndex(current), true, tmpval );
432 350 : for(unsigned j=0; j<tmpval.getNumberActive(); ++j) {
433 300 : unsigned jder=tmpval.getActiveIndex(j); myvals.addDerivative( 0, jder, tmpval.getDerivative(0, jder) );
434 1200 : for(unsigned i=2; i<myvessels[0]->getNumberOfQuantities(); ++i) myvals.addDerivative( i-1, jder, tmpval.getDerivative(i, jder) );
435 : }
436 : myvals.updateDynamicList();
437 : }
438 15264 : } else if( myvessels.size()>0 ) {
439 356 : std::vector<double> cvals( myvessels[0]->getNumberOfQuantities() );
440 356 : stashes[0]->retrieveSequentialValue( current, false, cvals );
441 356 : unsigned derbase=0; double totweight=cvals[0], tnorm = cvals[0]; myvals.setValue( 1, cvals[1] );
442 : // Get the derivatives as well if we are in apply
443 356 : if( in_apply ) {
444 : // This bit gets the total weight
445 150 : double weight0 = cvals[0]; // Store the current weight
446 250 : for(unsigned j=1; j<myvessels.size(); ++j) {
447 100 : stashes[j]->retrieveSequentialValue( current, false, cvals ); totweight *= cvals[0];
448 : }
449 : // And this bit the derivatives
450 150 : MultiValue& tmpval = stashes[0]->getTemporyMultiValue(0);
451 297 : if( tmpval.getNumberOfValues()!=myvessels[0]->getNumberOfQuantities() ||
452 147 : tmpval.getNumberOfDerivatives()!=myvessels[0]->getNumberOfDerivatives() )
453 3 : tmpval.resize( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
454 150 : stashes[0]->retrieveDerivatives( stashes[0]->getTrueIndex(current), false, tmpval );
455 31590 : for(unsigned j=0; j<tmpval.getNumberActive(); ++j) {
456 31440 : unsigned jder=tmpval.getActiveIndex(j);
457 31440 : myvals.addDerivative( 1, jder, tmpval.getDerivative(1,jder) );
458 31440 : myvals.addDerivative( 0, jder, (totweight/weight0)*tmpval.getDerivative(0,jder) );
459 : }
460 150 : derbase = myvessels[0]->getNumberOfDerivatives();
461 : }
462 604 : for(unsigned i=1; i<myvessels.size(); ++i) {
463 248 : if( cvals.size()!=myvessels[i]->getNumberOfQuantities() ) cvals.resize( myvessels[i]->getNumberOfQuantities() );
464 248 : stashes[i]->retrieveSequentialValue( current, false, cvals );
465 248 : tnorm *= cvals[0]; myvals.setValue( 1+i, cvals[1] );
466 : // Get the derivatives as well if we are in apply
467 248 : if( in_apply ) {
468 100 : MultiValue& tmpval = stashes[0]->getTemporyMultiValue(0);
469 200 : if( tmpval.getNumberOfValues()!=myvessels[0]->getNumberOfQuantities() ||
470 100 : tmpval.getNumberOfDerivatives()!=myvessels[0]->getNumberOfDerivatives() )
471 0 : tmpval.resize( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
472 100 : stashes[i]->retrieveDerivatives( stashes[i]->getTrueIndex(current), false, tmpval );
473 1090 : for(unsigned j=0; j<tmpval.getNumberActive(); ++j) {
474 : unsigned jder=tmpval.getActiveIndex(j);
475 990 : myvals.addDerivative( 1+i, derbase+jder, tmpval.getDerivative(1,jder) );
476 990 : myvals.addDerivative( 0, derbase+jder, (totweight/cvals[0])*tmpval.getDerivative(0,jder) );
477 : }
478 100 : derbase += myvessels[i]->getNumberOfDerivatives();
479 : }
480 : }
481 356 : myvals.setValue( 0, tnorm ); myvals.setValue( 1+myvessels.size(), ww );
482 356 : if( in_apply ) myvals.updateDynamicList();
483 : } else {
484 14908 : plumed_assert( !in_apply );
485 14908 : std::vector<std::unique_ptr<Value>> vv( myhist->getVectorOfValues() );
486 14908 : std::vector<double> val( getNumberOfArguments() ), der( getNumberOfArguments() );
487 : // Retrieve the location of the grid point at which we are evaluating the kernel
488 14908 : mygrid->getGridPointCoordinates( current, val );
489 14908 : if( kernel ) {
490 54492 : for(unsigned i=0; i<getNumberOfArguments(); ++i) vv[i]->set( val[i] );
491 : // Evaluate the histogram at the relevant grid point and set the values
492 29816 : double vvh = kernel->evaluate( Tools::unique2raw(vv), der,true); myvals.setValue( 1, vvh );
493 : } else {
494 0 : plumed_merror("normalisation of vectors does not work with arguments and spherical grids");
495 : // Evalulate dot product
496 : double dot=0; for(unsigned j=0; j<getNumberOfArguments(); ++j) { dot+=val[j]*getArgument(j); der[j]=val[j]; }
497 : // Von misses distribution for concentration parameter
498 : double newval = (myhist->von_misses_norm)*std::exp( (myhist->von_misses_concentration)*dot ); myvals.setValue( 1, newval );
499 : // And final derivatives
500 : for(unsigned j=0; j<getNumberOfArguments(); ++j) der[j] *= (myhist->von_misses_concentration)*newval;
501 : }
502 : // Set the derivatives and delete the vector of values
503 54492 : for(unsigned i=0; i<getNumberOfArguments(); ++i) { myvals.setDerivative( 1, i, der[i] ); }
504 14908 : }
505 15404 : }
506 :
507 133 : void Histogram::apply() {
508 133 : if( !myhist->wasForced() ) return ;
509 20 : in_apply=true;
510 : // Run the loop to calculate the forces
511 20 : runAllTasks(); finishAveraging();
512 : // We now need to retrieve the buffer and set the forces on the atoms
513 20 : myhist->applyForce( forcesToApply );
514 : // Now make the forces make sense for the virial
515 20 : unsigned fbase=0, tbase=0, vbase = getNumberOfDerivatives() - myvessels.size()*9;
516 200 : for(unsigned i=vbase; i<vbase+9; ++i) finalForces[i]=0.0;
517 50 : for(unsigned i=0; i<myvessels.size(); ++i) {
518 3555 : for(unsigned j=0; j<myvessels[i]->getNumberOfDerivatives()-9; ++j) {
519 3525 : finalForces[fbase + j] = forcesToApply[tbase + j];
520 : }
521 : unsigned k=0;
522 300 : for(unsigned j=myvessels[i]->getNumberOfDerivatives()-9; j<myvessels[i]->getNumberOfDerivatives(); ++j) {
523 270 : finalForces[vbase + k] += forcesToApply[tbase + j]; k++;
524 : }
525 30 : fbase += myvessels[i]->getNumberOfDerivatives() - 9;
526 30 : tbase += myvessels[i]->getNumberOfDerivatives();
527 : }
528 : // And set the final forces on the atoms
529 20 : setForcesOnAtoms( finalForces );
530 : // Reset everything for next regular loop
531 20 : in_apply=false;
532 : }
533 :
534 : }
535 : }
|
