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 <cstring>
24 : #include "ActionWithValue.h"
25 : #include "ActionAtomistic.h"
26 : #include "tools/Exception.h"
27 : #include "tools/OpenMP.h"
28 : #include "tools/OFile.h"
29 : #include "tools/Communicator.h"
30 : #include "PlumedMain.h"
31 :
32 : namespace PLMD {
33 :
34 2 : Value::Value():
35 4 : data(1,0.0),
36 2 : inputForce(1,0.0),
37 4 : shape() {}
38 :
39 18 : Value::Value(const std::string& valname):
40 36 : data(1,0.0),
41 18 : inputForce(1,0.0),
42 18 : name(valname),
43 36 : shape() {
44 18 : }
45 :
46 135780 : Value::Value(ActionWithValue* av,
47 : const std::string& valname,
48 : const bool withderiv,
49 135780 : const std::vector<std::size_t>&ss):
50 135780 : action(av),
51 135780 : name(valname),
52 135780 : hasDeriv(withderiv) {
53 135780 : if( action ) {
54 135324 : if( action->getName()=="ACCUMULATE" || action->getName()=="COLLECT" ) {
55 193 : valtype=average;
56 : }
57 : }
58 135780 : setShape( ss );
59 135780 : }
60 :
61 468 : void Value::setValType( const std::string& vtype ) {
62 468 : if( vtype=="normal" ) {
63 0 : valtype=normal;
64 468 : } else if( vtype=="constant" ) {
65 0 : valtype=constant;
66 468 : } else if( vtype=="average" ) {
67 0 : valtype=average;
68 468 : } else if( vtype=="calcFromAverage" ) {
69 468 : valtype=calcFromAverage;
70 : } else {
71 0 : plumed_merror("invalid valtype " + vtype );
72 : }
73 468 : }
74 :
75 144381 : void Value::setShape( const std::vector<std::size_t>&ss ) {
76 : std::size_t tot=1;
77 144381 : shape.resize( ss.size() );
78 176035 : for(unsigned i=0; i<shape.size(); ++i) {
79 31654 : tot = tot*ss[i];
80 31654 : shape[i]=ss[i];
81 : }
82 :
83 144381 : if( shape.size()>0 && hasDeriv ) {
84 : // This is for grids
85 660 : ngrid_der = shape.size();
86 660 : if( action ) {
87 659 : ngrid_der = action->getNumberOfDerivatives();
88 : }
89 660 : std::size_t ndata = tot*(1+ngrid_der);
90 660 : data.resize( ndata );
91 660 : inputForce.resize( tot );
92 143721 : } else if( shape.size()==0 ) {
93 : // This is for scalars
94 117885 : data.resize(1);
95 117885 : inputForce.resize(1);
96 25836 : } else if( shape.size()<2 ) {
97 : // This is for vectors (matrices have special version because we have sparse storage)
98 20954 : data.resize( tot );
99 20954 : inputForce.resize( tot );
100 : }
101 144381 : }
102 :
103 124386 : void Value::setupPeriodicity() {
104 124386 : if( min==0 && max==0 ) {
105 113258 : periodicity=notperiodic;
106 : } else {
107 11128 : periodicity=periodic;
108 11128 : max_minus_min=max-min;
109 11128 : plumed_massert(max_minus_min>0, "your function has a very strange domain?");
110 11128 : inv_max_minus_min=1.0/max_minus_min;
111 : }
112 124386 : }
113 :
114 2064398 : bool Value::isPeriodic()const {
115 2064398 : plumed_massert(periodicity!=unset,"periodicity should be set");
116 2064398 : return periodicity==periodic;
117 : }
118 :
119 863845 : bool Value::applyForce(std::vector<double>& forces ) const {
120 863845 : if( !hasForce || valtype!=normal ) {
121 : return false;
122 : }
123 : plumed_dbg_massert( data.size()-1==forces.size()," forces array has wrong size" );
124 12164 : const unsigned N=data.size()-1;
125 42530163 : for(unsigned i=0; i<N; ++i) {
126 42517999 : forces[i]=inputForce[0]*data[1+i];
127 : }
128 : return true;
129 : }
130 :
131 1109725 : void Value::setNotPeriodic() {
132 1109725 : min=0;
133 1109725 : max=0;
134 1109725 : periodicity=notperiodic;
135 1109725 : }
136 :
137 11128 : void Value::setDomain(const std::string& pmin,const std::string& pmax) {
138 11128 : str_min=pmin;
139 11128 : if( !Tools::convertNoexcept(str_min,min) ) {
140 0 : action->error("could not convert period string " + str_min + " to real");
141 : }
142 11128 : str_max=pmax;
143 11128 : if( !Tools::convertNoexcept(str_max,max) ) {
144 0 : action->error("could not convert period string " + str_max + " to read");
145 : }
146 11128 : setupPeriodicity();
147 11128 : }
148 :
149 37068 : void Value::getDomain(std::string&minout,std::string&maxout) const {
150 37068 : plumed_massert(periodicity==periodic,"function should be periodic");
151 37068 : minout=str_min;
152 37068 : maxout=str_max;
153 37068 : }
154 :
155 72 : void Value::getDomain(double&minout,double&maxout) const {
156 72 : plumed_massert(periodicity==periodic,"function should be periodic");
157 72 : minout=min;
158 72 : maxout=max;
159 72 : }
160 :
161 1395 : void Value::setGradients( ActionAtomistic* aa, unsigned& start ) {
162 : // Can't do gradients if we don't have derivatives
163 1395 : if( !hasDeriv ) {
164 : return;
165 : }
166 1359 : plumed_assert( shape.size()==0 );
167 16361 : for(unsigned j=0; j<aa->getNumberOfAtoms(); ++j) {
168 15002 : Vector der(data[1+start+3*j],data[1+start+3*j+1],data[1+start+3*j+2]);
169 15002 : aa->getGradient( j, der, gradients );
170 : }
171 1359 : start += aa->getNumberOfAtoms();
172 : }
173 :
174 2250 : void Value::passGradients( const double& der, std::map<AtomNumber,Vector>& g ) const {
175 18538 : for(const auto & p : gradients) {
176 16288 : AtomNumber iatom=p.first;
177 16288 : g[iatom] += p.second*der;
178 : }
179 2250 : }
180 :
181 3048 : double Value::projection(const Value& v1,const Value&v2) {
182 : double proj=0.0;
183 : const std::map<AtomNumber,Vector> & grad1(v1.gradients);
184 : const std::map<AtomNumber,Vector> & grad2(v2.gradients);
185 29884 : for(const auto & p1 : grad1) {
186 26836 : AtomNumber a=p1.first;
187 : const auto p2=grad2.find(a);
188 26836 : if(p2!=grad2.end()) {
189 18334 : proj+=dotProduct(p1.second,(*p2).second);
190 : }
191 : }
192 3048 : return proj;
193 : }
194 :
195 15231654 : ActionWithValue* Value::getPntrToAction() {
196 15231654 : plumed_assert( action!=nullptr );
197 15231654 : return action;
198 : }
199 :
200 28680809 : void Value::set(const std::size_t& n, const double& v ) {
201 28680809 : value_set=true;
202 28680809 : if( getRank()==0 ) {
203 25716 : plumed_assert( n==0 );
204 25716 : data[n]=v;
205 25716 : applyPeriodicity(n);
206 28655093 : } else if( !hasDeriv ) {
207 : plumed_dbg_massert( n<data.size(), "failing in " + getName() );
208 26767620 : data[n]=v;
209 26767620 : applyPeriodicity(n);
210 : } else {
211 1887473 : data[n*(1+ngrid_der)] = v;
212 : }
213 28680809 : }
214 :
215 74939 : void Value::push_back( const double& v ) {
216 74939 : value_set=true;
217 74939 : if( shape.size()==1 ) {
218 36380 : data.push_back(v);
219 36380 : shape[0]++;
220 38559 : } else if( shape.size()==2 ) {
221 38559 : data.push_back(v);
222 38559 : shape[0] = ( data.size() + shape[1] - 1 ) / shape[1];
223 : }
224 74939 : }
225 :
226 0 : std::size_t Value::getIndexInStore( const std::size_t& ival ) const {
227 0 : if( shape.size()==2 && ncols<shape[1] ) {
228 0 : unsigned irow = std::floor( ival / shape[1] ), jcol = ival%shape[1];
229 0 : for(unsigned i=0; i<getRowLength(irow); ++i) {
230 0 : if( getRowIndex(irow,i)==jcol ) {
231 0 : return irow*ncols+i;
232 : }
233 : }
234 0 : plumed_merror("cannot get store index");
235 : }
236 0 : return ival;
237 : }
238 :
239 1827201 : bool Value::checkValueIsActiveForMMul(const std::size_t task) const {
240 : const auto ncol = getRowLength(task);
241 1827201 : const auto base = task * getNumberOfColumns();
242 1827201 : if (hasDeriv) {
243 0 : for(std::size_t k=base; k<base+ncol; ++k) {
244 0 : if(std::fabs(data[k*(1+ngrid_der)])>0.0) {
245 : return true;
246 : }
247 : }
248 : return false;
249 : }
250 1827201 : return std::any_of(&data[base],&data[base]+ncol,[](double x) {
251 1274781 : return std::fabs(x)>0.0;
252 1827201 : });
253 :
254 : }
255 :
256 472179667 : double Value::get(const std::size_t ival, const bool trueind) const {
257 472179667 : if( hasDeriv ) {
258 11004253 : return data[ival*(1+ngrid_der)];
259 : }
260 : // plumed_dbg_massert is only evaluated in debug builds (NDEBUG not defined)
261 : if( action ) {
262 : plumed_dbg_massert( ival<getNumberOfValues(), "could not get value from " + name );
263 : }
264 461175414 : if( shape.size()==2 && ncols<shape[1] && trueind ) {
265 8620118 : const unsigned irow = std::floor( ival / shape[1] );
266 : const unsigned jcol = ival%shape[1];
267 : // This is a special treatment for the lower triangular matrices that are used when
268 : // we do ITRE with COLLECT_FRAMES
269 8620118 : if( ncols==0 ) {
270 0 : if( jcol<=irow ) {
271 0 : return data[0.5*irow*(irow+1) + jcol];
272 : }
273 : return 0.0;
274 : }
275 : /* I have to work on this
276 : auto begin = matrix_bookeeping.begin()+(1+ncols)*irow+1;
277 : auto end = matrix_bookeeping.begin()+(1+ncols)*irow+1+getRowLength(irow);
278 : auto i=std::find(begin,end,jcol);
279 : if (i!=end){
280 : return data[irow*ncols+end-i];
281 : }*/
282 2302090094 : for(unsigned i=0; i<getRowLength(irow); ++i) {
283 1150406138 : if( getRowIndex(irow,i)==jcol ) {
284 7981209 : return data[irow*ncols+i];
285 : }
286 : }
287 : return 0.0;
288 : }
289 452555296 : plumed_massert( ival<data.size(), "cannot get value from " + name );
290 452555296 : return data[ival];
291 : }
292 :
293 122065 : size_t Value::assignValues(View<double> target) {
294 122065 : const auto nvals=getNumberOfStoredValues ();
295 122065 : if( hasDeriv ) {
296 6490 : for(std::size_t j=0; j<nvals; ++j) {
297 6484 : target[j] = data[j*(1+ngrid_der)];
298 : }
299 : } else {
300 122059 : plumed_massert( data.size()>=nvals, "cannot get value from " + name );
301 122059 : std::memcpy(target.data(),data.data(),nvals*sizeof(double));
302 : }
303 122065 : return nvals;
304 : }
305 :
306 0 : size_t Value::assignValues(View<float> target) {
307 0 : const auto nvals=getNumberOfStoredValues ();
308 0 : if( hasDeriv ) {
309 0 : for(std::size_t j=0; j<nvals; ++j) {
310 0 : target[j] = data[j*(1+ngrid_der)];
311 : }
312 : } else {
313 0 : plumed_massert( data.size()>=nvals, "cannot get value from " + name );
314 0 : for(std::size_t j=0; j<nvals; ++j) {
315 0 : target[j] = data[j];
316 : }
317 : }
318 0 : return nvals;
319 : }
320 :
321 0 : void Value::addForceOutOfRange( const std::size_t& iforce ) const {
322 0 : plumed_merror( "can't add force to " + name );
323 : }
324 :
325 1425 : void Value::addForceToSparseMatrix(const std::size_t& iforce, double f) {
326 1425 : unsigned irow = iforce / shape[1], jcol = iforce%shape[1];
327 15498 : for(unsigned i=0; i<getRowLength(irow); ++i) {
328 6393 : if( getRowIndex(irow,i)==jcol ) {
329 69 : inputForce[irow*ncols+i]+=f;
330 69 : return;
331 : }
332 : }
333 : // elements that are not stored are zero; when their derivatives are zero too
334 : // (e.g. switching functions beyond the cutoff) the force on them has no effect
335 1356 : plumed_massert( derivativeIsZeroWhenValueIsZero || fabs(f)<epsilon, "force on an element that is not stored in sparse matrix " + name );
336 : }
337 :
338 312183 : size_t Value::addForces(View<const double> const forces) {
339 312183 : hasForce=true;
340 312183 : const auto nvals = getNumberOfStoredValues();
341 312183 : plumed_massert( inputForce.size()>=nvals, "can't add force to " + name );
342 : //I need at least nvals elements in forces
343 312183 : plumed_massert( forces.size()>=nvals, "can't add force to " + name );
344 : //is this daxpy?
345 71261727 : for(auto i=0u; i<nvals; ++i) {
346 70949544 : inputForce[i]+=forces[i];
347 : }
348 312183 : return nvals;
349 : }
350 :
351 8067 : void Value::reshapeMatrixStore( const unsigned& n ) {
352 : plumed_dbg_assert( shape.size()==2 && !hasDeriv );
353 8067 : ncols=n;
354 8067 : if( shape[1]>0 && ncols>shape[1] ) {
355 0 : ncols=shape[1];
356 : }
357 : // ncols==0 is the flag of the lower-triangular storage (see get): a sparse
358 : // matrix whose rows are all empty keeps one column of storage
359 8067 : if( shape[1]>0 && ncols==0 ) {
360 16 : ncols=1;
361 : }
362 8067 : unsigned size=shape[0]*ncols;
363 8067 : bookeeping_root=nullptr;
364 8067 : if( matrix_bookeeping.size()!=(size+shape[0]) ) {
365 4407 : data.resize( size, 0 );
366 4407 : inputForce.resize( size, 0 );
367 4407 : matrix_bookeeping.resize( size + shape[0], 0 );
368 4407 : if( ncols>=shape[1] ) {
369 271819 : for(unsigned i=0; i<shape[0]; ++i) {
370 268377 : matrix_bookeeping[(1+ncols)*i] = shape[1];
371 16500065 : for(unsigned j=0; j<shape[1]; ++j) {
372 16231688 : matrix_bookeeping[(1+ncols)*i+1+j]=j;
373 : }
374 : }
375 : }
376 : }
377 8067 : }
378 :
379 7076 : void Value::copyBookeepingArrayFromArgument( Value* myarg ) {
380 : plumed_dbg_assert( shape.size()==2 && !hasDeriv );
381 7076 : ncols = myarg->getNumberOfColumns();
382 7076 : matrix_bookeeping.resize( myarg->matrix_bookeeping.size() );
383 7076 : if( !matrix_bookeeping.empty() ) {
384 7076 : std::memcpy( matrix_bookeeping.data(), myarg->matrix_bookeeping.data(), matrix_bookeeping.size()*sizeof(matrix_bookeeping[0]) );
385 : }
386 7076 : bookeeping_root = myarg->getBookeepingRoot();
387 7076 : data.resize( shape[0]*ncols );
388 7076 : inputForce.resize( shape[0]*ncols );
389 7076 : }
390 :
391 0 : bool Value::ignoreStoredValue(const std::string& c) const {
392 0 : return false;
393 : }
394 :
395 8085 : void Value::setConstant() {
396 8085 : valtype=constant;
397 8085 : setShape( shape );
398 8085 : if( getRank()==2 && !hasDeriv ) {
399 149 : reshapeMatrixStore( shape[1] );
400 : }
401 8085 : derivativeIsZeroWhenValueIsZero=true;
402 8085 : }
403 :
404 1 : void Value::reshapeConstantValue( const std::vector<std::size_t>& sh ) {
405 1 : plumed_assert( valtype==constant && getNumberOfValues()==1 );
406 1 : double val = get(0);
407 1 : setShape( sh );
408 1 : if( getRank()==2 && !hasDeriv ) {
409 0 : reshapeMatrixStore( sh[1] );
410 : }
411 3 : for(unsigned i=0; i<getNumberOfValues(); ++i) {
412 2 : set( i, val );
413 : }
414 1 : }
415 :
416 504 : void Value::writeBinary(std::ostream&o) const {
417 504 : o.write(reinterpret_cast<const char*>(&data[0]),data.size()*sizeof(double));
418 504 : }
419 :
420 1355 : void Value::setSymmetric( const bool& sym ) {
421 1355 : plumed_assert( shape.size()==2 && !hasDeriv );
422 1355 : if( sym && shape[0]!=shape[1] ) {
423 0 : plumed_merror("non-square matrix cannot be symmetric");
424 : }
425 1355 : symmetric=sym;
426 1355 : }
427 :
428 1157 : void Value::retrieveEdgeList( unsigned& nedge, std::vector<std::pair<unsigned,unsigned> >& active, std::vector<double>& elems ) {
429 1157 : nedge=0;
430 : plumed_dbg_assert( shape.size()==2 && !hasDeriv );
431 : // Check we have enough space to store the edge list
432 1157 : if( elems.size()<shape[0]*ncols ) {
433 372 : elems.resize( shape[0]*ncols );
434 372 : active.resize( shape[0]*ncols );
435 : }
436 :
437 34035 : for(unsigned i=0; i<shape[0]; ++i) {
438 : unsigned ncol = getRowLength(i);
439 8717348 : for(unsigned j=0; j<ncol; ++j) {
440 8684470 : if( fabs(get(i*ncols+j,false))<epsilon ) {
441 8150688 : continue;
442 : }
443 533782 : if( symmetric && getRowIndex(i,j)>i ) {
444 58922 : continue;
445 : }
446 474860 : active[nedge].first = i;
447 474860 : active[nedge].second = getRowIndex(i,j);
448 474860 : elems[nedge] = get(i*ncols+j,false);
449 474860 : nedge++;
450 : }
451 : }
452 1157 : }
453 :
454 504 : void Value::readBinary(std::istream&i) {
455 504 : i.read(reinterpret_cast<char*>(&data[0]),data.size()*sizeof(double));
456 504 : }
457 :
458 1090311 : void Value::convertIndexToindices(const std::size_t& index, std::vector<unsigned>& indices ) const {
459 1090311 : if( hasDeriv || getRank()==1 ) {
460 113508 : std::size_t kk=index;
461 113508 : indices[0]=index%shape[0];
462 113508 : for(unsigned i=1; i<shape.size()-1; ++i) {
463 0 : kk=(kk-indices[i-1])/shape[i-1];
464 0 : indices[i]=kk%shape[i];
465 : }
466 113508 : if(shape.size()>=2) {
467 0 : indices[shape.size()-1]=(kk-indices[shape.size()-2])/shape[shape.size()-2];
468 : }
469 976803 : } else if( getRank()==2 ) {
470 976803 : indices[0]=std::floor( index/shape[1] );
471 976803 : indices[1] = index%shape[1];
472 : }
473 1090311 : }
474 :
475 634487 : void Value::print( OFile& ofile ) const {
476 634487 : if( isPeriodic() ) {
477 24462 : ofile.printField( "min_" + name, str_min );
478 48924 : ofile.printField("max_" + name, str_max );
479 : }
480 634487 : if( shape.size()==0 || getNumberOfValues()==1 ) {
481 632809 : ofile.printField( name, get(0) );
482 : } else {
483 1678 : std::vector<unsigned> indices( shape.size() );
484 1091769 : for(unsigned i=0; i<getNumberOfValues(); ++i) {
485 1090091 : convertIndexToindices( i, indices );
486 1090091 : std::string num, fname = name;
487 3156985 : for(unsigned ii=0; ii<shape.size(); ++ii) {
488 2066894 : Tools::convert( indices[ii]+1, num );
489 2066894 : fname += "." + num;
490 : }
491 1090091 : ofile.printField( fname, get(i) );
492 : }
493 1678 : }
494 2701381 : }
495 :
496 10454 : void Value::printForce( OFile& ofile ) const {
497 10454 : if( shape.size()==0 || getNumberOfValues()==1 ) {
498 10389 : ofile.printField( name, getForce(0) );
499 : } else {
500 65 : std::vector<unsigned> indices( shape.size() );
501 285 : for(unsigned i=0; i<getNumberOfValues(); ++i) {
502 220 : convertIndexToindices( i, indices );
503 220 : std::string num, fname = name;
504 440 : for(unsigned ii=0; ii<shape.size(); ++ii) {
505 220 : Tools::convert( indices[ii]+1, num );
506 440 : fname += "." + num;
507 : }
508 220 : plumed_assert( i<inputForce.size() );
509 220 : ofile.printField( fname, getForce(i) );
510 : }
511 65 : }
512 10454 : }
513 :
514 0 : unsigned Value::getGoodNumThreads( const unsigned& j, const unsigned& k ) const {
515 0 : return OpenMP::getGoodNumThreads( &data[j], (k-j) );
516 : }
517 :
518 9026476 : bool Value::calculateOnUpdate() const {
519 9026476 : return (valtype==average || valtype==calcFromAverage);
520 : }
521 :
522 129 : std::string Value::getValueType() const {
523 129 : if( getRank()==0 ) {
524 31 : return "scalar";
525 : }
526 98 : if( getRank()>0 && hasDerivatives() ) {
527 5 : return "grid";
528 : }
529 93 : if( getRank()==1 ) {
530 72 : return "vector";
531 : }
532 21 : if( getRank()==2 ) {
533 21 : return "matrix";
534 : }
535 0 : plumed_merror("unknown type for value " + getName() );
536 : return "";
537 : }
538 :
539 119 : bool Value::allElementsEqual() const {
540 119 : if( getNumberOfValues()==0 ) {
541 : return false;
542 : }
543 :
544 119 : double refval = get(0);
545 1411783 : for(unsigned i=1; i<getNumberOfValues(); ++i) {
546 1411664 : if( fabs(get(i)-refval)>epsilon ) {
547 : return false;
548 : }
549 : }
550 : return true;
551 : }
552 :
553 3463 : std::size_t Value::getLengthOfLongestColumn() const {
554 : plumed_dbg_assert( shape.size()==2 && !hasDeriv );
555 3463 : std::vector<unsigned> column_lengths( shape[1], 0 );
556 356126 : for(unsigned i=0; i<shape[0]; ++i) {
557 : unsigned nr = getRowLength(i);
558 10598689 : for(unsigned j=0; j<nr; ++j) {
559 10246026 : unsigned ind = getRowIndex( i, j );
560 10246026 : column_lengths[ind]++;
561 : }
562 : }
563 : // Find the longest column
564 3463 : std::size_t maxcol = column_lengths[0];
565 293781 : for(unsigned i=1; i<column_lengths.size(); ++i) {
566 290318 : if( column_lengths[i]>maxcol ) {
567 : maxcol = column_lengths[i];
568 : }
569 : }
570 3463 : return maxcol;
571 3463 : }
572 :
573 : }
|
