Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2018 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-arXiv-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 1.
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 rigurously.
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 :
78 : Much like in the Wang-Landau algorithm \cite wanglandau or in the multicanonical ensemble \cite Berg-PRL-1992 , a flat histogram is targeted.
79 : 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.
80 :
81 : 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:
82 : \f[
83 : \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
84 : \f]
85 : 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.
86 : \f$ \beta F_{\beta,P}(E,\mathcal{V},s) \f$ is not know from the start and is instead found during the simulation.
87 : Therefore \f$ p(E,\mathcal{V},s) \f$ is determined iteratively as done in the well tempered target distribution \cite Valsson-JCTC-2015.
88 :
89 : The output of these simulations can be reweighted in order to obtain information at all temperatures and pressures in the targeted region of TP plane.
90 : The reweighting can be performed using the action \ref REWEIGHT_TEMP_PRESS.
91 :
92 : The multicanonical ensemble (fixed volume) can be targeted using \ref TD_MULTICANONICAL.
93 :
94 : \par Examples
95 :
96 : The following input can be used to run a simulation in the multithermal-multibaric ensemble.
97 : The region of the temperature-pressure plane that will be explored is 260-350 K and 1 bar- 300 MPa.
98 : The energy and the volume are used as collective variables.
99 : Legendre polynomials are used to construct the two dimensional bias potential.
100 : The averaged stochastic gradient descent algorithm is chosen to optimize the VES functional.
101 : The target distribution is updated every 100 optimization steps (200 ps here) using the last estimation of the free energy.
102 :
103 : \plumedfile
104 : # Use energy and volume as CVs
105 : energy: ENERGY
106 : vol: VOLUME
107 :
108 : # Basis functions
109 : bf1: BF_LEGENDRE ORDER=10 MINIMUM=-14750 MAXIMUM=-12250
110 : bf2: BF_LEGENDRE ORDER=10 MINIMUM=6.5 MAXIMUM=8.25
111 :
112 : # Target distribution
113 : TD_MULTITHERMAL_MULTIBARIC ...
114 : MIN_TEMP=260
115 : MAX_TEMP=350
116 : MAX_PRESSURE=180.66422571 # 300 MPa
117 : MIN_PRESSURE=0.06022140857 # 1 bar
118 : PRESSURE=0.06022140857 # 1 bar
119 : STEPS_PRESSURE=20
120 : STEPS_TEMP=20
121 : SIGMA=50.,0.05
122 : THRESHOLD=1
123 : LABEL=td_multi
124 : ... TD_MULTITHERMAL_MULTIBARIC
125 :
126 : # Bias expansion
127 : VES_LINEAR_EXPANSION ...
128 : ARG=energy,vol
129 : BASIS_FUNCTIONS=bf1,bf2
130 : TEMP=300.0
131 : GRID_BINS=200,200
132 : TARGET_DISTRIBUTION=td_multi
133 : LABEL=b1
134 : ... VES_LINEAR_EXPANSION
135 :
136 : # Optimization algorithm
137 : OPT_AVERAGED_SGD ...
138 : BIAS=b1
139 : STRIDE=500
140 : LABEL=o1
141 : STEPSIZE=1.0
142 : FES_OUTPUT=500
143 : BIAS_OUTPUT=500
144 : TARGETDIST_OUTPUT=500
145 : COEFFS_OUTPUT=100
146 : TARGETDIST_STRIDE=100
147 : ... OPT_AVERAGED_SGD
148 :
149 : \endplumedfile
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 \ref FCCUBIC.
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 : \plumedfile
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 : \endplumedfile
212 :
213 : */
214 : //+ENDPLUMEDOC
215 :
216 : class TD_MultithermalMultibaric: public TargetDistribution {
217 : private:
218 : double threshold_, min_temp_, max_temp_;
219 : double min_press_, max_press_, press_;
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();
226 : double getValue(const std::vector<double>&) const;
227 6 : ~TD_MultithermalMultibaric() {}
228 : double GaussianSwitchingFunc(const double, const double, const double) const;
229 : };
230 :
231 :
232 7836 : PLUMED_REGISTER_ACTION(TD_MultithermalMultibaric,"TD_MULTITHERMAL_MULTIBARIC")
233 :
234 :
235 3 : void TD_MultithermalMultibaric::registerKeywords(Keywords& keys) {
236 3 : TargetDistribution::registerKeywords(keys);
237 15 : keys.add("compulsory","THRESHOLD","1","Maximum exploration free energy in kT.");
238 12 : keys.add("compulsory","MIN_TEMP","Minimum energy.");
239 12 : keys.add("compulsory","MAX_TEMP","Maximum energy.");
240 12 : keys.add("compulsory","MIN_PRESSURE","Minimum pressure.");
241 12 : keys.add("compulsory","MAX_PRESSURE","Maximum pressure.");
242 12 : keys.add("compulsory","PRESSURE","Target pressure of the barostat used in the MD engine.");
243 15 : keys.add("compulsory","STEPS_TEMP","20","Number of temperature steps.");
244 15 : keys.add("compulsory","STEPS_PRESSURE","20","Number of pressure steps.");
245 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 smooting is performed, this is the default behavior.");
246 3 : }
247 :
248 :
249 2 : TD_MultithermalMultibaric::TD_MultithermalMultibaric(const ActionOptions& ao):
250 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
251 : threshold_(1.0),
252 : min_temp_(0.0),
253 : max_temp_(1000.0),
254 : min_press_(0.0),
255 : max_press_(1000.0),
256 : sigma_(0.0),
257 : steps_temp_(20),
258 2 : steps_pressure_(20)
259 : {
260 2 : log.printf(" Multithermal-multibaric target distribution");
261 2 : log.printf("\n");
262 :
263 2 : log.printf(" Please read and cite ");
264 6 : log << plumed.cite("Piaggi and Parrinello, Phys. Rev. Lett. 122 (5), 050601 (2019)");
265 2 : log.printf(" and ");
266 6 : log << plumed.cite("Piaggi and Parrinello, arXiv preprint arXiv:1904.05624 (2019)");
267 2 : log.printf("\n");
268 :
269 :
270 4 : parse("THRESHOLD",threshold_);
271 2 : if(threshold_<=0.0) {
272 0 : plumed_merror("TD_MULTITHERMAL_MULTIBARIC target distribution: the value of the threshold should be positive.");
273 : }
274 4 : parse("MIN_TEMP",min_temp_);
275 4 : parse("MAX_TEMP",max_temp_);
276 4 : parse("MIN_PRESSURE",min_press_);
277 4 : parse("MAX_PRESSURE",max_press_);
278 4 : parse("PRESSURE",press_);
279 4 : parseVector("SIGMA",sigma_);
280 2 : if(sigma_.size()<2 || sigma_.size()>3) plumed_merror(getName()+": SIGMA takes 2 or 3 values as input.");
281 4 : parse("STEPS_TEMP",steps_temp_);
282 4 : parse("STEPS_PRESSURE",steps_pressure_);
283 2 : steps_temp_ += 1;
284 2 : steps_pressure_ += 1;
285 :
286 : setDynamic();
287 : setFesGridNeeded();
288 2 : checkRead();
289 2 : }
290 :
291 :
292 0 : double TD_MultithermalMultibaric::getValue(const std::vector<double>& argument) const {
293 0 : plumed_merror("getValue not implemented for TD_MultithermalMultibaric");
294 : return 0.0;
295 : }
296 :
297 :
298 6 : void TD_MultithermalMultibaric::updateGrid() {
299 6 : if (getStep() == 0) {
300 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");
301 2 : if(sigma_.size()!=targetDistGrid().getDimension()) plumed_merror(getName()+": mismatch between SIGMA dimension and number of arguments");
302 : // Use uniform TD
303 8 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
304 : double norm = 0.0;
305 23728 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
306 : double value = 1.0;
307 11862 : norm += integration_weights[l]*value;
308 11862 : targetDistGrid().setValue(l,value);
309 : }
310 4 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
311 2 : logTargetDistGrid().setMinToZero();
312 : } else {
313 4 : double beta = getBeta();
314 8 : double beta_prime_min = 1./(plumed.getAtoms().getKBoltzmann()*min_temp_);
315 8 : double beta_prime_max = 1./(plumed.getAtoms().getKBoltzmann()*max_temp_);
316 4 : plumed_massert(getFesGridPntr()!=NULL,"the FES grid has to be linked to use TD_MultithermalMultibaric!");
317 : // Set all to zero
318 47452 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
319 : double value = 0.0;
320 23724 : targetDistGrid().setValue(l,value);
321 : }
322 : // Loop over pressures and temperatures
323 84 : for(unsigned i=0; i<steps_temp_; i++) {
324 84 : double beta_prime=beta_prime_min + (beta_prime_max-beta_prime_min)*i/(steps_temp_-1);
325 1848 : for(unsigned j=0; j<steps_pressure_; j++) {
326 1764 : double pressure_prime=min_press_ + (max_press_-min_press_)*j/(steps_pressure_-1);
327 : // Find minimum for this pressure and temperature
328 : double minval=DBL_MAX;
329 20928096 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
330 31386852 : double energy = targetDistGrid().getPoint(l)[0];
331 31386852 : double volume = targetDistGrid().getPoint(l)[1];
332 10462284 : double value = getFesGridPntr()->getValue(l);
333 10462284 : value = beta*value + (beta_prime-beta)*energy + (beta_prime*pressure_prime-beta*press_)*volume;
334 10462284 : if(value<minval) {
335 : minval=value;
336 : }
337 : }
338 : // Now check which energies and volumes are below X kt
339 20926332 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
340 31386852 : double energy = targetDistGrid().getPoint(l)[0];
341 31386852 : double volume = targetDistGrid().getPoint(l)[1];
342 10462284 : double value = getFesGridPntr()->getValue(l);
343 10462284 : value = beta*value + (beta_prime-beta)*energy + (beta_prime*pressure_prime-beta*press_)*volume - minval;
344 10462284 : if (value<threshold_) {
345 : double value = 1.0;
346 113833 : targetDistGrid().setValue(l,value);
347 : }
348 : }
349 : }
350 : }
351 4 : std::vector<unsigned> nbin=targetDistGrid().getNbin();
352 4 : std::vector<double> dx=targetDistGrid().getDx();
353 4 : unsigned dim=targetDistGrid().getDimension();
354 : // Smoothening
355 47452 : for(Grid::index_t index=0; index<targetDistGrid().getSize(); index++) {
356 23724 : std::vector<unsigned> indices = targetDistGrid().getIndices(index);
357 23724 : std::vector<double> point = targetDistGrid().getPoint(index);
358 23724 : double value = targetDistGrid().getValue(index);
359 23724 : if (value>(1-1.e-5)) { // Apply only if this grid point was 1.
360 : // Apply gaussians around
361 4400 : std::vector<int> minBin(dim), maxBin(dim); // These cannot be unsigned
362 : // Only consider contributions less than n*sigma bins apart from the actual distance
363 15612 : for(unsigned k=0; k<dim; k++) {
364 33636 : int deltaBin=std::floor(5*sigma_[k]/dx[k]);
365 11212 : minBin[k]=indices[k] - deltaBin;
366 11212 : if (minBin[k] < 0) minBin[k]=0;
367 22424 : if (minBin[k] > (nbin[k]-1)) minBin[k]=nbin[k]-1;
368 11212 : maxBin[k]=indices[k] + deltaBin;
369 22424 : if (maxBin[k] > (nbin[k]-1)) maxBin[k]=nbin[k]-1;
370 : }
371 4400 : if (dim==2) {
372 38852 : for(unsigned l=minBin[0]; l<maxBin[0]+1; l++) {
373 546960 : for(unsigned m=minBin[1]; m<maxBin[1]+1; m++) {
374 256042 : std::vector<unsigned> indices_prime(dim);
375 256042 : indices_prime[0]=l;
376 256042 : indices_prime[1]=m;
377 256042 : Grid::index_t index_prime = targetDistGrid().getIndex(indices_prime);
378 256042 : std::vector<double> point_prime = targetDistGrid().getPoint(index_prime);
379 256042 : double value_prime = targetDistGrid().getValue(index_prime);
380 : // Apply gaussian
381 : double gaussian_value = 1;
382 512084 : for(unsigned k=0; k<dim; k++) {
383 2560420 : gaussian_value *= GaussianSwitchingFunc(point_prime[k],point[k],sigma_[k]);
384 : }
385 256042 : if (value_prime<gaussian_value) {
386 14938 : targetDistGrid().setValue(index_prime,gaussian_value);
387 : }
388 : }
389 : }
390 2412 : } else if (dim==3) {
391 76516 : for(unsigned l=minBin[0]; l<maxBin[0]+1; l++) {
392 402736 : for(unsigned m=minBin[1]; m<maxBin[1]+1; m++) {
393 2568380 : for(unsigned n=minBin[2]; n<maxBin[2]+1; n++) {
394 1118668 : std::vector<unsigned> indices_prime(dim);
395 1118668 : indices_prime[0]=l;
396 1118668 : indices_prime[1]=m;
397 1118668 : indices_prime[2]=n;
398 1118668 : Grid::index_t index_prime = targetDistGrid().getIndex(indices_prime);
399 1118668 : std::vector<double> point_prime = targetDistGrid().getPoint(index_prime);
400 1118668 : double value_prime = targetDistGrid().getValue(index_prime);
401 : // Apply gaussian
402 : double gaussian_value = 1;
403 3356004 : for(unsigned k=0; k<dim; k++) {
404 16780020 : gaussian_value *= GaussianSwitchingFunc(point_prime[k],point[k],sigma_[k]);
405 : }
406 1118668 : if (value_prime<gaussian_value) {
407 101881 : targetDistGrid().setValue(index_prime,gaussian_value);
408 : }
409 : }
410 : }
411 : }
412 : }
413 : }
414 : }
415 : // Normalize
416 16 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
417 : double norm = 0.0;
418 47456 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
419 23724 : double value = targetDistGrid().getValue(l);
420 23724 : norm += integration_weights[l]*value;
421 : }
422 8 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
423 4 : logTargetDistGrid().setMinToZero();
424 : }
425 6 : }
426 :
427 : inline
428 : double TD_MultithermalMultibaric::GaussianSwitchingFunc(const double argument, const double center, const double sigma) const {
429 3868088 : if(sigma>0.0) {
430 3868088 : double arg=(argument-center)/sigma;
431 3868088 : return exp(-0.5*arg*arg);
432 : }
433 : else {
434 : return 0.0;
435 : }
436 : }
437 :
438 :
439 : }
440 5874 : }
|
