Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 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 "tools/PDB.h"
25 : #include "PlumedMain.h"
26 : #include "ActionSet.h"
27 : #include <iostream>
28 : #ifdef __PLUMED_HAS_CREGEX
29 : #include <cstring>
30 : #include "regex.h"
31 : #endif
32 :
33 : using namespace std;
34 : namespace PLMD {
35 :
36 991 : void ActionWithArguments::registerKeywords(Keywords& keys) {
37 : 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 "
38 : "are referenced using the label of the action. If the label appears on its own then it is assumed that the Action calculates "
39 : "a single scalar value. The value of this scalar is thus used as the input to this new action. If * or *.* appears the "
40 : "scalars calculated by all the proceding actions in the input file are taken. Some actions have multi-component outputs and "
41 : "each component of the output has a specific label. For example a \\ref DISTANCE action labelled dist may have three componets "
42 : "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.*."
43 : "More information on the referencing of Actions can be found in the section of the manual on the PLUMED \\ref Syntax. "
44 : "Scalar values can also be "
45 : "referenced using POSIX regular expressions as detailed in the section on \\ref Regex. To use this feature you you must compile "
46 991 : "PLUMED with the appropriate flag.");
47 991 : }
48 :
49 948 : void ActionWithArguments::parseArgumentList(const std::string&key,std::vector<Value*>&arg) {
50 948 : vector<string> c; arg.clear(); parseVector(key,c);
51 948 : if( c.size()==0 && (keywords.style(key,"compulsory") || keywords.style(key,"hidden")) ) {
52 2 : std::string def; if( keywords.getDefaultValue(key,def) ) c.push_back( def );
53 : }
54 948 : interpretArgumentList(c,arg);
55 948 : }
56 :
57 6 : bool ActionWithArguments::parseArgumentList(const std::string&key,int i,std::vector<Value*>&arg) {
58 6 : vector<string> c;
59 6 : arg.clear();
60 6 : if(parseNumberedVector(key,i,c)) {
61 6 : interpretArgumentList(c,arg);
62 6 : return true;
63 0 : } else return false;
64 : }
65 :
66 965 : void ActionWithArguments::interpretArgumentList(const std::vector<std::string>& c, std::vector<Value*>&arg) {
67 2970 : for(unsigned i=0; i<c.size(); i++) {
68 : // is a regex? then just interpret it. The signal is ()
69 2005 : if(!c[i].compare(0,1,"(")) {
70 136 : unsigned l=c[i].length();
71 136 : if(!c[i].compare(l-1,1,")")) {
72 : // start regex parsing
73 : #ifdef __PLUMED_HAS_CREGEX
74 : // take the string enclosed in quotes and put in round brackets
75 136 : std::string myregex=c[i];
76 136 : log.printf(" Evaluating regexp for this action: %s \n",myregex.c_str());
77 : int errcode;
78 136 : regex_t *preg = (regex_t*)malloc(sizeof(regex_t)); // pointer to the regular expression
79 : regmatch_t *pmatch;
80 136 : if ((errcode=regcomp(preg, myregex.c_str(),REG_EXTENDED|REG_NEWLINE))) { // compile the regular expression
81 : char* errbuf;
82 : size_t errbuf_size;
83 : // one can check the errors asking to regerror
84 0 : errbuf_size = regerror(errcode, preg, NULL, 0);
85 0 : if (!(errbuf=(char*)malloc(errbuf_size))) {
86 0 : plumed_merror("cannot allocate the buffer for error detection in regexp!");
87 : };
88 0 : regerror(errcode, preg, errbuf, errbuf_size);
89 0 : error(errbuf);
90 : }
91 136 : plumed_massert(preg->re_nsub==1,"I can parse with only one subexpression");
92 136 : pmatch = (regmatch_t*)malloc(sizeof(regmatch_t)*preg->re_nsub);
93 : // select all the actions that have a value
94 272 : std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>();
95 136 : if( all.empty() ) error("your input file is not telling plumed to calculate anything");
96 621 : for(unsigned j=0; j<all.size(); j++) {
97 485 : std::vector<std::string> ss=all[j]->getComponentsVector();
98 173218 : for(unsigned k=0; k<ss.size(); ++k) {
99 172733 : unsigned ll=strlen(ss[k].c_str())+1;
100 : char*str;
101 172733 : str=new char [ll];
102 172733 : strcpy(str,ss[k].c_str());
103 172733 : char *ppstr=str;
104 172733 : if(!regexec(preg, ppstr, preg->re_nsub, pmatch, 0)) {
105 21899 : log.printf(" Something matched with \"%s\" : ",ss[k].c_str());
106 21899 : do {
107 21899 : if (pmatch[0].rm_so != -1) { /* The regex is matching part of a string */
108 : char *submatch;
109 21899 : size_t matchlen = pmatch[0].rm_eo - pmatch[0].rm_so;
110 21899 : submatch = (char*)malloc(matchlen+1);
111 21899 : strncpy(submatch, ppstr+pmatch[0].rm_so, matchlen+1);
112 21899 : submatch[matchlen]='\0';
113 21899 : log.printf(" subpattern %s\n", submatch);
114 : // this is the match: try to see if it is a valid action
115 21899 : std::string putativeVal(submatch);
116 21899 : if( all[j]->exists(putativeVal) ) {
117 15398 : arg.push_back(all[j]->copyOutput(putativeVal));
118 15398 : log.printf(" Action %s added! \n",putativeVal.c_str());
119 : }
120 21899 : free(submatch);
121 : };
122 21899 : ppstr += pmatch[0].rm_eo; /* Restart from last match */
123 21899 : } while(!regexec(preg,ppstr,preg->re_nsub,pmatch,0));
124 : }
125 172733 : delete [] str;
126 : }
127 485 : };
128 136 : regfree(preg);
129 136 : free(preg);
130 272 : free(pmatch);
131 : #else
132 : plumed_merror("Regexp support not compiled!");
133 : #endif
134 : } else {
135 0 : plumed_merror("did you want to use regexp to input arguments? enclose it between two round braces (...) with no spaces!");
136 : }
137 : } else {
138 1869 : std::size_t dot=c[i].find_first_of('.');
139 1869 : string a=c[i].substr(0,dot);
140 3738 : string name=c[i].substr(dot+1);
141 1869 : if(c[i].find(".")!=string::npos) { // if it contains a dot:
142 870 : if(a=="*" && name=="*") {
143 : // Take all values from all actions
144 1 : std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>();
145 1 : if( all.empty() ) error("your input file is not telling plumed to calculate anything");
146 9 : for(unsigned j=0; j<all.size(); j++) {
147 8 : for(int k=0; k<all[j]->getNumberOfComponents(); ++k) arg.push_back(all[j]->copyOutput(k));
148 1 : }
149 869 : } else if ( name=="*") {
150 : // Take all the values from an action with a specific name
151 356 : ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(a);
152 356 : if(!action) {
153 0 : std::string str=" (hint! the actions in this ActionSet are: ";
154 0 : str+=plumed.getActionSet().getLabelList()+")";
155 0 : error("cannot find action named " + a + str);
156 : }
157 356 : if( action->getNumberOfComponents()==0 ) error("found " + a +".* indicating use all components calculated by action with label " + a + " but this action has no components");
158 356 : for(int k=0; k<action->getNumberOfComponents(); ++k) arg.push_back(action->copyOutput(k));
159 513 : } else if ( a=="*" ) {
160 : // Take components from all actions with a specific name
161 7 : std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>();
162 7 : if( all.empty() ) error("your input file is not telling plumed to calculate anything");
163 7 : unsigned nval=0;
164 31 : for(unsigned j=0; j<all.size(); j++) {
165 24 : std::string flab; flab=all[j]->getLabel() + "." + name;
166 24 : if( all[j]->exists(flab) ) { arg.push_back(all[j]->copyOutput(flab)); nval++; }
167 24 : }
168 7 : if(nval==0) error("found no actions with a component called " + name );
169 : } else {
170 : // Take values with a specific name
171 506 : ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(a);
172 506 : if(!action) {
173 0 : std::string str=" (hint! the actions in this ActionSet are: ";
174 0 : str+=plumed.getActionSet().getLabelList()+")";
175 0 : error("cannot find action named " + a +str);
176 : }
177 506 : if( !(action->exists(c[i])) ) {
178 0 : std::string str=" (hint! the components in this actions are: ";
179 0 : str+=action->getComponentsList()+")";
180 0 : error("action " + a + " has no component named " + name + str);
181 : } ;
182 506 : arg.push_back(action->copyOutput(c[i]));
183 : }
184 : } else { // if it doesn't contain a dot
185 999 : if(c[i]=="*") {
186 : // Take all values from all actions
187 51 : std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>();
188 51 : if( all.empty() ) error("your input file is not telling plumed to calculate anything");
189 397 : for(unsigned j=0; j<all.size(); j++) {
190 346 : for(int k=0; k<all[j]->getNumberOfComponents(); ++k) arg.push_back(all[j]->copyOutput(k));
191 51 : }
192 : } else {
193 948 : ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(c[i]);
194 948 : if(!action) {
195 0 : std::string str=" (hint! the actions in this ActionSet are: ";
196 0 : str+=plumed.getActionSet().getLabelList()+")";
197 0 : error("cannot find action named " + c[i] + str );
198 : }
199 948 : if( !(action->exists(c[i])) ) {
200 0 : std::string str=" (hint! the components in this actions are: ";
201 0 : str+=action->getComponentsList()+")";
202 0 : error("action " + c[i] + " has no component named " + c[i] +str);
203 : };
204 948 : arg.push_back(action->copyOutput(c[i]));
205 : }
206 1869 : }
207 : }
208 : }
209 965 : }
210 :
211 317 : void ActionWithArguments::expandArgKeywordInPDB( PDB& pdb ) {
212 317 : std::vector<std::string> pdb_remark=pdb.getRemark();
213 634 : std::vector<std::string> arg_names;
214 317 : bool found=Tools::parseVector(pdb_remark,"ARG",arg_names);
215 317 : if( found ) {
216 0 : std::vector<Value*> arg_vals;
217 0 : interpretArgumentList( arg_names, arg_vals );
218 0 : std::string new_args="ARG=" + arg_vals[0]->getName();
219 0 : for(unsigned i=1; i<arg_vals.size(); ++i) new_args = new_args + "," + arg_vals[i]->getName();
220 0 : pdb.setArgKeyword( new_args );
221 317 : }
222 317 : }
223 :
224 2052 : void ActionWithArguments::requestArguments(const vector<Value*> &arg) {
225 2052 : plumed_massert(!lockRequestArguments,"requested argument list can only be changed in the prepare() method");
226 2052 : arguments=arg;
227 2052 : clearDependencies();
228 4104 : std::string fullname,name;
229 20608 : for(unsigned i=0; i<arguments.size(); i++) {
230 18556 : fullname=arguments[i]->getName();
231 18556 : if(fullname.find(".")!=string::npos) {
232 14400 : std::size_t dot=fullname.find_first_of('.');
233 14400 : name=fullname.substr(0,dot);
234 : } else {
235 4156 : name=fullname;
236 : }
237 18556 : ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(name);
238 18556 : plumed_massert(action,"cannot find action named (in requestArguments - this is weird)" + name);
239 18556 : addDependency(action);
240 2052 : }
241 2052 : }
242 :
243 947 : ActionWithArguments::ActionWithArguments(const ActionOptions&ao):
244 : Action(ao),
245 947 : lockRequestArguments(false)
246 : {
247 947 : if( keywords.exists("ARG") && !keywords.exists("DATA") ) {
248 897 : vector<Value*> arg;
249 897 : parseArgumentList("ARG",arg);
250 :
251 897 : if(!arg.empty()) {
252 894 : log.printf(" with arguments");
253 894 : for(unsigned i=0; i<arg.size(); i++) log.printf(" %s",arg[i]->getName().c_str());
254 894 : log.printf("\n");
255 : }
256 897 : requestArguments(arg);
257 : }
258 947 : }
259 :
260 10 : void ActionWithArguments::calculateNumericalDerivatives( ActionWithValue* a ) {
261 10 : if(!a) {
262 10 : a=dynamic_cast<ActionWithValue*>(this);
263 10 : plumed_massert(a,"cannot compute numerical derivatives for an action without values");
264 : }
265 :
266 10 : const int nval=a->getNumberOfComponents();
267 10 : const int npar=arguments.size();
268 10 : std::vector<double> value (nval*npar);
269 65 : for(int i=0; i<npar; i++) {
270 55 : double arg0=arguments[i]->get();
271 55 : arguments[i]->set(arg0+sqrt(epsilon));
272 55 : a->calculate();
273 55 : arguments[i]->set(arg0);
274 215 : for(int j=0; j<nval; j++) {
275 160 : value[i*nval+j]=a->getOutputQuantity(j);
276 : }
277 : }
278 10 : a->calculate();
279 10 : a->clearDerivatives();
280 40 : for(int j=0; j<nval; j++) {
281 30 : Value* v=a->copyOutput(j);
282 30 : if( v->hasDerivatives() ) for(int i=0; i<npar; i++) v->addDerivative(i,(value[i*nval+j]-a->getOutputQuantity(j))/sqrt(epsilon));
283 10 : }
284 10 : }
285 :
286 261 : double ActionWithArguments::getProjection(unsigned i,unsigned j)const {
287 261 : plumed_massert(i<arguments.size()," making projections with an index which is too large");
288 261 : plumed_massert(j<arguments.size()," making projections with an index which is too large");
289 261 : const Value* v1=arguments[i];
290 261 : const Value* v2=arguments[j];
291 261 : return Value::projection(*v1,*v2);
292 : }
293 :
294 0 : void ActionWithArguments::addForcesOnArguments( const std::vector<double>& forces ) {
295 0 : for(unsigned i=0; i<arguments.size(); ++i) arguments[i]->addForce( forces[i] );
296 0 : }
297 :
298 2523 : }
|
