Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "BiasRepresentation.h"
23 : #include "core/Value.h"
24 : #include "Communicator.h"
25 : #include <iostream>
26 : #include "KernelFunctions.h"
27 : #include "File.h"
28 : #include "Grid.h"
29 : #include "Tools.h"
30 : #include <algorithm>
31 : #include <cmath>
32 : #include <limits>
33 :
34 : namespace PLMD {
35 :
36 : namespace {
37 : /// factor converting the heights written in a hills file into the bias:
38 : /// with well-tempered metadynamics (biasf>1) the heights are scaled by biasf/(biasf-1),
39 : /// otherwise they are the bias itself
40 606 : double biasFactorToBias(double biasf) {
41 606 : return biasf>1.0 ? (biasf-1.)/biasf : 1.0;
42 : }
43 : }
44 :
45 : /// the constructor here
46 7 : BiasRepresentation::BiasRepresentation(const std::vector<Value*> & tmpvalues, Communicator &cc ):hasgrid(false),rescaledToBias(false),mycomm(cc) {
47 7 : lowI_=0.0;
48 7 : uppI_=0.0;
49 7 : doInt_=false;
50 7 : ndim=tmpvalues.size();
51 21 : for(int i=0; i<ndim; i++) {
52 14 : values.push_back(tmpvalues[i]);
53 14 : names.push_back(values[i]->getName());
54 : }
55 7 : resetSupport();
56 7 : }
57 :
58 : /// overload the constructor: add the sigma at constructor time
59 1 : BiasRepresentation::BiasRepresentation(const std::vector<Value*> & tmpvalues, Communicator &cc, const std::vector<double> & sigma ):
60 1 : hasgrid(false), rescaledToBias(false), histosigma(sigma),mycomm(cc) {
61 1 : lowI_=0.0;
62 1 : uppI_=0.0;
63 1 : doInt_=false;
64 1 : ndim=tmpvalues.size();
65 3 : for(int i=0; i<ndim; i++) {
66 2 : values.push_back(tmpvalues[i]);
67 2 : names.push_back(values[i]->getName());
68 : }
69 1 : resetSupport();
70 1 : }
71 :
72 : /// overload the constructor: add the grid at constructor time
73 14 : BiasRepresentation::BiasRepresentation(const std::vector<Value*> & tmpvalues, Communicator &cc, const std::vector<std::string> & gmin, const std::vector<std::string> & gmax,
74 14 : const std::vector<unsigned> & nbin, bool doInt, double lowI, double uppI):
75 14 : hasgrid(false), rescaledToBias(false), mycomm(cc) {
76 14 : ndim=tmpvalues.size();
77 38 : for(int i=0; i<ndim; i++) {
78 24 : values.push_back(tmpvalues[i]);
79 24 : names.push_back(values[i]->getName());
80 : }
81 14 : resetSupport();
82 14 : doInt_=doInt;
83 14 : lowI_=lowI;
84 14 : uppI_=uppI;
85 : // initialize the grid
86 14 : addGrid(gmin,gmax,nbin);
87 14 : }
88 :
89 : /// overload the constructor with some external sigmas: needed for histogram
90 1 : BiasRepresentation::BiasRepresentation(const std::vector<Value*> & tmpvalues, Communicator &cc, const std::vector<std::string> & gmin, const std::vector<std::string> & gmax,
91 1 : const std::vector<unsigned> & nbin, const std::vector<double> & sigma):
92 1 : hasgrid(false), rescaledToBias(false),histosigma(sigma),mycomm(cc) {
93 1 : lowI_=0.0;
94 1 : uppI_=0.0;
95 1 : doInt_=false;
96 1 : ndim=tmpvalues.size();
97 3 : for(int i=0; i<ndim; i++) {
98 2 : values.push_back(tmpvalues[i]);
99 2 : names.push_back(values[i]->getName());
100 : }
101 1 : resetSupport();
102 : // initialize the grid
103 1 : addGrid(gmin,gmax,nbin);
104 1 : }
105 :
106 15 : void BiasRepresentation::addGrid(const std::vector<std::string> & gmin, const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin ) {
107 15 : plumed_massert(nkernels==0,"you can set the grid before loading the hills");
108 15 : plumed_massert(hasgrid==false,"to build the grid you should not having the grid in this bias representation");
109 : std::string ss;
110 : ss="file.free";
111 15 : std::vector<Value*> vv;
112 41 : for(unsigned i=0; i<values.size(); i++) {
113 26 : vv.push_back(values[i]);
114 : }
115 30 : BiasGrid_=Tools::make_unique<Grid>(ss,vv,gmin,gmax,nbin,false,true);
116 15 : hasgrid=true;
117 : // geometry of the grid, used in addSeparableGaussian
118 30 : gridnbin=BiasGrid_->getNbin();
119 30 : gridpbc=BiasGrid_->getIsPeriodic();
120 30 : griddx=BiasGrid_->getDx();
121 15 : gridmin=BiasGrid_->getPoint(std::vector<unsigned>(ndim,0));
122 15 : gridstride.assign(ndim,1);
123 26 : for(int d=1; d<ndim; d++) {
124 11 : gridstride[d]=gridstride[d-1]*gridnbin[d-1];
125 : }
126 15 : sepIndex.resize(ndim);
127 15 : sepT.resize(ndim);
128 15 : sepT2.resize(ndim);
129 15 : sepExp.resize(ndim);
130 30 : }
131 :
132 6143 : bool BiasRepresentation::hasSigmaInInput() {
133 6143 : if(histosigma.size()==0) {
134 : return false;
135 : } else {
136 1092 : return true;
137 : }
138 : }
139 :
140 2 : void BiasRepresentation::setRescaledToBias(bool rescaled) {
141 2 : plumed_massert(nkernels==0,"you can set the rescaling function only before loading hills");
142 2 : rescaledToBias=rescaled;
143 2 : }
144 :
145 0 : const bool & BiasRepresentation::isRescaledToBias() {
146 0 : return rescaledToBias;
147 : }
148 :
149 0 : unsigned BiasRepresentation::getNumberOfDimensions() {
150 0 : return values.size();
151 : }
152 :
153 0 : std::vector<std::string> BiasRepresentation::getNames() {
154 0 : return names;
155 : }
156 :
157 0 : const std::string & BiasRepresentation::getName(unsigned i) {
158 0 : return names[i];
159 : }
160 :
161 0 : const std::vector<Value*>& BiasRepresentation::getPtrToValues() {
162 0 : return values;
163 : }
164 :
165 0 : Value* BiasRepresentation::getPtrToValue(unsigned i) {
166 0 : return values[i];
167 : }
168 :
169 1092 : std::unique_ptr<KernelFunctions> BiasRepresentation::readFromPoint(IFile *ifile) {
170 1092 : std::vector<double> cc( names.size() );
171 3276 : for(unsigned i=0; i<names.size(); ++i) {
172 2184 : ifile->scanField(names[i],cc[i]);
173 : }
174 1092 : double h=1.0;
175 2184 : return Tools::make_unique<KernelFunctions>(cc,histosigma,"stretched-gaussian","DIAGONAL",h);
176 : }
177 :
178 6023 : void BiasRepresentation::pushKernel( IFile *ifile, double time ) {
179 : double bf;
180 6023 : auto kk=readKernel(ifile,bf);
181 6023 : addKernel(std::move(kk),bf,time);
182 6023 : }
183 :
184 6143 : std::unique_ptr<KernelFunctions> BiasRepresentation::readKernel( IFile *ifile, double & bf ) {
185 6143 : std::unique_ptr<KernelFunctions> kk;
186 : // here below the reading of the kernel is completely hidden
187 6143 : if(histosigma.size()==0) {
188 5051 : ifile->allowIgnoredFields();
189 10102 : kk=KernelFunctions::read(ifile,true,names);
190 : } else {
191 : // when doing histogram assume gaussian with a given diagonal sigma
192 : // and neglect all the rest
193 2184 : kk=readFromPoint(ifile);
194 : }
195 : // the bias factor is not something about the kernels but
196 : // must be stored to keep the bias/free energy duality
197 12286 : if(ifile->FieldExist("biasf")) {
198 : std::string dummy;
199 6083 : ifile->scanField("biasf",dummy);
200 6083 : Tools::convert(dummy,bf);
201 : } else {
202 60 : bf=1.0;
203 : }
204 : // the domain does not pertain to the kernel but to the values here defined
205 : std::string mins,maxs,minv,maxv,mini,maxi;
206 : mins="min_";
207 : maxs="max_";
208 18208 : for(int i=0 ; i<ndim; i++) {
209 12065 : if(values[i]->isPeriodic()) {
210 22536 : ifile->scanField(mins+names[i],minv);
211 22536 : ifile->scanField(maxs+names[i],maxv);
212 : // verify that the domain is correct
213 11268 : values[i]->getDomain(mini,maxi);
214 11268 : plumed_massert(mini==minv,"the input periodicity in hills and in value definition does not match" );
215 11268 : plumed_massert(maxi==maxv,"the input periodicity in hills and in value definition does not match" );
216 : }
217 : }
218 6143 : return kk;
219 0 : }
220 :
221 6143 : void BiasRepresentation::addKernel( std::unique_ptr<KernelFunctions> kk, double bf, double time ) {
222 6143 : nkernels++;
223 6143 : lasttime=std::max(lasttime,time);
224 6143 : if(!hasgrid) {
225 : // without a grid only the support of the kernels is needed (see getMinMaxBin)
226 2368 : std::vector<double> ss;
227 2368 : if(histosigma.size()!=0) {
228 546 : ss=histosigma;
229 : } else {
230 3644 : ss=kk->getContinuousSupport();
231 : }
232 : const std::vector<double> cc=kk->getCenter();
233 : const int ndiv=10; // adjustable parameter: division per support
234 7104 : for(int j=0; j<ndim; j++) {
235 5064 : supmin[j]=std::min(supmin[j],cc[j]-ss[j]);
236 5006 : supmax[j]=std::max(supmax[j],cc[j]+ss[j]);
237 4776 : minbinsize[j]=std::min(minbinsize[j],ss[j]/double(ndiv));
238 : }
239 : return;
240 : }
241 : // with MPI each process adds a subset of the kernels to its own grid,
242 : // the grids are summed in getGridPtr
243 3775 : if(mycomm.Get_size()>1 && (nkernels-1)%mycomm.Get_size()!=mycomm.Get_rank()) {
244 : return;
245 : }
246 3775 : const double f=(rescaledToBias ? biasFactorToBias(bf) : 1.0);
247 3775 : if(!timeaverage || time<avt0) {
248 3745 : addToGrid(*BiasGrid_,*kk,f);
249 30 : } else if(!avhasend || time<=avt1) {
250 : // time average: the hills after t0 are accumulated with weights 1 and t,
251 : // and combined in getGridPtr with the weight (t1-t)/(t1-t0)
252 21 : addToGrid(*avS1_,*kk,f);
253 21 : addToGrid(*avS2_,*kk,f*time);
254 : }
255 : }
256 :
257 3787 : void BiasRepresentation::addToGrid( Grid & g, const KernelFunctions & kk, double f ) {
258 : double h,a,b,r2cut;
259 3787 : std::vector<double> w;
260 3787 : if(!doInt_ && kk.getSeparableGaussian(h,a,b,r2cut,w)) {
261 1057 : addSeparableGaussian(g,kk,f*h,a,b,r2cut,w);
262 : } else {
263 2730 : addGenericKernel(g,kk,f);
264 : }
265 3787 : }
266 :
267 2730 : void BiasRepresentation::addGenericKernel( Grid & g, const KernelFunctions & kk, double f ) {
268 2730 : std::vector<unsigned> nneighb;
269 2730 : if(doInt_&&(kk.getCenter()[0]+kk.getContinuousSupport()[0] > uppI_ || kk.getCenter()[0]-kk.getContinuousSupport()[0] < lowI_ )) {
270 0 : nneighb=g.getNbin();
271 : } else {
272 5460 : nneighb=kk.getSupport(g.getDx());
273 : }
274 2730 : std::vector<Grid::index_t> neighbors=g.getNeighbors(kk.getCenter(),nneighb);
275 2730 : std::vector<double> der(ndim);
276 2730 : std::vector<double> xx(ndim);
277 973234 : for(unsigned i=0; i<neighbors.size(); ++i) {
278 970504 : Grid::index_t ineigh=neighbors[i];
279 970504 : g.getPoint(ineigh,xx);
280 2911512 : for(int j=0; j<ndim; ++j) {
281 1941008 : values[j]->set(xx[j]);
282 : }
283 : double bias;
284 970504 : if(doInt_) {
285 0 : bias=kk.evaluate(values,der,true,doInt_,lowI_,uppI_);
286 : } else {
287 970504 : bias=kk.evaluate(values,der,true);
288 : }
289 970504 : bias*=f;
290 2911512 : for(int j=0; j<ndim; ++j) {
291 1941008 : der[j]*=f;
292 : }
293 970504 : g.addValueAndDerivatives(ineigh,bias,der);
294 : }
295 5460 : }
296 :
297 1057 : void BiasRepresentation::addSeparableGaussian( Grid & g, const KernelFunctions & kk, double h, double a, double b, double r2cut, const std::vector<double> & w ) {
298 : // The kernel is evaluated on the same grid points as in addGenericKernel, but the
299 : // Gaussian factors are computed once per dimension and then multiplied, instead of
300 : // computing an exponential at each point of the grid.
301 : const std::vector<double> c=kk.getCenter();
302 1057 : const std::vector<unsigned> nneigh=kk.getSupport(griddx);
303 2938 : for(int d=0; d<ndim; d++) {
304 1881 : auto & idx=sepIndex[d];
305 : auto & t=sepT[d];
306 : auto & q=sepT2[d];
307 : auto & e=sepExp[d];
308 1881 : idx.clear();
309 1881 : t.clear();
310 1881 : q.clear();
311 1881 : e.clear();
312 1881 : const int nb=gridnbin[d];
313 1881 : const int ic=static_cast<int>(std::floor((c[d]-gridmin[d])/griddx[d]));
314 1881 : int span=2*nneigh[d]+1;
315 : // along a periodic direction visit each point of the grid only once
316 1881 : if(gridpbc[d] && span>nb) {
317 : span=nb;
318 : }
319 23984 : for(int k=0; k<span; k++) {
320 22103 : int i=ic-static_cast<int>(nneigh[d])+k;
321 22103 : if(!gridpbc[d] && (i<0 || i>=nb)) {
322 70 : continue;
323 : }
324 22033 : if(gridpbc[d]) {
325 15476 : i=((i%nb)+nb)%nb;
326 : }
327 22033 : const double x=gridmin[d]+static_cast<double>(i)*griddx[d];
328 22033 : const double tt=-values[d]->difference(x,c[d])/w[d];
329 22033 : idx.push_back(i);
330 22033 : t.push_back(tt);
331 22033 : q.push_back(tt*tt);
332 22033 : e.push_back(std::exp(-0.5*tt*tt));
333 : }
334 1881 : if(idx.empty()) {
335 : return;
336 : }
337 : }
338 : // prefix (over the dimensions from the last to the second) of r2, of the product of the exponentials and of the index
339 1057 : std::vector<unsigned> pos(ndim,0);
340 1057 : std::vector<double> pr2(ndim+1,0.0), pe(ndim+1,1.0);
341 1057 : std::vector<Grid::index_t> pidx(ndim+1,0);
342 8300 : auto setLevel=[&](int l) {
343 8300 : pr2[l]=pr2[l+1]+sepT2[l][pos[l]];
344 8300 : pe[l]=pe[l+1]*sepExp[l][pos[l]];
345 8300 : pidx[l]=pidx[l+1]+static_cast<Grid::index_t>(sepIndex[l][pos[l]])*gridstride[l];
346 9357 : };
347 1881 : for(int l=ndim-1; l>=1; --l) {
348 824 : setLevel(l);
349 : }
350 1057 : std::vector<double> der(ndim), dd(ndim);
351 : const auto & idx0=sepIndex[0];
352 : const auto & t0=sepT[0];
353 : const auto & q0=sepT2[0];
354 : const auto & e0=sepExp[0];
355 1057 : const double ha=h*a;
356 1057 : const double hb=h*b;
357 : while(true) {
358 8533 : const double r2u=pr2[1];
359 8533 : const double eu=pe[1];
360 8533 : const Grid::index_t iu=pidx[1];
361 16833 : for(int l=1; l<ndim; ++l) {
362 8300 : der[l]=sepT[l][pos[l]]/w[l];
363 : }
364 112150 : for(unsigned k=0; k<idx0.size(); ++k) {
365 103617 : if(r2u+q0[k]>=r2cut) {
366 47641 : continue;
367 : }
368 55976 : const double ee=eu*e0[k];
369 55976 : const double kderiv=-ha*ee;
370 55976 : dd[0]=t0[k]/w[0]*kderiv;
371 107532 : for(int l=1; l<ndim; ++l) {
372 51556 : dd[l]=der[l]*kderiv;
373 : }
374 55976 : g.addValueAndDerivatives(iu+idx0[k],ha*ee+hb,dd.data(),ndim);
375 : }
376 : // next combination of the indices of the other dimensions
377 : int l=1;
378 9357 : while(l<ndim) {
379 8300 : pos[l]++;
380 8300 : if(pos[l]<sepIndex[l].size()) {
381 : break;
382 : }
383 824 : pos[l]=0;
384 824 : l++;
385 : }
386 8533 : if(l>=ndim) {
387 : break;
388 : }
389 14952 : for(int m=l; m>=1; --m) {
390 7476 : setLevel(m);
391 : }
392 : }
393 1057 : }
394 :
395 6198 : int BiasRepresentation::getNumberOfKernels() {
396 6198 : return nkernels;
397 : }
398 :
399 23 : Grid* BiasRepresentation::getGridPtr() {
400 23 : plumed_massert(hasgrid,"if you want the grid pointer then you should have defined a grid before");
401 23 : if(mycomm.Get_size()==1 && !timeaverage) {
402 21 : return BiasGrid_.get();
403 : }
404 : // sum the grids of all the processes and combine the grids of the time average
405 2 : const auto n=BiasGrid_->getSize();
406 6 : auto getData=[&](const Grid & g, std::vector<double> & val, std::vector<double> & der) {
407 6 : val.resize(n);
408 6 : der.resize(n*ndim);
409 366 : for(Grid::index_t i=0; i<n; i++) {
410 360 : val[i]=g.getValueAndDerivatives(i,&der[i*ndim],ndim);
411 : }
412 6 : if(mycomm.Get_size()>1) {
413 0 : mycomm.Sum(val);
414 0 : mycomm.Sum(der);
415 : }
416 8 : };
417 2 : std::vector<double> val,der;
418 2 : getData(*BiasGrid_,val,der);
419 2 : if(timeaverage) {
420 2 : std::vector<double> val1,der1,val2,der2;
421 2 : getData(*avS1_,val1,der1);
422 2 : getData(*avS2_,val2,der2);
423 2 : const double t1=getTimeAverageEnd();
424 : // weight (t1-t)/(t1-t0) for the hills between t0 and t1
425 2 : const double c1=(t1>avt0 ? t1/(t1-avt0) : 1.0);
426 2 : const double c2=(t1>avt0 ? -1.0/(t1-avt0) : 0.0);
427 122 : for(Grid::index_t i=0; i<n; i++) {
428 120 : val[i]+=c1*val1[i]+c2*val2[i];
429 : }
430 122 : for(Grid::index_t i=0; i<n*ndim; i++) {
431 120 : der[i]+=c1*der1[i]+c2*der2[i];
432 : }
433 : }
434 2 : reducedGrid_=Tools::make_unique<Grid>(*BiasGrid_);
435 2 : std::vector<double> dd(ndim);
436 122 : for(Grid::index_t i=0; i<n; i++) {
437 240 : for(int j=0; j<ndim; j++) {
438 120 : dd[j]=der[i*ndim+j];
439 : }
440 120 : reducedGrid_->setValueAndDerivatives(i,val[i],dd);
441 : }
442 : return reducedGrid_.get();
443 : }
444 :
445 1 : void BiasRepresentation::setTimeAverage(double t0, bool hasend, double t1) {
446 1 : plumed_massert(hasgrid,"the time average requires a grid");
447 1 : plumed_massert(nkernels==0,"the time average should be set before adding the kernels");
448 1 : timeaverage=true;
449 1 : avt0=t0;
450 1 : avhasend=hasend;
451 1 : avt1=t1;
452 1 : plumed_massert(!hasend || t1>t0,"the end of the time average should be after its beginning");
453 1 : avS1_=Tools::make_unique<Grid>(*BiasGrid_);
454 1 : avS2_=Tools::make_unique<Grid>(*BiasGrid_);
455 1 : }
456 :
457 2 : double BiasRepresentation::getTimeAverageEnd() const {
458 : // without an explicit end (or before reaching it) the average ends at the last hill
459 2 : if(!avhasend) {
460 0 : return lasttime;
461 : }
462 3 : return std::min(avt1,lasttime);
463 : }
464 :
465 3 : double BiasRepresentation::getLastTime() const {
466 3 : return lasttime;
467 : }
468 :
469 8 : void BiasRepresentation::getMinMaxBin(std::vector<double> &vmin, std::vector<double> &vmax, std::vector<unsigned> &vbin) {
470 8 : vmin=supmin;
471 8 : vmax=supmax;
472 8 : vbin.assign(ndim,0);
473 24 : for(int j=0; j<ndim; j++) {
474 : // reset to periodicity
475 16 : if(values[j]->isPeriodic()) {
476 : double minv,maxv;
477 8 : values[j]->getDomain(minv,maxv);
478 8 : if(minv>vmin[j]) {
479 0 : vmin[j]=minv;
480 : }
481 8 : if(maxv<vmax[j]) {
482 0 : vmax[j]=maxv;
483 : }
484 : }
485 16 : vbin[j]=static_cast<unsigned>(std::ceil((vmax[j]-vmin[j])/minbinsize[j]) );
486 : }
487 8 : }
488 :
489 0 : void BiasRepresentation::clear() {
490 0 : nkernels=0;
491 0 : lasttime=-std::numeric_limits<double>::max();
492 0 : resetSupport();
493 : // clear the grid
494 0 : if(hasgrid) {
495 0 : BiasGrid_->clear();
496 : }
497 0 : if(timeaverage) {
498 0 : avS1_->clear();
499 0 : avS2_->clear();
500 : }
501 0 : }
502 :
503 23 : void BiasRepresentation::resetSupport() {
504 23 : supmin.assign(ndim,10.e20);
505 23 : supmax.assign(ndim,-10.e20);
506 23 : minbinsize.assign(ndim,10.e20);
507 23 : }
508 :
509 :
510 : }
|
