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 : #include "tools/View.h"
33 :
34 : namespace PLMD {
35 :
36 : class OFile;
37 : class Communicator;
38 : class ActionWithValue;
39 : class ActionAtomistic;
40 :
41 : /// \ingroup TOOLBOX
42 : /// A class for holding the value of a function together with its derivatives.
43 : /// Typically, an object of type PLMD::ActionWithValue will contain one
44 : /// object of type PLUMD::Value that will be named after the label. If the
45 : /// PLMD::ActionWithValue is part of a class that calculates multiple components
46 : /// then the class will contain multiple that will be called label.component-name
47 : /// This class is used to pass information between different PLMD::Action
48 : /// objects. However, if you find a use for a temporary PLMD::Value in some method
49 : /// you are implementing please feel free to use it.
50 : class Value {
51 : friend struct ArgumentsBookkeeping;
52 : friend class ActionWithValue;
53 : friend class ActionWithVector;
54 : friend class ActionWithMatrix;
55 : friend class ActionAtomistic;
56 : friend class ActionWithArguments;
57 : friend class ActionWithVirtualAtom;
58 : friend class DomainDecomposition;
59 : template<typename T>
60 : friend class DataPassingObjectTyped;
61 : private:
62 : /// The action in which this quantity is calculated
63 : ActionWithValue* action=nullptr;
64 : /// Had the value been set
65 : bool value_set=false;
66 : /// The value of the quantity
67 : std::vector<double> data;
68 : /// The force acting on this quantity
69 : std::vector<double> inputForce;
70 : /// The value whose matrix bookkeeping this one was copied from (nullptr: its own).
71 : /// Lets the task manager pack the bookkeeping of arguments that share it once.
72 : const Value* bookeeping_root=nullptr;
73 : /// A flag telling us we have a force acting on this quantity
74 : bool hasForce=false;
75 : /// The way this value is used in the code
76 : /// normal = regular value that is determined during calculate
77 : /// constant = constnt value that is determined during startup and that doesn't change during simulation
78 : /// average = value that is averaged/collected over multiple steps of trajectory
79 : /// calcFromAverage = value that is calculated from an average value
80 : enum {normal,constant,average,calcFromAverage} valtype=normal;
81 : /// This is used by ActionWithValue to set the valtype
82 : void setValType( const std::string& vtype );
83 : /// This is used by ActionWithValue to determine if we need to calculate on update
84 : bool calculateOnUpdate() const ;
85 : /// The derivatives of the quantity stored in value
86 : std::map<AtomNumber,Vector> gradients;
87 : /// The name of this quantiy
88 : std::string name;
89 : /// 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)
90 : std::vector<std::size_t> shape;
91 : /// Does this quanity have derivatives
92 : bool hasDeriv=true;
93 : /// Variables for storing data
94 : unsigned bufstart=0;
95 : unsigned ngrid_der=0;
96 : std::size_t ncols=0;
97 : /// If we are storing a matrix is it symmetric?
98 : bool symmetric=false;
99 : /// This is a bookeeping array that holds the non-zero elements of the "sparse" matrix
100 : std::vector<unsigned> matrix_bookeeping;
101 : /// Is this quantity periodic
102 : enum {unset,periodic,notperiodic} periodicity=unset;
103 : /// Various quantities that describe the domain of this value
104 : std::string str_min;
105 : std::string str_max;
106 : double min=0.0;
107 : double max=0.0;
108 : double max_minus_min=0.0;
109 : double inv_max_minus_min=0.0;
110 : /// Is the derivative of this quantity zero when the value is zero
111 : bool derivativeIsZeroWhenValueIsZero=false;
112 : /// Complete the setup of the periodicity
113 : void setupPeriodicity();
114 : // bring value within PBCs
115 : void applyPeriodicity( const unsigned& ival );
116 : /// Add a force on an element of a sparse matrix (the uncommon case of addForce)
117 : void addForceToSparseMatrix( const std::size_t& iforce, double f );
118 : /// Error raised by addForce, out of line so that addForce can be inlined
119 : [[noreturn]] void addForceOutOfRange( const std::size_t& iforce ) const ;
120 : public:
121 : /// A constructor that can be used to make Vectors of values
122 : Value();
123 : /// A constructor that can be used to make Vectors of named values
124 : explicit Value(const std::string& valname);
125 : /// A constructor that is used throughout the code to setup the value poiters
126 : Value(ActionWithValue* av, const std::string& valname, const bool withderiv,const std::vector<std::size_t>&ss=std::vector<std::size_t>());
127 : /// Set the shape of the Value
128 : void setShape( const std::vector<std::size_t>&ss );
129 : /// Set the value of the function
130 : void set(double);
131 : /// Set the value of the stored data
132 : void set(const std::size_t& n, const double& v );
133 : /// Add something to the value of the function
134 : void add(double);
135 : /// Add something to the ith element of the data array
136 : void add(const std::size_t& n, const double& v );
137 : /// Get the location of this element of in the store
138 : std::size_t getIndexInStore( const std::size_t& ival ) const ;
139 : /// Get the value of the function
140 : double get( const std::size_t ival=0, const bool trueind=true ) const;
141 : /// A variant of get() for checking that at least one of the values on the row is !=0
142 : bool checkValueIsActiveForMMul(std::size_t task) const;
143 : /// A variant of get() for assigning data to an external view (assuems trueind=false), returns the number of arguments assigned
144 : std::size_t assignValues(View<double> target);
145 : /// Copy the stored values converted to single precision
146 : std::size_t assignValues(View<float> target);
147 : /// Find out if the value has been set
148 : bool valueHasBeenSet() const;
149 : /// Check if the value is periodic
150 : bool isPeriodic() const;
151 : /// Set the function not periodic
152 : void setNotPeriodic();
153 : /// Set the domain of the function
154 : void setDomain(const std::string&, const std::string&);
155 : /// Get the domain of the quantity
156 : void getDomain(std::string&,std::string&) const;
157 : /// Get the domain of the quantity
158 : void getDomain(double&,double&) const;
159 : /// Get the name of the quantity
160 : const std::string& getName() const;
161 : /// Check whether or not this particular quantity has derivatives
162 : bool hasDerivatives()const;
163 : /// Get the number of derivatives that this particular value has
164 : unsigned getNumberOfDerivatives() const;
165 : /// Set the number of derivatives
166 : void resizeDerivatives(int n);
167 : /// Set all the derivatives to zero
168 : void clearDerivatives( const bool force=false );
169 : /// Add some derivative to the ith component of the derivatives array
170 : void addDerivative(unsigned i,double d);
171 : /// Set the value of the ith component of the derivatives array
172 : void setDerivative(unsigned i, double d);
173 : /// Get the derivative with respect to component n
174 : double getDerivative(const unsigned n) const;
175 : /// Clear the input force on the variable
176 : void clearInputForce();
177 : /// Special method for clearing forces on variables used by DataPassingObject
178 : void clearInputForce( const std::vector<AtomNumber>& index );
179 : /// Set hasForce equal to true
180 : void addForce();
181 : /// Add some force on this value
182 : void addForce(double f);
183 : /// Add some force on the ival th component of this value
184 : void addForce( const std::size_t& ival, double f, const bool trueind=true );
185 : ///Add forces from a vector, imples trueInd=false and retunrs the number of forces assigned
186 : std::size_t addForces(View<const double> f);
187 : /// Get the value of the force on this colvar
188 : double getForce( const std::size_t& ival=0 ) const ;
189 : /// Apply the forces to the derivatives using the chain rule (if there are no forces this routine returns false)
190 : bool applyForce( std::vector<double>& forces ) const ;
191 : /// Calculate the difference between the instantaneous value of the function and some other point: other_point-inst_val
192 : double difference(double)const;
193 : /// Calculate the difference between two values of this function: d2 -d1
194 : double difference(double d1,double d2)const;
195 : /// This returns the pointer to the action where this value is calculated
196 : ActionWithValue* getPntrToAction();
197 : /// Bring back one value into the correct pbc if needed, else give back the value
198 : double bringBackInPbc(double d1)const;
199 : /// Get the difference between max and minimum of domain
200 : double getMaxMinusMin()const;
201 : /// This sets up the gradients
202 : void setGradients( ActionAtomistic* aa, unsigned& start );
203 : /// This passes gradients from one action to another
204 : void passGradients( const double& der, std::map<AtomNumber,Vector>& g ) const ;
205 : static double projection(const Value&,const Value&);
206 : /// Get the rank of the object that is contained in this value
207 : unsigned getRank() const ;
208 : /// Get the shape of the object that is contained in this value
209 : const std::vector<std::size_t>& getShape() const ;
210 : /// Reshape the storage for sparse matrices
211 : void reshapeMatrixStore( const unsigned& n );
212 : /// Copy the matrix bookeeping stuff
213 : void copyBookeepingArrayFromArgument( Value* myarg );
214 : /// Set the symmetric flag equal true for this matrix
215 : void setSymmetric( const bool& sym );
216 : /// Get the total number of scalars that are stored here
217 : std::size_t getNumberOfValues() const ;
218 : /// Get the number of values that are actually stored here once sparse matrices are taken into account
219 : std::size_t getNumberOfStoredValues() const ;
220 : /// Get the number of threads to use when assigning this value
221 : unsigned getGoodNumThreads( const unsigned& j, const unsigned& k ) const ;
222 : /// These are used for passing around the data in this value when we are doing replica exchange
223 : void writeBinary(std::ostream&o) const ;
224 : void readBinary(std::istream&i);
225 : /// These are used for making constant values
226 : bool isConstant() const ;
227 : void setConstant();
228 : void reshapeConstantValue( const std::vector<std::size_t>& sh );
229 : /// Check if forces have been added on this value
230 : bool forcesWereAdded() const ;
231 : /// Set a bool that tells us if the derivative is zero when the value is zero true
232 : void setDerivativeIsZeroWhenValueIsZero();
233 : /// Return a bool that tells us if the derivative is zero when the value is zero
234 : bool isDerivativeZeroWhenValueIsZero() const ;
235 : /// Convert the input index to its corresponding indices
236 : void convertIndexToindices(const std::size_t& index, std::vector<unsigned>& indices ) const ;
237 : /// Print out all the values in this Value
238 : void print( OFile& ofile ) const ;
239 : /// Print out all the forces in this Value
240 : void printForce( OFile& ofile ) const ;
241 : /// Are we to ignore the stored value
242 : bool ignoreStoredValue(const std::string& n) const ;
243 : /// Set a matrix element to be non zero
244 : void setMatrixBookeepingElement( const unsigned& i, const unsigned& n );
245 : /// The value that holds the original of the matrix bookkeeping of this value (see bookeeping_root)
246 : const Value* getBookeepingRoot() const ;
247 : ///
248 : unsigned getRowLength( const std::size_t& irow ) const ;
249 : ///
250 : unsigned getRowIndex( std::size_t irow, std::size_t jind ) const ;
251 : ///
252 : void setRowIndices( const std::size_t& irow, const std::vector<std::size_t>& ind );
253 : ///
254 : std::size_t getLengthOfLongestColumn() const ;
255 : ///
256 : std::size_t getNumberOfColumns() const ;
257 : ///
258 : bool isSymmetric() const ;
259 : /// Retrieve the non-zero edges in a matrix
260 : void retrieveEdgeList( unsigned& nedge, std::vector<std::pair<unsigned,unsigned> >& active, std::vector<double>& elems );
261 : /// Get the number of derivatives that the grid has
262 : unsigned getNumberOfGridDerivatives() const ;
263 : /// get the derivative of a grid at a point n with resepct to argument j
264 : double getGridDerivative(const unsigned& n, const unsigned& j ) const ;
265 : /// Add the derivatives of the grid to the corner
266 : void addGridDerivatives( const unsigned& n, const unsigned& j, const double& val );
267 : ///
268 : void setGridDerivatives( const unsigned& n, const unsigned& j, const double& val );
269 : /// Add another value to the end of the data vector held by this value. This is used in COLLECT
270 : void push_back( const double& val );
271 : /// Get the type of value that is stored here
272 : std::string getValueType() const ;
273 : /// Check if all the elements in the value have the same value
274 : bool allElementsEqual() const ;
275 : };
276 :
277 : inline
278 34009417 : void Value::applyPeriodicity(const unsigned& ival) {
279 34009417 : if(periodicity==periodic) {
280 3415604 : data[ival]=min+difference(min,data[ival]);
281 3415604 : if(data[ival]<min) {
282 1292954 : data[ival]+=max_minus_min;
283 : }
284 : }
285 34009417 : }
286 :
287 : inline
288 : void Value::set(double v) {
289 5595271 : value_set=true;
290 5663732 : data[0]=v;
291 5696522 : applyPeriodicity(0);
292 84148 : }
293 :
294 : inline
295 2377 : void Value::add(double v) {
296 2377 : value_set=true;
297 2377 : data[0]+=v;
298 2377 : applyPeriodicity(0);
299 2377 : }
300 :
301 : inline
302 10 : void Value::add(const std::size_t& n, const double& v ) {
303 10 : value_set=true;
304 10 : data[n]+=v;
305 10 : applyPeriodicity(n);
306 10 : }
307 :
308 : inline
309 : bool Value::valueHasBeenSet() const {
310 92 : return value_set;
311 : }
312 :
313 : inline
314 : const std::string& Value::getName()const {
315 5222687 : return name;
316 : }
317 :
318 : inline
319 69080 : unsigned Value::getNumberOfDerivatives() const {
320 69080 : plumed_massert(hasDeriv,"the derivatives array for this value has zero size");
321 69080 : if( shape.size()>0 ) {
322 2762 : return shape.size();
323 : }
324 66318 : return data.size() - 1;
325 : }
326 :
327 : inline
328 : double Value::getDerivative(const unsigned n) const {
329 : plumed_dbg_massert(n<getNumberOfDerivatives(),"you are asking for a derivative that is out of bounds");
330 20474195 : return data[1+n];
331 : }
332 :
333 : inline
334 : bool Value::hasDerivatives() const {
335 15730375 : return hasDeriv;
336 : }
337 :
338 : inline
339 4504546 : void Value::resizeDerivatives(int n) {
340 4504546 : if( shape.size()>0 ) {
341 : return;
342 : }
343 3581891 : if(hasDeriv) {
344 1888317 : data.resize(1+n);
345 : }
346 : }
347 :
348 : inline
349 : void Value::addDerivative(unsigned i,double d) {
350 : plumed_dbg_massert(i<getNumberOfDerivatives(),"derivative is out of bounds");
351 21578244 : data[1+i]+=d;
352 123 : }
353 :
354 : inline
355 : void Value::setDerivative(unsigned i, double d) {
356 : plumed_dbg_massert(i<getNumberOfDerivatives(),"derivative is out of bounds");
357 1545444 : data[1+i]=d;
358 0 : }
359 :
360 : inline
361 : void Value::clearInputForce() {
362 4489255 : if( !hasForce ) {
363 : return;
364 : }
365 1142383 : hasForce=false;
366 : std::fill(inputForce.begin(),inputForce.end(),0);
367 : }
368 :
369 : inline
370 182290 : void Value::clearInputForce( const std::vector<AtomNumber>& index ) {
371 182290 : if( !hasForce ) {
372 : return;
373 : }
374 101253 : hasForce=false;
375 2228361 : for(const auto & p : index) {
376 2127108 : inputForce[p.index()]=0;
377 : }
378 : }
379 :
380 : inline
381 3448791 : void Value::clearDerivatives( const bool force ) {
382 3448791 : if( !force && (valtype==constant || valtype==average) ) {
383 : return;
384 : }
385 :
386 3428165 : value_set=false;
387 3428165 : if( shape.size()>0 ) {
388 : std::fill(data.begin(), data.end(), 0);
389 3147382 : } else if( data.size()>1 ) {
390 : std::fill(data.begin()+1, data.end(), 0);
391 : }
392 : }
393 :
394 : inline
395 : void Value::addForce() {
396 37 : hasForce=true;
397 : }
398 :
399 : inline
400 : void Value::addForce(double f) {
401 414364 : hasForce=true;
402 414364 : inputForce[0]+=f;
403 3231 : }
404 :
405 : inline
406 11044601 : void Value::addForce(const std::size_t& iforce, double f, const bool trueind) {
407 11044601 : hasForce=true;
408 11044601 : if( shape.size()==2 && !hasDeriv && ncols<shape[1] && trueind ) {
409 1425 : addForceToSparseMatrix( iforce, f );
410 1425 : return;
411 : }
412 11043176 : if( iforce>=inputForce.size() ) {
413 0 : addForceOutOfRange( iforce );
414 : }
415 11043176 : inputForce[iforce]+=f;
416 : }
417 :
418 : inline
419 : bool Value::forcesWereAdded() const {
420 2181500 : return hasForce;
421 : }
422 :
423 : inline
424 : double Value::getForce( const std::size_t& ival ) const {
425 : plumed_dbg_assert( ival<inputForce.size() );
426 9025175 : return inputForce[ival];
427 : }
428 :
429 : /// d2-d1
430 : inline
431 10973605 : double Value::difference(double d1,double d2)const {
432 10973605 : if(periodicity==notperiodic) {
433 3403939 : return d2-d1;
434 7569666 : } else if(periodicity==periodic) {
435 7569666 : double s=(d2-d1)*inv_max_minus_min;
436 : // remember: pbc brings the difference in a range of -0.5:0.5
437 7569666 : s=Tools::pbc(s);
438 7569666 : return s*max_minus_min;
439 : } else {
440 0 : plumed_merror("periodicity should be set to compute differences");
441 : }
442 : }
443 :
444 : inline
445 19574 : double Value::bringBackInPbc(double d1)const {
446 19574 : return min+max_minus_min/2.+difference(min+max_minus_min/2., d1);
447 : }
448 :
449 : inline
450 3882016 : double Value::difference(double d)const {
451 3882016 : return difference(get(),d);
452 : }
453 :
454 : inline
455 : double Value::getMaxMinusMin()const {
456 : plumed_dbg_assert( periodicity==periodic );
457 238 : return max_minus_min;
458 : }
459 :
460 : inline
461 62498406 : unsigned Value::getRank() const {
462 62498406 : if( valtype==constant && shape.size()==1 && shape[0]==1 ) {
463 : return 0;
464 : }
465 62497827 : return shape.size();
466 : }
467 :
468 : inline
469 : const std::vector<std::size_t>& Value::getShape() const {
470 1158208 : return shape;
471 : }
472 :
473 : inline
474 52144693 : std::size_t Value::getNumberOfValues() const {
475 : std::size_t size=1;
476 70737377 : for(unsigned i=0; i<shape.size(); ++i) {
477 18592684 : size *= shape[i];
478 : }
479 52144693 : return size;
480 : }
481 :
482 : inline
483 11760202 : std::size_t Value::getNumberOfStoredValues() const {
484 11760202 : if( getRank()==2 && !hasDeriv ) {
485 1061853 : return shape[0]*ncols;
486 : }
487 10698349 : return getNumberOfValues();
488 : }
489 :
490 : inline
491 : bool Value::isConstant() const {
492 592438 : return valtype==constant;
493 : }
494 :
495 : inline
496 : void Value::setDerivativeIsZeroWhenValueIsZero() {
497 3724 : derivativeIsZeroWhenValueIsZero=true;
498 432 : }
499 :
500 : inline
501 : bool Value::isDerivativeZeroWhenValueIsZero() const {
502 11231932 : return derivativeIsZeroWhenValueIsZero;
503 : }
504 :
505 : inline
506 : const Value* Value::getBookeepingRoot() const {
507 16250 : return bookeeping_root ? bookeeping_root : this;
508 : }
509 :
510 : inline
511 : void Value::setMatrixBookeepingElement( const unsigned& i, const unsigned& n ) {
512 : plumed_dbg_assert( i<matrix_bookeeping.size() );
513 39465390 : bookeeping_root=nullptr;
514 39465390 : matrix_bookeeping[i]=n;
515 : }
516 :
517 : inline
518 : unsigned Value::getRowLength( const std::size_t& irow ) const {
519 1246568566 : if( matrix_bookeeping.size()==0 ) {
520 : return 0;
521 : }
522 : plumed_dbg_assert( (1+ncols)*irow<matrix_bookeeping.size() );
523 1246568566 : return matrix_bookeeping[(1+ncols)*irow];
524 : }
525 :
526 : inline
527 : unsigned Value::getRowIndex(const std::size_t irow, const std::size_t jind ) const {
528 : plumed_dbg_massert( (1+ncols)*irow+1+jind<matrix_bookeeping.size() && jind<matrix_bookeeping[(1+ncols)*irow], "failing in value " + name );
529 1258126231 : return matrix_bookeeping[(1+ncols)*irow+1+jind];
530 : }
531 :
532 : inline
533 3548 : void Value::setRowIndices( const std::size_t& irow, const std::vector<std::size_t>& ind ) {
534 : plumed_dbg_massert( (1+ncols)*irow+1+ind.size()<=matrix_bookeeping.size(), "problem in " + name );
535 3548 : std::size_t istart = (1+ncols)*irow;
536 3548 : matrix_bookeeping[istart] = ind.size();
537 3548 : ++istart;
538 18632 : for(unsigned i=0; i<ind.size(); ++i) {
539 15084 : matrix_bookeeping[istart] = ind[i];
540 15084 : ++istart;
541 : }
542 3548 : }
543 :
544 : inline
545 : std::size_t Value::getNumberOfColumns() const {
546 6381045 : return ncols;
547 : }
548 :
549 : inline
550 : bool Value::isSymmetric() const {
551 2302 : return symmetric;
552 : }
553 :
554 : inline
555 : unsigned Value::getNumberOfGridDerivatives() const {
556 7534 : return ngrid_der;
557 : }
558 :
559 : inline
560 : double Value::getGridDerivative(const unsigned& n, const unsigned& j ) const {
561 : plumed_dbg_assert( hasDeriv && n*(1+ngrid_der) + 1 + j < data.size() );
562 7981321 : return data[n*(1+ngrid_der) + 1 + j];
563 : }
564 :
565 : inline
566 563312 : void Value::addGridDerivatives( const unsigned& n, const unsigned& j, const double& val ) {
567 : plumed_dbg_assert( hasDeriv && n*(1+ngrid_der) + 1 + j < data.size() );
568 563312 : data[n*(1+ngrid_der) + 1 + j] += val;
569 563312 : }
570 :
571 : inline
572 : void Value::setGridDerivatives( const unsigned& n, const unsigned& j, const double& val ) {
573 : plumed_dbg_assert( hasDeriv && n*(1+ngrid_der) + 1 + j < data.size() );
574 3470841 : data[n*(1+ngrid_der) + 1 + j] = val;
575 : }
576 :
577 : }
578 : #endif
579 :
|
