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 15487 : 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 15487 : }
49 :
50 11205 : void ActionWithArguments::parseArgumentList(const std::string&key,std::vector<Value*>&arg) {
51 22410 : 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 11205 : parseVector(key,c);
58 16690 : 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 11198 : interpretArgumentList(c,plumed.getActionSet(),this,arg);
66 11205 : }
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 15095 : void ActionWithArguments::interpretArgumentList(const std::vector<std::string>& c, const ActionSet& as, Action* readact, std::vector<Value*>&arg) {
83 38278 : for(unsigned i=0; i<c.size(); i++) {
84 : // is a regex? then just interpret it. The signal is ()
85 23186 : 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 22967 : std::string a=c[i].substr(0,dot);
121 22967 : std::string name=c[i].substr(dot+1);
122 22967 : if(c[i].find(".")!=std::string::npos) { // if it contains a dot:
123 6689 : 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(unsigned k=0; k<all[j]->getNumberOfComponents(); ++k) {
136 10 : arg.push_back(all[j]->copyOutput(k));
137 : }
138 : }
139 6670 : } else if ( name=="*") {
140 : unsigned carg=arg.size();
141 : // Take all the values from an action with a specific name
142 676 : ActionShortcut* shortcut=as.getShortcutActionWithLabel(a);
143 676 : if( shortcut ) {
144 872 : shortcut->interpretDataLabel( a + "." + name, readact, arg );
145 : }
146 676 : if( arg.size()==carg ) {
147 : // Take all the values from an action with a specific name
148 442 : ActionWithValue* action=as.selectWithLabel<ActionWithValue*>(a);
149 442 : 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 442 : 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 7051 : for(unsigned k=0; k<action->getNumberOfComponents(); ++k) {
158 6609 : arg.push_back(action->copyOutput(k));
159 : }
160 : }
161 5994 : } 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 167 : for(unsigned j=0; j<shortcuts.size(); ++j) {
170 300 : shortcuts[j]->interpretDataLabel( shortcuts[j]->getShortcutLabel() + "." + name, readact, arg );
171 : }
172 : unsigned nval=0;
173 333 : for(unsigned j=0; j<all.size(); j++) {
174 : std::string flab;
175 632 : flab=all[j]->getLabel() + "." + name;
176 316 : 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 5977 : ActionWithValue* action=as.selectWithLabel<ActionWithValue*>(a);
187 5977 : ActionShortcut* shortcut=as.getShortcutActionWithLabel(a);
188 5977 : 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 5977 : } else if( action && action->exists(c[i]) ) {
193 5890 : arg.push_back(action->copyOutput(c[i]));
194 87 : } else if( shortcut ) {
195 : unsigned narg=arg.size();
196 174 : shortcut->interpretDataLabel( a + "." + name, readact, arg );
197 87 : 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 16296 : 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(unsigned k=0; k<all[j]->getNumberOfComponents(); ++k) {
224 703 : arg.push_back(all[j]->copyOutput(k));
225 : }
226 : }
227 : } else {
228 16189 : ActionWithValue* action=as.selectWithLabel<ActionWithValue*>(c[i]);
229 16189 : 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 16188 : 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 16189 : arg.push_back(action->copyOutput(c[i]));
240 : }
241 : }
242 : }
243 : }
244 66517 : for(unsigned i=0; i<arg.size(); ++i) {
245 51425 : 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 15092 : if( readact->keywords.exists("MASKED_INPUT_ALLOWED") || readact->keywords.exists("IS_SHORTCUT") || readact->keywords.exists("MASK") ) {
250 10678 : return;
251 : }
252 39308 : for(unsigned i=0; i<arg.size(); ++i) {
253 34894 : if( arg[i]->getRank()==0 ) {
254 32819 : continue;
255 : }
256 2075 : ActionWithVector* av=dynamic_cast<ActionWithVector*>( arg[i]->getPntrToAction() );
257 2075 : if( av && av->getNumberOfMasks()>=0 ) {
258 0 : readact->error("cannot use argument " + arg[i]->getName() + " in input as not all elements are computed");
259 : }
260 : }
261 : }
262 :
263 0 : void ActionWithArguments::expandArgKeywordInPDB( const PDB& pdb ) {
264 0 : std::vector<std::string> arg_names = pdb.getArgumentNames();
265 0 : if( arg_names.size()>0 ) {
266 : std::vector<Value*> arg_vals;
267 0 : interpretArgumentList( arg_names, plumed.getActionSet(), this, arg_vals );
268 : }
269 0 : }
270 :
271 185051 : void ActionWithArguments::requestArguments(const std::vector<Value*> &arg) {
272 185051 : plumed_massert(!lockRequestArguments,"requested argument list can only be changed in the prepare() method");
273 185051 : arguments=arg;
274 185051 : clearDependencies();
275 : std::string fullname;
276 : std::string argName;
277 1274304 : for(unsigned i=0; i<arguments.size(); i++) {
278 1089253 : fullname=arguments[i]->getName();
279 1089253 : if(fullname.find(".")!=std::string::npos) {
280 : std::size_t dot=fullname.find_first_of('.');
281 746726 : argName=fullname.substr(0,dot);
282 : } else {
283 : argName=fullname;
284 : }
285 1089253 : ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(argName);
286 1089253 : plumed_massert(action,"cannot find action named (in requestArguments - this is weird)" + argName);
287 1089253 : addDependency(action);
288 : }
289 185051 : ActionWithValue* av=dynamic_cast<ActionWithValue*>(this);
290 185051 : if(av) {
291 176655 : av->firststep=true;
292 : }
293 185051 : }
294 :
295 4 : void ActionWithArguments::requestExtraDependencies(const std::vector<Value*> &extra) {
296 4 : plumed_massert(!lockRequestArguments,"requested argument list can only be changed in the prepare() method");
297 : std::string fullname;
298 : std::string argName;
299 9 : for(unsigned i=0; i<extra.size(); i++) {
300 5 : fullname=extra[i]->getName();
301 5 : if(fullname.find(".")!=std::string::npos) {
302 : std::size_t dot=fullname.find_first_of('.');
303 0 : argName=fullname.substr(0,dot);
304 : } else {
305 : argName=fullname;
306 : }
307 5 : ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(argName);
308 5 : plumed_massert(action,"cannot find action named (in requestArguments - this is weird)" + argName);
309 5 : addDependency(action);
310 : }
311 4 : }
312 :
313 9050 : ActionWithArguments::ActionWithArguments(const ActionOptions&ao):
314 : Action(ao),
315 9050 : lockRequestArguments(false) {
316 9050 : if( keywords.exists("ARG") ) {
317 : std::vector<Value*> arg;
318 16770 : parseArgumentList("ARG",arg);
319 :
320 8385 : if(!arg.empty()) {
321 8047 : log.printf(" with arguments : \n");
322 42777 : for(unsigned i=0; i<arg.size(); i++) {
323 34730 : if( arg[i]->hasDerivatives() && arg[i]->getRank()>0 ) {
324 345 : log.printf(" function on grid with label %s \n",arg[i]->getName().c_str());
325 34385 : } else if( arg[i]->getRank()==2 ) {
326 2475 : log.printf(" matrix with label %s \n",arg[i]->getName().c_str());
327 31910 : } else if( arg[i]->getRank()==1 ) {
328 6093 : log.printf(" vector with label %s \n",arg[i]->getName().c_str());
329 25817 : } else if( arg[i]->getRank()==0 ) {
330 25817 : log.printf(" scalar with label %s \n",arg[i]->getName().c_str());
331 : } else {
332 0 : error("type of argument does not make sense");
333 : }
334 : }
335 : }
336 8385 : requestArguments(arg);
337 : }
338 9050 : }
339 :
340 58 : void ActionWithArguments::calculateNumericalDerivatives( ActionWithValue* a ) {
341 58 : if(!a) {
342 58 : a=castToActionWithValue();
343 58 : plumed_massert(a,"cannot compute numerical derivatives for an action without values");
344 : }
345 :
346 : const size_t nval=a->getNumberOfComponents();
347 : const size_t npar=arguments.size();
348 58 : std::vector<double> value (nval*npar);
349 161 : for(unsigned i=0; i<npar; i++) {
350 103 : double arg0=arguments[i]->get();
351 103 : arguments[i]->set(arg0+std::sqrt(epsilon));
352 103 : a->calculate();
353 103 : arguments[i]->set(arg0);
354 1367 : for(unsigned j=0; j<nval; j++) {
355 1264 : value[i*nval+j]=a->getOutputQuantity(j);
356 : }
357 : }
358 58 : a->calculate();
359 58 : a->clearDerivatives();
360 1192 : for(unsigned j=0; j<nval; j++) {
361 1134 : Value* v=a->copyOutput(j);
362 1134 : if( v->hasDerivatives() )
363 670 : for(unsigned i=0; i<npar; i++) {
364 400 : v->addDerivative(i,(value[i*nval+j]-a->getOutputQuantity(j))/std::sqrt(epsilon));
365 : }
366 : }
367 58 : }
368 :
369 261 : double ActionWithArguments::getProjection(unsigned i,unsigned j)const {
370 261 : plumed_massert(i<arguments.size()," making projections with an index which is too large");
371 261 : plumed_massert(j<arguments.size()," making projections with an index which is too large");
372 261 : const Value* v1=arguments[i];
373 261 : const Value* v2=arguments[j];
374 261 : return Value::projection(*v1,*v2);
375 : }
376 :
377 257403 : void ActionWithArguments::addForcesOnArguments( const unsigned& argstart, const std::vector<double>& forces, unsigned& ind ) {
378 257403 : unsigned nargs=arguments.size();
379 257403 : const ActionWithVector* av=dynamic_cast<const ActionWithVector*>( this );
380 257403 : if( av && av->getNumberOfMasks()>0 ) {
381 796 : nargs=nargs-av->getNumberOfMasks();
382 : }
383 595431 : for(unsigned i=0; i<nargs; ++i) {
384 338028 : ind += arguments[i]->addForces(View(&forces[ind],forces.size()-ind));
385 : }
386 257403 : }
387 :
388 83 : void ActionWithArguments::setGradients( Value* myval, unsigned& start ) const {
389 83 : if( !myval->hasDeriv ) {
390 : return;
391 : }
392 83 : plumed_assert( myval->getRank()==0 );
393 :
394 : bool scalar=true;
395 249 : for(unsigned i=0; i<arguments.size(); ++i ) {
396 166 : if( arguments[i]->getRank()!=0 ) {
397 : scalar=false;
398 : break;
399 : }
400 : }
401 83 : if( !scalar ) {
402 : bool constant=true;
403 0 : for(unsigned i=0; i<arguments.size(); ++i ) {
404 0 : if( !arguments[i]->isConstant() ) {
405 : constant=false;
406 : break;
407 : } else {
408 0 : start += arguments[i]->getNumberOfValues();
409 : }
410 : }
411 0 : if( !constant ) {
412 0 : error("cannot set gradient as unable to handle non-constant actions that take vectors/matrices/grids in input");
413 : }
414 : }
415 : // Now pass the gradients
416 249 : for(unsigned i=0; i<arguments.size(); ++i ) {
417 166 : arguments[i]->passGradients( myval->getDerivative(i), myval->gradients );
418 : }
419 : }
420 :
421 16537 : bool ActionWithArguments::calculateConstantValues( const bool& haveatoms ) {
422 16537 : ActionWithValue* awval = castToActionWithValue();
423 16537 : if( !awval || arguments.size()==0 ) {
424 : return false;
425 : }
426 : bool constant = true, atoms=false;
427 12522 : for(unsigned i=0; i<arguments.size(); ++i) {
428 12185 : auto * ptr=arguments[i]->getPntrToAction();
429 12185 : plumed_assert(ptr); // needed for following calls, see #1046
430 12185 : ActionAtomistic* aa=ptr->castToActionAtomistic();
431 12185 : if( aa ) {
432 8829 : ActionWithVector* awvec=dynamic_cast<ActionWithVector*>( arguments[i]->getPntrToAction() );
433 8829 : if( !awvec || aa->getNumberOfAtoms()>0 ) {
434 : atoms=true;
435 : }
436 : }
437 12185 : if( !arguments[i]->isConstant() ) {
438 : constant=false;
439 : break;
440 : }
441 : }
442 11713 : if( constant ) {
443 : // Set everything constant first as we need to set the shape
444 698 : for(unsigned i=0; i<awval->getNumberOfComponents(); ++i) {
445 361 : (awval->copyOutput(i))->setConstant();
446 : }
447 337 : if( !haveatoms ) {
448 319 : log.printf(" values stored by this action are computed during startup and stay fixed during the simulation\n");
449 : }
450 337 : if( atoms ) {
451 36 : return haveatoms;
452 : }
453 : }
454 : // Now do the calculation and store the values if we don't need anything from the atoms
455 11677 : if( constant && !haveatoms ) {
456 301 : plumed_assert( !atoms );
457 301 : activate();
458 301 : calculate();
459 301 : deactivate();
460 626 : for(unsigned i=0; i<awval->getNumberOfComponents(); ++i) {
461 325 : unsigned nv = awval->copyOutput(i)->getNumberOfValues();
462 325 : log.printf(" %d values stored in component labelled %s are : ", nv, (awval->copyOutput(i))->getName().c_str() );
463 6827 : for(unsigned j=0; j<nv; ++j) {
464 6502 : log.printf(" %f", (awval->copyOutput(i))->get(j) );
465 : }
466 325 : log.printf("\n");
467 : }
468 : }
469 : return constant;
470 : }
471 :
472 : }
|
