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 :
26 : #include "core/ActionRegister.h"
27 : #include "tools/Grid.h"
28 : #include "core/PlumedMain.h"
29 :
30 :
31 :
32 : namespace PLMD {
33 : namespace ves {
34 :
35 : //+PLUMEDOC VES_TARGETDIST TD_WELLTEMPERED
36 : /*
37 : Well-tempered target distribution (dynamic).
38 :
39 : Use as a target distribution the well-tempered distribution \cite Barducci:2008
40 : given by
41 : \f[
42 : p(\mathbf{s}) =
43 : \frac{e^{-(\beta/\gamma) F(\mathbf{s})}}
44 : {\int d\mathbf{s}\, e^{-(\beta/\gamma) F(\mathbf{s})}} =
45 : \frac{[P_{0}(\mathbf{s})]^{1/\gamma}}
46 : {\int d\mathbf{s}\, [P_{0}(\mathbf{s})]^{1/\gamma}}
47 : \f]
48 : where \f$\gamma\f$ is a so-called bias factor and \f$P_{0}(\mathbf{s})\f$ is the
49 : unbiased canonical distribution of the CVs. This target distribution thus
50 : correponds to a biased ensemble where, as compared to the unbiased one,
51 : the probability peaks have been broaden and the fluctations of the CVs are
52 : enhanced.
53 : The value of the bias factor \f$\gamma\f$ determines by how much the fluctations
54 : are enhanced.
55 :
56 : The well-tempered distribution can be view as sampling on
57 : an effective free energy surface \f$\tilde{F}(\mathbf{s}) = (1/\gamma) F(\mathbf{s})\f$
58 : which has largely the same metastable states as the original \f$F(\mathbf{s})\f$
59 : but with barriers that have been reduced by a factor of \f$\gamma\f$.
60 : Generally one should use a value of \f$\gamma\f$ that results in
61 : effective barriers on the order of few \f$k_{\mathrm{B}}T\f$
62 : such that thermal fluctuations can easily induce transitions
63 : between different metastable states.
64 :
65 : At convergence the relationship between the bias potential and the free
66 : energy surface is given by
67 : \f[
68 : F(\mathbf{s}) = - \left(\frac{1}{1-\gamma^{-1}} \right) V(\mathbf{s})
69 : \f]
70 :
71 : This target distribution depends directly on the free energy surface
72 : \f$F(\mathbf{s})\f$ which is quantity that we do not know a-priori and
73 : want to obtain. Therefore, this target distribution
74 : is iteratively updated \cite Valsson-JCTC-2015 according to
75 : \f[
76 : p^{(m+1)}(\mathbf{s}) =
77 : \frac{e^{-(\beta/\gamma) F^{(m+1)}(\mathbf{s})}}
78 : {\int d\mathbf{s}\, e^{-(\beta/\gamma) F^{(m+1)}(\mathbf{s})}}
79 : \f]
80 : where \f$F^{(m+1)}(\mathbf{s})\f$ is the current best estimate of the
81 : free energy surface obtained according to
82 : \f[
83 : F^{(m+1)}(\mathbf{s}) =
84 : - V^{(m+1)}(\mathbf{s}) - \frac{1}{\beta} \log p^{(m)}(\mathbf{s}) =
85 : - V^{(m+1)}(\mathbf{s}) + \frac{1}{\gamma} F^{(m)}(\mathbf{s})
86 : \f]
87 : The frequency of performing this update needs to be set in the
88 : optimizer used in the calculation. Normally it is sufficient
89 : to do it every 100-1000 bias update iterations.
90 :
91 : \par Examples
92 :
93 : Employ a well-tempered target distribution with a bias factor of 10
94 : \plumedfile
95 : td_welltemp: TD_WELLTEMPERED BIASFACTOR=10
96 : \endplumedfile
97 :
98 : */
99 : //+ENDPLUMEDOC
100 :
101 : class TD_WellTempered: public TargetDistribution {
102 : private:
103 : double bias_factor_;
104 : public:
105 : static void registerKeywords(Keywords&);
106 : explicit TD_WellTempered(const ActionOptions& ao);
107 : void updateGrid();
108 : double getValue(const std::vector<double>&) const;
109 28 : ~TD_WellTempered() {}
110 : };
111 :
112 :
113 6480 : PLUMED_REGISTER_ACTION(TD_WellTempered,"TD_WELLTEMPERED")
114 :
115 :
116 29 : void TD_WellTempered::registerKeywords(Keywords& keys) {
117 29 : TargetDistribution::registerKeywords(keys);
118 116 : keys.add("compulsory","BIASFACTOR","The bias factor used for the well-tempered distribution.");
119 29 : }
120 :
121 :
122 28 : TD_WellTempered::TD_WellTempered(const ActionOptions& ao):
123 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
124 28 : bias_factor_(0.0)
125 : {
126 28 : log.printf(" Well-tempered target distribution, see and cite ");
127 84 : log << plumed.cite("Valsson and Parrinello, J. Chem. Theory Comput. 11, 1996-2002 (2015)");
128 84 : log << plumed.cite("Barducci, Bussi, and Parrinello, Phys. Rev. Lett. 100, 020603 (2008)");
129 28 : log.printf("\n");
130 56 : parse("BIASFACTOR",bias_factor_);
131 28 : if(bias_factor_<=1.0) {
132 0 : plumed_merror("TD_WELLTEMPERED target distribution: the value of the bias factor doesn't make sense, it should be larger than 1.0");
133 : }
134 : setDynamic();
135 : setFesGridNeeded();
136 28 : checkRead();
137 28 : }
138 :
139 :
140 0 : double TD_WellTempered::getValue(const std::vector<double>& argument) const {
141 0 : plumed_merror("getValue not implemented for TD_WellTempered");
142 : return 0.0;
143 : }
144 :
145 :
146 308 : void TD_WellTempered::updateGrid() {
147 308 : double beta_prime = getBeta()/bias_factor_;
148 308 : plumed_massert(getFesGridPntr()!=NULL,"the FES grid has to be linked to use TD_WellTempered!");
149 1232 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
150 : double norm = 0.0;
151 2253240 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
152 1126466 : double value = beta_prime * getFesGridPntr()->getValue(l);
153 1126466 : logTargetDistGrid().setValue(l,value);
154 1126466 : value = exp(-value);
155 1126466 : norm += integration_weights[l]*value;
156 1126466 : targetDistGrid().setValue(l,value);
157 : }
158 308 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
159 308 : logTargetDistGrid().setMinToZero();
160 308 : }
161 :
162 :
163 : }
164 4839 : }
|
