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 6839559 : 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 56822909 : void Value::applyPeriodicity(const unsigned& ival) {
255 56822909 : if(periodicity==periodic) {
256 2386213 : data[ival]=min+difference(min,data[ival]);
257 2386213 : if(data[ival]<min)data[ival]+=max_minus_min;
258 : }
259 56822909 : }
260 :
261 : inline
262 : void Value::set(double v) {
263 5312386 : value_set=true;
264 5535434 : data[0]=v;
265 5491272 : applyPeriodicity(0);
266 64142 : }
267 :
268 : inline
269 2144 : void Value::add(double v) {
270 2144 : value_set=true;
271 2144 : data[0]+=v;
272 2144 : applyPeriodicity(0);
273 2144 : }
274 :
275 : inline
276 46452643 : void Value::add(const std::size_t& n, const double& v ) {
277 46452643 : value_set=true; data[n]+=v; applyPeriodicity(n);
278 46452643 : }
279 :
280 : inline
281 : bool Value::valueHasBeenSet() const {
282 231208647 : return value_set;
283 : }
284 :
285 : inline
286 : const std::string& Value::getName()const {
287 5378457 : return name;
288 : }
289 :
290 : inline
291 16294106 : unsigned Value::getNumberOfDerivatives() const {
292 16294106 : plumed_massert(hasDeriv,"the derivatives array for this value has zero size");
293 16294106 : if( shape.size()>0 ) return shape.size();
294 16294106 : return data.size() - 1;
295 : }
296 :
297 : inline
298 : double Value::getDerivative(const unsigned n) const {
299 : plumed_dbg_massert(n<getNumberOfDerivatives(),"you are asking for a derivative that is out of bounds");
300 19322011 : return data[1+n];
301 : }
302 :
303 : inline
304 : bool Value::hasDerivatives() const {
305 225345145 : return hasDeriv;
306 : }
307 :
308 : inline
309 4448078 : void Value::resizeDerivatives(int n) {
310 4448078 : if( shape.size()>0 ) return;
311 3551835 : if(hasDeriv) data.resize(1+n);
312 : }
313 :
314 : inline
315 : void Value::addDerivative(unsigned i,double d) {
316 : plumed_dbg_massert(i<getNumberOfDerivatives(),"derivative is out of bounds");
317 20064107 : data[1+i]+=d;
318 : }
319 :
320 : inline
321 : void Value::setDerivative(unsigned i, double d) {
322 : plumed_dbg_massert(i<getNumberOfDerivatives(),"derivative is out of bounds");
323 17649718 : data[1+i]=d;
324 : }
325 :
326 : inline
327 : void Value::clearInputForce() {
328 3256152 : if( !hasForce ) return;
329 548999 : hasForce=false; std::fill(inputForce.begin(),inputForce.end(),0);
330 : }
331 :
332 : inline
333 116735 : void Value::clearInputForce( const std::vector<AtomNumber>& index ) {
334 116735 : if( !hasForce ) return;
335 2003241 : hasForce=false; for(const auto & p : index) inputForce[p.index()]=0;
336 : }
337 :
338 : inline
339 2810205 : void Value::clearDerivatives( const bool force ) {
340 2810205 : if( !force && (valtype==constant || valtype==average) ) return;
341 :
342 2790102 : value_set=false;
343 2790102 : if( data.size()>1 ) std::fill(data.begin()+1, data.end(), 0);
344 : }
345 :
346 : inline
347 : void Value::addForce(double f) {
348 179001 : hasForce=true;
349 179001 : inputForce[0]+=f;
350 2623 : }
351 :
352 : inline
353 : bool Value::forcesWereAdded() const {
354 19678787 : return hasForce;
355 : }
356 :
357 : inline
358 : double Value::getForce( const std::size_t& ival ) const {
359 : plumed_dbg_assert( ival<inputForce.size() );
360 24554279 : return inputForce[ival];
361 : }
362 :
363 : /// d2-d1
364 : inline
365 135679951 : double Value::difference(double d1,double d2)const {
366 135679951 : if(periodicity==notperiodic) {
367 33852260 : return d2-d1;
368 101827691 : } else if(periodicity==periodic) {
369 101827691 : double s=(d2-d1)*inv_max_minus_min;
370 : // remember: pbc brings the difference in a range of -0.5:0.5
371 101827691 : s=Tools::pbc(s);
372 101827691 : return s*max_minus_min;
373 0 : } else plumed_merror("periodicity should be set to compute differences");
374 : }
375 :
376 : inline
377 3173 : double Value::bringBackInPbc(double d1)const {
378 3173 : return min+max_minus_min/2.+difference(min+max_minus_min/2., d1);
379 : }
380 :
381 : inline
382 3989526 : double Value::difference(double d)const {
383 3989526 : return difference(get(),d);
384 : }
385 :
386 : inline
387 : double Value::getMaxMinusMin()const {
388 : plumed_dbg_assert( periodicity==periodic );
389 0 : return max_minus_min;
390 : }
391 :
392 : inline
393 : unsigned Value::getRank() const {
394 708195681 : return shape.size();
395 : }
396 :
397 : inline
398 : const std::vector<unsigned>& Value::getShape() const {
399 471251 : return shape;
400 : }
401 :
402 : inline
403 45659110 : unsigned Value::getNumberOfValues() const {
404 55648774 : unsigned size=1; for(unsigned i=0; i<shape.size(); ++i) size *= shape[i];
405 45659110 : return size;
406 : }
407 :
408 : inline
409 2713975 : unsigned Value::getNumberOfStoredValues() const {
410 2713975 : if( getRank()==2 && !hasDeriv ) return shape[0]*ncols;
411 449604 : return getNumberOfValues();
412 : }
413 :
414 : inline
415 : bool Value::isConstant() const {
416 397480 : return valtype==constant;
417 : }
418 :
419 : inline
420 : void Value::setDerivativeIsZeroWhenValueIsZero() {
421 1199 : derivativeIsZeroWhenValueIsZero=true;
422 15 : }
423 :
424 : inline
425 : bool Value::isDerivativeZeroWhenValueIsZero() const {
426 1798 : return derivativeIsZeroWhenValueIsZero;
427 : }
428 :
429 : inline
430 : unsigned Value::getPositionInStream() const {
431 1052590383 : return streampos;
432 : }
433 :
434 : inline
435 : unsigned Value::getPositionInMatrixStash() const {
436 28657011 : return matpos;
437 : }
438 :
439 : inline
440 : void Value::setMatrixBookeepingStart( const unsigned& b ) {
441 7794 : book_start = b;
442 : }
443 :
444 : inline
445 : unsigned Value::getMatrixBookeepingStart() const {
446 14169812 : return book_start;
447 : }
448 :
449 : inline
450 : void Value::setMatrixBookeepingElement( const unsigned& i, const unsigned& n ) {
451 : plumed_dbg_assert( i<matrix_bookeeping.size() );
452 5060232 : matrix_bookeeping[i]=n;
453 : }
454 :
455 : inline
456 : bool Value::valueIsStored() const {
457 117898873 : return storedata;
458 : }
459 :
460 : inline
461 : unsigned Value::getRowLength( const unsigned& irow ) const {
462 : plumed_dbg_assert( (1+ncols)*irow<matrix_bookeeping.size() );
463 34407153 : return matrix_bookeeping[(1+ncols)*irow];
464 : }
465 :
466 : inline
467 : unsigned Value::getRowIndex( const unsigned& irow, const unsigned& jind ) const {
468 : plumed_dbg_assert( (1+ncols)*irow+1+jind<matrix_bookeeping.size() && jind<matrix_bookeeping[(1+ncols)*irow] );
469 54232896 : return matrix_bookeeping[(1+ncols)*irow+1+jind];
470 : }
471 :
472 : inline
473 : unsigned Value::getNumberOfColumns() const {
474 17498698 : return ncols;
475 : }
476 :
477 : inline
478 : bool Value::isSymmetric() const {
479 1456 : return symmetric;
480 : }
481 :
482 : inline
483 : unsigned Value::getNumberOfGridDerivatives() const {
484 1911185 : return ngrid_der;
485 : }
486 :
487 : inline
488 : double Value::getGridDerivative(const unsigned& n, const unsigned& j ) const {
489 : plumed_dbg_assert( hasDeriv && n*(1+ngrid_der) + 1 + j < data.size() );
490 8329578 : return data[n*(1+ngrid_der) + 1 + j];
491 : }
492 :
493 : inline
494 1089247 : void Value::addGridDerivatives( const unsigned& n, const unsigned& j, const double& val ) {
495 : plumed_dbg_assert( hasDeriv && n*(1+ngrid_der) + 1 + j < data.size() );
496 1089247 : data[n*(1+ngrid_der) + 1 + j] += val;
497 1089247 : }
498 :
499 : inline
500 : void Value::setGridDerivatives( const unsigned& n, const unsigned& j, const double& val ) {
501 : plumed_dbg_assert( hasDeriv && n*(1+ngrid_der) + 1 + j < data.size() );
502 48400 : data[n*(1+ngrid_der) + 1 + j] = val;
503 : }
504 :
505 : }
506 : #endif
507 :
|
