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