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 "ActionWithArguments.h"
23 : #include "ActionWithValue.h"
24 : #include "ActionAtomistic.h"
25 : #include "ActionForInterface.h"
26 : #include "ActionWithVector.h"
27 : #include "ActionWithVirtualAtom.h"
28 : #include "ActionShortcut.h"
29 : #include "tools/PDB.h"
30 : #include "PlumedMain.h"
31 : #include "ActionSet.h"
32 : #include <iostream>
33 : #include <regex>
34 :
35 : namespace PLMD {
36 :
37 18411 : void ActionWithArguments::registerKeywords(Keywords& keys) {
38 36822 : keys.reserve("numbered","ARG","the input for this action is the scalar output from one or more other actions. The particular scalars that you will use "
39 : "are referenced using the label of the action. If the label appears on its own then it is assumed that the Action calculates "
40 : "a single scalar value. The value of this scalar is thus used as the input to this new action. If * or *.* appears the "
41 : "scalars calculated by all the proceeding actions in the input file are taken. Some actions have multi-component outputs and "
42 : "each component of the output has a specific label. For example a \\ref DISTANCE action labelled dist may have three components "
43 : "x, y and z. To take just the x component you should use dist.x, if you wish to take all three components then use dist.*."
44 : "More information on the referencing of Actions can be found in the section of the manual on the PLUMED \\ref Syntax. "
45 : "Scalar values can also be "
46 : "referenced using POSIX regular expressions as detailed in the section on \\ref Regex. To use this feature you you must compile "
47 : "PLUMED with the appropriate flag.");
48 18411 : }
49 :
50 10021 : void ActionWithArguments::parseArgumentList(const std::string&key,std::vector<Value*>&arg) {
51 : std::string def;
52 : std::vector<std::string> c;
53 : arg.clear();
54 10021 : parseVector(key,c);
55 10970 : if( c.size()==0 && (keywords.style(key,"compulsory") || keywords.style(key,"hidden")) ) {
56 7 : if( keywords.getDefaultValue(key,def) ) {
57 0 : c.push_back( def );
58 : } else {
59 : return;
60 : }
61 : }
62 10014 : interpretArgumentList(c,plumed.getActionSet(),this,arg);
63 10021 : }
64 :
65 104 : bool ActionWithArguments::parseArgumentList(const std::string&key,int i,std::vector<Value*>&arg) {
66 : std::vector<std::string> c;
67 : arg.clear();
68 104 : if(parseNumberedVector(key,i,c)) {
69 44 : interpretArgumentList(c,plumed.getActionSet(),this,arg);
70 : return true;
71 : } else {
72 : return false;
73 : }
74 104 : }
75 :
76 15387 : void ActionWithArguments::interpretArgumentList(const std::vector<std::string>& c, const ActionSet& as, Action* readact, std::vector<Value*>&arg) {
77 41261 : for(unsigned i=0; i<c.size(); i++) {
78 : // is a regex? then just interpret it. The signal is ()
79 25877 : if(!c[i].compare(0,1,"(")) {
80 219 : unsigned l=c[i].length();
81 219 : if(!c[i].compare(l-1,1,")")) {
82 : // start regex parsing
83 : bool found_something=false;
84 : // take the string enclosed in quotes and put in round brackets
85 218 : std::string myregex=c[i];
86 218 : std::vector<ActionWithValue*> all=as.select<ActionWithValue*>();
87 218 : if( all.empty() ) {
88 1 : readact->error("your input file is not telling plumed to calculate anything");
89 : }
90 :
91 : try {
92 218 : std::regex txt_regex(myregex,std::regex::extended);
93 217 : plumed_massert(txt_regex.mark_count()==1,"I can parse with only one subexpression");
94 24239 : for(unsigned j=0; j<all.size(); j++) {
95 24022 : std::vector<std::string> ss=all[j]->getComponentsVector();
96 349707 : for(unsigned k=0; k<ss.size(); ++k) {
97 325685 : if(std::regex_match(ss[k],txt_regex)) {
98 21601 : arg.push_back(all[j]->copyOutput(ss[k]));
99 : found_something=true;
100 : }
101 : }
102 24022 : }
103 218 : } catch(std::regex_error & e) {
104 3 : plumed_error()<<"Error parsing regular expression: "<<e.what();
105 1 : }
106 217 : if(!found_something) {
107 0 : plumed_error()<<"There isn't any action matching your regex " << myregex;
108 : }
109 : } else {
110 2 : plumed_merror("did you want to use regexp to input arguments? enclose it between two round braces (...) with no spaces!");
111 : }
112 : } else {
113 : std::size_t dot=c[i].find_first_of('.');
114 25658 : std::string a=c[i].substr(0,dot);
115 25658 : std::string name=c[i].substr(dot+1);
116 25658 : if(c[i].find(".")!=std::string::npos) { // if it contains a dot:
117 6534 : if(a=="*" && name=="*") {
118 : // Take all values from all actions
119 1 : std::vector<ActionWithValue*> all=as.select<ActionWithValue*>();
120 1 : if( all.empty() ) {
121 0 : readact->error("your input file is not telling plumed to calculate anything");
122 : }
123 17 : for(unsigned j=0; j<all.size(); j++) {
124 16 : plumed_assert(all[j]); // needed for following calls, see #1046
125 16 : ActionForInterface* ap=all[j]->castToActionForInterface();
126 16 : if( ap ) {
127 8 : continue;
128 : }
129 18 : for(int k=0; k<all[j]->getNumberOfComponents(); ++k) {
130 10 : arg.push_back(all[j]->copyOutput(k));
131 : }
132 : }
133 6515 : } else if ( name=="*") {
134 : unsigned carg=arg.size();
135 : // Take all the values from an action with a specific name
136 650 : ActionShortcut* shortcut=as.getShortcutActionWithLabel(a);
137 650 : if( shortcut ) {
138 822 : shortcut->interpretDataLabel( a + "." + name, readact, arg );
139 : }
140 650 : if( arg.size()==carg ) {
141 : // Take all the values from an action with a specific name
142 395 : ActionWithValue* action=as.selectWithLabel<ActionWithValue*>(a);
143 395 : if(!action) {
144 0 : std::string str=" (hint! the actions with value in this ActionSet are: ";
145 0 : str+=as.getLabelList<ActionWithValue*>()+")";
146 0 : readact->error("cannot find action named " + a + str);
147 : }
148 395 : if( action->getNumberOfComponents()==0 ) {
149 0 : readact->error("found " + a +".* indicating use all components calculated by action with label " + a + " but this action has no components");
150 : }
151 6306 : for(int k=0; k<action->getNumberOfComponents(); ++k) {
152 5911 : arg.push_back(action->copyOutput(k));
153 : }
154 : }
155 5865 : } else if ( a=="*" ) {
156 17 : std::vector<ActionShortcut*> shortcuts=as.select<ActionShortcut*>();
157 : // Take components from all actions with a specific name
158 17 : std::vector<ActionWithValue*> all=as.select<ActionWithValue*>();
159 17 : if( all.empty() ) {
160 0 : readact->error("your input file is not telling plumed to calculate anything");
161 : }
162 : unsigned carg=arg.size();
163 161 : for(unsigned j=0; j<shortcuts.size(); ++j) {
164 288 : shortcuts[j]->interpretDataLabel( shortcuts[j]->getShortcutLabel() + "." + name, readact, arg );
165 : }
166 : unsigned nval=0;
167 320 : for(unsigned j=0; j<all.size(); j++) {
168 : std::string flab;
169 606 : flab=all[j]->getLabel() + "." + name;
170 303 : if( all[j]->exists(flab) ) {
171 44 : arg.push_back(all[j]->copyOutput(flab));
172 44 : nval++;
173 : }
174 : }
175 17 : if(nval==0 && arg.size()==carg) {
176 0 : readact->error("found no actions with a component called " + name );
177 : }
178 : } else {
179 : // Take values with a specific name
180 5848 : ActionWithValue* action=as.selectWithLabel<ActionWithValue*>(a);
181 5848 : ActionShortcut* shortcut=as.getShortcutActionWithLabel(a);
182 5848 : if( !shortcut && !action ) {
183 0 : std::string str=" (hint! the actions with value in this ActionSet are: ";
184 0 : str+=as.getLabelList<ActionWithValue*>()+")";
185 0 : readact->error("cannot find action named " + a +str);
186 5848 : } else if( action && action->exists(c[i]) ) {
187 5760 : arg.push_back(action->copyOutput(c[i]));
188 88 : } else if( shortcut ) {
189 : unsigned narg=arg.size();
190 176 : shortcut->interpretDataLabel( a + "." + name, readact, arg );
191 88 : if( arg.size()==narg ) {
192 0 : readact->error("found no element in " + a + " with label " + name );
193 : }
194 : } else {
195 0 : std::string str=" (hint! the components in this actions are: ";
196 0 : str+=action->getComponentsList()+")";
197 0 : readact->error("action " + a + " has no component named " + name + str);
198 : }
199 : }
200 : } else { // if it doesn't contain a dot
201 19142 : if(c[i]=="*") {
202 : // Take all values from all actions
203 107 : std::vector<ActionWithValue*> all=as.select<ActionWithValue*>();
204 107 : if( all.empty() ) {
205 0 : readact->error("your input file is not telling plumed to calculate anything");
206 : }
207 1615 : for(unsigned j=0; j<all.size(); j++) {
208 1508 : plumed_assert(all[j]); // needed for following calls, see #1046
209 1508 : ActionWithVirtualAtom* av=all[j]->castToActionWithVirtualAtom();
210 1508 : if( av ) {
211 57 : continue;
212 : }
213 1451 : ActionForInterface* ap=all[j]->castToActionForInterface();
214 1451 : if( ap && all[j]->getName()!="ENERGY" ) {
215 834 : continue;
216 : }
217 1320 : for(int k=0; k<all[j]->getNumberOfComponents(); ++k) {
218 703 : arg.push_back(all[j]->copyOutput(k));
219 : }
220 : }
221 : } else {
222 19035 : ActionWithValue* action=as.selectWithLabel<ActionWithValue*>(c[i]);
223 19035 : if(!action) {
224 1 : std::string str=" (hint! the actions with value in this ActionSet are: ";
225 2 : str+=as.getLabelList<ActionWithValue*>()+")";
226 3 : readact->error("cannot find action named " + c[i] + str );
227 : }
228 19034 : if( !(action->exists(c[i])) ) {
229 0 : std::string str=" (hint! the components in this actions are: ";
230 0 : str+=action->getComponentsList()+")";
231 0 : readact->error("action " + c[i] + " has no component named " + c[i] +str);
232 : };
233 19035 : arg.push_back(action->copyOutput(c[i]));
234 : }
235 : }
236 : }
237 : }
238 15384 : }
239 :
240 0 : void ActionWithArguments::expandArgKeywordInPDB( const PDB& pdb ) {
241 0 : std::vector<std::string> arg_names = pdb.getArgumentNames();
242 0 : if( arg_names.size()>0 ) {
243 : std::vector<Value*> arg_vals;
244 0 : interpretArgumentList( arg_names, plumed.getActionSet(), this, arg_vals );
245 : }
246 0 : }
247 :
248 186077 : void ActionWithArguments::requestArguments(const std::vector<Value*> &arg) {
249 186077 : plumed_massert(!lockRequestArguments,"requested argument list can only be changed in the prepare() method");
250 186077 : arguments=arg;
251 186077 : clearDependencies();
252 : std::string fullname;
253 : std::string name;
254 1276136 : for(unsigned i=0; i<arguments.size(); i++) {
255 1090059 : fullname=arguments[i]->getName();
256 1090059 : if(fullname.find(".")!=std::string::npos) {
257 : std::size_t dot=fullname.find_first_of('.');
258 746184 : name=fullname.substr(0,dot);
259 : } else {
260 : name=fullname;
261 : }
262 1090059 : ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(name);
263 1090059 : plumed_massert(action,"cannot find action named (in requestArguments - this is weird)" + name);
264 1090059 : addDependency(action);
265 : }
266 186077 : ActionWithValue* av=dynamic_cast<ActionWithValue*>(this);
267 186077 : if(av) {
268 176400 : av->firststep=true;
269 : }
270 186077 : }
271 :
272 4 : void ActionWithArguments::requestExtraDependencies(const std::vector<Value*> &extra) {
273 4 : plumed_massert(!lockRequestArguments,"requested argument list can only be changed in the prepare() method");
274 : std::string fullname;
275 : std::string name;
276 9 : for(unsigned i=0; i<extra.size(); i++) {
277 5 : fullname=extra[i]->getName();
278 5 : if(fullname.find(".")!=std::string::npos) {
279 : std::size_t dot=fullname.find_first_of('.');
280 0 : name=fullname.substr(0,dot);
281 : } else {
282 : name=fullname;
283 : }
284 5 : ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(name);
285 5 : plumed_massert(action,"cannot find action named (in requestArguments - this is weird)" + name);
286 5 : addDependency(action);
287 : }
288 4 : }
289 :
290 10345 : ActionWithArguments::ActionWithArguments(const ActionOptions&ao):
291 : Action(ao),
292 10345 : lockRequestArguments(false) {
293 20690 : if( keywords.exists("ARG") ) {
294 : std::vector<Value*> arg;
295 19332 : parseArgumentList("ARG",arg);
296 :
297 9666 : if(!arg.empty()) {
298 9328 : log.printf(" with arguments : \n");
299 45663 : for(unsigned i=0; i<arg.size(); i++) {
300 36335 : if( arg[i]->hasDerivatives() && arg[i]->getRank()>0 ) {
301 1009 : log.printf(" function on grid with label %s \n",arg[i]->getName().c_str());
302 35326 : } else if( arg[i]->getRank()==2 ) {
303 2611 : log.printf(" matrix with label %s \n",arg[i]->getName().c_str());
304 32715 : } else if( arg[i]->getRank()==1 ) {
305 6135 : log.printf(" vector with label %s \n",arg[i]->getName().c_str());
306 26580 : } else if( arg[i]->getRank()==0 ) {
307 26580 : log.printf(" scalar with label %s \n",arg[i]->getName().c_str());
308 : } else {
309 0 : error("type of argument does not make sense");
310 : }
311 : }
312 : }
313 9666 : requestArguments(arg);
314 : }
315 10345 : }
316 :
317 58 : void ActionWithArguments::calculateNumericalDerivatives( ActionWithValue* a ) {
318 58 : if(!a) {
319 58 : a=castToActionWithValue();
320 58 : plumed_massert(a,"cannot compute numerical derivatives for an action without values");
321 : }
322 :
323 58 : const size_t nval=a->getNumberOfComponents();
324 : const size_t npar=arguments.size();
325 58 : std::vector<double> value (nval*npar);
326 161 : for(int i=0; i<npar; i++) {
327 103 : double arg0=arguments[i]->get();
328 103 : arguments[i]->set(arg0+std::sqrt(epsilon));
329 103 : a->calculate();
330 103 : arguments[i]->set(arg0);
331 1367 : for(int j=0; j<nval; j++) {
332 1264 : value[i*nval+j]=a->getOutputQuantity(j);
333 : }
334 : }
335 58 : a->calculate();
336 58 : a->clearDerivatives();
337 1192 : for(int j=0; j<nval; j++) {
338 1134 : Value* v=a->copyOutput(j);
339 1134 : if( v->hasDerivatives() )
340 670 : for(int i=0; i<npar; i++) {
341 400 : v->addDerivative(i,(value[i*nval+j]-a->getOutputQuantity(j))/std::sqrt(epsilon));
342 : }
343 : }
344 58 : }
345 :
346 261 : double ActionWithArguments::getProjection(unsigned i,unsigned j)const {
347 261 : plumed_massert(i<arguments.size()," making projections with an index which is too large");
348 261 : plumed_massert(j<arguments.size()," making projections with an index which is too large");
349 261 : const Value* v1=arguments[i];
350 261 : const Value* v2=arguments[j];
351 261 : return Value::projection(*v1,*v2);
352 : }
353 :
354 201798 : void ActionWithArguments::addForcesOnArguments( const unsigned& argstart, const std::vector<double>& forces, unsigned& ind, const std::string& c ) {
355 529572 : for(unsigned i=0; i<arguments.size(); ++i) {
356 327774 : if( i==0 && getName().find("EVALUATE_FUNCTION_FROM_GRID")!=std::string::npos ) {
357 3350 : continue ;
358 : }
359 324424 : if( !arguments[i]->ignoreStoredValue(c) || arguments[i]->getRank()==0 || (arguments[i]->getRank()>0 && arguments[i]->hasDerivatives()) ) {
360 302069 : unsigned nvals = arguments[i]->getNumberOfStoredValues();
361 33595635 : for(unsigned j=0; j<nvals; ++j) {
362 33293566 : arguments[i]->addForce( j, forces[ind], false );
363 33293566 : ind++;
364 : }
365 : }
366 : }
367 201798 : }
368 :
369 83 : void ActionWithArguments::setGradients( Value* myval, unsigned& start ) const {
370 83 : if( !myval->hasDeriv ) {
371 : return;
372 : }
373 83 : plumed_assert( myval->getRank()==0 );
374 :
375 : bool scalar=true;
376 249 : for(unsigned i=0; i<arguments.size(); ++i ) {
377 166 : if( arguments[i]->getRank()!=0 ) {
378 : scalar=false;
379 : break;
380 : }
381 : }
382 83 : if( !scalar ) {
383 : bool constant=true;
384 0 : for(unsigned i=0; i<arguments.size(); ++i ) {
385 0 : if( !arguments[i]->isConstant() ) {
386 : constant=false;
387 : break;
388 : } else {
389 0 : start += arguments[i]->getNumberOfValues();
390 : }
391 : }
392 0 : if( !constant ) {
393 0 : error("cannot set gradient as unable to handle non-constant actions that take vectors/matrices/grids in input");
394 : }
395 : }
396 : // Now pass the gradients
397 249 : for(unsigned i=0; i<arguments.size(); ++i ) {
398 166 : arguments[i]->passGradients( myval->getDerivative(i), myval->gradients );
399 : }
400 : }
401 :
402 18273 : bool ActionWithArguments::calculateConstantValues( const bool& haveatoms ) {
403 18273 : ActionWithValue* av = castToActionWithValue();
404 18273 : if( !av || arguments.size()==0 ) {
405 : return false;
406 : }
407 : bool constant = true, atoms=false;
408 14792 : for(unsigned i=0; i<arguments.size(); ++i) {
409 13680 : auto * ptr=arguments[i]->getPntrToAction();
410 13680 : plumed_assert(ptr); // needed for following calls, see #1046
411 13680 : ActionAtomistic* aa=ptr->castToActionAtomistic();
412 13680 : if( aa ) {
413 10796 : ActionWithVector* av=dynamic_cast<ActionWithVector*>( arguments[i]->getPntrToAction() );
414 10796 : if( !av || aa->getNumberOfAtoms()>0 ) {
415 : atoms=true;
416 : }
417 : }
418 13680 : if( !arguments[i]->isConstant() ) {
419 : constant=false;
420 : break;
421 : }
422 : }
423 13432 : if( constant ) {
424 : // Set everything constant first as we need to set the shape
425 2248 : for(unsigned i=0; i<av->getNumberOfComponents(); ++i) {
426 1136 : (av->copyOutput(i))->setConstant();
427 : }
428 1112 : if( !haveatoms ) {
429 1094 : log.printf(" values stored by this action are computed during startup and stay fixed during the simulation\n");
430 : }
431 1112 : if( atoms ) {
432 36 : return haveatoms;
433 : }
434 : }
435 : // Now do the calculation and store the values if we don't need anything from the atoms
436 13396 : if( constant && !haveatoms ) {
437 1076 : plumed_assert( !atoms );
438 1076 : activate();
439 1076 : calculate();
440 1076 : deactivate();
441 2176 : for(unsigned i=0; i<av->getNumberOfComponents(); ++i) {
442 1100 : unsigned nv = av->copyOutput(i)->getNumberOfValues();
443 1100 : log.printf(" %d values stored in component labelled %s are : ", nv, (av->copyOutput(i))->getName().c_str() );
444 3939 : for(unsigned j=0; j<nv; ++j) {
445 2839 : log.printf(" %f", (av->copyOutput(i))->get(j) );
446 : }
447 1100 : log.printf("\n");
448 : }
449 : }
450 : return constant;
451 : }
452 :
453 : }
|
