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 "Value.h"
23 : #include "ActionWithValue.h"
24 : #include "ActionAtomistic.h"
25 : #include "ActionWithArguments.h"
26 : #include "ActionWithVector.h"
27 : #include "ActionWithVirtualAtom.h"
28 : #include "tools/Exception.h"
29 : #include "tools/OpenMP.h"
30 : #include "tools/OFile.h"
31 : #include "PlumedMain.h"
32 :
33 : namespace PLMD {
34 :
35 164 : Value::Value():
36 164 : action(NULL),
37 164 : value_set(false),
38 164 : hasForce(false),
39 164 : storedata(false),
40 : shape(std::vector<unsigned>()),
41 164 : hasDeriv(true),
42 164 : bufstart(0),
43 164 : streampos(0),
44 164 : matpos(0),
45 164 : ngrid_der(0),
46 164 : ncols(0),
47 164 : book_start(0),
48 164 : symmetric(false),
49 164 : periodicity(unset),
50 164 : min(0.0),
51 164 : max(0.0),
52 164 : max_minus_min(0.0),
53 164 : inv_max_minus_min(0.0),
54 164 : derivativeIsZeroWhenValueIsZero(false) {
55 164 : data.resize(1);
56 164 : inputForce.resize(1);
57 164 : }
58 :
59 18 : Value::Value(const std::string& name):
60 18 : action(NULL),
61 18 : value_set(false),
62 18 : hasForce(false),
63 18 : name(name),
64 18 : storedata(false),
65 : shape(std::vector<unsigned>()),
66 18 : hasDeriv(true),
67 18 : bufstart(0),
68 18 : streampos(0),
69 18 : ngrid_der(0),
70 18 : matpos(0),
71 18 : ncols(0),
72 18 : book_start(0),
73 18 : symmetric(false),
74 18 : periodicity(unset),
75 18 : min(0.0),
76 18 : max(0.0),
77 18 : max_minus_min(0.0),
78 18 : inv_max_minus_min(0.0),
79 18 : derivativeIsZeroWhenValueIsZero(false) {
80 18 : data.resize(1);
81 18 : inputForce.resize(1);
82 18 : data[0]=inputForce[0]=0;
83 18 : }
84 :
85 117896 : Value::Value(ActionWithValue* av, const std::string& name, const bool withderiv, const std::vector<unsigned>&ss):
86 117896 : action(av),
87 117896 : value_set(false),
88 117896 : hasForce(false),
89 117896 : name(name),
90 117896 : storedata(false),
91 117896 : hasDeriv(withderiv),
92 117896 : bufstart(0),
93 117896 : streampos(0),
94 117896 : ngrid_der(0),
95 117896 : matpos(0),
96 117896 : ncols(0),
97 117896 : book_start(0),
98 117896 : symmetric(false),
99 117896 : periodicity(unset),
100 117896 : min(0.0),
101 117896 : max(0.0),
102 117896 : max_minus_min(0.0),
103 117896 : inv_max_minus_min(0.0),
104 117896 : derivativeIsZeroWhenValueIsZero(false) {
105 117896 : if( action ) {
106 117440 : if( action->getName()=="ACCUMULATE" || action->getName()=="COLLECT" ) {
107 157 : valtype=average;
108 : }
109 : }
110 117896 : if( action ) {
111 227514 : storedata=action->getName()=="PUT" || valtype==average;
112 : }
113 117896 : if( ss.size() && withderiv ) {
114 716 : storedata=true;
115 : }
116 117896 : setShape( ss );
117 117896 : }
118 :
119 446 : void Value::setValType( const std::string& vtype ) {
120 446 : if( vtype=="normal" ) {
121 0 : valtype=normal;
122 446 : } else if( vtype=="constant" ) {
123 0 : valtype=constant;
124 446 : } else if( vtype=="average" ) {
125 0 : valtype=average;
126 446 : } else if( vtype=="calcFromAverage" ) {
127 446 : valtype=calcFromAverage;
128 : } else {
129 0 : plumed_merror("invalid valtype " + vtype );
130 : }
131 446 : }
132 :
133 130302 : void Value::setShape( const std::vector<unsigned>&ss ) {
134 : std::size_t tot=1;
135 130302 : shape.resize( ss.size() );
136 161355 : for(unsigned i=0; i<shape.size(); ++i) {
137 31053 : tot = tot*ss[i];
138 31053 : shape[i]=ss[i];
139 : }
140 :
141 130302 : if( shape.size()>0 && hasDeriv ) {
142 : // This is for grids
143 2005 : ngrid_der = shape.size();
144 2005 : if( action ) {
145 2004 : ngrid_der = action->getNumberOfDerivatives();
146 : }
147 2005 : std::size_t ndata = tot*(1+ngrid_der);
148 2005 : data.resize( ndata );
149 2005 : inputForce.resize( tot );
150 128297 : } else if( shape.size()==0 ) {
151 : // This is for scalars
152 105569 : data.resize(1);
153 105569 : inputForce.resize(1);
154 22728 : } else if( storedata && shape.size()<2 ) {
155 : // This is for vectors (matrices have special version because we have sparse storage)
156 12566 : data.resize( tot );
157 12566 : inputForce.resize( tot );
158 : }
159 130302 : }
160 :
161 108980 : void Value::setupPeriodicity() {
162 108980 : if( min==0 && max==0 ) {
163 101189 : periodicity=notperiodic;
164 : } else {
165 7791 : periodicity=periodic;
166 7791 : max_minus_min=max-min;
167 7791 : plumed_massert(max_minus_min>0, "your function has a very strange domain?");
168 7791 : inv_max_minus_min=1.0/max_minus_min;
169 : }
170 108980 : }
171 :
172 1694017 : bool Value::isPeriodic()const {
173 1694017 : plumed_massert(periodicity!=unset,"periodicity should be set");
174 1694017 : return periodicity==periodic;
175 : }
176 :
177 422407 : bool Value::applyForce(std::vector<double>& forces ) const {
178 422407 : if( !hasForce || valtype!=normal ) {
179 : return false;
180 : }
181 : plumed_dbg_massert( data.size()-1==forces.size()," forces array has wrong size" );
182 9645 : const unsigned N=data.size()-1;
183 41732577 : for(unsigned i=0; i<N; ++i) {
184 41722932 : forces[i]=inputForce[0]*data[1+i];
185 : }
186 : return true;
187 : }
188 :
189 1103162 : void Value::setNotPeriodic() {
190 1103162 : min=0;
191 1103162 : max=0;
192 1103162 : periodicity=notperiodic;
193 1103162 : }
194 :
195 7791 : void Value::setDomain(const std::string& pmin,const std::string& pmax) {
196 7791 : str_min=pmin;
197 7791 : if( !Tools::convertNoexcept(str_min,min) ) {
198 0 : action->error("could not convert period string " + str_min + " to real");
199 : }
200 7791 : str_max=pmax;
201 7791 : if( !Tools::convertNoexcept(str_max,max) ) {
202 0 : action->error("could not convert period string " + str_max + " to read");
203 : }
204 7791 : setupPeriodicity();
205 7791 : }
206 :
207 16783 : void Value::getDomain(std::string&minout,std::string&maxout) const {
208 16783 : plumed_massert(periodicity==periodic,"function should be periodic");
209 16783 : minout=str_min;
210 16783 : maxout=str_max;
211 16783 : }
212 :
213 34 : void Value::getDomain(double&minout,double&maxout) const {
214 34 : plumed_massert(periodicity==periodic,"function should be periodic");
215 34 : minout=min;
216 34 : maxout=max;
217 34 : }
218 :
219 190 : void Value::setGradients( ActionAtomistic* aa, unsigned& start ) {
220 : // Can't do gradients if we don't have derivatives
221 190 : if( !hasDeriv ) {
222 : return;
223 : }
224 154 : plumed_assert( shape.size()==0 );
225 8362 : for(unsigned j=0; j<aa->getNumberOfAtoms(); ++j) {
226 8208 : Vector der(data[1+start+3*j],data[1+start+3*j+1],data[1+start+3*j+2]);
227 8208 : aa->getGradient( j, der, gradients );
228 : }
229 154 : start += aa->getNumberOfAtoms();
230 : }
231 :
232 283 : void Value::passGradients( const double& der, std::map<AtomNumber,Vector>& g ) const {
233 8921 : for(const auto & p : gradients) {
234 8638 : AtomNumber iatom=p.first;
235 8638 : g[iatom] += p.second*der;
236 : }
237 283 : }
238 :
239 261 : double Value::projection(const Value& v1,const Value&v2) {
240 : double proj=0.0;
241 : const std::map<AtomNumber,Vector> & grad1(v1.gradients);
242 : const std::map<AtomNumber,Vector> & grad2(v2.gradients);
243 16167 : for(const auto & p1 : grad1) {
244 15906 : AtomNumber a=p1.first;
245 : const auto p2=grad2.find(a);
246 15906 : if(p2!=grad2.end()) {
247 8224 : proj+=dotProduct(p1.second,(*p2).second);
248 : }
249 : }
250 261 : return proj;
251 : }
252 :
253 18740351 : ActionWithValue* Value::getPntrToAction() {
254 18740351 : plumed_assert( action!=NULL );
255 18740351 : return action;
256 : }
257 :
258 4948751 : void Value::set(const std::size_t& n, const double& v ) {
259 4948751 : value_set=true;
260 4948751 : if( getRank()==0 ) {
261 25675 : plumed_assert( n==0 );
262 25675 : data[n]=v;
263 25675 : applyPeriodicity(n);
264 4923076 : } else if( !hasDeriv ) {
265 : plumed_dbg_massert( n<data.size(), "failing in " + getName() );
266 4605760 : data[n]=v;
267 4605760 : applyPeriodicity(n);
268 : } else {
269 317316 : data[n*(1+ngrid_der)] = v;
270 : }
271 4948751 : }
272 :
273 74915 : void Value::push_back( const double& v ) {
274 74915 : value_set=true;
275 74915 : if( shape.size()==1 ) {
276 36356 : data.push_back(v);
277 36356 : shape[0]++;
278 38559 : } else if( shape.size()==2 ) {
279 38559 : data.push_back(v);
280 38559 : shape[0] = std::ceil( data.size() / shape[1] );
281 : }
282 74915 : }
283 :
284 3223564 : std::size_t Value::getIndexInStore( const std::size_t& ival ) const {
285 3223564 : if( shape.size()==2 && ncols<shape[1] ) {
286 2712800 : unsigned irow = std::floor( ival / shape[1] ), jcol = ival%shape[1];
287 11801200 : for(unsigned i=0; i<getRowLength(irow); ++i) {
288 11801200 : if( getRowIndex(irow,i)==jcol ) {
289 2712800 : return irow*ncols+i;
290 : }
291 : }
292 0 : plumed_merror("cannot get store index");
293 : }
294 510764 : return ival;
295 : }
296 :
297 557678237 : double Value::get(const std::size_t& ival, const bool trueind) const {
298 557678237 : if( hasDeriv ) {
299 14924522 : return data[ival*(1+ngrid_der)];
300 : }
301 : #ifdef DNDEBUG
302 : if( action ) {
303 : plumed_dbg_massert( ival<getNumberOfValues(), "could not get value from " + name );
304 : }
305 : #endif
306 542753715 : if( shape.size()==2 && ncols<shape[1] && trueind ) {
307 2000962 : unsigned irow = std::floor( ival / shape[1] ), jcol = ival%shape[1];
308 : // This is a special treatment for the lower triangular matrices that are used when
309 : // we do ITRE with COLLECT_FRAMES
310 2000962 : if( ncols==0 ) {
311 0 : if( jcol<=irow ) {
312 0 : return data[0.5*irow*(irow+1) + jcol];
313 : }
314 : return 0;
315 : }
316 20655470 : for(unsigned i=0; i<getRowLength(irow); ++i) {
317 19057356 : if( getRowIndex(irow,i)==jcol ) {
318 402848 : return data[irow*ncols+i];
319 : }
320 : }
321 : return 0.0;
322 : }
323 540752753 : plumed_massert( ival<data.size(), "cannot get value from " + name );
324 540752753 : return data[ival];
325 : }
326 :
327 39324858 : void Value::addForce(const std::size_t& iforce, double f, const bool trueind) {
328 39324858 : hasForce=true;
329 39324858 : if( shape.size()==2 && !hasDeriv && ncols<shape[1] && trueind ) {
330 0 : unsigned irow = std::floor( iforce / shape[0] ), jcol = iforce%shape[0];
331 0 : for(unsigned i=0; i<getRowLength(irow); ++i) {
332 0 : if( getRowIndex(irow,i)==jcol ) {
333 0 : inputForce[irow*ncols+i]+=f;
334 : return;
335 : }
336 : }
337 0 : plumed_assert( fabs(f)<epsilon );
338 : return;
339 : }
340 39324858 : plumed_massert( iforce<inputForce.size(), "can't add force to " + name );
341 39324858 : inputForce[iforce]+=f;
342 : }
343 :
344 :
345 12206 : void Value::buildDataStore( const bool forprint ) {
346 12206 : if( getRank()==0 ) {
347 : return;
348 : }
349 5310 : storedata=true;
350 5310 : setShape( shape );
351 5310 : if( !forprint ) {
352 : return ;
353 : }
354 192 : ActionWithVector* av=dynamic_cast<ActionWithVector*>( action );
355 192 : if( av ) {
356 166 : (av->getFirstActionInChain())->never_reduce_tasks=true;
357 : }
358 : }
359 :
360 7102 : void Value::reshapeMatrixStore( const unsigned& n ) {
361 : plumed_dbg_assert( shape.size()==2 && !hasDeriv );
362 7102 : if( !storedata ) {
363 : return ;
364 : }
365 7102 : ncols=n;
366 7102 : if( ncols>shape[1] ) {
367 750 : ncols=shape[1];
368 : }
369 7102 : unsigned size=shape[0]*ncols;
370 7102 : if( matrix_bookeeping.size()!=(size+shape[0]) ) {
371 2584 : data.resize( size );
372 2584 : inputForce.resize( size );
373 2584 : matrix_bookeeping.resize( size + shape[0], 0 );
374 2584 : if( ncols>=shape[1] ) {
375 248634 : for(unsigned i=0; i<shape[0]; ++i) {
376 246113 : matrix_bookeeping[(1+ncols)*i] = shape[1];
377 23653381 : for(unsigned j=0; j<shape[1]; ++j) {
378 23407268 : matrix_bookeeping[(1+ncols)*i+1+j]=j;
379 : }
380 : }
381 : }
382 : }
383 7102 : if( ncols<shape[1] ) {
384 : std::fill(matrix_bookeeping.begin(), matrix_bookeeping.end(), 0);
385 : }
386 : }
387 :
388 35389 : void Value::setPositionInMatrixStash( const unsigned& p ) {
389 : plumed_dbg_assert( shape.size()==2 && !hasDeriv );
390 35389 : matpos=p;
391 35389 : }
392 :
393 15992448 : bool Value::ignoreStoredValue(const std::string& c) const {
394 15992448 : if( !storedata && shape.size()>0 ) {
395 : return true;
396 : }
397 6832530 : ActionWithVector* av=dynamic_cast<ActionWithVector*>(action);
398 6832530 : if( av ) {
399 6075155 : return (av->getFirstActionInChain())->getLabel()==c;
400 : }
401 : return false;
402 : }
403 :
404 5573 : void Value::setConstant() {
405 5573 : valtype=constant;
406 5573 : storedata=true;
407 5573 : setShape( shape );
408 5573 : if( getRank()==2 && !hasDeriv ) {
409 125 : reshapeMatrixStore( shape[1] );
410 : }
411 5573 : }
412 :
413 456 : void Value::writeBinary(std::ostream&o) const {
414 456 : o.write(reinterpret_cast<const char*>(&data[0]),data.size()*sizeof(double));
415 456 : }
416 :
417 1274 : void Value::setSymmetric( const bool& sym ) {
418 1274 : plumed_assert( shape.size()==2 && !hasDeriv );
419 1274 : if( sym && shape[0]!=shape[1] ) {
420 0 : plumed_merror("non-square matrix cannot be symmetric");
421 : }
422 1274 : symmetric=sym;
423 1274 : }
424 :
425 2220 : void Value::retrieveEdgeList( unsigned& nedge, std::vector<std::pair<unsigned,unsigned> >& active, std::vector<double>& elems ) {
426 2220 : nedge=0;
427 : plumed_dbg_assert( shape.size()==2 && !hasDeriv );
428 : // Check we have enough space to store the edge list
429 2220 : if( elems.size()<shape[0]*ncols ) {
430 2120 : elems.resize( shape[0]*ncols );
431 2120 : active.resize( shape[0]*ncols );
432 : }
433 :
434 269579 : for(unsigned i=0; i<shape[0]; ++i) {
435 : unsigned ncol = getRowLength(i);
436 15002116 : for(unsigned j=0; j<ncol; ++j) {
437 14734757 : if( fabs(get(i*ncols+j,false))<epsilon ) {
438 11469720 : continue;
439 : }
440 3265037 : if( symmetric && getRowIndex(i,j)>i ) {
441 51540 : continue;
442 : }
443 3213497 : active[nedge].first = i;
444 3213497 : active[nedge].second = getRowIndex(i,j);
445 3213497 : elems[nedge] = get(i*ncols+j,false);
446 3213497 : nedge++;
447 : }
448 : }
449 2220 : }
450 :
451 456 : void Value::readBinary(std::istream&i) {
452 456 : i.read(reinterpret_cast<char*>(&data[0]),data.size()*sizeof(double));
453 456 : }
454 :
455 775666 : void Value::convertIndexToindices(const std::size_t& index, std::vector<unsigned>& indices ) const {
456 775666 : if( hasDeriv || getRank()==1 ) {
457 12302 : std::size_t kk=index;
458 12302 : indices[0]=index%shape[0];
459 12302 : for(unsigned i=1; i<shape.size()-1; ++i) {
460 0 : kk=(kk-indices[i-1])/shape[i-1];
461 0 : indices[i]=kk%shape[i];
462 : }
463 12302 : if(shape.size()>=2) {
464 0 : indices[shape.size()-1]=(kk-indices[shape.size()-2])/shape[shape.size()-2];
465 : }
466 763364 : } else if( getRank()==2 ) {
467 763364 : indices[0]=std::floor( index/shape[1] );
468 763364 : indices[1] = index%shape[1];
469 : }
470 775666 : }
471 :
472 554900 : void Value::print( OFile& ofile ) const {
473 554900 : if( isPeriodic() ) {
474 11537 : ofile.printField( "min_" + name, str_min );
475 23074 : ofile.printField("max_" + name, str_max );
476 : }
477 554900 : if( shape.size()==0 || getNumberOfValues()==1 ) {
478 553581 : ofile.printField( name, get(0) );
479 : } else {
480 1319 : std::vector<unsigned> indices( shape.size() );
481 776985 : for(unsigned i=0; i<getNumberOfValues(); ++i) {
482 775666 : convertIndexToindices( i, indices );
483 775666 : std::string num, fname = name;
484 2314696 : for(unsigned i=0; i<shape.size(); ++i) {
485 1539030 : Tools::convert( indices[i]+1, num );
486 3078060 : fname += "." + num;
487 : }
488 775666 : ofile.printField( fname, get(i) );
489 : }
490 : }
491 554900 : }
492 :
493 240958 : unsigned Value::getGoodNumThreads( const unsigned& j, const unsigned& k ) const {
494 240958 : return OpenMP::getGoodNumThreads( &data[j], (k-j) );
495 : }
496 :
497 7580559 : bool Value::calculateOnUpdate() const {
498 7580559 : return (valtype==average || valtype==calcFromAverage);
499 : }
500 :
501 129 : std::string Value::getValueType() const {
502 129 : if( getRank()==0 ) {
503 33 : return "scalar";
504 : }
505 96 : if( getRank()>0 && hasDerivatives() ) {
506 5 : return "grid";
507 : }
508 91 : if( getRank()==1 ) {
509 70 : return "vector";
510 : }
511 21 : if( getRank()==2 ) {
512 21 : return "matrix";
513 : }
514 0 : plumed_merror("unknown type for value " + getName() );
515 : return "";
516 : }
517 :
518 : }
|
