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 7226691 : 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 : /// Derivative writes to constant scalars are ignored.
152 : void addDerivative(unsigned i,double d);
153 : /// Set the value of the ith component of the derivatives array
154 : /// Derivative writes to constant scalars are ignored.
155 : void setDerivative(unsigned i, double d);
156 : /// Get the derivative with respect to component n
157 : double getDerivative(const unsigned n) const;
158 : /// Clear the input force on the variable
159 : void clearInputForce();
160 : /// Special method for clearing forces on variables used by DataPassingObject
161 : void clearInputForce( const std::vector<AtomNumber>& index );
162 : /// Add some force on this value
163 : void addForce(double f);
164 : /// Add some force on the ival th component of this value
165 : void addForce( const std::size_t& ival, double f, const bool trueind=true );
166 : /// Get the value of the force on this colvar
167 : double getForce( const std::size_t& ival=0 ) const ;
168 : /// Apply the forces to the derivatives using the chain rule (if there are no forces this routine returns false)
169 : bool applyForce( std::vector<double>& forces ) const ;
170 : /// Calculate the difference between the instantaneous value of the function and some other point: other_point-inst_val
171 : double difference(double)const;
172 : /// Calculate the difference between two values of this function: d2 -d1
173 : double difference(double d1,double d2)const;
174 : /// This returns the pointer to the action where this value is calculated
175 : ActionWithValue* getPntrToAction();
176 : /// Bring back one value into the correct pbc if needed, else give back the value
177 : double bringBackInPbc(double d1)const;
178 : /// Get the difference between max and minimum of domain
179 : double getMaxMinusMin()const;
180 : /// This sets up the gradients
181 : void setGradients( ActionAtomistic* aa, unsigned& start );
182 : /// This passes gradients from one action to another
183 : void passGradients( const double& der, std::map<AtomNumber,Vector>& g ) const ;
184 : static double projection(const Value&,const Value&);
185 : /// Get the rank of the object that is contained in this value
186 : unsigned getRank() const ;
187 : /// Get the shape of the object that is contained in this value
188 : const std::vector<unsigned>& getShape() const ;
189 : /// This turns on storing of vectors/matrices
190 : void buildDataStore( const bool forprint=false );
191 : /// Reshape the storage for sparse matrices
192 : void reshapeMatrixStore( const unsigned& n );
193 : /// Set the symmetric flag equal true for this matrix
194 : void setSymmetric( const bool& sym );
195 : /// Get the total number of scalars that are stored here
196 : unsigned getNumberOfValues() const ;
197 : /// Get the number of values that are actually stored here once sparse matrices are taken into account
198 : unsigned getNumberOfStoredValues() const ;
199 : /// Get the number of threads to use when assigning this value
200 : unsigned getGoodNumThreads( const unsigned& j, const unsigned& k ) const ;
201 : /// These are used for passing around the data in this value when we are doing replica exchange
202 : void writeBinary(std::ostream&o) const ;
203 : void readBinary(std::istream&i);
204 : /// These are used for making constant values
205 : bool isConstant() const ;
206 : void setConstant();
207 : /// Check if forces have been added on this value
208 : bool forcesWereAdded() const ;
209 : /// Set a bool that tells us if the derivative is zero when the value is zero true
210 : void setDerivativeIsZeroWhenValueIsZero();
211 : /// Return a bool that tells us if the derivative is zero when the value is zero
212 : bool isDerivativeZeroWhenValueIsZero() const ;
213 : ///
214 : unsigned getPositionInStream() const ;
215 : /// This stuff handles where to look for the start of the row that contains the row of the matrix
216 : void setPositionInMatrixStash( const unsigned& p );
217 : unsigned getPositionInMatrixStash() const ;
218 : /// This stuff handles where to keep the bookeeping stuff for storing the sparse matrix
219 : void setMatrixBookeepingStart( const unsigned& b );
220 : unsigned getMatrixBookeepingStart() const ;
221 : /// Convert the input index to its corresponding indices
222 : void convertIndexToindices(const std::size_t& index, std::vector<unsigned>& indices ) const ;
223 : /// Print out all the values in this Value
224 : void print( OFile& ofile ) const ;
225 : /// Are we to ignore the stored value
226 : bool ignoreStoredValue(const std::string& n) const ;
227 : /// Set a matrix element to be non zero
228 : void setMatrixBookeepingElement( const unsigned& i, const unsigned& n );
229 : ///
230 : unsigned getRowLength( const unsigned& irow ) const ;
231 : ///
232 : unsigned getRowIndex( const unsigned& irow, const unsigned& jind ) const ;
233 : /// Are we storing this value
234 : bool valueIsStored() const ;
235 : ///
236 : unsigned getNumberOfColumns() const ;
237 : ///
238 : bool isSymmetric() const ;
239 : /// Retrieve the non-zero edges in a matrix
240 : void retrieveEdgeList( unsigned& nedge, std::vector<std::pair<unsigned,unsigned> >& active, std::vector<double>& elems );
241 : /// Get the number of derivatives that the grid has
242 : unsigned getNumberOfGridDerivatives() const ;
243 : /// get the derivative of a grid at a point n with resepct to argument j
244 : double getGridDerivative(const unsigned& n, const unsigned& j ) const ;
245 : /// Add the derivatives of the grid to the corner
246 : void addGridDerivatives( const unsigned& n, const unsigned& j, const double& val );
247 : ///
248 : void setGridDerivatives( const unsigned& n, const unsigned& j, const double& val );
249 : /// Add another value to the end of the data vector held by this value. This is used in COLLECT
250 : void push_back( const double& val );
251 : /// Get the type of value that is stored here
252 : std::string getValueType() const ;
253 : };
254 :
255 : inline
256 60743030 : void Value::applyPeriodicity(const unsigned& ival) {
257 60743030 : if(periodicity==periodic) {
258 2386914 : data[ival]=min+difference(min,data[ival]);
259 2386914 : if(data[ival]<min) {
260 790797 : data[ival]+=max_minus_min;
261 : }
262 : }
263 60743030 : }
264 :
265 : inline
266 : void Value::set(double v) {
267 5581241 : value_set=true;
268 5803384 : data[0]=v;
269 5763273 : applyPeriodicity(0);
270 68188 : }
271 :
272 : inline
273 2152 : void Value::add(double v) {
274 2152 : value_set=true;
275 2152 : data[0]+=v;
276 2152 : applyPeriodicity(0);
277 2152 : }
278 :
279 : inline
280 47255742 : void Value::add(const std::size_t& n, const double& v ) {
281 47255742 : value_set=true;
282 47255742 : data[n]+=v;
283 47255742 : applyPeriodicity(n);
284 47255742 : }
285 :
286 : inline
287 : bool Value::valueHasBeenSet() const {
288 714925640 : return value_set;
289 : }
290 :
291 : inline
292 : const std::string& Value::getName()const {
293 5424339 : return name;
294 : }
295 :
296 : inline
297 10256256 : unsigned Value::getNumberOfDerivatives() const {
298 10256256 : plumed_massert(hasDeriv,"the derivatives array for this value has zero size");
299 10256256 : if( shape.size()>0 ) {
300 0 : return shape.size();
301 : }
302 10256256 : return data.size() - 1;
303 : }
304 :
305 : inline
306 : double Value::getDerivative(const unsigned n) const {
307 : plumed_dbg_massert(n<getNumberOfDerivatives(),"you are asking for a derivative that is out of bounds");
308 19337925 : return data[1+n];
309 : }
310 :
311 : inline
312 : bool Value::hasDerivatives() const {
313 446869417 : return hasDeriv;
314 : }
315 :
316 : inline
317 4455975 : void Value::resizeDerivatives(int n) {
318 4455975 : if( shape.size()>0 ) {
319 : return;
320 : }
321 3554056 : if(hasDeriv) {
322 1868949 : data.resize(1+n);
323 : }
324 : }
325 :
326 : inline
327 : void Value::addDerivative(unsigned i,double d) {
328 : plumed_dbg_assert( shape.empty() );
329 : // Constant scalars store only their value, without derivative storage.
330 71444088 : if(valtype==constant) {
331 : return;
332 : }
333 : plumed_dbg_massert(i<getNumberOfDerivatives(),"derivative is out of bounds");
334 71444056 : data[1+i]+=d;
335 : }
336 :
337 : inline
338 : void Value::setDerivative(unsigned i, double d) {
339 : plumed_dbg_assert( shape.empty() );
340 : // Constant scalars store only their value, without derivative storage.
341 13943876 : if(valtype==constant) {
342 : return;
343 : }
344 : plumed_dbg_massert(i<getNumberOfDerivatives(),"derivative is out of bounds");
345 13940778 : data[1+i]=d;
346 : }
347 :
348 : inline
349 : void Value::clearInputForce() {
350 4258249 : if( !hasForce ) {
351 : return;
352 : }
353 951693 : hasForce=false;
354 : std::fill(inputForce.begin(),inputForce.end(),0);
355 : }
356 :
357 : inline
358 118825 : void Value::clearInputForce( const std::vector<AtomNumber>& index ) {
359 118825 : if( !hasForce ) {
360 : return;
361 : }
362 63762 : hasForce=false;
363 2019957 : for(const auto & p : index) {
364 1956195 : inputForce[p.index()]=0;
365 : }
366 : }
367 :
368 : inline
369 3379131 : void Value::clearDerivatives( const bool force ) {
370 3379131 : if( !force && (valtype==constant || valtype==average) ) {
371 : return;
372 : }
373 :
374 3358968 : value_set=false;
375 3358968 : if( data.size()>1 ) {
376 : std::fill(data.begin()+1, data.end(), 0);
377 : }
378 : }
379 :
380 : inline
381 : void Value::addForce(double f) {
382 312488 : hasForce=true;
383 312488 : inputForce[0]+=f;
384 3231 : }
385 :
386 : inline
387 : bool Value::forcesWereAdded() const {
388 24790098 : return hasForce;
389 : }
390 :
391 : inline
392 : double Value::getForce( const std::size_t& ival ) const {
393 : plumed_dbg_assert( ival<inputForce.size() );
394 24934718 : return inputForce[ival];
395 : }
396 :
397 : /// d2-d1
398 : inline
399 141153121 : double Value::difference(double d1,double d2)const {
400 141153121 : if(periodicity==notperiodic) {
401 39331655 : return d2-d1;
402 101821466 : } else if(periodicity==periodic) {
403 101821466 : double s=(d2-d1)*inv_max_minus_min;
404 : // remember: pbc brings the difference in a range of -0.5:0.5
405 101821466 : s=Tools::pbc(s);
406 101821466 : return s*max_minus_min;
407 : } else {
408 0 : plumed_merror("periodicity should be set to compute differences");
409 : }
410 : }
411 :
412 : inline
413 3319 : double Value::bringBackInPbc(double d1)const {
414 3319 : return min+max_minus_min/2.+difference(min+max_minus_min/2., d1);
415 : }
416 :
417 : inline
418 3989526 : double Value::difference(double d)const {
419 3989526 : return difference(get(),d);
420 : }
421 :
422 : inline
423 : double Value::getMaxMinusMin()const {
424 : plumed_dbg_assert( periodicity==periodic );
425 44 : return max_minus_min;
426 : }
427 :
428 : inline
429 : unsigned Value::getRank() const {
430 1802680863 : return shape.size();
431 : }
432 :
433 : inline
434 : const std::vector<unsigned>& Value::getShape() const {
435 1003994 : return shape;
436 : }
437 :
438 : inline
439 46149144 : unsigned Value::getNumberOfValues() const {
440 : unsigned size=1;
441 56852471 : for(unsigned i=0; i<shape.size(); ++i) {
442 10703327 : size *= shape[i];
443 : }
444 46149144 : return size;
445 : }
446 :
447 : inline
448 2716252 : unsigned Value::getNumberOfStoredValues() const {
449 2716252 : if( getRank()==2 && !hasDeriv ) {
450 2264373 : return shape[0]*ncols;
451 : }
452 451879 : return getNumberOfValues();
453 : }
454 :
455 : inline
456 : bool Value::isConstant() const {
457 543093 : return valtype==constant;
458 : }
459 :
460 : inline
461 : void Value::setDerivativeIsZeroWhenValueIsZero() {
462 1215 : derivativeIsZeroWhenValueIsZero=true;
463 15 : }
464 :
465 : inline
466 : bool Value::isDerivativeZeroWhenValueIsZero() const {
467 2001 : return derivativeIsZeroWhenValueIsZero;
468 : }
469 :
470 : inline
471 : unsigned Value::getPositionInStream() const {
472 2535309574 : return streampos;
473 : }
474 :
475 : inline
476 : unsigned Value::getPositionInMatrixStash() const {
477 45376201 : return matpos;
478 : }
479 :
480 : inline
481 : void Value::setMatrixBookeepingStart( const unsigned& b ) {
482 7798 : book_start = b;
483 : }
484 :
485 : inline
486 : unsigned Value::getMatrixBookeepingStart() const {
487 18488832 : return book_start;
488 : }
489 :
490 : inline
491 : void Value::setMatrixBookeepingElement( const unsigned& i, const unsigned& n ) {
492 : plumed_dbg_assert( i<matrix_bookeeping.size() );
493 5060232 : matrix_bookeeping[i]=n;
494 : }
495 :
496 : inline
497 : bool Value::valueIsStored() const {
498 219642692 : return storedata;
499 : }
500 :
501 : inline
502 : unsigned Value::getRowLength( const unsigned& irow ) const {
503 : plumed_dbg_assert( (1+ncols)*irow<matrix_bookeeping.size() );
504 37846135 : return matrix_bookeeping[(1+ncols)*irow];
505 : }
506 :
507 : inline
508 : unsigned Value::getRowIndex( const unsigned& irow, const unsigned& jind ) const {
509 : plumed_dbg_assert( (1+ncols)*irow+1+jind<matrix_bookeeping.size() && jind<matrix_bookeeping[(1+ncols)*irow] );
510 145642238 : return matrix_bookeeping[(1+ncols)*irow+1+jind];
511 : }
512 :
513 : inline
514 : unsigned Value::getNumberOfColumns() const {
515 200649953 : return ncols;
516 : }
517 :
518 : inline
519 : bool Value::isSymmetric() const {
520 1458 : return symmetric;
521 : }
522 :
523 : inline
524 : unsigned Value::getNumberOfGridDerivatives() const {
525 1911721 : return ngrid_der;
526 : }
527 :
528 : inline
529 : double Value::getGridDerivative(const unsigned& n, const unsigned& j ) const {
530 : plumed_dbg_assert( hasDeriv && n*(1+ngrid_der) + 1 + j < data.size() );
531 8330694 : return data[n*(1+ngrid_der) + 1 + j];
532 : }
533 :
534 : inline
535 1089991 : void Value::addGridDerivatives( const unsigned& n, const unsigned& j, const double& val ) {
536 : plumed_dbg_assert( hasDeriv && n*(1+ngrid_der) + 1 + j < data.size() );
537 1089991 : data[n*(1+ngrid_der) + 1 + j] += val;
538 1089991 : }
539 :
540 : inline
541 : void Value::setGridDerivatives( const unsigned& n, const unsigned& j, const double& val ) {
542 : plumed_dbg_assert( hasDeriv && n*(1+ngrid_der) + 1 + j < data.size() );
543 48400 : data[n*(1+ngrid_der) + 1 + j] = val;
544 : }
545 :
546 : }
547 : #endif
548 :
|
