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_MULTITHERMAL_MULTIBARIC
36 : /*
37 : Multithermal-multibaric target distribution (dynamic).
38 :
39 : Use the target distribution to sample the multithermal-multibaric ensemble \cite Piaggi-PRL-2019 \cite Okumura-CPL-2004.
40 : In this way, in a single molecular dynamics simulation one can obtain information
41 : about the simulated system in a range of temperatures and pressures.
42 : This range is determined through the keywords MIN_TEMP, MAX_TEMP, MIN_PRESSURE, and MAX_PRESSURE.
43 : One should also specified the target pressure of the barostat with the keyword PRESSURE.
44 :
45 : The collective variables (CVs) used to construct the bias potential must be:
46 : 1. the potential energy and the volume or,
47 : 2. the potential energy, the volume, and an order parameter.
48 :
49 : Other choices of CVs or a different order of the above mentioned CVs are nonsensical.
50 : The third CV, the order parameter, must be used when the region of the phase diagram under study is crossed by a first order phase transition \cite Piaggi-JCP-2019 .
51 :
52 : The algorithm will explore the free energy at each temperature and pressure up to a predefined free
53 : energy threshold \f$\epsilon\f$ specified through the keyword THRESHOLD (in kT units).
54 : If only the energy and the volume are being biased, i.e. no phase transition is considered, then THRESHOLD can be set to around 5.
55 : If also an order parameter is used then the THRESHOLD should be greater than the barrier for the transformation in kT.
56 : For small systems undergoing a freezing transition THRESHOLD is typically between 20 and 50.
57 :
58 : It is also important to specify the number of intermediate temperatures and pressures to consider.
59 : This is done through the keywords STEPS_TEMP and STEPS_PRESSURE.
60 : If the number of intermediate temperature and pressures is too small, then holes might appear in the target distribution.
61 : If it is too large, the performance will deteriorate with no additional advantage.
62 :
63 : We now describe the algorithm more rigorously.
64 : The target distribution is given by
65 : \f[
66 : p(E,\mathcal{V},s)=
67 : \begin{cases}
68 : 1/\Omega_{E,\mathcal{V},s} & \text{if there is at least one } \beta',P' \text{ such} \\
69 : & \text{that } \beta' F_{\beta',P'}(E,\mathcal{V},s)<\epsilon \text{ with} \\
70 : & \beta_1>\beta'>\beta_2 \text{ and } P_1<P'<P_2 \\
71 : 0 & \text{otherwise}
72 : \end{cases}
73 : \f]
74 : with \f$F_{\beta',P'}(E,\mathcal{V},s)\f$ the free energy as a function of energy \f$E\f$ and volume \f$\mathcal{V}\f$ (and optionally the order parameter \f$s\f$) at temperature \f$\beta'\f$ and pressure \f$P'\f$, \f$\Omega_{E,\mathcal{V},s}\f$ is a normalization constant, and \f$\epsilon\f$ is the THRESHOLD.
75 : In practice the condition \f$\beta' F_{\beta',P'}(E,\mathcal{V},s)<\epsilon\f$ is checked in equally spaced points in each dimension \f$\beta'\f$ and \f$P'\f$.
76 : The number of points is determined with the keywords STEPS_TEMP and STEPS_PRESSURE.
77 : In practice the target distribution is never set to zero but rather to a small value controlled by the keyword EPSILON.
78 : The small value is e^-EPSILON.
79 :
80 : Much like in the Wang-Landau algorithm \cite wanglandau or in the multicanonical ensemble \cite Berg-PRL-1992 , a flat histogram is targeted.
81 : The idea behind this choice of target distribution is that all regions of potential energy and volume (and optionally order parameter) that are relevant at all temperatures \f$\beta_1<\beta'<\beta_2\f$ and pressure \f$P_1<P'<P_2\f$ are included in the distribution.
82 :
83 : The free energy at temperature \f$\beta'\f$ and pressure \f$P'\f$ is calculated from the free energy at \f$\beta\f$ and \f$P\f$ using:
84 : \f[
85 : \beta' F_{\beta',P'}(E,\mathcal{V},s) = \beta F_{\beta,P}(E,\mathcal{V},s) + (\beta' - \beta) E + (\beta' P' - \beta P ) \mathcal{V} + C
86 : \f]
87 : with \f$C\f$ such that \f$F_{\beta',P'}(E_m,\mathcal{V}_m,s_m)=0\f$ with \f$E_{m},\mathcal{V}_m,s_m\f$ the position of the free energy minimum.
88 : \f$ \beta F_{\beta,P}(E,\mathcal{V},s) \f$ is not know from the start and is instead found during the simulation.
89 : Therefore \f$ p(E,\mathcal{V},s) \f$ is determined iteratively as done in the well tempered target distribution \cite Valsson-JCTC-2015.
90 :
91 : The output of these simulations can be reweighted in order to obtain information at all temperatures and pressures in the targeted region of Temperature-Pressure plane.
92 : The reweighting can be performed using the action \ref REWEIGHT_TEMP_PRESS.
93 :
94 : The multicanonical ensemble (fixed volume) can be targeted using \ref TD_MULTICANONICAL.
95 :
96 : \par Examples
97 :
98 : The following input can be used to run a simulation in the multithermal-multibaric ensemble.
99 : The region of the temperature-pressure plane that will be explored is 260-350 K and 1 bar- 300 MPa.
100 : The energy and the volume are used as collective variables.
101 : Legendre polynomials are used to construct the two dimensional bias potential.
102 : The averaged stochastic gradient descent algorithm is chosen to optimize the VES functional.
103 : The target distribution is updated every 100 optimization steps (200 ps here) using the last estimation of the free energy.
104 :
105 : \plumedfile
106 : # Use energy and volume as CVs
107 : energy: ENERGY
108 : vol: VOLUME
109 :
110 : # Basis functions
111 : bf1: BF_LEGENDRE ORDER=10 MINIMUM=-14750 MAXIMUM=-12250
112 : bf2: BF_LEGENDRE ORDER=10 MINIMUM=6.5 MAXIMUM=8.25
113 :
114 : # Target distribution - 1 bar = 0.06022140857 and 300 MPa = 180.66422571
115 : TD_MULTITHERMAL_MULTIBARIC ...
116 : MIN_TEMP=260
117 : MAX_TEMP=350
118 : MAX_PRESSURE=180.66422571
119 : MIN_PRESSURE=0.06022140857
120 : PRESSURE=0.06022140857
121 : LABEL=td_multi
122 : ... TD_MULTITHERMAL_MULTIBARIC
123 :
124 : # Bias expansion
125 : VES_LINEAR_EXPANSION ...
126 : ARG=energy,vol
127 : BASIS_FUNCTIONS=bf1,bf2
128 : TEMP=300.0
129 : GRID_BINS=200,200
130 : TARGET_DISTRIBUTION=td_multi
131 : LABEL=b1
132 : ... VES_LINEAR_EXPANSION
133 :
134 : # Optimization algorithm
135 : OPT_AVERAGED_SGD ...
136 : BIAS=b1
137 : STRIDE=500
138 : LABEL=o1
139 : STEPSIZE=1.0
140 : FES_OUTPUT=500
141 : BIAS_OUTPUT=500
142 : TARGETDIST_OUTPUT=500
143 : COEFFS_OUTPUT=100
144 : TARGETDIST_STRIDE=100
145 : ... OPT_AVERAGED_SGD
146 :
147 : \endplumedfile
148 :
149 :
150 : The multithermal-multibaric target distribution can also be used to explore regions of the phase diagram crossed by first order phase transitions.
151 : Consider a system of 250 atoms that crystallizes in the FCC crystal structure.
152 : The region of the temperature-pressure plane that will be explored is 350-450 K and 1bar-1GPa.
153 : We assume that inside this region we can find the liquid-FCC coexistence line that we would like to obtain.
154 : In this case in addition to the energy and volume, an order parameter must also be biased.
155 : The energy, volume, and an order parameter are used as collective variables to construct the bias potential.
156 : We choose as order parameter the \ref FCCUBIC.
157 : Legendre polynomials are used to construct the three dimensional bias potential.
158 : The averaged stochastic gradient descent algorithm is chosen to optimize the VES functional.
159 : The target distribution is updated every 100 optimization steps (200 ps here) using the last estimation of the free energy.
160 :
161 : \plumedfile
162 : # Use energy, volume and FCCUBIC as CVs
163 : energy: ENERGY
164 : vol: VOLUME
165 : fcc: FCCUBIC SPECIES=1-256 SWITCH={CUBIC D_0=0.4 D_MAX=0.5} MORE_THAN={RATIONAL R_0=0.45 NN=12 MM=24}
166 :
167 : # Basis functions
168 : bf1: BF_LEGENDRE ORDER=8 MINIMUM=-26500 MAXIMUM=-23500
169 : bf2: BF_LEGENDRE ORDER=8 MINIMUM=8.0 MAXIMUM=11.5
170 : bf3: BF_LEGENDRE ORDER=8 MINIMUM=0.0 MAXIMUM=256.0
171 :
172 : # Target distribution
173 : TD_MULTITHERMAL_MULTIBARIC ...
174 : LABEL=td_multitp
175 : MIN_TEMP=350.0
176 : MAX_TEMP=450.0
177 : MIN_PRESSURE=0.06022140857
178 : MAX_PRESSURE=602.2140857
179 : PRESSURE=301.10704285
180 : SIGMA=250.0,0.1,10.0
181 : THRESHOLD=15
182 : STEPS_TEMP=20
183 : STEPS_PRESSURE=20
184 : ... TD_MULTITHERMAL_MULTIBARIC
185 :
186 : # Expansion
187 : VES_LINEAR_EXPANSION ...
188 : ARG=energy,vol,fcc.morethan
189 : BASIS_FUNCTIONS=bf1,bf2,bf3
190 : TEMP=400.0
191 : GRID_BINS=40,40,40
192 : TARGET_DISTRIBUTION=td_multitp
193 : LABEL=b1
194 : ... VES_LINEAR_EXPANSION
195 :
196 : # Optimization algorithm
197 : OPT_AVERAGED_SGD ...
198 : BIAS=b1
199 : STRIDE=500
200 : LABEL=o1
201 : STEPSIZE=1.0
202 : FES_OUTPUT=500
203 : BIAS_OUTPUT=500
204 : TARGETDIST_OUTPUT=500
205 : COEFFS_OUTPUT=100
206 : TARGETDIST_STRIDE=500
207 : ... OPT_AVERAGED_SGD
208 :
209 : \endplumedfile
210 :
211 : */
212 : //+ENDPLUMEDOC
213 :
214 : class TD_MultithermalMultibaric: public TargetDistribution {
215 : private:
216 : double threshold_, min_temp_, max_temp_;
217 : double min_press_, max_press_, press_;
218 : double epsilon_;
219 : bool smoothening_;
220 : std::vector<double> sigma_;
221 : unsigned steps_temp_, steps_pressure_;
222 : public:
223 : static void registerKeywords(Keywords&);
224 : explicit TD_MultithermalMultibaric(const ActionOptions& ao);
225 : void updateGrid() override;
226 : double getValue(const std::vector<double>&) const override;
227 6 : ~TD_MultithermalMultibaric() {}
228 : double GaussianSwitchingFunc(const double, const double, const double) const;
229 : };
230 :
231 :
232 10421 : PLUMED_REGISTER_ACTION(TD_MultithermalMultibaric,"TD_MULTITHERMAL_MULTIBARIC")
233 :
234 :
235 3 : void TD_MultithermalMultibaric::registerKeywords(Keywords& keys) {
236 3 : TargetDistribution::registerKeywords(keys);
237 6 : keys.add("compulsory","THRESHOLD","5","Maximum exploration free energy in kT.");
238 6 : keys.add("compulsory","EPSILON","10","The zeros of the target distribution are changed to e^-EPSILON.");
239 6 : keys.add("compulsory","MIN_TEMP","Minimum energy.");
240 6 : keys.add("compulsory","MAX_TEMP","Maximum energy.");
241 6 : keys.add("compulsory","MIN_PRESSURE","Minimum pressure.");
242 6 : keys.add("compulsory","MAX_PRESSURE","Maximum pressure.");
243 6 : keys.add("compulsory","PRESSURE","Target pressure of the barostat used in the MD engine.");
244 6 : keys.add("compulsory","STEPS_TEMP","20","Number of temperature steps.");
245 6 : keys.add("compulsory","STEPS_PRESSURE","20","Number of pressure steps.");
246 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.");
247 3 : }
248 :
249 :
250 2 : TD_MultithermalMultibaric::TD_MultithermalMultibaric(const ActionOptions& ao):
251 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
252 2 : threshold_(5.0),
253 2 : min_temp_(0.0),
254 2 : max_temp_(1000.0),
255 2 : min_press_(0.0),
256 2 : max_press_(1000.0),
257 2 : epsilon_(10.0),
258 2 : smoothening_(true),
259 2 : steps_temp_(20),
260 2 : steps_pressure_(20)
261 : {
262 2 : log.printf(" Multithermal-multibaric target distribution");
263 2 : log.printf("\n");
264 :
265 2 : log.printf(" Please read and cite ");
266 4 : log << plumed.cite("Piaggi and Parrinello, Phys. Rev. Lett. 122 (5), 050601 (2019)");
267 2 : log.printf(" and ");
268 4 : log << plumed.cite("Piaggi and Parrinello, J. Chem. Phys. 150 (24), 244119 (2019)");
269 2 : log.printf("\n");
270 :
271 :
272 2 : parse("THRESHOLD",threshold_);
273 2 : if(threshold_<=0.0) {
274 0 : plumed_merror(getName()+": the value of the threshold should be positive.");
275 : }
276 2 : log.printf(" exploring free energy up to %f kT for each temperature and pressure\n",threshold_);
277 2 : parse("MIN_TEMP",min_temp_);
278 2 : parse("MAX_TEMP",max_temp_);
279 2 : log.printf(" temperatures between %f and %f will be explored \n",min_temp_,max_temp_);
280 2 : parse("MIN_PRESSURE",min_press_);
281 2 : parse("MAX_PRESSURE",max_press_);
282 2 : log.printf(" pressures between %f and %f will be explored \n",min_press_,max_press_);
283 2 : parse("PRESSURE",press_);
284 2 : log.printf(" pressure of the barostat should have been set to %f. Please check this in the MD engine \n",press_);
285 4 : parseVector("SIGMA",sigma_);
286 2 : if(sigma_.size()==0) smoothening_=false;
287 2 : if(smoothening_ && (sigma_.size()<2 || sigma_.size()>3) ) plumed_merror(getName()+": SIGMA takes 2 or 3 values as input.");
288 2 : if (smoothening_) {
289 2 : log.printf(" the target distribution will be smoothed using sigma values");
290 7 : for(unsigned i=0; i<sigma_.size(); ++i) log.printf(" %f",sigma_[i]);
291 2 : log.printf("\n");
292 : }
293 2 : parse("STEPS_TEMP",steps_temp_);
294 2 : parse("STEPS_PRESSURE",steps_pressure_);
295 2 : log.printf(" %d steps in temperatures and %d steps in pressure will be employed \n",steps_temp_,steps_pressure_);
296 2 : steps_temp_ += 1;
297 2 : steps_pressure_ += 1;
298 2 : parse("EPSILON",epsilon_);
299 2 : if(epsilon_<=1.0) {
300 0 : plumed_merror(getName()+": the value of epsilon should be greater than 1.");
301 : }
302 2 : log.printf(" the non relevant regions of the target distribution are set to e^-%f \n",epsilon_);
303 : setDynamic();
304 : setFesGridNeeded();
305 2 : checkRead();
306 2 : }
307 :
308 :
309 0 : double TD_MultithermalMultibaric::getValue(const std::vector<double>& argument) const {
310 0 : plumed_merror("getValue not implemented for TD_MultithermalMultibaric");
311 : return 0.0;
312 : }
313 :
314 :
315 4 : void TD_MultithermalMultibaric::updateGrid() {
316 4 : if (getStep() == 0) {
317 2 : if(targetDistGrid().getDimension()>3 || targetDistGrid().getDimension()<2) plumed_merror(getName()+" works only with 2 or 3 arguments, i.e. energy and volume, or energy, volume, and CV");
318 2 : if(smoothening_ && sigma_.size()!=targetDistGrid().getDimension()) plumed_merror(getName()+": mismatch between SIGMA dimension and number of arguments");
319 : // Use uniform TD
320 4 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
321 : double norm = 0.0;
322 11864 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
323 : double value = 1.0;
324 11862 : norm += integration_weights[l]*value;
325 11862 : targetDistGrid().setValue(l,value);
326 : }
327 2 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
328 : } else {
329 2 : double beta = getBeta();
330 2 : double beta_prime_min = 1./(plumed.getAtoms().getKBoltzmann()*min_temp_);
331 2 : double beta_prime_max = 1./(plumed.getAtoms().getKBoltzmann()*max_temp_);
332 2 : plumed_massert(getFesGridPntr()!=NULL,"the FES grid has to be linked to use TD_MultithermalMultibaric!");
333 : // Set all to current epsilon value
334 11864 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
335 11862 : double value = exp(-1.0*epsilon_);
336 11862 : targetDistGrid().setValue(l,value);
337 : }
338 : // Loop over pressures and temperatures
339 44 : for(unsigned i=0; i<steps_temp_; i++) {
340 42 : double beta_prime=beta_prime_min + (beta_prime_max-beta_prime_min)*i/(steps_temp_-1);
341 924 : for(unsigned j=0; j<steps_pressure_; j++) {
342 882 : double pressure_prime=min_press_ + (max_press_-min_press_)*j/(steps_pressure_-1);
343 : // Find minimum for this pressure and temperature
344 : double minval=DBL_MAX;
345 5232024 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
346 10462284 : double energy = targetDistGrid().getPoint(l)[0];
347 10462284 : double volume = targetDistGrid().getPoint(l)[1];
348 5231142 : double value = getFesGridPntr()->getValue(l);
349 5231142 : value = beta*value + (beta_prime-beta)*energy + (beta_prime*pressure_prime-beta*press_)*volume;
350 5231142 : if(value<minval) {
351 : minval=value;
352 : }
353 : }
354 : // Now check which energies and volumes are below X kt
355 5232024 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
356 10462284 : double energy = targetDistGrid().getPoint(l)[0];
357 10462284 : double volume = targetDistGrid().getPoint(l)[1];
358 5231142 : double value = getFesGridPntr()->getValue(l);
359 5231142 : value = beta*value + (beta_prime-beta)*energy + (beta_prime*pressure_prime-beta*press_)*volume - minval;
360 5231142 : if (value<threshold_) {
361 : double value = 1.0;
362 65261 : targetDistGrid().setValue(l,value);
363 : }
364 : }
365 : }
366 : }
367 2 : if (smoothening_) {
368 2 : std::vector<unsigned> nbin=targetDistGrid().getNbin();
369 2 : std::vector<double> dx=targetDistGrid().getDx();
370 2 : unsigned dim=targetDistGrid().getDimension();
371 : // Smoothening
372 11864 : for(Grid::index_t index=0; index<targetDistGrid().getSize(); index++) {
373 11862 : std::vector<unsigned> indices = targetDistGrid().getIndices(index);
374 11862 : std::vector<double> point = targetDistGrid().getPoint(index);
375 11862 : double value = targetDistGrid().getValue(index);
376 11862 : if (value>(1-1.e-5)) { // Apply only if this grid point was 1.
377 : // Apply gaussians around
378 1242 : std::vector<int> minBin(dim), maxBin(dim); // These cannot be unsigned
379 : // Only consider contributions less than n*sigma bins apart from the actual distance
380 4384 : for(unsigned k=0; k<dim; k++) {
381 3142 : int deltaBin=std::floor(6*sigma_[k]/dx[k]);
382 3142 : minBin[k]=indices[k] - deltaBin;
383 3142 : if (minBin[k] < 0) minBin[k]=0;
384 3142 : if (minBin[k] > (nbin[k]-1)) minBin[k]=nbin[k]-1;
385 3142 : maxBin[k]=indices[k] + deltaBin;
386 3142 : if (maxBin[k] > (nbin[k]-1)) maxBin[k]=nbin[k]-1;
387 : }
388 1242 : if (dim==2) {
389 7008 : for(unsigned l=minBin[0]; l<maxBin[0]+1; l++) {
390 115632 : for(unsigned m=minBin[1]; m<maxBin[1]+1; m++) {
391 109208 : std::vector<unsigned> indices_prime(dim);
392 109208 : indices_prime[0]=l;
393 109208 : indices_prime[1]=m;
394 109208 : Grid::index_t index_prime = targetDistGrid().getIndex(indices_prime);
395 109208 : std::vector<double> point_prime = targetDistGrid().getPoint(index_prime);
396 109208 : double value_prime = targetDistGrid().getValue(index_prime);
397 : // Apply gaussian
398 : double gaussian_value = 1;
399 327624 : for(unsigned k=0; k<dim; k++) {
400 436832 : gaussian_value *= GaussianSwitchingFunc(point_prime[k],point[k],sigma_[k]);
401 : }
402 109208 : if (value_prime<gaussian_value) {
403 2761 : targetDistGrid().setValue(index_prime,gaussian_value);
404 : }
405 : }
406 : }
407 658 : } else if (dim==3) {
408 12352 : for(unsigned l=minBin[0]; l<maxBin[0]+1; l++) {
409 84952 : for(unsigned m=minBin[1]; m<maxBin[1]+1; m++) {
410 509608 : for(unsigned n=minBin[2]; n<maxBin[2]+1; n++) {
411 436350 : std::vector<unsigned> indices_prime(dim);
412 436350 : indices_prime[0]=l;
413 436350 : indices_prime[1]=m;
414 436350 : indices_prime[2]=n;
415 436350 : Grid::index_t index_prime = targetDistGrid().getIndex(indices_prime);
416 436350 : std::vector<double> point_prime = targetDistGrid().getPoint(index_prime);
417 436350 : double value_prime = targetDistGrid().getValue(index_prime);
418 : // Apply gaussian
419 : double gaussian_value = 1;
420 1745400 : for(unsigned k=0; k<dim; k++) {
421 2618100 : gaussian_value *= GaussianSwitchingFunc(point_prime[k],point[k],sigma_[k]);
422 : }
423 436350 : if (value_prime<gaussian_value) {
424 13789 : targetDistGrid().setValue(index_prime,gaussian_value);
425 : }
426 : }
427 : }
428 : }
429 : }
430 : }
431 : }
432 : }
433 : // Normalize
434 4 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
435 : double norm = 0.0;
436 11864 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
437 11862 : double value = targetDistGrid().getValue(l);
438 11862 : norm += integration_weights[l]*value;
439 : }
440 2 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
441 : }
442 4 : updateLogTargetDistGrid();
443 4 : }
444 :
445 : inline
446 : double TD_MultithermalMultibaric::GaussianSwitchingFunc(const double argument, const double center, const double sigma) const {
447 1527466 : if(sigma>0.0) {
448 1527466 : double arg=(argument-center)/sigma;
449 1527466 : return exp(-0.5*arg*arg);
450 : }
451 : else {
452 : return 0.0;
453 : }
454 : }
455 :
456 :
457 : }
458 : }
|
