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 : #include "ActionRegister.h"
23 : #include "Function.h"
24 : #include "lepton/Lepton.h"
25 :
26 : namespace PLMD {
27 : namespace function {
28 :
29 : //+PLUMEDOC FUNCTION CUSTOM
30 : /*
31 : Calculate a combination of variables using a custom expression.
32 :
33 : This action computes an arbitrary function of one or more
34 : collective variables. Arguments are chosen with the ARG keyword,
35 : and the function is provided with the FUNC string. Notice that this
36 : string should contain no space. Within FUNC, one can refer to the
37 : arguments as x,y,z, and t (up to four variables provided as ARG).
38 : This names can be customized using the VAR keyword (see examples below).
39 :
40 : This function is implemented using the Lepton library, that allows to evaluate
41 : algebraic expressions and to automatically differentiate them.
42 :
43 : If you want a function that depends not only on collective variables
44 : but also on time you can use the \subpage TIME action.
45 :
46 : \par Examples
47 :
48 : The following input tells plumed to perform a metadynamics
49 : using as a CV the difference between two distances.
50 : \plumedfile
51 : dAB: DISTANCE ATOMS=10,12
52 : dAC: DISTANCE ATOMS=10,15
53 : diff: CUSTOM ARG=dAB,dAC FUNC=y-x PERIODIC=NO
54 : # notice: the previous line could be replaced with the following
55 : # diff: COMBINE ARG=dAB,dAC COEFFICIENTS=-1,1
56 : METAD ARG=diff SIGMA=0.1 HEIGHT=0.5 BIASFACTOR=10 PACE=100
57 : \endplumedfile
58 : (see also \ref DISTANCE, \ref COMBINE, and \ref METAD).
59 : Notice that forces applied to diff will be correctly propagated
60 : to atoms 10, 12, and 15.
61 : Also notice that since CUSTOM is used without the VAR option
62 : the two arguments should be referred to as x and y in the expression FUNC.
63 : For simple functions
64 : such as this one it is possible to use \ref COMBINE.
65 :
66 : The following input tells plumed to print the angle between vectors
67 : identified by atoms 1,2 and atoms 2,3
68 : its square (as computed from the x,y,z components) and the distance
69 : again as computed from the square root of the square.
70 : \plumedfile
71 : DISTANCE LABEL=d1 ATOMS=1,2 COMPONENTS
72 : DISTANCE LABEL=d2 ATOMS=2,3 COMPONENTS
73 : CUSTOM ...
74 : LABEL=theta
75 : ARG=d1.x,d1.y,d1.z,d2.x,d2.y,d2.z
76 : VAR=ax,ay,az,bx,by,bz
77 : FUNC=acos((ax*bx+ay*by+az*bz)/sqrt((ax*ax+ay*ay+az*az)*(bx*bx+by*by+bz*bz)))
78 : PERIODIC=NO
79 : ... CUSTOM
80 : PRINT ARG=theta
81 : \endplumedfile
82 : (See also \ref PRINT and \ref DISTANCE).
83 :
84 : Notice that this action implements a large number of functions (trigonometric, exp, log, etc).
85 : Among the useful functions, have a look at the step function (that is the Heaviside function).
86 : `step(x)` is defined as 1 when `x` is positive and `0` when x is negative. This allows for
87 : a straightforward implementation of if clauses.
88 :
89 : For example, imagine that you want to implement a restraint that only acts when a
90 : distance is larger than 0.5. You can do it with
91 : \plumedfile
92 : d: DISTANCE ATOMS=10,15
93 : m: CUSTOM ARG=d FUNC=0.5*step(0.5-x)+x*step(x-0.5) PERIODIC=NO
94 : # check the function you are applying:
95 : PRINT ARG=d,m FILE=checkme
96 : RESTRAINT ARG=d AT=0.5 KAPPA=10.0
97 : \endplumedfile
98 : (see also \ref DISTANCE, \ref PRINT, and \ref RESTRAINT)
99 :
100 : The meaning of the function `0.5*step(0.5-x)+x*step(x-0.5)` is:
101 : - If x<0.5 (step(0.5-x)!=0) use 0.5
102 : - If x>0.5 (step(x-0.5)!=0) use x
103 : Notice that the same could have been obtained using an \ref UPPER_WALLS
104 : However, with CUSTOM you can create way more complex definitions.
105 :
106 : \warning If you apply forces on the variable (as in the previous example) you should
107 : make sure that the variable is continuous!
108 : Conversely, if you are just analyzing a trajectory you can safely use
109 : discontinuous variables.
110 :
111 : A possible continuity check with gnuplot is
112 : \verbatim
113 : # this allow to step function to be used in gnuplot:
114 : gnuplot> step(x)=0.5*(erf(x*10000000)+1)
115 : # here you can test your function
116 : gnuplot> p 0.5*step(0.5-x)+x*step(x-0.5)
117 : \endverbatim
118 :
119 : Also notice that you can easily make logical operations on the conditions that you
120 : create. The equivalent of the AND operator is the product: `step(1.0-x)*step(x-0.5)` is
121 : only equal to 1 when x is between 0.5 and 1.0. By combining negation and AND you can obtain an OR. That is,
122 : `1-step(1.0-x)*step(x-0.5)` is only equal to 1 when x is outside the 0.5-1.0 interval.
123 :
124 : CUSTOM can be used in combination with \ref DISTANCE to implement variants of the
125 : DISTANCE keyword that were present in PLUMED 1.3 and that allowed to compute
126 : the distance of a point from a line defined by two other points, or the progression
127 : along that line.
128 : \plumedfile
129 : # take center of atoms 1 to 10 as reference point 1
130 : p1: CENTER ATOMS=1-10
131 : # take center of atoms 11 to 20 as reference point 2
132 : p2: CENTER ATOMS=11-20
133 : # take center of atoms 21 to 30 as reference point 3
134 : p3: CENTER ATOMS=21-30
135 :
136 : # compute distances
137 : d12: DISTANCE ATOMS=p1,p2
138 : d13: DISTANCE ATOMS=p1,p3
139 : d23: DISTANCE ATOMS=p2,p3
140 :
141 : # compute progress variable of the projection of point p3
142 : # along the vector joining p1 and p2
143 : # notice that progress is measured from the middle point
144 : onaxis: CUSTOM ARG=d13,d23,d12 FUNC=(0.5*(y^2-x^2)/z) PERIODIC=NO
145 :
146 : # compute between point p3 and the vector joining p1 and p2
147 : fromaxis: CUSTOM ARG=d13,d23,d12,onaxis VAR=x,y,z,o FUNC=(0.5*(y^2+x^2)-o^2-0.25*z^2) PERIODIC=NO
148 :
149 : PRINT ARG=onaxis,fromaxis
150 :
151 : \endplumedfile
152 :
153 : Notice that these equations have been used to combine \ref RMSD
154 : from different snapshots of a protein so as to define
155 : progression (S) and distance (Z) variables \cite perez2015atp.
156 :
157 :
158 : */
159 : //+ENDPLUMEDOC
160 :
161 :
162 : class Custom :
163 : public Function {
164 : lepton::CompiledExpression expression;
165 : std::vector<lepton::CompiledExpression> expression_deriv;
166 : std::vector<std::string> var;
167 : std::string func;
168 : std::vector<double> values;
169 : std::vector<char*> names;
170 : std::vector<double*> lepton_ref;
171 : std::vector<double*> lepton_ref_deriv;
172 : public:
173 : explicit Custom(const ActionOptions&);
174 : void calculate() override;
175 : static void registerKeywords(Keywords& keys);
176 : };
177 :
178 13879 : PLUMED_REGISTER_ACTION(Custom,"CUSTOM")
179 :
180 : //+PLUMEDOC FUNCTION MATHEVAL
181 : /*
182 : An alias to the CUSTOM function that can also be used to calaculate combinations of variables using a custom expression.
183 :
184 : Documentation for this action is identical to that for \ref CUSTOM
185 :
186 : This alias is kept in order to maintain compatibility with previous PLUMED versions.
187 : However, notice that as of PLUMED 2.5 the libmatheval library is not linked anymore,
188 : and the \ref MATHEVAL function is implemented using the Lepton library.
189 :
190 : \par Examples
191 :
192 : Just replace \ref CUSTOM with \ref MATHEVAL.
193 :
194 : \plumedfile
195 : d: DISTANCE ATOMS=10,15
196 : m: MATHEVAL ARG=d FUNC=0.5*step(0.5-x)+x*step(x-0.5) PERIODIC=NO
197 : # check the function you are applying:
198 : PRINT ARG=d,m FILE=checkme
199 : RESTRAINT ARG=d AT=0.5 KAPPA=10.0
200 : \endplumedfile
201 : (see also \ref DISTANCE, \ref PRINT, and \ref RESTRAINT)
202 :
203 : */
204 : //+ENDPLUMEDOC
205 :
206 : class Matheval :
207 : public Custom {
208 : };
209 :
210 14587 : PLUMED_REGISTER_ACTION(Custom,"MATHEVAL")
211 :
212 456 : void Custom::registerKeywords(Keywords& keys) {
213 456 : Function::registerKeywords(keys);
214 456 : keys.use("ARG");
215 456 : keys.use("PERIODIC");
216 912 : keys.add("compulsory","FUNC","the function you wish to evaluate");
217 912 : keys.add("optional","VAR","the names to give each of the arguments in the function. If you have up to three arguments in your function you can use x, y and z to refer to them. Otherwise you must use this flag to give your variables names.");
218 456 : }
219 :
220 448 : Custom::Custom(const ActionOptions&ao):
221 : Action(ao),
222 : Function(ao),
223 448 : expression_deriv(getNumberOfArguments()),
224 448 : values(getNumberOfArguments()),
225 448 : names(getNumberOfArguments()),
226 448 : lepton_ref(getNumberOfArguments(),nullptr),
227 1344 : lepton_ref_deriv(getNumberOfArguments()*getNumberOfArguments(),nullptr) {
228 896 : parseVector("VAR",var);
229 448 : if(var.size()==0) {
230 435 : var.resize(getNumberOfArguments());
231 435 : if(getNumberOfArguments()>3) {
232 0 : error("Using more than 3 arguments you should explicitly write their names with VAR");
233 : }
234 435 : if(var.size()>0) {
235 : var[0]="x";
236 : }
237 435 : if(var.size()>1) {
238 : var[1]="y";
239 : }
240 435 : if(var.size()>2) {
241 : var[2]="z";
242 : }
243 : }
244 448 : if(var.size()!=getNumberOfArguments()) {
245 0 : error("Size of VAR array should be the same as number of arguments");
246 : }
247 448 : parse("FUNC",func);
248 448 : addValueWithDerivatives();
249 448 : checkRead();
250 :
251 448 : log.printf(" with function : %s\n",func.c_str());
252 448 : log.printf(" with variables :");
253 996 : for(unsigned i=0; i<var.size(); i++) {
254 548 : log.printf(" %s",var[i].c_str());
255 : }
256 448 : log.printf("\n");
257 :
258 448 : lepton::ParsedExpression pe=lepton::Parser::parse(func).optimize(lepton::Constants());
259 448 : log<<" function as parsed by lepton: "<<pe<<"\n";
260 448 : expression=pe.createCompiledExpression();
261 970 : for(auto &p: expression.getVariables()) {
262 522 : if(std::find(var.begin(),var.end(),p)==var.end()) {
263 0 : error("variable " + p + " is not defined");
264 : }
265 : }
266 448 : log<<" derivatives as computed by lepton:\n";
267 996 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
268 1096 : lepton::ParsedExpression pe=lepton::Parser::parse(func).differentiate(var[i]).optimize(lepton::Constants());
269 548 : log<<" "<<pe<<"\n";
270 548 : expression_deriv[i]=pe.createCompiledExpression();
271 : }
272 :
273 996 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
274 : try {
275 548 : lepton_ref[i]=&expression.getVariableReference(var[i]);
276 26 : } catch(const PLMD::lepton::Exception& exc) {
277 : // this is necessary since in some cases lepton things a variable is not present even though it is present
278 : // e.g. func=0*x
279 26 : }
280 : }
281 996 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
282 1554 : for(unsigned j=0; j<getNumberOfArguments(); j++) {
283 : try {
284 1006 : lepton_ref_deriv[i*getNumberOfArguments()+j]=&expression_deriv[i].getVariableReference(var[j]);
285 513 : } catch(const PLMD::lepton::Exception& exc) {
286 : // this is necessary since in some cases lepton things a variable is not present even though it is present
287 : // e.g. func=0*x
288 513 : }
289 : }
290 : }
291 448 : }
292 :
293 40568 : void Custom::calculate() {
294 81777 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
295 41209 : if(lepton_ref[i]) {
296 40119 : *lepton_ref[i]=getArgument(i);
297 : }
298 : }
299 40568 : setValue(expression.evaluate());
300 81777 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
301 85176 : for(unsigned j=0; j<getNumberOfArguments(); j++) {
302 43967 : if(lepton_ref_deriv[i*getNumberOfArguments()+j]) {
303 34220 : *lepton_ref_deriv[i*getNumberOfArguments()+j]=getArgument(j);
304 : }
305 : }
306 41209 : setDerivative(i,expression_deriv[i].evaluate());
307 : }
308 40568 : }
309 :
310 : }
311 : }
312 :
313 :
|
