1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2011-2021 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
23 #include "Grid.h"
24 #include "Tools.h"
25 #include "core/Value.h"
26 #include "File.h"
27 #include "Exception.h"
28 #include "KernelFunctions.h"
29 #include "RootFindingBase.h"
30 #include "Communicator.h"
31
32 #include <vector>
33 #include <cmath>
34 #include <iostream>
35 #include <sstream>
36 #include <cstdio>
37 #include <cfloat>
38 #include <array>
39
40 using namespace std;
41 namespace PLMD {
42
43 constexpr size_t GridBase::maxdim;
44
GridBase(const std::string & funcl,const std::vector<Value * > & args,const vector<std::string> & gmin,const vector<std::string> & gmax,const vector<unsigned> & nbin,bool dospline,bool usederiv)45 GridBase::GridBase(const std::string& funcl, const std::vector<Value*> & args, const vector<std::string> & gmin,
46 const vector<std::string> & gmax, const vector<unsigned> & nbin, bool dospline, bool usederiv) {
47 // various checks
48 plumed_assert(args.size()<=maxdim) << "grid dim cannot exceed "<<maxdim;
49 plumed_massert(args.size()==gmin.size(),"grid min dimensions in input do not match number of arguments");
50 plumed_massert(args.size()==nbin.size(),"number of bins on input do not match number of arguments");
51 plumed_massert(args.size()==gmax.size(),"grid max dimensions in input do not match number of arguments");
52 unsigned dim=gmax.size();
53 std::vector<std::string> names;
54 std::vector<bool> isperiodic;
55 std::vector<string> pmin,pmax;
56 names.resize( dim );
57 isperiodic.resize( dim );
58 pmin.resize( dim );
59 pmax.resize( dim );
60 for(unsigned int i=0; i<dim; ++i) {
61 names[i]=args[i]->getName();
62 if( args[i]->isPeriodic() ) {
63 isperiodic[i]=true;
64 args[i]->getDomain( pmin[i], pmax[i] );
65 } else {
66 isperiodic[i]=false;
67 pmin[i]="0.";
68 pmax[i]="0.";
69 }
70 }
71 // this is a value-independent initializator
72 Init(funcl,names,gmin,gmax,nbin,dospline,usederiv,isperiodic,pmin,pmax);
73 }
74
GridBase(const std::string & funcl,const std::vector<string> & names,const std::vector<std::string> & gmin,const vector<std::string> & gmax,const std::vector<unsigned> & nbin,bool dospline,bool usederiv,const std::vector<bool> & isperiodic,const std::vector<std::string> & pmin,const std::vector<std::string> & pmax)75 GridBase::GridBase(const std::string& funcl, const std::vector<string> &names, const std::vector<std::string> & gmin,
76 const vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline, bool usederiv, const std::vector<bool> &isperiodic, const std::vector<std::string> &pmin, const std::vector<std::string> &pmax ) {
77 // this calls the initializator
78 Init(funcl,names,gmin,gmax,nbin,dospline,usederiv,isperiodic,pmin,pmax);
79 }
80
Init(const std::string & funcl,const std::vector<std::string> & names,const vector<std::string> & gmin,const std::vector<std::string> & gmax,const std::vector<unsigned> & nbin,bool dospline,bool usederiv,const std::vector<bool> & isperiodic,const std::vector<std::string> & pmin,const std::vector<std::string> & pmax)81 void GridBase::Init(const std::string& funcl, const std::vector<std::string> &names, const vector<std::string> & gmin,
82 const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline, bool usederiv,
83 const std::vector<bool> &isperiodic, const std::vector<std::string> &pmin, const std::vector<std::string> &pmax ) {
84 fmt_="%14.9f";
85 // various checks
86 plumed_assert(names.size()<=maxdim) << "grid size cannot exceed "<<maxdim;
87 plumed_massert(names.size()==gmin.size(),"grid dimensions in input do not match number of arguments");
88 plumed_massert(names.size()==nbin.size(),"grid dimensions in input do not match number of arguments");
89 plumed_massert(names.size()==gmax.size(),"grid dimensions in input do not match number of arguments");
90 dimension_=gmax.size();
91 str_min_=gmin; str_max_=gmax;
92 argnames.resize( dimension_ );
93 min_.resize( dimension_ );
94 max_.resize( dimension_ );
95 pbc_.resize( dimension_ );
96 for(unsigned int i=0; i<dimension_; ++i) {
97 argnames[i]=names[i];
98 if( isperiodic[i] ) {
99 pbc_[i]=true;
100 str_min_[i]=pmin[i];
101 str_max_[i]=pmax[i];
102 } else {
103 pbc_[i]=false;
104 }
105 Tools::convert(str_min_[i],min_[i]);
106 Tools::convert(str_max_[i],max_[i]);
107 funcname=funcl;
108 plumed_massert(max_[i]>min_[i],"maximum in grid must be larger than minimum");
109 plumed_massert(nbin[i]>0,"number of grid points must be greater than zero");
110 }
111 nbin_=nbin;
112 dospline_=dospline;
113 usederiv_=usederiv;
114 if(dospline_) plumed_assert(dospline_==usederiv_);
115 maxsize_=1;
116 for(unsigned int i=0; i<dimension_; ++i) {
117 dx_.push_back( (max_[i]-min_[i])/static_cast<double>( nbin_[i] ) );
118 if( !pbc_[i] ) { max_[i] += dx_[i]; nbin_[i] += 1; }
119 maxsize_*=nbin_[i];
120 }
121 }
122
getMin() const123 vector<std::string> GridBase::getMin() const {
124 return str_min_;
125 }
126
getMax() const127 vector<std::string> GridBase::getMax() const {
128 return str_max_;
129 }
130
getDx() const131 vector<double> GridBase::getDx() const {
132 return dx_;
133 }
134
getDx(index_t j) const135 double GridBase::getDx(index_t j) const {
136 return dx_[j];
137 }
138
getBinVolume() const139 double GridBase::getBinVolume() const {
140 double vol=1.;
141 for(unsigned i=0; i<dx_.size(); ++i) vol*=dx_[i];
142 return vol;
143 }
144
getIsPeriodic() const145 vector<bool> GridBase::getIsPeriodic() const {
146 return pbc_;
147 }
148
getNbin() const149 vector<unsigned> GridBase::getNbin() const {
150 return nbin_;
151 }
152
getArgNames() const153 vector<string> GridBase::getArgNames() const {
154 return argnames;
155 }
156
157
getDimension() const158 unsigned GridBase::getDimension() const {
159 return dimension_;
160 }
161
162 // we are flattening arrays using a column-major order
getIndex(const vector<unsigned> & indices) const163 GridBase::index_t GridBase::getIndex(const vector<unsigned> & indices) const {
164 plumed_dbg_assert(indices.size()==dimension_);
165 for(unsigned int i=0; i<dimension_; i++)
166 if(indices[i]>=nbin_[i]) {
167 std::string is;
168 Tools::convert(i,is);
169 std::string msg="ERROR: the system is looking for a value outside the grid along the " + is + " ("+getArgNames()[i]+")";
170 plumed_merror(msg+" index!");
171 }
172 index_t index=indices[dimension_-1];
173 for(unsigned int i=dimension_-1; i>0; --i) {
174 index=index*nbin_[i-1]+indices[i-1];
175 }
176 return index;
177 }
178
getIndex(const vector<double> & x) const179 GridBase::index_t GridBase::getIndex(const vector<double> & x) const {
180 plumed_dbg_assert(x.size()==dimension_);
181 return getIndex(getIndices(x));
182 }
183
184 // we are flattening arrays using a column-major order
getIndices(index_t index) const185 vector<unsigned> GridBase::getIndices(index_t index) const {
186 vector<unsigned> indices(dimension_);
187 index_t kk=index;
188 indices[0]=(index%nbin_[0]);
189 for(unsigned int i=1; i<dimension_-1; ++i) {
190 kk=(kk-indices[i-1])/nbin_[i-1];
191 indices[i]=(kk%nbin_[i]);
192 }
193 if(dimension_>=2) {
194 indices[dimension_-1]=((kk-indices[dimension_-2])/nbin_[dimension_-2]);
195 }
196 return indices;
197 }
198
getIndices(index_t index,std::vector<unsigned> & indices) const199 void GridBase::getIndices(index_t index, std::vector<unsigned>& indices) const {
200 if (indices.size()!=dimension_) indices.resize(dimension_);
201 index_t kk=index;
202 indices[0]=(index%nbin_[0]);
203 for(unsigned int i=1; i<dimension_-1; ++i) {
204 kk=(kk-indices[i-1])/nbin_[i-1];
205 indices[i]=(kk%nbin_[i]);
206 }
207 if(dimension_>=2) {
208 indices[dimension_-1]=((kk-indices[dimension_-2])/nbin_[dimension_-2]);
209 }
210 }
211
getIndices(const vector<double> & x) const212 vector<unsigned> GridBase::getIndices(const vector<double> & x) const {
213 plumed_dbg_assert(x.size()==dimension_);
214 vector<unsigned> indices(dimension_);
215 for(unsigned int i=0; i<dimension_; ++i) {
216 indices[i] = unsigned(floor((x[i]-min_[i])/dx_[i]));
217 }
218 return indices;
219 }
220
getIndices(const vector<double> & x,std::vector<unsigned> & indices) const221 void GridBase::getIndices(const vector<double> & x, std::vector<unsigned>& indices) const {
222 plumed_dbg_assert(x.size()==dimension_);
223 if (indices.size()!=dimension_) indices.resize(dimension_);
224 for(unsigned int i=0; i<dimension_; ++i) {
225 indices[i] = unsigned(floor((x[i]-min_[i])/dx_[i]));
226 }
227 }
228
getPoint(const vector<unsigned> & indices) const229 vector<double> GridBase::getPoint(const vector<unsigned> & indices) const {
230 plumed_dbg_assert(indices.size()==dimension_);
231 vector<double> x(dimension_);
232 for(unsigned int i=0; i<dimension_; ++i) {
233 x[i]=min_[i]+(double)(indices[i])*dx_[i];
234 }
235 return x;
236 }
237
getPoint(index_t index) const238 vector<double> GridBase::getPoint(index_t index) const {
239 plumed_dbg_assert(index<maxsize_);
240 return getPoint(getIndices(index));
241 }
242
getPoint(const vector<double> & x) const243 vector<double> GridBase::getPoint(const vector<double> & x) const {
244 plumed_dbg_assert(x.size()==dimension_);
245 return getPoint(getIndices(x));
246 }
247
getPoint(index_t index,std::vector<double> & point) const248 void GridBase::getPoint(index_t index,std::vector<double> & point) const {
249 plumed_dbg_assert(index<maxsize_);
250 getPoint(getIndices(index),point);
251 }
252
getPoint(const std::vector<unsigned> & indices,std::vector<double> & point) const253 void GridBase::getPoint(const std::vector<unsigned> & indices,std::vector<double> & point) const {
254 plumed_dbg_assert(indices.size()==dimension_);
255 plumed_dbg_assert(point.size()==dimension_);
256 for(unsigned int i=0; i<dimension_; ++i) {
257 point[i]=min_[i]+(double)(indices[i])*dx_[i];
258 }
259 }
260
getPoint(const std::vector<double> & x,std::vector<double> & point) const261 void GridBase::getPoint(const std::vector<double> & x,std::vector<double> & point) const {
262 plumed_dbg_assert(x.size()==dimension_);
263 getPoint(getIndices(x),point);
264 }
265
getNeighbors(const vector<unsigned> & indices,const vector<unsigned> & nneigh) const266 vector<GridBase::index_t> GridBase::getNeighbors
267 (const vector<unsigned> &indices,const vector<unsigned> &nneigh)const {
268 plumed_dbg_assert(indices.size()==dimension_ && nneigh.size()==dimension_);
269
270 vector<index_t> neighbors;
271 vector<unsigned> small_bin(dimension_);
272
273 unsigned small_nbin=1;
274 for(unsigned j=0; j<dimension_; ++j) {
275 small_bin[j]=(2*nneigh[j]+1);
276 small_nbin*=small_bin[j];
277 }
278
279 vector<unsigned> small_indices(dimension_);
280 vector<unsigned> tmp_indices;
281 for(unsigned index=0; index<small_nbin; ++index) {
282 tmp_indices.resize(dimension_);
283 unsigned kk=index;
284 small_indices[0]=(index%small_bin[0]);
285 for(unsigned i=1; i<dimension_-1; ++i) {
286 kk=(kk-small_indices[i-1])/small_bin[i-1];
287 small_indices[i]=(kk%small_bin[i]);
288 }
289 if(dimension_>=2) {
290 small_indices[dimension_-1]=((kk-small_indices[dimension_-2])/small_bin[dimension_-2]);
291 }
292 unsigned ll=0;
293 for(unsigned i=0; i<dimension_; ++i) {
294 int i0=small_indices[i]-nneigh[i]+indices[i];
295 if(!pbc_[i] && i0<0) continue;
296 if(!pbc_[i] && i0>=static_cast<int>(nbin_[i])) continue;
297 if( pbc_[i] && i0<0) i0=nbin_[i]-(-i0)%nbin_[i];
298 if( pbc_[i] && i0>=static_cast<int>(nbin_[i])) i0%=nbin_[i];
299 tmp_indices[ll]=static_cast<unsigned>(i0);
300 ll++;
301 }
302 tmp_indices.resize(ll);
303 if(tmp_indices.size()==dimension_) {neighbors.push_back(getIndex(tmp_indices));}
304 }
305 return neighbors;
306 }
307
getNeighbors(const vector<double> & x,const vector<unsigned> & nneigh) const308 vector<GridBase::index_t> GridBase::getNeighbors
309 (const vector<double> & x,const vector<unsigned> & nneigh)const {
310 plumed_dbg_assert(x.size()==dimension_ && nneigh.size()==dimension_);
311 return getNeighbors(getIndices(x),nneigh);
312 }
313
getNeighbors(index_t index,const vector<unsigned> & nneigh) const314 vector<GridBase::index_t> GridBase::getNeighbors
315 (index_t index,const vector<unsigned> & nneigh)const {
316 plumed_dbg_assert(index<maxsize_ && nneigh.size()==dimension_);
317 return getNeighbors(getIndices(index),nneigh);
318 }
319
getSplineNeighbors(const vector<unsigned> & indices,vector<GridBase::index_t> & neighbors,unsigned & nneighbors) const320 void GridBase::getSplineNeighbors(const vector<unsigned> & indices, vector<GridBase::index_t>& neighbors, unsigned& nneighbors)const {
321 plumed_dbg_assert(indices.size()==dimension_);
322 unsigned nneigh=unsigned(pow(2.0,int(dimension_)));
323 if (neighbors.size()!=nneigh) neighbors.resize(nneigh);
324
325 vector<unsigned> nindices(dimension_);
326 unsigned inind; nneighbors = 0;
327 for(unsigned int i=0; i<nneigh; ++i) {
328 unsigned tmp=i; inind=0;
329 for(unsigned int j=0; j<dimension_; ++j) {
330 unsigned i0=tmp%2+indices[j];
331 tmp/=2;
332 if(!pbc_[j] && i0==nbin_[j]) continue;
333 if( pbc_[j] && i0==nbin_[j]) i0=0;
334 nindices[inind++]=i0;
335 }
336 if(inind==dimension_) neighbors[nneighbors++]=getIndex(nindices);
337 }
338 }
339
getNearestNeighbors(const index_t index) const340 vector<GridBase::index_t> GridBase::getNearestNeighbors(const index_t index) const {
341 vector<index_t> nearest_neighs = vector<index_t>();
342 for (unsigned i = 0; i < dimension_; i++) {
343 vector<unsigned> neighsneeded = vector<unsigned>(dimension_, 0);
344 neighsneeded[i] = 1;
345 vector<index_t> singledim_nearest_neighs = getNeighbors(index, neighsneeded);
346 for (unsigned j = 0; j < singledim_nearest_neighs.size(); j++) {
347 index_t neigh = singledim_nearest_neighs[j];
348 if (neigh != index) {
349 nearest_neighs.push_back(neigh);
350 }
351 }
352 }
353 return nearest_neighs;
354 }
355
getNearestNeighbors(const vector<unsigned> & indices) const356 vector<GridBase::index_t> GridBase::getNearestNeighbors(const vector<unsigned> &indices) const {
357 plumed_dbg_assert(indices.size() == dimension_);
358 return getNearestNeighbors(getIndex(indices));
359 }
360
361
addKernel(const KernelFunctions & kernel)362 void GridBase::addKernel( const KernelFunctions& kernel ) {
363 plumed_dbg_assert( kernel.ndim()==dimension_ );
364 std::vector<unsigned> nneighb=kernel.getSupport( dx_ );
365 std::vector<index_t> neighbors=getNeighbors( kernel.getCenter(), nneighb );
366 std::vector<double> xx( dimension_ );
367 std::vector<std::unique_ptr<Value>> vv( dimension_ );
368 std::string str_min, str_max;
369 for(unsigned i=0; i<dimension_; ++i) {
370 vv[i].reset(new Value());
371 if( pbc_[i] ) {
372 Tools::convert(min_[i],str_min);
373 Tools::convert(max_[i],str_max);
374 vv[i]->setDomain( str_min, str_max );
375 } else {
376 vv[i]->setNotPeriodic();
377 }
378 }
379
380 // vv_ptr contains plain pointers obtained from vv.
381 // this is the simplest way to replace a unique_ptr here.
382 // perhaps the interface of kernel.evaluate() should be changed
383 // in order to accept a std::vector<std::unique_ptr<Value>>
384 auto vv_ptr=Tools::unique2raw(vv);
385
386 std::vector<double> der( dimension_ );
387 for(unsigned i=0; i<neighbors.size(); ++i) {
388 index_t ineigh=neighbors[i];
389 getPoint( ineigh, xx );
390 for(unsigned j=0; j<dimension_; ++j) vv[j]->set(xx[j]);
391 double newval = kernel.evaluate( vv_ptr, der, usederiv_ );
392 if( usederiv_ ) addValueAndDerivatives( ineigh, newval, der );
393 else addValue( ineigh, newval );
394 }
395 }
396
397
getValue(const vector<unsigned> & indices) const398 double GridBase::getValue(const vector<unsigned> & indices) const {
399 return getValue(getIndex(indices));
400 }
401
getValue(const vector<double> & x) const402 double GridBase::getValue(const vector<double> & x) const {
403 if(!dospline_) {
404 return getValue(getIndex(x));
405 } else {
406 vector<double> der(dimension_);
407 return getValueAndDerivatives(x,der);
408 }
409 }
410
getValueAndDerivatives(const vector<unsigned> & indices,vector<double> & der) const411 double GridBase::getValueAndDerivatives
412 (const vector<unsigned> & indices, vector<double>& der) const {
413 return getValueAndDerivatives(getIndex(indices),der);
414 }
415
getValueAndDerivatives(const vector<double> & x,vector<double> & der) const416 double GridBase::getValueAndDerivatives
417 (const vector<double> & x, vector<double>& der) const {
418 plumed_dbg_assert(der.size()==dimension_ && usederiv_);
419
420 if(dospline_) {
421 double X,X2,X3,value;
422 std::array<double,maxdim> fd, C, D;
423 std::vector<double> dder(dimension_);
424 // reset
425 value=0.0;
426 for(unsigned int i=0; i<dimension_; ++i) der[i]=0.0;
427
428 vector<unsigned> indices(dimension_);
429 getIndices(x, indices);
430 vector<double> xfloor(dimension_);
431 getPoint(indices, xfloor);
432 vector<index_t> neigh; unsigned nneigh; getSplineNeighbors(indices, neigh, nneigh);
433
434 // loop over neighbors
435 vector<unsigned> nindices;
436 for(unsigned int ipoint=0; ipoint<nneigh; ++ipoint) {
437 double grid=getValueAndDerivatives(neigh[ipoint],dder);
438 getIndices(neigh[ipoint], nindices);
439 double ff=1.0;
440
441 for(unsigned j=0; j<dimension_; ++j) {
442 int x0=1;
443 if(nindices[j]==indices[j]) x0=0;
444 double dx=getDx(j);
445 X=fabs((x[j]-xfloor[j])/dx-(double)x0);
446 X2=X*X;
447 X3=X2*X;
448 double yy;
449 if(fabs(grid)<0.0000001) yy=0.0;
450 else yy=-dder[j]/grid;
451 C[j]=(1.0-3.0*X2+2.0*X3) - (x0?-1.0:1.0)*yy*(X-2.0*X2+X3)*dx;
452 D[j]=( -6.0*X +6.0*X2) - (x0?-1.0:1.0)*yy*(1.0-4.0*X +3.0*X2)*dx;
453 D[j]*=(x0?-1.0:1.0)/dx;
454 ff*=C[j];
455 }
456 for(unsigned j=0; j<dimension_; ++j) {
457 fd[j]=D[j];
458 for(unsigned i=0; i<dimension_; ++i) if(i!=j) fd[j]*=C[i];
459 }
460 value+=grid*ff;
461 for(unsigned j=0; j<dimension_; ++j) der[j]+=grid*fd[j];
462 }
463 return value;
464 } else {
465 return getValueAndDerivatives(getIndex(x),der);
466 }
467 }
468
setValue(const vector<unsigned> & indices,double value)469 void GridBase::setValue(const vector<unsigned> & indices, double value) {
470 setValue(getIndex(indices),value);
471 }
472
setValueAndDerivatives(const vector<unsigned> & indices,double value,vector<double> & der)473 void GridBase::setValueAndDerivatives
474 (const vector<unsigned> & indices, double value, vector<double>& der) {
475 setValueAndDerivatives(getIndex(indices),value,der);
476 }
477
addValue(const vector<unsigned> & indices,double value)478 void GridBase::addValue(const vector<unsigned> & indices, double value) {
479 addValue(getIndex(indices),value);
480 }
481
addValueAndDerivatives(const vector<unsigned> & indices,double value,vector<double> & der)482 void GridBase::addValueAndDerivatives
483 (const vector<unsigned> & indices, double value, vector<double>& der) {
484 addValueAndDerivatives(getIndex(indices),value,der);
485 }
486
writeHeader(OFile & ofile)487 void GridBase::writeHeader(OFile& ofile) {
488 for(unsigned i=0; i<dimension_; ++i) {
489 ofile.addConstantField("min_" + argnames[i]);
490 ofile.addConstantField("max_" + argnames[i]);
491 ofile.addConstantField("nbins_" + argnames[i]);
492 ofile.addConstantField("periodic_" + argnames[i]);
493 }
494 }
495
clear()496 void Grid::clear() {
497 grid_.assign(maxsize_,0.0);
498 if(usederiv_) der_.assign(maxsize_*dimension_,0.0);
499 }
500
writeToFile(OFile & ofile)501 void Grid::writeToFile(OFile& ofile) {
502 vector<double> xx(dimension_);
503 vector<double> der(dimension_);
504 double f;
505 writeHeader(ofile);
506 for(index_t i=0; i<getSize(); ++i) {
507 xx=getPoint(i);
508 if(usederiv_) {f=getValueAndDerivatives(i,der);}
509 else {f=getValue(i);}
510 if(i>0 && dimension_>1 && getIndices(i)[dimension_-2]==0) ofile.printf("\n");
511 for(unsigned j=0; j<dimension_; ++j) {
512 ofile.printField("min_" + argnames[j], str_min_[j] );
513 ofile.printField("max_" + argnames[j], str_max_[j] );
514 ofile.printField("nbins_" + argnames[j], static_cast<int>(nbin_[j]) );
515 if( pbc_[j] ) ofile.printField("periodic_" + argnames[j], "true" );
516 else ofile.printField("periodic_" + argnames[j], "false" );
517 }
518 for(unsigned j=0; j<dimension_; ++j) { ofile.fmtField(" "+fmt_); ofile.printField(argnames[j],xx[j]); }
519 ofile.fmtField(" "+fmt_); ofile.printField(funcname,f);
520 if(usederiv_) for(unsigned j=0; j<dimension_; ++j) { ofile.fmtField(" "+fmt_); ofile.printField("der_" + argnames[j],der[j]); }
521 ofile.printField();
522 }
523 }
524
writeCubeFile(OFile & ofile,const double & lunit)525 void GridBase::writeCubeFile(OFile& ofile, const double& lunit) {
526 plumed_assert( dimension_==3 );
527 ofile.printf("PLUMED CUBE FILE\n");
528 ofile.printf("OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z\n");
529 // Number of atoms followed by position of origin (origin set so that center of grid is in center of cell)
530 ofile.printf("%d %f %f %f\n",1,-0.5*lunit*(max_[0]-min_[0]),-0.5*lunit*(max_[1]-min_[1]),-0.5*lunit*(max_[2]-min_[2]));
531 ofile.printf("%u %f %f %f\n",nbin_[0],lunit*dx_[0],0.0,0.0); // Number of bins in each direction followed by
532 ofile.printf("%u %f %f %f\n",nbin_[1],0.0,lunit*dx_[1],0.0); // shape of voxel
533 ofile.printf("%u %f %f %f\n",nbin_[2],0.0,0.0,lunit*dx_[2]);
534 ofile.printf("%d %f %f %f\n",1,0.0,0.0,0.0); // Fake atom otherwise VMD doesn't work
535 std::vector<unsigned> pp(3);
536 for(pp[0]=0; pp[0]<nbin_[0]; ++pp[0]) {
537 for(pp[1]=0; pp[1]<nbin_[1]; ++pp[1]) {
538 for(pp[2]=0; pp[2]<nbin_[2]; ++pp[2]) {
539 ofile.printf("%f ",getValue(pp) );
540 if(pp[2]%6==5) ofile.printf("\n");
541 }
542 ofile.printf("\n");
543 }
544 }
545 }
546
create(const std::string & funcl,const std::vector<Value * > & args,IFile & ifile,const vector<std::string> & gmin,const vector<std::string> & gmax,const vector<unsigned> & nbin,bool dosparse,bool dospline,bool doder)547 std::unique_ptr<GridBase> GridBase::create(const std::string& funcl, const std::vector<Value*> & args, IFile& ifile,
548 const vector<std::string> & gmin,const vector<std::string> & gmax,
549 const vector<unsigned> & nbin,bool dosparse, bool dospline, bool doder) {
550 std::unique_ptr<GridBase> grid=GridBase::create(funcl,args,ifile,dosparse,dospline,doder);
551 std::vector<unsigned> cbin( grid->getNbin() );
552 std::vector<std::string> cmin( grid->getMin() ), cmax( grid->getMax() );
553 for(unsigned i=0; i<args.size(); ++i) {
554 plumed_massert( cmin[i]==gmin[i], "mismatched grid min" );
555 plumed_massert( cmax[i]==gmax[i], "mismatched grid max" );
556 if( args[i]->isPeriodic() ) {
557 plumed_massert( cbin[i]==nbin[i], "mismatched grid nbins" );
558 } else {
559 plumed_massert( (cbin[i]-1)==nbin[i], "mismatched grid nbins");
560 }
561 }
562 return grid;
563 }
564
create(const std::string & funcl,const std::vector<Value * > & args,IFile & ifile,bool dosparse,bool dospline,bool doder)565 std::unique_ptr<GridBase> GridBase::create(const std::string& funcl, const std::vector<Value*> & args, IFile& ifile, bool dosparse, bool dospline, bool doder)
566 {
567 std::unique_ptr<GridBase> grid;
568 unsigned nvar=args.size(); bool hasder=false; std::string pstring;
569 std::vector<int> gbin1(nvar); std::vector<unsigned> gbin(nvar);
570 std::vector<std::string> labels(nvar),gmin(nvar),gmax(nvar);
571 std::vector<std::string> fieldnames; ifile.scanFieldList( fieldnames );
572 // Retrieve names for fields
573 for(unsigned i=0; i<args.size(); ++i) labels[i]=args[i]->getName();
574 // And read the stuff from the header
575 plumed_massert( ifile.FieldExist( funcl ), "no column labelled " + funcl + " in in grid input");
576 for(unsigned i=0; i<args.size(); ++i) {
577 ifile.scanField( "min_" + labels[i], gmin[i]);
578 ifile.scanField( "max_" + labels[i], gmax[i]);
579 ifile.scanField( "periodic_" + labels[i], pstring );
580 ifile.scanField( "nbins_" + labels[i], gbin1[i]);
581 plumed_assert( gbin1[i]>0 );
582 if( args[i]->isPeriodic() ) {
583 plumed_massert( pstring=="true", "input value is periodic but grid is not");
584 std::string pmin, pmax;
585 args[i]->getDomain( pmin, pmax ); gbin[i]=gbin1[i];
586 if( pmin!=gmin[i] || pmax!=gmax[i] ) plumed_merror("mismatch between grid boundaries and periods of values");
587 } else {
588 gbin[i]=gbin1[i]-1; // Note header in grid file indicates one more bin that there should be when data is not periodic
589 plumed_massert( pstring=="false", "input value is not periodic but grid is");
590 }
591 hasder=ifile.FieldExist( "der_" + args[i]->getName() );
592 if( doder && !hasder ) plumed_merror("missing derivatives from grid file");
593 for(unsigned j=0; j<fieldnames.size(); ++j) {
594 for(unsigned k=i+1; k<args.size(); ++k) {
595 if( fieldnames[j]==labels[k] ) plumed_merror("arguments in input are not in same order as in grid file");
596 }
597 if( fieldnames[j]==labels[i] ) break;
598 }
599 }
600
601 if(!dosparse) {grid.reset(new Grid(funcl,args,gmin,gmax,gbin,dospline,doder));}
602 else {grid.reset(new SparseGrid(funcl,args,gmin,gmax,gbin,dospline,doder));}
603
604 vector<double> xx(nvar),dder(nvar);
605 vector<double> dx=grid->getDx();
606 double f,x;
607 while( ifile.scanField(funcl,f) ) {
608 for(unsigned i=0; i<nvar; ++i) {
609 ifile.scanField(labels[i],x); xx[i]=x+dx[i]/2.0;
610 ifile.scanField( "min_" + labels[i], gmin[i]);
611 ifile.scanField( "max_" + labels[i], gmax[i]);
612 ifile.scanField( "nbins_" + labels[i], gbin1[i]);
613 ifile.scanField( "periodic_" + labels[i], pstring );
614 }
615 if(hasder) { for(unsigned i=0; i<nvar; ++i) { ifile.scanField( "der_" + args[i]->getName(), dder[i] ); } }
616 index_t index=grid->getIndex(xx);
617 if(doder) {grid->setValueAndDerivatives(index,f,dder);}
618 else {grid->setValue(index,f);}
619 ifile.scanField();
620 }
621 return grid;
622 }
623
getMinValue() const624 double Grid::getMinValue() const {
625 double minval;
626 minval=DBL_MAX;
627 for(index_t i=0; i<grid_.size(); ++i) {
628 if(grid_[i]<minval)minval=grid_[i];
629 }
630 return minval;
631 }
632
getMaxValue() const633 double Grid::getMaxValue() const {
634 double maxval;
635 maxval=DBL_MIN;
636 for(index_t i=0; i<grid_.size(); ++i) {
637 if(grid_[i]>maxval)maxval=grid_[i];
638 }
639 return maxval;
640 }
641
scaleAllValuesAndDerivatives(const double & scalef)642 void Grid::scaleAllValuesAndDerivatives( const double& scalef ) {
643 if(usederiv_) {
644 for(index_t i=0; i<grid_.size(); ++i) {
645 grid_[i]*=scalef;
646 for(unsigned j=0; j<dimension_; ++j) der_[i*dimension_+j]*=scalef;
647 }
648 } else {
649 for(index_t i=0; i<grid_.size(); ++i) grid_[i]*=scalef;
650 }
651 }
652
logAllValuesAndDerivatives(const double & scalef)653 void Grid::logAllValuesAndDerivatives( const double& scalef ) {
654 if(usederiv_) {
655 for(index_t i=0; i<grid_.size(); ++i) {
656 grid_[i] = scalef*log(grid_[i]);
657 for(unsigned j=0; j<dimension_; ++j) der_[i*dimension_+j] = scalef/der_[i*dimension_+j];
658 }
659 } else {
660 for(index_t i=0; i<grid_.size(); ++i) grid_[i] = scalef*log(grid_[i]);
661 }
662 }
663
setMinToZero()664 void Grid::setMinToZero() {
665 double min=grid_[0];
666 for(index_t i=1; i<grid_.size(); ++i) if(grid_[i]<min) min=grid_[i];
667 for(index_t i=0; i<grid_.size(); ++i) grid_[i] -= min;
668 }
669
applyFunctionAllValuesAndDerivatives(double (* func)(double val),double (* funcder)(double valder))670 void Grid::applyFunctionAllValuesAndDerivatives( double (*func)(double val), double (*funcder)(double valder) ) {
671 if(usederiv_) {
672 for(index_t i=0; i<grid_.size(); ++i) {
673 grid_[i]=func(grid_[i]);
674 for(unsigned j=0; j<dimension_; ++j) der_[i*dimension_+j]=funcder(der_[i*dimension_+j]);
675 }
676 } else {
677 for(index_t i=0; i<grid_.size(); ++i) grid_[i]=func(grid_[i]);
678 }
679 }
680
getDifferenceFromContour(const std::vector<double> & x,std::vector<double> & der) const681 double Grid::getDifferenceFromContour( const std::vector<double>& x, std::vector<double>& der ) const {
682 return getValueAndDerivatives( x, der ) - contour_location;
683 }
684
findSetOfPointsOnContour(const double & target,const std::vector<bool> & nosearch,unsigned & npoints,std::vector<std::vector<double>> & points)685 void Grid::findSetOfPointsOnContour(const double& target, const std::vector<bool>& nosearch,
686 unsigned& npoints, std::vector<std::vector<double> >& points ) {
687 // Set contour location for function
688 contour_location=target;
689 // Resize points to maximum possible value
690 points.resize( dimension_*maxsize_ );
691
692 // Two points for search
693 std::vector<unsigned> ind(dimension_);
694 std::vector<double> direction( dimension_, 0 );
695
696 // Run over whole grid
697 npoints=0; RootFindingBase<Grid> mymin( this );
698 for(unsigned i=0; i<maxsize_; ++i) {
699 for(unsigned j=0; j<dimension_; ++j) ind[j]=getIndices(i)[j];
700
701 // Get the value of a point on the grid
702 double val1=getValue(i) - target;
703
704 // Now search for contour in each direction
705 bool edge=false;
706 for(unsigned j=0; j<dimension_; ++j) {
707 if( nosearch[j] ) continue ;
708 // Make sure we don't search at the edge of the grid
709 if( !pbc_[j] && (ind[j]+1)==nbin_[j] ) continue;
710 else if( (ind[j]+1)==nbin_[j] ) { edge=true; ind[j]=0; }
711 else ind[j]+=1;
712 double val2=getValue(ind) - target;
713 if( val1*val2<0 ) {
714 // Use initial point location as first guess for search
715 points[npoints].resize(dimension_); for(unsigned k=0; k<dimension_; ++k) points[npoints][k]=getPoint(i)[k];
716 // Setup direction vector
717 direction[j]=0.999999999*dx_[j];
718 // And do proper search for contour point
719 mymin.linesearch( direction, points[npoints], &Grid::getDifferenceFromContour );
720 direction[j]=0.0; npoints++;
721 }
722 if( pbc_[j] && edge ) { edge=false; ind[j]=nbin_[j]-1; }
723 else ind[j]-=1;
724 }
725 }
726 }
727
728 /// OVERRIDES ARE BELOW
729
getSize() const730 Grid::index_t Grid::getSize() const {
731 return maxsize_;
732 }
733
getValue(index_t index) const734 double Grid::getValue(index_t index) const {
735 plumed_dbg_assert(index<maxsize_);
736 return grid_[index];
737 }
738
getValueAndDerivatives(index_t index,vector<double> & der) const739 double Grid::getValueAndDerivatives
740 (index_t index, vector<double>& der) const {
741 plumed_dbg_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
742 der.resize(dimension_);
743 for(unsigned i=0; i<dimension_; i++) der[i]=der_[dimension_*index+i];
744 return grid_[index];
745 }
746
setValue(index_t index,double value)747 void Grid::setValue(index_t index, double value) {
748 plumed_dbg_assert(index<maxsize_ && !usederiv_);
749 grid_[index]=value;
750 }
751
setValueAndDerivatives(index_t index,double value,vector<double> & der)752 void Grid::setValueAndDerivatives
753 (index_t index, double value, vector<double>& der) {
754 plumed_dbg_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
755 grid_[index]=value;
756 for(unsigned i=0; i<dimension_; i++) der_[dimension_*index+i]=der[i];
757 }
758
addValue(index_t index,double value)759 void Grid::addValue(index_t index, double value) {
760 plumed_dbg_assert(index<maxsize_ && !usederiv_);
761 grid_[index]+=value;
762 }
763
addValueAndDerivatives(index_t index,double value,vector<double> & der)764 void Grid::addValueAndDerivatives
765 (index_t index, double value, vector<double>& der) {
766 plumed_dbg_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
767 grid_[index]+=value;
768 for(unsigned int i=0; i<dimension_; ++i) der_[index*dimension_+i]+=der[i];
769 }
770
getSize() const771 Grid::index_t SparseGrid::getSize() const {
772 return map_.size();
773 }
774
getMaxSize() const775 Grid::index_t SparseGrid::getMaxSize() const {
776 return maxsize_;
777 }
778
getValue(index_t index) const779 double SparseGrid::getValue(index_t index)const {
780 plumed_assert(index<maxsize_);
781 double value=0.0;
782 const auto it=map_.find(index);
783 if(it!=map_.end()) value=it->second;
784 return value;
785 }
786
getValueAndDerivatives(index_t index,vector<double> & der) const787 double SparseGrid::getValueAndDerivatives
788 (index_t index, vector<double>& der)const {
789 plumed_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
790 double value=0.0;
791 for(unsigned int i=0; i<dimension_; ++i) der[i]=0.0;
792 const auto it=map_.find(index);
793 if(it!=map_.end()) value=it->second;
794 const auto itder=der_.find(index);
795 if(itder!=der_.end()) der=itder->second;
796 return value;
797 }
798
setValue(index_t index,double value)799 void SparseGrid::setValue(index_t index, double value) {
800 plumed_assert(index<maxsize_ && !usederiv_);
801 map_[index]=value;
802 }
803
setValueAndDerivatives(index_t index,double value,vector<double> & der)804 void SparseGrid::setValueAndDerivatives
805 (index_t index, double value, vector<double>& der) {
806 plumed_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
807 map_[index]=value;
808 der_[index]=der;
809 }
810
addValue(index_t index,double value)811 void SparseGrid::addValue(index_t index, double value) {
812 plumed_assert(index<maxsize_ && !usederiv_);
813 map_[index]+=value;
814 }
815
addValueAndDerivatives(index_t index,double value,vector<double> & der)816 void SparseGrid::addValueAndDerivatives
817 (index_t index, double value, vector<double>& der) {
818 plumed_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
819 map_[index]+=value;
820 der_[index].resize(dimension_);
821 for(unsigned int i=0; i<dimension_; ++i) der_[index][i]+=der[i];
822 }
823
writeToFile(OFile & ofile)824 void SparseGrid::writeToFile(OFile& ofile) {
825 vector<double> xx(dimension_);
826 vector<double> der(dimension_);
827 double f;
828 writeHeader(ofile);
829 ofile.fmtField(" "+fmt_);
830 for(const auto & it : map_) {
831 index_t i=it.first;
832 xx=getPoint(i);
833 if(usederiv_) {f=getValueAndDerivatives(i,der);}
834 else {f=getValue(i);}
835 if(i>0 && dimension_>1 && getIndices(i)[dimension_-2]==0) ofile.printf("\n");
836 for(unsigned j=0; j<dimension_; ++j) {
837 ofile.printField("min_" + argnames[j], str_min_[j] );
838 ofile.printField("max_" + argnames[j], str_max_[j] );
839 ofile.printField("nbins_" + argnames[j], static_cast<int>(nbin_[j]) );
840 if( pbc_[j] ) ofile.printField("periodic_" + argnames[j], "true" );
841 else ofile.printField("periodic_" + argnames[j], "false" );
842 }
843 for(unsigned j=0; j<dimension_; ++j) ofile.printField(argnames[j],xx[j]);
844 ofile.printField(funcname, f);
845 if(usederiv_) { for(unsigned j=0; j<dimension_; ++j) ofile.printField("der_" + argnames[j],der[j]); }
846 ofile.printField();
847 }
848 }
849
getMinValue() const850 double SparseGrid::getMinValue() const {
851 double minval;
852 minval=0.0;
853 for(auto const & i : map_) {
854 if(i.second<minval) minval=i.second;
855 }
856 return minval;
857 }
858
getMaxValue() const859 double SparseGrid::getMaxValue() const {
860 double maxval;
861 maxval=0.0;
862 for(auto const & i : map_) {
863 if(i.second>maxval) maxval=i.second;
864 }
865 return maxval;
866 }
867
projectOnLowDimension(double & val,std::vector<int> & vHigh,WeightBase * ptr2obj)868 void Grid::projectOnLowDimension(double &val, std::vector<int> &vHigh, WeightBase * ptr2obj ) {
869 unsigned i=0;
870 for(i=0; i<vHigh.size(); i++) {
871 if(vHigh[i]<0) { // this bin needs to be integrated out
872 // parallelize here???
873 for(unsigned j=0; j<(getNbin())[i]; j++) {
874 vHigh[i]=int(j);
875 projectOnLowDimension(val,vHigh,ptr2obj); // recursive function: this is the core of the mechanism
876 vHigh[i]=-1;
877 }
878 return; //
879 }
880 }
881 // when there are no more bin to dig in then retrieve the value
882 if(i==vHigh.size()) {
883 //std::cerr<<"POINT: ";
884 //for(unsigned j=0;j<vHigh.size();j++){
885 // std::cerr<<vHigh[j]<<" ";
886 //}
887 std::vector<unsigned> vv(vHigh.size());
888 for(unsigned j=0; j<vHigh.size(); j++)vv[j]=unsigned(vHigh[j]);
889 //
890 // this is the real assignment !!!!! (hack this to have bias or other stuff)
891 //
892
893 // this case: produce fes
894 //val+=exp(beta*getValue(vv)) ;
895 double myv=getValue(vv);
896 val=ptr2obj->projectInnerLoop(val,myv) ;
897 // to be added: bias (same as before without negative sign)
898 //std::cerr<<" VAL: "<<val <<endl;
899 }
900 }
901
project(const std::vector<std::string> & proj,WeightBase * ptr2obj)902 Grid Grid::project(const std::vector<std::string> & proj, WeightBase *ptr2obj ) {
903 // find extrema only for the projection
904 vector<string> smallMin,smallMax;
905 vector<unsigned> smallBin;
906 vector<unsigned> dimMapping;
907 vector<bool> smallIsPeriodic;
908 vector<string> smallName;
909
910 // check if the two key methods are there
911 WeightBase* pp = dynamic_cast<WeightBase*>(ptr2obj);
912 if (!pp)plumed_merror("This WeightBase is not complete: you need a projectInnerLoop and projectOuterLoop ");
913
914 for(unsigned j=0; j<proj.size(); j++) {
915 for(unsigned i=0; i<getArgNames().size(); i++) {
916 if(proj[j]==getArgNames()[i]) {
917 unsigned offset;
918 // note that at sizetime the non periodic dimension get a bin more
919 if(getIsPeriodic()[i]) {offset=0;} else {offset=1;}
920 smallMax.push_back(getMax()[i]);
921 smallMin.push_back(getMin()[i]);
922 smallBin.push_back(getNbin()[i]-offset);
923 smallIsPeriodic.push_back(getIsPeriodic()[i]);
924 dimMapping.push_back(i);
925 smallName.push_back(getArgNames()[i]);
926 break;
927 }
928 }
929 }
930 Grid smallgrid("projection",smallName,smallMin,smallMax,smallBin,false,false,smallIsPeriodic,smallMin,smallMax);
931 // check that the two grids are commensurate
932 for(unsigned i=0; i<dimMapping.size(); i++) {
933 plumed_massert( (smallgrid.getMax())[i] == (getMax())[dimMapping[i]], "the two input grids are not compatible in max" );
934 plumed_massert( (smallgrid.getMin())[i] == (getMin())[dimMapping[i]], "the two input grids are not compatible in min" );
935 plumed_massert( (smallgrid.getNbin())[i]== (getNbin())[dimMapping[i]], "the two input grids are not compatible in bin" );
936 }
937 vector<unsigned> toBeIntegrated;
938 for(unsigned i=0; i<getArgNames().size(); i++) {
939 bool doappend=true;
940 for(unsigned j=0; j<dimMapping.size(); j++) {
941 if(dimMapping[j]==i) {doappend=false; break;}
942 }
943 if(doappend)toBeIntegrated.push_back(i);
944 }
945 //for(unsigned i=0;i<dimMapping.size();i++ ){
946 // cerr<<"Dimension to preserve "<<dimMapping[i]<<endl;
947 //}
948 //for(unsigned i=0;i<toBeIntegrated.size();i++ ){
949 // cerr<<"Dimension to integrate "<<toBeIntegrated[i]<<endl;
950 //}
951
952 // loop over all the points in the Grid, find the corresponding fixed index, rotate over all the other ones
953 for(unsigned i=0; i<smallgrid.getSize(); i++) {
954 std::vector<unsigned> v;
955 v=smallgrid.getIndices(i);
956 std::vector<int> vHigh((getArgNames()).size(),-1);
957 for(unsigned j=0; j<dimMapping.size(); j++)vHigh[dimMapping[j]]=int(v[j]);
958 // the vector vhigh now contains at the beginning the index of the low dimension and -1 for the values that need to be integrated
959 double initval=0.;
960 projectOnLowDimension(initval,vHigh, ptr2obj);
961 smallgrid.setValue(i,initval);
962 }
963 // reset to zero just for biasing (this option can be evtl enabled in a future...)
964 //double vmin;vmin=-smallgrid.getMinValue()+1;
965 for(unsigned i=0; i<smallgrid.getSize(); i++) {
966 // //if(dynamic_cast<BiasWeight*>(ptr2obj)){
967 // // smallgrid.addValue(i,vmin);// go to 1
968 // //}
969 double vv=smallgrid.getValue(i);
970 smallgrid.setValue(i,ptr2obj->projectOuterLoop(vv));
971 // //if(dynamic_cast<BiasWeight*>(ptr2obj)){
972 // // smallgrid.addValue(i,-vmin);// bring back to the value
973 // //}
974 }
975
976 return smallgrid;
977 }
978
integrate(std::vector<unsigned> & npoints)979 double Grid::integrate( std::vector<unsigned>& npoints ) {
980 plumed_dbg_assert( npoints.size()==dimension_ ); plumed_assert( dospline_ );
981
982 unsigned ntotgrid=1; double box_vol=1.0;
983 std::vector<double> ispacing( npoints.size() );
984 for(unsigned j=0; j<dimension_; ++j) {
985 if( !pbc_[j] ) {
986 ispacing[j] = ( max_[j] - dx_[j] - min_[j] ) / static_cast<double>( npoints[j] );
987 npoints[j]+=1;
988 } else {
989 ispacing[j] = ( max_[j] - min_[j] ) / static_cast<double>( npoints[j] );
990 }
991 ntotgrid*=npoints[j]; box_vol*=ispacing[j];
992 }
993
994 std::vector<double> vals( dimension_ );
995 std::vector<unsigned> t_index( dimension_ ); double integral=0.0;
996 for(unsigned i=0; i<ntotgrid; ++i) {
997 t_index[0]=(i%npoints[0]);
998 unsigned kk=i;
999 for(unsigned j=1; j<dimension_-1; ++j) { kk=(kk-t_index[j-1])/npoints[j-1]; t_index[j]=(kk%npoints[j]); }
1000 if( dimension_>=2 ) t_index[dimension_-1]=((kk-t_index[dimension_-2])/npoints[dimension_-2]);
1001
1002 for(unsigned j=0; j<dimension_; ++j) vals[j]=min_[j] + t_index[j]*ispacing[j];
1003
1004 integral += getValue( vals );
1005 }
1006
1007 return box_vol*integral;
1008 }
1009
mpiSumValuesAndDerivatives(Communicator & comm)1010 void Grid::mpiSumValuesAndDerivatives( Communicator& comm ) {
1011 comm.Sum( grid_ ); for(unsigned i=0; i<der_.size(); ++i) comm.Sum( der_[i] );
1012 }
1013
1014
indexed_lt(pair<Grid::index_t,double> const & x,pair<Grid::index_t,double> const & y)1015 bool indexed_lt(pair<Grid::index_t, double> const &x, pair<Grid::index_t, double> const &y) {
1016 return x.second < y.second;
1017 }
1018
findMaximalPathMinimum(const std::vector<double> & source,const std::vector<double> & sink)1019 double GridBase::findMaximalPathMinimum(const std::vector<double> &source, const std::vector<double> &sink) {
1020 plumed_dbg_assert(source.size() == dimension_);
1021 plumed_dbg_assert(sink.size() == dimension_);
1022 // Start and end indices
1023 index_t source_idx = getIndex(source);
1024 index_t sink_idx = getIndex(sink);
1025 // Path cost
1026 double maximal_minimum = 0;
1027 // In one dimension, path searching is very easy--either go one way if it's not periodic,
1028 // or go both ways if it is periodic. There's no reason to pay the cost of Dijkstra.
1029 if (dimension_ == 1) {
1030 // Do a search from the grid source to grid sink that does not
1031 // cross the grid boundary.
1032 double curr_min_bias = getValue(source_idx);
1033 // Either search from a high source to a low sink.
1034 if (source_idx > sink_idx) {
1035 for (index_t i = source_idx; i >= sink_idx; i--) {
1036 if (curr_min_bias == 0.0) {
1037 break;
1038 }
1039 curr_min_bias = fmin(curr_min_bias, getValue(i));
1040 }
1041 // Or search from a low source to a high sink.
1042 } else if (source_idx < sink_idx) {
1043 for (index_t i = source_idx; i <= sink_idx; i++) {
1044 if (curr_min_bias == 0.0) {
1045 break;
1046 }
1047 curr_min_bias = fmin(curr_min_bias, getValue(i));
1048 }
1049 }
1050 maximal_minimum = curr_min_bias;
1051 // If the grid is periodic, also do the search that crosses
1052 // the grid boundary.
1053 if (pbc_[0]) {
1054 double curr_min_bias = getValue(source_idx);
1055 // Either go from a high source to the upper boundary and
1056 // then from the bottom boundary to the sink
1057 if (source_idx > sink_idx) {
1058 for (index_t i = source_idx; i < maxsize_; i++) {
1059 if (curr_min_bias == 0.0) {
1060 break;
1061 }
1062 curr_min_bias = fmin(curr_min_bias, getValue(i));
1063 }
1064 for (index_t i = 0; i <= sink_idx; i++) {
1065 if (curr_min_bias == 0.0) {
1066 break;
1067 }
1068 curr_min_bias = fmin(curr_min_bias, getValue(i));
1069 }
1070 // Or go from a low source to the bottom boundary and
1071 // then from the high boundary to the sink
1072 } else if (source_idx < sink_idx) {
1073 for (index_t i = source_idx; i > 0; i--) {
1074 if (curr_min_bias == 0.0) {
1075 break;
1076 }
1077 curr_min_bias = fmin(curr_min_bias, getValue(i));
1078 }
1079 curr_min_bias = fmin(curr_min_bias, getValue(0));
1080 for (index_t i = maxsize_ - 1; i <= sink_idx; i--) {
1081 if (curr_min_bias == 0.0) {
1082 break;
1083 }
1084 curr_min_bias = fmin(curr_min_bias, getValue(i));
1085 }
1086 }
1087 // If the boundary crossing paths was more biased, it's
1088 // minimal bias replaces the non-boundary-crossing path's
1089 // minimum.
1090 maximal_minimum = fmax(maximal_minimum, curr_min_bias);
1091 }
1092 // The one dimensional path search is complete.
1093 return maximal_minimum;
1094 // In two or more dimensions, path searching isn't trivial and we really
1095 // do need to use a path search algorithm. Dijkstra is the simplest decent
1096 // one. Using it we've never found the path search to be performance
1097 // limiting in any solvated biomolecule test system, but faster options are
1098 // easy to imagine if they become necessary. NB-In this case, we're actually
1099 // using a greedy variant of Dijkstra's algorithm where the first possible
1100 // path to a point always controls the path cost to that point. The structure
1101 // of the cost function in this case guarantees that the calculated costs will
1102 // be correct using this variant even though fine details of the paths may not
1103 // match a normal Dijkstra search.
1104 } else if (dimension_ > 1) {
1105 // Prepare calculation temporaries for Dijkstra's algorithm.
1106 // Minimal path costs from source to a given grid point
1107 vector<double> mins_from_source = vector<double>(maxsize_, -1.0);
1108 // Heap for tracking available steps, steps are recorded as std::pairs of
1109 // an index and a value.
1110 vector< pair<index_t, double> > next_steps;
1111 pair<index_t, double> curr_indexed_val;
1112 make_heap(next_steps.begin(), next_steps.end(), indexed_lt);
1113 // The search begins at the source index.
1114 next_steps.push_back(pair<index_t, double>(source_idx, getValue(source_idx)));
1115 push_heap(next_steps.begin(), next_steps.end(), indexed_lt);
1116 // At first no points have been examined and the optimal path has not been found.
1117 index_t n_examined = 0;
1118 bool path_not_found = true;
1119 // Until a path is found,
1120 while (path_not_found) {
1121 // Examine the grid point currently most accessible from
1122 // the set of all previously explored grid points by popping
1123 // it from the top of the heap.
1124 pop_heap(next_steps.begin(), next_steps.end(), indexed_lt);
1125 curr_indexed_val = next_steps.back();
1126 next_steps.pop_back();
1127 n_examined++;
1128 // Check if this point is the sink point, and if so
1129 // finish the loop.
1130 if (curr_indexed_val.first == sink_idx) {
1131 path_not_found = false;
1132 maximal_minimum = curr_indexed_val.second;
1133 break;
1134 // Check if this point has reached the worst possible
1135 // value, and if so stop looking for paths.
1136 } else if (curr_indexed_val.second == 0.0) {
1137 maximal_minimum = 0.0;
1138 break;
1139 }
1140 // If the search is not over, add this grid point's neighbors to the
1141 // possible next points to search for the sink.
1142 vector<index_t> neighs = getNearestNeighbors(curr_indexed_val.first);
1143 for (unsigned k = 0; k < neighs.size(); k++) {
1144 index_t i = neighs[k];
1145 // If the neighbor has not already been added to the list of possible next steps,
1146 if (mins_from_source[i] == -1.0) {
1147 // Set the cost to reach it via a path through the current point being examined.
1148 mins_from_source[i] = fmin(curr_indexed_val.second, getValue(i));
1149 // Add the neighboring point to the heap of potential next steps.
1150 next_steps.push_back(pair<index_t, double>(i, mins_from_source[i]));
1151 push_heap(next_steps.begin(), next_steps.end(), indexed_lt);
1152 }
1153 }
1154 // Move on to the next best looking step along any of the paths
1155 // growing from the source.
1156 }
1157 // The multidimensional path search is now complete.
1158 return maximal_minimum;
1159 }
1160 return 0.0;
1161 }
1162
1163 }
1164