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 37574 : void ActionWithValue::registerKeywords(Keywords& keys) {
32 37574 : 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 75148 : keys.addFlag("NUMERICAL_DERIVATIVES", false, "calculate the derivatives for these quantities numerically");
38 75148 : keys.add("hidden","HAS_VALUES","this is used in json output to determine those actions that have values");
39 37574 : }
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 1026 : void ActionWithValue::useCustomisableComponents(Keywords& keys) {
47 2052 : if( !keys.outputComponentExists(".#!custom") ) {
48 2052 : keys.addOutputComponent(".#!custom","default","the names of the output components for this action depend on the actions input file see the example inputs below for details");
49 : }
50 1026 : 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 1026 : }
54 :
55 35143 : ActionWithValue::ActionWithValue(const ActionOptions&ao):
56 : Action(ao),
57 35143 : firststep(true),
58 35143 : noderiv(true),
59 35143 : numericalDerivatives(false) {
60 70286 : if( keywords.exists("NUMERICAL_DERIVATIVES") ) {
61 39424 : parseFlag("NUMERICAL_DERIVATIVES",numericalDerivatives);
62 : }
63 70286 : if(!keywords.exists("NO_ACTION_LOG") && numericalDerivatives) {
64 55 : log.printf(" using numerical derivatives\n");
65 : }
66 35143 : }
67 :
68 35143 : ActionWithValue::~ActionWithValue() {
69 : // empty destructor to delete unique_ptr
70 35143 : }
71 :
72 3129941 : void ActionWithValue::clearInputForces( const bool& force ) {
73 7388113 : for(unsigned i=0; i<values.size(); i++) {
74 : values[i]->clearInputForce();
75 : }
76 3129941 : }
77 :
78 2107801 : void ActionWithValue::clearDerivatives( const bool& force ) {
79 2107801 : unsigned nt = OpenMP::getNumThreads();
80 2107801 : #pragma omp parallel num_threads(nt)
81 : {
82 : #pragma omp for
83 : for(unsigned i=0; i<values.size(); i++) {
84 : values[i]->clearDerivatives();
85 : }
86 : }
87 2107801 : }
88 :
89 : // -- These are the routine for copying the value pointers to other classes -- //
90 :
91 935076 : bool ActionWithValue::exists( const std::string& name ) const {
92 59788468 : for(unsigned i=0; i<values.size(); ++i) {
93 58879307 : if (values[i]->name==name) {
94 : return true;
95 : }
96 : }
97 : return false;
98 : }
99 :
100 19 : void ActionWithValue::getMatrixColumnTitles( std::vector<std::string>& argnames ) const {
101 19 : plumed_assert( getNumberOfComponents()==1 && getConstPntrToComponent(0)->getRank()==2 );
102 19 : unsigned nargs = getConstPntrToComponent(0)->getShape()[1];
103 19 : std::string aname = getConstPntrToComponent(0)->getName();
104 206 : for(unsigned j=0; j<nargs; ++j) {
105 : std::string nn;
106 187 : Tools::convert( j+1, nn );
107 374 : argnames.push_back( aname + "." + nn );
108 : }
109 19 : }
110 :
111 33373939 : Value* ActionWithValue::copyOutput( const std::string& name ) const {
112 111340296 : for(unsigned i=0; i<values.size(); ++i) {
113 111340296 : if (values[i]->name==name) {
114 33373939 : return values[i].get();
115 : }
116 : }
117 0 : plumed_merror("there is no pointer with name " + name);
118 : }
119 :
120 1563091 : Value* ActionWithValue::copyOutput( const unsigned& n ) const {
121 1563091 : plumed_massert(n<values.size(),"you have requested a pointer that is out of bounds");
122 1563091 : return values[n].get();
123 : }
124 :
125 : // -- HERE WE HAVE THE STUFF FOR THE DEFAULT VALUE -- //
126 :
127 15364 : void ActionWithValue::addValue( const std::vector<unsigned>& shape ) {
128 30728 : if( !keywords.outputComponentExists(".#!value") ) {
129 138 : warning("documentation for the value calculated by this action has not been included");
130 : }
131 15364 : plumed_massert(values.empty(),"You have already added the default value for this action");
132 15364 : values.emplace_back(Tools::make_unique<Value>(this,getLabel(), false, shape ) );
133 15364 : }
134 :
135 5683 : void ActionWithValue::addValueWithDerivatives( const std::vector<unsigned>& shape ) {
136 11366 : if( !keywords.outputComponentExists(".#!value") ) {
137 106 : warning("documentation for the value calculated by this action has not been included");
138 : }
139 5683 : plumed_massert(values.empty(),"You have already added the default value for this action");
140 5683 : values.emplace_back(Tools::make_unique<Value>(this,getLabel(), true, shape ) );
141 5683 : }
142 :
143 20204 : void ActionWithValue::setNotPeriodic() {
144 20204 : plumed_massert(values.size()==1,"The number of components is not equal to one");
145 20204 : plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default");
146 20204 : values[0]->min=0;
147 20204 : values[0]->max=0;
148 20204 : values[0]->setupPeriodicity();
149 20204 : }
150 :
151 842 : void ActionWithValue::setPeriodic( const std::string& min, const std::string& max ) {
152 842 : plumed_massert(values.size()==1,"The number of components is not equal to one");
153 842 : plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default");
154 842 : values[0]->setDomain( min, max );
155 842 : }
156 :
157 : // -- HERE WE HAVE THE STUFF FOR NAMED VALUES / COMPONENTS -- //
158 :
159 46188 : void ActionWithValue::addComponent( const std::string& name, const std::vector<unsigned>& shape ) {
160 46188 : if( !keywords.outputComponentExists(name) ) {
161 0 : plumed_merror("a description of component " + name + " has not been added to the manual. Components should be registered like keywords in "
162 : "registerKeywords as described in the developer docs.");
163 : }
164 : std::string thename;
165 92376 : thename=getLabel() + "." + name;
166 18799781 : for(unsigned i=0; i<values.size(); ++i) {
167 18753593 : plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components");
168 18753593 : plumed_massert(values[i]->name!=thename,"there is already a value with this name: "+thename);
169 18753593 : 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"
170 : "Remove the line addComponent(\"bias\") from your bias.");
171 : }
172 46188 : values.emplace_back(Tools::make_unique<Value>(this,thename, false, shape ) );
173 46188 : std::string msg=" added component to this action: "+thename+" \n";
174 46188 : log.printf(msg.c_str());
175 46188 : }
176 :
177 42017 : void ActionWithValue::addComponentWithDerivatives( const std::string& name, const std::vector<unsigned>& shape ) {
178 42017 : if( !keywords.outputComponentExists(name) ) {
179 0 : plumed_merror("a description of component " + name + " has not been added to the manual. Components should be registered like keywords in "
180 : "registerKeywords as described in the developer doc.");
181 : }
182 : std::string thename;
183 84034 : thename=getLabel() + "." + name;
184 2504770 : for(unsigned i=0; i<values.size(); ++i) {
185 2462753 : plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components");
186 2462753 : plumed_massert(values[i]->name!=thename,"there is already a value with this name: "+thename);
187 2462753 : 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"
188 : "Remove the line addComponentWithDerivatives(\"bias\") from your bias.");
189 : }
190 42017 : values.emplace_back(Tools::make_unique<Value>(this,thename, true, shape ) );
191 42017 : std::string msg=" added component to this action: "+thename+" \n";
192 42017 : log.printf(msg.c_str());
193 42017 : }
194 :
195 23 : std::string ActionWithValue::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
196 46 : if( keys.outputComponentExists(".#!custom") ) {
197 0 : return "a quantity calculated by the action " + getName() + " with label " + getLabel();
198 : }
199 23 : std::size_t und=cname.find_last_of("_");
200 23 : std::size_t hyph=cname.find_first_of("-");
201 23 : if( und!=std::string::npos ) {
202 0 : return keys.getOutputComponentDescription(cname.substr(und)) + " This particular component measures this quantity for the input CV named " + cname.substr(0,und);
203 : }
204 23 : if( hyph!=std::string::npos ) {
205 14 : return keys.getOutputComponentDescription(cname.substr(0,hyph)) + " This is the " + cname.substr(hyph+1) + "th of these quantities";
206 : }
207 16 : 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" );
208 16 : return keys.getOutputComponentDescription( cname );
209 : }
210 :
211 908548 : int ActionWithValue::getComponent( const std::string& name ) const {
212 908548 : plumed_massert( !exists( getLabel() ), "You should not be calling this routine if you are using a value");
213 : std::string thename;
214 1817096 : thename=getLabel() + "." + name;
215 54929886 : for(unsigned i=0; i<values.size(); ++i) {
216 54929886 : if (values[i]->name==thename) {
217 908548 : return i;
218 : }
219 : }
220 0 : plumed_merror("there is no component with name " + name);
221 : }
222 :
223 0 : std::string ActionWithValue::getComponentsList( ) const {
224 : std::string complist;
225 0 : for(unsigned i=0; i<values.size(); ++i) {
226 0 : complist+=values[i]->name+" ";
227 : }
228 0 : return complist;
229 : }
230 :
231 24087 : std::vector<std::string> ActionWithValue::getComponentsVector( ) const {
232 : std::vector<std::string> complist;
233 349877 : for(unsigned i=0; i<values.size(); ++i) {
234 325790 : complist.push_back(values[i]->name);
235 : }
236 24087 : return complist;
237 0 : }
238 :
239 84707 : void ActionWithValue::componentIsNotPeriodic( const std::string& name ) {
240 84707 : int kk=getComponent(name);
241 84707 : values[kk]->min=0;
242 84707 : values[kk]->max=0;
243 84707 : values[kk]->setupPeriodicity();
244 84707 : }
245 :
246 139 : void ActionWithValue::componentIsPeriodic( const std::string& name, const std::string& min, const std::string& max ) {
247 139 : int kk=getComponent(name);
248 139 : values[kk]->setDomain(min,max);
249 139 : }
250 :
251 2074283 : void ActionWithValue::setGradientsIfNeeded() {
252 4148566 : if(isOptionOn("GRADIENTS")) {
253 201 : ActionAtomistic* aa=castToActionAtomistic();
254 201 : if(aa) {
255 308 : for(unsigned i=0; i<values.size(); i++) {
256 190 : unsigned start=0;
257 : values[i]->gradients.clear();
258 190 : values[i]->setGradients( aa, start );
259 : }
260 : } else {
261 83 : ActionWithArguments* aarg = castToActionWithArguments();
262 83 : if( !aarg ) {
263 0 : plumed_merror( "failing in " + getLabel() );
264 : }
265 166 : for(unsigned i=0; i<values.size(); i++) {
266 83 : unsigned start=0;
267 : values[i]->gradients.clear();
268 83 : aarg->setGradients( values[i].get(), start );
269 : }
270 : }
271 : }
272 2074283 : }
273 :
274 1652107 : void ActionWithValue::turnOnDerivatives() {
275 : // Turn on the derivatives
276 1652107 : noderiv=false;
277 : // Resize the derivatives
278 6014265 : for(unsigned i=0; i<values.size(); ++i) {
279 4362158 : values[i]->resizeDerivatives( getNumberOfDerivatives() );
280 : }
281 : // And turn on the derivatives in all actions on which we are dependent
282 3297495 : for(unsigned i=0; i<getDependencies().size(); ++i) {
283 1645388 : ActionWithValue* vv=getDependencies()[i]->castToActionWithValue();
284 1645388 : if(vv) {
285 1645388 : vv->turnOnDerivatives();
286 : }
287 : }
288 1652107 : }
289 :
290 823702 : Value* ActionWithValue::getPntrToComponent( const std::string& name ) {
291 823702 : int kk=getComponent(name);
292 823702 : return values[kk].get();
293 : }
294 :
295 2455688966 : const Value* ActionWithValue::getConstPntrToComponent(int n) const {
296 : plumed_dbg_massert(n<values.size(),"you have requested a pointer that is out of bounds");
297 2455688966 : return values[n].get();
298 : }
299 :
300 77022328 : Value* ActionWithValue::getPntrToComponent( int n ) {
301 : plumed_dbg_massert(n<values.size(),"you have requested a pointer that is out of bounds");
302 77022328 : return values[n].get();
303 : }
304 :
305 6058301 : bool ActionWithValue::calculateOnUpdate() {
306 6058301 : if( firststep ) {
307 209037 : ActionWithArguments* aa=dynamic_cast<ActionWithArguments*>(this);
308 209037 : if(aa) {
309 184595 : const std::vector<Value*> & args(aa->getArguments());
310 1262985 : for(const auto & p : args ) {
311 1078812 : if( p->calculateOnUpdate() ) {
312 885 : for(unsigned i=0; i<values.size(); ++i) {
313 926 : values[i]->setValType("calcFromAverage");
314 : }
315 : break;
316 : }
317 : }
318 : }
319 209037 : firststep=false;
320 : }
321 14500748 : for(unsigned i=0; i<values.size(); ++i) {
322 8451321 : if( values[i]->calculateOnUpdate() ) {
323 : return true;
324 : }
325 : }
326 : return false;
327 : }
328 :
329 636149 : bool ActionWithValue::checkForForces() {
330 636149 : const unsigned ncp=getNumberOfComponents();
331 636149 : unsigned nder=getNumberOfDerivatives();
332 636149 : if( ncp==0 || nder==0 ) {
333 : return false;
334 : }
335 :
336 : unsigned nvalsWithForce=0;
337 635003 : valsToForce.resize(ncp);
338 1670946 : for(unsigned i=0; i<ncp; ++i) {
339 1035943 : if( values[i]->hasForce && !values[i]->isConstant() ) {
340 424647 : valsToForce[nvalsWithForce]=i;
341 424647 : nvalsWithForce++;
342 : }
343 : }
344 635003 : if( nvalsWithForce==0 ) {
345 : return false;
346 : }
347 :
348 : // Make sure forces to apply is empty of forces
349 373908 : if( forcesForApply.size()!=nder ) {
350 4189 : forcesForApply.resize( nder );
351 : }
352 : std::fill(forcesForApply.begin(),forcesForApply.end(),0);
353 :
354 : unsigned stride=1;
355 : unsigned rank=0;
356 373908 : if(ncp>4*comm.Get_size()) {
357 22970 : stride=comm.Get_size();
358 22970 : rank=comm.Get_rank();
359 : }
360 :
361 373908 : unsigned nt=OpenMP::getNumThreads();
362 373908 : if(nt>ncp/(4*stride)) {
363 : nt=1;
364 : }
365 :
366 : // Note: plain loops are used here on purpose instead of BLAS daxpy.
367 : // An external multithreaded BLAS (e.g. OpenBLAS, MKL) would spawn its own
368 : // threads for long vectors (nder>~10000), ignoring the number of threads
369 : // set by the MD engine and oversubscribing the cores (see issue #1121).
370 373908 : #pragma omp parallel num_threads(nt)
371 : {
372 : std::vector<double> omp_f;
373 : if( nt>1 ) {
374 : omp_f.resize(nder,0);
375 : }
376 : double* const dest = ( nt>1 ? omp_f.data() : forcesForApply.data() );
377 : #pragma omp for
378 : for(unsigned i=rank; i<nvalsWithForce; i+=stride) {
379 : const double ff=values[valsToForce[i]]->inputForce[0];
380 : const double* const thisderiv = values[valsToForce[i]]->data.data()+1;
381 : #pragma omp simd
382 : for(unsigned j=0; j<nder; ++j) {
383 : dest[j] += ff*thisderiv[j];
384 : }
385 : }
386 : if( nt>1 ) {
387 : #pragma omp critical
388 : {
389 : double* const f = forcesForApply.data();
390 : const double* const src = omp_f.data();
391 : #pragma omp simd
392 : for(unsigned j=0; j<nder; ++j) {
393 : f[j] += src[j];
394 : }
395 : }
396 : }
397 : }
398 :
399 373908 : if(ncp>4*comm.Get_size()) {
400 22970 : comm.Sum(&forcesForApply[0],nder);
401 : }
402 : return true;
403 : }
404 :
405 : }
|
