Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2021 The VES code team
3 : (see the PEOPLE-VES file at the root of this folder for a list of names)
4 :
5 : See http://www.ves-code.org for more information.
6 :
7 : This file is part of VES code module.
8 :
9 : The VES code module 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 : The VES code module 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 the VES code module. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 :
23 : #include "TargetDistribution.h"
24 : #include "GridIntegrationWeights.h"
25 : #include "core/ActionRegister.h"
26 : #include "tools/Grid.h"
27 : #include "core/PlumedMain.h"
28 : #include <cfloat>
29 :
30 :
31 : namespace PLMD {
32 : namespace ves {
33 :
34 : //+PLUMEDOC VES_TARGETDIST TD_MULTICANONICAL
35 : /*
36 : Multicanonical target distribution (dynamic).
37 :
38 : Use the target distribution to sample the multicanonical ensemble \cite Berg-PRL-1992 \cite Piaggi-PRL-2019.
39 : In this way, in a single molecular dynamics simulation one can obtain information about the system in a range of temperatures.
40 : This range is determined through the keywords MIN_TEMP and MAX_TEMP.
41 :
42 : The collective variables (CVs) used to construct the bias potential must be:
43 : 1. the energy or,
44 : 2. the energy and an order parameter.
45 :
46 : Other choices of CVs or a different order of the above mentioned CVs are nonsensical.
47 : The second CV, the order parameter, must be used when one aims at studying a first order phase transition in the chosen temperature interval \cite Piaggi-JCP-2019.
48 :
49 : The algorithm will explore the free energy at each temperature up to a predefined free
50 : energy threshold \f$\epsilon\f$ specified through the keyword THRESHOLD (in kT units).
51 : If only the energy is biased, i.e. no phase transition is considered, then THRESHOLD can be set to around 5.
52 : If also an order parameter is used then the THRESHOLD should be greater than the barrier for the transformation in kT.
53 : For small systems undergoing a freezing transition THRESHOLD is typically between 20 and 50.
54 :
55 : When only the potential energy is used as CV the method is equivalent to the Wang-Landau algorithm \cite wanglandau.
56 : The advantage with respect to Wang-Landau is that instead of sampling the potential energy indiscriminately, an interval is chosen on the fly based on the minimum and maximum targeted temperatures.
57 :
58 : The algorithm works as follows.
59 : The target distribution for the potential energy is chosen to be:
60 :
61 : \f[
62 : p(E)= \left\{\begin{array}{ll}
63 : \frac{1}{E_2-E_1} & \mathrm{if} \quad E_1<E<E_2 \\
64 : 0 & \mathrm{otherwise}
65 : \end{array}\right.
66 : \f]
67 :
68 : where the energy limits \f$E_1\f$ and \f$E_2\f$ are yet to be determined.
69 : Clearly the interval \f$E_1–E_2\f$ chosen is related to the interval of temperatures \f$T_1-T_2\f$.
70 : To link these two intervals we make use of the following relation:
71 : \f[
72 : \beta' F_{\beta'}(E) = \beta F_{\beta}(E) + (\beta' - \beta) E + C,
73 : \f]
74 : where \f$F_{\beta}(E)\f$ is determined during the optimization and we shall choose \f$C\f$ such that \f$F_{\beta'}(E_{m})=0\f$ with \f$E_{m}\f$ the position of the free energy minimum.
75 : Using this relation we employ an iterative procedure to find the energy interval.
76 : At iteration \f$k\f$ we have the estimates \f$E_1^k\f$ and \f$E_2^k\f$ for \f$E_1\f$ and \f$E_2\f$, and the target distribution is:
77 : \f[
78 : p^k(E)=\frac{1}{E_2^k-E_1^k} \quad \mathrm{for} \quad E_1^k<E<E_2^k.
79 : \f]
80 : \f$E_1^k\f$ and \f$E_2^k\f$ are obtained from the leftmost solution of \f$\beta_2 F_{\beta_2}^{k-1}(E_1^k)=\epsilon\f$ and the rightmost solution of \f$\beta_1 F_{\beta_1}^{k-1}(E_2^k)=\epsilon\f$.
81 : The procedure is repeated until convergence.
82 : This iterative approach is similar to that in \ref TD_WELLTEMPERED.
83 :
84 : The version of this algorithm in which the energy and an order parameter are biased is similar to the one described in \ref TD_MULTITHERMAL_MULTIBARIC.
85 :
86 : The output of these simulations can be reweighted in order to obtain information at all temperatures in the targeted temperature interval.
87 : The reweighting can be performed using the action \ref REWEIGHT_TEMP_PRESS.
88 :
89 : \par Examples
90 :
91 : The following input can be used to run a simulation in the multicanonical ensemble.
92 : The temperature interval to be explored is 400-600 K.
93 : The energy is used as collective variable.
94 : Legendre polynomials are used to construct the bias potential.
95 : The averaged stochastic gradient descent algorithm is chosen to optimize the VES functional.
96 : The target distribution is updated every 100 optimization steps (200 ps here) using the last estimation of the free energy.
97 :
98 : \plumedfile
99 : # Use energy and volume as CVs
100 : energy: ENERGY
101 :
102 : # Basis functions
103 : bf1: BF_LEGENDRE ORDER=20 MINIMUM=-25000 MAXIMUM=-23500
104 :
105 : # Target distributions
106 : TD_MULTICANONICAL ...
107 : LABEL=td_multi
108 : MIN_TEMP=400
109 : MAX_TEMP=600
110 : ... TD_MULTICANONICAL
111 :
112 : # Expansion
113 : VES_LINEAR_EXPANSION ...
114 : ARG=energy
115 : BASIS_FUNCTIONS=bf1
116 : TEMP=500.0
117 : GRID_BINS=1000
118 : TARGET_DISTRIBUTION=td_multi
119 : LABEL=b1
120 : ... VES_LINEAR_EXPANSION
121 :
122 : # Optimization algorithm
123 : OPT_AVERAGED_SGD ...
124 : BIAS=b1
125 : STRIDE=500
126 : LABEL=o1
127 : STEPSIZE=1.0
128 : FES_OUTPUT=500
129 : BIAS_OUTPUT=500
130 : TARGETDIST_OUTPUT=500
131 : COEFFS_OUTPUT=10
132 : TARGETDIST_STRIDE=100
133 : ... OPT_AVERAGED_SGD
134 :
135 : \endplumedfile
136 :
137 : The multicanonical target distribution can also be used to explore a temperature interval in which a first order phase transitions is observed.
138 :
139 : */
140 : //+ENDPLUMEDOC
141 :
142 : class TD_Multicanonical: public TargetDistribution {
143 : private:
144 : double threshold_, min_temp_, max_temp_;
145 : std::vector<double> sigma_;
146 : unsigned steps_temp_;
147 : double epsilon_;
148 : bool smoothening_;
149 : public:
150 : static void registerKeywords(Keywords&);
151 : explicit TD_Multicanonical(const ActionOptions& ao);
152 : void updateGrid() override;
153 : double getValue(const std::vector<double>&) const override;
154 4 : ~TD_Multicanonical() {}
155 : double GaussianSwitchingFunc(const double, const double, const double) const;
156 : };
157 :
158 :
159 : PLUMED_REGISTER_ACTION(TD_Multicanonical,"TD_MULTICANONICAL")
160 :
161 :
162 4 : void TD_Multicanonical::registerKeywords(Keywords& keys) {
163 4 : TargetDistribution::registerKeywords(keys);
164 4 : keys.add("compulsory","THRESHOLD","5","Maximum exploration free energy in kT.");
165 4 : keys.add("compulsory","EPSILON","10","The zeros of the target distribution are changed to e^-EPSILON.");
166 4 : keys.add("compulsory","MIN_TEMP","Minimum temperature.");
167 4 : keys.add("compulsory","MAX_TEMP","Maximum temperature.");
168 4 : keys.add("optional","STEPS_TEMP","Number of temperature steps. Only for the 2D version, i.e. energy and order parameter.");
169 4 : keys.add("optional","SIGMA","The standard deviation parameters of the Gaussian kernels used for smoothing the target distribution. One value must be specified for each argument, i.e. one value per CV. A value of 0.0 means that no smoothing is performed, this is the default behavior.");
170 4 : }
171 :
172 :
173 2 : TD_Multicanonical::TD_Multicanonical(const ActionOptions& ao):
174 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
175 2 : threshold_(5.0),
176 2 : min_temp_(0.0),
177 2 : max_temp_(1000.0),
178 4 : sigma_(0.0),
179 2 : steps_temp_(20),
180 2 : epsilon_(10.0),
181 2 : smoothening_(true) {
182 2 : log.printf(" Multicanonical target distribution");
183 2 : log.printf("\n");
184 2 : log.printf(" Please read and cite ");
185 4 : log << plumed.cite("Piaggi and Parrinello, Phys. Rev. Lett. 122 (5), 050601 (2019)");
186 2 : log.printf(" and ");
187 4 : log << plumed.cite("Piaggi and Parrinello, J. Chem. Phys. 150 (24), 244119 (2019)");
188 2 : log.printf("\n");
189 2 : parse("THRESHOLD",threshold_);
190 2 : if(threshold_<=0.0) {
191 0 : plumed_merror(getName()+": the value of the threshold should be positive.");
192 : }
193 2 : log.printf(" exploring free energy up to %f kT for each temperature \n",threshold_);
194 :
195 2 : parse("MIN_TEMP",min_temp_);
196 2 : parse("MAX_TEMP",max_temp_);
197 2 : log.printf(" temperatures between %f and %f will be explored \n",min_temp_,max_temp_);
198 4 : parseVector("SIGMA",sigma_);
199 2 : if(sigma_.size()==0) {
200 0 : smoothening_=false;
201 : }
202 2 : if(smoothening_ && (sigma_.size()<1 || sigma_.size()>2) ) {
203 0 : plumed_merror(getName()+": SIGMA takes 1 or 2 values as input.");
204 : }
205 2 : if (smoothening_) {
206 2 : log.printf(" the target distribution will be smoothed using sigma values");
207 5 : for(unsigned i=0; i<sigma_.size(); ++i) {
208 3 : log.printf(" %f",sigma_[i]);
209 : }
210 2 : log.printf("\n");
211 : }
212 :
213 2 : parse("STEPS_TEMP",steps_temp_); // Only used in the 2D version
214 2 : steps_temp_ += 1;
215 2 : log.printf(" %d steps in temperatures will be employed (if TD is two-dimensional) \n",steps_temp_);
216 :
217 2 : parse("EPSILON",epsilon_);
218 2 : if(epsilon_<=1.0) {
219 0 : plumed_merror(getName()+": the value of epsilon should be greater than 1.");
220 : }
221 2 : log.printf(" the non relevant regions of the target distribution are set to e^-%f \n",epsilon_);
222 :
223 : setDynamic();
224 : setFesGridNeeded();
225 2 : checkRead();
226 2 : }
227 :
228 :
229 0 : double TD_Multicanonical::getValue(const std::vector<double>& argument) const {
230 0 : plumed_merror("getValue not implemented for TD_Multicanonical");
231 : return 0.0;
232 : }
233 :
234 :
235 14 : void TD_Multicanonical::updateGrid() {
236 14 : if (getStep() == 0) {
237 2 : if(targetDistGrid().getDimension()>2 || targetDistGrid().getDimension()<1) {
238 0 : plumed_merror(getName()+" works only with 1 or 2 arguments, i.e. energy, or energy and CV");
239 : }
240 2 : if(smoothening_ && sigma_.size()!=targetDistGrid().getDimension()) {
241 0 : plumed_merror(getName()+": mismatch between SIGMA dimension and number of arguments");
242 : }
243 : // Use uniform TD
244 4 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
245 : double norm = 0.0;
246 2704 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
247 : double value = 1.0;
248 2702 : norm += integration_weights[l]*value;
249 2702 : targetDistGrid().setValue(l,value);
250 : }
251 2 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
252 : } else {
253 : // Two variants: 1D and 2D
254 12 : if(targetDistGrid().getDimension()==1) {
255 : // 1D variant: Multicanonical without order parameter
256 : // In this variant we find the minimum and maximum relevant potential energies.
257 : // Using this information we construct a uniform target distribution in between these two.
258 10 : double beta = getBeta();
259 10 : double beta_prime_min = 1./(getKBoltzmann()*min_temp_);
260 10 : double beta_prime_max = 1./(getKBoltzmann()*max_temp_);
261 10 : plumed_massert(getFesGridPntr()!=NULL,"the FES grid has to be linked to use TD_Multicanonical!");
262 : // Find minimum of F(U) at temperature min
263 : double minval=DBL_MAX;
264 10 : Grid::index_t minindex = (targetDistGrid().getSize())/2;
265 10 : double minpos = targetDistGrid().getPoint(minindex)[0];
266 1020 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
267 1010 : double value = getFesGridPntr()->getValue(l);
268 1010 : double argument = targetDistGrid().getPoint(l)[0];
269 1010 : value = beta*value + (beta_prime_min-beta)*argument;
270 1010 : if(value<minval) {
271 : minval=value;
272 : minpos=argument;
273 : minindex=l;
274 : }
275 : }
276 : // Find minimum energy at low temperature
277 10 : double minimum_low = minpos;
278 11 : for(Grid::index_t l=minindex; l>1; l-=1) {
279 11 : double argument = targetDistGrid().getPoint(l)[0];
280 11 : double argument_next = targetDistGrid().getPoint(l-1)[0];
281 11 : double value = getFesGridPntr()->getValue(l);
282 11 : double value_next = getFesGridPntr()->getValue(l-1);
283 11 : value = beta*value + (beta_prime_min-beta)*argument - minval;
284 11 : value_next = beta*value_next + (beta_prime_min-beta)*argument_next - minval;
285 11 : if (value<threshold_ && value_next>threshold_) {
286 10 : minimum_low = argument_next;
287 10 : break;
288 : }
289 : }
290 : // Find maximum energy at low temperature
291 10 : double maximum_low = minpos;
292 12 : for(Grid::index_t l=minindex; l<(targetDistGrid().getSize()-1); l++) {
293 12 : double argument = targetDistGrid().getPoint(l)[0];
294 12 : double argument_next = targetDistGrid().getPoint(l+1)[0];
295 12 : double value = getFesGridPntr()->getValue(l);
296 12 : double value_next = getFesGridPntr()->getValue(l+1);
297 12 : value = beta*value + (beta_prime_min-beta)*argument - minval;
298 12 : value_next = beta*value_next + (beta_prime_min-beta)*argument_next - minval;
299 12 : if (value<threshold_ && value_next>threshold_) {
300 10 : maximum_low = argument_next;
301 10 : break;
302 : }
303 : }
304 : // Find minimum of F(U) at temperature max
305 : minval=DBL_MAX;
306 1020 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
307 1010 : double value = getFesGridPntr()->getValue(l);
308 1010 : double argument = targetDistGrid().getPoint(l)[0];
309 1010 : value = beta*value + (beta_prime_max-beta)*argument;
310 1010 : if(value<minval) {
311 : minval=value;
312 : minpos=argument;
313 : minindex=l;
314 : }
315 : }
316 : // Find minimum energy at high temperature
317 10 : double minimum_high = minpos;
318 13 : for(Grid::index_t l=minindex; l>1; l-=1) {
319 13 : double argument = targetDistGrid().getPoint(l)[0];
320 13 : double argument_next = targetDistGrid().getPoint(l-1)[0];
321 13 : double value = getFesGridPntr()->getValue(l);
322 13 : double value_next = getFesGridPntr()->getValue(l-1);
323 13 : value = beta*value + (beta_prime_max-beta)*argument - minval;
324 13 : value_next = beta*value_next + (beta_prime_max-beta)*argument_next - minval;
325 13 : if (value<threshold_ && value_next>threshold_) {
326 10 : minimum_high = argument_next;
327 10 : break;
328 : }
329 : }
330 : // Find maximum energy at high temperature
331 10 : double maximum_high = minpos;
332 11 : for(Grid::index_t l=minindex; l<(targetDistGrid().getSize()-1); l++) {
333 11 : double argument = targetDistGrid().getPoint(l)[0];
334 11 : double argument_next = targetDistGrid().getPoint(l+1)[0];
335 11 : double value = getFesGridPntr()->getValue(l);
336 11 : double value_next = getFesGridPntr()->getValue(l+1);
337 11 : value = beta*value + (beta_prime_max-beta)*argument - minval;
338 11 : value_next = beta*value_next + (beta_prime_max-beta)*argument_next - minval;
339 11 : if (value<threshold_ && value_next>threshold_) {
340 10 : maximum_high = argument_next;
341 10 : break;
342 : }
343 : }
344 10 : double minimum = std::min(minimum_low,minimum_high);
345 10 : double maximum = std::max(maximum_low,maximum_high);
346 : // Construct uniform TD in the interval between minimum and maximum
347 20 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
348 : double norm = 0.0;
349 1020 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
350 1010 : double argument = targetDistGrid().getPoint(l)[0];
351 : double value = 1.0;
352 : double tmp;
353 1010 : if(argument < minimum) {
354 217 : if (smoothening_) {
355 217 : tmp = GaussianSwitchingFunc(argument,minimum,sigma_[0]);
356 : } else {
357 0 : tmp = exp(-1.0*epsilon_);
358 : }
359 793 : } else if(argument > maximum) {
360 199 : if (smoothening_) {
361 199 : tmp = GaussianSwitchingFunc(argument,maximum,sigma_[0]);
362 : } else {
363 0 : tmp = exp(-1.0*epsilon_);
364 : }
365 : } else {
366 : tmp = 1.0;
367 : }
368 : value *= tmp;
369 1010 : norm += integration_weights[l]*value;
370 1010 : targetDistGrid().setValue(l,value);
371 : }
372 10 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
373 2 : } else if(targetDistGrid().getDimension()==2) {
374 : // 2D variant: Multicanonical with order parameter
375 : // In this variant we find for each temperature the relevant region of potential energy and order parameter.
376 : // The target distribution will be the union of the relevant regions at all temperatures in the temperature interval.
377 2 : double beta = getBeta();
378 2 : double beta_prime_min = 1./(getKBoltzmann()*min_temp_);
379 2 : double beta_prime_max = 1./(getKBoltzmann()*max_temp_);
380 2 : plumed_massert(getFesGridPntr()!=NULL,"the FES grid has to be linked to use TD_Multicanonical!");
381 : // Set all to zero
382 5204 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
383 5202 : double value = exp(-1.0*epsilon_);
384 5202 : targetDistGrid().setValue(l,value);
385 : }
386 : // Loop over temperatures
387 44 : for(unsigned i=0; i<steps_temp_; i++) {
388 42 : double beta_prime=beta_prime_min + (beta_prime_max-beta_prime_min)*i/(steps_temp_-1);
389 : // Find minimum for this temperature
390 : double minval=DBL_MAX;
391 109284 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
392 109242 : double energy = targetDistGrid().getPoint(l)[0];
393 109242 : double value = getFesGridPntr()->getValue(l);
394 109242 : value = beta*value + (beta_prime-beta)*energy;
395 109242 : if(value<minval) {
396 : minval=value;
397 : }
398 : }
399 : // Now check which energies and volumes are below X kt
400 109284 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
401 109242 : double energy = targetDistGrid().getPoint(l)[0];
402 109242 : double value = getFesGridPntr()->getValue(l);
403 109242 : value = beta*value + (beta_prime-beta)*energy - minval;
404 109242 : if (value<threshold_) {
405 : double value = 1.0;
406 7076 : targetDistGrid().setValue(l,value);
407 : }
408 : }
409 : }
410 2 : if (smoothening_) {
411 2 : std::vector<unsigned> nbin=targetDistGrid().getNbin();
412 2 : std::vector<double> dx=targetDistGrid().getDx();
413 : // Smoothening
414 104 : for(unsigned i=0; i<nbin[0]; i++) {
415 5304 : for(unsigned j=0; j<nbin[1]; j++) {
416 5202 : std::vector<unsigned> indices(2);
417 5202 : indices[0]=i;
418 5202 : indices[1]=j;
419 5202 : Grid::index_t index = targetDistGrid().getIndex(indices);
420 5202 : double energy = targetDistGrid().getPoint(index)[0];
421 5202 : double volume = targetDistGrid().getPoint(index)[1];
422 5202 : double value = targetDistGrid().getValue(index);
423 5202 : if (value>(1-1.e-5)) { // Apply only if this grid point was 1.
424 : // Apply gaussians around
425 773 : std::vector<int> minBin(2), maxBin(2), deltaBin(2); // These cannot be unsigned
426 : // Only consider contributions less than n*sigma bins apart from the actual distance
427 773 : deltaBin[0]=std::floor(6*sigma_[0]/dx[0]);;
428 773 : deltaBin[1]=std::floor(6*sigma_[1]/dx[1]);;
429 : // For energy
430 773 : minBin[0]=i - deltaBin[0];
431 773 : if (minBin[0] < 0) {
432 406 : minBin[0]=0;
433 : }
434 773 : if (minBin[0] > (nbin[0]-1)) {
435 0 : minBin[0]=nbin[0]-1;
436 : }
437 773 : maxBin[0]=i + deltaBin[0];
438 773 : if (maxBin[0] > (nbin[0]-1)) {
439 349 : maxBin[0]=nbin[0]-1;
440 : }
441 : // For volume
442 773 : minBin[1]=j - deltaBin[1];
443 773 : if (minBin[1] < 0) {
444 655 : minBin[1]=0;
445 : }
446 773 : if (minBin[1] > (nbin[1]-1)) {
447 0 : minBin[1]=nbin[1]-1;
448 : }
449 773 : maxBin[1]=j + deltaBin[1];
450 773 : if (maxBin[1] > (nbin[1]-1)) {
451 86 : maxBin[1]=nbin[1]-1;
452 : }
453 31273 : for(unsigned l=minBin[0]; l<maxBin[0]+1; l++) {
454 549973 : for(unsigned m=minBin[1]; m<maxBin[1]+1; m++) {
455 519473 : std::vector<unsigned> indices_prime(2);
456 519473 : indices_prime[0]=l;
457 519473 : indices_prime[1]=m;
458 519473 : Grid::index_t index_prime = targetDistGrid().getIndex(indices_prime);
459 519473 : double energy_prime = targetDistGrid().getPoint(index_prime)[0];
460 519473 : double volume_prime = targetDistGrid().getPoint(index_prime)[1];
461 519473 : double value_prime = targetDistGrid().getValue(index_prime);
462 : // Apply gaussian
463 1558419 : double gaussian_value = GaussianSwitchingFunc(energy_prime,energy,sigma_[0])*GaussianSwitchingFunc(volume_prime,volume,sigma_[1]);
464 519473 : if (value_prime<gaussian_value) {
465 19817 : targetDistGrid().setValue(index_prime,gaussian_value);
466 : }
467 : }
468 : }
469 : }
470 : }
471 : }
472 : }
473 : // Normalize
474 4 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
475 : double norm = 0.0;
476 5204 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
477 5202 : double value = targetDistGrid().getValue(l);
478 5202 : norm += integration_weights[l]*value;
479 : }
480 2 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
481 : } else {
482 0 : plumed_merror(getName()+": Number of arguments for this target distribution must be 1 or 2");
483 : }
484 : }
485 14 : updateLogTargetDistGrid();
486 14 : }
487 :
488 : inline
489 : double TD_Multicanonical::GaussianSwitchingFunc(const double argument, const double center, const double sigma) const {
490 1039362 : if(sigma>0.0) {
491 1039362 : double arg=(argument-center)/sigma;
492 1039362 : return exp(-0.5*arg*arg);
493 : } else {
494 : return 0.0;
495 : }
496 : }
497 :
498 :
499 : }
500 : }
|
