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_MULTICANONICAL
36 : /*
37 : Multicanonical target distribution (dynamic).
38 :
39 : Use the target distribution to sample the multicanonical ensemble \cite Berg-PRL-1992 \cite Piaggi-PRL-2019.
40 : In this way, in a single molecular dynamics simulation one can obtain information about the system in a range of temperatures.
41 : This range is determined through the keywords MIN_TEMP and MAX_TEMP.
42 :
43 : The collective variables (CVs) used to construct the bias potential must be:
44 : 1. the energy or,
45 : 2. the energy and an order parameter.
46 :
47 : Other choices of CVs or a different order of the above mentioned CVs are nonsensical.
48 : The second CV, the order parameter, must be used when one aims at studying a first order phase transition in the chosen temperature interval \cite Piaggi-arXiv-2019.
49 :
50 : The algorithm will explore the free energy at each temperature up to a predefined free
51 : energy threshold \f$\epsilon\f$ specified through the keyword THRESHOLD (in kT units).
52 : If only the energy is biased, i.e. no phase transition is considered, then TRESHOLD can be set to 1.
53 : If also an order parameter is used then the THRESHOLD should be greater than the barrier for the transformation in kT.
54 : For small systems undergoing a freezing transition THRESHOLD is typically between 20 and 50.
55 :
56 : When only the potential energy is used as CV the method is equivalent to the Wang-Landau algorithm \cite wanglandau.
57 : The advantage with respect to Wang-Landau is that instead of sampling the potential energy indiscriminately, an interval is chosen on the fly based on the minimum and maximum targeted temperatures.
58 :
59 : The algorithm works as follows.
60 : The target distribution for the potential energy is chosen to be:
61 : \f[
62 : p(E)= \begin{cases}
63 : \frac{1}{E_2-E_1} & \mathrm{if} \quad E_1<E<E_2 \\
64 : 0 & \mathrm{otherwise}
65 : \end{cases}
66 : \f]
67 : where the energy limits \f$E_1\f$ and \f$E_2\f$ are yet to be determined.
68 : Clearly the interval \f$E_1–E_2\f$ chosen is related to the interval of temperatures \f$T_1-T_2\f$.
69 : To link these two intervals we make use of the following relation:
70 : \f[
71 : \beta' F_{\beta'}(E) = \beta F_{\beta}(E) + (\beta' - \beta) E + C,
72 : \f]
73 : where \f$F_{\beta}(E)\f$ is determined during the optimization and we shall choose \f$C\f$ such that \f$F_{\beta'}(E_{m})=0\f$ with \f$E_{m}\f$ the position of the free energy minimum.
74 : Using this relation we employ an iterative procedure to find the energy interval.
75 : At iteration \f$k\f$ we have the estimates \f$E_1^k\f$ and \f$E_2^k\f$ for \f$E_1\f$ and \f$E_2\f$, and the target distribution is:
76 : \f[
77 : p^k(E)=\frac{1}{E_2^k-E_1^k} \quad \mathrm{for} \quad E_1^k<E<E_2^k.
78 : \f]
79 : \f$E_1^k\f$ and \f$E_2^k\f$ are obtained from the leftmost solution of \f$\beta_2 F_{\beta_2}^{k-1}(E_1^k)=\epsilon\f$ and the rightmost solution of \f$\beta_1 F_{\beta_1}^{k-1}(E_2^k)=\epsilon\f$.
80 : The procedure is repeated until convergence.
81 : This iterative approach is similar to that in \ref TD_WELLTEMPERED.
82 :
83 : The version of this algorithm in which the energy and an order parameter are biased is similar to the one described in \ref TD_MULTITHERMAL_MULTIBARIC.
84 :
85 : The output of these simulations can be reweighted in order to obtain information at all temperatures in the targeted temperature interval.
86 : The reweighting can be performed using the action \ref REWEIGHT_TEMP_PRESS.
87 :
88 : \par Examples
89 :
90 : The following input can be used to run a simulation in the multicanonical ensemble.
91 : The temperature interval to be explored is 400-600 K.
92 : The energy is used as collective variable.
93 : Legendre polynomials are used to construct the bias potential.
94 : The averaged stochastic gradient descent algorithm is chosen to optimize the VES functional.
95 : The target distribution is updated every 100 optimization steps (200 ps here) using the last estimation of the free energy.
96 :
97 : \plumedfile
98 : # Use energy and volume as CVs
99 : energy: ENERGY
100 :
101 : # Basis functions
102 : bf1: BF_LEGENDRE ORDER=20 MINIMUM=-25000 MAXIMUM=-23500
103 :
104 : # Target distributions
105 : TD_MULTICANONICAL ...
106 : LABEL=td_multi
107 : SIGMA=50.0
108 : MIN_TEMP=400
109 : MAX_TEMP=600
110 : THRESHOLD=1
111 : ... TD_MULTICANONICAL
112 :
113 : # Expansion
114 : VES_LINEAR_EXPANSION ...
115 : ARG=energy
116 : BASIS_FUNCTIONS=bf1
117 : TEMP=500.0
118 : GRID_BINS=1000
119 : TARGET_DISTRIBUTION=td_multi
120 : LABEL=b1
121 : ... VES_LINEAR_EXPANSION
122 :
123 : # Optimization algorithm
124 : OPT_AVERAGED_SGD ...
125 : BIAS=b1
126 : STRIDE=500
127 : LABEL=o1
128 : STEPSIZE=1.0
129 : FES_OUTPUT=500
130 : BIAS_OUTPUT=500
131 : TARGETDIST_OUTPUT=500
132 : COEFFS_OUTPUT=10
133 : TARGETDIST_STRIDE=100
134 : ... OPT_AVERAGED_SGD
135 :
136 : \endplumedfile
137 :
138 : The multicanonical target distribution can also be used to explore a temperature interval in which a first order phase transitions is observed.
139 :
140 : */
141 : //+ENDPLUMEDOC
142 :
143 : class TD_Multicanonical: public TargetDistribution {
144 : private:
145 : double threshold_, min_temp_, max_temp_;
146 : std::vector<double> sigma_;
147 : unsigned steps_temp_;
148 : public:
149 : static void registerKeywords(Keywords&);
150 : explicit TD_Multicanonical(const ActionOptions& ao);
151 : void updateGrid();
152 : double getValue(const std::vector<double>&) const;
153 6 : ~TD_Multicanonical() {}
154 : double GaussianSwitchingFunc(const double, const double, const double) const;
155 : };
156 :
157 :
158 7836 : PLUMED_REGISTER_ACTION(TD_Multicanonical,"TD_MULTICANONICAL")
159 :
160 :
161 3 : void TD_Multicanonical::registerKeywords(Keywords& keys) {
162 3 : TargetDistribution::registerKeywords(keys);
163 15 : keys.add("compulsory","THRESHOLD","1","Maximum exploration free energy in kT.");
164 12 : keys.add("compulsory","MIN_TEMP","Minimum temperature.");
165 12 : keys.add("compulsory","MAX_TEMP","Maximum temperature.");
166 12 : keys.add("optional","STEPS_TEMP","Number of temperature steps. Only for the 2D version, i.e. energy and order parameter.");
167 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.");
168 3 : }
169 :
170 :
171 2 : TD_Multicanonical::TD_Multicanonical(const ActionOptions& ao):
172 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
173 : threshold_(1.0),
174 : min_temp_(0.0),
175 : max_temp_(1000.0),
176 : sigma_(0.0),
177 2 : steps_temp_(20)
178 : {
179 2 : log.printf(" Multicanonical target distribution");
180 2 : log.printf("\n");
181 2 : log.printf(" Please read and cite ");
182 6 : log << plumed.cite("Piaggi and Parrinello, Phys. Rev. Lett. 122 (5), 050601 (2019)");
183 2 : log.printf(" and ");
184 6 : log << plumed.cite("Piaggi and Parrinello, arXiv preprint arXiv:1904.05624 (2019)");
185 2 : log.printf("\n");
186 4 : parse("THRESHOLD",threshold_);
187 2 : if(threshold_<=0.0) {
188 0 : plumed_merror("TD_MULTICANONICAL target distribution: the value of the threshold should be positive.");
189 : }
190 :
191 4 : parse("MIN_TEMP",min_temp_);
192 4 : parse("MAX_TEMP",max_temp_);
193 4 : parseVector("SIGMA",sigma_);
194 4 : if(sigma_.size()<1 || sigma_.size()>2) plumed_merror(getName()+": SIGMA takes 1 or 2 values as input.");
195 4 : parse("STEPS_TEMP",steps_temp_); // Only used in the 2D version
196 2 : steps_temp_ += 1;
197 :
198 : setDynamic();
199 : setFesGridNeeded();
200 2 : checkRead();
201 2 : }
202 :
203 :
204 0 : double TD_Multicanonical::getValue(const std::vector<double>& argument) const {
205 0 : plumed_merror("getValue not implemented for TD_Multicanonical");
206 : return 0.0;
207 : }
208 :
209 :
210 14 : void TD_Multicanonical::updateGrid() {
211 14 : if (getStep() == 0) {
212 2 : if(targetDistGrid().getDimension()>2 && targetDistGrid().getDimension()<1) plumed_merror(getName()+" works only with 1 or 2 arguments, i.e. energy, or energy and CV");
213 2 : if(sigma_.size()!=targetDistGrid().getDimension()) plumed_merror(getName()+": mismatch between SIGMA dimension and number of arguments");
214 : // Use uniform TD
215 8 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
216 : double norm = 0.0;
217 5408 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
218 : double value = 1.0;
219 2702 : norm += integration_weights[l]*value;
220 2702 : targetDistGrid().setValue(l,value);
221 : }
222 4 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
223 2 : logTargetDistGrid().setMinToZero();
224 : } else {
225 : // Two variants: 1D and 2D
226 : // This could be done with one variant but the 1D variant is useful for pedagogical purposes.
227 12 : if(targetDistGrid().getDimension()==1) {
228 : // 1D variant: Multicanonical without order parameter
229 : // In this variant we find the minimum and maximum relevant potential energies.
230 : // Using this information we construct a uniform target distribution inbetween these two.
231 10 : double beta = getBeta();
232 20 : double beta_prime_min = 1./(plumed.getAtoms().getKBoltzmann()*min_temp_);
233 20 : double beta_prime_max = 1./(plumed.getAtoms().getKBoltzmann()*max_temp_);
234 10 : plumed_massert(getFesGridPntr()!=NULL,"the FES grid has to be linked to use TD_Multicanonical!");
235 : // Find minimum of F(U) at temperature min
236 : double minval=DBL_MAX;
237 10 : Grid::index_t minindex = (targetDistGrid().getSize())/2;
238 30 : double minpos = targetDistGrid().getPoint(minindex)[0];
239 2040 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
240 1010 : double value = getFesGridPntr()->getValue(l);
241 3030 : double argument = targetDistGrid().getPoint(l)[0];
242 1010 : value = beta*value + (beta_prime_min-beta)*argument;
243 1010 : if(value<minval) {
244 : minval=value;
245 : minpos=argument;
246 : minindex=l;
247 : }
248 : }
249 : // Find minimum energy at low temperature
250 : double minimum_low = minpos;
251 12 : for(Grid::index_t l=minindex; l>1; l-=1) {
252 36 : double argument = targetDistGrid().getPoint(l)[0];
253 48 : double argument_next = targetDistGrid().getPoint(l-1)[0];
254 12 : double value = getFesGridPntr()->getValue(l);
255 12 : double value_next = getFesGridPntr()->getValue(l-1);
256 12 : value = beta*value + (beta_prime_min-beta)*argument - minval;
257 12 : value_next = beta*value_next + (beta_prime_min-beta)*argument_next - minval;
258 12 : if (value<threshold_ && value_next>threshold_) {
259 : minimum_low = argument_next;
260 : break;
261 : }
262 : }
263 : // Find maximum energy at low temperature
264 : double maximum_low = minpos;
265 22 : for(Grid::index_t l=minindex; l<(targetDistGrid().getSize()-1); l++) {
266 36 : double argument = targetDistGrid().getPoint(l)[0];
267 48 : double argument_next = targetDistGrid().getPoint(l+1)[0];
268 12 : double value = getFesGridPntr()->getValue(l);
269 12 : double value_next = getFesGridPntr()->getValue(l+1);
270 12 : value = beta*value + (beta_prime_min-beta)*argument - minval;
271 12 : value_next = beta*value_next + (beta_prime_min-beta)*argument_next - minval;
272 12 : if (value<threshold_ && value_next>threshold_) {
273 : maximum_low = argument_next;
274 : break;
275 : }
276 : }
277 : // Find minimum of F(U) at temperature max
278 : minval=DBL_MAX;
279 2040 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
280 1010 : double value = getFesGridPntr()->getValue(l);
281 3030 : double argument = targetDistGrid().getPoint(l)[0];
282 1010 : value = beta*value + (beta_prime_max-beta)*argument;
283 1010 : if(value<minval) {
284 : minval=value;
285 : minpos=argument;
286 : minindex=l;
287 : }
288 : }
289 : // Find minimum energy at high temperature
290 : double minimum_high = minpos;
291 13 : for(Grid::index_t l=minindex; l>1; l-=1) {
292 39 : double argument = targetDistGrid().getPoint(l)[0];
293 52 : double argument_next = targetDistGrid().getPoint(l-1)[0];
294 13 : double value = getFesGridPntr()->getValue(l);
295 13 : double value_next = getFesGridPntr()->getValue(l-1);
296 13 : value = beta*value + (beta_prime_max-beta)*argument - minval;
297 13 : value_next = beta*value_next + (beta_prime_max-beta)*argument_next - minval;
298 13 : if (value<threshold_ && value_next>threshold_) {
299 : minimum_high = argument_next;
300 : break;
301 : }
302 : }
303 : // Find maximum energy at high temperature
304 : double maximum_high = minpos;
305 21 : for(Grid::index_t l=minindex; l<(targetDistGrid().getSize()-1); l++) {
306 33 : double argument = targetDistGrid().getPoint(l)[0];
307 44 : double argument_next = targetDistGrid().getPoint(l+1)[0];
308 11 : double value = getFesGridPntr()->getValue(l);
309 11 : double value_next = getFesGridPntr()->getValue(l+1);
310 11 : value = beta*value + (beta_prime_max-beta)*argument - minval;
311 11 : value_next = beta*value_next + (beta_prime_max-beta)*argument_next - minval;
312 11 : if (value<threshold_ && value_next>threshold_) {
313 : maximum_high = argument_next;
314 : break;
315 : }
316 : }
317 : double minimum = minimum_low;
318 10 : if (minimum_high<minimum_low) minimum=minimum_high;
319 : double maximum = maximum_low;
320 10 : if (maximum_high>maximum_low) maximum=maximum_high;
321 : // Construct uniform TD in the interval between minimum and maximum
322 40 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
323 : double norm = 0.0;
324 2040 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
325 3030 : double argument = targetDistGrid().getPoint(l)[0];
326 : double value = 1.0;
327 : double tmp;
328 1010 : if(argument < minimum) {
329 209 : tmp = GaussianSwitchingFunc(argument,minimum,sigma_[0]);
330 : }
331 801 : else if(argument > maximum) {
332 188 : tmp = GaussianSwitchingFunc(argument,maximum,sigma_[0]);
333 : }
334 : else {
335 : tmp = 1.0;
336 : }
337 : value *= tmp;
338 1010 : norm += integration_weights[l]*value;
339 1010 : targetDistGrid().setValue(l,value);
340 : }
341 20 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
342 10 : logTargetDistGrid().setMinToZero();
343 2 : } else if(targetDistGrid().getDimension()==2) {
344 : // 2D variant: Multicanonical with order parameter
345 : // In this variant we find for each temperature the relevant region of potential energy and order parameter.
346 : // The target distribution will be the union of the relevant regions at all temperatures in the temperature interval.
347 2 : double beta = getBeta();
348 4 : double beta_prime_min = 1./(plumed.getAtoms().getKBoltzmann()*min_temp_);
349 4 : double beta_prime_max = 1./(plumed.getAtoms().getKBoltzmann()*max_temp_);
350 2 : plumed_massert(getFesGridPntr()!=NULL,"the FES grid has to be linked to use TD_MulticanonicalWithCV!");
351 : // Set all to zero
352 10406 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
353 : double value = 0.0;
354 5202 : targetDistGrid().setValue(l,value);
355 : }
356 : // Loop over temperatures
357 42 : for(unsigned i=0; i<steps_temp_; i++) {
358 42 : double beta_prime=beta_prime_min + (beta_prime_max-beta_prime_min)*i/(steps_temp_-1);
359 : // Find minimum for this temperature
360 : double minval=DBL_MAX;
361 218568 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
362 327726 : double energy = targetDistGrid().getPoint(l)[0];
363 109242 : double value = getFesGridPntr()->getValue(l);
364 109242 : value = beta*value + (beta_prime-beta)*energy;
365 109242 : if(value<minval) {
366 : minval=value;
367 : }
368 : }
369 : // Now check which energies and volumes are below X kt
370 218526 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
371 327726 : double energy = targetDistGrid().getPoint(l)[0];
372 109242 : double value = getFesGridPntr()->getValue(l);
373 109242 : value = beta*value + (beta_prime-beta)*energy - minval;
374 109242 : if (value<threshold_) {
375 : double value = 1.0;
376 7151 : targetDistGrid().setValue(l,value);
377 : }
378 : }
379 : }
380 2 : std::vector<unsigned> nbin=targetDistGrid().getNbin();
381 2 : std::vector<double> dx=targetDistGrid().getDx();
382 : // Smoothening
383 206 : for(unsigned i=0; i<nbin[0]; i++) {
384 10506 : for(unsigned j=0; j<nbin[1]; j++) {
385 5202 : std::vector<unsigned> indices(2);
386 5202 : indices[0]=i;
387 5202 : indices[1]=j;
388 5202 : Grid::index_t index = targetDistGrid().getIndex(indices);
389 15606 : double energy = targetDistGrid().getPoint(index)[0];
390 15606 : double volume = targetDistGrid().getPoint(index)[1];
391 5202 : double value = targetDistGrid().getValue(index);
392 5202 : if (value>(1-1.e-5)) { // Apply only if this grid point was 1.
393 : // Apply gaussians around
394 777 : std::vector<int> minBin(2), maxBin(2), deltaBin(2); // These cannot be unsigned
395 : // Only consider contributions less than n*sigma bins apart from the actual distance
396 1554 : deltaBin[0]=std::floor(5*sigma_[0]/dx[0]);;
397 1554 : deltaBin[1]=std::floor(5*sigma_[1]/dx[1]);;
398 : // For energy
399 777 : minBin[0]=i - deltaBin[0];
400 777 : if (minBin[0] < 0) minBin[0]=0;
401 1554 : if (minBin[0] > (nbin[0]-1)) minBin[0]=nbin[0]-1;
402 777 : maxBin[0]=i + deltaBin[0];
403 1554 : if (maxBin[0] > (nbin[0]-1)) maxBin[0]=nbin[0]-1;
404 : // For volume
405 777 : minBin[1]=j - deltaBin[1];
406 777 : if (minBin[1] < 0) minBin[1]=0;
407 1554 : if (minBin[1] > (nbin[1]-1)) minBin[1]=nbin[1]-1;
408 777 : maxBin[1]=j + deltaBin[1];
409 1554 : if (maxBin[1] > (nbin[1]-1)) maxBin[1]=nbin[1]-1;
410 58472 : for(unsigned l=minBin[0]; l<maxBin[0]+1; l++) {
411 966356 : for(unsigned m=minBin[1]; m<maxBin[1]+1; m++) {
412 454719 : std::vector<unsigned> indices_prime(2);
413 454719 : indices_prime[0]=l;
414 454719 : indices_prime[1]=m;
415 454719 : Grid::index_t index_prime = targetDistGrid().getIndex(indices_prime);
416 1364157 : double energy_prime = targetDistGrid().getPoint(index_prime)[0];
417 1364157 : double volume_prime = targetDistGrid().getPoint(index_prime)[1];
418 454719 : double value_prime = targetDistGrid().getValue(index_prime);
419 : // Apply gaussian
420 1364157 : double gaussian_value = GaussianSwitchingFunc(energy_prime,energy,sigma_[0])*GaussianSwitchingFunc(volume_prime,volume,sigma_[1]);
421 454719 : if (value_prime<gaussian_value) {
422 39745 : targetDistGrid().setValue(index_prime,gaussian_value);
423 : }
424 : }
425 : }
426 : }
427 : }
428 : }
429 : // Normalize
430 8 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
431 : double norm = 0.0;
432 10408 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
433 5202 : double value = targetDistGrid().getValue(l);
434 5202 : norm += integration_weights[l]*value;
435 : }
436 4 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
437 2 : logTargetDistGrid().setMinToZero();
438 0 : } else plumed_merror(getName()+": Number of arguments for this target distribution must be 1 or 2");
439 : }
440 14 : }
441 :
442 : inline
443 : double TD_Multicanonical::GaussianSwitchingFunc(const double argument, const double center, const double sigma) const {
444 909835 : if(sigma>0.0) {
445 909835 : double arg=(argument-center)/sigma;
446 909835 : return exp(-0.5*arg*arg);
447 : }
448 : else {
449 : return 0.0;
450 : }
451 : }
452 :
453 :
454 : }
455 5874 : }
|
