Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 : #ifndef __PLUMED_core_Value_h
23 : #define __PLUMED_core_Value_h
24 :
25 : #include <vector>
26 : #include <string>
27 : #include <map>
28 : #include "tools/Exception.h"
29 : #include "tools/Tools.h"
30 : #include "tools/AtomNumber.h"
31 : #include "tools/Vector.h"
32 :
33 : namespace PLMD {
34 :
35 : class OFile;
36 : class ActionWithValue;
37 : class ActionAtomistic;
38 :
39 : /// \ingroup TOOLBOX
40 : /// A class for holding the value of a function together with its derivatives.
41 : /// Typically, an object of type PLMD::ActionWithValue will contain one
42 : /// object of type PLUMD::Value that will be named after the label. If the
43 : /// PLMD::ActionWithValue is part of a class that calculates multiple components
44 : /// then the class will contain multiple that will be called label.component-name
45 : /// This class is used to pass information between different PLMD::Action
46 : /// objects. However, if you find a use for a tempory PLMD::Value in some method
47 : /// you are implementing please feel free to use it.
48 : class Value {
49 : friend class ActionWithValue;
50 : friend class ActionWithVector;
51 : friend class ActionAtomistic;
52 : friend class ActionWithArguments;
53 : friend class ActionWithVirtualAtom;
54 : friend class DomainDecomposition;
55 : template<typename T>
56 : friend class DataPassingObjectTyped;
57 : private:
58 : /// The action in which this quantity is calculated
59 : ActionWithValue* action;
60 : /// Had the value been set
61 : bool value_set;
62 : /// The value of the quantity
63 : std::vector<double> data;
64 : /// The force acting on this quantity
65 : std::vector<double> inputForce;
66 : /// A flag telling us we have a force acting on this quantity
67 : bool hasForce;
68 : /// The way this value is used in the code
69 : /// normal = regular value that is determined during calculate
70 : /// constant = constnt value that is determined during startup and that doesn't change during simulation
71 : /// average = value that is averaged/collected over multiple steps of trajectory
72 : /// calcFromAverage = value that is calculated from an average value
73 : enum {normal,constant,average,calcFromAverage} valtype=normal;
74 : /// This is used by ActionWithValue to set the valtype
75 : void setValType( const std::string& vtype );
76 : /// This is used by ActionWithValue to determine if we need to calculate on update
77 : bool calculateOnUpdate() const ;
78 : /// The derivatives of the quantity stored in value
79 : std::map<AtomNumber,Vector> gradients;
80 : /// The name of this quantiy
81 : std::string name;
82 : /// Are we storing the data for this value if it is vector or matrix
83 : bool storedata;
84 : /// What is the shape of the value (0 dimensional=scalar, n dimensional with derivatives=grid, 1 dimensional no derivatives=vector, 2 dimensional no derivatives=matrix)
85 : std::vector<unsigned> shape;
86 : /// Does this quanity have derivatives
87 : bool hasDeriv;
88 : /// Variables for storing data
89 : unsigned bufstart, streampos, matpos, ngrid_der, ncols, book_start;
90 : /// If we are storing a matrix is it symmetric?
91 : bool symmetric;
92 : /// This is a bookeeping array that holds the non-zero elements of the "sparse" matrix
93 : std::vector<unsigned> matrix_bookeeping;
94 : /// Is this quantity periodic
95 : enum {unset,periodic,notperiodic} periodicity;
96 : /// Various quantities that describe the domain of this value
97 : std::string str_min, str_max;
98 : double min,max;
99 : double max_minus_min;
100 : double inv_max_minus_min;
101 : /// Is the derivative of this quantity zero when the value is zero
102 : bool derivativeIsZeroWhenValueIsZero;
103 : /// Complete the setup of the periodicity
104 : void setupPeriodicity();
105 : // bring value within PBCs
106 : void applyPeriodicity( const unsigned& ival );
107 : public:
108 : /// A constructor that can be used to make Vectors of values
109 : Value();
110 : /// A constructor that can be used to make Vectors of named values
111 : explicit Value(const std::string& name);
112 : /// A constructor that is used throughout the code to setup the value poiters
113 : Value(ActionWithValue* av, const std::string& name, const bool withderiv,const std::vector<unsigned>&ss=std::vector<unsigned>());
114 : /// Set the shape of the Value
115 : void setShape( const std::vector<unsigned>&ss );
116 : /// Set the value of the function
117 : void set(double);
118 : /// Set the value of the stored data
119 : void set(const std::size_t& n, const double& v );
120 : /// Add something to the value of the function
121 : void add(double);
122 : /// Add something to the ith element of the data array
123 : void add(const std::size_t& n, const double& v );
124 : /// Get the location of this element of in the store
125 : std::size_t getIndexInStore( const std::size_t& ival ) const ;
126 : /// Get the value of the function
127 6852532 : double get( const std::size_t& ival=0, const bool trueind=true ) const;
128 : /// Find out if the value has been set
129 : bool valueHasBeenSet() const;
130 : /// Check if the value is periodic
131 : bool isPeriodic() const;
132 : /// Set the function not periodic
133 : void setNotPeriodic();
134 : /// Set the domain of the function
135 : void setDomain(const std::string&, const std::string&);
136 : /// Get the domain of the quantity
137 : void getDomain(std::string&,std::string&) const;
138 : /// Get the domain of the quantity
139 : void getDomain(double&,double&) const;
140 : /// Get the name of the quantity
141 : const std::string& getName() const;
142 : /// Check whether or not this particular quantity has derivatives
143 : bool hasDerivatives()const;
144 : /// Get the number of derivatives that this particular value has
145 : unsigned getNumberOfDerivatives() const;
146 : /// Set the number of derivatives
147 : void resizeDerivatives(int n);
148 : /// Set all the derivatives to zero
149 : void clearDerivatives( const bool force=false );
150 : /// Add some derivative to the ith component of the derivatives array
151 : void addDerivative(unsigned i,double d);
152 : /// Set the value of the ith component of the derivatives array
153 : void setDerivative(unsigned i, double d);
154 : /// Get the derivative with respect to component n
155 : double getDerivative(const unsigned n) const;
156 : /// Clear the input force on the variable
157 : void clearInputForce();
158 : /// Special method for clearing forces on variables used by DataPassingObject
159 : void clearInputForce( const std::vector<AtomNumber>& index );
160 : /// Add some force on this value
161 : void addForce(double f);
162 : /// Add some force on the ival th component of this value
163 : void addForce( const std::size_t& ival, double f, const bool trueind=true );
164 : /// Get the value of the force on this colvar
165 : double getForce( const std::size_t& ival=0 ) const ;
166 : /// Apply the forces to the derivatives using the chain rule (if there are no forces this routine returns false)
167 : bool applyForce( std::vector<double>& forces ) const ;
168 : /// Calculate the difference between the instantaneous value of the function and some other point: other_point-inst_val
169 : double difference(double)const;
170 : /// Calculate the difference between two values of this function: d2 -d1
171 : double difference(double d1,double d2)const;
172 : /// This returns the pointer to the action where this value is calculated
173 : ActionWithValue* getPntrToAction();
174 : /// Bring back one value into the correct pbc if needed, else give back the value
175 : double bringBackInPbc(double d1)const;
176 : /// Get the difference between max and minimum of domain
177 : double getMaxMinusMin()const;
178 : /// This sets up the gradients
179 : void setGradients( ActionAtomistic* aa, unsigned& start );
180 : /// This passes gradients from one action to another
181 : void passGradients( const double& der, std::map<AtomNumber,Vector>& g ) const ;
182 : static double projection(const Value&,const Value&);
183 : /// Get the rank of the object that is contained in this value
184 : unsigned getRank() const ;
185 : /// Get the shape of the object that is contained in this value
186 : const std::vector<unsigned>& getShape() const ;
187 : /// This turns on storing of vectors/matrices
188 : void buildDataStore( const bool forprint=false );
189 : /// Reshape the storage for sparse matrices
190 : void reshapeMatrixStore( const unsigned& n );
191 : /// Set the symmetric flag equal true for this matrix
192 : void setSymmetric( const bool& sym );
193 : /// Get the total number of scalars that are stored here
194 : unsigned getNumberOfValues() const ;
195 : /// Get the number of values that are actually stored here once sparse matrices are taken into account
196 : unsigned getNumberOfStoredValues() const ;
197 : /// Get the number of threads to use when assigning this value
198 : unsigned getGoodNumThreads( const unsigned& j, const unsigned& k ) const ;
199 : /// These are used for passing around the data in this value when we are doing replica exchange
200 : void writeBinary(std::ostream&o) const ;
201 : void readBinary(std::istream&i);
202 : /// These are used for making constant values
203 : bool isConstant() const ;
204 : void setConstant();
205 : /// Check if forces have been added on this value
206 : bool forcesWereAdded() const ;
207 : /// Set a bool that tells us if the derivative is zero when the value is zero true
208 : void setDerivativeIsZeroWhenValueIsZero();
209 : /// Return a bool that tells us if the derivative is zero when the value is zero
210 : bool isDerivativeZeroWhenValueIsZero() const ;
211 : ///
212 : unsigned getPositionInStream() const ;
213 : /// This stuff handles where to look for the start of the row that contains the row of the matrix
214 : void setPositionInMatrixStash( const unsigned& p );
215 : unsigned getPositionInMatrixStash() const ;
216 : /// This stuff handles where to keep the bookeeping stuff for storing the sparse matrix
217 : void setMatrixBookeepingStart( const unsigned& b );
218 : unsigned getMatrixBookeepingStart() const ;
219 : /// Convert the input index to its corresponding indices
220 : void convertIndexToindices(const std::size_t& index, std::vector<unsigned>& indices ) const ;
221 : /// Print out all the values in this Value
222 : void print( OFile& ofile ) const ;
223 : /// Are we to ignore the stored value
224 : bool ignoreStoredValue(const std::string& n) const ;
225 : /// Set a matrix element to be non zero
226 : void setMatrixBookeepingElement( const unsigned& i, const unsigned& n );
227 : ///
228 : unsigned getRowLength( const unsigned& irow ) const ;
229 : ///
230 : unsigned getRowIndex( const unsigned& irow, const unsigned& jind ) const ;
231 : /// Are we storing this value
232 : bool valueIsStored() const ;
233 : ///
234 : unsigned getNumberOfColumns() const ;
235 : ///
236 : bool isSymmetric() const ;
237 : /// Retrieve the non-zero edges in a matrix
238 : void retrieveEdgeList( unsigned& nedge, std::vector<std::pair<unsigned,unsigned> >& active, std::vector<double>& elems );
239 : /// Get the number of derivatives that the grid has
240 : unsigned getNumberOfGridDerivatives() const ;
241 : /// get the derivative of a grid at a point n with resepct to argument j
242 : double getGridDerivative(const unsigned& n, const unsigned& j ) const ;
243 : /// Add the derivatives of the grid to the corner
244 : void addGridDerivatives( const unsigned& n, const unsigned& j, const double& val );
245 : ///
246 : void setGridDerivatives( const unsigned& n, const unsigned& j, const double& val );
247 : /// Add another value to the end of the data vector held by this value. This is used in COLLECT
248 : void push_back( const double& val );
249 : /// Get the type of value that is stored here
250 : std::string getValueType() const ;
251 : };
252 :
253 : inline
254 56843038 : void Value::applyPeriodicity(const unsigned& ival) {
255 56843038 : if(periodicity==periodic) {
256 2457001 : data[ival]=min+difference(min,data[ival]);
257 2457001 : if(data[ival]<min) {
258 828453 : data[ival]+=max_minus_min;
259 : }
260 : }
261 56843038 : }
262 :
263 : inline
264 : void Value::set(double v) {
265 5325177 : value_set=true;
266 5547340 : data[0]=v;
267 5507186 : applyPeriodicity(0);
268 66605 : }
269 :
270 : inline
271 2144 : void Value::add(double v) {
272 2144 : value_set=true;
273 2144 : data[0]+=v;
274 2144 : applyPeriodicity(0);
275 2144 : }
276 :
277 : inline
278 46452643 : void Value::add(const std::size_t& n, const double& v ) {
279 46452643 : value_set=true;
280 46452643 : data[n]+=v;
281 46452643 : applyPeriodicity(n);
282 46452643 : }
283 :
284 : inline
285 : bool Value::valueHasBeenSet() const {
286 235211657 : return value_set;
287 : }
288 :
289 : inline
290 : const std::string& Value::getName()const {
291 5392349 : return name;
292 : }
293 :
294 : inline
295 16372320 : unsigned Value::getNumberOfDerivatives() const {
296 16372320 : plumed_massert(hasDeriv,"the derivatives array for this value has zero size");
297 16372320 : if( shape.size()>0 ) {
298 0 : return shape.size();
299 : }
300 16372320 : return data.size() - 1;
301 : }
302 :
303 : inline
304 : double Value::getDerivative(const unsigned n) const {
305 : plumed_dbg_massert(n<getNumberOfDerivatives(),"you are asking for a derivative that is out of bounds");
306 19334473 : return data[1+n];
307 : }
308 :
309 : inline
310 : bool Value::hasDerivatives() const {
311 225379799 : return hasDeriv;
312 : }
313 :
314 : inline
315 4448719 : void Value::resizeDerivatives(int n) {
316 4448719 : if( shape.size()>0 ) {
317 : return;
318 : }
319 3551582 : if(hasDeriv) {
320 1866457 : data.resize(1+n);
321 : }
322 : }
323 :
324 : inline
325 : void Value::addDerivative(unsigned i,double d) {
326 : plumed_dbg_massert(i<getNumberOfDerivatives(),"derivative is out of bounds");
327 20363344 : data[1+i]+=d;
328 : }
329 :
330 : inline
331 : void Value::setDerivative(unsigned i, double d) {
332 : plumed_dbg_massert(i<getNumberOfDerivatives(),"derivative is out of bounds");
333 17727802 : data[1+i]=d;
334 : }
335 :
336 : inline
337 : void Value::clearInputForce() {
338 3284108 : if( !hasForce ) {
339 : return;
340 : }
341 559654 : hasForce=false;
342 : std::fill(inputForce.begin(),inputForce.end(),0);
343 : }
344 :
345 : inline
346 116940 : void Value::clearInputForce( const std::vector<AtomNumber>& index ) {
347 116940 : if( !hasForce ) {
348 : return;
349 : }
350 62736 : hasForce=false;
351 2014386 : for(const auto & p : index) {
352 1951650 : inputForce[p.index()]=0;
353 : }
354 : }
355 :
356 : inline
357 2825523 : void Value::clearDerivatives( const bool force ) {
358 2825523 : if( !force && (valtype==constant || valtype==average) ) {
359 : return;
360 : }
361 :
362 2805420 : value_set=false;
363 2805420 : if( data.size()>1 ) {
364 : std::fill(data.begin()+1, data.end(), 0);
365 : }
366 : }
367 :
368 : inline
369 : void Value::addForce(double f) {
370 181575 : hasForce=true;
371 181575 : inputForce[0]+=f;
372 2623 : }
373 :
374 : inline
375 : bool Value::forcesWereAdded() const {
376 19702072 : return hasForce;
377 : }
378 :
379 : inline
380 : double Value::getForce( const std::size_t& ival ) const {
381 : plumed_dbg_assert( ival<inputForce.size() );
382 24565385 : return inputForce[ival];
383 : }
384 :
385 : /// d2-d1
386 : inline
387 135752510 : double Value::difference(double d1,double d2)const {
388 135752510 : if(periodicity==notperiodic) {
389 33853863 : return d2-d1;
390 101898647 : } else if(periodicity==periodic) {
391 101898647 : double s=(d2-d1)*inv_max_minus_min;
392 : // remember: pbc brings the difference in a range of -0.5:0.5
393 101898647 : s=Tools::pbc(s);
394 101898647 : return s*max_minus_min;
395 : } else {
396 0 : plumed_merror("periodicity should be set to compute differences");
397 : }
398 : }
399 :
400 : inline
401 3173 : double Value::bringBackInPbc(double d1)const {
402 3173 : return min+max_minus_min/2.+difference(min+max_minus_min/2., d1);
403 : }
404 :
405 : inline
406 3989526 : double Value::difference(double d)const {
407 3989526 : return difference(get(),d);
408 : }
409 :
410 : inline
411 : double Value::getMaxMinusMin()const {
412 : plumed_dbg_assert( periodicity==periodic );
413 0 : return max_minus_min;
414 : }
415 :
416 : inline
417 : unsigned Value::getRank() const {
418 708331592 : return shape.size();
419 : }
420 :
421 : inline
422 : const std::vector<unsigned>& Value::getShape() const {
423 485900 : return shape;
424 : }
425 :
426 : inline
427 45599091 : unsigned Value::getNumberOfValues() const {
428 : unsigned size=1;
429 55624521 : for(unsigned i=0; i<shape.size(); ++i) {
430 10025430 : size *= shape[i];
431 : }
432 45599091 : return size;
433 : }
434 :
435 : inline
436 2713987 : unsigned Value::getNumberOfStoredValues() const {
437 2713987 : if( getRank()==2 && !hasDeriv ) {
438 2264371 : return shape[0]*ncols;
439 : }
440 449616 : return getNumberOfValues();
441 : }
442 :
443 : inline
444 : bool Value::isConstant() const {
445 400198 : return valtype==constant;
446 : }
447 :
448 : inline
449 : void Value::setDerivativeIsZeroWhenValueIsZero() {
450 1205 : derivativeIsZeroWhenValueIsZero=true;
451 15 : }
452 :
453 : inline
454 : bool Value::isDerivativeZeroWhenValueIsZero() const {
455 1798 : return derivativeIsZeroWhenValueIsZero;
456 : }
457 :
458 : inline
459 : unsigned Value::getPositionInStream() const {
460 1061459775 : return streampos;
461 : }
462 :
463 : inline
464 : unsigned Value::getPositionInMatrixStash() const {
465 28657011 : return matpos;
466 : }
467 :
468 : inline
469 : void Value::setMatrixBookeepingStart( const unsigned& b ) {
470 7794 : book_start = b;
471 : }
472 :
473 : inline
474 : unsigned Value::getMatrixBookeepingStart() const {
475 14169812 : return book_start;
476 : }
477 :
478 : inline
479 : void Value::setMatrixBookeepingElement( const unsigned& i, const unsigned& n ) {
480 : plumed_dbg_assert( i<matrix_bookeeping.size() );
481 5060232 : matrix_bookeeping[i]=n;
482 : }
483 :
484 : inline
485 : bool Value::valueIsStored() const {
486 117898873 : return storedata;
487 : }
488 :
489 : inline
490 : unsigned Value::getRowLength( const unsigned& irow ) const {
491 : plumed_dbg_assert( (1+ncols)*irow<matrix_bookeeping.size() );
492 34407153 : return matrix_bookeeping[(1+ncols)*irow];
493 : }
494 :
495 : inline
496 : unsigned Value::getRowIndex( const unsigned& irow, const unsigned& jind ) const {
497 : plumed_dbg_assert( (1+ncols)*irow+1+jind<matrix_bookeeping.size() && jind<matrix_bookeeping[(1+ncols)*irow] );
498 54232896 : return matrix_bookeeping[(1+ncols)*irow+1+jind];
499 : }
500 :
501 : inline
502 : unsigned Value::getNumberOfColumns() const {
503 17498698 : return ncols;
504 : }
505 :
506 : inline
507 : bool Value::isSymmetric() const {
508 1456 : return symmetric;
509 : }
510 :
511 : inline
512 : unsigned Value::getNumberOfGridDerivatives() const {
513 1911185 : return ngrid_der;
514 : }
515 :
516 : inline
517 : double Value::getGridDerivative(const unsigned& n, const unsigned& j ) const {
518 : plumed_dbg_assert( hasDeriv && n*(1+ngrid_der) + 1 + j < data.size() );
519 8329578 : return data[n*(1+ngrid_der) + 1 + j];
520 : }
521 :
522 : inline
523 1089247 : void Value::addGridDerivatives( const unsigned& n, const unsigned& j, const double& val ) {
524 : plumed_dbg_assert( hasDeriv && n*(1+ngrid_der) + 1 + j < data.size() );
525 1089247 : data[n*(1+ngrid_der) + 1 + j] += val;
526 1089247 : }
527 :
528 : inline
529 : void Value::setGridDerivatives( const unsigned& n, const unsigned& j, const double& val ) {
530 : plumed_dbg_assert( hasDeriv && n*(1+ngrid_der) + 1 + j < data.size() );
531 48400 : data[n*(1+ngrid_der) + 1 + j] = val;
532 : }
533 :
534 : }
535 : #endif
536 :
|
