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 8 : keys.add("compulsory","THRESHOLD","5","Maximum exploration free energy in kT.");
165 8 : keys.add("compulsory","EPSILON","10","The zeros of the target distribution are changed to e^-EPSILON.");
166 8 : keys.add("compulsory","MIN_TEMP","Minimum temperature.");
167 8 : keys.add("compulsory","MAX_TEMP","Maximum temperature.");
168 8 : keys.add("optional","STEPS_TEMP","Number of temperature steps. Only for the 2D version, i.e. energy and order parameter.");
169 8 : 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 : {
183 2 : log.printf(" Multicanonical target distribution");
184 2 : log.printf("\n");
185 2 : log.printf(" Please read and cite ");
186 4 : log << plumed.cite("Piaggi and Parrinello, Phys. Rev. Lett. 122 (5), 050601 (2019)");
187 2 : log.printf(" and ");
188 4 : log << plumed.cite("Piaggi and Parrinello, J. Chem. Phys. 150 (24), 244119 (2019)");
189 2 : log.printf("\n");
190 2 : parse("THRESHOLD",threshold_);
191 2 : if(threshold_<=0.0) {
192 0 : plumed_merror(getName()+": the value of the threshold should be positive.");
193 : }
194 2 : log.printf(" exploring free energy up to %f kT for each temperature \n",threshold_);
195 :
196 2 : parse("MIN_TEMP",min_temp_);
197 2 : parse("MAX_TEMP",max_temp_);
198 2 : log.printf(" temperatures between %f and %f will be explored \n",min_temp_,max_temp_);
199 4 : parseVector("SIGMA",sigma_);
200 2 : if(sigma_.size()==0) smoothening_=false;
201 2 : if(smoothening_ && (sigma_.size()<1 || sigma_.size()>2) ) plumed_merror(getName()+": SIGMA takes 1 or 2 values as input.");
202 2 : if (smoothening_) {
203 2 : log.printf(" the target distribution will be smoothed using sigma values");
204 5 : for(unsigned i=0; i<sigma_.size(); ++i) log.printf(" %f",sigma_[i]);
205 2 : log.printf("\n");
206 : }
207 :
208 2 : parse("STEPS_TEMP",steps_temp_); // Only used in the 2D version
209 2 : steps_temp_ += 1;
210 2 : log.printf(" %d steps in temperatures will be employed (if TD is two-dimensional) \n",steps_temp_);
211 :
212 2 : parse("EPSILON",epsilon_);
213 2 : if(epsilon_<=1.0) {
214 0 : plumed_merror(getName()+": the value of epsilon should be greater than 1.");
215 : }
216 2 : log.printf(" the non relevant regions of the target distribution are set to e^-%f \n",epsilon_);
217 :
218 : setDynamic();
219 : setFesGridNeeded();
220 2 : checkRead();
221 2 : }
222 :
223 :
224 0 : double TD_Multicanonical::getValue(const std::vector<double>& argument) const {
225 0 : plumed_merror("getValue not implemented for TD_Multicanonical");
226 : return 0.0;
227 : }
228 :
229 :
230 14 : void TD_Multicanonical::updateGrid() {
231 14 : if (getStep() == 0) {
232 2 : if(targetDistGrid().getDimension()>2 || targetDistGrid().getDimension()<1) plumed_merror(getName()+" works only with 1 or 2 arguments, i.e. energy, or energy and CV");
233 2 : if(smoothening_ && sigma_.size()!=targetDistGrid().getDimension()) plumed_merror(getName()+": mismatch between SIGMA dimension and number of arguments");
234 : // Use uniform TD
235 4 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
236 : double norm = 0.0;
237 2704 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
238 : double value = 1.0;
239 2702 : norm += integration_weights[l]*value;
240 2702 : targetDistGrid().setValue(l,value);
241 : }
242 2 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
243 : } else {
244 : // Two variants: 1D and 2D
245 12 : if(targetDistGrid().getDimension()==1) {
246 : // 1D variant: Multicanonical without order parameter
247 : // In this variant we find the minimum and maximum relevant potential energies.
248 : // Using this information we construct a uniform target distribution in between these two.
249 10 : double beta = getBeta();
250 10 : double beta_prime_min = 1./(getKBoltzmann()*min_temp_);
251 10 : double beta_prime_max = 1./(getKBoltzmann()*max_temp_);
252 10 : plumed_massert(getFesGridPntr()!=NULL,"the FES grid has to be linked to use TD_Multicanonical!");
253 : // Find minimum of F(U) at temperature min
254 : double minval=DBL_MAX;
255 10 : Grid::index_t minindex = (targetDistGrid().getSize())/2;
256 10 : double minpos = targetDistGrid().getPoint(minindex)[0];
257 1020 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
258 1010 : double value = getFesGridPntr()->getValue(l);
259 1010 : double argument = targetDistGrid().getPoint(l)[0];
260 1010 : value = beta*value + (beta_prime_min-beta)*argument;
261 1010 : if(value<minval) {
262 : minval=value;
263 : minpos=argument;
264 : minindex=l;
265 : }
266 : }
267 : // Find minimum energy at low temperature
268 10 : double minimum_low = minpos;
269 11 : for(Grid::index_t l=minindex; l>1; l-=1) {
270 11 : double argument = targetDistGrid().getPoint(l)[0];
271 11 : double argument_next = targetDistGrid().getPoint(l-1)[0];
272 11 : double value = getFesGridPntr()->getValue(l);
273 11 : double value_next = getFesGridPntr()->getValue(l-1);
274 11 : value = beta*value + (beta_prime_min-beta)*argument - minval;
275 11 : value_next = beta*value_next + (beta_prime_min-beta)*argument_next - minval;
276 11 : if (value<threshold_ && value_next>threshold_) {
277 10 : minimum_low = argument_next;
278 10 : break;
279 : }
280 : }
281 : // Find maximum energy at low temperature
282 10 : double maximum_low = minpos;
283 12 : for(Grid::index_t l=minindex; l<(targetDistGrid().getSize()-1); l++) {
284 12 : double argument = targetDistGrid().getPoint(l)[0];
285 12 : double argument_next = targetDistGrid().getPoint(l+1)[0];
286 12 : double value = getFesGridPntr()->getValue(l);
287 12 : double value_next = getFesGridPntr()->getValue(l+1);
288 12 : value = beta*value + (beta_prime_min-beta)*argument - minval;
289 12 : value_next = beta*value_next + (beta_prime_min-beta)*argument_next - minval;
290 12 : if (value<threshold_ && value_next>threshold_) {
291 10 : maximum_low = argument_next;
292 10 : break;
293 : }
294 : }
295 : // Find minimum of F(U) at temperature max
296 : minval=DBL_MAX;
297 1020 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
298 1010 : double value = getFesGridPntr()->getValue(l);
299 1010 : double argument = targetDistGrid().getPoint(l)[0];
300 1010 : value = beta*value + (beta_prime_max-beta)*argument;
301 1010 : if(value<minval) {
302 : minval=value;
303 : minpos=argument;
304 : minindex=l;
305 : }
306 : }
307 : // Find minimum energy at high temperature
308 10 : double minimum_high = minpos;
309 13 : for(Grid::index_t l=minindex; l>1; l-=1) {
310 13 : double argument = targetDistGrid().getPoint(l)[0];
311 13 : double argument_next = targetDistGrid().getPoint(l-1)[0];
312 13 : double value = getFesGridPntr()->getValue(l);
313 13 : double value_next = getFesGridPntr()->getValue(l-1);
314 13 : value = beta*value + (beta_prime_max-beta)*argument - minval;
315 13 : value_next = beta*value_next + (beta_prime_max-beta)*argument_next - minval;
316 13 : if (value<threshold_ && value_next>threshold_) {
317 10 : minimum_high = argument_next;
318 10 : break;
319 : }
320 : }
321 : // Find maximum energy at high temperature
322 10 : double maximum_high = minpos;
323 11 : for(Grid::index_t l=minindex; l<(targetDistGrid().getSize()-1); l++) {
324 11 : double argument = targetDistGrid().getPoint(l)[0];
325 11 : double argument_next = targetDistGrid().getPoint(l+1)[0];
326 11 : double value = getFesGridPntr()->getValue(l);
327 11 : double value_next = getFesGridPntr()->getValue(l+1);
328 11 : value = beta*value + (beta_prime_max-beta)*argument - minval;
329 11 : value_next = beta*value_next + (beta_prime_max-beta)*argument_next - minval;
330 11 : if (value<threshold_ && value_next>threshold_) {
331 10 : maximum_high = argument_next;
332 10 : break;
333 : }
334 : }
335 10 : double minimum = std::min(minimum_low,minimum_high);
336 10 : double maximum = std::max(maximum_low,maximum_high);
337 : // Construct uniform TD in the interval between minimum and maximum
338 20 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
339 : double norm = 0.0;
340 1020 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
341 1010 : double argument = targetDistGrid().getPoint(l)[0];
342 : double value = 1.0;
343 : double tmp;
344 1010 : if(argument < minimum) {
345 217 : if (smoothening_) tmp = GaussianSwitchingFunc(argument,minimum,sigma_[0]);
346 0 : else tmp = exp(-1.0*epsilon_);
347 : }
348 793 : else if(argument > maximum) {
349 199 : if (smoothening_) tmp = GaussianSwitchingFunc(argument,maximum,sigma_[0]);
350 0 : else tmp = exp(-1.0*epsilon_);
351 : }
352 : else {
353 : tmp = 1.0;
354 : }
355 : value *= tmp;
356 1010 : norm += integration_weights[l]*value;
357 1010 : targetDistGrid().setValue(l,value);
358 : }
359 10 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
360 2 : } else if(targetDistGrid().getDimension()==2) {
361 : // 2D variant: Multicanonical with order parameter
362 : // In this variant we find for each temperature the relevant region of potential energy and order parameter.
363 : // The target distribution will be the union of the relevant regions at all temperatures in the temperature interval.
364 2 : double beta = getBeta();
365 2 : double beta_prime_min = 1./(getKBoltzmann()*min_temp_);
366 2 : double beta_prime_max = 1./(getKBoltzmann()*max_temp_);
367 2 : plumed_massert(getFesGridPntr()!=NULL,"the FES grid has to be linked to use TD_Multicanonical!");
368 : // Set all to zero
369 5204 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
370 5202 : double value = exp(-1.0*epsilon_);
371 5202 : targetDistGrid().setValue(l,value);
372 : }
373 : // Loop over temperatures
374 44 : for(unsigned i=0; i<steps_temp_; i++) {
375 42 : double beta_prime=beta_prime_min + (beta_prime_max-beta_prime_min)*i/(steps_temp_-1);
376 : // Find minimum for this temperature
377 : double minval=DBL_MAX;
378 109284 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
379 109242 : double energy = targetDistGrid().getPoint(l)[0];
380 109242 : double value = getFesGridPntr()->getValue(l);
381 109242 : value = beta*value + (beta_prime-beta)*energy;
382 109242 : if(value<minval) {
383 : minval=value;
384 : }
385 : }
386 : // Now check which energies and volumes are below X kt
387 109284 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
388 109242 : double energy = targetDistGrid().getPoint(l)[0];
389 109242 : double value = getFesGridPntr()->getValue(l);
390 109242 : value = beta*value + (beta_prime-beta)*energy - minval;
391 109242 : if (value<threshold_) {
392 : double value = 1.0;
393 7076 : targetDistGrid().setValue(l,value);
394 : }
395 : }
396 : }
397 2 : if (smoothening_) {
398 2 : std::vector<unsigned> nbin=targetDistGrid().getNbin();
399 2 : std::vector<double> dx=targetDistGrid().getDx();
400 : // Smoothening
401 104 : for(unsigned i=0; i<nbin[0]; i++) {
402 5304 : for(unsigned j=0; j<nbin[1]; j++) {
403 5202 : std::vector<unsigned> indices(2);
404 5202 : indices[0]=i;
405 5202 : indices[1]=j;
406 5202 : Grid::index_t index = targetDistGrid().getIndex(indices);
407 5202 : double energy = targetDistGrid().getPoint(index)[0];
408 5202 : double volume = targetDistGrid().getPoint(index)[1];
409 5202 : double value = targetDistGrid().getValue(index);
410 5202 : if (value>(1-1.e-5)) { // Apply only if this grid point was 1.
411 : // Apply gaussians around
412 773 : std::vector<int> minBin(2), maxBin(2), deltaBin(2); // These cannot be unsigned
413 : // Only consider contributions less than n*sigma bins apart from the actual distance
414 773 : deltaBin[0]=std::floor(6*sigma_[0]/dx[0]);;
415 773 : deltaBin[1]=std::floor(6*sigma_[1]/dx[1]);;
416 : // For energy
417 773 : minBin[0]=i - deltaBin[0];
418 773 : if (minBin[0] < 0) minBin[0]=0;
419 773 : if (minBin[0] > (nbin[0]-1)) minBin[0]=nbin[0]-1;
420 773 : maxBin[0]=i + deltaBin[0];
421 773 : if (maxBin[0] > (nbin[0]-1)) maxBin[0]=nbin[0]-1;
422 : // For volume
423 773 : minBin[1]=j - deltaBin[1];
424 773 : if (minBin[1] < 0) minBin[1]=0;
425 773 : if (minBin[1] > (nbin[1]-1)) minBin[1]=nbin[1]-1;
426 773 : maxBin[1]=j + deltaBin[1];
427 773 : if (maxBin[1] > (nbin[1]-1)) maxBin[1]=nbin[1]-1;
428 31273 : for(unsigned l=minBin[0]; l<maxBin[0]+1; l++) {
429 549973 : for(unsigned m=minBin[1]; m<maxBin[1]+1; m++) {
430 519473 : std::vector<unsigned> indices_prime(2);
431 519473 : indices_prime[0]=l;
432 519473 : indices_prime[1]=m;
433 519473 : Grid::index_t index_prime = targetDistGrid().getIndex(indices_prime);
434 519473 : double energy_prime = targetDistGrid().getPoint(index_prime)[0];
435 519473 : double volume_prime = targetDistGrid().getPoint(index_prime)[1];
436 519473 : double value_prime = targetDistGrid().getValue(index_prime);
437 : // Apply gaussian
438 1558419 : double gaussian_value = GaussianSwitchingFunc(energy_prime,energy,sigma_[0])*GaussianSwitchingFunc(volume_prime,volume,sigma_[1]);
439 519473 : if (value_prime<gaussian_value) {
440 19817 : targetDistGrid().setValue(index_prime,gaussian_value);
441 : }
442 : }
443 : }
444 : }
445 : }
446 : }
447 : }
448 : // Normalize
449 4 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
450 : double norm = 0.0;
451 5204 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
452 5202 : double value = targetDistGrid().getValue(l);
453 5202 : norm += integration_weights[l]*value;
454 : }
455 2 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
456 0 : } else plumed_merror(getName()+": Number of arguments for this target distribution must be 1 or 2");
457 : }
458 14 : updateLogTargetDistGrid();
459 14 : }
460 :
461 : inline
462 : double TD_Multicanonical::GaussianSwitchingFunc(const double argument, const double center, const double sigma) const {
463 1039362 : if(sigma>0.0) {
464 1039362 : double arg=(argument-center)/sigma;
465 1039362 : return exp(-0.5*arg*arg);
466 : }
467 : else {
468 : return 0.0;
469 : }
470 : }
471 :
472 :
473 : }
474 : }
|
