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 17846 : 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 17846 : }
49 :
50 9869 : void ActionWithArguments::parseArgumentList(const std::string&key,std::vector<Value*>&arg) {
51 19762 : if( keywords.getArgumentType(key).length()==0 ) warning("keyword " + key + " for reading arguments should is registered using Keyword::add rather than Keyword::addInputKeyword. The keyword will thus not appear in the correct place in the manual");
52 9869 : std::string def; std::vector<std::string> c; arg.clear(); parseVector(key,c);
53 10720 : if( c.size()==0 && (keywords.style(key,"compulsory") || keywords.style(key,"hidden")) ) {
54 7 : if( keywords.getDefaultValue(key,def) ) c.push_back( def );
55 : else return;
56 : }
57 9862 : interpretArgumentList(c,plumed.getActionSet(),this,arg);
58 9869 : }
59 :
60 104 : bool ActionWithArguments::parseArgumentList(const std::string&key,int i,std::vector<Value*>&arg) {
61 208 : if( keywords.getArgumentType(key).length()==0 ) warning("keyword " + key + " for reading argument should is registered using Keyword::add rather than Keyword::addInputKeyword. The keyword will thus not appear in the correct place in the manual");
62 : std::vector<std::string> c;
63 : arg.clear();
64 104 : if(parseNumberedVector(key,i,c)) {
65 44 : interpretArgumentList(c,plumed.getActionSet(),this,arg);
66 : return true;
67 : } else return false;
68 104 : }
69 :
70 15213 : void ActionWithArguments::interpretArgumentList(const std::vector<std::string>& c, const ActionSet& as, Action* readact, std::vector<Value*>&arg) {
71 40875 : for(unsigned i=0; i<c.size(); i++) {
72 : // is a regex? then just interpret it. The signal is ()
73 25665 : if(!c[i].compare(0,1,"(")) {
74 219 : unsigned l=c[i].length();
75 219 : if(!c[i].compare(l-1,1,")")) {
76 : // start regex parsing
77 : bool found_something=false;
78 : // take the string enclosed in quotes and put in round brackets
79 218 : std::string myregex=c[i];
80 218 : std::vector<ActionWithValue*> all=as.select<ActionWithValue*>();
81 219 : if( all.empty() ) readact->error("your input file is not telling plumed to calculate anything");
82 :
83 : try {
84 218 : std::regex txt_regex(myregex,std::regex::extended);
85 217 : plumed_massert(txt_regex.mark_count()==1,"I can parse with only one subexpression");
86 24239 : for(unsigned j=0; j<all.size(); j++) {
87 24022 : std::vector<std::string> ss=all[j]->getComponentsVector();
88 349707 : for(unsigned k=0; k<ss.size(); ++k) {
89 325685 : if(std::regex_match(ss[k],txt_regex)) {
90 21601 : arg.push_back(all[j]->copyOutput(ss[k]));
91 : found_something=true;
92 : }
93 : }
94 24022 : }
95 218 : } catch(std::regex_error & e) {
96 3 : plumed_error()<<"Error parsing regular expression: "<<e.what();
97 1 : }
98 217 : if(!found_something) plumed_error()<<"There isn't any action matching your regex " << myregex;
99 : } else {
100 2 : plumed_merror("did you want to use regexp to input arguments? enclose it between two round braces (...) with no spaces!");
101 : }
102 : } else {
103 : std::size_t dot=c[i].find_first_of('.');
104 25446 : std::string a=c[i].substr(0,dot);
105 25446 : std::string name=c[i].substr(dot+1);
106 25446 : if(c[i].find(".")!=std::string::npos) { // if it contains a dot:
107 6380 : if(a=="*" && name=="*") {
108 : // Take all values from all actions
109 1 : std::vector<ActionWithValue*> all=as.select<ActionWithValue*>();
110 1 : if( all.empty() ) readact->error("your input file is not telling plumed to calculate anything");
111 17 : for(unsigned j=0; j<all.size(); j++) {
112 16 : plumed_assert(all[j]); // needed for following calls, see #1046
113 16 : ActionForInterface* ap=all[j]->castToActionForInterface(); if( ap ) continue;
114 18 : for(int k=0; k<all[j]->getNumberOfComponents(); ++k) arg.push_back(all[j]->copyOutput(k));
115 : }
116 6365 : } else if ( name=="*") {
117 : unsigned carg=arg.size();
118 : // Take all the values from an action with a specific name
119 606 : ActionShortcut* shortcut=as.getShortcutActionWithLabel(a);
120 991 : if( shortcut ) shortcut->interpretDataLabel( a + "." + name, readact, arg );
121 606 : if( arg.size()==carg ) {
122 : // Take all the values from an action with a specific name
123 375 : ActionWithValue* action=as.selectWithLabel<ActionWithValue*>(a);
124 375 : if(!action) {
125 0 : std::string str=" (hint! the actions with value in this ActionSet are: ";
126 0 : str+=as.getLabelList<ActionWithValue*>()+")";
127 0 : readact->error("cannot find action named " + a + str);
128 : }
129 375 : if( action->getNumberOfComponents()==0 ) readact->error("found " + a +".* indicating use all components calculated by action with label " + a + " but this action has no components");
130 6171 : for(int k=0; k<action->getNumberOfComponents(); ++k) arg.push_back(action->copyOutput(k));
131 : }
132 5759 : } else if ( a=="*" ) {
133 13 : std::vector<ActionShortcut*> shortcuts=as.select<ActionShortcut*>();
134 : // Take components from all actions with a specific name
135 13 : std::vector<ActionWithValue*> all=as.select<ActionWithValue*>();
136 13 : if( all.empty() ) readact->error("your input file is not telling plumed to calculate anything");
137 : unsigned carg=arg.size();
138 157 : for(unsigned j=0; j<shortcuts.size(); ++j) {
139 288 : shortcuts[j]->interpretDataLabel( shortcuts[j]->getShortcutLabel() + "." + name, readact, arg );
140 : }
141 : unsigned nval=0;
142 276 : for(unsigned j=0; j<all.size(); j++) {
143 526 : std::string flab; flab=all[j]->getLabel() + "." + name;
144 263 : if( all[j]->exists(flab) ) { arg.push_back(all[j]->copyOutput(flab)); nval++; }
145 : }
146 13 : if(nval==0 && arg.size()==carg) readact->error("found no actions with a component called " + name );
147 : } else {
148 : // Take values with a specific name
149 5746 : ActionWithValue* action=as.selectWithLabel<ActionWithValue*>(a);
150 5746 : ActionShortcut* shortcut=as.getShortcutActionWithLabel(a);
151 5746 : if( !shortcut && !action ) {
152 0 : std::string str=" (hint! the actions with value in this ActionSet are: ";
153 0 : str+=as.getLabelList<ActionWithValue*>()+")";
154 0 : readact->error("cannot find action named " + a +str);
155 5746 : } else if( action && action->exists(c[i]) ) {
156 5662 : arg.push_back(action->copyOutput(c[i]));
157 84 : } else if( shortcut ) {
158 168 : unsigned narg=arg.size(); shortcut->interpretDataLabel( a + "." + name, readact, arg );
159 84 : if( arg.size()==narg ) readact->error("found no element in " + a + " with label " + name );
160 : } else {
161 0 : std::string str=" (hint! the components in this actions are: ";
162 0 : str+=action->getComponentsList()+")";
163 0 : readact->error("action " + a + " has no component named " + name + str);
164 : }
165 : }
166 : } else { // if it doesn't contain a dot
167 19080 : if(c[i]=="*") {
168 : // Take all values from all actions
169 107 : std::vector<ActionWithValue*> all=as.select<ActionWithValue*>();
170 107 : if( all.empty() ) readact->error("your input file is not telling plumed to calculate anything");
171 1615 : for(unsigned j=0; j<all.size(); j++) {
172 1508 : plumed_assert(all[j]); // needed for following calls, see #1046
173 1508 : ActionWithVirtualAtom* av=all[j]->castToActionWithVirtualAtom(); if( av ) continue;
174 1451 : ActionForInterface* ap=all[j]->castToActionForInterface(); if( ap && all[j]->getName()!="ENERGY" ) continue;
175 1320 : for(int k=0; k<all[j]->getNumberOfComponents(); ++k) arg.push_back(all[j]->copyOutput(k));
176 : }
177 : } else {
178 18973 : ActionWithValue* action=as.selectWithLabel<ActionWithValue*>(c[i]);
179 18973 : if(!action) {
180 1 : std::string str=" (hint! the actions with value in this ActionSet are: ";
181 2 : str+=as.getLabelList<ActionWithValue*>()+")";
182 3 : readact->error("cannot find action named " + c[i] + str );
183 : }
184 18972 : if( !(action->exists(c[i])) ) {
185 0 : std::string str=" (hint! the components in this actions are: ";
186 0 : str+=action->getComponentsList()+")";
187 0 : readact->error("action " + c[i] + " has no component named " + c[i] +str);
188 : };
189 18973 : arg.push_back(action->copyOutput(c[i]));
190 : }
191 : }
192 : }
193 : }
194 68378 : for(unsigned i=0; i<arg.size(); ++i) {
195 53168 : if( !readact->keywords.checkArgumentType( arg[i]->getRank(), arg[i]->hasDerivatives() ) ) readact->warning("documentation for input type is not provided in " + readact->getName() );
196 : }
197 15210 : }
198 :
199 0 : void ActionWithArguments::expandArgKeywordInPDB( const PDB& pdb ) {
200 0 : std::vector<std::string> arg_names = pdb.getArgumentNames();
201 0 : if( arg_names.size()>0 ) {
202 : std::vector<Value*> arg_vals;
203 0 : interpretArgumentList( arg_names, plumed.getActionSet(), this, arg_vals );
204 : }
205 0 : }
206 :
207 185938 : void ActionWithArguments::requestArguments(const std::vector<Value*> &arg) {
208 185938 : plumed_massert(!lockRequestArguments,"requested argument list can only be changed in the prepare() method");
209 185938 : arguments=arg;
210 185938 : clearDependencies();
211 : std::string fullname;
212 : std::string name;
213 1275721 : for(unsigned i=0; i<arguments.size(); i++) {
214 1089783 : fullname=arguments[i]->getName();
215 1089783 : if(fullname.find(".")!=std::string::npos) {
216 : std::size_t dot=fullname.find_first_of('.');
217 745754 : name=fullname.substr(0,dot);
218 : } else {
219 : name=fullname;
220 : }
221 1089783 : ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(name);
222 1089783 : plumed_massert(action,"cannot find action named (in requestArguments - this is weird)" + name);
223 1089783 : addDependency(action);
224 : }
225 185938 : ActionWithValue* av=dynamic_cast<ActionWithValue*>(this);
226 185938 : if(av) av->firststep=true;
227 185938 : }
228 :
229 4 : void ActionWithArguments::requestExtraDependencies(const std::vector<Value*> &extra) {
230 4 : plumed_massert(!lockRequestArguments,"requested argument list can only be changed in the prepare() method");
231 : std::string fullname;
232 : std::string name;
233 9 : for(unsigned i=0; i<extra.size(); i++) {
234 5 : fullname=extra[i]->getName();
235 5 : if(fullname.find(".")!=std::string::npos) {
236 : std::size_t dot=fullname.find_first_of('.');
237 0 : name=fullname.substr(0,dot);
238 : } else {
239 : name=fullname;
240 : }
241 5 : ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(name);
242 5 : plumed_massert(action,"cannot find action named (in requestArguments - this is weird)" + name);
243 5 : addDependency(action);
244 : }
245 4 : }
246 :
247 10200 : ActionWithArguments::ActionWithArguments(const ActionOptions&ao):
248 : Action(ao),
249 10200 : lockRequestArguments(false)
250 : {
251 20400 : if( keywords.exists("ARG") ) {
252 : std::vector<Value*> arg;
253 19054 : parseArgumentList("ARG",arg);
254 :
255 9527 : if(!arg.empty()) {
256 9225 : log.printf(" with arguments : \n");
257 45284 : for(unsigned i=0; i<arg.size(); i++) {
258 36059 : if( arg[i]->hasDerivatives() && arg[i]->getRank()>0 ) log.printf(" function on grid with label %s \n",arg[i]->getName().c_str());
259 35050 : else if( arg[i]->getRank()==2 ) log.printf(" matrix with label %s \n",arg[i]->getName().c_str());
260 32444 : else if( arg[i]->getRank()==1 ) log.printf(" vector with label %s \n",arg[i]->getName().c_str());
261 26331 : else if( arg[i]->getRank()==0 ) log.printf(" scalar with label %s \n",arg[i]->getName().c_str());
262 0 : else error("type of argument does not make sense");
263 : }
264 : }
265 9527 : requestArguments(arg);
266 : }
267 10200 : }
268 :
269 58 : void ActionWithArguments::calculateNumericalDerivatives( ActionWithValue* a ) {
270 58 : if(!a) {
271 58 : a=castToActionWithValue();
272 58 : plumed_massert(a,"cannot compute numerical derivatives for an action without values");
273 : }
274 :
275 58 : const size_t nval=a->getNumberOfComponents();
276 : const size_t npar=arguments.size();
277 58 : std::vector<double> value (nval*npar);
278 161 : for(int i=0; i<npar; i++) {
279 103 : double arg0=arguments[i]->get();
280 103 : arguments[i]->set(arg0+std::sqrt(epsilon));
281 103 : a->calculate();
282 103 : arguments[i]->set(arg0);
283 1367 : for(int j=0; j<nval; j++) {
284 1264 : value[i*nval+j]=a->getOutputQuantity(j);
285 : }
286 : }
287 58 : a->calculate();
288 58 : a->clearDerivatives();
289 1192 : for(int j=0; j<nval; j++) {
290 1134 : Value* v=a->copyOutput(j);
291 1804 : if( v->hasDerivatives() ) for(int i=0; i<npar; i++) v->addDerivative(i,(value[i*nval+j]-a->getOutputQuantity(j))/std::sqrt(epsilon));
292 : }
293 58 : }
294 :
295 261 : double ActionWithArguments::getProjection(unsigned i,unsigned j)const {
296 261 : plumed_massert(i<arguments.size()," making projections with an index which is too large");
297 261 : plumed_massert(j<arguments.size()," making projections with an index which is too large");
298 261 : const Value* v1=arguments[i];
299 261 : const Value* v2=arguments[j];
300 261 : return Value::projection(*v1,*v2);
301 : }
302 :
303 201619 : void ActionWithArguments::addForcesOnArguments( const unsigned& argstart, const std::vector<double>& forces, unsigned& ind, const std::string& c ) {
304 529373 : for(unsigned i=0; i<arguments.size(); ++i) {
305 327754 : if( i==0 && getName().find("EVALUATE_FUNCTION_FROM_GRID")!=std::string::npos ) continue ;
306 324404 : if( !arguments[i]->ignoreStoredValue(c) || arguments[i]->getRank()==0 || (arguments[i]->getRank()>0 && arguments[i]->hasDerivatives()) ) {
307 302049 : unsigned nvals = arguments[i]->getNumberOfStoredValues();
308 33595347 : for(unsigned j=0; j<nvals; ++j) { arguments[i]->addForce( j, forces[ind], false ); ind++; }
309 : }
310 : }
311 201619 : }
312 :
313 83 : void ActionWithArguments::setGradients( Value* myval, unsigned& start ) const {
314 83 : if( !myval->hasDeriv ) return; plumed_assert( myval->getRank()==0 );
315 :
316 : bool scalar=true;
317 249 : for(unsigned i=0; i<arguments.size(); ++i ) {
318 166 : if( arguments[i]->getRank()!=0 ) { scalar=false; break; }
319 : }
320 83 : if( !scalar ) {
321 : bool constant=true;
322 0 : for(unsigned i=0; i<arguments.size(); ++i ) {
323 0 : if( !arguments[i]->isConstant() ) { constant=false; break; }
324 0 : else start += arguments[i]->getNumberOfValues();
325 : }
326 0 : if( !constant ) error("cannot set gradient as unable to handle non-constant actions that take vectors/matrices/grids in input");
327 : }
328 : // Now pass the gradients
329 249 : for(unsigned i=0; i<arguments.size(); ++i ) arguments[i]->passGradients( myval->getDerivative(i), myval->gradients );
330 : }
331 :
332 18111 : bool ActionWithArguments::calculateConstantValues( const bool& haveatoms ) {
333 18111 : ActionWithValue* av = castToActionWithValue();
334 18111 : if( !av || arguments.size()==0 ) return false;
335 : bool constant = true, atoms=false;
336 14786 : for(unsigned i=0; i<arguments.size(); ++i) {
337 13674 : auto * ptr=arguments[i]->getPntrToAction();
338 13674 : plumed_assert(ptr); // needed for following calls, see #1046
339 13674 : ActionAtomistic* aa=ptr->castToActionAtomistic();
340 13674 : if( aa ) {
341 10770 : ActionWithVector* av=dynamic_cast<ActionWithVector*>( arguments[i]->getPntrToAction() );
342 10770 : if( !av || aa->getNumberOfAtoms()>0 ) atoms=true;
343 : }
344 13674 : if( !arguments[i]->isConstant() ) { constant=false; break; }
345 : }
346 13426 : if( constant ) {
347 : // Set everything constant first as we need to set the shape
348 2248 : for(unsigned i=0; i<av->getNumberOfComponents(); ++i) (av->copyOutput(i))->setConstant();
349 1112 : if( !haveatoms ) log.printf(" values stored by this action are computed during startup and stay fixed during the simulation\n");
350 1112 : if( atoms ) return haveatoms;
351 : }
352 : // Now do the calculation and store the values if we don't need anything from the atoms
353 13390 : if( constant && !haveatoms ) {
354 1076 : plumed_assert( !atoms ); activate(); calculate(); deactivate();
355 2176 : for(unsigned i=0; i<av->getNumberOfComponents(); ++i) {
356 1100 : unsigned nv = av->copyOutput(i)->getNumberOfValues();
357 1100 : log.printf(" %d values stored in component labelled %s are : ", nv, (av->copyOutput(i))->getName().c_str() );
358 3939 : for(unsigned j=0; j<nv; ++j) log.printf(" %f", (av->copyOutput(i))->get(j) );
359 1100 : log.printf("\n");
360 : }
361 : }
362 : return constant;
363 : }
364 :
365 : }
|
