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 : #ifndef __PLUMED_bias_MetaDKernel_h
23 : #define __PLUMED_bias_MetaDKernel_h
24 :
25 : #include "core/Value.h"
26 : #include "tools/Grid.h"
27 : #include "tools/Matrix.h"
28 : #include "tools/Tools.h"
29 :
30 : #include <cmath>
31 : #include <vector>
32 :
33 : namespace PLMD {
34 : namespace bias {
35 :
36 : /// Work arrays of the OpenMP threads, written at every hill or grid point. When two threads write in the
37 : /// same cache line (false sharing), or one writes where another one reads, both become much slower, and
38 : /// using more threads can take longer than using one: this happened with small vectors allocated by
39 : /// each thread, depending on where the heap placed them. Here the arrays of all the threads are in a
40 : /// single allocation, separated from each other and from the ends of the allocation by 128 bytes that
41 : /// are never used.
42 : template<typename T>
43 26025 : class ThreadWorkArrays {
44 : static constexpr unsigned gap=128/sizeof(T);
45 : unsigned stride;
46 : std::vector<T> buffer;
47 : public:
48 63742 : ThreadWorkArrays(unsigned nthreads, unsigned size):
49 63742 : stride(size+gap),
50 63742 : buffer(gap+nthreads*stride,T()) {}
51 : T* get(unsigned thread) {
52 : return buffer.data()+gap+thread*stride;
53 : }
54 : };
55 :
56 : /// A Gaussian hill deposited by METAD or PBMETAD.
57 : /// For multivariate hills (ADAPTIVE) sigma contains the upper triangle (row-major) of the
58 : /// inverse covariance matrix, otherwise it contains the widths along each CV.
59 : struct MetaDGaussian {
60 : bool multivariate; // this is required to discriminate the one dimensional case
61 : double height;
62 : std::vector<double> center;
63 : std::vector<double> sigma;
64 : std::vector<double> invsigma;
65 : /// multivariate hills: the full inverse covariance matrix (row-major), built from sigma
66 : std::vector<double> metric;
67 24368 : MetaDGaussian(const bool m,
68 : const double h,
69 : const std::vector<double>& c,
70 24368 : const std::vector<double>& s):
71 24368 : multivariate(m),
72 24368 : height(h),
73 24368 : center(c),
74 24368 : sigma(s),
75 24368 : invsigma(s) {
76 : // to avoid troubles from zero element in flexible hills
77 62532 : for(unsigned i=0; i<invsigma.size(); ++i) {
78 38164 : if(std::abs(invsigma[i])>1.e-20) {
79 38082 : invsigma[i]=1.0/invsigma[i] ;
80 : } else {
81 82 : invsigma[i]=0.0;
82 : }
83 : }
84 24368 : if(multivariate) {
85 : const std::size_t ncv=center.size();
86 4603 : metric.resize(ncv*ncv);
87 : unsigned k=0;
88 11459 : for(unsigned i=0; i<ncv; i++) {
89 15965 : for(unsigned j=i; j<ncv; j++) {
90 9109 : metric[i*ncv+j]=metric[j*ncv+i]=sigma[k];
91 9109 : k++;
92 : }
93 : }
94 : }
95 24368 : }
96 : };
97 :
98 : /// Evaluation of the (possibly stretched) Gaussian hills of METAD as a function of its arguments.
99 : /// Gaussians are truncated at dp2cutoff. With INTERVAL (one dimensional only) the CV is
100 : /// clamped to [lower,upper] and the force outside the interval is zero.
101 255 : class MetaDKernel {
102 : /// periodicity of the arguments, cached to compute the differences inline
103 : /// (same arithmetic as Value::difference)
104 : std::vector<char> periodic_;
105 : std::vector<double> period_;
106 : std::vector<double> invperiod_;
107 : double stretchA_=1.0;
108 : double stretchB_=0.0;
109 : bool interval_=false;
110 : double lower_=-1;
111 : double upper_=-1;
112 : public:
113 : /// d2-d1, with the periodicity of argument i (as Value::difference)
114 7703574 : double difference(unsigned i,double d1,double d2) const {
115 7703574 : if(!periodic_[i]) {
116 5249373 : return d2-d1;
117 : }
118 2454201 : double s=(d2-d1)*invperiod_[i];
119 2454201 : s=Tools::pbc(s);
120 2454201 : return s*period_[i];
121 : }
122 450 : void setArguments(const std::vector<Value*>& args) {
123 450 : periodic_.resize(args.size());
124 450 : period_.assign(args.size(),0.0);
125 450 : invperiod_.assign(args.size(),0.0);
126 1094 : for(unsigned i=0; i<args.size(); i++) {
127 644 : periodic_[i]=args[i]->isPeriodic();
128 644 : if(periodic_[i]) {
129 194 : period_[i]=args[i]->getMaxMinusMin();
130 194 : invperiod_[i]=1.0/period_[i];
131 : }
132 : }
133 450 : }
134 : /// stretched Gaussians: A*exp(-dp2)+B, so that they go to zero at the cutoff
135 : void setStretch(double A,double B) {
136 450 : stretchA_=A;
137 450 : stretchB_=B;
138 450 : }
139 : void setInterval(double lower,double upper) {
140 19 : interval_=true;
141 19 : lower_=lower;
142 19 : upper_=upper;
143 19 : }
144 : /// value of hill in cv
145 : double evaluate(const std::vector<double>& cv, const MetaDGaussian& hill) const;
146 : /// value of hill in cv; the derivatives are subtracted to der; dp_ is a work array of size ncv
147 : double evaluateWithDerivatives(const std::vector<double>& cv, const MetaDGaussian& hill, double* der, double* dp_) const;
148 : /// value of hill at a point, given the differences dp of the point from the center of the hill
149 : /// (multiplied by invsigma for diagonal hills, plain differences for multivariate hills);
150 : /// the derivatives are subtracted to der. Outside the cutoff it returns 0, leaves der unchanged and sets
151 : /// inside=false. Same arithmetic as evaluateWithDerivatives, INTERVAL is not taken into account.
152 : double evaluateFromDifferences(const double* dp, const MetaDGaussian& hill, double* der, bool& inside) const;
153 : /// true if INTERVAL is used
154 : bool hasInterval() const {
155 2565073 : return interval_;
156 : }
157 : /// as evaluate, for a diagonal hill given by pointers to its center and inverse widths (no INTERVAL)
158 : double evaluateDiagonal(const std::vector<double>& cv, const double* center, const double* invsigma, double height) const;
159 : /// as evaluateWithDerivatives, for a diagonal hill given by pointers to its center and inverse widths (no INTERVAL)
160 : double evaluateDiagonalWithDerivatives(const std::vector<double>& cv, const double* center, const double* invsigma, double height,
161 : double* der, double* dp_) const;
162 : /// half-width of the box containing the hill (truncated at dp2cutoff) along each CV; it can be infinite for a degenerate hill
163 : std::vector<double> getCutoff(const MetaDGaussian& hill) const;
164 : /// number of grid points along each direction that are within the cutoff of hill
165 : std::vector<unsigned> getSupport(const MetaDGaussian& hill, const GridBase& grid) const;
166 : /// squared distance of cv from the center of hill, in units of the hill widths (d^T*S^-1*d, S covariance);
167 : /// the hill is truncated where this is 2*dp2cutoff. INTERVAL is not taken into account.
168 : double getDistance2(const std::vector<double>& cv, const MetaDGaussian& hill) const;
169 : /// variance of hill along CV i (diagonal of the covariance matrix for multivariate hills)
170 : double getVariance(const MetaDGaussian& hill, unsigned i) const;
171 : };
172 :
173 : inline
174 95663 : double MetaDKernel::evaluate(const std::vector<double>& cv, const MetaDGaussian& hill) const {
175 95663 : unsigned ncv=cv.size();
176 :
177 : // with INTERVAL the CV is replaced by the limit when it is beyond it (cv is const)
178 : double tmpcv[1];
179 : const double *pcv=nullptr; // pointer to cv
180 95663 : if(ncv>0) {
181 : pcv=&cv[0];
182 : }
183 95663 : if(interval_) {
184 24086 : plumed_assert(ncv==1);
185 24086 : tmpcv[0]=cv[0];
186 24086 : if(cv[0]<lower_) {
187 611 : tmpcv[0]=lower_;
188 : }
189 24086 : if(cv[0]>upper_) {
190 17776 : tmpcv[0]=upper_;
191 : }
192 : pcv=&(tmpcv[0]);
193 : }
194 :
195 : double dp2=0.0;
196 95663 : if(hill.multivariate) {
197 : // full inverse covariance matrix
198 : const double* mymatrix=hill.metric.data();
199 150693 : for(unsigned i=0; i<ncv; i++) {
200 79116 : double dp_i=difference(i,hill.center[i],pcv[i]);
201 165771 : for(unsigned j=i; j<ncv; j++) {
202 86655 : if(i==j) {
203 79116 : dp2+=dp_i*dp_i*mymatrix[i*ncv+j]*0.5;
204 : } else {
205 7539 : double dp_j=difference(j,hill.center[j],pcv[j]);
206 7539 : dp2+=dp_i*dp_j*mymatrix[i*ncv+j];
207 : }
208 : }
209 : }
210 : } else {
211 48172 : for(unsigned i=0; i<ncv; i++) {
212 24086 : double dp=difference(i,hill.center[i],pcv[i])*hill.invsigma[i];
213 24086 : dp2+=dp*dp;
214 : }
215 24086 : dp2*=0.5;
216 : }
217 :
218 : double bias=0.0;
219 95663 : if(dp2<dp2cutoff) {
220 55993 : bias=hill.height*(stretchA_*std::exp(-dp2)+stretchB_);
221 : }
222 :
223 95663 : return bias;
224 : }
225 :
226 : inline
227 911295 : double MetaDKernel::evaluateWithDerivatives(const std::vector<double>& cv, const MetaDGaussian& hill, double* der, double* dp_) const {
228 911295 : unsigned ncv=cv.size();
229 :
230 : // with INTERVAL the CV is replaced by the limit when it is beyond it (cv is const)
231 : const double *pcv=nullptr; // pointer to cv
232 : double tmpcv[1]; // tmp array with cv (to be used with INTERVAL)
233 911295 : if(ncv>0) {
234 : pcv=&cv[0];
235 : }
236 911295 : if(interval_) {
237 188110 : plumed_assert(ncv==1);
238 188110 : tmpcv[0]=cv[0];
239 188110 : if(cv[0]<lower_) {
240 29663 : tmpcv[0]=lower_;
241 : }
242 188110 : if(cv[0]>upper_) {
243 119374 : tmpcv[0]=upper_;
244 : }
245 : pcv=&(tmpcv[0]);
246 : }
247 :
248 : bool int_der=false;
249 911295 : if(interval_) {
250 188110 : if(cv[0]<lower_ || cv[0]>upper_) {
251 : int_der=true;
252 : }
253 : }
254 :
255 : double dp2=0.0;
256 : double bias=0.0;
257 911295 : if(hill.multivariate) {
258 : // full inverse covariance matrix
259 : const double* mymatrix=hill.metric.data();
260 1741523 : for(unsigned i=0; i<ncv; i++) {
261 1013010 : dp_[i]=difference(i,hill.center[i],pcv[i]);
262 2310517 : for(unsigned j=i; j<ncv; j++) {
263 1297507 : if(i==j) {
264 1013010 : dp2+=dp_[i]*dp_[i]*mymatrix[i*ncv+j]*0.5;
265 : } else {
266 284497 : double dp_j=difference(j,hill.center[j],pcv[j]);
267 284497 : dp2+=dp_[i]*dp_j*mymatrix[i*ncv+j];
268 : }
269 : }
270 : }
271 728513 : if(dp2<dp2cutoff) {
272 459804 : bias=hill.height*std::exp(-dp2);
273 459804 : if(!int_der) {
274 1133968 : for(unsigned i=0; i<ncv; i++) {
275 : double tmp=0.0;
276 1777048 : for(unsigned j=0; j<ncv; j++) {
277 1102884 : tmp += dp_[j]*mymatrix[i*ncv+j]*bias;
278 : }
279 674164 : der[i]-=tmp*stretchA_;
280 : }
281 : } else {
282 0 : for(unsigned i=0; i<ncv; i++) {
283 0 : der[i]=0.;
284 : }
285 : }
286 459804 : bias=stretchA_*bias+hill.height*stretchB_;
287 : }
288 : } else {
289 365564 : for(unsigned i=0; i<ncv; i++) {
290 182782 : dp_[i]=difference(i,hill.center[i],pcv[i])*hill.invsigma[i];
291 182782 : dp2+=dp_[i]*dp_[i];
292 : }
293 182782 : dp2*=0.5;
294 182782 : if(dp2<dp2cutoff) {
295 159844 : bias=hill.height*std::exp(-dp2);
296 159844 : if(!int_der) {
297 63898 : for(unsigned i=0; i<ncv; i++) {
298 31949 : der[i]-=bias*dp_[i]*hill.invsigma[i]*stretchA_;
299 : }
300 : } else {
301 255790 : for(unsigned i=0; i<ncv; i++) {
302 127895 : der[i]=0.;
303 : }
304 : }
305 159844 : bias=stretchA_*bias+hill.height*stretchB_;
306 : }
307 : }
308 :
309 911295 : return bias;
310 : }
311 :
312 : inline
313 24826 : std::vector<double> MetaDKernel::getCutoff(const MetaDGaussian& hill) const {
314 24826 : std::vector<double> cutoff;
315 24826 : unsigned ncv=periodic_.size();
316 :
317 : // traditional or flexible hill?
318 24826 : if(hill.multivariate) {
319 : unsigned k=0;
320 1900 : Matrix<double> mymatrix(ncv,ncv);
321 4644 : for(unsigned i=0; i<ncv; i++) {
322 6332 : for(unsigned j=i; j<ncv; j++) {
323 : // recompose the full inverse matrix
324 3588 : mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k];
325 3588 : k++;
326 : }
327 : }
328 : // the covariance matrix is the inverse of the matrix stored in the hill
329 1900 : Matrix<double> mycov(ncv,ncv);
330 1900 : Invert(mymatrix,mycov);
331 : // the hill is truncated on the ellipsoid 0.5*d^T*mymatrix*d = dp2cutoff,
332 : // whose bounding box has half-width sqrt(2*dp2cutoff*cov_ii) along CV i
333 4644 : for(unsigned i=0; i<ncv; i++) {
334 2744 : cutoff.push_back(std::sqrt(2.0*dp2cutoff)*std::sqrt(mycov(i,i)));
335 : }
336 : } else {
337 53986 : for(unsigned i=0; i<ncv; ++i) {
338 31060 : cutoff.push_back(std::sqrt(2.0*dp2cutoff)*hill.sigma[i]);
339 : }
340 : }
341 24826 : return cutoff;
342 : }
343 :
344 : inline
345 12413 : std::vector<unsigned> MetaDKernel::getSupport(const MetaDGaussian& hill, const GridBase& grid) const {
346 12413 : std::vector<unsigned> nneigh;
347 12413 : const std::vector<double> cutoff=getCutoff(hill);
348 12413 : unsigned ncv=periodic_.size();
349 :
350 : // number of grid points within the cutoff along CV i; it is never larger than the grid
351 : // (a support larger than the grid, or not finite for a degenerate hill, covers the whole grid)
352 12413 : const std::vector<unsigned> nbin=grid.getNbin();
353 12413 : const std::vector<double> dx=grid.getDx();
354 16325 : auto npoints=[&](unsigned i) {
355 16325 : const double n=std::ceil(cutoff[i]/dx[i]);
356 16325 : return (n>=0.0 && n<nbin[i]) ? static_cast<unsigned>(n) : nbin[i];
357 12413 : };
358 12413 : if(interval_) {
359 721 : if(hill.center[0]+cutoff[0] > upper_ || hill.center[0]-cutoff[0] < lower_) {
360 : // in this case, we updated the entire grid to avoid problems
361 577 : return nbin;
362 : } else {
363 144 : nneigh.push_back(npoints(0));
364 : return nneigh;
365 : }
366 : } else {
367 27873 : for(unsigned i=0; i<ncv; i++) {
368 16181 : nneigh.push_back(npoints(i));
369 : }
370 : }
371 :
372 : return nneigh;
373 : }
374 :
375 : inline
376 92905 : double MetaDKernel::getDistance2(const std::vector<double>& cv, const MetaDGaussian& hill) const {
377 92905 : const unsigned ncv=periodic_.size();
378 : double dist2=0.0;
379 92905 : if(hill.multivariate) {
380 : // hill.sigma holds the upper triangle of the inverse covariance matrix
381 41015 : std::vector<double> d(ncv);
382 123045 : for(unsigned i=0; i<ncv; i++) {
383 82030 : d[i]=difference(i,cv[i],hill.center[i]);
384 : }
385 : unsigned k=0;
386 123045 : for(unsigned i=0; i<ncv; i++) {
387 205075 : for(unsigned j=i; j<ncv; j++) {
388 164060 : dist2+=(i==j ? 1.0 : 2.0)*d[i]*d[j]*hill.sigma[k];
389 123045 : k++;
390 : }
391 : }
392 : } else {
393 146074 : for(unsigned i=0; i<ncv; i++) {
394 94184 : const double d=difference(i,cv[i],hill.center[i])/hill.sigma[i];
395 94184 : dist2+=d*d;
396 : }
397 : }
398 92905 : return dist2;
399 : }
400 :
401 : inline
402 0 : double MetaDKernel::getVariance(const MetaDGaussian& hill, unsigned i) const {
403 0 : if(!hill.multivariate) {
404 0 : return hill.sigma[i]*hill.sigma[i];
405 : }
406 0 : const unsigned ncv=periodic_.size();
407 0 : Matrix<double> mymatrix(ncv,ncv);
408 : unsigned k=0;
409 0 : for(unsigned j=0; j<ncv; j++) {
410 0 : for(unsigned l=j; l<ncv; l++) {
411 0 : mymatrix(j,l)=mymatrix(l,j)=hill.sigma[k];
412 0 : k++;
413 : }
414 : }
415 0 : Matrix<double> mycov(ncv,ncv);
416 0 : Invert(mymatrix,mycov);
417 0 : return mycov(i,i);
418 : }
419 :
420 : inline
421 6973135 : double MetaDKernel::evaluateFromDifferences(const double* dp, const MetaDGaussian& hill, double* der, bool& inside) const {
422 6973135 : const unsigned ncv=hill.center.size();
423 : double dp2=0.0;
424 : double bias=0.0;
425 6973135 : if(hill.multivariate) {
426 : const double* mymatrix=hill.metric.data();
427 3724496 : for(unsigned i=0; i<ncv; i++) {
428 6143408 : for(unsigned j=i; j<ncv; j++) {
429 3676432 : if(i==j) {
430 2466976 : dp2+=dp[i]*dp[i]*mymatrix[i*ncv+j]*0.5;
431 : } else {
432 1209456 : dp2+=dp[i]*dp[j]*mymatrix[i*ncv+j];
433 : }
434 : }
435 : }
436 1257520 : inside=(dp2<dp2cutoff);
437 1257520 : if(inside) {
438 924970 : bias=hill.height*std::exp(-dp2);
439 2727648 : for(unsigned i=0; i<ncv; i++) {
440 : double tmp=0.0;
441 5360772 : for(unsigned j=0; j<ncv; j++) {
442 3558094 : tmp += dp[j]*mymatrix[i*ncv+j]*bias;
443 : }
444 1802678 : der[i]-=tmp*stretchA_;
445 : }
446 924970 : bias=stretchA_*bias+hill.height*stretchB_;
447 : }
448 : } else {
449 16899376 : for(unsigned i=0; i<ncv; i++) {
450 11183761 : dp2+=dp[i]*dp[i];
451 : }
452 5715615 : dp2*=0.5;
453 5715615 : inside=(dp2<dp2cutoff);
454 5715615 : if(inside) {
455 4016953 : bias=hill.height*std::exp(-dp2);
456 11820652 : for(unsigned i=0; i<ncv; i++) {
457 7803699 : der[i]-=bias*dp[i]*hill.invsigma[i]*stretchA_;
458 : }
459 4016953 : bias=stretchA_*bias+hill.height*stretchB_;
460 : }
461 : }
462 6973135 : return bias;
463 : }
464 :
465 : inline
466 69370 : double MetaDKernel::evaluateDiagonal(const std::vector<double>& cv, const double* center, const double* invsigma, double height) const {
467 69370 : const unsigned ncv=cv.size();
468 : double dp2=0.0;
469 141252 : for(unsigned i=0; i<ncv; i++) {
470 71882 : double dp=difference(i,center[i],cv[i])*invsigma[i];
471 71882 : dp2+=dp*dp;
472 : }
473 69370 : dp2*=0.5;
474 : double bias=0.0;
475 69370 : if(dp2<dp2cutoff) {
476 52353 : bias=height*(stretchA_*std::exp(-dp2)+stretchB_);
477 : }
478 69370 : return bias;
479 : }
480 :
481 : inline
482 2769520 : double MetaDKernel::evaluateDiagonalWithDerivatives(const std::vector<double>& cv, const double* center, const double* invsigma, double height,
483 : double* der, double* dp_) const {
484 2769520 : const unsigned ncv=cv.size();
485 : double dp2=0.0;
486 8006110 : for(unsigned i=0; i<ncv; i++) {
487 5236590 : dp_[i]=difference(i,center[i],cv[i])*invsigma[i];
488 5236590 : dp2+=dp_[i]*dp_[i];
489 : }
490 2769520 : dp2*=0.5;
491 : double bias=0.0;
492 2769520 : if(dp2<dp2cutoff) {
493 1623497 : bias=height*std::exp(-dp2);
494 4646800 : for(unsigned i=0; i<ncv; i++) {
495 3023303 : der[i]-=bias*dp_[i]*invsigma[i]*stretchA_;
496 : }
497 1623497 : bias=stretchA_*bias+height*stretchB_;
498 : }
499 2769520 : return bias;
500 : }
501 :
502 : }
503 : }
504 :
505 : #endif
|
