Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2019 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 7293113 : 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 is used throughout the code to setup the value poiters
91 : Value(ActionWithValue* av, const std::string& name, const bool withderiv);
92 : /// Set the value of the function
93 : void set(double);
94 : /// Add something to the value of the function
95 : void add(double);
96 : /// Get the value of the function
97 : double get() const;
98 : /// Find out if the value has been set
99 : bool valueHasBeenSet() const;
100 : /// Check if the value is periodic
101 : bool isPeriodic() const;
102 : /// Set the function not periodic
103 : void setNotPeriodic();
104 : /// Set the domain of the function
105 : void setDomain(const std::string&, const std::string&);
106 : /// Get the domain of the quantity
107 : void getDomain(std::string&,std::string&) const;
108 : /// Get the domain of the quantity
109 : void getDomain(double&,double&) const;
110 : /// Get the name of the quantity
111 : const std::string& getName() const;
112 : /// Check whether or not this particular quantity has derivatives
113 : bool hasDerivatives()const;
114 : /// Get the number of derivatives that this particular value has
115 : unsigned getNumberOfDerivatives() const;
116 : /// Set the number of derivatives
117 : void resizeDerivatives(int n);
118 : /// Set all the derivatives to zero
119 : void clearDerivatives();
120 : /// Add some derivative to the ith component of the derivatives array
121 : void addDerivative(unsigned i,double d);
122 : /// Set the value of the ith component of the derivatives array
123 : void setDerivative(unsigned i, double d);
124 : /// Apply the chain rule to the derivatives
125 : void chainRule(double df);
126 : /// Get the derivative with respect to component n
127 : double getDerivative(const unsigned n) const;
128 : /// Clear the input force on the variable
129 : void clearInputForce();
130 : /// Add some force on this value
131 : void addForce(double f);
132 : /// Get the value of the force on this colvar
133 : double getForce() const ;
134 : /// Apply the forces to the derivatives using the chain rule (if there are no forces this routine returns false)
135 : bool applyForce( std::vector<double>& forces ) const ;
136 : /// Calculate the difference between the instantaneous value of the function and some other point: other_point-inst_val
137 : double difference(double)const;
138 : /// Calculate the difference between two values of this function: d2 -d1
139 : double difference(double d1,double d2)const;
140 : /// This returns the pointer to the action where this value is calculated
141 : ActionWithValue* getPntrToAction();
142 : /// Bring back one value into the correct pbc if needed, else give back the value
143 : double bringBackInPbc(double d1)const;
144 : /// Get the difference between max and minimum of domain
145 : double getMaxMinusMin()const;
146 : /// This sets up the gradients
147 : void setGradients();
148 : static double projection(const Value&,const Value&);
149 : };
150 :
151 : void copy( const Value& val1, Value& val2 );
152 : void copy( const Value& val1, Value* val2 );
153 : void add( const Value& val1, Value* valout );
154 :
155 : inline
156 41957264 : void Value::applyPeriodicity() {
157 41957264 : if(periodicity==periodic) {
158 30689730 : value=min+difference(min,value);
159 30688249 : if(value<min)value+=max_minus_min;
160 : }
161 41955783 : }
162 :
163 : inline
164 : void product( const Value& val1, const Value& val2, Value& valout ) {
165 : plumed_assert( val1.derivatives.size()==val2.derivatives.size() );
166 : if( valout.derivatives.size()!=val1.derivatives.size() ) valout.resizeDerivatives( val1.derivatives.size() );
167 : valout.value_set=false;
168 : valout.clearDerivatives();
169 : double u=val1.value;
170 : double v=val2.value;
171 : for(unsigned i=0; i<val1.derivatives.size(); ++i) {
172 : valout.addDerivative(i, u*val2.derivatives[i] + v*val1.derivatives[i] );
173 : }
174 : valout.set( u*v );
175 : }
176 :
177 : inline
178 : void quotient( const Value& val1, const Value& val2, Value* valout ) {
179 : plumed_assert( val1.derivatives.size()==val2.derivatives.size() );
180 : if( valout->derivatives.size()!=val1.derivatives.size() ) valout->resizeDerivatives( val1.derivatives.size() );
181 : valout->value_set=false;
182 : valout->clearDerivatives();
183 : double u=val1.get();
184 : double v=val2.get();
185 : for(unsigned i=0; i<val1.getNumberOfDerivatives(); ++i) {
186 : valout->addDerivative(i, v*val1.getDerivative(i) - u*val2.getDerivative(i) );
187 : }
188 : valout->chainRule( 1/(v*v) ); valout->set( u / v );
189 : }
190 :
191 : inline
192 : void Value::set(double v) {
193 41949529 : value_set=true;
194 41949529 : value=v;
195 41949529 : applyPeriodicity();
196 : }
197 :
198 : inline
199 : void Value::add(double v) {
200 : value_set=true;
201 : value+=v;
202 : applyPeriodicity();
203 : }
204 :
205 : inline
206 : double Value::get()const {
207 120040463 : return value;
208 : }
209 :
210 : inline
211 : bool Value::valueHasBeenSet() const {
212 0 : return value_set;
213 : }
214 :
215 : inline
216 : const std::string& Value::getName()const {
217 : return name;
218 : }
219 :
220 : inline
221 346896 : unsigned Value::getNumberOfDerivatives() const {
222 346896 : plumed_massert(hasDeriv,"the derivatives array for this value has zero size");
223 346896 : return derivatives.size();
224 : }
225 :
226 : inline
227 : double Value::getDerivative(const unsigned n) const {
228 : plumed_dbg_massert(n<derivatives.size(),"you are asking for a derivative that is out of bounds");
229 84545728 : return derivatives[n];
230 : }
231 :
232 : inline
233 : bool Value::hasDerivatives() const {
234 13040 : return hasDeriv;
235 : }
236 :
237 : inline
238 : void Value::resizeDerivatives(int n) {
239 4102750286 : if(hasDeriv) derivatives.resize(n);
240 : }
241 :
242 : inline
243 : void Value::addDerivative(unsigned i,double d) {
244 : plumed_dbg_massert(i<derivatives.size(),"derivative is out of bounds");
245 179306246 : derivatives[i]+=d;
246 : }
247 :
248 : inline
249 : void Value::setDerivative(unsigned i, double d) {
250 : plumed_dbg_massert(i<derivatives.size(),"derivative is out of bounds");
251 111698 : derivatives[i]=d;
252 : }
253 :
254 : inline
255 : void Value::chainRule(double df) {
256 36302 : for(unsigned i=0; i<derivatives.size(); ++i) derivatives[i]*=df;
257 : }
258 :
259 : inline
260 : void Value::clearInputForce() {
261 723853 : hasForce=false;
262 723853 : inputForce=0.0;
263 : }
264 :
265 : inline
266 : void Value::clearDerivatives() {
267 716132 : value_set=false;
268 : std::fill(derivatives.begin(), derivatives.end(), 0);
269 : }
270 :
271 : inline
272 : void Value::addForce(double f) {
273 : plumed_dbg_massert(hasDerivatives(),"forces can only be added to values with derivatives");
274 223499 : hasForce=true;
275 223499 : inputForce+=f;
276 : }
277 :
278 : inline
279 : double Value::getForce() const {
280 4794 : return inputForce;
281 : }
282 : /// d2-d1
283 : inline
284 164581527 : double Value::difference(double d1,double d2)const {
285 164581527 : if(periodicity==notperiodic) {
286 99183529 : return d2-d1;
287 65397998 : } else if(periodicity==periodic) {
288 65397998 : double s=(d2-d1)*inv_max_minus_min;
289 : // remember: pbc brings the difference in a range of -0.5:0.5
290 : s=Tools::pbc(s);
291 65397998 : return s*max_minus_min;
292 0 : } else plumed_merror("periodicity should be set to compute differences");
293 : }
294 :
295 : inline
296 : double Value::bringBackInPbc(double d1)const {
297 2373 : return min+max_minus_min/2.+difference(min+max_minus_min/2., d1);
298 : }
299 :
300 : inline
301 : double Value::difference(double d)const {
302 119233727 : return difference(get(),d);
303 : }
304 :
305 : inline
306 : double Value::getMaxMinusMin()const {
307 : plumed_dbg_assert( periodicity==periodic );
308 770 : return max_minus_min;
309 : }
310 :
311 : }
312 :
313 : #endif
314 :
|
