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 32340 : void ActionWithValue::registerKeywords(Keywords& keys) {
33 32340 : 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 64680 : keys.addFlag("NUMERICAL_DERIVATIVES", false, "calculate the derivatives for these quantities numerically");
39 64680 : keys.add("hidden","HAS_VALUES","this is used in json output to determine those actions that have values");
40 32340 : }
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 2562 : if( !keys.outputComponentExists(".#!custom") ) 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 854 : keys.setComponentsIntroduction("The names of the components in this action can be customized by the user in the "
50 : "actions input file. However, in addition to the components that can be customized the "
51 : "following quantities will always be output");
52 854 : }
53 :
54 31378 : ActionWithValue::ActionWithValue(const ActionOptions&ao):
55 : Action(ao),
56 31378 : firststep(true),
57 31378 : noderiv(true),
58 31378 : numericalDerivatives(false)
59 : {
60 78821 : if( keywords.exists("NUMERICAL_DERIVATIVES") ) parseFlag("NUMERICAL_DERIVATIVES",numericalDerivatives);
61 62756 : if(!keywords.exists("NO_ACTION_LOG") && numericalDerivatives) log.printf(" using numerical derivatives\n");
62 31378 : }
63 :
64 31378 : ActionWithValue::~ActionWithValue() {
65 : // empty destructor to delete unique_ptr
66 31378 : }
67 :
68 2409500 : void ActionWithValue::clearInputForces( const bool& force ) {
69 5665575 : for(unsigned i=0; i<values.size(); i++) values[i]->clearInputForce();
70 2409500 : }
71 :
72 1888715 : void ActionWithValue::clearDerivatives( const bool& force ) {
73 1888715 : unsigned nt = OpenMP::getNumThreads();
74 1888715 : #pragma omp parallel num_threads(nt)
75 : {
76 : #pragma omp for
77 : for(unsigned i=0; i<values.size(); i++) values[i]->clearDerivatives();
78 : }
79 1888715 : }
80 :
81 : // -- These are the routine for copying the value pointers to other classes -- //
82 :
83 10717649 : bool ActionWithValue::exists( const std::string& name ) const {
84 101878095 : for(unsigned i=0; i<values.size(); ++i) {
85 91185159 : if (values[i]->name==name) return true;
86 : }
87 : return false;
88 : }
89 :
90 19 : void ActionWithValue::getMatrixColumnTitles( std::vector<std::string>& argnames ) const {
91 19 : plumed_assert( getNumberOfComponents()==1 && getConstPntrToComponent(0)->getRank()==2 );
92 19 : unsigned nargs = getConstPntrToComponent(0)->getShape()[1]; std::string aname = getConstPntrToComponent(0)->getName();
93 393 : for(unsigned j=0; j<nargs; ++j) { std::string nn; Tools::convert( j+1, nn ); argnames.push_back( aname + "." + nn ); }
94 19 : }
95 :
96 33605533 : Value* ActionWithValue::copyOutput( const std::string& name ) const {
97 112036508 : for(unsigned i=0; i<values.size(); ++i) {
98 112036508 : if (values[i]->name==name) return values[i].get();
99 : }
100 0 : plumed_merror("there is no pointer with name " + name);
101 : }
102 :
103 1143728 : Value* ActionWithValue::copyOutput( const unsigned& n ) const {
104 1143728 : plumed_massert(n<values.size(),"you have requested a pointer that is out of bounds");
105 1143728 : return values[n].get();
106 : }
107 :
108 : // -- HERE WE HAVE THE STUFF FOR THE DEFAULT VALUE -- //
109 :
110 12714 : void ActionWithValue::addValue( const std::vector<unsigned>& shape ) {
111 25428 : if( !keywords.outputComponentExists(".#!value") ) warning("documentation for the value calculated by this action has not been included");
112 25428 : else plumed_massert( keywords.componentHasCorrectType(".#!value",shape.size(),false), "documentation for type of value is incorrect");
113 12714 : plumed_massert(values.empty(),"You have already added the default value for this action");
114 12714 : values.emplace_back(Tools::make_unique<Value>(this,getLabel(), false, shape ) );
115 12714 : }
116 :
117 5142 : void ActionWithValue::addValueWithDerivatives( const std::vector<unsigned>& shape ) {
118 10286 : if( !keywords.outputComponentExists(".#!value") ) warning("documentation for the value calculated by this action has not been included");
119 10280 : else plumed_massert( keywords.componentHasCorrectType(".#!value",shape.size(),true), "documentation for type of value is incorrect");
120 5142 : plumed_massert(values.empty(),"You have already added the default value for this action");
121 5142 : values.emplace_back(Tools::make_unique<Value>(this,getLabel(), true, shape ) );
122 5142 : }
123 :
124 17030 : void ActionWithValue::setNotPeriodic() {
125 17030 : plumed_massert(values.size()==1,"The number of components is not equal to one");
126 17030 : plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default");
127 17030 : values[0]->min=0; values[0]->max=0;
128 17030 : values[0]->setupPeriodicity();
129 17030 : }
130 :
131 825 : void ActionWithValue::setPeriodic( const std::string& min, const std::string& max ) {
132 825 : plumed_massert(values.size()==1,"The number of components is not equal to one");
133 825 : plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default");
134 825 : values[0]->setDomain( min, max );
135 825 : }
136 :
137 : // -- HERE WE HAVE THE STUFF FOR NAMED VALUES / COMPONENTS -- //
138 :
139 46056 : void ActionWithValue::addComponent( const std::string& name, const std::vector<unsigned>& shape ) {
140 46056 : if( !keywords.outputComponentExists(name) ) {
141 0 : plumed_merror("a description of component " + name + " has not been added to the manual. Components should be registered like keywords in "
142 : "registerKeywords as described in the developer docs.");
143 : }
144 46056 : plumed_massert( keywords.componentHasCorrectType(name,shape.size(),false), "documentation for type of component " + name + " is incorrect");
145 92112 : std::string thename; thename=getLabel() + "." + name;
146 18798986 : for(unsigned i=0; i<values.size(); ++i) {
147 18752930 : plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components");
148 18752930 : plumed_massert(values[i]->name!=thename,"there is already a value with this name: "+thename);
149 18752930 : 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"
150 : "Remove the line addComponent(\"bias\") from your bias.");
151 : }
152 46056 : values.emplace_back(Tools::make_unique<Value>(this,thename, false, shape ) );
153 46056 : std::string msg=" added component to this action: "+thename+" \n";
154 46056 : log.printf(msg.c_str());
155 46056 : }
156 :
157 42173 : void ActionWithValue::addComponentWithDerivatives( const std::string& name, const std::vector<unsigned>& shape ) {
158 42173 : if( !keywords.outputComponentExists(name) ) {
159 0 : plumed_merror("a description of component " + name + " has not been added to the manual. Components should be registered like keywords in "
160 : "registerKeywords as described in the developer doc.");
161 : }
162 42173 : plumed_massert( keywords.componentHasCorrectType(name,shape.size(),true), "documentation for type of component " + name + " is incorrect");
163 84346 : std::string thename; thename=getLabel() + "." + name;
164 2505321 : for(unsigned i=0; i<values.size(); ++i) {
165 2463148 : plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components");
166 2463148 : plumed_massert(values[i]->name!=thename,"there is already a value with this name: "+thename);
167 2463148 : 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"
168 : "Remove the line addComponentWithDerivatives(\"bias\") from your bias.");
169 : }
170 42173 : values.emplace_back(Tools::make_unique<Value>(this,thename, true, shape ) );
171 42173 : std::string msg=" added component to this action: "+thename+" \n";
172 42173 : log.printf(msg.c_str());
173 42173 : }
174 :
175 20 : std::string ActionWithValue::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
176 40 : if( keys.outputComponentExists(".#!custom") ) return "a quantity calculated by the action " + getName() + " with label " + getLabel();
177 20 : std::size_t und=cname.find_last_of("_"); std::size_t hyph=cname.find_first_of("-");
178 20 : if( und!=std::string::npos ) return keys.getOutputComponentDescription(cname.substr(und)) + " This particular component measures this quantity for the input CV named " + cname.substr(0,und);
179 27 : if( hyph!=std::string::npos ) return keys.getOutputComponentDescription(cname.substr(0,hyph)) + " This is the " + cname.substr(hyph+1) + "th of these quantities";
180 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" );
181 13 : return keys.getOutputComponentDescription( cname );
182 : }
183 :
184 10692363 : int ActionWithValue::getComponent( const std::string& name ) const {
185 10692363 : plumed_massert( !exists( getLabel() ), "You should not be calling this routine if you are using a value");
186 21384726 : std::string thename; thename=getLabel() + "." + name;
187 64408643 : for(unsigned i=0; i<values.size(); ++i) {
188 64408643 : if (values[i]->name==thename) return i;
189 : }
190 0 : plumed_merror("there is no component with name " + name);
191 : }
192 :
193 0 : std::string ActionWithValue::getComponentsList( ) const {
194 : std::string complist;
195 0 : for(unsigned i=0; i<values.size(); ++i) {
196 0 : complist+=values[i]->name+" ";
197 : }
198 0 : return complist;
199 : }
200 :
201 24075 : std::vector<std::string> ActionWithValue::getComponentsVector( ) const {
202 : std::vector<std::string> complist;
203 349845 : for(unsigned i=0; i<values.size(); ++i) {
204 325770 : complist.push_back(values[i]->name);
205 : }
206 24075 : return complist;
207 0 : }
208 :
209 85009 : void ActionWithValue::componentIsNotPeriodic( const std::string& name ) {
210 85009 : int kk=getComponent(name);
211 85009 : values[kk]->min=0; values[kk]->max=0;
212 85009 : values[kk]->setupPeriodicity();
213 85009 : }
214 :
215 139 : void ActionWithValue::componentIsPeriodic( const std::string& name, const std::string& min, const std::string& max ) {
216 139 : int kk=getComponent(name);
217 139 : values[kk]->setDomain(min,max);
218 139 : }
219 :
220 1854053 : void ActionWithValue::setGradientsIfNeeded() {
221 3708106 : if(isOptionOn("GRADIENTS")) {
222 201 : ActionAtomistic* aa=castToActionAtomistic();
223 201 : if(aa) {
224 308 : for(unsigned i=0; i<values.size(); i++) { unsigned start=0; values[i]->gradients.clear(); values[i]->setGradients( aa, start ); }
225 : } else {
226 83 : ActionWithArguments* aarg = castToActionWithArguments();
227 83 : if( !aarg ) plumed_merror( "failing in " + getLabel() );
228 166 : for(unsigned i=0; i<values.size(); i++) { unsigned start=0; values[i]->gradients.clear(); aarg->setGradients( values[i].get(), start ); }
229 : }
230 : }
231 1854053 : }
232 :
233 1640136 : void ActionWithValue::turnOnDerivatives() {
234 : // Turn on the derivatives
235 1640136 : noderiv=false;
236 : // Resize the derivatives
237 5993745 : for(unsigned i=0; i<values.size(); ++i) values[i]->resizeDerivatives( getNumberOfDerivatives() );
238 : // And turn on the derivatives in all actions on which we are dependent
239 3274319 : for(unsigned i=0; i<getDependencies().size(); ++i) {
240 1634183 : ActionWithValue* vv=getDependencies()[i]->castToActionWithValue();
241 1634183 : if(vv) vv->turnOnDerivatives();
242 : }
243 1640136 : }
244 :
245 10607215 : Value* ActionWithValue::getPntrToComponent( const std::string& name ) {
246 10607215 : int kk=getComponent(name);
247 10607215 : return values[kk].get();
248 : }
249 :
250 1190410193 : const Value* ActionWithValue::getConstPntrToComponent(int n) const {
251 : plumed_dbg_massert(n<values.size(),"you have requested a pointer that is out of bounds");
252 1190410193 : return values[n].get();
253 : }
254 :
255 87336097 : Value* ActionWithValue::getPntrToComponent( int n ) {
256 : plumed_dbg_massert(n<values.size(),"you have requested a pointer that is out of bounds");
257 87336097 : return values[n].get();
258 : }
259 :
260 4616827 : bool ActionWithValue::calculateOnUpdate() {
261 4616827 : if( firststep ) {
262 205953 : ActionWithArguments* aa=dynamic_cast<ActionWithArguments*>(this);
263 205953 : if(aa) {
264 183959 : const std::vector<Value*> & args(aa->getArguments());
265 1261557 : for(const auto & p : args ) {
266 1078012 : if( p->calculateOnUpdate() ) {
267 869 : for(unsigned i=0; i<values.size(); ++i) values[i]->setValType("calcFromAverage");
268 : break;
269 : }
270 : }
271 : }
272 205953 : firststep=false;
273 : }
274 11056154 : for(unsigned i=0; i<values.size(); ++i) {
275 6448185 : if( values[i]->calculateOnUpdate() ) return true;
276 : }
277 : return false;
278 : }
279 :
280 500897 : bool ActionWithValue::checkForForces() {
281 500897 : const unsigned ncp=getNumberOfComponents();
282 500897 : unsigned nder=getNumberOfDerivatives();
283 500897 : if( ncp==0 || nder==0 ) return false;
284 :
285 : unsigned nvalsWithForce=0;
286 499751 : valsToForce.resize(ncp);
287 1398434 : for(unsigned i=0; i<ncp; ++i) {
288 898683 : if( values[i]->hasForce && !values[i]->isConstant() ) {
289 290687 : valsToForce[nvalsWithForce]=i; nvalsWithForce++;
290 : }
291 : }
292 499751 : if( nvalsWithForce==0 ) return false;
293 :
294 : // Make sure forces to apply is empty of forces
295 239948 : if( forcesForApply.size()!=nder ) forcesForApply.resize( nder );
296 : std::fill(forcesForApply.begin(),forcesForApply.end(),0);
297 :
298 : unsigned stride=1;
299 : unsigned rank=0;
300 239948 : if(ncp>4*comm.Get_size()) {
301 22970 : stride=comm.Get_size();
302 22970 : rank=comm.Get_rank();
303 : }
304 :
305 239948 : unsigned nt=OpenMP::getNumThreads();
306 239948 : if(nt>ncp/(4*stride)) nt=1;
307 :
308 239948 : #pragma omp parallel num_threads(nt)
309 : {
310 : std::vector<double> omp_f;
311 : if( nt>1 ) omp_f.resize(nder,0);
312 : #pragma omp for
313 : for(unsigned i=rank; i<nvalsWithForce; i+=stride) {
314 : double ff=values[valsToForce[i]]->inputForce[0];
315 : std::vector<double> & thisderiv( values[valsToForce[i]]->data );
316 : int nn=nder;
317 : int one1=1;
318 : int one2=1;
319 : if( nt>1 ) plumed_blas_daxpy(&nn,&ff,thisderiv.data()+1,&one1,omp_f.data(),&one2);
320 : else plumed_blas_daxpy(&nn,&ff,thisderiv.data()+1,&one1,forcesForApply.data(),&one2);
321 : // if( nt>1 ) for(unsigned j=0; j<nder; ++j) omp_f[j] += ff*thisderiv[1+j];
322 : //else for(unsigned j=0; j<nder; ++j) forcesForApply[j] += ff*thisderiv[1+j];
323 : }
324 : #pragma omp critical
325 : {
326 : if( nt>1 ) {
327 : int nn=forcesForApply.size();
328 : double one0=1.0;
329 : int one1=1;
330 : int one2=1;
331 : plumed_blas_daxpy(&nn,&one0,omp_f.data(),&one1,forcesForApply.data(),&one2);
332 : }
333 : // for(unsigned j=0; j<forcesForApply.size(); ++j) {
334 : // forcesForApply[j]+=omp_f[j];
335 : // }
336 : }
337 : }
338 :
339 239948 : if(ncp>4*comm.Get_size()) comm.Sum(&forcesForApply[0],nder);
340 : return true;
341 : }
342 :
343 : }
|
