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 :
29 : #include "lepton/Lepton.h"
30 :
31 :
32 : namespace PLMD {
33 : namespace ves {
34 :
35 17677 : static std::map<std::string, double> leptonConstants= {
36 : {"e", std::exp(1.0)},
37 : {"log2e", 1.0/std::log(2.0)},
38 : {"log10e", 1.0/std::log(10.0)},
39 : {"ln2", std::log(2.0)},
40 : {"ln10", std::log(10.0)},
41 : {"pi", pi},
42 : {"pi_2", pi*0.5},
43 : {"pi_4", pi*0.25},
44 : // {"1_pi", 1.0/pi},
45 : // {"2_pi", 2.0/pi},
46 : // {"2_sqrtpi", 2.0/std::sqrt(pi)},
47 : {"sqrt2", std::sqrt(2.0)},
48 : {"sqrt1_2", std::sqrt(0.5)}
49 16070 : };
50 :
51 :
52 : //+PLUMEDOC VES_TARGETDIST TD_CUSTOM
53 : /*
54 : Target distribution given by an arbitrary mathematical expression (static or dynamic).
55 :
56 : Use as a target distribution the distribution defined by
57 : \f[
58 : p(\mathbf{s}) =
59 : \frac{f(\mathbf{s})}{\int d\mathbf{s} \, f(\mathbf{s})}
60 : \f]
61 : where \f$f(\mathbf{s})\f$ is some arbitrary mathematical function that
62 : is parsed by the lepton library.
63 :
64 : The function \f$f(\mathbf{s})\f$ is given by the FUNCTION keywords by
65 : using _s1_,_s2_,..., as variables for the arguments
66 : \f$\mathbf{s}=(s_1,s_2,\ldots,s_d)\f$.
67 : If one variable is not given the target distribution will be
68 : taken as uniform in that argument.
69 :
70 : It is also possible to include the free energy surface \f$F(\mathbf{s})\f$
71 : in the target distribution by using the _FE_ variable. In this case the
72 : target distribution is dynamic and needs to be updated with current
73 : best estimate of \f$F(\mathbf{s})\f$, similarly as for the
74 : \ref TD_WELLTEMPERED "well-tempered target distribution".
75 : Furthermore, the inverse temperature \f$\beta = (k_{\mathrm{B}}T)^{-1}\f$ and
76 : the thermal energy \f$k_{\mathrm{B}}T\f$ can be included
77 : by using the _beta_ and _kBT_ variables.
78 :
79 : The target distribution will be automatically normalized over the region on
80 : which it is defined on. Therefore, the function given in
81 : FUNCTION needs to be non-negative and normalizable. The
82 : code will perform checks to make sure that this is indeed the case.
83 :
84 :
85 : \par Examples
86 :
87 : Here we use as shifted [Maxwell-Boltzmann distribution](https://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_distribution)
88 : as a target distribution in one-dimension.
89 : Note that it is not need to include the normalization factor as the distribution will be
90 : automatically normalized.
91 : \plumedfile
92 : TD_CUSTOM ...
93 : FUNCTION=(s1+20)^2*exp(-(s1+20)^2/(2*10.0^2))
94 : LABEL=td
95 : ... TD_CUSTOM
96 : \endplumedfile
97 :
98 : Here we have a two dimensional target distribution where we
99 : use a [generalized normal distribution](https://en.wikipedia.org/wiki/Generalized_normal_distribution)
100 : for argument \f$s_2\f$ while the distribution for \f$s_1\f$ is taken as
101 : uniform as the variable _s1_ is not included in the function.
102 : \plumedfile
103 : TD_CUSTOM ...
104 : FUNCTION=exp(-(abs(s2-20.0)/5.0)^4.0)
105 : LABEL=td
106 : ... TD_CUSTOM
107 : \endplumedfile
108 :
109 : By using the _FE_ variable the target distribution can depend on
110 : the free energy surface \f$F(\mathbf{s})\f$. For example,
111 : the following input is identical to using \ref TD_WELLTEMPERED with
112 : BIASFACTOR=10.
113 : \plumedfile
114 : TD_CUSTOM ...
115 : FUNCTION=exp(-(beta/10.0)*FE)
116 : LABEL=td
117 : ... TD_CUSTOM
118 : \endplumedfile
119 : Here the inverse temperature is automatically obtained by using the _beta_
120 : variable. It is also possible to use the _kBT_ variable. The following
121 : syntax will give the exact same results as the syntax above
122 : \plumedfile
123 : TD_CUSTOM ...
124 : FUNCTION=exp(-(1.0/(kBT*10.0))*FE)}
125 : LABEL=td
126 : ... TD_CUSTOM
127 : \endplumedfile
128 :
129 :
130 : */
131 : //+ENDPLUMEDOC
132 :
133 : class TD_Custom : public TargetDistribution {
134 : private:
135 : void setupAdditionalGrids(const std::vector<Value*>&, const std::vector<std::string>&, const std::vector<std::string>&, const std::vector<unsigned int>&);
136 : //
137 : lepton::CompiledExpression expression;
138 : //
139 : std::vector<unsigned int> cv_var_idx_;
140 : std::vector<std::string> cv_var_str_;
141 : //
142 : std::string cv_var_prefix_str_;
143 : std::string fes_var_str_;
144 : std::string kbt_var_str_;
145 : std::string beta_var_str_;
146 : //
147 : bool use_fes_;
148 : bool use_kbt_;
149 : bool use_beta_;
150 : public:
151 : static void registerKeywords( Keywords&);
152 : explicit TD_Custom(const ActionOptions& ao);
153 : void updateGrid();
154 : double getValue(const std::vector<double>&) const;
155 32 : ~TD_Custom() {};
156 : };
157 :
158 4837 : PLUMED_REGISTER_ACTION(TD_Custom,"TD_CUSTOM")
159 :
160 :
161 17 : void TD_Custom::registerKeywords(Keywords& keys) {
162 17 : TargetDistribution::registerKeywords(keys);
163 17 : keys.add("compulsory","FUNCTION","The function you wish to use for the target distribution where you should use the variables _s1_,_s2_,... for the arguments. You can also use the current estimate of the FES by using the variable _FE_ and the temperature by using the _kBT_ and _beta_ variables.");
164 17 : keys.use("WELLTEMPERED_FACTOR");
165 17 : keys.use("SHIFT_TO_ZERO");
166 17 : }
167 :
168 :
169 16 : TD_Custom::TD_Custom(const ActionOptions& ao):
170 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
171 : //
172 : cv_var_idx_(0),
173 : cv_var_str_(0),
174 : //
175 : cv_var_prefix_str_("s"),
176 : fes_var_str_("FE"),
177 : kbt_var_str_("kBT"),
178 : beta_var_str_("beta"),
179 : //
180 : use_fes_(false),
181 : use_kbt_(false),
182 16 : use_beta_(false)
183 : {
184 16 : std::string func_str;
185 16 : parse("FUNCTION",func_str);
186 16 : checkRead();
187 : //
188 : try {
189 16 : lepton::ParsedExpression pe=lepton::Parser::parse(func_str).optimize(leptonConstants);
190 16 : log<<" function as parsed by lepton: "<<pe<<"\n";
191 16 : expression=pe.createCompiledExpression();
192 : }
193 0 : catch(PLMD::lepton::Exception& exc) {
194 0 : plumed_merror("There was some problem in parsing the function "+func_str+" given in FUNCTION with lepton");
195 : }
196 :
197 39 : for(auto &p: expression.getVariables()) {
198 23 : std::string curr_var = p;
199 : unsigned int cv_idx;
200 23 : if(curr_var.substr(0,cv_var_prefix_str_.size())==cv_var_prefix_str_ && Tools::convert(curr_var.substr(cv_var_prefix_str_.size()),cv_idx) && cv_idx>0) {
201 16 : cv_var_idx_.push_back(cv_idx-1);
202 : }
203 7 : else if(curr_var==fes_var_str_) {
204 3 : use_fes_=true;
205 3 : setDynamic();
206 3 : setFesGridNeeded();
207 : }
208 4 : else if(curr_var==kbt_var_str_) {
209 2 : use_kbt_=true;
210 : }
211 2 : else if(curr_var==beta_var_str_) {
212 2 : use_beta_=true;
213 : }
214 : else {
215 0 : plumed_merror(getName()+": problem with parsing formula with lepton, cannot recognise the variable "+curr_var);
216 : }
217 23 : }
218 : //
219 16 : std::sort(cv_var_idx_.begin(),cv_var_idx_.end());
220 16 : cv_var_str_.resize(cv_var_idx_.size());
221 32 : for(unsigned int j=0; j<cv_var_idx_.size(); j++) {
222 16 : std::string str1; Tools::convert(cv_var_idx_[j]+1,str1);
223 16 : cv_var_str_[j] = cv_var_prefix_str_+str1;
224 32 : }
225 16 : }
226 :
227 :
228 16 : void TD_Custom::setupAdditionalGrids(const std::vector<Value*>& arguments, const std::vector<std::string>& min, const std::vector<std::string>& max, const std::vector<unsigned int>& nbins) {
229 16 : if(cv_var_idx_.size()>0 && cv_var_idx_[cv_var_idx_.size()-1]>getDimension()) {
230 0 : plumed_merror(getName()+": mismatch between CVs given in FUNC and the dimension of the target distribution");
231 : }
232 16 : }
233 :
234 :
235 0 : double TD_Custom::getValue(const std::vector<double>& argument) const {
236 0 : plumed_merror("getValue not implemented for TD_Custom");
237 : return 0.0;
238 : }
239 :
240 :
241 46 : void TD_Custom::updateGrid() {
242 46 : if(use_fes_) {
243 33 : plumed_massert(getFesGridPntr()!=NULL,"the FES grid has to be linked to the free energy in the target distribution");
244 : }
245 46 : if(use_kbt_) {
246 : try {
247 22 : expression.getVariableReference(kbt_var_str_) = 1.0/getBeta();
248 0 : } catch(PLMD::lepton::Exception& exc) {}
249 : }
250 46 : if(use_beta_) {
251 : try {
252 22 : expression.getVariableReference(beta_var_str_) = getBeta();
253 0 : } catch(PLMD::lepton::Exception& exc) {}
254 : }
255 : //
256 46 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
257 46 : double norm = 0.0;
258 : //
259 40909 : for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
260 40863 : std::vector<double> point = targetDistGrid().getPoint(l);
261 99928 : for(unsigned int k=0; k<cv_var_str_.size() ; k++) {
262 : try {
263 59065 : expression.getVariableReference(cv_var_str_[k]) = point[cv_var_idx_[k]];
264 0 : } catch(PLMD::lepton::Exception& exc) {}
265 : }
266 40863 : if(use_fes_) {
267 : try {
268 3300 : expression.getVariableReference(fes_var_str_) = getFesGridPntr()->getValue(l);
269 0 : } catch(PLMD::lepton::Exception& exc) {}
270 : }
271 40863 : double value = expression.evaluate();
272 :
273 40863 : if(value<0.0 && !isTargetDistGridShiftedToZero()) {plumed_merror(getName()+": The target distribution function gives negative values. You should change the definition of the function used for the target distribution to avoid this. You can also use the SHIFT_TO_ZERO keyword to avoid this problem.");}
274 40863 : targetDistGrid().setValue(l,value);
275 40863 : norm += integration_weights[l]*value;
276 40863 : logTargetDistGrid().setValue(l,-std::log(value));
277 40863 : }
278 46 : if(norm>0.0) {
279 45 : targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
280 : }
281 1 : else if(!isTargetDistGridShiftedToZero()) {
282 0 : plumed_merror(getName()+": The target distribution function cannot be normalized proberly. You should change the definition of the function used for the target distribution to avoid this. You can also use the SHIFT_TO_ZERO keyword to avoid this problem.");
283 : }
284 46 : logTargetDistGrid().setMinToZero();
285 46 : }
286 :
287 :
288 : }
289 4821 : }
|
