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 "ActionWithValue.h"
23 : #include "ActionWithArguments.h"
24 : #include "ActionAtomistic.h"
25 : #include "tools/Exception.h"
26 : #include "tools/OpenMP.h"
27 : #include "tools/Communicator.h"
28 :
29 : namespace PLMD {
30 :
31 47811 : void ActionWithValue::registerKeywords(Keywords& keys) {
32 95622 : keys.setComponentsIntroduction("By default the value of the calculated quantity can be referenced elsewhere in the "
33 : "input file by using the label of the action. Alternatively this Action can be used "
34 : "to calculate the following quantities by employing the keywords listed "
35 : "below. These quantities can be referenced elsewhere in the input by using this Action's "
36 : "label followed by a dot and the name of the quantity required from the list below.");
37 47811 : keys.addFlag("NUMERICAL_DERIVATIVES", false, "calculate the derivatives for these quantities numerically");
38 47811 : keys.add("hidden","HAS_VALUES","this is used in json output to determine those actions that have values");
39 47811 : }
40 :
41 0 : void ActionWithValue::noAnalyticalDerivatives(Keywords& keys) {
42 0 : keys.remove("NUMERICAL_DERIVATIVES");
43 0 : keys.addFlag("NUMERICAL_DERIVATIVES",false,"analytical derivatives are not implemented for this keyword so numerical derivatives are always used");
44 0 : }
45 :
46 2936 : void ActionWithValue::useCustomisableComponents(Keywords& keys) {
47 5872 : if( !keys.outputComponentExists(".#!custom") ) {
48 4268 : keys.addOutputComponent(".#!custom","default","scalar","the names of the output components for this action depend on the actions input file see the example inputs below for details");
49 : }
50 2936 : keys.setComponentsIntroduction("The names of the components in this action can be customized by the user in the "
51 : "actions input file. However, in addition to the components that can be customized the "
52 : "following quantities will always be output");
53 2936 : }
54 :
55 43453 : ActionWithValue::ActionWithValue(const ActionOptions&ao):
56 : Action(ao),
57 43453 : firststep(true),
58 43453 : noderiv(true),
59 43453 : numericalDerivatives(false) {
60 43453 : if( keywords.exists("NUMERICAL_DERIVATIVES") ) {
61 22778 : parseFlag("NUMERICAL_DERIVATIVES",numericalDerivatives);
62 : }
63 43453 : if(!keywords.exists("NO_ACTION_LOG") && numericalDerivatives) {
64 43 : log.printf(" using numerical derivatives\n");
65 : }
66 43453 : }
67 :
68 43453 : ActionWithValue::~ActionWithValue() {
69 : // empty destructor to delete unique_ptr
70 43453 : }
71 :
72 3354040 : void ActionWithValue::clearInputForces( const bool& force ) {
73 7843218 : for(unsigned i=0; i<values.size(); i++) {
74 : values[i]->clearInputForce();
75 : }
76 3354040 : }
77 :
78 2146636 : void ActionWithValue::clearDerivatives( const bool& force ) {
79 : // Opening an OpenMP region costs microseconds and this is called for every
80 : // active action at every step: use threads only when there is enough to clear.
81 : #ifdef _OPENMP
82 2146636 : unsigned nt=OpenMP::getNumThreads();
83 2146636 : if( nt>1 ) {
84 : std::size_t tot=0;
85 5584099 : for(const auto & v : values) {
86 3441303 : tot+=v->data.size();
87 : }
88 2142796 : if( values.size()<2 || tot<32768 ) {
89 : nt=1;
90 342 : } else if( nt>values.size() ) {
91 0 : nt=values.size();
92 : }
93 : }
94 2146636 : if( nt>1 ) {
95 342 : #pragma omp parallel for num_threads(nt)
96 : for(unsigned i=0; i<values.size(); i++) {
97 : values[i]->clearDerivatives();
98 : }
99 342 : return;
100 : }
101 : #endif //_OPENMP
102 5575727 : for(unsigned i=0; i<values.size(); i++) {
103 3429433 : values[i]->clearDerivatives();
104 : }
105 : }
106 :
107 : // -- These are the routine for copying the value pointers to other classes -- //
108 :
109 722866 : bool ActionWithValue::exists( const std::string& valname ) const {
110 45430713 : for(unsigned i=0; i<values.size(); ++i) {
111 44738297 : if (values[i]->name==valname) {
112 : return true;
113 : }
114 : }
115 : return false;
116 : }
117 :
118 19 : void ActionWithValue::getMatrixColumnTitles( std::vector<std::string>& argnames ) const {
119 19 : plumed_assert( getNumberOfComponents()==1 && getConstPntrToComponent(0)->getRank()==2 );
120 19 : unsigned nargs = getConstPntrToComponent(0)->getShape()[1];
121 19 : std::string aname = getConstPntrToComponent(0)->getName();
122 206 : for(unsigned j=0; j<nargs; ++j) {
123 : std::string nn;
124 187 : Tools::convert( j+1, nn );
125 374 : argnames.push_back( aname + "." + nn );
126 : }
127 19 : }
128 :
129 125508 : Value* ActionWithValue::copyOutput( const std::string& valname ) const {
130 11634934 : for(unsigned i=0; i<values.size(); ++i) {
131 11634934 : if (values[i]->name==valname) {
132 125508 : return values[i].get();
133 : }
134 : }
135 0 : plumed_merror("there is no pointer with name " + valname);
136 : }
137 :
138 1838456 : Value* ActionWithValue::copyOutput( const unsigned& n ) const {
139 1838456 : plumed_massert(n<values.size(),"you have requested a pointer that is out of bounds");
140 1838456 : return values[n].get();
141 : }
142 :
143 : // -- HERE WE HAVE THE STUFF FOR THE DEFAULT VALUE -- //
144 :
145 21273 : void ActionWithValue::addValue( const std::vector<std::size_t>& shape ) {
146 42546 : if( !keywords.outputComponentExists(".#!value") ) {
147 74 : warning("documentation for the value calculated by this action has not been included");
148 : } else {
149 42472 : plumed_massert( keywords.componentHasCorrectType(".#!value",shape.size(),false), "documentation for type of value of " + getName() + " is missing or incorrect");
150 : }
151 21273 : plumed_massert(values.empty(),"You have already added the default value for this action");
152 21273 : values.emplace_back(Tools::make_unique<Value>(this,getLabel(), false, shape ) );
153 21273 : }
154 :
155 5725 : void ActionWithValue::addValueWithDerivatives( const std::vector<std::size_t>& shape ) {
156 11450 : if( !keywords.outputComponentExists(".#!value") ) {
157 24 : warning("documentation for the value calculated by this action has not been included");
158 : } else {
159 11426 : plumed_massert( keywords.componentHasCorrectType(".#!value",shape.size(),true), "documentation for type of value of " + getName() + " is missing or incorrect");
160 : }
161 5725 : plumed_massert(values.empty(),"You have already added the default value for this action");
162 5725 : values.emplace_back(Tools::make_unique<Value>(this,getLabel(), true, shape ) );
163 5725 : }
164 :
165 21894 : void ActionWithValue::setNotPeriodic() {
166 21894 : plumed_massert(values.size()==1,"The number of components is not equal to one");
167 21894 : plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default");
168 21894 : values[0]->min=0;
169 21894 : values[0]->max=0;
170 21894 : values[0]->setupPeriodicity();
171 21894 : }
172 :
173 969 : void ActionWithValue::setPeriodic( const std::string& min, const std::string& max ) {
174 969 : plumed_massert(values.size()==1,"The number of components is not equal to one");
175 969 : plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default");
176 969 : values[0]->setDomain( min, max );
177 969 : }
178 :
179 : // -- HERE WE HAVE THE STUFF FOR NAMED VALUES / COMPONENTS -- //
180 :
181 50668 : void ActionWithValue::addComponent( const std::string& valname, const std::vector<std::size_t>& shape ) {
182 50668 : if( !keywords.outputComponentExists(valname) ) {
183 0 : plumed_merror("a description of component " + valname + " has not been added to the manual. Components should be registered like keywords in "
184 : "registerKeywords as described in the developer docs.");
185 : }
186 50668 : plumed_massert( keywords.componentHasCorrectType(valname,shape.size(),false), "documentation for type of component " + valname + " of " + getName() + " is missing or incorrect");
187 : std::string thename;
188 50668 : thename=getLabel() + "." + valname;
189 18880865 : for(unsigned i=0; i<values.size(); ++i) {
190 18830197 : plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components");
191 18830197 : plumed_massert(values[i]->name!=thename,"there is already a value with this name: "+thename);
192 18830197 : plumed_massert(values[i]->name!=thename&&valname!="bias","Since PLUMED 2.3 the component 'bias' is automatically added to all biases by the general constructor!\n"
193 : "Remove the line addComponent(\"bias\") from your bias.");
194 : }
195 50668 : values.emplace_back(Tools::make_unique<Value>(this,thename, false, shape ) );
196 50668 : log.printf(" added component to this action: %s \n", thename.c_str() );
197 50668 : }
198 :
199 45204 : void ActionWithValue::addComponentWithDerivatives( const std::string& valname, const std::vector<std::size_t>& shape ) {
200 45204 : if( !keywords.outputComponentExists(valname) ) {
201 0 : plumed_merror("a description of component " + valname + " has not been added to the manual. Components should be registered like keywords in "
202 : "registerKeywords as described in the developer doc.");
203 : }
204 45204 : plumed_massert( keywords.componentHasCorrectType(valname,shape.size(),true), "documentation for type of component " + valname + " of " + getName() + " is missing or incorrect");
205 : std::string thename;
206 45204 : thename=getLabel() + "." + valname;
207 2516940 : for(unsigned i=0; i<values.size(); ++i) {
208 2471736 : plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components");
209 2471736 : plumed_massert(values[i]->name!=thename,"there is already a value with this name: "+thename);
210 2471736 : plumed_massert(values[i]->name!=thename&&valname!="bias","Since PLUMED 2.3 the component 'bias' is automatically added to all biases by the general constructor!\n"
211 : "Remove the line addComponentWithDerivatives(\"bias\") from your bias.");
212 : }
213 45204 : values.emplace_back(Tools::make_unique<Value>(this,thename, true, shape ) );
214 45204 : log.printf(" added component to this action: %s \n", thename.c_str() );
215 45204 : }
216 :
217 30 : std::string ActionWithValue::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
218 60 : if( keys.outputComponentExists(".#!custom") ) {
219 0 : return "a quantity calculated by the action " + getName() + " with label " + getLabel();
220 : }
221 30 : std::size_t und=cname.find_last_of("_");
222 30 : std::size_t hyph=cname.find_first_of("-");
223 30 : if( und!=std::string::npos ) {
224 0 : return keys.getOutputComponentDescription(cname.substr(und)) + " This particular component measures this quantity for the input CV named " + cname.substr(0,und);
225 : }
226 30 : if( hyph!=std::string::npos ) {
227 21 : return keys.getOutputComponentDescription(cname.substr(0,hyph)) + " This is the " + cname.substr(hyph+1) + "th of these quantities";
228 : }
229 23 : plumed_massert( keys.outputComponentExists(cname), "component " + cname + " does not exist in " + keys.getDisplayName() + " if the component names are customizable then you should override this function" );
230 23 : return keys.getOutputComponentDescription( cname );
231 : }
232 :
233 691724 : int ActionWithValue::getComponent( const std::string& valname ) const {
234 691724 : plumed_massert( !exists( getLabel() ), "You should not be calling this routine if you are using a value");
235 : std::string thename;
236 691724 : thename=getLabel() + "." + valname;
237 43499207 : for(unsigned i=0; i<values.size(); ++i) {
238 43499207 : if (values[i]->name==thename) {
239 691724 : return i;
240 : }
241 : }
242 0 : plumed_merror("there is no component with name " + valname);
243 : }
244 :
245 0 : std::string ActionWithValue::getComponentsList( ) const {
246 : std::string complist;
247 0 : for(unsigned i=0; i<values.size(); ++i) {
248 0 : complist+=values[i]->name+" ";
249 : }
250 0 : return complist;
251 : }
252 :
253 24434 : std::vector<std::string> ActionWithValue::getComponentsVector( ) const {
254 24434 : std::vector<std::string> complist;
255 353041 : for(unsigned i=0; i<values.size(); ++i) {
256 328607 : complist.push_back(values[i]->name);
257 : }
258 24434 : return complist;
259 0 : }
260 :
261 91364 : void ActionWithValue::componentIsNotPeriodic( const std::string& valname ) {
262 91364 : int kk=getComponent(valname);
263 91364 : values[kk]->min=0;
264 91364 : values[kk]->max=0;
265 91364 : values[kk]->setupPeriodicity();
266 91364 : }
267 :
268 139 : void ActionWithValue::componentIsPeriodic( const std::string& valname, const std::string& min, const std::string& max ) {
269 139 : int kk=getComponent(valname);
270 139 : values[kk]->setDomain(min,max);
271 139 : }
272 :
273 2183188 : void ActionWithValue::setGradientsIfNeeded() {
274 4366376 : if(isOptionOn("GRADIENTS")) {
275 2335 : ActionAtomistic* aa=castToActionAtomistic();
276 2335 : if(aa) {
277 2718 : for(unsigned i=0; i<values.size(); i++) {
278 1395 : unsigned start=0;
279 : values[i]->gradients.clear();
280 1395 : values[i]->setGradients( aa, start );
281 : }
282 : } else {
283 1012 : ActionWithArguments* aarg = castToActionWithArguments();
284 1012 : if( !aarg ) {
285 0 : plumed_merror( "failing in " + getLabel() );
286 : }
287 2684 : for(unsigned i=0; i<values.size(); i++) {
288 1672 : unsigned start=0;
289 : values[i]->gradients.clear();
290 1672 : aarg->setGradients( values[i].get(), start );
291 : }
292 : }
293 : }
294 2183188 : }
295 :
296 1706259 : void ActionWithValue::turnOnDerivatives() {
297 : // Turn on the derivatives
298 1706259 : noderiv=false;
299 : // Resize the derivatives
300 6112953 : for(unsigned i=0; i<values.size(); ++i) {
301 4406694 : values[i]->resizeDerivatives( getNumberOfDerivatives() );
302 : }
303 : // And turn on the derivatives in all actions on which we are dependent
304 3403745 : for(unsigned i=0; i<getDependencies().size(); ++i) {
305 1697486 : ActionWithValue* vv=getDependencies()[i]->castToActionWithValue();
306 1697486 : if(vv) {
307 1697486 : vv->turnOnDerivatives();
308 : }
309 : }
310 1706259 : }
311 :
312 600221 : Value* ActionWithValue::getPntrToComponent( const std::string& valname ) {
313 600221 : int kk=getComponent(valname);
314 600221 : return values[kk].get();
315 : }
316 :
317 2455066 : const Value* ActionWithValue::getConstPntrToComponent(unsigned n) const {
318 : plumed_dbg_massert(n<values.size(),"you have requested a pointer that is out of bounds");
319 2455066 : return values[n].get();
320 : }
321 :
322 7746536 : Value* ActionWithValue::getPntrToComponent( unsigned n ) {
323 : plumed_dbg_massert(n<values.size(),"you have requested a pointer that is out of bounds");
324 7746536 : return values[n].get();
325 : }
326 :
327 6587325 : bool ActionWithValue::calculateOnUpdate() {
328 6587325 : if( firststep ) {
329 39572 : ActionWithArguments* aa=dynamic_cast<ActionWithArguments*>(this);
330 39572 : if(aa) {
331 10461 : const std::vector<Value*> & args(aa->getArguments());
332 39426 : for(const auto & p : args ) {
333 29392 : if( p->calculateOnUpdate() ) {
334 895 : for(unsigned i=0; i<values.size(); ++i) {
335 936 : values[i]->setValType("calcFromAverage");
336 : }
337 : break;
338 : }
339 : }
340 : }
341 39572 : firststep=false;
342 : }
343 15575547 : for(unsigned i=0; i<values.size(); ++i) {
344 8997084 : if( values[i]->calculateOnUpdate() ) {
345 : return true;
346 : }
347 : }
348 : return false;
349 : }
350 :
351 606214 : bool ActionWithValue::checkForForces() {
352 606214 : const unsigned ncp=getNumberOfComponents();
353 606214 : unsigned nder=getNumberOfDerivatives();
354 606214 : if( ncp==0 || nder==0 ) {
355 : return false;
356 : }
357 :
358 : unsigned nvalsWithForce=0;
359 605068 : valsToForce.resize(ncp);
360 1475760 : for(unsigned i=0; i<ncp; ++i) {
361 870692 : if( values[i]->hasForce && !values[i]->isConstant() ) {
362 363401 : valsToForce[nvalsWithForce]=i;
363 363401 : nvalsWithForce++;
364 : }
365 : }
366 605068 : if( nvalsWithForce==0 ) {
367 : return false;
368 : }
369 :
370 : // Make sure forces to apply is empty of forces
371 357614 : if( forcesForApply.size()!=nder ) {
372 3766 : forcesForApply.resize( nder );
373 : }
374 : std::fill(forcesForApply.begin(),forcesForApply.end(),0);
375 :
376 : unsigned stride=1;
377 : unsigned rank=0;
378 357614 : if(ncp>static_cast<unsigned>(4*comm.Get_size())) {
379 65 : stride=comm.Get_size();
380 65 : rank=comm.Get_rank();
381 : }
382 :
383 357614 : unsigned nt=OpenMP::getNumThreads();
384 357614 : if(nt>ncp/(4*stride)) {
385 : nt=1;
386 : }
387 :
388 : // Note: plain loops are used here on purpose instead of BLAS daxpy.
389 : // An external multithreaded BLAS (e.g. OpenBLAS, MKL) would spawn its own
390 : // threads for long vectors (nder>~10000), ignoring the number of threads
391 : // set by the MD engine and oversubscribing the cores (see issue #1121).
392 23 : if( nt==1 ) {
393 : // serial path: opening an OpenMP region would cost more than the loop for small actions
394 : double* const dest=forcesForApply.data();
395 717957 : for(unsigned i=rank; i<nvalsWithForce; i+=stride) {
396 : // a plain loop is vectorized by the compiler and is as fast as an external
397 : // blas daxpy, without the risk of calling a threaded blas within an OpenMP region
398 360366 : const double ff=values[valsToForce[i]]->inputForce[0];
399 : const double* thisderiv=values[valsToForce[i]]->data.data()+1;
400 15008514 : for(unsigned j=0; j<nder; ++j) {
401 14648148 : dest[j]+=ff*thisderiv[j];
402 : }
403 : }
404 : } else {
405 23 : #pragma omp parallel num_threads(nt)
406 : {
407 : std::vector<double> omp_f(nder,0);
408 : double* const dest=omp_f.data();
409 : #pragma omp for
410 : for(unsigned i=rank; i<nvalsWithForce; i+=stride) {
411 : const double ff=values[valsToForce[i]]->inputForce[0];
412 : const double* thisderiv=values[valsToForce[i]]->data.data()+1;
413 : for(unsigned j=0; j<nder; ++j) {
414 : dest[j]+=ff*thisderiv[j];
415 : }
416 : }
417 : #pragma omp critical
418 : {
419 : for(unsigned j=0; j<nder; ++j) {
420 : forcesForApply[j]+=omp_f[j];
421 : }
422 : }
423 : }
424 : }
425 :
426 357614 : if(ncp>static_cast<unsigned>(4*comm.Get_size())) {
427 65 : comm.Sum(&forcesForApply[0],nder);
428 : }
429 : return true;
430 : }
431 :
432 : }
|
