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 44392773 : void Value::applyPeriodicity() {
159 44392773 : if(periodicity==periodic) {
160 30707608 : value=min+difference(min,value);
161 30707608 : if(value<min) {
162 14872392 : value+=max_minus_min;
163 : }
164 : }
165 44392773 : }
166 :
167 : inline
168 : void product( const Value& val1, const Value& val2, Value& valout ) {
169 : plumed_assert( val1.derivatives.size()==val2.derivatives.size() );
170 : if( valout.derivatives.size()!=val1.derivatives.size() ) {
171 : valout.resizeDerivatives( val1.derivatives.size() );
172 : }
173 : valout.value_set=false;
174 : valout.clearDerivatives();
175 : double u=val1.value;
176 : double v=val2.value;
177 : for(unsigned i=0; i<val1.derivatives.size(); ++i) {
178 : valout.addDerivative(i, u*val2.derivatives[i] + v*val1.derivatives[i] );
179 : }
180 : valout.set( u*v );
181 : }
182 :
183 : inline
184 : void quotient( const Value& val1, const Value& val2, Value* valout ) {
185 : plumed_assert( val1.derivatives.size()==val2.derivatives.size() );
186 : if( valout->derivatives.size()!=val1.derivatives.size() ) {
187 : valout->resizeDerivatives( val1.derivatives.size() );
188 : }
189 : valout->value_set=false;
190 : valout->clearDerivatives();
191 : double u=val1.get();
192 : double v=val2.get();
193 : for(unsigned i=0; i<val1.getNumberOfDerivatives(); ++i) {
194 : valout->addDerivative(i, v*val1.getDerivative(i) - u*val2.getDerivative(i) );
195 : }
196 : valout->chainRule( 1/(v*v) );
197 : valout->set( u / v );
198 : }
199 :
200 : inline
201 : void Value::set(double v) {
202 44381919 : value_set=true;
203 44381919 : value=v;
204 44192936 : applyPeriodicity();
205 42557058 : }
206 :
207 : inline
208 : void Value::add(double v) {
209 : value_set=true;
210 : value+=v;
211 : applyPeriodicity();
212 : }
213 :
214 : inline
215 : double Value::get()const {
216 3248012 : return value;
217 : }
218 :
219 : inline
220 : bool Value::valueHasBeenSet() const {
221 0 : return value_set;
222 : }
223 :
224 : inline
225 : const std::string& Value::getName()const {
226 8676977 : return name;
227 : }
228 :
229 : inline
230 378648 : unsigned Value::getNumberOfDerivatives() const {
231 378648 : plumed_massert(hasDeriv,"the derivatives array for this value has zero size");
232 378648 : return derivatives.size();
233 : }
234 :
235 : inline
236 : double Value::getDerivative(const unsigned n) const {
237 : plumed_dbg_massert(n<derivatives.size(),"you are asking for a derivative that is out of bounds");
238 20376386 : return derivatives[n];
239 : }
240 :
241 : inline
242 : bool Value::hasDerivatives() const {
243 13276 : return hasDeriv;
244 : }
245 :
246 : inline
247 : void Value::resizeDerivatives(int n) {
248 4102753188 : if(hasDeriv) {
249 2053106339 : derivatives.resize(n);
250 : }
251 : }
252 :
253 : inline
254 : void Value::addDerivative(unsigned i,double d) {
255 : plumed_dbg_massert(i<derivatives.size(),"derivative is out of bounds");
256 23202491 : derivatives[i]+=d;
257 10 : }
258 :
259 : inline
260 : void Value::setDerivative(unsigned i, double d) {
261 : plumed_dbg_massert(i<derivatives.size(),"derivative is out of bounds");
262 63116 : derivatives[i]=d;
263 : }
264 :
265 : inline
266 434 : void Value::chainRule(double df) {
267 18368 : for(unsigned i=0; i<derivatives.size(); ++i) {
268 17934 : derivatives[i]*=df;
269 : }
270 434 : }
271 :
272 : inline
273 : void Value::clearInputForce() {
274 2104624 : hasForce=false;
275 2104624 : inputForce=0.0;
276 : }
277 :
278 : inline
279 : void Value::clearDerivatives() {
280 2523 : value_set=false;
281 : std::fill(derivatives.begin(), derivatives.end(), 0);
282 : }
283 :
284 : inline
285 : void Value::addForce(double f) {
286 : plumed_dbg_massert(hasDerivatives(),"forces can only be added to values with derivatives");
287 194130 : hasForce=true;
288 194130 : inputForce+=f;
289 2623 : }
290 :
291 : inline
292 : double Value::getForce() const {
293 4985 : return inputForce;
294 : }
295 : /// d2-d1
296 : inline
297 155961574 : double Value::difference(double d1,double d2)const {
298 155961574 : if(periodicity==notperiodic) {
299 90591338 : return d2-d1;
300 65370236 : } else if(periodicity==periodic) {
301 65370236 : double s=(d2-d1)*inv_max_minus_min;
302 : // remember: pbc brings the difference in a range of -0.5:0.5
303 65370236 : s=Tools::pbc(s);
304 65370236 : return s*max_minus_min;
305 : } else {
306 0 : plumed_merror("periodicity should be set to compute differences");
307 : }
308 : }
309 :
310 : inline
311 3173 : double Value::bringBackInPbc(double d1)const {
312 3173 : return min+max_minus_min/2.+difference(min+max_minus_min/2., d1);
313 : }
314 :
315 : inline
316 : double Value::difference(double d)const {
317 119526971 : return difference(get(),d);
318 : }
319 :
320 : inline
321 : double Value::getMaxMinusMin()const {
322 : plumed_dbg_assert( periodicity==periodic );
323 518 : return max_minus_min;
324 : }
325 :
326 : }
327 :
328 : #endif
329 :
|
