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 2455190 : 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 : /// This sets up the gradients
145 : void setGradients();
146 : static double projection(const Value&,const Value&);
147 : };
148 :
149 : void copy( const Value& val1, Value& val2 );
150 : void copy( const Value& val1, Value* val2 );
151 : void add( const Value& val1, Value* valout );
152 :
153 : inline
154 39215596 : void Value::applyPeriodicity() {
155 39215596 : if(periodicity==periodic) {
156 30688010 : value=min+difference(min,value);
157 30688010 : if(value<min)value+=max_minus_min;
158 : }
159 39215596 : }
160 :
161 : inline
162 : void product( const Value& val1, const Value& val2, Value& valout ) {
163 : plumed_assert( val1.derivatives.size()==val2.derivatives.size() );
164 : if( valout.derivatives.size()!=val1.derivatives.size() ) valout.resizeDerivatives( val1.derivatives.size() );
165 : valout.value_set=false;
166 : valout.clearDerivatives();
167 : double u=val1.value;
168 : double v=val2.value;
169 : for(unsigned i=0; i<val1.derivatives.size(); ++i) {
170 : valout.addDerivative(i, u*val2.derivatives[i] + v*val1.derivatives[i] );
171 : }
172 : valout.set( u*v );
173 : }
174 :
175 : inline
176 : void quotient( const Value& val1, const Value& val2, Value* valout ) {
177 : plumed_assert( val1.derivatives.size()==val2.derivatives.size() );
178 : if( valout->derivatives.size()!=val1.derivatives.size() ) valout->resizeDerivatives( val1.derivatives.size() );
179 : valout->value_set=false;
180 : valout->clearDerivatives();
181 : double u=val1.get();
182 : double v=val2.get();
183 : for(unsigned i=0; i<val1.getNumberOfDerivatives(); ++i) {
184 : valout->addDerivative(i, v*val1.getDerivative(i) - u*val2.getDerivative(i) );
185 : }
186 : valout->chainRule( 1/(v*v) ); valout->set( u / v );
187 : }
188 :
189 : inline
190 : void Value::set(double v) {
191 39089215 : value_set=true;
192 39215596 : value=v;
193 39215596 : applyPeriodicity();
194 : }
195 :
196 : inline
197 : void Value::add(double v) {
198 : value_set=true;
199 : value+=v;
200 : applyPeriodicity();
201 : }
202 :
203 : inline
204 : double Value::get()const {
205 42061876 : return value;
206 : }
207 :
208 : inline
209 : bool Value::valueHasBeenSet() const {
210 0 : return value_set;
211 : }
212 :
213 : inline
214 : const std::string& Value::getName()const {
215 3245018 : return name;
216 : }
217 :
218 : inline
219 338586 : unsigned Value::getNumberOfDerivatives() const {
220 338586 : plumed_massert(hasDeriv,"the derivatives array for this value has zero size");
221 338586 : return derivatives.size();
222 : }
223 :
224 : inline
225 : double Value::getDerivative(const unsigned n) const {
226 : plumed_dbg_massert(n<derivatives.size(),"you are asking for a derivative that is out of bounds");
227 82785980 : return derivatives[n];
228 : }
229 :
230 : inline
231 : bool Value::hasDerivatives() const {
232 13030 : return hasDeriv;
233 : }
234 :
235 : inline
236 : void Value::resizeDerivatives(int n) {
237 4102745976 : if(hasDeriv) derivatives.resize(n);
238 : }
239 :
240 : inline
241 : void Value::addDerivative(unsigned i,double d) {
242 : plumed_dbg_massert(i<derivatives.size(),"derivative is out of bounds");
243 129695818 : derivatives[i]+=d;
244 : }
245 :
246 : inline
247 : void Value::setDerivative(unsigned i, double d) {
248 : plumed_dbg_massert(i<derivatives.size(),"derivative is out of bounds");
249 110981 : derivatives[i]=d;
250 : }
251 :
252 : inline
253 : void Value::chainRule(double df) {
254 48790 : for(unsigned i=0; i<derivatives.size(); ++i) derivatives[i]*=df;
255 : }
256 :
257 : inline
258 : void Value::clearInputForce() {
259 426184 : hasForce=false;
260 426184 : inputForce=0.0;
261 : }
262 :
263 : inline
264 : void Value::clearDerivatives() {
265 433601 : value_set=false;
266 867202 : std::fill(derivatives.begin(), derivatives.end(), 0);
267 : }
268 :
269 : inline
270 : void Value::addForce(double f) {
271 : plumed_dbg_massert(hasDerivatives(),"forces can only be added to values with derivatives");
272 118645 : hasForce=true;
273 118645 : inputForce+=f;
274 : }
275 :
276 : inline
277 : double Value::getForce() const {
278 4698 : return inputForce;
279 : }
280 : /// d2-d1
281 : inline
282 86277149 : double Value::difference(double d1,double d2)const {
283 86277149 : if(periodicity==notperiodic) {
284 21544706 : return d2-d1;
285 64732443 : } else if(periodicity==periodic) {
286 64732443 : double s=(d2-d1)*inv_max_minus_min;
287 : // remember: pbc brings the difference in a range of -0.5:0.5
288 : s=Tools::pbc(s);
289 64732443 : return s*max_minus_min;
290 0 : } else plumed_merror("periodicity should be set to compute differences");
291 : }
292 :
293 : inline
294 : double Value::bringBackInPbc(double d1)const {
295 690 : return min+max_minus_min/2.+difference(min+max_minus_min/2., d1);
296 : }
297 :
298 : inline
299 : double Value::difference(double d)const {
300 41605998 : return difference(get(),d);
301 : }
302 :
303 : }
304 :
305 : #endif
306 :
|
