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