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 tempory 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 : /// A flag telling us we have a force acting on this quantity
71 : bool hasForce=false;
72 : /// The way this value is used in the code
73 : /// normal = regular value that is determined during calculate
74 : /// constant = constnt value that is determined during startup and that doesn't change during simulation
75 : /// average = value that is averaged/collected over multiple steps of trajectory
76 : /// calcFromAverage = value that is calculated from an average value
77 : enum {normal,constant,average,calcFromAverage} valtype=normal;
78 : /// This is used by ActionWithValue to set the valtype
79 : void setValType( const std::string& vtype );
80 : /// This is used by ActionWithValue to determine if we need to calculate on update
81 : bool calculateOnUpdate() const ;
82 : /// The derivatives of the quantity stored in value
83 : std::map<AtomNumber,Vector> gradients;
84 : /// The name of this quantiy
85 : std::string name;
86 : /// 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)
87 : std::vector<std::size_t> shape;
88 : /// Does this quanity have derivatives
89 : bool hasDeriv=true;
90 : /// Variables for storing data
91 : unsigned bufstart=0;
92 : unsigned ngrid_der=0;
93 : std::size_t ncols=0;
94 : /// If we are storing a matrix is it symmetric?
95 : bool symmetric=false;
96 : /// This is a bookeeping array that holds the non-zero elements of the "sparse" matrix
97 : std::vector<unsigned> matrix_bookeeping;
98 : /// Is this quantity periodic
99 : enum {unset,periodic,notperiodic} periodicity=unset;
100 : /// Various quantities that describe the domain of this value
101 : std::string str_min;
102 : std::string str_max;
103 : double min=0.0;
104 : double max=0.0;
105 : double max_minus_min=0.0;
106 : double inv_max_minus_min=0.0;
107 : /// Is the derivative of this quantity zero when the value is zero
108 : bool derivativeIsZeroWhenValueIsZero=false;
109 : /// Complete the setup of the periodicity
110 : void setupPeriodicity();
111 : // bring value within PBCs
112 : void applyPeriodicity( const unsigned& ival );
113 : public:
114 : /// A constructor that can be used to make Vectors of values
115 : Value();
116 : /// A constructor that can be used to make Vectors of named values
117 : explicit Value(const std::string& valname);
118 : /// A constructor that is used throughout the code to setup the value poiters
119 : Value(ActionWithValue* av, const std::string& valname, const bool withderiv,const std::vector<std::size_t>&ss=std::vector<std::size_t>());
120 : /// Set the shape of the Value
121 : void setShape( const std::vector<std::size_t>&ss );
122 : /// Set the value of the function
123 : void set(double);
124 : /// Set the value of the stored data
125 : void set(const std::size_t& n, const double& v );
126 : /// Add something to the value of the function
127 : void add(double);
128 : /// Add something to the ith element of the data array
129 : void add(const std::size_t& n, const double& v );
130 : /// Get the location of this element of in the store
131 : std::size_t getIndexInStore( const std::size_t& ival ) const ;
132 : /// Get the value of the function
133 : double get( const std::size_t ival=0, const bool trueind=true ) const;
134 : /// A variant of get() for checking that at least one of the values on the row is !=0
135 : bool checkValueIsActiveForMMul(std::size_t task) const;
136 : /// A variant of get() for assigning data to an external view (assuems trueind=false), returns the number of arguments assigned
137 : std::size_t assignValues(View<double> target);
138 : /// Find out if the value has been set
139 : bool valueHasBeenSet() const;
140 : /// Check if the value is periodic
141 : bool isPeriodic() const;
142 : /// Set the function not periodic
143 : void setNotPeriodic();
144 : /// Set the domain of the function
145 : void setDomain(const std::string&, const std::string&);
146 : /// Get the domain of the quantity
147 : void getDomain(std::string&,std::string&) const;
148 : /// Get the domain of the quantity
149 : void getDomain(double&,double&) const;
150 : /// Get the name of the quantity
151 : const std::string& getName() const;
152 : /// Check whether or not this particular quantity has derivatives
153 : bool hasDerivatives()const;
154 : /// Get the number of derivatives that this particular value has
155 : unsigned getNumberOfDerivatives() const;
156 : /// Set the number of derivatives
157 : void resizeDerivatives(int n);
158 : /// Set all the derivatives to zero
159 : void clearDerivatives( const bool force=false );
160 : /// Add some derivative to the ith component of the derivatives array
161 : void addDerivative(unsigned i,double d);
162 : /// Set the value of the ith component of the derivatives array
163 : void setDerivative(unsigned i, double d);
164 : /// Get the derivative with respect to component n
165 : double getDerivative(const unsigned n) const;
166 : /// Clear the input force on the variable
167 : void clearInputForce();
168 : /// Special method for clearing forces on variables used by DataPassingObject
169 : void clearInputForce( const std::vector<AtomNumber>& index );
170 : /// Set hasForce equal to true
171 : void addForce();
172 : /// Add some force on this value
173 : void addForce(double f);
174 : /// Add some force on the ival th component of this value
175 : void addForce( const std::size_t& ival, double f, const bool trueind=true );
176 : ///Add forces from a vector, imples trueInd=false and retunrs the number of forces assigned
177 : std::size_t addForces(View<const double> f);
178 : /// Get the value of the force on this colvar
179 : double getForce( const std::size_t& ival=0 ) const ;
180 : /// Apply the forces to the derivatives using the chain rule (if there are no forces this routine returns false)
181 : bool applyForce( std::vector<double>& forces ) const ;
182 : /// Calculate the difference between the instantaneous value of the function and some other point: other_point-inst_val
183 : double difference(double)const;
184 : /// Calculate the difference between two values of this function: d2 -d1
185 : double difference(double d1,double d2)const;
186 : /// This returns the pointer to the action where this value is calculated
187 : ActionWithValue* getPntrToAction();
188 : /// Bring back one value into the correct pbc if needed, else give back the value
189 : double bringBackInPbc(double d1)const;
190 : /// Get the difference between max and minimum of domain
191 : double getMaxMinusMin()const;
192 : /// This sets up the gradients
193 : void setGradients( ActionAtomistic* aa, unsigned& start );
194 : /// This passes gradients from one action to another
195 : void passGradients( const double& der, std::map<AtomNumber,Vector>& g ) const ;
196 : static double projection(const Value&,const Value&);
197 : /// Get the rank of the object that is contained in this value
198 : unsigned getRank() const ;
199 : /// Get the shape of the object that is contained in this value
200 : const std::vector<std::size_t>& getShape() const ;
201 : /// Reshape the storage for sparse matrices
202 : void reshapeMatrixStore( const unsigned& n );
203 : /// Copy the matrix bookeeping stuff
204 : void copyBookeepingArrayFromArgument( Value* myarg );
205 : /// Set the symmetric flag equal true for this matrix
206 : void setSymmetric( const bool& sym );
207 : /// Get the total number of scalars that are stored here
208 : std::size_t getNumberOfValues() const ;
209 : /// Get the number of values that are actually stored here once sparse matrices are taken into account
210 : std::size_t getNumberOfStoredValues() const ;
211 : /// Get the number of threads to use when assigning this value
212 : unsigned getGoodNumThreads( const unsigned& j, const unsigned& k ) const ;
213 : /// These are used for passing around the data in this value when we are doing replica exchange
214 : void writeBinary(std::ostream&o) const ;
215 : void readBinary(std::istream&i);
216 : /// These are used for making constant values
217 : bool isConstant() const ;
218 : void setConstant();
219 : void reshapeConstantValue( const std::vector<std::size_t>& sh );
220 : /// Check if forces have been added on this value
221 : bool forcesWereAdded() const ;
222 : /// Set a bool that tells us if the derivative is zero when the value is zero true
223 : void setDerivativeIsZeroWhenValueIsZero();
224 : /// Return a bool that tells us if the derivative is zero when the value is zero
225 : bool isDerivativeZeroWhenValueIsZero() const ;
226 : /// Convert the input index to its corresponding indices
227 : void convertIndexToindices(const std::size_t& index, std::vector<unsigned>& indices ) const ;
228 : /// Print out all the values in this Value
229 : void print( OFile& ofile ) const ;
230 : /// Print out all the forces in this Value
231 : void printForce( OFile& ofile ) const ;
232 : /// Are we to ignore the stored value
233 : bool ignoreStoredValue(const std::string& n) const ;
234 : /// Set a matrix element to be non zero
235 : void setMatrixBookeepingElement( const unsigned& i, const unsigned& n );
236 : ///
237 : unsigned getRowLength( const std::size_t& irow ) const ;
238 : ///
239 : unsigned getRowIndex( std::size_t irow, std::size_t jind ) const ;
240 : ///
241 : void setRowIndices( const std::size_t& irow, const std::vector<std::size_t>& ind );
242 : ///
243 : std::size_t getNumberOfColumns() const ;
244 : ///
245 : bool isSymmetric() const ;
246 : /// Retrieve the non-zero edges in a matrix
247 : void retrieveEdgeList( unsigned& nedge, std::vector<std::pair<unsigned,unsigned> >& active, std::vector<double>& elems );
248 : /// Get the number of derivatives that the grid has
249 : unsigned getNumberOfGridDerivatives() const ;
250 : /// get the derivative of a grid at a point n with resepct to argument j
251 : double getGridDerivative(const unsigned& n, const unsigned& j ) const ;
252 : /// Add the derivatives of the grid to the corner
253 : void addGridDerivatives( const unsigned& n, const unsigned& j, const double& val );
254 : ///
255 : void setGridDerivatives( const unsigned& n, const unsigned& j, const double& val );
256 : /// Add another value to the end of the data vector held by this value. This is used in COLLECT
257 : void push_back( const double& val );
258 : /// Get the type of value that is stored here
259 : std::string getValueType() const ;
260 : /// Check if all the elements in the value have the same value
261 : bool allElementsEqual() const ;
262 : };
263 :
264 : inline
265 603494694 : void Value::applyPeriodicity(const unsigned& ival) {
266 603494694 : if(periodicity==periodic) {
267 3701550 : data[ival]=min+difference(min,data[ival]);
268 3701550 : if(data[ival]<min) {
269 1423738 : data[ival]+=max_minus_min;
270 : }
271 : }
272 603494694 : }
273 :
274 : inline
275 : void Value::set(double v) {
276 5355836 : value_set=true;
277 5577999 : data[0]=v;
278 5537845 : applyPeriodicity(0);
279 69304 : }
280 :
281 : inline
282 2148 : void Value::add(double v) {
283 2148 : value_set=true;
284 2148 : data[0]+=v;
285 2148 : applyPeriodicity(0);
286 2148 : }
287 :
288 : inline
289 10 : void Value::add(const std::size_t& n, const double& v ) {
290 10 : value_set=true;
291 10 : data[n]+=v;
292 10 : applyPeriodicity(n);
293 10 : }
294 :
295 : inline
296 : bool Value::valueHasBeenSet() const {
297 54 : return value_set;
298 : }
299 :
300 : inline
301 : const std::string& Value::getName()const {
302 5296613 : return name;
303 : }
304 :
305 : inline
306 67299 : unsigned Value::getNumberOfDerivatives() const {
307 67299 : plumed_massert(hasDeriv,"the derivatives array for this value has zero size");
308 67299 : if( shape.size()>0 ) {
309 2456 : return shape.size();
310 : }
311 64843 : return data.size() - 1;
312 : }
313 :
314 : inline
315 : double Value::getDerivative(const unsigned n) const {
316 : plumed_dbg_massert(n<getNumberOfDerivatives(),"you are asking for a derivative that is out of bounds");
317 19467571 : return data[1+n];
318 : }
319 :
320 : inline
321 : bool Value::hasDerivatives() const {
322 14823863 : return hasDeriv;
323 : }
324 :
325 : inline
326 4452615 : void Value::resizeDerivatives(int n) {
327 4452615 : if( shape.size()>0 ) {
328 : return;
329 : }
330 3555999 : if(hasDeriv) {
331 1868316 : data.resize(1+n);
332 : }
333 : }
334 :
335 : inline
336 : void Value::addDerivative(unsigned i,double d) {
337 : plumed_dbg_massert(i<getNumberOfDerivatives(),"derivative is out of bounds");
338 20475446 : data[1+i]+=d;
339 315 : }
340 :
341 : inline
342 : void Value::setDerivative(unsigned i, double d) {
343 : plumed_dbg_massert(i<getNumberOfDerivatives(),"derivative is out of bounds");
344 1512682 : data[1+i]=d;
345 0 : }
346 :
347 : inline
348 : void Value::clearInputForce() {
349 3294356 : if( !hasForce ) {
350 : return;
351 : }
352 697645 : hasForce=false;
353 : std::fill(inputForce.begin(),inputForce.end(),0);
354 : }
355 :
356 : inline
357 117270 : void Value::clearInputForce( const std::vector<AtomNumber>& index ) {
358 117270 : if( !hasForce ) {
359 : return;
360 : }
361 62802 : hasForce=false;
362 2019756 : for(const auto & p : index) {
363 1956954 : inputForce[p.index()]=0;
364 : }
365 : }
366 :
367 : inline
368 2812598 : void Value::clearDerivatives( const bool force ) {
369 2812598 : if( !force && (valtype==constant || valtype==average) ) {
370 : return;
371 : }
372 :
373 2792490 : value_set=false;
374 2792490 : if( shape.size()>0 ) {
375 : std::fill(data.begin(), data.end(), 0);
376 2505707 : } else if( data.size()>1 ) {
377 : std::fill(data.begin()+1, data.end(), 0);
378 : }
379 : }
380 :
381 : inline
382 : void Value::addForce() {
383 37 : hasForce=true;
384 : }
385 :
386 : inline
387 : void Value::addForce(double f) {
388 185785 : hasForce=true;
389 185785 : inputForce[0]+=f;
390 2623 : }
391 :
392 : inline
393 : bool Value::forcesWereAdded() const {
394 1255382 : return hasForce;
395 : }
396 :
397 : inline
398 : double Value::getForce( const std::size_t& ival ) const {
399 : plumed_dbg_assert( ival<inputForce.size() );
400 125611113 : return inputForce[ival];
401 : }
402 :
403 : /// d2-d1
404 : inline
405 12552985 : double Value::difference(double d1,double d2)const {
406 12552985 : if(periodicity==notperiodic) {
407 4583138 : return d2-d1;
408 7969847 : } else if(periodicity==periodic) {
409 7969847 : double s=(d2-d1)*inv_max_minus_min;
410 : // remember: pbc brings the difference in a range of -0.5:0.5
411 7969847 : s=Tools::pbc(s);
412 7969847 : return s*max_minus_min;
413 : } else {
414 0 : plumed_merror("periodicity should be set to compute differences");
415 : }
416 : }
417 :
418 : inline
419 3173 : double Value::bringBackInPbc(double d1)const {
420 3173 : return min+max_minus_min/2.+difference(min+max_minus_min/2., d1);
421 : }
422 :
423 : inline
424 3989526 : double Value::difference(double d)const {
425 3989526 : return difference(get(),d);
426 : }
427 :
428 : inline
429 : double Value::getMaxMinusMin()const {
430 : plumed_dbg_assert( periodicity==periodic );
431 0 : return max_minus_min;
432 : }
433 :
434 : inline
435 630639273 : unsigned Value::getRank() const {
436 630639273 : if( valtype==constant && shape.size()==1 && shape[0]==1 ) {
437 : return 0;
438 : }
439 630638690 : return shape.size();
440 : }
441 :
442 : inline
443 : const std::vector<std::size_t>& Value::getShape() const {
444 518048 : return shape;
445 : }
446 :
447 : inline
448 49828276 : std::size_t Value::getNumberOfValues() const {
449 : std::size_t size=1;
450 64400850 : for(unsigned i=0; i<shape.size(); ++i) {
451 14572574 : size *= shape[i];
452 : }
453 49828276 : return size;
454 : }
455 :
456 : inline
457 11176658 : std::size_t Value::getNumberOfStoredValues() const {
458 11176658 : if( getRank()==2 && !hasDeriv ) {
459 1032546 : return shape[0]*ncols;
460 : }
461 10144112 : return getNumberOfValues();
462 : }
463 :
464 : inline
465 : bool Value::isConstant() const {
466 353986 : return valtype==constant;
467 : }
468 :
469 : inline
470 : void Value::setDerivativeIsZeroWhenValueIsZero() {
471 3167 : derivativeIsZeroWhenValueIsZero=true;
472 305 : }
473 :
474 : inline
475 : bool Value::isDerivativeZeroWhenValueIsZero() const {
476 12430341 : return derivativeIsZeroWhenValueIsZero;
477 : }
478 :
479 : inline
480 : void Value::setMatrixBookeepingElement( const unsigned& i, const unsigned& n ) {
481 : plumed_dbg_assert( i<matrix_bookeeping.size() );
482 65856363 : matrix_bookeeping[i]=n;
483 : }
484 :
485 : inline
486 : unsigned Value::getRowLength( const std::size_t& irow ) const {
487 16845576659 : if( matrix_bookeeping.size()==0 ) {
488 : return 0;
489 : }
490 : plumed_dbg_assert( (1+ncols)*irow<matrix_bookeeping.size() );
491 16845576659 : return matrix_bookeeping[(1+ncols)*irow];
492 : }
493 :
494 : inline
495 : unsigned Value::getRowIndex(const std::size_t irow, const std::size_t jind ) const {
496 : plumed_dbg_massert( (1+ncols)*irow+1+jind<matrix_bookeeping.size() && jind<matrix_bookeeping[(1+ncols)*irow], "failing in value " + name );
497 16841377931 : return matrix_bookeeping[(1+ncols)*irow+1+jind];
498 : }
499 :
500 : inline
501 3548 : void Value::setRowIndices( const std::size_t& irow, const std::vector<std::size_t>& ind ) {
502 : plumed_dbg_massert( (1+ncols)*irow+1+ind.size()<=matrix_bookeeping.size(), "problem in " + name );
503 3548 : std::size_t istart = (1+ncols)*irow;
504 3548 : matrix_bookeeping[istart] = ind.size();
505 3548 : ++istart;
506 18632 : for(unsigned i=0; i<ind.size(); ++i) {
507 15084 : matrix_bookeeping[istart] = ind[i];
508 15084 : ++istart;
509 : }
510 3548 : }
511 :
512 : inline
513 : std::size_t Value::getNumberOfColumns() const {
514 6242169 : return ncols;
515 : }
516 :
517 : inline
518 : bool Value::isSymmetric() const {
519 1944 : return symmetric;
520 : }
521 :
522 : inline
523 : unsigned Value::getNumberOfGridDerivatives() const {
524 6724 : return ngrid_der;
525 : }
526 :
527 : inline
528 : double Value::getGridDerivative(const unsigned& n, const unsigned& j ) const {
529 : plumed_dbg_assert( hasDeriv && n*(1+ngrid_der) + 1 + j < data.size() );
530 7955705 : return data[n*(1+ngrid_der) + 1 + j];
531 : }
532 :
533 : inline
534 554747 : void Value::addGridDerivatives( const unsigned& n, const unsigned& j, const double& val ) {
535 : plumed_dbg_assert( hasDeriv && n*(1+ngrid_der) + 1 + j < data.size() );
536 554747 : data[n*(1+ngrid_der) + 1 + j] += val;
537 554747 : }
538 :
539 : inline
540 : void Value::setGridDerivatives( const unsigned& n, const unsigned& j, const double& val ) {
541 : plumed_dbg_assert( hasDeriv && n*(1+ngrid_der) + 1 + j < data.size() );
542 3442324 : data[n*(1+ngrid_der) + 1 + j] = val;
543 : }
544 :
545 : }
546 : #endif
547 :
|
