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