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 4 : ~TD_MultithermalMultibaric() {}
228 : double GaussianSwitchingFunc(const double, const double, const double) const;
229 : };
230 :
231 :
232 13787 : PLUMED_REGISTER_ACTION(TD_MultithermalMultibaric,"TD_MULTITHERMAL_MULTIBARIC")
233 :
234 :
235 6 : void TD_MultithermalMultibaric::registerKeywords(Keywords& keys) {
236 6 : TargetDistribution::registerKeywords(keys);
237 12 : keys.add("compulsory","THRESHOLD","5","Maximum exploration free energy in kT.");
238 12 : keys.add("compulsory","EPSILON","10","The zeros of the target distribution are changed to e^-EPSILON.");
239 12 : keys.add("compulsory","MIN_TEMP","Minimum energy.");
240 12 : keys.add("compulsory","MAX_TEMP","Maximum energy.");
241 12 : keys.add("compulsory","MIN_PRESSURE","Minimum pressure.");
242 12 : keys.add("compulsory","MAX_PRESSURE","Maximum pressure.");
243 12 : keys.add("compulsory","PRESSURE","Target pressure of the barostat used in the MD engine.");
244 12 : keys.add("compulsory","STEPS_TEMP","20","Number of temperature steps.");
245 12 : keys.add("compulsory","STEPS_PRESSURE","20","Number of pressure steps.");
246 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.");
247 6 : }
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 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) {
286 0 : smoothening_=false;
287 : }
288 2 : if(smoothening_ && (sigma_.size()<2 || sigma_.size()>3) ) {
289 0 : plumed_merror(getName()+": SIGMA takes 2 or 3 values as input.");
290 : }
291 2 : if (smoothening_) {
292 2 : log.printf(" the target distribution will be smoothed using sigma values");
293 7 : for(unsigned i=0; i<sigma_.size(); ++i) {
294 5 : log.printf(" %f",sigma_[i]);
295 : }
296 2 : log.printf("\n");
297 : }
298 2 : parse("STEPS_TEMP",steps_temp_);
299 2 : parse("STEPS_PRESSURE",steps_pressure_);
300 2 : log.printf(" %d steps in temperatures and %d steps in pressure will be employed \n",steps_temp_,steps_pressure_);
301 2 : steps_temp_ += 1;
302 2 : steps_pressure_ += 1;
303 2 : parse("EPSILON",epsilon_);
304 2 : if(epsilon_<=1.0) {
305 0 : plumed_merror(getName()+": the value of epsilon should be greater than 1.");
306 : }
307 2 : log.printf(" the non relevant regions of the target distribution are set to e^-%f \n",epsilon_);
308 : setDynamic();
309 : setFesGridNeeded();
310 2 : checkRead();
311 2 : }
312 :
313 :
314 0 : double TD_MultithermalMultibaric::getValue(const std::vector<double>& argument) const {
315 0 : plumed_merror("getValue not implemented for TD_MultithermalMultibaric");
316 : return 0.0;
317 : }
318 :
319 :
320 4 : void TD_MultithermalMultibaric::updateGrid() {
321 4 : if (getStep() == 0) {
322 2 : if(targetDistGrid().getDimension()>3 || targetDistGrid().getDimension()<2) {
323 0 : plumed_merror(getName()+" works only with 2 or 3 arguments, i.e. energy and volume, or energy, volume, and CV");
324 : }
325 2 : if(smoothening_ && sigma_.size()!=targetDistGrid().getDimension()) {
326 0 : plumed_merror(getName()+": mismatch between SIGMA dimension and number of arguments");
327 : }
328 : // Use uniform TD
329 4 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
330 : double norm = 0.0;
331 11864 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
332 : double value = 1.0;
333 11862 : norm += integration_weights[l]*value;
334 11862 : targetDistGrid().setValue(l,value);
335 : }
336 2 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
337 : } else {
338 2 : double beta = getBeta();
339 2 : double beta_prime_min = 1./(plumed.getAtoms().getKBoltzmann()*min_temp_);
340 2 : double beta_prime_max = 1./(plumed.getAtoms().getKBoltzmann()*max_temp_);
341 2 : plumed_massert(getFesGridPntr()!=NULL,"the FES grid has to be linked to use TD_MultithermalMultibaric!");
342 : // Set all to current epsilon value
343 11864 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
344 11862 : double value = exp(-1.0*epsilon_);
345 11862 : targetDistGrid().setValue(l,value);
346 : }
347 : // Loop over pressures and temperatures
348 44 : for(unsigned i=0; i<steps_temp_; i++) {
349 42 : double beta_prime=beta_prime_min + (beta_prime_max-beta_prime_min)*i/(steps_temp_-1);
350 924 : for(unsigned j=0; j<steps_pressure_; j++) {
351 882 : double pressure_prime=min_press_ + (max_press_-min_press_)*j/(steps_pressure_-1);
352 : // Find minimum for this pressure and temperature
353 : double minval=DBL_MAX;
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;
359 5231142 : if(value<minval) {
360 : minval=value;
361 : }
362 : }
363 : // Now check which energies and volumes are below X kt
364 5232024 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
365 5231142 : double energy = targetDistGrid().getPoint(l)[0];
366 5231142 : double volume = targetDistGrid().getPoint(l)[1];
367 5231142 : double value = getFesGridPntr()->getValue(l);
368 5231142 : value = beta*value + (beta_prime-beta)*energy + (beta_prime*pressure_prime-beta*press_)*volume - minval;
369 5231142 : if (value<threshold_) {
370 : double value = 1.0;
371 65261 : targetDistGrid().setValue(l,value);
372 : }
373 : }
374 : }
375 : }
376 2 : if (smoothening_) {
377 2 : std::vector<unsigned> nbin=targetDistGrid().getNbin();
378 2 : std::vector<double> dx=targetDistGrid().getDx();
379 2 : unsigned dim=targetDistGrid().getDimension();
380 : // Smoothening
381 11864 : for(Grid::index_t index=0; index<targetDistGrid().getSize(); index++) {
382 11862 : std::vector<unsigned> indices = targetDistGrid().getIndices(index);
383 11862 : std::vector<double> point = targetDistGrid().getPoint(index);
384 11862 : double value = targetDistGrid().getValue(index);
385 11862 : if (value>(1-1.e-5)) { // Apply only if this grid point was 1.
386 : // Apply gaussians around
387 1242 : std::vector<int> minBin(dim), maxBin(dim); // These cannot be unsigned
388 : // Only consider contributions less than n*sigma bins apart from the actual distance
389 4384 : for(unsigned k=0; k<dim; k++) {
390 3142 : int deltaBin=std::floor(6*sigma_[k]/dx[k]);
391 3142 : minBin[k]=indices[k] - deltaBin;
392 3142 : if (minBin[k] < 0) {
393 947 : minBin[k]=0;
394 : }
395 3142 : if (minBin[k] > (nbin[k]-1)) {
396 0 : minBin[k]=nbin[k]-1;
397 : }
398 3142 : maxBin[k]=indices[k] + deltaBin;
399 3142 : if (maxBin[k] > (nbin[k]-1)) {
400 541 : maxBin[k]=nbin[k]-1;
401 : }
402 : }
403 1242 : if (dim==2) {
404 7008 : for(unsigned l=minBin[0]; l<maxBin[0]+1; l++) {
405 115632 : for(unsigned m=minBin[1]; m<maxBin[1]+1; m++) {
406 109208 : std::vector<unsigned> indices_prime(dim);
407 109208 : indices_prime[0]=l;
408 109208 : indices_prime[1]=m;
409 109208 : Grid::index_t index_prime = targetDistGrid().getIndex(indices_prime);
410 109208 : std::vector<double> point_prime = targetDistGrid().getPoint(index_prime);
411 109208 : double value_prime = targetDistGrid().getValue(index_prime);
412 : // Apply gaussian
413 : double gaussian_value = 1;
414 327624 : for(unsigned k=0; k<dim; k++) {
415 436832 : gaussian_value *= GaussianSwitchingFunc(point_prime[k],point[k],sigma_[k]);
416 : }
417 109208 : if (value_prime<gaussian_value) {
418 2761 : targetDistGrid().setValue(index_prime,gaussian_value);
419 : }
420 : }
421 : }
422 658 : } else if (dim==3) {
423 12352 : for(unsigned l=minBin[0]; l<maxBin[0]+1; l++) {
424 84952 : for(unsigned m=minBin[1]; m<maxBin[1]+1; m++) {
425 509608 : for(unsigned n=minBin[2]; n<maxBin[2]+1; n++) {
426 436350 : std::vector<unsigned> indices_prime(dim);
427 436350 : indices_prime[0]=l;
428 436350 : indices_prime[1]=m;
429 436350 : indices_prime[2]=n;
430 436350 : Grid::index_t index_prime = targetDistGrid().getIndex(indices_prime);
431 436350 : std::vector<double> point_prime = targetDistGrid().getPoint(index_prime);
432 436350 : double value_prime = targetDistGrid().getValue(index_prime);
433 : // Apply gaussian
434 : double gaussian_value = 1;
435 1745400 : for(unsigned k=0; k<dim; k++) {
436 2618100 : gaussian_value *= GaussianSwitchingFunc(point_prime[k],point[k],sigma_[k]);
437 : }
438 436350 : if (value_prime<gaussian_value) {
439 13789 : targetDistGrid().setValue(index_prime,gaussian_value);
440 : }
441 : }
442 : }
443 : }
444 : }
445 : }
446 : }
447 : }
448 : // Normalize
449 4 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
450 : double norm = 0.0;
451 11864 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
452 11862 : double value = targetDistGrid().getValue(l);
453 11862 : norm += integration_weights[l]*value;
454 : }
455 2 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
456 : }
457 4 : updateLogTargetDistGrid();
458 4 : }
459 :
460 : inline
461 : double TD_MultithermalMultibaric::GaussianSwitchingFunc(const double argument, const double center, const double sigma) const {
462 1527466 : if(sigma>0.0) {
463 1527466 : double arg=(argument-center)/sigma;
464 1527466 : return exp(-0.5*arg*arg);
465 : } else {
466 : return 0.0;
467 : }
468 : }
469 :
470 :
471 : }
472 : }
|
