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