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 18017 : void ActionWithArguments::registerKeywords(Keywords& keys) {
38 : // 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 18017 : }
49 :
50 9993 : void ActionWithArguments::parseArgumentList(const std::string&key,std::vector<Value*>&arg) {
51 19986 : if( keywords.getArgumentType(key).length()==0 ) {
52 48 : warning("keyword " + key + " for reading arguments is registered using Keyword::add rather than Keyword::addInputKeyword. The keyword will thus not appear in the correct place in the manual");
53 : }
54 : std::string def;
55 : std::vector<std::string> c;
56 : arg.clear();
57 9993 : parseVector(key,c);
58 10938 : if( c.size()==0 && (keywords.style(key,"compulsory") || keywords.style(key,"hidden")) ) {
59 7 : if( keywords.getDefaultValue(key,def) ) {
60 0 : c.push_back( def );
61 : } else {
62 : return;
63 : }
64 : }
65 9986 : interpretArgumentList(c,plumed.getActionSet(),this,arg);
66 9993 : }
67 :
68 104 : bool ActionWithArguments::parseArgumentList(const std::string&key,int i,std::vector<Value*>&arg) {
69 208 : if( keywords.getArgumentType(key).length()==0 ) {
70 0 : warning("keyword " + key + " for reading argument is registered using Keyword::add rather than Keyword::addInputKeyword. The keyword will thus not appear in the correct place in the manual");
71 : }
72 : std::vector<std::string> c;
73 : arg.clear();
74 104 : if(parseNumberedVector(key,i,c)) {
75 44 : interpretArgumentList(c,plumed.getActionSet(),this,arg);
76 : return true;
77 : } else {
78 : return false;
79 : }
80 104 : }
81 :
82 15345 : void ActionWithArguments::interpretArgumentList(const std::vector<std::string>& c, const ActionSet& as, Action* readact, std::vector<Value*>&arg) {
83 41173 : for(unsigned i=0; i<c.size(); i++) {
84 : // is a regex? then just interpret it. The signal is ()
85 25831 : if(!c[i].compare(0,1,"(")) {
86 219 : unsigned l=c[i].length();
87 219 : if(!c[i].compare(l-1,1,")")) {
88 : // start regex parsing
89 : bool found_something=false;
90 : // take the string enclosed in quotes and put in round brackets
91 218 : std::string myregex=c[i];
92 218 : std::vector<ActionWithValue*> all=as.select<ActionWithValue*>();
93 218 : if( all.empty() ) {
94 1 : readact->error("your input file is not telling plumed to calculate anything");
95 : }
96 :
97 : try {
98 218 : std::regex txt_regex(myregex,std::regex::extended);
99 217 : plumed_massert(txt_regex.mark_count()==1,"I can parse with only one subexpression");
100 24240 : for(unsigned j=0; j<all.size(); j++) {
101 24023 : std::vector<std::string> ss=all[j]->getComponentsVector();
102 349709 : for(unsigned k=0; k<ss.size(); ++k) {
103 325686 : if(std::regex_match(ss[k],txt_regex)) {
104 21602 : arg.push_back(all[j]->copyOutput(ss[k]));
105 : found_something=true;
106 : }
107 : }
108 24023 : }
109 218 : } catch(std::regex_error & e) {
110 3 : plumed_error()<<"Error parsing regular expression: "<<e.what();
111 1 : }
112 217 : if(!found_something) {
113 0 : plumed_error()<<"There isn't any action matching your regex " << myregex;
114 : }
115 : } else {
116 2 : plumed_merror("did you want to use regexp to input arguments? enclose it between two round braces (...) with no spaces!");
117 : }
118 : } else {
119 : std::size_t dot=c[i].find_first_of('.');
120 25612 : std::string a=c[i].substr(0,dot);
121 25612 : std::string name=c[i].substr(dot+1);
122 25612 : if(c[i].find(".")!=std::string::npos) { // if it contains a dot:
123 6532 : if(a=="*" && name=="*") {
124 : // Take all values from all actions
125 1 : std::vector<ActionWithValue*> all=as.select<ActionWithValue*>();
126 1 : if( all.empty() ) {
127 0 : readact->error("your input file is not telling plumed to calculate anything");
128 : }
129 17 : for(unsigned j=0; j<all.size(); j++) {
130 16 : plumed_assert(all[j]); // needed for following calls, see #1046
131 16 : ActionForInterface* ap=all[j]->castToActionForInterface();
132 16 : if( ap ) {
133 8 : continue;
134 : }
135 18 : for(int k=0; k<all[j]->getNumberOfComponents(); ++k) {
136 10 : arg.push_back(all[j]->copyOutput(k));
137 : }
138 : }
139 6513 : } else if ( name=="*") {
140 : unsigned carg=arg.size();
141 : // Take all the values from an action with a specific name
142 648 : ActionShortcut* shortcut=as.getShortcutActionWithLabel(a);
143 648 : if( shortcut ) {
144 818 : shortcut->interpretDataLabel( a + "." + name, readact, arg );
145 : }
146 648 : if( arg.size()==carg ) {
147 : // Take all the values from an action with a specific name
148 393 : ActionWithValue* action=as.selectWithLabel<ActionWithValue*>(a);
149 393 : if(!action) {
150 0 : std::string str=" (hint! the actions with value in this ActionSet are: ";
151 0 : str+=as.getLabelList<ActionWithValue*>()+")";
152 0 : readact->error("cannot find action named " + a + str);
153 : }
154 393 : if( action->getNumberOfComponents()==0 ) {
155 0 : readact->error("found " + a +".* indicating use all components calculated by action with label " + a + " but this action has no components");
156 : }
157 6298 : for(int k=0; k<action->getNumberOfComponents(); ++k) {
158 5905 : arg.push_back(action->copyOutput(k));
159 : }
160 : }
161 5865 : } else if ( a=="*" ) {
162 17 : std::vector<ActionShortcut*> shortcuts=as.select<ActionShortcut*>();
163 : // Take components from all actions with a specific name
164 17 : std::vector<ActionWithValue*> all=as.select<ActionWithValue*>();
165 17 : if( all.empty() ) {
166 0 : readact->error("your input file is not telling plumed to calculate anything");
167 : }
168 : unsigned carg=arg.size();
169 161 : for(unsigned j=0; j<shortcuts.size(); ++j) {
170 288 : shortcuts[j]->interpretDataLabel( shortcuts[j]->getShortcutLabel() + "." + name, readact, arg );
171 : }
172 : unsigned nval=0;
173 320 : for(unsigned j=0; j<all.size(); j++) {
174 : std::string flab;
175 606 : flab=all[j]->getLabel() + "." + name;
176 303 : if( all[j]->exists(flab) ) {
177 44 : arg.push_back(all[j]->copyOutput(flab));
178 44 : nval++;
179 : }
180 : }
181 17 : if(nval==0 && arg.size()==carg) {
182 0 : readact->error("found no actions with a component called " + name );
183 : }
184 : } else {
185 : // Take values with a specific name
186 5848 : ActionWithValue* action=as.selectWithLabel<ActionWithValue*>(a);
187 5848 : ActionShortcut* shortcut=as.getShortcutActionWithLabel(a);
188 5848 : if( !shortcut && !action ) {
189 0 : std::string str=" (hint! the actions with value in this ActionSet are: ";
190 0 : str+=as.getLabelList<ActionWithValue*>()+")";
191 0 : readact->error("cannot find action named " + a +str);
192 5848 : } else if( action && action->exists(c[i]) ) {
193 5760 : arg.push_back(action->copyOutput(c[i]));
194 88 : } else if( shortcut ) {
195 : unsigned narg=arg.size();
196 176 : shortcut->interpretDataLabel( a + "." + name, readact, arg );
197 88 : if( arg.size()==narg ) {
198 0 : readact->error("found no element in " + a + " with label " + name );
199 : }
200 : } else {
201 0 : std::string str=" (hint! the components in this actions are: ";
202 0 : str+=action->getComponentsList()+")";
203 0 : readact->error("action " + a + " has no component named " + name + str);
204 : }
205 : }
206 : } else { // if it doesn't contain a dot
207 19098 : if(c[i]=="*") {
208 : // Take all values from all actions
209 107 : std::vector<ActionWithValue*> all=as.select<ActionWithValue*>();
210 107 : if( all.empty() ) {
211 0 : readact->error("your input file is not telling plumed to calculate anything");
212 : }
213 1615 : for(unsigned j=0; j<all.size(); j++) {
214 1508 : plumed_assert(all[j]); // needed for following calls, see #1046
215 1508 : ActionWithVirtualAtom* av=all[j]->castToActionWithVirtualAtom();
216 1508 : if( av ) {
217 57 : continue;
218 : }
219 1451 : ActionForInterface* ap=all[j]->castToActionForInterface();
220 1451 : if( ap && all[j]->getName()!="ENERGY" ) {
221 834 : continue;
222 : }
223 1320 : for(int k=0; k<all[j]->getNumberOfComponents(); ++k) {
224 703 : arg.push_back(all[j]->copyOutput(k));
225 : }
226 : }
227 : } else {
228 18991 : ActionWithValue* action=as.selectWithLabel<ActionWithValue*>(c[i]);
229 18991 : if(!action) {
230 1 : std::string str=" (hint! the actions with value in this ActionSet are: ";
231 2 : str+=as.getLabelList<ActionWithValue*>()+")";
232 3 : readact->error("cannot find action named " + c[i] + str );
233 : }
234 18990 : if( !(action->exists(c[i])) ) {
235 0 : std::string str=" (hint! the components in this actions are: ";
236 0 : str+=action->getComponentsList()+")";
237 0 : readact->error("action " + c[i] + " has no component named " + c[i] +str);
238 : };
239 18991 : arg.push_back(action->copyOutput(c[i]));
240 : }
241 : }
242 : }
243 : }
244 68772 : for(unsigned i=0; i<arg.size(); ++i) {
245 53430 : if( !readact->keywords.checkArgumentType( arg[i]->getRank(), arg[i]->hasDerivatives() ) ) {
246 0 : readact->error("documentation for input type is not provided in " + readact->getName() );
247 : }
248 : }
249 15342 : }
250 :
251 0 : void ActionWithArguments::expandArgKeywordInPDB( const PDB& pdb ) {
252 0 : std::vector<std::string> arg_names = pdb.getArgumentNames();
253 0 : if( arg_names.size()>0 ) {
254 : std::vector<Value*> arg_vals;
255 0 : interpretArgumentList( arg_names, plumed.getActionSet(), this, arg_vals );
256 : }
257 0 : }
258 :
259 186049 : void ActionWithArguments::requestArguments(const std::vector<Value*> &arg) {
260 186049 : plumed_massert(!lockRequestArguments,"requested argument list can only be changed in the prepare() method");
261 186049 : arguments=arg;
262 186049 : clearDependencies();
263 : std::string fullname;
264 : std::string name;
265 1276076 : for(unsigned i=0; i<arguments.size(); i++) {
266 1090027 : fullname=arguments[i]->getName();
267 1090027 : if(fullname.find(".")!=std::string::npos) {
268 : std::size_t dot=fullname.find_first_of('.');
269 746172 : name=fullname.substr(0,dot);
270 : } else {
271 : name=fullname;
272 : }
273 1090027 : ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(name);
274 1090027 : plumed_massert(action,"cannot find action named (in requestArguments - this is weird)" + name);
275 1090027 : addDependency(action);
276 : }
277 186049 : ActionWithValue* av=dynamic_cast<ActionWithValue*>(this);
278 186049 : if(av) {
279 176400 : av->firststep=true;
280 : }
281 186049 : }
282 :
283 4 : void ActionWithArguments::requestExtraDependencies(const std::vector<Value*> &extra) {
284 4 : plumed_massert(!lockRequestArguments,"requested argument list can only be changed in the prepare() method");
285 : std::string fullname;
286 : std::string name;
287 9 : for(unsigned i=0; i<extra.size(); i++) {
288 5 : fullname=extra[i]->getName();
289 5 : if(fullname.find(".")!=std::string::npos) {
290 : std::size_t dot=fullname.find_first_of('.');
291 0 : name=fullname.substr(0,dot);
292 : } else {
293 : name=fullname;
294 : }
295 5 : ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(name);
296 5 : plumed_massert(action,"cannot find action named (in requestArguments - this is weird)" + name);
297 5 : addDependency(action);
298 : }
299 4 : }
300 :
301 10317 : ActionWithArguments::ActionWithArguments(const ActionOptions&ao):
302 : Action(ao),
303 10317 : lockRequestArguments(false) {
304 10317 : if( keywords.exists("ARG") ) {
305 : std::vector<Value*> arg;
306 19276 : parseArgumentList("ARG",arg);
307 :
308 9638 : if(!arg.empty()) {
309 9302 : log.printf(" with arguments : \n");
310 45605 : for(unsigned i=0; i<arg.size(); i++) {
311 72606 : if( arg[i]->hasDerivatives() && arg[i]->getRank()>0 ) {
312 1009 : log.printf(" function on grid with label %s \n",arg[i]->getName().c_str());
313 35294 : } else if( arg[i]->getRank()==2 ) {
314 2606 : log.printf(" matrix with label %s \n",arg[i]->getName().c_str());
315 32688 : } else if( arg[i]->getRank()==1 ) {
316 6131 : log.printf(" vector with label %s \n",arg[i]->getName().c_str());
317 26557 : } else if( arg[i]->getRank()==0 ) {
318 26557 : log.printf(" scalar with label %s \n",arg[i]->getName().c_str());
319 : } else {
320 0 : error("type of argument does not make sense");
321 : }
322 : }
323 : }
324 9638 : requestArguments(arg);
325 : }
326 10317 : }
327 :
328 58 : void ActionWithArguments::calculateNumericalDerivatives( ActionWithValue* a ) {
329 58 : if(!a) {
330 58 : a=castToActionWithValue();
331 58 : plumed_massert(a,"cannot compute numerical derivatives for an action without values");
332 : }
333 :
334 58 : const size_t nval=a->getNumberOfComponents();
335 : const size_t npar=arguments.size();
336 58 : std::vector<double> value (nval*npar);
337 161 : for(int i=0; i<npar; i++) {
338 103 : double arg0=arguments[i]->get();
339 103 : arguments[i]->set(arg0+std::sqrt(epsilon));
340 103 : a->calculate();
341 103 : arguments[i]->set(arg0);
342 1367 : for(int j=0; j<nval; j++) {
343 1264 : value[i*nval+j]=a->getOutputQuantity(j);
344 : }
345 : }
346 58 : a->calculate();
347 58 : a->clearDerivatives();
348 1192 : for(int j=0; j<nval; j++) {
349 1134 : Value* v=a->copyOutput(j);
350 1134 : if( v->hasDerivatives() )
351 670 : for(int i=0; i<npar; i++) {
352 400 : v->addDerivative(i,(value[i*nval+j]-a->getOutputQuantity(j))/std::sqrt(epsilon));
353 : }
354 : }
355 58 : }
356 :
357 261 : double ActionWithArguments::getProjection(unsigned i,unsigned j)const {
358 261 : plumed_massert(i<arguments.size()," making projections with an index which is too large");
359 261 : plumed_massert(j<arguments.size()," making projections with an index which is too large");
360 261 : const Value* v1=arguments[i];
361 261 : const Value* v2=arguments[j];
362 261 : return Value::projection(*v1,*v2);
363 : }
364 :
365 201793 : void ActionWithArguments::addForcesOnArguments( const unsigned& argstart, const std::vector<double>& forces, unsigned& ind, const std::string& c ) {
366 529559 : for(unsigned i=0; i<arguments.size(); ++i) {
367 327766 : if( i==0 && getName().find("EVALUATE_FUNCTION_FROM_GRID")!=std::string::npos ) {
368 3350 : continue ;
369 : }
370 324416 : if( !arguments[i]->ignoreStoredValue(c) || arguments[i]->getRank()==0 || (arguments[i]->getRank()>0 && arguments[i]->hasDerivatives()) ) {
371 302061 : unsigned nvals = arguments[i]->getNumberOfStoredValues();
372 33595371 : for(unsigned j=0; j<nvals; ++j) {
373 33293310 : arguments[i]->addForce( j, forces[ind], false );
374 33293310 : ind++;
375 : }
376 : }
377 : }
378 201793 : }
379 :
380 83 : void ActionWithArguments::setGradients( Value* myval, unsigned& start ) const {
381 83 : if( !myval->hasDeriv ) {
382 : return;
383 : }
384 83 : plumed_assert( myval->getRank()==0 );
385 :
386 : bool scalar=true;
387 249 : for(unsigned i=0; i<arguments.size(); ++i ) {
388 166 : if( arguments[i]->getRank()!=0 ) {
389 : scalar=false;
390 : break;
391 : }
392 : }
393 83 : if( !scalar ) {
394 : bool constant=true;
395 0 : for(unsigned i=0; i<arguments.size(); ++i ) {
396 0 : if( !arguments[i]->isConstant() ) {
397 : constant=false;
398 : break;
399 : } else {
400 0 : start += arguments[i]->getNumberOfValues();
401 : }
402 : }
403 0 : if( !constant ) {
404 0 : error("cannot set gradient as unable to handle non-constant actions that take vectors/matrices/grids in input");
405 : }
406 : }
407 : // Now pass the gradients
408 249 : for(unsigned i=0; i<arguments.size(); ++i ) {
409 166 : arguments[i]->passGradients( myval->getDerivative(i), myval->gradients );
410 : }
411 : }
412 :
413 18227 : bool ActionWithArguments::calculateConstantValues( const bool& haveatoms ) {
414 18227 : ActionWithValue* av = castToActionWithValue();
415 18227 : if( !av || arguments.size()==0 ) {
416 : return false;
417 : }
418 : bool constant = true, atoms=false;
419 14773 : for(unsigned i=0; i<arguments.size(); ++i) {
420 13661 : auto * ptr=arguments[i]->getPntrToAction();
421 13661 : plumed_assert(ptr); // needed for following calls, see #1046
422 13661 : ActionAtomistic* aa=ptr->castToActionAtomistic();
423 13661 : if( aa ) {
424 10777 : ActionWithVector* av=dynamic_cast<ActionWithVector*>( arguments[i]->getPntrToAction() );
425 10777 : if( !av || aa->getNumberOfAtoms()>0 ) {
426 : atoms=true;
427 : }
428 : }
429 13661 : if( !arguments[i]->isConstant() ) {
430 : constant=false;
431 : break;
432 : }
433 : }
434 13413 : if( constant ) {
435 : // Set everything constant first as we need to set the shape
436 2248 : for(unsigned i=0; i<av->getNumberOfComponents(); ++i) {
437 1136 : (av->copyOutput(i))->setConstant();
438 : }
439 1112 : if( !haveatoms ) {
440 1094 : log.printf(" values stored by this action are computed during startup and stay fixed during the simulation\n");
441 : }
442 1112 : if( atoms ) {
443 36 : return haveatoms;
444 : }
445 : }
446 : // Now do the calculation and store the values if we don't need anything from the atoms
447 13377 : if( constant && !haveatoms ) {
448 1076 : plumed_assert( !atoms );
449 1076 : activate();
450 1076 : calculate();
451 1076 : deactivate();
452 2176 : for(unsigned i=0; i<av->getNumberOfComponents(); ++i) {
453 1100 : unsigned nv = av->copyOutput(i)->getNumberOfValues();
454 1100 : log.printf(" %d values stored in component labelled %s are : ", nv, (av->copyOutput(i))->getName().c_str() );
455 3939 : for(unsigned j=0; j<nv; ++j) {
456 2839 : log.printf(" %f", (av->copyOutput(i))->get(j) );
457 : }
458 1100 : log.printf("\n");
459 : }
460 : }
461 : return constant;
462 : }
463 :
464 : }
|
