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 : #include "blas/blas.h"
29 :
30 : namespace PLMD {
31 :
32 30576 : void ActionWithValue::registerKeywords(Keywords& keys) {
33 61152 : keys.setComponentsIntroduction("By default the value of the calculated quantity can be referenced elsewhere in the "
34 : "input file by using the label of the action. Alternatively this Action can be used "
35 : "to calculate the following quantities by employing the keywords listed "
36 : "below. These quantities can be referenced elsewhere in the input by using this Action's "
37 : "label followed by a dot and the name of the quantity required from the list below.");
38 30576 : keys.addFlag("NUMERICAL_DERIVATIVES", false, "calculate the derivatives for these quantities numerically");
39 30576 : keys.add("hidden","HAS_VALUES","this is used in json output to determine those actions that have values");
40 30576 : }
41 :
42 0 : void ActionWithValue::noAnalyticalDerivatives(Keywords& keys) {
43 0 : keys.remove("NUMERICAL_DERIVATIVES");
44 0 : keys.addFlag("NUMERICAL_DERIVATIVES",false,"analytical derivatives are not implemented for this keyword so numerical derivatives are always used");
45 0 : }
46 :
47 2347 : void ActionWithValue::useCustomisableComponents(Keywords& keys) {
48 4694 : if( !keys.outputComponentExists(".#!custom") ) {
49 3218 : 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");
50 : }
51 2347 : keys.setComponentsIntroduction("The names of the components in this action can be customized by the user in the "
52 : "actions input file. However, in addition to the components that can be customized the "
53 : "following quantities will always be output");
54 2347 : }
55 :
56 30583 : ActionWithValue::ActionWithValue(const ActionOptions&ao):
57 : Action(ao),
58 30583 : firststep(true),
59 30583 : noderiv(true),
60 30583 : numericalDerivatives(false) {
61 30583 : if( keywords.exists("NUMERICAL_DERIVATIVES") ) {
62 13782 : parseFlag("NUMERICAL_DERIVATIVES",numericalDerivatives);
63 : }
64 30583 : if(!keywords.exists("NO_ACTION_LOG") && numericalDerivatives) {
65 66 : log.printf(" using numerical derivatives\n");
66 : }
67 30583 : }
68 :
69 30583 : ActionWithValue::~ActionWithValue() {
70 : // empty destructor to delete unique_ptr
71 30583 : }
72 :
73 2413401 : void ActionWithValue::clearInputForces( const bool& force ) {
74 5707680 : for(unsigned i=0; i<values.size(); i++) {
75 : values[i]->clearInputForce();
76 : }
77 2413401 : }
78 :
79 1851182 : void ActionWithValue::clearDerivatives( const bool& force ) {
80 : #ifdef _OPENMP
81 : //nt is unused if openmp is not declared
82 1851182 : const unsigned nt=OpenMP::getNumThreads();
83 : #endif //_OPENMP
84 1851182 : #pragma omp parallel num_threads(nt)
85 : {
86 : #pragma omp for
87 : for(unsigned i=0; i<values.size(); i++) {
88 : values[i]->clearDerivatives();
89 : }
90 : }
91 1851182 : }
92 :
93 : // -- These are the routine for copying the value pointers to other classes -- //
94 :
95 11840099 : bool ActionWithValue::exists( const std::string& valname ) const {
96 107508540 : for(unsigned i=0; i<values.size(); ++i) {
97 95690604 : if (values[i]->name==valname) {
98 : return true;
99 : }
100 : }
101 : return false;
102 : }
103 :
104 19 : void ActionWithValue::getMatrixColumnTitles( std::vector<std::string>& argnames ) const {
105 19 : plumed_assert( getNumberOfComponents()==1 && getConstPntrToComponent(0)->getRank()==2 );
106 19 : unsigned nargs = getConstPntrToComponent(0)->getShape()[1];
107 19 : std::string aname = getConstPntrToComponent(0)->getName();
108 206 : for(unsigned j=0; j<nargs; ++j) {
109 : std::string nn;
110 187 : Tools::convert( j+1, nn );
111 374 : argnames.push_back( aname + "." + nn );
112 : }
113 19 : }
114 :
115 33412553 : Value* ActionWithValue::copyOutput( const std::string& valname ) const {
116 111471179 : for(unsigned i=0; i<values.size(); ++i) {
117 111471179 : if (values[i]->name==valname) {
118 33412553 : return values[i].get();
119 : }
120 : }
121 0 : plumed_merror("there is no pointer with name " + valname);
122 : }
123 :
124 1338559 : Value* ActionWithValue::copyOutput( const unsigned& n ) const {
125 1338559 : plumed_massert(n<values.size(),"you have requested a pointer that is out of bounds");
126 1338559 : return values[n].get();
127 : }
128 :
129 : // -- HERE WE HAVE THE STUFF FOR THE DEFAULT VALUE -- //
130 :
131 13643 : void ActionWithValue::addValue( const std::vector<std::size_t>& shape ) {
132 27286 : if( !keywords.outputComponentExists(".#!value") ) {
133 40 : warning("documentation for the value calculated by this action has not been included");
134 : } else {
135 27246 : plumed_massert( keywords.componentHasCorrectType(".#!value",shape.size(),false), "documentation for type of value is incorrect");
136 : }
137 13643 : plumed_massert(values.empty(),"You have already added the default value for this action");
138 13643 : values.emplace_back(Tools::make_unique<Value>(this,getLabel(), false, shape ) );
139 13643 : }
140 :
141 3341 : void ActionWithValue::addValueWithDerivatives( const std::vector<std::size_t>& shape ) {
142 6682 : if( !keywords.outputComponentExists(".#!value") ) {
143 4 : warning("documentation for the value calculated by this action has not been included");
144 : } else {
145 6678 : plumed_massert( keywords.componentHasCorrectType(".#!value",shape.size(),true), "documentation for type of value is incorrect");
146 : }
147 3341 : plumed_massert(values.empty(),"You have already added the default value for this action");
148 3341 : values.emplace_back(Tools::make_unique<Value>(this,getLabel(), true, shape ) );
149 3341 : }
150 :
151 15870 : void ActionWithValue::setNotPeriodic() {
152 15870 : plumed_massert(values.size()==1,"The number of components is not equal to one");
153 15870 : plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default");
154 15870 : values[0]->min=0;
155 15870 : values[0]->max=0;
156 15870 : values[0]->setupPeriodicity();
157 15870 : }
158 :
159 843 : void ActionWithValue::setPeriodic( const std::string& min, const std::string& max ) {
160 843 : plumed_massert(values.size()==1,"The number of components is not equal to one");
161 843 : plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default");
162 843 : values[0]->setDomain( min, max );
163 843 : }
164 :
165 : // -- HERE WE HAVE THE STUFF FOR NAMED VALUES / COMPONENTS -- //
166 :
167 46326 : void ActionWithValue::addComponent( const std::string& valname, const std::vector<std::size_t>& shape ) {
168 46326 : if( !keywords.outputComponentExists(valname) ) {
169 0 : plumed_merror("a description of component " + valname + " has not been added to the manual. Components should be registered like keywords in "
170 : "registerKeywords as described in the developer docs.");
171 : }
172 46326 : plumed_massert( keywords.componentHasCorrectType(valname,shape.size(),false), "documentation for type of component " + valname + " is incorrect");
173 : std::string thename;
174 92652 : thename=getLabel() + "." + valname;
175 18800954 : for(unsigned i=0; i<values.size(); ++i) {
176 18754628 : plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components");
177 18754628 : plumed_massert(values[i]->name!=thename,"there is already a value with this name: "+thename);
178 18754628 : 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"
179 : "Remove the line addComponent(\"bias\") from your bias.");
180 : }
181 46326 : values.emplace_back(Tools::make_unique<Value>(this,thename, false, shape ) );
182 46326 : log.printf(" added component to this action: %s \n", thename.c_str() );
183 46326 : }
184 :
185 42057 : void ActionWithValue::addComponentWithDerivatives( const std::string& valname, const std::vector<std::size_t>& shape ) {
186 42057 : if( !keywords.outputComponentExists(valname) ) {
187 0 : plumed_merror("a description of component " + valname + " has not been added to the manual. Components should be registered like keywords in "
188 : "registerKeywords as described in the developer doc.");
189 : }
190 42057 : plumed_massert( keywords.componentHasCorrectType(valname,shape.size(),true), "documentation for type of component " + valname + " is incorrect");
191 : std::string thename;
192 84114 : thename=getLabel() + "." + valname;
193 2505241 : for(unsigned i=0; i<values.size(); ++i) {
194 2463184 : plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components");
195 2463184 : plumed_massert(values[i]->name!=thename,"there is already a value with this name: "+thename);
196 2463184 : 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"
197 : "Remove the line addComponentWithDerivatives(\"bias\") from your bias.");
198 : }
199 42057 : values.emplace_back(Tools::make_unique<Value>(this,thename, true, shape ) );
200 42057 : log.printf(" added component to this action: %s \n", thename.c_str() );
201 42057 : }
202 :
203 20 : std::string ActionWithValue::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
204 40 : if( keys.outputComponentExists(".#!custom") ) {
205 0 : return "a quantity calculated by the action " + getName() + " with label " + getLabel();
206 : }
207 20 : std::size_t und=cname.find_last_of("_");
208 20 : std::size_t hyph=cname.find_first_of("-");
209 20 : if( und!=std::string::npos ) {
210 0 : return keys.getOutputComponentDescription(cname.substr(und)) + " This particular component measures this quantity for the input CV named " + cname.substr(0,und);
211 : }
212 20 : if( hyph!=std::string::npos ) {
213 14 : return keys.getOutputComponentDescription(cname.substr(0,hyph)) + " This is the " + cname.substr(hyph+1) + "th of these quantities";
214 : }
215 13 : 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" );
216 13 : return keys.getOutputComponentDescription( cname );
217 : }
218 :
219 11817306 : int ActionWithValue::getComponent( const std::string& valname ) const {
220 11817306 : plumed_massert( !exists( getLabel() ), "You should not be calling this routine if you are using a value");
221 : std::string thename;
222 23634612 : thename=getLabel() + "." + valname;
223 65539434 : for(unsigned i=0; i<values.size(); ++i) {
224 65539434 : if (values[i]->name==thename) {
225 11817306 : return i;
226 : }
227 : }
228 0 : plumed_merror("there is no component with name " + valname);
229 : }
230 :
231 0 : std::string ActionWithValue::getComponentsList( ) const {
232 : std::string complist;
233 0 : for(unsigned i=0; i<values.size(); ++i) {
234 0 : complist+=values[i]->name+" ";
235 : }
236 0 : return complist;
237 : }
238 :
239 24076 : std::vector<std::string> ActionWithValue::getComponentsVector( ) const {
240 : std::vector<std::string> complist;
241 349847 : for(unsigned i=0; i<values.size(); ++i) {
242 325771 : complist.push_back(values[i]->name);
243 : }
244 24076 : return complist;
245 0 : }
246 :
247 84450 : void ActionWithValue::componentIsNotPeriodic( const std::string& valname ) {
248 84450 : int kk=getComponent(valname);
249 84450 : values[kk]->min=0;
250 84450 : values[kk]->max=0;
251 84450 : values[kk]->setupPeriodicity();
252 84450 : }
253 :
254 139 : void ActionWithValue::componentIsPeriodic( const std::string& valname, const std::string& min, const std::string& max ) {
255 139 : int kk=getComponent(valname);
256 139 : values[kk]->setDomain(min,max);
257 139 : }
258 :
259 1886605 : void ActionWithValue::setGradientsIfNeeded() {
260 3773210 : if(isOptionOn("GRADIENTS")) {
261 201 : ActionAtomistic* aa=castToActionAtomistic();
262 201 : if(aa) {
263 308 : for(unsigned i=0; i<values.size(); i++) {
264 190 : unsigned start=0;
265 : values[i]->gradients.clear();
266 190 : values[i]->setGradients( aa, start );
267 : }
268 : } else {
269 83 : ActionWithArguments* aarg = castToActionWithArguments();
270 83 : if( !aarg ) {
271 0 : plumed_merror( "failing in " + getLabel() );
272 : }
273 166 : for(unsigned i=0; i<values.size(); i++) {
274 83 : unsigned start=0;
275 : values[i]->gradients.clear();
276 83 : aarg->setGradients( values[i].get(), start );
277 : }
278 : }
279 : }
280 1886605 : }
281 :
282 1652653 : void ActionWithValue::turnOnDerivatives() {
283 : // Turn on the derivatives
284 1652653 : noderiv=false;
285 : // Resize the derivatives
286 6011471 : for(unsigned i=0; i<values.size(); ++i) {
287 4358818 : values[i]->resizeDerivatives( getNumberOfDerivatives() );
288 : }
289 : // And turn on the derivatives in all actions on which we are dependent
290 3299196 : for(unsigned i=0; i<getDependencies().size(); ++i) {
291 1646543 : ActionWithValue* vv=getDependencies()[i]->castToActionWithValue();
292 1646543 : if(vv) {
293 1646543 : vv->turnOnDerivatives();
294 : }
295 : }
296 1652653 : }
297 :
298 11732717 : Value* ActionWithValue::getPntrToComponent( const std::string& valname ) {
299 11732717 : int kk=getComponent(valname);
300 11732717 : return values[kk].get();
301 : }
302 :
303 39657375 : const Value* ActionWithValue::getConstPntrToComponent(unsigned n) const {
304 : plumed_dbg_massert(n<values.size(),"you have requested a pointer that is out of bounds");
305 39657375 : return values[n].get();
306 : }
307 :
308 139963366 : Value* ActionWithValue::getPntrToComponent( unsigned n ) {
309 : plumed_dbg_massert(n<values.size(),"you have requested a pointer that is out of bounds");
310 139963366 : return values[n].get();
311 : }
312 :
313 4719781 : bool ActionWithValue::calculateOnUpdate() {
314 4719781 : if( firststep ) {
315 205126 : ActionWithArguments* aa=dynamic_cast<ActionWithArguments*>(this);
316 205126 : if(aa) {
317 182752 : const std::vector<Value*> & args(aa->getArguments());
318 1259016 : for(const auto & p : args ) {
319 1076671 : if( p->calculateOnUpdate() ) {
320 855 : for(unsigned i=0; i<values.size(); ++i) {
321 896 : values[i]->setValType("calcFromAverage");
322 : }
323 : break;
324 : }
325 : }
326 : }
327 205126 : firststep=false;
328 : }
329 11325582 : for(unsigned i=0; i<values.size(); ++i) {
330 6614619 : if( values[i]->calculateOnUpdate() ) {
331 : return true;
332 : }
333 : }
334 : return false;
335 : }
336 :
337 462923 : bool ActionWithValue::checkForForces() {
338 462923 : const unsigned ncp=getNumberOfComponents();
339 462923 : unsigned nder=getNumberOfDerivatives();
340 462923 : if( ncp==0 || nder==0 ) {
341 : return false;
342 : }
343 :
344 : unsigned nvalsWithForce=0;
345 461777 : valsToForce.resize(ncp);
346 1338888 : for(unsigned i=0; i<ncp; ++i) {
347 877111 : if( values[i]->hasForce && !values[i]->isConstant() ) {
348 256383 : valsToForce[nvalsWithForce]=i;
349 256383 : nvalsWithForce++;
350 : }
351 : }
352 461777 : if( nvalsWithForce==0 ) {
353 : return false;
354 : }
355 :
356 : // Make sure forces to apply is empty of forces
357 205216 : if( forcesForApply.size()!=nder ) {
358 3428 : forcesForApply.resize( nder );
359 : }
360 : std::fill(forcesForApply.begin(),forcesForApply.end(),0);
361 :
362 : unsigned stride=1;
363 : unsigned rank=0;
364 205216 : if(ncp>static_cast<unsigned>(4*comm.Get_size())) {
365 22980 : stride=comm.Get_size();
366 22980 : rank=comm.Get_rank();
367 : }
368 :
369 205216 : unsigned nt=OpenMP::getNumThreads();
370 205216 : if(nt>ncp/(4*stride)) {
371 : nt=1;
372 : }
373 :
374 205216 : #pragma omp parallel num_threads(nt)
375 : {
376 : std::vector<double> omp_f;
377 : if( nt>1 ) {
378 : omp_f.resize(nder,0);
379 : }
380 : #pragma omp for
381 : for(unsigned i=rank; i<nvalsWithForce; i+=stride) {
382 : double ff=values[valsToForce[i]]->inputForce[0];
383 : std::vector<double> & thisderiv( values[valsToForce[i]]->data );
384 : int nn=nder;
385 : int one1=1;
386 : int one2=1;
387 : if( nt>1 ) {
388 : plumed_blas_daxpy(&nn,&ff,thisderiv.data()+1,&one1,omp_f.data(),&one2);
389 : } else {
390 : plumed_blas_daxpy(&nn,&ff,thisderiv.data()+1,&one1,forcesForApply.data(),&one2);
391 : }
392 : // if( nt>1 ) for(unsigned j=0; j<nder; ++j) omp_f[j] += ff*thisderiv[1+j];
393 : //else for(unsigned j=0; j<nder; ++j) forcesForApply[j] += ff*thisderiv[1+j];
394 : }
395 : #pragma omp critical
396 : {
397 : if( nt>1 ) {
398 : int nn=forcesForApply.size();
399 : double one0=1.0;
400 : int one1=1;
401 : int one2=1;
402 : plumed_blas_daxpy(&nn,&one0,omp_f.data(),&one1,forcesForApply.data(),&one2);
403 : }
404 : // for(unsigned j=0; j<forcesForApply.size(); ++j) {
405 : // forcesForApply[j]+=omp_f[j];
406 : // }
407 : }
408 : }
409 :
410 205216 : if(ncp>static_cast<unsigned>(4*comm.Get_size())) {
411 22980 : comm.Sum(&forcesForApply[0],nder);
412 : }
413 : return true;
414 : }
415 :
416 : }
|
