Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2026 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 "MetaDGrid.h"
23 : #include "MetaDHills.h"
24 : #include "core/Action.h"
25 : #include "core/Value.h"
26 : #include "tools/File.h"
27 : #include "tools/OpenMP.h"
28 : #include "tools/Tools.h"
29 : #include <algorithm>
30 : #include <cmath>
31 : #include <limits>
32 :
33 : namespace PLMD {
34 : namespace bias {
35 :
36 339 : bool MetaDBiasGrid::setupBins(Action& action, const std::vector<Value*>& args, std::vector<std::string>& gmin, std::vector<std::string>& gmax,
37 : std::vector<unsigned>& gbin, std::vector<double>& gspacing, bool adaptive, const std::vector<double>& sigma0,
38 : const std::vector<double>& sigma0min) {
39 : bool nobounds=false;
40 339 : if(gmin.size()!=gmax.size()) {
41 0 : action.error("GRID_MAX and GRID_MIN should be either present or absent");
42 : }
43 339 : if(gbin.size()!=0 && gmin.size()==0) {
44 0 : action.error("If GRID_BIN is present also GRID_MIN and GRID_MAX should be present");
45 : }
46 489 : for(const auto s : gspacing) {
47 150 : if(!(s>0.0)) {
48 0 : action.error("GRID_SPACING should be positive");
49 : }
50 : }
51 339 : if(gmin.size()!=0) {
52 144 : if(gbin.size()==0 && gspacing.size()==0) {
53 17 : if(!adaptive) {
54 13 : action.log<<" Binsize not specified, 1/5 of sigma will be used\n";
55 13 : plumed_assert(sigma0.size()==args.size());
56 13 : gspacing.resize(args.size());
57 33 : for(unsigned i=0; i<gspacing.size(); i++) {
58 20 : gspacing[i]=0.2*sigma0[i];
59 : }
60 : } else {
61 : // with adaptive hills and grid a sigma min must be specified
62 12 : for(unsigned i=0; i<sigma0min.size(); i++)
63 8 : if(sigma0min[i]<=0) {
64 0 : action.error("When using ADAPTIVE Gaussians on a grid SIGMA_MIN must be specified");
65 : }
66 4 : action.log<<" Binsize not specified, 1/5 of sigma_min will be used\n";
67 4 : gspacing.resize(args.size());
68 12 : for(unsigned i=0; i<gspacing.size(); i++) {
69 8 : gspacing[i]=0.2*sigma0min[i];
70 : }
71 : }
72 127 : } else if(gspacing.size()!=0 && gbin.size()==0) {
73 60 : action.log<<" The number of bins will be estimated from GRID_SPACING\n";
74 67 : } else if(gspacing.size()!=0 && gbin.size()!=0) {
75 1 : action.log<<" You specified both GRID_BIN and GRID_SPACING\n";
76 1 : action.log<<" The more conservative (highest) number of bins will be used for each variable\n";
77 : }
78 144 : if(gbin.size()==0) {
79 77 : gbin.assign(args.size(),1);
80 : }
81 144 : if(gspacing.size()!=0)
82 224 : for(unsigned i=0; i<args.size(); i++) {
83 : double a,b;
84 147 : Tools::convert(gmin[i],a);
85 146 : Tools::convert(gmax[i],b);
86 146 : unsigned n=std::ceil(((b-a)/gspacing[i]));
87 146 : if(gbin[i]<n) {
88 145 : gbin[i]=n;
89 : }
90 : }
91 : }
92 338 : if(gspacing.size()!=0 && gmin.size()==0) {
93 : // no range: the grid is placed around the first hill, on the lattice with a point at 0
94 : // (a first range of one bin is used until then); periodic variables are covered over the whole period
95 : nobounds=true;
96 18 : gmin.resize(args.size());
97 18 : gmax.resize(args.size());
98 18 : gbin.resize(args.size());
99 49 : for(unsigned i=0; i<args.size(); i++) {
100 31 : if(!(gspacing[i]>0.0)) {
101 0 : action.error("GRID_SPACING should be positive");
102 : }
103 31 : if(args[i]->isPeriodic()) {
104 13 : args[i]->getDomain(gmin[i],gmax[i]);
105 : double a,b;
106 13 : Tools::convert(gmin[i],a);
107 13 : Tools::convert(gmax[i],b);
108 13 : gbin[i]=std::ceil((b-a)/gspacing[i]);
109 : } else {
110 : gmin[i]="0";
111 36 : gmax[i]=exactString(gspacing[i]);
112 18 : gbin[i]=1;
113 : }
114 : }
115 : }
116 338 : return nobounds;
117 : }
118 :
119 205 : void MetaDBiasGrid::checkSpacing(Action& action, const std::vector<std::string>& gmin, const std::vector<std::string>& gmax,
120 : const std::vector<unsigned>& gbin, bool adaptive, const std::vector<double>& sigma0,
121 : const std::vector<double>& sigma0min) {
122 497 : for(unsigned i=0; i<gbin.size(); i++) {
123 : double a,b;
124 292 : Tools::convert(gmin[i],a);
125 292 : Tools::convert(gmax[i],b);
126 292 : const double mesh=(b-a)/static_cast<double>(gbin[i]);
127 292 : if(!adaptive) {
128 263 : if(mesh>0.5*sigma0[i]) {
129 67 : action.log<<" WARNING: Using a "<<action.getName()<<" with a Grid Spacing larger than half of the Gaussians width (SIGMA) can produce artifacts\n";
130 : }
131 : } else {
132 29 : if(!(sigma0min[i]>0.0)) {
133 0 : action.error("When using ADAPTIVE Gaussians on a grid SIGMA_MIN must be specified");
134 : }
135 29 : if(mesh>0.5*sigma0min[i]) {
136 12 : action.log<<" WARNING: to use a "<<action.getName()<<" with a GRID and ADAPTIVE you need to set a Grid Spacing lower than half of the Gaussians (SIGMA_MIN) \n";
137 : }
138 : }
139 : }
140 205 : }
141 :
142 338 : void MetaDBiasGrid::checkSizeLimit(Action& action, double limit) {
143 338 : if(!(limit>=1.0)) {
144 0 : action.error("GRID_SIZE_LIMIT should be at least 1");
145 : }
146 : // grid points are indexed with unsigned integers (see GridBase)
147 338 : if(limit>static_cast<double>(std::numeric_limits<unsigned>::max())) {
148 0 : action.error("GRID_SIZE_LIMIT cannot be larger than "+exactString(std::numeric_limits<unsigned>::max()));
149 : }
150 338 : }
151 :
152 205 : MetaDBiasGrid::MetaDBiasGrid(Action& action):
153 205 : action_(action) {
154 205 : }
155 :
156 162 : void MetaDBiasGrid::create(const std::string& funcl, const std::vector<Value*>& args, const std::vector<std::string>& gmin,
157 : const std::vector<std::string>& gmax, const std::vector<unsigned>& gbin, bool sparse, bool spline, bool nobounds) {
158 162 : sparse_=sparse;
159 162 : if(!sparse) {
160 148 : grid_=Tools::make_unique<Grid>(funcl,args,gmin,gmax,gbin,spline,true);
161 : } else {
162 14 : grid_=Tools::make_unique<SparseGrid>(funcl,args,gmin,gmax,gbin,spline,true);
163 : }
164 162 : pristine_=nobounds;
165 162 : std::vector<std::string> actualmin=grid_->getMin();
166 162 : std::vector<std::string> actualmax=grid_->getMax();
167 392 : for(unsigned i=0; i<args.size(); i++) {
168 : std::string is;
169 230 : Tools::convert(i,is);
170 230 : if(gmin[i]!=actualmin[i]) {
171 0 : action_.error("GRID_MIN["+is+"] must be adjusted to "+actualmin[i]+" to fit periodicity");
172 : }
173 230 : if(gmax[i]!=actualmax[i]) {
174 0 : action_.error("GRID_MAX["+is+"] must be adjusted to "+actualmax[i]+" to fit periodicity");
175 : }
176 : }
177 162 : }
178 :
179 43 : void MetaDBiasGrid::read(IFile& file, const std::string& filename, const std::string& funcl, const std::vector<Value*>& args,
180 : const std::vector<std::string>& gmin, const std::vector<std::string>& gmax, const std::vector<unsigned>& gbin,
181 : bool sparse, bool spline, bool nobounds) {
182 43 : sparse_=sparse;
183 85 : grid_=GridBase::create(funcl, args, file, sparse, spline, true);
184 42 : if(grid_->getDimension()!=args.size()) {
185 0 : action_.error("mismatch between dimensionality of input grid and number of arguments");
186 : }
187 103 : for(unsigned i=0; i<args.size(); ++i) {
188 122 : if( args[i]->isPeriodic()!=grid_->getIsPeriodic()[i] ) {
189 0 : action_.error("periodicity mismatch between arguments and input bias");
190 : }
191 : // the grid in the file can have a different range (e.g. it was extended), but it must be on the same
192 : // lattice as the grid in input: same spacing, points at the same positions
193 : std::string is;
194 61 : Tools::convert(i,is);
195 : double min,max,fmin;
196 61 : Tools::convert(gmin[i],min);
197 61 : Tools::convert(gmax[i],max);
198 61 : Tools::convert(grid_->getMin()[i],fmin);
199 61 : const double dx=(max-min)/gbin[i];
200 61 : const double fdx=grid_->getDx(i);
201 61 : if(args[i]->isPeriodic()) {
202 27 : if(grid_->getNbin()[i]!=gbin[i]) {
203 0 : action_.error("the number of bins of the grid in "+filename+" for the periodic variable "+is+" does not match the input");
204 : }
205 : } else {
206 34 : const double shift=(fmin-min)/dx;
207 : // (the bounds of an extended grid are written with 12 significant digits)
208 34 : if(std::fabs(fdx-dx)>1e-6*dx || std::fabs(shift-std::round(shift))>1e-4) {
209 0 : action_.error("the grid in "+filename+" for variable "+is+" is not on the lattice of the grid in input: "
210 0 : "its spacing should be "+exactString(dx)+" and its points should be at "+gmin[i]+" plus multiples of the spacing");
211 : }
212 : }
213 : }
214 : // a grid without initial range read from a file written before the first hill: it is placed around the first hill
215 42 : pristine_=false;
216 42 : if(nobounds) {
217 2 : bool empty=true;
218 2 : grid_->forEachStoredPoint([&](GridBase::index_t,double value) {
219 55 : empty=empty && value==0.0;
220 : });
221 2 : pristine_=empty;
222 : }
223 42 : }
224 :
225 204 : void MetaDBiasGrid::setSizeLimit(double limit) {
226 204 : sizelimit_=limit;
227 : // an initial grid larger than GRID_SIZE_LIMIT (e.g. given with GRID_MIN, GRID_MAX, GRID_BIN) can still be extended
228 : double size=1.0;
229 495 : for(const auto n : grid_->getNbin()) {
230 291 : size*=n;
231 : }
232 204 : if(size>sizelimit_) {
233 1 : sizelimit_=std::min(2.0*size,static_cast<double>(std::numeric_limits<unsigned>::max()));
234 2 : action_.log.printf(" The initial grid has more points than GRID_SIZE_LIMIT: the limit is set to %s points\n",exactString(sizelimit_).c_str());
235 : }
236 204 : }
237 :
238 15 : void MetaDBiasGrid::setInterval(double lower, double upper, double time) {
239 15 : interval_=true;
240 15 : lower_=lower;
241 15 : upper_=upper;
242 : // the bias beyond the limits is constant: the grid contains the interval and at least one point
243 : // beyond each limit, which holds this constant value (see extend)
244 15 : const std::vector<double> center(1,0.5*(lower_+upper_));
245 15 : const std::vector<double> halfwidth(1,0.5*(upper_-lower_)+grid_->getDx(0));
246 15 : extend(center,halfwidth,time);
247 15 : }
248 :
249 12458 : void MetaDBiasGrid::extend(const std::vector<double>& x, const std::vector<double>& halfwidth, double time) {
250 12458 : const unsigned ncv=grid_->getDimension();
251 12458 : const std::vector<bool> pbc=grid_->getIsPeriodic();
252 12458 : const std::vector<unsigned> nbin=grid_->getNbin();
253 : // points to add at the lower and upper end along each CV (for a pristine grid: the change of the range, also negative)
254 12458 : std::vector<long long> lower(ncv,0), upper(ncv,0);
255 : // as above, with some more points so that the grid is not extended at every hill
256 12458 : std::vector<long long> lowergrow(ncv,0), uppergrow(ncv,0);
257 : bool extend=false;
258 29421 : for(unsigned i=0; i<ncv; i++) {
259 16963 : if(pbc[i]) {
260 8493 : continue;
261 : }
262 : // points within halfwidth on each side (as in MetaDKernel::getSupport, which adds the hill to the points from
263 : // the one below x minus this number to the one below x plus this number), plus one more on each side,
264 : // so that the bias and its derivatives are zero at the boundary of the grid
265 8470 : double n=std::ceil(halfwidth[i]/grid_->getDx(i));
266 8470 : if(!std::isfinite(n) || n<0.0) {
267 : // support not finite (degenerate hill): the hill is added to the whole grid, x only must be in it
268 : n=0.0;
269 : }
270 8470 : const double k=std::floor(grid_->getLatticeCoordinate(i,x[i]));
271 8470 : const double lo=k-n-1.0;
272 8470 : const double hi=k+n+1.0;
273 8470 : if(!std::isfinite(lo) || !std::isfinite(hi) || std::fabs(lo)>1e15 || std::fabs(hi)>1e15) {
274 : std::string is;
275 0 : Tools::convert(i,is);
276 0 : action_.error("cannot extend the grid along variable "+is+" to contain the point "+exactString(x[i]));
277 : }
278 8470 : const double last=static_cast<double>(nbin[i])-1.0;
279 8470 : if(pristine_) {
280 17 : lower[i]=lowergrow[i]=static_cast<long long>(-lo);
281 17 : upper[i]=uppergrow[i]=static_cast<long long>(hi-last);
282 : extend=true;
283 : } else {
284 : // geometric growth: at least a quarter of the current number of points
285 8453 : const long long grow=static_cast<long long>(std::ceil(0.25*nbin[i]));
286 8453 : if(lo<0.0) {
287 72 : lower[i]=static_cast<long long>(-lo);
288 72 : lowergrow[i]=std::max(lower[i],grow);
289 : extend=true;
290 : }
291 8453 : if(hi>last) {
292 57 : upper[i]=static_cast<long long>(hi-last);
293 57 : uppergrow[i]=std::max(upper[i],grow);
294 : extend=true;
295 : }
296 : }
297 : }
298 12458 : const bool pristine=pristine_;
299 12458 : pristine_=false;
300 12458 : if(!extend) {
301 : return;
302 : }
303 132 : auto size=[&](const std::vector<long long>& lo,const std::vector<long long>& up) {
304 : double n=1.0;
305 318 : for(unsigned i=0; i<ncv; i++) {
306 186 : n*=static_cast<double>(nbin[i])+lo[i]+up[i];
307 : }
308 132 : return n;
309 132 : };
310 132 : if(size(lowergrow,uppergrow)>sizelimit_) {
311 0 : lowergrow=lower;
312 0 : uppergrow=upper;
313 0 : if(size(lower,upper)>sizelimit_) {
314 : std::string point;
315 0 : for(unsigned i=0; i<ncv; i++) {
316 0 : point+=(i>0?",":"")+exactString(x[i]);
317 : }
318 0 : action_.error("a Gaussian at "+point+" requires a grid with "+exactString(size(lower,upper))+" points, more than GRID_SIZE_LIMIT="
319 0 : +exactString(sizelimit_)+": increase GRID_SIZE_LIMIT, or check the variables (e.g. a distance that grows without limit) "
320 : "and the width of the Gaussians (e.g. SIGMA_MAX with ADAPTIVE)");
321 : }
322 : }
323 132 : grid_->resize(lowergrow,uppergrow);
324 132 : if(interval_ && !pristine) {
325 : // With INTERVAL, all the points beyond a limit have the same value (the hills are computed at the limit there) and
326 : // zero derivative: the new points beyond a limit take the value of the old point at that end of the grid.
327 24 : const unsigned n=grid_->getNbin()[0];
328 24 : std::vector<double> der(1);
329 48 : auto fill=[&](unsigned from,unsigned first,unsigned last,bool beyond) {
330 48 : if(!beyond) {
331 : return;
332 : }
333 24 : const double value=grid_->getValueAndDerivatives(from,der);
334 24 : der[0]=0.0;
335 1427 : for(unsigned k=first; k<last; k++) {
336 1403 : grid_->setValueAndDerivatives(k,value,der);
337 : }
338 24 : };
339 24 : const unsigned first=static_cast<unsigned>(lowergrow[0]);
340 24 : const unsigned last=n-1-static_cast<unsigned>(uppergrow[0]);
341 : // (a point at the limit, e.g. in a grid file written before grids were extended, has the same value)
342 42 : fill(first,0,first,lowergrow[0]>0 && grid_->getPoint(first)[0]<=lower_);
343 32 : fill(last,last+1,n,uppergrow[0]>0 && grid_->getPoint(last)[0]>=upper_);
344 : }
345 132 : const std::vector<std::string> gmin=grid_->getMin();
346 132 : const std::vector<std::string> gmax=grid_->getMax();
347 132 : const std::vector<unsigned> newnbin=grid_->getNbin();
348 132 : Log& log=action_.log;
349 132 : if(midline_) {
350 6 : log.printf("\n");
351 6 : midline_=false;
352 : }
353 247 : log.printf(" %s grid: min",pristine ? "Placing the" : "Extending the");
354 318 : for(unsigned i=0; i<ncv; i++) {
355 186 : log.printf(" %s",gmin[i].c_str());
356 : }
357 132 : log.printf(" max");
358 318 : for(unsigned i=0; i<ncv; i++) {
359 186 : log.printf(" %s",gmax[i].c_str());
360 : }
361 132 : log.printf(" points");
362 318 : for(unsigned i=0; i<ncv; i++) {
363 186 : log.printf(" %u",newnbin[i]);
364 : }
365 132 : log.printf(" at time %f\n",time);
366 12458 : }
367 :
368 12413 : void MetaDBiasGrid::addHill(const MetaDGaussian& hill, const MetaDKernel& kernel, double time) {
369 : // the grid must contain all the points where the hill is not zero
370 12413 : extend(hill.center,kernel.getCutoff(hill),time);
371 12413 : if(!interval_) {
372 11692 : addHillToPoints(hill,kernel);
373 : } else {
374 721 : size_t ncv=grid_->getDimension();
375 721 : std::vector<unsigned> nneighb=kernel.getSupport(hill,*grid_);
376 721 : std::vector<Grid::index_t> neighbors=grid_->getNeighbors(hill.center,nneighb);
377 721 : std::vector<double> der(ncv);
378 721 : std::vector<double> xx(ncv);
379 : // every MPI process adds the whole hill to its copy of the grid (see addHillToPoints)
380 721 : std::vector<double> dp(ncv);
381 85934 : for(size_t i=0; i<neighbors.size(); ++i) {
382 85213 : Grid::index_t ineigh=neighbors[i];
383 170426 : for(unsigned j=0; j<ncv; ++j) {
384 85213 : der[j]=0.0;
385 : }
386 85213 : grid_->getPoint(ineigh,xx);
387 85213 : double bias=kernel.evaluateWithDerivatives(xx,hill,der.data(),dp.data());
388 85213 : grid_->addValueAndDerivatives(ineigh,bias,der);
389 : }
390 721 : }
391 12413 : }
392 :
393 11692 : void MetaDBiasGrid::addHillToPoints(const MetaDGaussian& hill, const MetaDKernel& kernel) {
394 : // The hill is evaluated on the grid points within its support. The difference of a grid point from the
395 : // center of the hill along CV d only depends on the index of the point along d, so the differences are
396 : // computed once per direction; the value at each point is then computed from them
397 : // (MetaDKernel::evaluateFromDifferences, same arithmetic as MetaDKernel::evaluateWithDerivatives).
398 11692 : const unsigned ncv=grid_->getDimension();
399 11692 : const std::vector<unsigned> nneigh=kernel.getSupport(hill,*grid_);
400 11692 : const std::vector<unsigned> nbin=grid_->getNbin();
401 11692 : const std::vector<bool> pbc=grid_->getIsPeriodic();
402 11692 : std::vector<unsigned> center(ncv);
403 11692 : grid_->getIndices(hill.center,center);
404 : // grid indices and differences along each direction (as in GridBase::getNeighbors)
405 11692 : std::vector<std::vector<unsigned>> indices(ncv);
406 11692 : std::vector<std::vector<double>> differences(ncv);
407 11692 : std::vector<unsigned> tmpindices(center);
408 11692 : std::vector<double> point(ncv);
409 27873 : for(unsigned d=0; d<ncv; d++) {
410 16181 : unsigned nb=2*nneigh[d]+1;
411 16181 : if(pbc[d] && nb>nbin[d]) {
412 : nb=nbin[d];
413 : }
414 644039 : for(unsigned k=0; k<nb; k++) {
415 627858 : int i0=static_cast<int>(k)-static_cast<int>(nneigh[d])+static_cast<int>(center[d]);
416 627858 : const int n=nbin[d];
417 627858 : if(pbc[d]) {
418 374794 : i0=((i0%n)+n)%n;
419 253064 : } else if(i0<0 || i0>=n) {
420 0 : continue;
421 : }
422 627858 : tmpindices[d]=i0;
423 627858 : grid_->getPoint(tmpindices.data(),ncv,point.data(),ncv);
424 627858 : double dd=kernel.difference(d,hill.center[d],point[d]);
425 627858 : if(!hill.multivariate) {
426 545084 : dd*=hill.invsigma[d];
427 : }
428 627858 : indices[d].push_back(i0);
429 627858 : differences[d].push_back(dd);
430 : }
431 16181 : tmpindices[d]=center[d];
432 : }
433 : // linear index of a point (column-major order, as in GridBase::getIndex)
434 11692 : std::vector<Grid::index_t> stride(ncv);
435 : Grid::index_t npoints=1;
436 27873 : for(unsigned d=0; d<ncv; d++) {
437 16181 : stride[d]=(d==0 ? 1 : stride[d-1]*nbin[d-1]);
438 16181 : npoints*=indices[d].size();
439 : }
440 11692 : if(npoints==0) {
441 : return;
442 : }
443 : // Loop over all the points of the box, the first direction being the fastest.
444 : // With MPI every process adds the whole hill to its copy of the grid: the hills are the same on all
445 : // processes, so the grids stay identical without communication (splitting the points among the processes
446 : // requires summing values and derivatives of the whole box at every hill, which costs more than computing them).
447 : // With OpenMP the box is split in blocks of consecutive points; each point is computed and added by one
448 : // thread only, with the same arithmetic as the serial loop, so the result does not depend on the number of
449 : // threads. Sparse grids are filled serially, as new points are inserted.
450 11692 : unsigned nt=OpenMP::getNumThreads();
451 11692 : if(sparse_ || npoints<static_cast<Grid::index_t>(1024*nt)) {
452 : nt=1;
453 : }
454 11692 : ThreadWorkArrays<unsigned> kwork(nt,ncv);
455 11692 : ThreadWorkArrays<double> dpwork(nt,ncv);
456 11692 : ThreadWorkArrays<double> derwork(nt,ncv);
457 : GridBase* grid=grid_.get();
458 11692 : const bool sparse=sparse_;
459 11692 : #pragma omp parallel for num_threads(nt) schedule(static,1)
460 : for(unsigned block=0; block<nt; block++) {
461 : const Grid::index_t first=(npoints*block)/nt;
462 : const Grid::index_t last=(npoints*(block+1))/nt;
463 : // indices of the first point of the block along each direction
464 : unsigned* k=kwork.get(block);
465 : Grid::index_t m=first;
466 : for(unsigned d=0; d<ncv; d++) {
467 : k[d]=m%indices[d].size();
468 : m/=indices[d].size();
469 : }
470 : double* dp=dpwork.get(block);
471 : double* der=derwork.get(block);
472 : for(Grid::index_t n=first; n<last; n++) {
473 : for(unsigned d=0; d<ncv; d++) {
474 : dp[d]=differences[d][k[d]];
475 : der[d]=0.0;
476 : }
477 : bool inside;
478 : const double bias=kernel.evaluateFromDifferences(dp,hill,der,inside);
479 : if(inside || sparse) {
480 : // adding zero outside the cutoff only matters for sparse grids, where it creates the point
481 : Grid::index_t index=0;
482 : for(unsigned d=0; d<ncv; d++) {
483 : index+=indices[d][k[d]]*stride[d];
484 : }
485 : grid->addValueAndDerivatives(index,bias,der,ncv);
486 : }
487 : // next point
488 : for(unsigned d=0; d<ncv; d++) {
489 : if(++k[d]<indices[d].size()) {
490 : break;
491 : }
492 : k[d]=0;
493 : }
494 : }
495 : }
496 23384 : }
497 :
498 9261 : double MetaDBiasGrid::getBias(const std::vector<double>& cv) const {
499 : // with INTERVAL, beyond the limits the bias is the one at the limit (the grid contains the interval and at least
500 : // one point beyond each limit)
501 9261 : if(interval_ && (cv[0]<lower_ || cv[0]>upper_)) {
502 444 : return grid_->getValue(std::vector<double>(1,std::min(std::max(cv[0],lower_),upper_)));
503 : }
504 : // the grid contains all the points where the hills are not zero
505 9113 : if(grid_->isInside(cv)) {
506 9084 : return grid_->getValue(cv);
507 : }
508 : return 0.0;
509 : }
510 :
511 38261 : double MetaDBiasGrid::getBiasAndDerivatives(const std::vector<double>& cv, std::vector<double>& der) const {
512 38261 : const unsigned ncv=grid_->getDimension();
513 : // with INTERVAL, beyond the limits the bias is the one at the limit and the force is zero
514 38261 : if(interval_ && (cv[0]<lower_ || cv[0]>upper_)) {
515 1395 : der[0]=0.0;
516 4185 : return grid_->getValue(std::vector<double>(1,std::min(std::max(cv[0],lower_),upper_)));
517 : }
518 : // the grid contains all the points where the hills are not zero
519 36866 : if(grid_->isInside(cv)) {
520 36567 : std::vector<double> vder(ncv);
521 36567 : const double bias=grid_->getValueAndDerivatives(cv,vder);
522 87519 : for(unsigned i=0; i<ncv; i++) {
523 50952 : der[i]=vder[i];
524 : }
525 : return bias;
526 : }
527 665 : for(unsigned i=0; i<ncv; i++) {
528 366 : der[i]=0.0;
529 : }
530 : return 0.0;
531 : }
532 :
533 134 : void MetaDBiasGrid::write(OFile& file, bool append, const std::string& comment) const {
534 134 : if(append) {
535 : // this results in a repetition of the header
536 46 : file.clearFields();
537 : } else {
538 : // this overwrites the grid written before
539 88 : file.rewind();
540 : }
541 134 : grid_->writeToFile(file);
542 134 : if(comment.length()>0) {
543 : // (a comment before the header would hide it to IFile)
544 0 : file.printf("%s",comment.c_str());
545 : }
546 134 : if(!append) {
547 88 : file.flush();
548 : }
549 134 : }
550 :
551 : }
552 : }
|
