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 32621 : void ActionWithValue::registerKeywords(Keywords& keys) {
33 65242 : 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 32621 : keys.addFlag("NUMERICAL_DERIVATIVES", false, "calculate the derivatives for these quantities numerically");
39 32621 : keys.add("hidden","HAS_VALUES","this is used in json output to determine those actions that have values");
40 32621 : }
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 854 : void ActionWithValue::useCustomisableComponents(Keywords& keys) {
48 1708 : if( !keys.outputComponentExists(".#!custom") ) {
49 1708 : 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 854 : 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 854 : }
55 :
56 31334 : ActionWithValue::ActionWithValue(const ActionOptions&ao):
57 : Action(ao),
58 31334 : firststep(true),
59 31334 : noderiv(true),
60 31334 : numericalDerivatives(false) {
61 31334 : if( keywords.exists("NUMERICAL_DERIVATIVES") ) {
62 32570 : parseFlag("NUMERICAL_DERIVATIVES",numericalDerivatives);
63 : }
64 31334 : if(!keywords.exists("NO_ACTION_LOG") && numericalDerivatives) {
65 66 : log.printf(" using numerical derivatives\n");
66 : }
67 31334 : }
68 :
69 31334 : ActionWithValue::~ActionWithValue() {
70 : // empty destructor to delete unique_ptr
71 31334 : }
72 :
73 2428773 : void ActionWithValue::clearInputForces( const bool& force ) {
74 5712804 : for(unsigned i=0; i<values.size(); i++) {
75 : values[i]->clearInputForce();
76 : }
77 2428773 : }
78 :
79 1893240 : void ActionWithValue::clearDerivatives( const bool& force ) {
80 1893240 : unsigned nt = OpenMP::getNumThreads();
81 1893240 : #pragma omp parallel num_threads(nt)
82 : {
83 : #pragma omp for
84 : for(unsigned i=0; i<values.size(); i++) {
85 : values[i]->clearDerivatives();
86 : }
87 : }
88 1893240 : }
89 :
90 : // -- These are the routine for copying the value pointers to other classes -- //
91 :
92 11841535 : bool ActionWithValue::exists( const std::string& name ) const {
93 107518560 : for(unsigned i=0; i<values.size(); ++i) {
94 95701862 : if (values[i]->name==name) {
95 : return true;
96 : }
97 : }
98 : return false;
99 : }
100 :
101 19 : void ActionWithValue::getMatrixColumnTitles( std::vector<std::string>& argnames ) const {
102 19 : plumed_assert( getNumberOfComponents()==1 && getConstPntrToComponent(0)->getRank()==2 );
103 19 : unsigned nargs = getConstPntrToComponent(0)->getShape()[1];
104 19 : std::string aname = getConstPntrToComponent(0)->getName();
105 206 : for(unsigned j=0; j<nargs; ++j) {
106 : std::string nn;
107 187 : Tools::convert( j+1, nn );
108 374 : argnames.push_back( aname + "." + nn );
109 : }
110 19 : }
111 :
112 33372862 : Value* ActionWithValue::copyOutput( const std::string& name ) const {
113 111339123 : for(unsigned i=0; i<values.size(); ++i) {
114 111339123 : if (values[i]->name==name) {
115 33372862 : return values[i].get();
116 : }
117 : }
118 0 : plumed_merror("there is no pointer with name " + name);
119 : }
120 :
121 1155329 : Value* ActionWithValue::copyOutput( const unsigned& n ) const {
122 1155329 : plumed_massert(n<values.size(),"you have requested a pointer that is out of bounds");
123 1155329 : return values[n].get();
124 : }
125 :
126 : // -- HERE WE HAVE THE STUFF FOR THE DEFAULT VALUE -- //
127 :
128 12873 : void ActionWithValue::addValue( const std::vector<unsigned>& shape ) {
129 25746 : if( !keywords.outputComponentExists(".#!value") ) {
130 0 : warning("documentation for the value calculated by this action has not been included");
131 : } else {
132 25746 : plumed_massert( keywords.componentHasCorrectType(".#!value",shape.size(),false), "documentation for type of value is incorrect");
133 : }
134 12873 : plumed_massert(values.empty(),"You have already added the default value for this action");
135 12873 : values.emplace_back(Tools::make_unique<Value>(this,getLabel(), false, shape ) );
136 12873 : }
137 :
138 5145 : void ActionWithValue::addValueWithDerivatives( const std::vector<unsigned>& shape ) {
139 10290 : if( !keywords.outputComponentExists(".#!value") ) {
140 4 : warning("documentation for the value calculated by this action has not been included");
141 : } else {
142 10286 : plumed_massert( keywords.componentHasCorrectType(".#!value",shape.size(),true), "documentation for type of value is incorrect");
143 : }
144 5145 : plumed_massert(values.empty(),"You have already added the default value for this action");
145 5145 : values.emplace_back(Tools::make_unique<Value>(this,getLabel(), true, shape ) );
146 5145 : }
147 :
148 17168 : void ActionWithValue::setNotPeriodic() {
149 17168 : plumed_massert(values.size()==1,"The number of components is not equal to one");
150 17168 : plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default");
151 17168 : values[0]->min=0;
152 17168 : values[0]->max=0;
153 17168 : values[0]->setupPeriodicity();
154 17168 : }
155 :
156 849 : void ActionWithValue::setPeriodic( const std::string& min, const std::string& max ) {
157 849 : plumed_massert(values.size()==1,"The number of components is not equal to one");
158 849 : plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default");
159 849 : values[0]->setDomain( min, max );
160 849 : }
161 :
162 : // -- HERE WE HAVE THE STUFF FOR NAMED VALUES / COMPONENTS -- //
163 :
164 45824 : void ActionWithValue::addComponent( const std::string& name, const std::vector<unsigned>& shape ) {
165 45824 : if( !keywords.outputComponentExists(name) ) {
166 0 : plumed_merror("a description of component " + name + " has not been added to the manual. Components should be registered like keywords in "
167 : "registerKeywords as described in the developer docs.");
168 : }
169 45824 : plumed_massert( keywords.componentHasCorrectType(name,shape.size(),false), "documentation for type of component " + name + " is incorrect");
170 : std::string thename;
171 91648 : thename=getLabel() + "." + name;
172 18798881 : for(unsigned i=0; i<values.size(); ++i) {
173 18753057 : plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components");
174 18753057 : plumed_massert(values[i]->name!=thename,"there is already a value with this name: "+thename);
175 18753057 : plumed_massert(values[i]->name!=thename&&name!="bias","Since PLUMED 2.3 the component 'bias' is automatically added to all biases by the general constructor!\n"
176 : "Remove the line addComponent(\"bias\") from your bias.");
177 : }
178 45824 : values.emplace_back(Tools::make_unique<Value>(this,thename, false, shape ) );
179 45824 : std::string msg=" added component to this action: "+thename+" \n";
180 45824 : log.printf(msg.c_str());
181 45824 : }
182 :
183 41417 : void ActionWithValue::addComponentWithDerivatives( const std::string& name, const std::vector<unsigned>& shape ) {
184 41417 : if( !keywords.outputComponentExists(name) ) {
185 0 : plumed_merror("a description of component " + name + " has not been added to the manual. Components should be registered like keywords in "
186 : "registerKeywords as described in the developer doc.");
187 : }
188 41417 : plumed_massert( keywords.componentHasCorrectType(name,shape.size(),true), "documentation for type of component " + name + " is incorrect");
189 : std::string thename;
190 82834 : thename=getLabel() + "." + name;
191 2504025 : for(unsigned i=0; i<values.size(); ++i) {
192 2462608 : plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components");
193 2462608 : plumed_massert(values[i]->name!=thename,"there is already a value with this name: "+thename);
194 2462608 : plumed_massert(values[i]->name!=thename&&name!="bias","Since PLUMED 2.3 the component 'bias' is automatically added to all biases by the general constructor!\n"
195 : "Remove the line addComponentWithDerivatives(\"bias\") from your bias.");
196 : }
197 41417 : values.emplace_back(Tools::make_unique<Value>(this,thename, true, shape ) );
198 41417 : std::string msg=" added component to this action: "+thename+" \n";
199 41417 : log.printf(msg.c_str());
200 41417 : }
201 :
202 20 : std::string ActionWithValue::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
203 40 : if( keys.outputComponentExists(".#!custom") ) {
204 0 : return "a quantity calculated by the action " + getName() + " with label " + getLabel();
205 : }
206 20 : std::size_t und=cname.find_last_of("_");
207 20 : std::size_t hyph=cname.find_first_of("-");
208 20 : if( und!=std::string::npos ) {
209 0 : return keys.getOutputComponentDescription(cname.substr(und)) + " This particular component measures this quantity for the input CV named " + cname.substr(0,und);
210 : }
211 20 : if( hyph!=std::string::npos ) {
212 14 : return keys.getOutputComponentDescription(cname.substr(0,hyph)) + " This is the " + cname.substr(hyph+1) + "th of these quantities";
213 : }
214 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" );
215 13 : return keys.getOutputComponentDescription( cname );
216 : }
217 :
218 11816085 : int ActionWithValue::getComponent( const std::string& name ) const {
219 11816085 : plumed_massert( !exists( getLabel() ), "You should not be calling this routine if you are using a value");
220 : std::string thename;
221 23632170 : thename=getLabel() + "." + name;
222 65541339 : for(unsigned i=0; i<values.size(); ++i) {
223 65541339 : if (values[i]->name==thename) {
224 11816085 : return i;
225 : }
226 : }
227 0 : plumed_merror("there is no component with name " + name);
228 : }
229 :
230 0 : std::string ActionWithValue::getComponentsList( ) const {
231 : std::string complist;
232 0 : for(unsigned i=0; i<values.size(); ++i) {
233 0 : complist+=values[i]->name+" ";
234 : }
235 0 : return complist;
236 : }
237 :
238 24076 : std::vector<std::string> ActionWithValue::getComponentsVector( ) const {
239 : std::vector<std::string> complist;
240 349847 : for(unsigned i=0; i<values.size(); ++i) {
241 325771 : complist.push_back(values[i]->name);
242 : }
243 24076 : return complist;
244 0 : }
245 :
246 84021 : void ActionWithValue::componentIsNotPeriodic( const std::string& name ) {
247 84021 : int kk=getComponent(name);
248 84021 : values[kk]->min=0;
249 84021 : values[kk]->max=0;
250 84021 : values[kk]->setupPeriodicity();
251 84021 : }
252 :
253 139 : void ActionWithValue::componentIsPeriodic( const std::string& name, const std::string& min, const std::string& max ) {
254 139 : int kk=getComponent(name);
255 139 : values[kk]->setDomain(min,max);
256 139 : }
257 :
258 1858577 : void ActionWithValue::setGradientsIfNeeded() {
259 3717154 : if(isOptionOn("GRADIENTS")) {
260 201 : ActionAtomistic* aa=castToActionAtomistic();
261 201 : if(aa) {
262 308 : for(unsigned i=0; i<values.size(); i++) {
263 190 : unsigned start=0;
264 : values[i]->gradients.clear();
265 190 : values[i]->setGradients( aa, start );
266 : }
267 : } else {
268 83 : ActionWithArguments* aarg = castToActionWithArguments();
269 83 : if( !aarg ) {
270 0 : plumed_merror( "failing in " + getLabel() );
271 : }
272 166 : for(unsigned i=0; i<values.size(); i++) {
273 83 : unsigned start=0;
274 : values[i]->gradients.clear();
275 83 : aarg->setGradients( values[i].get(), start );
276 : }
277 : }
278 : }
279 1858577 : }
280 :
281 1641897 : void ActionWithValue::turnOnDerivatives() {
282 : // Turn on the derivatives
283 1641897 : noderiv=false;
284 : // Resize the derivatives
285 5997559 : for(unsigned i=0; i<values.size(); ++i) {
286 4355662 : values[i]->resizeDerivatives( getNumberOfDerivatives() );
287 : }
288 : // And turn on the derivatives in all actions on which we are dependent
289 3277708 : for(unsigned i=0; i<getDependencies().size(); ++i) {
290 1635811 : ActionWithValue* vv=getDependencies()[i]->castToActionWithValue();
291 1635811 : if(vv) {
292 1635811 : vv->turnOnDerivatives();
293 : }
294 : }
295 1641897 : }
296 :
297 11731925 : Value* ActionWithValue::getPntrToComponent( const std::string& name ) {
298 11731925 : int kk=getComponent(name);
299 11731925 : return values[kk].get();
300 : }
301 :
302 1199389313 : const Value* ActionWithValue::getConstPntrToComponent(int n) const {
303 : plumed_dbg_massert(n<values.size(),"you have requested a pointer that is out of bounds");
304 1199389313 : return values[n].get();
305 : }
306 :
307 87512337 : Value* ActionWithValue::getPntrToComponent( int n ) {
308 : plumed_dbg_massert(n<values.size(),"you have requested a pointer that is out of bounds");
309 87512337 : return values[n].get();
310 : }
311 :
312 4655285 : bool ActionWithValue::calculateOnUpdate() {
313 4655285 : if( firststep ) {
314 205909 : ActionWithArguments* aa=dynamic_cast<ActionWithArguments*>(this);
315 205909 : if(aa) {
316 184033 : const std::vector<Value*> & args(aa->getArguments());
317 1261742 : for(const auto & p : args ) {
318 1078114 : if( p->calculateOnUpdate() ) {
319 851 : for(unsigned i=0; i<values.size(); ++i) {
320 892 : values[i]->setValType("calcFromAverage");
321 : }
322 : break;
323 : }
324 : }
325 : }
326 205909 : firststep=false;
327 : }
328 11148890 : for(unsigned i=0; i<values.size(); ++i) {
329 6502445 : if( values[i]->calculateOnUpdate() ) {
330 : return true;
331 : }
332 : }
333 : return false;
334 : }
335 :
336 502945 : bool ActionWithValue::checkForForces() {
337 502945 : const unsigned ncp=getNumberOfComponents();
338 502945 : unsigned nder=getNumberOfDerivatives();
339 502945 : if( ncp==0 || nder==0 ) {
340 : return false;
341 : }
342 :
343 : unsigned nvalsWithForce=0;
344 501799 : valsToForce.resize(ncp);
345 1402530 : for(unsigned i=0; i<ncp; ++i) {
346 900731 : if( values[i]->hasForce && !values[i]->isConstant() ) {
347 292733 : valsToForce[nvalsWithForce]=i;
348 292733 : nvalsWithForce++;
349 : }
350 : }
351 501799 : if( nvalsWithForce==0 ) {
352 : return false;
353 : }
354 :
355 : // Make sure forces to apply is empty of forces
356 241994 : if( forcesForApply.size()!=nder ) {
357 3674 : forcesForApply.resize( nder );
358 : }
359 : std::fill(forcesForApply.begin(),forcesForApply.end(),0);
360 :
361 : unsigned stride=1;
362 : unsigned rank=0;
363 241994 : if(ncp>4*comm.Get_size()) {
364 22970 : stride=comm.Get_size();
365 22970 : rank=comm.Get_rank();
366 : }
367 :
368 241994 : unsigned nt=OpenMP::getNumThreads();
369 241994 : if(nt>ncp/(4*stride)) {
370 : nt=1;
371 : }
372 :
373 241994 : #pragma omp parallel num_threads(nt)
374 : {
375 : std::vector<double> omp_f;
376 : if( nt>1 ) {
377 : omp_f.resize(nder,0);
378 : }
379 : #pragma omp for
380 : for(unsigned i=rank; i<nvalsWithForce; i+=stride) {
381 : double ff=values[valsToForce[i]]->inputForce[0];
382 : std::vector<double> & thisderiv( values[valsToForce[i]]->data );
383 : int nn=nder;
384 : int one1=1;
385 : int one2=1;
386 : if( nt>1 ) {
387 : plumed_blas_daxpy(&nn,&ff,thisderiv.data()+1,&one1,omp_f.data(),&one2);
388 : } else {
389 : plumed_blas_daxpy(&nn,&ff,thisderiv.data()+1,&one1,forcesForApply.data(),&one2);
390 : }
391 : // if( nt>1 ) for(unsigned j=0; j<nder; ++j) omp_f[j] += ff*thisderiv[1+j];
392 : //else for(unsigned j=0; j<nder; ++j) forcesForApply[j] += ff*thisderiv[1+j];
393 : }
394 : #pragma omp critical
395 : {
396 : if( nt>1 ) {
397 : int nn=forcesForApply.size();
398 : double one0=1.0;
399 : int one1=1;
400 : int one2=1;
401 : plumed_blas_daxpy(&nn,&one0,omp_f.data(),&one1,forcesForApply.data(),&one2);
402 : }
403 : // for(unsigned j=0; j<forcesForApply.size(); ++j) {
404 : // forcesForApply[j]+=omp_f[j];
405 : // }
406 : }
407 : }
408 :
409 241994 : if(ncp>4*comm.Get_size()) {
410 22970 : comm.Sum(&forcesForApply[0],nder);
411 : }
412 : return true;
413 : }
414 :
415 : }
|
