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 "tools/Exception.h"
26 : #include "tools/OpenMP.h"
27 : #include "tools/OFile.h"
28 : #include "tools/Communicator.h"
29 : #include "PlumedMain.h"
30 :
31 : namespace PLMD {
32 :
33 2 : Value::Value():
34 4 : data(1,0.0),
35 2 : inputForce(1,0.0),
36 4 : shape() {}
37 :
38 18 : Value::Value(const std::string& valname):
39 36 : data(1,0.0),
40 18 : inputForce(1,0.0),
41 18 : name(valname),
42 36 : shape() {
43 18 : }
44 :
45 118004 : Value::Value(ActionWithValue* av,
46 : const std::string& valname,
47 : const bool withderiv,
48 118004 : const std::vector<std::size_t>&ss):
49 118004 : action(av),
50 118004 : name(valname),
51 118004 : hasDeriv(withderiv) {
52 118004 : if( action ) {
53 117548 : if( action->getName()=="ACCUMULATE" || action->getName()=="COLLECT" ) {
54 159 : valtype=average;
55 : }
56 : }
57 118004 : setShape( ss );
58 118004 : }
59 :
60 448 : void Value::setValType( const std::string& vtype ) {
61 448 : if( vtype=="normal" ) {
62 0 : valtype=normal;
63 448 : } else if( vtype=="constant" ) {
64 0 : valtype=constant;
65 448 : } else if( vtype=="average" ) {
66 0 : valtype=average;
67 448 : } else if( vtype=="calcFromAverage" ) {
68 448 : valtype=calcFromAverage;
69 : } else {
70 0 : plumed_merror("invalid valtype " + vtype );
71 : }
72 448 : }
73 :
74 123397 : void Value::setShape( const std::vector<std::size_t>&ss ) {
75 : std::size_t tot=1;
76 123397 : shape.resize( ss.size() );
77 145167 : for(unsigned i=0; i<shape.size(); ++i) {
78 21770 : tot = tot*ss[i];
79 21770 : shape[i]=ss[i];
80 : }
81 :
82 123397 : if( shape.size()>0 && hasDeriv ) {
83 : // This is for grids
84 511 : ngrid_der = shape.size();
85 511 : if( action ) {
86 510 : ngrid_der = action->getNumberOfDerivatives();
87 : }
88 511 : std::size_t ndata = tot*(1+ngrid_der);
89 511 : data.resize( ndata );
90 511 : inputForce.resize( tot );
91 122886 : } else if( shape.size()==0 ) {
92 : // This is for scalars
93 105480 : data.resize(1);
94 105480 : inputForce.resize(1);
95 17406 : } else if( shape.size()<2 ) {
96 : // This is for vectors (matrices have special version because we have sparse storage)
97 13813 : data.resize( tot );
98 13813 : inputForce.resize( tot );
99 : }
100 123397 : }
101 :
102 108050 : void Value::setupPeriodicity() {
103 108050 : if( min==0 && max==0 ) {
104 100320 : periodicity=notperiodic;
105 : } else {
106 7730 : periodicity=periodic;
107 7730 : max_minus_min=max-min;
108 7730 : plumed_massert(max_minus_min>0, "your function has a very strange domain?");
109 7730 : inv_max_minus_min=1.0/max_minus_min;
110 : }
111 108050 : }
112 :
113 1698627 : bool Value::isPeriodic()const {
114 1698627 : plumed_massert(periodicity!=unset,"periodicity should be set");
115 1698627 : return periodicity==periodic;
116 : }
117 :
118 429627 : bool Value::applyForce(std::vector<double>& forces ) const {
119 429627 : if( !hasForce || valtype!=normal ) {
120 : return false;
121 : }
122 : plumed_dbg_massert( data.size()-1==forces.size()," forces array has wrong size" );
123 9579 : const unsigned N=data.size()-1;
124 41730663 : for(unsigned i=0; i<N; ++i) {
125 41721084 : forces[i]=inputForce[0]*data[1+i];
126 : }
127 : return true;
128 : }
129 :
130 1105012 : void Value::setNotPeriodic() {
131 1105012 : min=0;
132 1105012 : max=0;
133 1105012 : periodicity=notperiodic;
134 1105012 : }
135 :
136 7730 : void Value::setDomain(const std::string& pmin,const std::string& pmax) {
137 7730 : str_min=pmin;
138 7730 : if( !Tools::convertNoexcept(str_min,min) ) {
139 0 : action->error("could not convert period string " + str_min + " to real");
140 : }
141 7730 : str_max=pmax;
142 7730 : if( !Tools::convertNoexcept(str_max,max) ) {
143 0 : action->error("could not convert period string " + str_max + " to read");
144 : }
145 7730 : setupPeriodicity();
146 7730 : }
147 :
148 16878 : void Value::getDomain(std::string&minout,std::string&maxout) const {
149 16878 : plumed_massert(periodicity==periodic,"function should be periodic");
150 16878 : minout=str_min;
151 16878 : maxout=str_max;
152 16878 : }
153 :
154 34 : void Value::getDomain(double&minout,double&maxout) const {
155 34 : plumed_massert(periodicity==periodic,"function should be periodic");
156 34 : minout=min;
157 34 : maxout=max;
158 34 : }
159 :
160 190 : void Value::setGradients( ActionAtomistic* aa, unsigned& start ) {
161 : // Can't do gradients if we don't have derivatives
162 190 : if( !hasDeriv ) {
163 : return;
164 : }
165 154 : plumed_assert( shape.size()==0 );
166 8362 : for(unsigned j=0; j<aa->getNumberOfAtoms(); ++j) {
167 8208 : Vector der(data[1+start+3*j],data[1+start+3*j+1],data[1+start+3*j+2]);
168 8208 : aa->getGradient( j, der, gradients );
169 : }
170 154 : start += aa->getNumberOfAtoms();
171 : }
172 :
173 283 : void Value::passGradients( const double& der, std::map<AtomNumber,Vector>& g ) const {
174 8921 : for(const auto & p : gradients) {
175 8638 : AtomNumber iatom=p.first;
176 8638 : g[iatom] += p.second*der;
177 : }
178 283 : }
179 :
180 261 : double Value::projection(const Value& v1,const Value&v2) {
181 : double proj=0.0;
182 : const std::map<AtomNumber,Vector> & grad1(v1.gradients);
183 : const std::map<AtomNumber,Vector> & grad2(v2.gradients);
184 16167 : for(const auto & p1 : grad1) {
185 15906 : AtomNumber a=p1.first;
186 : const auto p2=grad2.find(a);
187 15906 : if(p2!=grad2.end()) {
188 8224 : proj+=dotProduct(p1.second,(*p2).second);
189 : }
190 : }
191 261 : return proj;
192 : }
193 :
194 14718824 : ActionWithValue* Value::getPntrToAction() {
195 14718824 : plumed_assert( action!=nullptr );
196 14718824 : return action;
197 : }
198 :
199 599565821 : void Value::set(const std::size_t& n, const double& v ) {
200 599565821 : value_set=true;
201 599565821 : if( getRank()==0 ) {
202 25636 : plumed_assert( n==0 );
203 25636 : data[n]=v;
204 25636 : applyPeriodicity(n);
205 599540185 : } else if( !hasDeriv ) {
206 : plumed_dbg_massert( n<data.size(), "failing in " + getName() );
207 597679425 : data[n]=v;
208 597679425 : applyPeriodicity(n);
209 : } else {
210 1860760 : data[n*(1+ngrid_der)] = v;
211 : }
212 599565821 : }
213 :
214 74915 : void Value::push_back( const double& v ) {
215 74915 : value_set=true;
216 74915 : if( shape.size()==1 ) {
217 36356 : data.push_back(v);
218 36356 : shape[0]++;
219 38559 : } else if( shape.size()==2 ) {
220 38559 : data.push_back(v);
221 38559 : shape[0] = std::ceil( data.size() / shape[1] );
222 : }
223 74915 : }
224 :
225 0 : std::size_t Value::getIndexInStore( const std::size_t& ival ) const {
226 0 : if( shape.size()==2 && ncols<shape[1] ) {
227 0 : unsigned irow = std::floor( ival / shape[1] ), jcol = ival%shape[1];
228 0 : for(unsigned i=0; i<getRowLength(irow); ++i) {
229 0 : if( getRowIndex(irow,i)==jcol ) {
230 0 : return irow*ncols+i;
231 : }
232 : }
233 0 : plumed_merror("cannot get store index");
234 : }
235 0 : return ival;
236 : }
237 :
238 1824673 : bool Value::checkValueIsActiveForMMul(const std::size_t task) const {
239 : const auto ncol = getRowLength(task);
240 1824673 : const auto base = task * getNumberOfColumns();
241 1824673 : if (hasDeriv) {
242 0 : for(std::size_t k=base; k<base+ncol; ++k) {
243 0 : if(std::fabs(data[k*(1+ngrid_der)])>0.0) {
244 : return true;
245 : }
246 : }
247 : return false;
248 : }
249 1824673 : return std::any_of(&data[base],&data[base]+ncol,[](double x) {
250 92409090 : return std::fabs(x)>0.0;
251 1824673 : });
252 :
253 : }
254 :
255 549231947 : double Value::get(const std::size_t ival, const bool trueind) const {
256 549231947 : if( hasDeriv ) {
257 10411340 : return data[ival*(1+ngrid_der)];
258 : }
259 : #ifdef DNDEBUG
260 : if( action ) {
261 : plumed_dbg_massert( ival<getNumberOfValues(), "could not get value from " + name );
262 : }
263 : #endif
264 538820607 : if( shape.size()==2 && trueind ) {
265 296697025 : const unsigned irow = std::floor( ival / shape[1] );
266 296697025 : 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 296697025 : 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 33322214998 : for(unsigned i=0; i<getRowLength(irow); ++i) {
283 16660907123 : if( getRowIndex(irow,i)==jcol ) {
284 296496649 : return data[irow*ncols+i];
285 : }
286 : }
287 : return 0.0;
288 : }
289 242123582 : plumed_massert( ival<data.size(), "cannot get value from " + name );
290 242123582 : return data[ival];
291 : }
292 :
293 86113 : size_t Value::assignValues(View<double> target) {
294 86113 : const auto nvals=getNumberOfStoredValues ();
295 86113 : 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 86107 : plumed_massert( data.size()>=nvals, "cannot get value from " + name );
301 86107 : std::memcpy(target.data(),data.data(),nvals*sizeof(double));
302 : }
303 86113 : return nvals;
304 : }
305 :
306 9351187 : void Value::addForce(const std::size_t& iforce, double f, const bool trueind) {
307 9351187 : hasForce=true;
308 9351187 : if( shape.size()==2 && !hasDeriv && ncols<shape[1] && trueind ) {
309 0 : unsigned irow = std::floor( iforce / shape[0] ), jcol = iforce%shape[0];
310 0 : for(unsigned i=0; i<getRowLength(irow); ++i) {
311 0 : if( getRowIndex(irow,i)==jcol ) {
312 0 : inputForce[irow*ncols+i]+=f;
313 0 : return;
314 : }
315 : }
316 0 : plumed_assert( fabs(f)<epsilon );
317 : return;
318 : }
319 9351187 : plumed_massert( iforce<inputForce.size(), "can't add force to " + name );
320 9351187 : inputForce[iforce]+=f;
321 : }
322 :
323 338028 : size_t Value::addForces(View<const double> const forces) {
324 338028 : hasForce=true;
325 338028 : const auto nvals = getNumberOfStoredValues();
326 338028 : plumed_massert( inputForce.size()>=nvals, "can't add force to " + name );
327 : //I need at least nvals elements in forces
328 338028 : plumed_massert( forces.size()>=nvals, "can't add force to " + name );
329 : /*
330 : {//this gives a very little speedup (+1 step in the 60s dragrace)
331 : auto f = forces.begin();
332 : const auto end=inputForce.begin()+nvals;
333 : for(auto iptf=inputForce.begin() ; iptf<end; ++iptf, ++f) {
334 : *iptf+=*f;
335 : }
336 : }
337 : */
338 : //is this daxpy?
339 180959438 : for(auto i=0u; i<nvals; ++i) {
340 180621410 : inputForce[i]+=forces[i];
341 : }
342 338028 : return nvals;
343 : }
344 :
345 4259 : void Value::reshapeMatrixStore( const unsigned& n ) {
346 : plumed_dbg_assert( shape.size()==2 && !hasDeriv );
347 4259 : ncols=n;
348 4259 : if( shape[1]>0 && ncols>shape[1] ) {
349 0 : ncols=shape[1];
350 : }
351 4259 : unsigned size=shape[0]*ncols;
352 4259 : if( matrix_bookeeping.size()!=(size+shape[0]) ) {
353 3224 : data.resize( size, 0 );
354 3224 : inputForce.resize( size, 0 );
355 3224 : matrix_bookeeping.resize( size + shape[0], 0 );
356 3224 : if( ncols>=shape[1] ) {
357 277813 : for(unsigned i=0; i<shape[0]; ++i) {
358 275312 : matrix_bookeeping[(1+ncols)*i] = shape[1];
359 87628346 : for(unsigned j=0; j<shape[1]; ++j) {
360 87353034 : matrix_bookeeping[(1+ncols)*i+1+j]=j;
361 : }
362 : }
363 : }
364 : }
365 4259 : }
366 :
367 6435 : void Value::copyBookeepingArrayFromArgument( Value* myarg ) {
368 : plumed_dbg_assert( shape.size()==2 && !hasDeriv );
369 6435 : ncols = myarg->getNumberOfColumns();
370 6435 : matrix_bookeeping.resize( myarg->matrix_bookeeping.size() );
371 602329564 : for(unsigned i=0; i<matrix_bookeeping.size(); ++i) {
372 602323129 : matrix_bookeeping[i] = myarg->matrix_bookeeping[i];
373 : }
374 6435 : data.resize( shape[0]*ncols );
375 6435 : inputForce.resize( shape[0]*ncols );
376 6435 : }
377 :
378 0 : bool Value::ignoreStoredValue(const std::string& c) const {
379 0 : return false;
380 : }
381 :
382 4937 : void Value::setConstant() {
383 4937 : valtype=constant;
384 4937 : setShape( shape );
385 4937 : if( getRank()==2 && !hasDeriv ) {
386 137 : reshapeMatrixStore( shape[1] );
387 : }
388 4937 : derivativeIsZeroWhenValueIsZero=true;
389 4937 : }
390 :
391 1 : void Value::reshapeConstantValue( const std::vector<std::size_t>& sh ) {
392 1 : plumed_assert( valtype==constant && getNumberOfValues()==1 );
393 1 : double val = get(0);
394 1 : setShape( sh );
395 1 : if( getRank()==2 && !hasDeriv ) {
396 0 : reshapeMatrixStore( sh[1] );
397 : }
398 3 : for(unsigned i=0; i<getNumberOfValues(); ++i) {
399 2 : set( i, val );
400 : }
401 1 : }
402 :
403 456 : void Value::writeBinary(std::ostream&o) const {
404 456 : o.write(reinterpret_cast<const char*>(&data[0]),data.size()*sizeof(double));
405 456 : }
406 :
407 1173 : void Value::setSymmetric( const bool& sym ) {
408 1173 : plumed_assert( shape.size()==2 && !hasDeriv );
409 1173 : if( sym && shape[0]!=shape[1] ) {
410 0 : plumed_merror("non-square matrix cannot be symmetric");
411 : }
412 1173 : symmetric=sym;
413 1173 : }
414 :
415 3412 : void Value::retrieveEdgeList( unsigned& nedge, std::vector<std::pair<unsigned,unsigned> >& active, std::vector<double>& elems ) {
416 3412 : nedge=0;
417 : plumed_dbg_assert( shape.size()==2 && !hasDeriv );
418 : // Check we have enough space to store the edge list
419 3412 : if( elems.size()<shape[0]*ncols ) {
420 2727 : elems.resize( shape[0]*ncols );
421 2727 : active.resize( shape[0]*ncols );
422 : }
423 :
424 361355 : for(unsigned i=0; i<shape[0]; ++i) {
425 : unsigned ncol = getRowLength(i);
426 17331451 : for(unsigned j=0; j<ncol; ++j) {
427 16973508 : if( fabs(get(i*ncols+j,false))<epsilon ) {
428 10342139 : continue;
429 : }
430 6631369 : if( symmetric && getRowIndex(i,j)>i ) {
431 51540 : continue;
432 : }
433 6579829 : active[nedge].first = i;
434 6579829 : active[nedge].second = getRowIndex(i,j);
435 6579829 : elems[nedge] = get(i*ncols+j,false);
436 6579829 : nedge++;
437 : }
438 : }
439 3412 : }
440 :
441 456 : void Value::readBinary(std::istream&i) {
442 456 : i.read(reinterpret_cast<char*>(&data[0]),data.size()*sizeof(double));
443 456 : }
444 :
445 731483 : void Value::convertIndexToindices(const std::size_t& index, std::vector<unsigned>& indices ) const {
446 731483 : if( hasDeriv || getRank()==1 ) {
447 13537 : std::size_t kk=index;
448 13537 : indices[0]=index%shape[0];
449 13537 : for(unsigned i=1; i<shape.size()-1; ++i) {
450 0 : kk=(kk-indices[i-1])/shape[i-1];
451 0 : indices[i]=kk%shape[i];
452 : }
453 13537 : if(shape.size()>=2) {
454 0 : indices[shape.size()-1]=(kk-indices[shape.size()-2])/shape[shape.size()-2];
455 : }
456 717946 : } else if( getRank()==2 ) {
457 717946 : indices[0]=std::floor( index/shape[1] );
458 717946 : indices[1] = index%shape[1];
459 : }
460 731483 : }
461 :
462 556218 : void Value::print( OFile& ofile ) const {
463 556218 : if( isPeriodic() ) {
464 11542 : ofile.printField( "min_" + name, str_min );
465 23084 : ofile.printField("max_" + name, str_max );
466 : }
467 556218 : if( shape.size()==0 || getNumberOfValues()==1 ) {
468 554817 : ofile.printField( name, get(0) );
469 : } else {
470 1401 : std::vector<unsigned> indices( shape.size() );
471 732784 : for(unsigned i=0; i<getNumberOfValues(); ++i) {
472 731383 : convertIndexToindices( i, indices );
473 731383 : std::string num, fname = name;
474 2180712 : for(unsigned ii=0; ii<shape.size(); ++ii) {
475 1449329 : Tools::convert( indices[ii]+1, num );
476 2898658 : fname += "." + num;
477 : }
478 731383 : ofile.printField( fname, get(i) );
479 : }
480 : }
481 556218 : }
482 :
483 4693 : void Value::printForce( OFile& ofile ) const {
484 4693 : if( shape.size()==0 || getNumberOfValues()==1 ) {
485 4688 : ofile.printField( name, getForce(0) );
486 : } else {
487 5 : std::vector<unsigned> indices( shape.size() );
488 105 : for(unsigned i=0; i<getNumberOfValues(); ++i) {
489 100 : convertIndexToindices( i, indices );
490 100 : std::string num, fname = name;
491 200 : for(unsigned ii=0; ii<shape.size(); ++ii) {
492 100 : Tools::convert( indices[ii]+1, num );
493 200 : fname += "." + num;
494 : }
495 100 : plumed_assert( i<inputForce.size() );
496 100 : ofile.printField( fname, getForce(i) );
497 : }
498 : }
499 4693 : }
500 :
501 256790 : unsigned Value::getGoodNumThreads( const unsigned& j, const unsigned& k ) const {
502 256790 : return OpenMP::getGoodNumThreads( &data[j], (k-j) );
503 : }
504 :
505 7691290 : bool Value::calculateOnUpdate() const {
506 7691290 : return (valtype==average || valtype==calcFromAverage);
507 : }
508 :
509 129 : std::string Value::getValueType() const {
510 129 : if( getRank()==0 ) {
511 31 : return "scalar";
512 : }
513 98 : if( getRank()>0 && hasDerivatives() ) {
514 5 : return "grid";
515 : }
516 93 : if( getRank()==1 ) {
517 72 : return "vector";
518 : }
519 21 : if( getRank()==2 ) {
520 21 : return "matrix";
521 : }
522 0 : plumed_merror("unknown type for value " + getName() );
523 : return "";
524 : }
525 :
526 102 : bool Value::allElementsEqual() const {
527 102 : if( getNumberOfValues()==0 ) {
528 : return false;
529 : }
530 :
531 102 : double refval = get(0);
532 757335 : for(unsigned i=1; i<getNumberOfValues(); ++i) {
533 757233 : if( fabs(get(i)-refval)>epsilon ) {
534 : return false;
535 : }
536 : }
537 : return true;
538 : }
539 :
540 : }
|
