Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2018 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 : #include "core/ActionRegister.h"
23 : #include "Function.h"
24 :
25 : #include <cmath>
26 :
27 : using namespace std;
28 :
29 : namespace PLMD {
30 : namespace function {
31 :
32 : //+PLUMEDOC FUNCTION PIECEWISE
33 : /*
34 : Compute a piecewise straight line through its arguments that passes through
35 : a set of ordered control points.
36 :
37 : For variables less than the first
38 : (greater than the last) point, the value of the first (last) point is used.
39 :
40 : \f[
41 : \frac{y_{i+1}-y_i}{x_{i+1}-x_i}(s-x_i)+y_i ; if x_i<s<x_{i+1}
42 : \f]
43 : \f[
44 : y_N ; if x>x_{N-1}
45 : \f]
46 : \f[
47 : y_1 ; if x<x_0
48 : \f]
49 :
50 : Control points are passed using the POINT0=... POINT1=... syntax as in the example below
51 :
52 : If one argument is supplied, it results in a scalar quantity.
53 : If multiple arguments are supplied, it results
54 : in a vector of values. Each value will be named as the name of the original
55 : argument with suffix _pfunc.
56 :
57 : \par Examples
58 :
59 : \verbatim
60 : dist1: DISTANCE ATOMS=1,10
61 : dist2: DISTANCE ATOMS=2,11
62 :
63 : pw: PIECEWISE POINT0=1,10 POINT1=1,PI POINT2=3,10 ARG=dist1
64 : ppww: PIECEWISE POINT0=1,10 POINT1=1,PI POINT2=3,10 ARG=dist1,dist2
65 : PRINT ARG=pw,ppww.dist1_pfunc,ppww.dist2_pfunc
66 : \endverbatim
67 : (See also \ref PRINT and \ref DISTANCE).
68 :
69 :
70 : */
71 : //+ENDPLUMEDOC
72 :
73 :
74 4 : class Piecewise :
75 : public Function
76 : {
77 : std::vector<std::pair<double,double> > points;
78 : public:
79 : explicit Piecewise(const ActionOptions&);
80 : void calculate();
81 : static void registerKeywords(Keywords& keys);
82 : };
83 :
84 :
85 2525 : PLUMED_REGISTER_ACTION(Piecewise,"PIECEWISE")
86 :
87 3 : void Piecewise::registerKeywords(Keywords& keys) {
88 3 : Function::registerKeywords(keys);
89 3 : keys.use("ARG");
90 3 : keys.add("numbered","POINT","This keyword is used to specify the various points in the function above.");
91 3 : keys.reset_style("POINT","compulsory");
92 3 : componentsAreNotOptional(keys);
93 : keys.addOutputComponent("_pfunc","default","one or multiple instances of this quantity will be referenceable elsewhere "
94 : "in the input file. These quantities will be named with the arguments of the "
95 : "function followed by the character string _pfunc. These quantities tell the "
96 3 : "user the values of the piecewise functions of each of the arguments.");
97 3 : }
98 :
99 2 : Piecewise::Piecewise(const ActionOptions&ao):
100 : Action(ao),
101 2 : Function(ao)
102 : {
103 8 : for(int i=0;; i++) {
104 8 : std::vector<double> pp;
105 8 : if(!parseNumberedVector("POINT",i,pp) ) break;
106 6 : if(pp.size()!=2) error("points should be in x,y format");
107 6 : points.push_back(std::pair<double,double>(pp[0],pp[1]));
108 6 : if(i>0 && points[i].first<=points[i-1].first) error("points abscissas should be monotonously increasing");
109 12 : }
110 :
111 5 : for(unsigned i=0; i<getNumberOfArguments(); i++)
112 3 : if(getPntrToArgument(i)->isPeriodic())
113 0 : error("Cannot use PIECEWISE on periodic arguments");
114 :
115 2 : if(getNumberOfArguments()==1) {
116 1 : addValueWithDerivatives();
117 1 : setNotPeriodic();
118 : } else {
119 3 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
120 2 : addComponentWithDerivatives( getPntrToArgument(i)->getName()+"_pfunc" );
121 2 : getPntrToComponent(i)->setNotPeriodic();
122 : }
123 : }
124 2 : checkRead();
125 :
126 2 : log.printf(" on points:");
127 2 : for(unsigned i=0; i<points.size(); i++) log.printf(" (%f,%f)",points[i].first,points[i].second);
128 2 : log.printf("\n");
129 2 : }
130 :
131 10 : void Piecewise::calculate() {
132 25 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
133 15 : double val=getArgument(i);
134 15 : unsigned p=0;
135 37 : for(; p<points.size(); p++) {
136 33 : if(val<points[p].first) break;
137 : }
138 : double f,d;
139 15 : if(p==0) {
140 5 : f=points[0].second;
141 5 : d=0.0;
142 10 : } else if(p==points.size()) {
143 4 : f=points[points.size()-1].second;
144 4 : d=0.0;
145 : } else {
146 6 : double m=(points[p].second-points[p-1].second) / (points[p].first-points[p-1].first);
147 6 : f=m*(val-points[p-1].first)+points[p-1].second;
148 6 : d=m;
149 : }
150 15 : if(getNumberOfArguments()==1) {
151 5 : setValue(f);
152 5 : setDerivative(i,d);
153 : } else {
154 10 : Value* v=getPntrToComponent(i);
155 10 : v->set(f);
156 10 : v->addDerivative(i,d);
157 : }
158 : }
159 10 : }
160 :
161 : }
162 2523 : }
163 :
164 :
|
