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 ActionWithValue;
36 :
37 : /// \ingroup TOOLBOX
38 : /// A class for holding the value of a function together with its derivatives.
39 : /// Typically, an object of type PLMD::ActionWithValue will contain one
40 : /// object of type PLUMD::Value that will be named after the label. If the
41 : /// PLMD::ActionWithValue is part of a class that calculates multiple components
42 : /// then the class will contain multiple that will be called label.component-name
43 : /// This class is used to pass information between different PLMD::Action
44 : /// objects. However, if you find a use for a tempory PLMD::Value in some method
45 : /// you are implementing please feel free to use it.
46 : class Value {
47 : friend class ActionWithValue;
48 : /// This copies the contents of a value into a second value (just the derivatives and value)
49 : friend void copy( const Value& val1, Value& val2 );
50 : /// This copies the contents of a value into a second value (but second value is a pointer)
51 : friend void copy( const Value& val, Value* val2 );
52 : /// This adds some derivatives onto the value
53 : friend void add( const Value& val1, Value* valout );
54 : /// This calculates val1*val2 and sorts out the derivatives
55 : friend void product( const Value& val1, const Value& val2, Value& valout );
56 : /// This calculates va1/val2 and sorts out the derivatives
57 : friend void quotient( const Value& val1, const Value& val2, Value* valout );
58 : private:
59 : /// The action in which this quantity is calculated
60 : ActionWithValue* action;
61 : /// Had the value been set
62 : bool value_set;
63 : /// The value of the quantity
64 : double value;
65 : /// The force acting on this quantity
66 : double inputForce;
67 : /// A flag telling us we have a force acting on this quantity
68 : bool hasForce;
69 : /// The derivatives of the quantity stored in value
70 : std::vector<double> derivatives;
71 : std::map<AtomNumber,Vector> gradients;
72 : /// The name of this quantiy
73 : std::string name;
74 : /// Does this quanity have derivatives
75 : bool hasDeriv;
76 : /// Is this quantity periodic
77 : enum {unset,periodic,notperiodic} periodicity;
78 : /// Various quantities that describe the domain of this value
79 : std::string str_min, str_max;
80 : double min,max;
81 : double max_minus_min;
82 : double inv_max_minus_min;
83 : /// Complete the setup of the periodicity
84 : void setupPeriodicity();
85 : // bring value within PBCs
86 : void applyPeriodicity();
87 : public:
88 : /// A constructor that can be used to make Vectors of values
89 : Value();
90 : /// A constructor that can be used to make Vectors of named values
91 : explicit Value(const std::string& name);
92 : /// A constructor that is used throughout the code to setup the value poiters
93 : Value(ActionWithValue* av, const std::string& name, const bool withderiv);
94 : /// Set the value of the function
95 : void set(double);
96 : /// Add something to the value of the function
97 : void add(double);
98 : /// Get the value of the function
99 : double get() const;
100 : /// Find out if the value has been set
101 : bool valueHasBeenSet() const;
102 : /// Check if the value is periodic
103 : bool isPeriodic() const;
104 : /// Set the function not periodic
105 : void setNotPeriodic();
106 : /// Set the domain of the function
107 : void setDomain(const std::string&, const std::string&);
108 : /// Get the domain of the quantity
109 : void getDomain(std::string&,std::string&) const;
110 : /// Get the domain of the quantity
111 : void getDomain(double&,double&) const;
112 : /// Get the name of the quantity
113 : const std::string& getName() const;
114 : /// Check whether or not this particular quantity has derivatives
115 : bool hasDerivatives()const;
116 : /// Get the number of derivatives that this particular value has
117 : unsigned getNumberOfDerivatives() const;
118 : /// Set the number of derivatives
119 : void resizeDerivatives(int n);
120 : /// Set all the derivatives to zero
121 : void clearDerivatives();
122 : /// Add some derivative to the ith component of the derivatives array
123 : void addDerivative(unsigned i,double d);
124 : /// Set the value of the ith component of the derivatives array
125 : void setDerivative(unsigned i, double d);
126 : /// Apply the chain rule to the derivatives
127 : void chainRule(double df);
128 : /// Get the derivative with respect to component n
129 : double getDerivative(const unsigned n) const;
130 : /// Clear the input force on the variable
131 : void clearInputForce();
132 : /// Add some force on this value
133 : void addForce(double f);
134 : /// Get the value of the force on this colvar
135 : double getForce() const ;
136 : /// Apply the forces to the derivatives using the chain rule (if there are no forces this routine returns false)
137 : bool applyForce( std::vector<double>& forces ) const ;
138 : /// Calculate the difference between the instantaneous value of the function and some other point: other_point-inst_val
139 : double difference(double)const;
140 : /// Calculate the difference between two values of this function: d2 -d1
141 : double difference(double d1,double d2)const;
142 : /// This returns the pointer to the action where this value is calculated
143 : ActionWithValue* getPntrToAction();
144 : /// Bring back one value into the correct pbc if needed, else give back the value
145 : double bringBackInPbc(double d1)const;
146 : /// Get the difference between max and minimum of domain
147 : double getMaxMinusMin()const;
148 : /// This sets up the gradients
149 : void setGradients();
150 : static double projection(const Value&,const Value&);
151 : };
152 :
153 : void copy( const Value& val1, Value& val2 );
154 : void copy( const Value& val1, Value* val2 );
155 : void add( const Value& val1, Value* valout );
156 :
157 : inline
158 44362470 : void Value::applyPeriodicity() {
159 44362470 : if(periodicity==periodic) {
160 30706139 : value=min+difference(min,value);
161 30706139 : if(value<min)value+=max_minus_min;
162 : }
163 44362470 : }
164 :
165 : inline
166 : void product( const Value& val1, const Value& val2, Value& valout ) {
167 : plumed_assert( val1.derivatives.size()==val2.derivatives.size() );
168 : if( valout.derivatives.size()!=val1.derivatives.size() ) valout.resizeDerivatives( val1.derivatives.size() );
169 : valout.value_set=false;
170 : valout.clearDerivatives();
171 : double u=val1.value;
172 : double v=val2.value;
173 : for(unsigned i=0; i<val1.derivatives.size(); ++i) {
174 : valout.addDerivative(i, u*val2.derivatives[i] + v*val1.derivatives[i] );
175 : }
176 : valout.set( u*v );
177 : }
178 :
179 : inline
180 : void quotient( const Value& val1, const Value& val2, Value* valout ) {
181 : plumed_assert( val1.derivatives.size()==val2.derivatives.size() );
182 : if( valout->derivatives.size()!=val1.derivatives.size() ) valout->resizeDerivatives( val1.derivatives.size() );
183 : valout->value_set=false;
184 : valout->clearDerivatives();
185 : double u=val1.get();
186 : double v=val2.get();
187 : for(unsigned i=0; i<val1.getNumberOfDerivatives(); ++i) {
188 : valout->addDerivative(i, v*val1.getDerivative(i) - u*val2.getDerivative(i) );
189 : }
190 : valout->chainRule( 1/(v*v) ); valout->set( u / v );
191 : }
192 :
193 : inline
194 : void Value::set(double v) {
195 44351760 : value_set=true;
196 44351760 : value=v;
197 44162997 : applyPeriodicity();
198 42537046 : }
199 :
200 : inline
201 : void Value::add(double v) {
202 : value_set=true;
203 : value+=v;
204 : applyPeriodicity();
205 : }
206 :
207 : inline
208 : double Value::get()const {
209 3224447 : return value;
210 : }
211 :
212 : inline
213 : bool Value::valueHasBeenSet() const {
214 0 : return value_set;
215 : }
216 :
217 : inline
218 : const std::string& Value::getName()const {
219 8263808 : return name;
220 : }
221 :
222 : inline
223 329575 : unsigned Value::getNumberOfDerivatives() const {
224 329575 : plumed_massert(hasDeriv,"the derivatives array for this value has zero size");
225 329575 : return derivatives.size();
226 : }
227 :
228 : inline
229 : double Value::getDerivative(const unsigned n) const {
230 : plumed_dbg_massert(n<derivatives.size(),"you are asking for a derivative that is out of bounds");
231 19967312 : return derivatives[n];
232 : }
233 :
234 : inline
235 : bool Value::hasDerivatives() const {
236 13272 : return hasDeriv;
237 : }
238 :
239 : inline
240 : void Value::resizeDerivatives(int n) {
241 4102750470 : if(hasDeriv) derivatives.resize(n);
242 : }
243 :
244 : inline
245 : void Value::addDerivative(unsigned i,double d) {
246 : plumed_dbg_massert(i<derivatives.size(),"derivative is out of bounds");
247 22368600 : derivatives[i]+=d;
248 10 : }
249 :
250 : inline
251 : void Value::setDerivative(unsigned i, double d) {
252 : plumed_dbg_massert(i<derivatives.size(),"derivative is out of bounds");
253 63116 : derivatives[i]=d;
254 : }
255 :
256 : inline
257 434 : void Value::chainRule(double df) {
258 18368 : for(unsigned i=0; i<derivatives.size(); ++i) derivatives[i]*=df;
259 434 : }
260 :
261 : inline
262 : void Value::clearInputForce() {
263 2081053 : hasForce=false;
264 2081053 : inputForce=0.0;
265 : }
266 :
267 : inline
268 : void Value::clearDerivatives() {
269 2523 : value_set=false;
270 : std::fill(derivatives.begin(), derivatives.end(), 0);
271 : }
272 :
273 : inline
274 : void Value::addForce(double f) {
275 : plumed_dbg_massert(hasDerivatives(),"forces can only be added to values with derivatives");
276 188450 : hasForce=true;
277 188450 : inputForce+=f;
278 2623 : }
279 :
280 : inline
281 : double Value::getForce() const {
282 4823 : return inputForce;
283 : }
284 : /// d2-d1
285 : inline
286 155957424 : double Value::difference(double d1,double d2)const {
287 155957424 : if(periodicity==notperiodic) {
288 90588656 : return d2-d1;
289 65368768 : } else if(periodicity==periodic) {
290 65368768 : double s=(d2-d1)*inv_max_minus_min;
291 : // remember: pbc brings the difference in a range of -0.5:0.5
292 65368768 : s=Tools::pbc(s);
293 65368768 : return s*max_minus_min;
294 0 : } else plumed_merror("periodicity should be set to compute differences");
295 : }
296 :
297 : inline
298 3173 : double Value::bringBackInPbc(double d1)const {
299 3173 : return min+max_minus_min/2.+difference(min+max_minus_min/2., d1);
300 : }
301 :
302 : inline
303 : double Value::difference(double d)const {
304 119526971 : return difference(get(),d);
305 : }
306 :
307 : inline
308 : double Value::getMaxMinusMin()const {
309 : plumed_dbg_assert( periodicity==periodic );
310 518 : return max_minus_min;
311 : }
312 :
313 : }
314 :
315 : #endif
316 :
|
