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