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