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 : {
56 164 : data.resize(1); 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 : {
81 18 : data.resize(1); inputForce.resize(1);
82 18 : data[0]=inputForce[0]=0;
83 18 : }
84 :
85 118726 : Value::Value(ActionWithValue* av, const std::string& name, const bool withderiv, const std::vector<unsigned>&ss):
86 118726 : action(av),
87 118726 : value_set(false),
88 118726 : hasForce(false),
89 118726 : name(name),
90 118726 : storedata(false),
91 118726 : hasDeriv(withderiv),
92 118726 : bufstart(0),
93 118726 : streampos(0),
94 118726 : ngrid_der(0),
95 118726 : matpos(0),
96 118726 : ncols(0),
97 118726 : book_start(0),
98 118726 : symmetric(false),
99 118726 : periodicity(unset),
100 118726 : min(0.0),
101 118726 : max(0.0),
102 118726 : max_minus_min(0.0),
103 118726 : inv_max_minus_min(0.0),
104 118726 : derivativeIsZeroWhenValueIsZero(false)
105 : {
106 118726 : if( action ) {
107 118270 : if( action->getName()=="ACCUMULATE" || action->getName()=="COLLECT" ) valtype=average;
108 : }
109 229756 : if( action ) storedata=action->getName()=="PUT" || valtype==average;
110 118726 : if( ss.size() && withderiv ) storedata=true;
111 118726 : setShape( ss );
112 118726 : }
113 :
114 455 : void Value::setValType( const std::string& vtype ) {
115 455 : if( vtype=="normal" ) valtype=normal;
116 455 : else if( vtype=="constant" ) valtype=constant;
117 455 : else if( vtype=="average" ) valtype=average;
118 455 : else if( vtype=="calcFromAverage" ) valtype=calcFromAverage;
119 0 : else plumed_merror("invalid valtype " + vtype );
120 455 : }
121 :
122 131049 : void Value::setShape( const std::vector<unsigned>&ss ) {
123 131049 : std::size_t tot=1; shape.resize( ss.size() );
124 161855 : for(unsigned i=0; i<shape.size(); ++i) { tot = tot*ss[i]; shape[i]=ss[i]; }
125 :
126 131049 : if( shape.size()>0 && hasDeriv ) {
127 : // This is for grids
128 2005 : ngrid_der = shape.size();
129 2005 : if( action ) ngrid_der = action->getNumberOfDerivatives();
130 2005 : std::size_t ndata = tot*(1+ngrid_der);
131 2005 : data.resize( ndata ); inputForce.resize( tot );
132 129044 : } else if( shape.size()==0 ) {
133 : // This is for scalars
134 106522 : data.resize(1); inputForce.resize(1);
135 22522 : } else if( storedata && shape.size()<2 ) {
136 : // This is for vectors (matrices have special version because we have sparse storage)
137 12419 : data.resize( tot ); inputForce.resize( tot );
138 : }
139 131049 : }
140 :
141 109806 : void Value::setupPeriodicity() {
142 109806 : if( min==0 && max==0 ) {
143 102039 : periodicity=notperiodic;
144 : } else {
145 7767 : periodicity=periodic;
146 7767 : max_minus_min=max-min;
147 7767 : plumed_massert(max_minus_min>0, "your function has a very strange domain?");
148 7767 : inv_max_minus_min=1.0/max_minus_min;
149 : }
150 109806 : }
151 :
152 1685651 : bool Value::isPeriodic()const {
153 1685651 : plumed_massert(periodicity!=unset,"periodicity should be set");
154 1685651 : return periodicity==periodic;
155 : }
156 :
157 409249 : bool Value::applyForce(std::vector<double>& forces ) const {
158 409249 : if( !hasForce || valtype!=normal ) return false;
159 : plumed_dbg_massert( data.size()-1==forces.size()," forces array has wrong size" );
160 9165 : const unsigned N=data.size()-1;
161 41672481 : for(unsigned i=0; i<N; ++i) forces[i]=inputForce[0]*data[1+i];
162 : return true;
163 : }
164 :
165 1102468 : void Value::setNotPeriodic() {
166 1102468 : min=0; max=0; periodicity=notperiodic;
167 1102468 : }
168 :
169 7767 : void Value::setDomain(const std::string& pmin,const std::string& pmax) {
170 7767 : str_min=pmin;
171 7767 : if( !Tools::convertNoexcept(str_min,min) ) action->error("could not convert period string " + str_min + " to real");
172 7767 : str_max=pmax;
173 7767 : if( !Tools::convertNoexcept(str_max,max) ) action->error("could not convert period string " + str_max + " to read");
174 7767 : setupPeriodicity();
175 7767 : }
176 :
177 16759 : void Value::getDomain(std::string&minout,std::string&maxout) const {
178 16759 : plumed_massert(periodicity==periodic,"function should be periodic");
179 16759 : minout=str_min;
180 16759 : maxout=str_max;
181 16759 : }
182 :
183 34 : void Value::getDomain(double&minout,double&maxout) const {
184 34 : plumed_massert(periodicity==periodic,"function should be periodic");
185 34 : minout=min;
186 34 : maxout=max;
187 34 : }
188 :
189 190 : void Value::setGradients( ActionAtomistic* aa, unsigned& start ) {
190 : // Can't do gradients if we don't have derivatives
191 190 : if( !hasDeriv ) return;
192 154 : plumed_assert( shape.size()==0 );
193 8362 : for(unsigned j=0; j<aa->getNumberOfAtoms(); ++j) {
194 8208 : Vector der(data[1+start+3*j],data[1+start+3*j+1],data[1+start+3*j+2]);
195 8208 : aa->getGradient( j, der, gradients );
196 : }
197 154 : start += aa->getNumberOfAtoms();
198 : }
199 :
200 283 : void Value::passGradients( const double& der, std::map<AtomNumber,Vector>& g ) const {
201 8921 : for(const auto & p : gradients) { AtomNumber iatom=p.first; g[iatom] += p.second*der; }
202 283 : }
203 :
204 261 : double Value::projection(const Value& v1,const Value&v2) {
205 : double proj=0.0;
206 : const std::map<AtomNumber,Vector> & grad1(v1.gradients);
207 : const std::map<AtomNumber,Vector> & grad2(v2.gradients);
208 16167 : for(const auto & p1 : grad1) {
209 15906 : AtomNumber a=p1.first;
210 : const auto p2=grad2.find(a);
211 15906 : if(p2!=grad2.end()) {
212 8224 : proj+=dotProduct(p1.second,(*p2).second);
213 : }
214 : }
215 261 : return proj;
216 : }
217 :
218 18718290 : ActionWithValue* Value::getPntrToAction() {
219 18718290 : plumed_assert( action!=NULL );
220 18718290 : return action;
221 : }
222 :
223 4948688 : void Value::set(const std::size_t& n, const double& v ) {
224 4948688 : value_set=true;
225 4948688 : if( getRank()==0 ) { plumed_assert( n==0 ); data[n]=v; applyPeriodicity(n); }
226 4923013 : else if( !hasDeriv ) { plumed_dbg_massert( n<data.size(), "failing in " + getName() ); data[n]=v; applyPeriodicity(n); }
227 317316 : else { data[n*(1+ngrid_der)] = v; }
228 4948688 : }
229 :
230 74915 : void Value::push_back( const double& v ) {
231 74915 : value_set=true;
232 74915 : if( shape.size()==1 ) {
233 36356 : data.push_back(v); shape[0]++;
234 38559 : } else if( shape.size()==2 ) {
235 38559 : data.push_back(v);
236 38559 : shape[0] = std::ceil( data.size() / shape[1] );
237 : }
238 74915 : }
239 :
240 3223564 : std::size_t Value::getIndexInStore( const std::size_t& ival ) const {
241 3223564 : if( shape.size()==2 && ncols<shape[1] ) {
242 2712800 : unsigned irow = std::floor( ival / shape[1] ), jcol = ival%shape[1];
243 11801200 : for(unsigned i=0; i<getRowLength(irow); ++i) {
244 11801200 : if( getRowIndex(irow,i)==jcol ) return irow*ncols+i;
245 : }
246 0 : plumed_merror("cannot get store index");
247 : }
248 510764 : return ival;
249 : }
250 :
251 557639572 : double Value::get(const std::size_t& ival, const bool trueind) const {
252 557639572 : if( hasDeriv ) return data[ival*(1+ngrid_der)];
253 : #ifdef DNDEBUG
254 : if( action ) plumed_dbg_massert( ival<getNumberOfValues(), "could not get value from " + name );
255 : #endif
256 542726586 : if( shape.size()==2 && ncols<shape[1] && trueind ) {
257 2000962 : unsigned irow = std::floor( ival / shape[1] ), jcol = ival%shape[1];
258 : // This is a special treatment for the lower triangular matrices that are used when
259 : // we do ITRE with COLLECT_FRAMES
260 2000962 : if( ncols==0 ) {
261 0 : if( jcol<=irow ) return data[0.5*irow*(irow+1) + jcol];
262 : return 0;
263 : }
264 20655470 : for(unsigned i=0; i<getRowLength(irow); ++i) {
265 19057356 : if( getRowIndex(irow,i)==jcol ) return data[irow*ncols+i];
266 : }
267 : return 0.0;
268 : }
269 540725624 : plumed_massert( ival<data.size(), "cannot get value from " + name );
270 540725624 : return data[ival];
271 : }
272 :
273 39305298 : void Value::addForce(const std::size_t& iforce, double f, const bool trueind) {
274 39305298 : hasForce=true;
275 39305298 : if( shape.size()==2 && !hasDeriv && ncols<shape[1] && trueind ) {
276 0 : unsigned irow = std::floor( iforce / shape[0] ), jcol = iforce%shape[0];
277 0 : for(unsigned i=0; i<getRowLength(irow); ++i) {
278 0 : if( getRowIndex(irow,i)==jcol ) { inputForce[irow*ncols+i]+=f; return; }
279 : }
280 0 : plumed_assert( fabs(f)<epsilon ); return;
281 : }
282 39305298 : plumed_massert( iforce<inputForce.size(), "can't add force to " + name );
283 39305298 : inputForce[iforce]+=f;
284 : }
285 :
286 :
287 12065 : void Value::buildDataStore( const bool forprint ) {
288 12065 : if( getRank()==0 ) return;
289 5289 : storedata=true; setShape( shape );
290 5289 : if( !forprint ) return ;
291 192 : ActionWithVector* av=dynamic_cast<ActionWithVector*>( action );
292 192 : if( av ) (av->getFirstActionInChain())->never_reduce_tasks=true;
293 : }
294 :
295 7082 : void Value::reshapeMatrixStore( const unsigned& n ) {
296 : plumed_dbg_assert( shape.size()==2 && !hasDeriv );
297 7082 : if( !storedata ) return ;
298 7082 : ncols=n; if( ncols>shape[1] ) ncols=shape[1];
299 7082 : unsigned size=shape[0]*ncols;
300 7082 : if( matrix_bookeeping.size()!=(size+shape[0]) ) {
301 2562 : data.resize( size ); inputForce.resize( size );
302 2562 : matrix_bookeeping.resize( size + shape[0], 0 );
303 2562 : if( ncols>=shape[1] ) {
304 248548 : for(unsigned i=0; i<shape[0]; ++i) {
305 246049 : matrix_bookeeping[(1+ncols)*i] = shape[1];
306 23653127 : for(unsigned j=0; j<shape[1]; ++j) matrix_bookeeping[(1+ncols)*i+1+j]=j;
307 : }
308 : }
309 : }
310 7082 : if( ncols<shape[1] ) std::fill(matrix_bookeeping.begin(), matrix_bookeeping.end(), 0);
311 : }
312 :
313 35389 : void Value::setPositionInMatrixStash( const unsigned& p ) {
314 : plumed_dbg_assert( shape.size()==2 && !hasDeriv );
315 35389 : matpos=p;
316 35389 : }
317 :
318 15992076 : bool Value::ignoreStoredValue(const std::string& c) const {
319 15992076 : if( !storedata && shape.size()>0 ) return true;
320 6832518 : ActionWithVector* av=dynamic_cast<ActionWithVector*>(action);
321 6832518 : if( av ) return (av->getFirstActionInChain())->getLabel()==c;
322 : return false;
323 : }
324 :
325 5510 : void Value::setConstant() {
326 5510 : valtype=constant; storedata=true; setShape( shape );
327 5510 : if( getRank()==2 && !hasDeriv ) reshapeMatrixStore( shape[1] );
328 5510 : }
329 :
330 456 : void Value::writeBinary(std::ostream&o) const {
331 456 : o.write(reinterpret_cast<const char*>(&data[0]),data.size()*sizeof(double));
332 456 : }
333 :
334 1274 : void Value::setSymmetric( const bool& sym ) {
335 1274 : plumed_assert( shape.size()==2 && !hasDeriv );
336 1274 : if( sym && shape[0]!=shape[1] ) plumed_merror("non-square matrix cannot be symmetric");
337 1274 : symmetric=sym;
338 1274 : }
339 :
340 2220 : void Value::retrieveEdgeList( unsigned& nedge, std::vector<std::pair<unsigned,unsigned> >& active, std::vector<double>& elems ) {
341 2220 : nedge=0; plumed_dbg_assert( shape.size()==2 && !hasDeriv );
342 : // Check we have enough space to store the edge list
343 2220 : if( elems.size()<shape[0]*ncols ) { elems.resize( shape[0]*ncols ); active.resize( shape[0]*ncols ); }
344 :
345 269579 : for(unsigned i=0; i<shape[0]; ++i) {
346 : unsigned ncol = getRowLength(i);
347 15002116 : for(unsigned j=0; j<ncol; ++j) {
348 14734757 : if( fabs(get(i*ncols+j,false))<epsilon ) continue;
349 3265037 : if( symmetric && getRowIndex(i,j)>i ) continue;
350 3213497 : active[nedge].first = i; active[nedge].second = getRowIndex(i,j);
351 3213497 : elems[nedge] = get(i*ncols+j,false); nedge++;
352 : }
353 : }
354 2220 : }
355 :
356 456 : void Value::readBinary(std::istream&i) {
357 456 : i.read(reinterpret_cast<char*>(&data[0]),data.size()*sizeof(double));
358 456 : }
359 :
360 775666 : void Value::convertIndexToindices(const std::size_t& index, std::vector<unsigned>& indices ) const {
361 775666 : if( hasDeriv || getRank()==1 ) {
362 12302 : std::size_t kk=index; indices[0]=index%shape[0];
363 12302 : for(unsigned i=1; i<shape.size()-1; ++i) {
364 0 : kk=(kk-indices[i-1])/shape[i-1];
365 0 : indices[i]=kk%shape[i];
366 : }
367 12302 : if(shape.size()>=2) indices[shape.size()-1]=(kk-indices[shape.size()-2])/shape[shape.size()-2];
368 763364 : } else if( getRank()==2 ) {
369 763364 : indices[0]=std::floor( index/shape[1] ); indices[1] = index%shape[1];
370 : }
371 775666 : }
372 :
373 548216 : void Value::print( OFile& ofile ) const {
374 571290 : if( isPeriodic() ) { ofile.printField( "min_" + name, str_min ); ofile.printField("max_" + name, str_max ); }
375 548216 : if( shape.size()==0 || getNumberOfValues()==1 ) {
376 546897 : ofile.printField( name, get(0) );
377 : } else {
378 1319 : std::vector<unsigned> indices( shape.size() );
379 776985 : for(unsigned i=0; i<getNumberOfValues(); ++i) {
380 775666 : convertIndexToindices( i, indices ); std::string num, fname = name;
381 2314696 : for(unsigned i=0; i<shape.size(); ++i) { Tools::convert( indices[i]+1, num ); fname += "." + num; }
382 775666 : ofile.printField( fname, get(i) );
383 : }
384 : }
385 548216 : }
386 :
387 232611 : unsigned Value::getGoodNumThreads( const unsigned& j, const unsigned& k ) const {
388 232611 : return OpenMP::getGoodNumThreads( &data[j], (k-j) );
389 : }
390 :
391 7526197 : bool Value::calculateOnUpdate() const {
392 7526197 : return (valtype==average || valtype==calcFromAverage);
393 : }
394 :
395 132 : std::string Value::getValueType() const {
396 132 : if( getRank()==0 ) return "scalar";
397 99 : if( getRank()>0 && hasDerivatives() ) return "grid";
398 94 : if( getRank()==1 ) return "vector";
399 21 : if( getRank()==2 ) return "matrix";
400 0 : plumed_merror("unknown type for value " + getName() );
401 : return "";
402 : }
403 :
404 : }
|
