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","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 1140751 : Value* ActionWithValue::copyOutput( const unsigned& n ) const {
104 1140751 : plumed_massert(n<values.size(),"you have requested a pointer that is out of bounds");
105 1140751 : 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 25494 : if( !keywords.outputComponentExists(".#!value") ) warning("documentation for the value calculated by this action has not been included");
112 12714 : plumed_massert(values.empty(),"You have already added the default value for this action");
113 12714 : values.emplace_back(Tools::make_unique<Value>(this,getLabel(), false, shape ) );
114 12714 : }
115 :
116 5142 : void ActionWithValue::addValueWithDerivatives( const std::vector<unsigned>& shape ) {
117 10332 : if( !keywords.outputComponentExists(".#!value") ) warning("documentation for the value calculated by this action has not been included");
118 5142 : plumed_massert(values.empty(),"You have already added the default value for this action");
119 5142 : values.emplace_back(Tools::make_unique<Value>(this,getLabel(), true, shape ) );
120 5142 : }
121 :
122 17030 : void ActionWithValue::setNotPeriodic() {
123 17030 : plumed_massert(values.size()==1,"The number of components is not equal to one");
124 17030 : plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default");
125 17030 : values[0]->min=0; values[0]->max=0;
126 17030 : values[0]->setupPeriodicity();
127 17030 : }
128 :
129 825 : void ActionWithValue::setPeriodic( const std::string& min, const std::string& max ) {
130 825 : plumed_massert(values.size()==1,"The number of components is not equal to one");
131 825 : plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default");
132 825 : values[0]->setDomain( min, max );
133 825 : }
134 :
135 : // -- HERE WE HAVE THE STUFF FOR NAMED VALUES / COMPONENTS -- //
136 :
137 46056 : void ActionWithValue::addComponent( const std::string& name, const std::vector<unsigned>& shape ) {
138 46056 : if( !keywords.outputComponentExists(name) ) {
139 0 : plumed_merror("a description of component " + name + " has not been added to the manual. Components should be registered like keywords in "
140 : "registerKeywords as described in the developer docs.");
141 : }
142 92112 : std::string thename; thename=getLabel() + "." + name;
143 18798986 : for(unsigned i=0; i<values.size(); ++i) {
144 18752930 : plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components");
145 18752930 : plumed_massert(values[i]->name!=thename,"there is already a value with this name: "+thename);
146 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"
147 : "Remove the line addComponent(\"bias\") from your bias.");
148 : }
149 46056 : values.emplace_back(Tools::make_unique<Value>(this,thename, false, shape ) );
150 46056 : std::string msg=" added component to this action: "+thename+" \n";
151 46056 : log.printf(msg.c_str());
152 46056 : }
153 :
154 42173 : void ActionWithValue::addComponentWithDerivatives( const std::string& name, const std::vector<unsigned>& shape ) {
155 42173 : if( !keywords.outputComponentExists(name) ) {
156 0 : plumed_merror("a description of component " + name + " has not been added to the manual. Components should be registered like keywords in "
157 : "registerKeywords as described in the developer doc.");
158 : }
159 84346 : std::string thename; thename=getLabel() + "." + name;
160 2505321 : for(unsigned i=0; i<values.size(); ++i) {
161 2463148 : plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components");
162 2463148 : plumed_massert(values[i]->name!=thename,"there is already a value with this name: "+thename);
163 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"
164 : "Remove the line addComponentWithDerivatives(\"bias\") from your bias.");
165 : }
166 42173 : values.emplace_back(Tools::make_unique<Value>(this,thename, true, shape ) );
167 42173 : std::string msg=" added component to this action: "+thename+" \n";
168 42173 : log.printf(msg.c_str());
169 42173 : }
170 :
171 23 : std::string ActionWithValue::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
172 46 : if( keys.outputComponentExists(".#!custom") ) return "a quantity calculated by the action " + getName() + " with label " + getLabel();
173 23 : std::size_t und=cname.find_last_of("_"); std::size_t hyph=cname.find_first_of("-");
174 23 : 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);
175 30 : if( hyph!=std::string::npos ) return keys.getOutputComponentDescription(cname.substr(0,hyph)) + " This is the " + cname.substr(hyph+1) + "th of these quantities";
176 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" );
177 16 : return keys.getOutputComponentDescription( cname );
178 : }
179 :
180 10692363 : int ActionWithValue::getComponent( const std::string& name ) const {
181 10692363 : plumed_massert( !exists( getLabel() ), "You should not be calling this routine if you are using a value");
182 21384726 : std::string thename; thename=getLabel() + "." + name;
183 64408643 : for(unsigned i=0; i<values.size(); ++i) {
184 64408643 : if (values[i]->name==thename) return i;
185 : }
186 0 : plumed_merror("there is no component with name " + name);
187 : }
188 :
189 0 : std::string ActionWithValue::getComponentsList( ) const {
190 : std::string complist;
191 0 : for(unsigned i=0; i<values.size(); ++i) {
192 0 : complist+=values[i]->name+" ";
193 : }
194 0 : return complist;
195 : }
196 :
197 24075 : std::vector<std::string> ActionWithValue::getComponentsVector( ) const {
198 : std::vector<std::string> complist;
199 349845 : for(unsigned i=0; i<values.size(); ++i) {
200 325770 : complist.push_back(values[i]->name);
201 : }
202 24075 : return complist;
203 0 : }
204 :
205 85009 : void ActionWithValue::componentIsNotPeriodic( const std::string& name ) {
206 85009 : int kk=getComponent(name);
207 85009 : values[kk]->min=0; values[kk]->max=0;
208 85009 : values[kk]->setupPeriodicity();
209 85009 : }
210 :
211 139 : void ActionWithValue::componentIsPeriodic( const std::string& name, const std::string& min, const std::string& max ) {
212 139 : int kk=getComponent(name);
213 139 : values[kk]->setDomain(min,max);
214 139 : }
215 :
216 1854053 : void ActionWithValue::setGradientsIfNeeded() {
217 3708106 : if(isOptionOn("GRADIENTS")) {
218 201 : ActionAtomistic* aa=castToActionAtomistic();
219 201 : if(aa) {
220 308 : for(unsigned i=0; i<values.size(); i++) { unsigned start=0; values[i]->gradients.clear(); values[i]->setGradients( aa, start ); }
221 : } else {
222 83 : ActionWithArguments* aarg = castToActionWithArguments();
223 83 : if( !aarg ) plumed_merror( "failing in " + getLabel() );
224 166 : for(unsigned i=0; i<values.size(); i++) { unsigned start=0; values[i]->gradients.clear(); aarg->setGradients( values[i].get(), start ); }
225 : }
226 : }
227 1854053 : }
228 :
229 1640136 : void ActionWithValue::turnOnDerivatives() {
230 : // Turn on the derivatives
231 1640136 : noderiv=false;
232 : // Resize the derivatives
233 5993745 : for(unsigned i=0; i<values.size(); ++i) values[i]->resizeDerivatives( getNumberOfDerivatives() );
234 : // And turn on the derivatives in all actions on which we are dependent
235 3274319 : for(unsigned i=0; i<getDependencies().size(); ++i) {
236 1634183 : ActionWithValue* vv=getDependencies()[i]->castToActionWithValue();
237 1634183 : if(vv) vv->turnOnDerivatives();
238 : }
239 1640136 : }
240 :
241 10607215 : Value* ActionWithValue::getPntrToComponent( const std::string& name ) {
242 10607215 : int kk=getComponent(name);
243 10607215 : return values[kk].get();
244 : }
245 :
246 1190410193 : const Value* ActionWithValue::getConstPntrToComponent(int n) const {
247 : plumed_dbg_massert(n<values.size(),"you have requested a pointer that is out of bounds");
248 1190410193 : return values[n].get();
249 : }
250 :
251 87336097 : Value* ActionWithValue::getPntrToComponent( int n ) {
252 : plumed_dbg_massert(n<values.size(),"you have requested a pointer that is out of bounds");
253 87336097 : return values[n].get();
254 : }
255 :
256 4616827 : bool ActionWithValue::calculateOnUpdate() {
257 4616827 : if( firststep ) {
258 205953 : ActionWithArguments* aa=dynamic_cast<ActionWithArguments*>(this);
259 205953 : if(aa) {
260 183959 : const std::vector<Value*> & args(aa->getArguments());
261 1261557 : for(const auto & p : args ) {
262 1078012 : if( p->calculateOnUpdate() ) {
263 869 : for(unsigned i=0; i<values.size(); ++i) values[i]->setValType("calcFromAverage");
264 : break;
265 : }
266 : }
267 : }
268 205953 : firststep=false;
269 : }
270 11056154 : for(unsigned i=0; i<values.size(); ++i) {
271 6448185 : if( values[i]->calculateOnUpdate() ) return true;
272 : }
273 : return false;
274 : }
275 :
276 500897 : bool ActionWithValue::checkForForces() {
277 500897 : const unsigned ncp=getNumberOfComponents();
278 500897 : unsigned nder=getNumberOfDerivatives();
279 500897 : if( ncp==0 || nder==0 ) return false;
280 :
281 : unsigned nvalsWithForce=0;
282 499751 : valsToForce.resize(ncp);
283 1398434 : for(unsigned i=0; i<ncp; ++i) {
284 898683 : if( values[i]->hasForce && !values[i]->isConstant() ) {
285 290687 : valsToForce[nvalsWithForce]=i; nvalsWithForce++;
286 : }
287 : }
288 499751 : if( nvalsWithForce==0 ) return false;
289 :
290 : // Make sure forces to apply is empty of forces
291 239948 : if( forcesForApply.size()!=nder ) forcesForApply.resize( nder );
292 : std::fill(forcesForApply.begin(),forcesForApply.end(),0);
293 :
294 : unsigned stride=1;
295 : unsigned rank=0;
296 239948 : if(ncp>4*comm.Get_size()) {
297 22970 : stride=comm.Get_size();
298 22970 : rank=comm.Get_rank();
299 : }
300 :
301 239948 : unsigned nt=OpenMP::getNumThreads();
302 239948 : if(nt>ncp/(4*stride)) nt=1;
303 :
304 239948 : #pragma omp parallel num_threads(nt)
305 : {
306 : std::vector<double> omp_f;
307 : if( nt>1 ) omp_f.resize(nder,0);
308 : #pragma omp for
309 : for(unsigned i=rank; i<nvalsWithForce; i+=stride) {
310 : double ff=values[valsToForce[i]]->inputForce[0];
311 : std::vector<double> & thisderiv( values[valsToForce[i]]->data );
312 : int nn=nder;
313 : int one1=1;
314 : int one2=1;
315 : if( nt>1 ) plumed_blas_daxpy(&nn,&ff,thisderiv.data()+1,&one1,omp_f.data(),&one2);
316 : else plumed_blas_daxpy(&nn,&ff,thisderiv.data()+1,&one1,forcesForApply.data(),&one2);
317 : // if( nt>1 ) for(unsigned j=0; j<nder; ++j) omp_f[j] += ff*thisderiv[1+j];
318 : //else for(unsigned j=0; j<nder; ++j) forcesForApply[j] += ff*thisderiv[1+j];
319 : }
320 : #pragma omp critical
321 : {
322 : if( nt>1 ) {
323 : int nn=forcesForApply.size();
324 : double one0=1.0;
325 : int one1=1;
326 : int one2=1;
327 : plumed_blas_daxpy(&nn,&one0,omp_f.data(),&one1,forcesForApply.data(),&one2);
328 : }
329 : // for(unsigned j=0; j<forcesForApply.size(); ++j) {
330 : // forcesForApply[j]+=omp_f[j];
331 : // }
332 : }
333 : }
334 :
335 239948 : if(ncp>4*comm.Get_size()) comm.Sum(&forcesForApply[0],nder);
336 : return true;
337 : }
338 :
339 : }
|
