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 33187 : void ActionWithValue::registerKeywords(Keywords& keys) {
33 33187 : 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 66374 : keys.addFlag("NUMERICAL_DERIVATIVES", false, "calculate the derivatives for these quantities numerically");
39 66374 : keys.add("hidden","HAS_VALUES","this is used in json output to determine those actions that have values");
40 33187 : }
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 882 : void ActionWithValue::useCustomisableComponents(Keywords& keys) {
48 1764 : if( !keys.outputComponentExists(".#!custom") ) {
49 1764 : 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");
50 : }
51 882 : 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 882 : }
55 :
56 31405 : ActionWithValue::ActionWithValue(const ActionOptions&ao):
57 : Action(ao),
58 31405 : firststep(true),
59 31405 : noderiv(true),
60 31405 : numericalDerivatives(false) {
61 62810 : if( keywords.exists("NUMERICAL_DERIVATIVES") ) {
62 32696 : parseFlag("NUMERICAL_DERIVATIVES",numericalDerivatives);
63 : }
64 62810 : if(!keywords.exists("NO_ACTION_LOG") && numericalDerivatives) {
65 66 : log.printf(" using numerical derivatives\n");
66 : }
67 31405 : }
68 :
69 31405 : ActionWithValue::~ActionWithValue() {
70 : // empty destructor to delete unique_ptr
71 31405 : }
72 :
73 2428875 : void ActionWithValue::clearInputForces( const bool& force ) {
74 5713021 : for(unsigned i=0; i<values.size(); i++) {
75 : values[i]->clearInputForce();
76 : }
77 2428875 : }
78 :
79 1893270 : void ActionWithValue::clearDerivatives( const bool& force ) {
80 1893270 : unsigned nt = OpenMP::getNumThreads();
81 1893270 : #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 1893270 : }
89 :
90 : // -- These are the routine for copying the value pointers to other classes -- //
91 :
92 11841939 : bool ActionWithValue::exists( const std::string& name ) const {
93 107520023 : for(unsigned i=0; i<values.size(); ++i) {
94 95702965 : 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 33372905 : Value* ActionWithValue::copyOutput( const std::string& name ) const {
113 111339166 : for(unsigned i=0; i<values.size(); ++i) {
114 111339166 : if (values[i]->name==name) {
115 33372905 : return values[i].get();
116 : }
117 : }
118 0 : plumed_merror("there is no pointer with name " + name);
119 : }
120 :
121 1152736 : Value* ActionWithValue::copyOutput( const unsigned& n ) const {
122 1152736 : plumed_massert(n<values.size(),"you have requested a pointer that is out of bounds");
123 1152736 : return values[n].get();
124 : }
125 :
126 : // -- HERE WE HAVE THE STUFF FOR THE DEFAULT VALUE -- //
127 :
128 12920 : void ActionWithValue::addValue( const std::vector<unsigned>& shape ) {
129 25840 : if( !keywords.outputComponentExists(".#!value") ) {
130 138 : warning("documentation for the value calculated by this action has not been included");
131 : }
132 12920 : plumed_massert(values.empty(),"You have already added the default value for this action");
133 12920 : values.emplace_back(Tools::make_unique<Value>(this,getLabel(), false, shape ) );
134 12920 : }
135 :
136 5158 : void ActionWithValue::addValueWithDerivatives( const std::vector<unsigned>& shape ) {
137 10316 : if( !keywords.outputComponentExists(".#!value") ) {
138 104 : warning("documentation for the value calculated by this action has not been included");
139 : }
140 5158 : plumed_massert(values.empty(),"You have already added the default value for this action");
141 5158 : values.emplace_back(Tools::make_unique<Value>(this,getLabel(), true, shape ) );
142 5158 : }
143 :
144 17252 : void ActionWithValue::setNotPeriodic() {
145 17252 : plumed_massert(values.size()==1,"The number of components is not equal to one");
146 17252 : plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default");
147 17252 : values[0]->min=0;
148 17252 : values[0]->max=0;
149 17252 : values[0]->setupPeriodicity();
150 17252 : }
151 :
152 825 : void ActionWithValue::setPeriodic( const std::string& min, const std::string& max ) {
153 825 : plumed_massert(values.size()==1,"The number of components is not equal to one");
154 825 : plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default");
155 825 : values[0]->setDomain( min, max );
156 825 : }
157 :
158 : // -- HERE WE HAVE THE STUFF FOR NAMED VALUES / COMPONENTS -- //
159 :
160 45827 : void ActionWithValue::addComponent( const std::string& name, const std::vector<unsigned>& shape ) {
161 45827 : if( !keywords.outputComponentExists(name) ) {
162 0 : plumed_merror("a description of component " + name + " has not been added to the manual. Components should be registered like keywords in "
163 : "registerKeywords as described in the developer docs.");
164 : }
165 : std::string thename;
166 91654 : thename=getLabel() + "." + name;
167 18798887 : for(unsigned i=0; i<values.size(); ++i) {
168 18753060 : plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components");
169 18753060 : plumed_massert(values[i]->name!=thename,"there is already a value with this name: "+thename);
170 18753060 : 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"
171 : "Remove the line addComponent(\"bias\") from your bias.");
172 : }
173 45827 : values.emplace_back(Tools::make_unique<Value>(this,thename, false, shape ) );
174 45827 : std::string msg=" added component to this action: "+thename+" \n";
175 45827 : log.printf(msg.c_str());
176 45827 : }
177 :
178 41426 : void ActionWithValue::addComponentWithDerivatives( const std::string& name, const std::vector<unsigned>& shape ) {
179 41426 : if( !keywords.outputComponentExists(name) ) {
180 0 : plumed_merror("a description of component " + name + " has not been added to the manual. Components should be registered like keywords in "
181 : "registerKeywords as described in the developer doc.");
182 : }
183 : std::string thename;
184 82852 : thename=getLabel() + "." + name;
185 2504040 : for(unsigned i=0; i<values.size(); ++i) {
186 2462614 : plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components");
187 2462614 : plumed_massert(values[i]->name!=thename,"there is already a value with this name: "+thename);
188 2462614 : 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"
189 : "Remove the line addComponentWithDerivatives(\"bias\") from your bias.");
190 : }
191 41426 : values.emplace_back(Tools::make_unique<Value>(this,thename, true, shape ) );
192 41426 : std::string msg=" added component to this action: "+thename+" \n";
193 41426 : log.printf(msg.c_str());
194 41426 : }
195 :
196 23 : std::string ActionWithValue::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
197 46 : if( keys.outputComponentExists(".#!custom") ) {
198 0 : return "a quantity calculated by the action " + getName() + " with label " + getLabel();
199 : }
200 23 : std::size_t und=cname.find_last_of("_");
201 23 : std::size_t hyph=cname.find_first_of("-");
202 23 : if( und!=std::string::npos ) {
203 0 : return keys.getOutputComponentDescription(cname.substr(und)) + " This particular component measures this quantity for the input CV named " + cname.substr(0,und);
204 : }
205 23 : if( hyph!=std::string::npos ) {
206 14 : return keys.getOutputComponentDescription(cname.substr(0,hyph)) + " This is the " + cname.substr(hyph+1) + "th of these quantities";
207 : }
208 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" );
209 16 : return keys.getOutputComponentDescription( cname );
210 : }
211 :
212 11816445 : int ActionWithValue::getComponent( const std::string& name ) const {
213 11816445 : plumed_massert( !exists( getLabel() ), "You should not be calling this routine if you are using a value");
214 : std::string thename;
215 23632890 : thename=getLabel() + "." + name;
216 65542053 : for(unsigned i=0; i<values.size(); ++i) {
217 65542053 : if (values[i]->name==thename) {
218 11816445 : return i;
219 : }
220 : }
221 0 : plumed_merror("there is no component with name " + name);
222 : }
223 :
224 0 : std::string ActionWithValue::getComponentsList( ) const {
225 : std::string complist;
226 0 : for(unsigned i=0; i<values.size(); ++i) {
227 0 : complist+=values[i]->name+" ";
228 : }
229 0 : return complist;
230 : }
231 :
232 24075 : std::vector<std::string> ActionWithValue::getComponentsVector( ) const {
233 : std::vector<std::string> complist;
234 349845 : for(unsigned i=0; i<values.size(); ++i) {
235 325770 : complist.push_back(values[i]->name);
236 : }
237 24075 : return complist;
238 0 : }
239 :
240 84033 : void ActionWithValue::componentIsNotPeriodic( const std::string& name ) {
241 84033 : int kk=getComponent(name);
242 84033 : values[kk]->min=0;
243 84033 : values[kk]->max=0;
244 84033 : values[kk]->setupPeriodicity();
245 84033 : }
246 :
247 139 : void ActionWithValue::componentIsPeriodic( const std::string& name, const std::string& min, const std::string& max ) {
248 139 : int kk=getComponent(name);
249 139 : values[kk]->setDomain(min,max);
250 139 : }
251 :
252 1858602 : void ActionWithValue::setGradientsIfNeeded() {
253 3717204 : if(isOptionOn("GRADIENTS")) {
254 201 : ActionAtomistic* aa=castToActionAtomistic();
255 201 : if(aa) {
256 308 : for(unsigned i=0; i<values.size(); i++) {
257 190 : unsigned start=0;
258 : values[i]->gradients.clear();
259 190 : values[i]->setGradients( aa, start );
260 : }
261 : } else {
262 83 : ActionWithArguments* aarg = castToActionWithArguments();
263 83 : if( !aarg ) {
264 0 : plumed_merror( "failing in " + getLabel() );
265 : }
266 166 : for(unsigned i=0; i<values.size(); i++) {
267 83 : unsigned start=0;
268 : values[i]->gradients.clear();
269 83 : aarg->setGradients( values[i].get(), start );
270 : }
271 : }
272 : }
273 1858602 : }
274 :
275 1642126 : void ActionWithValue::turnOnDerivatives() {
276 : // Turn on the derivatives
277 1642126 : noderiv=false;
278 : // Resize the derivatives
279 5997939 : for(unsigned i=0; i<values.size(); ++i) {
280 4355813 : values[i]->resizeDerivatives( getNumberOfDerivatives() );
281 : }
282 : // And turn on the derivatives in all actions on which we are dependent
283 3278154 : for(unsigned i=0; i<getDependencies().size(); ++i) {
284 1636028 : ActionWithValue* vv=getDependencies()[i]->castToActionWithValue();
285 1636028 : if(vv) {
286 1636028 : vv->turnOnDerivatives();
287 : }
288 : }
289 1642126 : }
290 :
291 11732273 : Value* ActionWithValue::getPntrToComponent( const std::string& name ) {
292 11732273 : int kk=getComponent(name);
293 11732273 : return values[kk].get();
294 : }
295 :
296 1199391010 : const Value* ActionWithValue::getConstPntrToComponent(int n) const {
297 : plumed_dbg_massert(n<values.size(),"you have requested a pointer that is out of bounds");
298 1199391010 : return values[n].get();
299 : }
300 :
301 87512997 : Value* ActionWithValue::getPntrToComponent( int n ) {
302 : plumed_dbg_massert(n<values.size(),"you have requested a pointer that is out of bounds");
303 87512997 : return values[n].get();
304 : }
305 :
306 4655603 : bool ActionWithValue::calculateOnUpdate() {
307 4655603 : if( firststep ) {
308 205980 : ActionWithArguments* aa=dynamic_cast<ActionWithArguments*>(this);
309 205980 : if(aa) {
310 184050 : const std::vector<Value*> & args(aa->getArguments());
311 1261770 : for(const auto & p : args ) {
312 1078134 : if( p->calculateOnUpdate() ) {
313 869 : for(unsigned i=0; i<values.size(); ++i) {
314 910 : values[i]->setValType("calcFromAverage");
315 : }
316 : break;
317 : }
318 : }
319 : }
320 205980 : firststep=false;
321 : }
322 11149540 : for(unsigned i=0; i<values.size(); ++i) {
323 6502795 : if( values[i]->calculateOnUpdate() ) {
324 : return true;
325 : }
326 : }
327 : return false;
328 : }
329 :
330 502958 : bool ActionWithValue::checkForForces() {
331 502958 : const unsigned ncp=getNumberOfComponents();
332 502958 : unsigned nder=getNumberOfDerivatives();
333 502958 : if( ncp==0 || nder==0 ) {
334 : return false;
335 : }
336 :
337 : unsigned nvalsWithForce=0;
338 501812 : valsToForce.resize(ncp);
339 1402576 : for(unsigned i=0; i<ncp; ++i) {
340 900764 : if( values[i]->hasForce && !values[i]->isConstant() ) {
341 292738 : valsToForce[nvalsWithForce]=i;
342 292738 : nvalsWithForce++;
343 : }
344 : }
345 501812 : if( nvalsWithForce==0 ) {
346 : return false;
347 : }
348 :
349 : // Make sure forces to apply is empty of forces
350 241999 : if( forcesForApply.size()!=nder ) {
351 3679 : forcesForApply.resize( nder );
352 : }
353 : std::fill(forcesForApply.begin(),forcesForApply.end(),0);
354 :
355 : unsigned stride=1;
356 : unsigned rank=0;
357 241999 : if(ncp>4*comm.Get_size()) {
358 22970 : stride=comm.Get_size();
359 22970 : rank=comm.Get_rank();
360 : }
361 :
362 241999 : unsigned nt=OpenMP::getNumThreads();
363 241999 : if(nt>ncp/(4*stride)) {
364 : nt=1;
365 : }
366 :
367 241999 : #pragma omp parallel num_threads(nt)
368 : {
369 : std::vector<double> omp_f;
370 : if( nt>1 ) {
371 : omp_f.resize(nder,0);
372 : }
373 : #pragma omp for
374 : for(unsigned i=rank; i<nvalsWithForce; i+=stride) {
375 : double ff=values[valsToForce[i]]->inputForce[0];
376 : std::vector<double> & thisderiv( values[valsToForce[i]]->data );
377 : int nn=nder;
378 : int one1=1;
379 : int one2=1;
380 : if( nt>1 ) {
381 : plumed_blas_daxpy(&nn,&ff,thisderiv.data()+1,&one1,omp_f.data(),&one2);
382 : } else {
383 : plumed_blas_daxpy(&nn,&ff,thisderiv.data()+1,&one1,forcesForApply.data(),&one2);
384 : }
385 : // if( nt>1 ) for(unsigned j=0; j<nder; ++j) omp_f[j] += ff*thisderiv[1+j];
386 : //else for(unsigned j=0; j<nder; ++j) forcesForApply[j] += ff*thisderiv[1+j];
387 : }
388 : #pragma omp critical
389 : {
390 : if( nt>1 ) {
391 : int nn=forcesForApply.size();
392 : double one0=1.0;
393 : int one1=1;
394 : int one2=1;
395 : plumed_blas_daxpy(&nn,&one0,omp_f.data(),&one1,forcesForApply.data(),&one2);
396 : }
397 : // for(unsigned j=0; j<forcesForApply.size(); ++j) {
398 : // forcesForApply[j]+=omp_f[j];
399 : // }
400 : }
401 : }
402 :
403 241999 : if(ncp>4*comm.Get_size()) {
404 22970 : comm.Sum(&forcesForApply[0],nder);
405 : }
406 : return true;
407 : }
408 :
409 : }
|
