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 "BasisFunctions.h"
24 : #include "TargetDistribution.h"
25 : #include "VesBias.h"
26 : #include "VesTools.h"
27 : #include "GridIntegrationWeights.h"
28 :
29 : #include "tools/Grid.h"
30 : #include "tools/Tools.h"
31 :
32 :
33 : namespace PLMD {
34 : namespace ves {
35 :
36 188 : void BasisFunctions::registerKeywords(Keywords& keys) {
37 188 : Action::registerKeywords(keys);
38 752 : keys.add("compulsory","ORDER","The order of the basis function expansion.");
39 752 : keys.add("compulsory","MINIMUM","The minimum of the interval on which the basis functions are defined.");
40 752 : keys.add("compulsory","MAXIMUM","The maximum of the interval on which the basis functions are defined.");
41 752 : keys.add("hidden","NGRID_POINTS","The number of grid points used for numerical integrals");
42 564 : keys.addFlag("DEBUG_INFO",false,"Print out more detailed information about the basis set. Useful for debugging.");
43 564 : keys.addFlag("NUMERICAL_INTEGRALS",false,"Calculate basis function integral for the uniform distribution numerically. Useful for debugging.");
44 188 : }
45 :
46 :
47 180 : BasisFunctions::BasisFunctions(const ActionOptions&ao):
48 : Action(ao),
49 : print_debug_info_(false),
50 : has_been_set(false),
51 : description_("Undefined"),
52 : type_("Undefined"),
53 : norder_(0),
54 : nbasis_(1),
55 : bf_label_prefix_("f"),
56 180 : bf_labels_(nbasis_,"f0"),
57 : periodic_(false),
58 : interval_bounded_(true),
59 : interval_intrinsic_min_str_("1.0"),
60 : interval_intrinsic_max_str_("-1.0"),
61 : interval_intrinsic_min_(1.0),
62 : interval_intrinsic_max_(-1.0),
63 : interval_intrinsic_range_(0.0),
64 : interval_intrinsic_mean_(0.0),
65 : interval_min_str_(""),
66 : interval_max_str_(""),
67 : interval_min_(0.0),
68 : interval_max_(0.0),
69 : interval_range_(0.0),
70 : interval_mean_(0.0),
71 : argT_derivf_(1.0),
72 : numerical_uniform_integrals_(false),
73 : nbins_(1001),
74 180 : uniform_integrals_(nbasis_,0.0),
75 : vesbias_pntr_(NULL),
76 900 : action_pntr_(NULL)
77 : {
78 180 : bf_keywords_.push_back(getName());
79 360 : if(keywords.exists("ORDER")) {
80 504 : parse("ORDER",norder_); addKeywordToList("ORDER",norder_);
81 : }
82 180 : nbasis_=norder_+1;
83 : //
84 : std::string str_imin; std::string str_imax;
85 540 : if(keywords.exists("MINIMUM") && keywords.exists("MAXIMUM")) {
86 890 : parse("MINIMUM",str_imin); addKeywordToList("MINIMUM",str_imin);
87 890 : parse("MAXIMUM",str_imax); addKeywordToList("MAXIMUM",str_imax);
88 : }
89 : else {
90 : str_imin = "-1.0";
91 : str_imax = "1.0";
92 : }
93 : interval_min_str_ = str_imin;
94 : interval_max_str_ = str_imax;
95 180 : if(!Tools::convert(str_imin,interval_min_)) {
96 0 : plumed_merror(getName()+": cannot convert the value given in MINIMUM to a double");
97 : }
98 180 : if(!Tools::convert(str_imax,interval_max_)) {
99 0 : plumed_merror(getName()+": cannot convert the value given in MAXIMUM to a double");
100 : }
101 180 : if(interval_min_>interval_max_) {plumed_merror(getName()+": MINIMUM and MAXIMUM are not correctly defined");}
102 : //
103 360 : parseFlag("DEBUG_INFO",print_debug_info_);
104 360 : if(keywords.exists("NUMERICAL_INTEGRALS")) {
105 312 : parseFlag("NUMERICAL_INTEGRALS",numerical_uniform_integrals_);
106 : }
107 360 : if(keywords.exists("NGRID_POINTS")) {
108 356 : parse("NGRID_POINTS",nbins_);
109 : }
110 : // log.printf(" %s \n",getKeywordString().c_str());
111 :
112 180 : }
113 :
114 :
115 22 : void BasisFunctions::setIntrinsicInterval(const double interval_intrinsic_min_in, const double interval_intrinsic_max_in) {
116 22 : interval_intrinsic_min_ = interval_intrinsic_min_in;
117 22 : interval_intrinsic_max_ = interval_intrinsic_max_in;
118 22 : VesTools::convertDbl2Str(interval_intrinsic_min_,interval_intrinsic_min_str_);
119 22 : VesTools::convertDbl2Str(interval_intrinsic_max_,interval_intrinsic_max_str_);
120 22 : plumed_massert(interval_intrinsic_min_<interval_intrinsic_max_,"setIntrinsicInterval: intrinsic intervals are not defined correctly");
121 22 : }
122 :
123 :
124 158 : void BasisFunctions::setIntrinsicInterval(const std::string& interval_intrinsic_min_str_in, const std::string& interval_intrinsic_max_str_in) {
125 158 : interval_intrinsic_min_str_ = interval_intrinsic_min_str_in;
126 158 : interval_intrinsic_max_str_ = interval_intrinsic_max_str_in;
127 158 : if(!Tools::convert(interval_intrinsic_min_str_,interval_intrinsic_min_)) {
128 0 : plumed_merror("setIntrinsicInterval: cannot convert string value given for the minimum of the intrinsic interval to a double");
129 : }
130 158 : if(!Tools::convert(interval_intrinsic_max_str_,interval_intrinsic_max_)) {
131 0 : plumed_merror("setIntrinsicInterval: cannot convert string value given for the maximum of the intrinsic interval to a double");
132 : }
133 158 : plumed_massert(interval_intrinsic_min_<interval_intrinsic_max_,"setIntrinsicInterval: intrinsic intervals are not defined correctly");
134 158 : }
135 :
136 :
137 0 : void BasisFunctions::setInterval(const double interval_min_in, const double interval_max_in) {
138 0 : interval_min_ = interval_min_in;
139 0 : interval_max_ = interval_max_in;
140 0 : VesTools::convertDbl2Str(interval_min_,interval_min_str_);
141 0 : VesTools::convertDbl2Str(interval_max_,interval_max_str_);
142 0 : plumed_massert(interval_min_<interval_max_,"setInterval: intervals are not defined correctly");
143 0 : }
144 :
145 :
146 2 : void BasisFunctions::setInterval(const std::string& interval_min_str_in, const std::string& interval_max_str_in) {
147 2 : interval_min_str_ = interval_min_str_in;
148 2 : interval_max_str_ = interval_max_str_in;
149 2 : if(!Tools::convert(interval_min_str_,interval_min_)) {
150 0 : plumed_merror("setInterval: cannot convert string value given for the minimum of the interval to a double");
151 : }
152 2 : if(!Tools::convert(interval_max_str_,interval_max_)) {
153 0 : plumed_merror("setInterval: cannot convert string value given for the maximum of the interval to a double");
154 : }
155 2 : plumed_massert(interval_min_<interval_max_,"setInterval: intervals are not defined correctly");
156 2 : }
157 :
158 :
159 180 : void BasisFunctions::setupInterval() {
160 : // if(!intervalBounded()){plumed_merror("setupInterval() only works for bounded interval");}
161 180 : interval_intrinsic_range_ = interval_intrinsic_max_-interval_intrinsic_min_;
162 180 : interval_intrinsic_mean_ = 0.5*(interval_intrinsic_max_+interval_intrinsic_min_);
163 180 : interval_range_ = interval_max_-interval_min_;
164 180 : interval_mean_ = 0.5*(interval_max_+interval_min_);
165 180 : argT_derivf_ = interval_intrinsic_range_/interval_range_;
166 180 : }
167 :
168 :
169 67 : void BasisFunctions::setupLabels() {
170 1385 : for(unsigned int i=0; i < nbasis_; i++) {
171 659 : std::string is; Tools::convert(i,is);
172 2636 : bf_labels_[i]=bf_label_prefix_+is+"(s)";
173 : }
174 67 : }
175 :
176 :
177 22 : void BasisFunctions::setupUniformIntegrals() {
178 22 : numerical_uniform_integrals_=true;
179 22 : numericalUniformIntegrals();
180 22 : }
181 :
182 :
183 180 : void BasisFunctions::setupBF() {
184 180 : if(interval_intrinsic_min_>interval_intrinsic_max_) {plumed_merror("setupBF: default intervals are not correctly set");}
185 180 : setupInterval();
186 180 : setupLabels();
187 180 : if(bf_labels_.size()==1) {plumed_merror("setupBF: the labels of the basis functions are not correct.");}
188 180 : if(!numerical_uniform_integrals_) {setupUniformIntegrals();}
189 6 : else {numericalUniformIntegrals();}
190 180 : if(uniform_integrals_.size()==1) {plumed_merror("setupBF: the integrals of the basis functions is not correct.");}
191 360 : if(type_=="Undefined") {plumed_merror("setupBF: the type of the basis function is not defined.");}
192 360 : if(description_=="Undefined") {plumed_merror("setupBF: the description of the basis function is not defined.");}
193 180 : has_been_set=true;
194 180 : printInfo();
195 180 : }
196 :
197 :
198 180 : void BasisFunctions::printInfo() const {
199 180 : if(!has_been_set) {plumed_merror("the basis set has not be setup correctly");}
200 180 : log.printf(" One-dimensional basis set\n");
201 360 : log.printf(" Description: %s\n",description_.c_str());
202 360 : log.printf(" Type: %s\n",type_.c_str());
203 180 : if(periodic_) {log.printf(" The basis functions are periodic\n");}
204 180 : log.printf(" Order of basis set: %u\n",norder_);
205 180 : log.printf(" Number of basis functions: %u\n",nbasis_);
206 : // log.printf(" Interval of basis set: %f to %f\n",interval_min_,interval_max_);
207 360 : log.printf(" Interval of basis set: %s to %s\n",interval_min_str_.c_str(),interval_max_str_.c_str());
208 180 : log.printf(" Description of basis functions:\n");
209 3892 : for(unsigned int i=0; i < nbasis_; i++) {log.printf(" %2u %10s\n",i,bf_labels_[i].c_str());}
210 : //
211 180 : if(print_debug_info_) {
212 34 : log.printf(" Debug information:\n");
213 : // log.printf(" Default interval of basis set: [%f,%f]\n",interval_intrinsic_min_,interval_intrinsic_max_);
214 68 : log.printf(" Intrinsic interval of basis set: [%s,%s]\n",interval_intrinsic_min_str_.c_str(),interval_intrinsic_max_str_.c_str());
215 34 : log.printf(" Intrinsic interval of basis set: range=%f, mean=%f\n",interval_intrinsic_range_,interval_intrinsic_mean_);
216 : // log.printf(" Defined interval of basis set: [%f,%f]\n",interval_min_,interval_max_);
217 68 : log.printf(" Defined interval of basis set: [%s,%s]\n",interval_min_str_.c_str(),interval_max_str_.c_str());
218 34 : log.printf(" Defined interval of basis set: range=%f, mean=%f\n",interval_range_,interval_mean_);
219 34 : log.printf(" Derivative factor due to interval translation: %f\n",argT_derivf_);
220 34 : log.printf(" Integral of basis functions over the interval:\n");
221 34 : if(numerical_uniform_integrals_) {log.printf(" Note: calculated numerically\n");}
222 882 : for(unsigned int i=0; i < nbasis_; i++) {log.printf(" %2u %16.10f\n",i,uniform_integrals_[i]);}
223 34 : log.printf(" --------------------------\n");
224 : }
225 180 : }
226 :
227 :
228 0 : void BasisFunctions::linkVesBias(VesBias* vesbias_pntr_in) {
229 0 : vesbias_pntr_ = vesbias_pntr_in;
230 0 : action_pntr_ = static_cast<Action*>(vesbias_pntr_in);
231 0 : }
232 :
233 :
234 0 : void BasisFunctions::linkAction(Action* action_pntr_in) {
235 0 : action_pntr_ = action_pntr_in;
236 0 : }
237 :
238 :
239 28 : void BasisFunctions::numericalUniformIntegrals() {
240 112 : std::vector<std::string> grid_min(1); grid_min[0]=intervalMinStr();
241 112 : std::vector<std::string> grid_max(1); grid_max[0]=intervalMaxStr();
242 56 : std::vector<unsigned int> grid_bins(1); grid_bins[0]=nbins_;
243 84 : std::vector<Value*> arguments(1); arguments[0]= new Value(NULL,"arg",false);
244 44 : if(arePeriodic()) {arguments[0]->setDomain(intervalMinStr(),intervalMaxStr());}
245 24 : else {arguments[0]->setNotPeriodic();}
246 56 : Grid* uniform_grid = new Grid("uniform",arguments,grid_min,grid_max,grid_bins,false,false);
247 : //
248 28 : double inverse_normalization = 1.0/(intervalMax()-intervalMin());
249 56132 : for(Grid::index_t l=0; l<uniform_grid->getSize(); l++) {
250 28052 : uniform_grid->setValue(l,inverse_normalization);
251 : }
252 56 : uniform_integrals_ = numericalTargetDistributionIntegralsFromGrid(uniform_grid);
253 28 : delete arguments[0]; arguments.clear();
254 28 : delete uniform_grid;
255 28 : }
256 :
257 :
258 144 : std::vector<double> BasisFunctions::numericalTargetDistributionIntegralsFromGrid(const Grid* grid_pntr) const {
259 144 : plumed_massert(grid_pntr!=NULL,"the grid is not defined");
260 144 : plumed_massert(grid_pntr->getDimension()==1,"the target distribution grid should be one dimensional");
261 : //
262 144 : std::vector<double> targetdist_integrals(nbasis_,0.0);
263 432 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(grid_pntr);
264 :
265 125998 : for(Grid::index_t k=0; k < grid_pntr->getSize(); k++) {
266 188781 : double arg = grid_pntr->getPoint(k)[0];
267 62927 : std::vector<double> bf_values(nbasis_);
268 62927 : std::vector<double> bf_derivs(nbasis_);
269 62927 : bool inside=true;
270 62927 : double argT=0.0;
271 62927 : getAllValues(arg,argT,inside,bf_values,bf_derivs);
272 1458925 : for(unsigned int i=0; i < nbasis_; i++) {
273 2093997 : targetdist_integrals[i] += (integration_weights[k] * grid_pntr->getValue(k)) * bf_values[i];
274 : }
275 : }
276 : // assume that the first function is the constant
277 144 : bool inside=true;
278 144 : double argT=0.0;
279 144 : targetdist_integrals[0] = getValue(0.0,0,argT,inside);
280 144 : return targetdist_integrals;
281 : }
282 :
283 :
284 175 : std::vector<double> BasisFunctions::getTargetDistributionIntegrals(const TargetDistribution* targetdist_pntr) const {
285 175 : if(targetdist_pntr==NULL) {
286 : return getUniformIntegrals();
287 : }
288 : else {
289 : Grid* targetdist_grid = targetdist_pntr->getTargetDistGridPntr();
290 116 : return numericalTargetDistributionIntegralsFromGrid(targetdist_grid);
291 : }
292 : }
293 :
294 :
295 118 : std::string BasisFunctions::getKeywordString() const {
296 : std::string str_keywords=bf_keywords_[0];
297 2300 : for(unsigned int i=1; i<bf_keywords_.size(); i++) {str_keywords+=" "+bf_keywords_[i];}
298 118 : return str_keywords;
299 : }
300 :
301 :
302 144 : double BasisFunctions::getValue(const double arg, const unsigned int n, double& argT, bool& inside_range) const {
303 144 : plumed_massert(n<numberOfBasisFunctions(),"getValue: n is outside range of the defined order of the basis set");
304 144 : inside_range=true;
305 144 : std::vector<double> tmp_values(numberOfBasisFunctions());
306 144 : std::vector<double> tmp_derivs(numberOfBasisFunctions());
307 144 : getAllValues(arg, argT, inside_range, tmp_values, tmp_derivs);
308 432 : return tmp_values[n];
309 : }
310 :
311 :
312 3007 : void BasisFunctions::getAllValuesNumericalDerivs(const double arg, double& argT, bool& inside_range, std::vector<double>& values, std::vector<double>& derivs) const {
313 : // use forward difference, unless very close to the boundary
314 : double delta = sqrt(epsilon);
315 6014 : if((arg+delta)>intervalMax()) {
316 : delta *= -1.0;
317 : }
318 3007 : inside_range=true;
319 3007 : std::vector<double> values_delta(numberOfBasisFunctions());
320 3007 : std::vector<double> derivs_dummy(numberOfBasisFunctions());
321 3007 : getAllValues(arg+delta, argT, inside_range, values_delta, derivs_dummy);
322 3007 : getAllValues(arg, argT, inside_range, values, derivs_dummy);
323 79977 : for(unsigned int i=0; i<numberOfBasisFunctions(); i++) {
324 153940 : derivs[i] = (values_delta[i]-values[i])/delta;
325 : }
326 3007 : }
327 :
328 :
329 59 : void BasisFunctions::getMultipleValue(const std::vector<double>& args, std::vector<double>& argsT, std::vector<std::vector<double> >& values, std::vector<std::vector<double> >& derivs, const bool numerical_deriv) const {
330 59 : argsT.resize(args.size());
331 59 : values.clear();
332 59 : derivs.clear();
333 56341 : for(unsigned int i=0; i<args.size(); i++) {
334 18741 : std::vector<double> tmp_values(getNumberOfBasisFunctions());
335 18741 : std::vector<double> tmp_derivs(getNumberOfBasisFunctions());
336 18741 : bool inside_interval=true;
337 18741 : if(!numerical_deriv) {
338 31468 : getAllValues(args[i],argsT[i],inside_interval,tmp_values,tmp_derivs);
339 : } else {
340 3007 : getAllValuesNumericalDerivs(args[i],argsT[i],inside_interval,tmp_values,tmp_derivs);
341 : }
342 18741 : values.push_back(tmp_values);
343 18741 : derivs.push_back(tmp_derivs);
344 : }
345 59 : }
346 :
347 :
348 59 : void BasisFunctions::writeBasisFunctionsToFile(OFile& ofile_values, OFile& ofile_derivs, const std::string& min_in, const std::string& max_in, unsigned int nbins_in, const bool ignore_periodicity, const std::string& output_fmt, const bool numerical_deriv) const {
349 118 : std::vector<std::string> min(1); min[0]=min_in;
350 118 : std::vector<std::string> max(1); max[0]=max_in;
351 118 : std::vector<unsigned int> nbins(1); nbins[0]=nbins_in;
352 59 : std::vector<Value*> value_pntr(1);
353 177 : value_pntr[0]= new Value(NULL,"arg",false);
354 131 : if(arePeriodic() && !ignore_periodicity) {value_pntr[0]->setDomain(intervalMinStr(),intervalMaxStr());}
355 41 : else {value_pntr[0]->setNotPeriodic();}
356 177 : Grid args_grid = Grid("grid",value_pntr,min,max,nbins,false,false);
357 :
358 59 : std::vector<double> args(args_grid.getSize(),0.0);
359 56341 : for(unsigned int i=0; i<args.size(); i++) {
360 56223 : args[i] = args_grid.getPoint(i)[0];
361 : }
362 : std::vector<double> argsT;
363 59 : std::vector<std::vector<double> > values;
364 59 : std::vector<std::vector<double> > derivs;
365 :
366 354 : ofile_values.addConstantField("bf_keywords").printField("bf_keywords","{"+getKeywordString()+"}");
367 354 : ofile_derivs.addConstantField("bf_keywords").printField("bf_keywords","{"+getKeywordString()+"}");
368 :
369 295 : ofile_values.addConstantField("min").printField("min",intervalMinStr());
370 295 : ofile_values.addConstantField("max").printField("max",intervalMaxStr());
371 :
372 295 : ofile_derivs.addConstantField("min").printField("min",intervalMinStr());
373 295 : ofile_derivs.addConstantField("max").printField("max",intervalMaxStr());
374 :
375 295 : ofile_values.addConstantField("nbins").printField("nbins",static_cast<int>(args_grid.getNbin()[0]));
376 295 : ofile_derivs.addConstantField("nbins").printField("nbins",static_cast<int>(args_grid.getNbin()[0]));
377 :
378 59 : if(arePeriodic()) {
379 84 : ofile_values.addConstantField("periodic").printField("periodic","true");
380 84 : ofile_derivs.addConstantField("periodic").printField("periodic","true");
381 : }
382 : else {
383 152 : ofile_values.addConstantField("periodic").printField("periodic","false");
384 152 : ofile_derivs.addConstantField("periodic").printField("periodic","false");
385 : }
386 :
387 59 : getMultipleValue(args,argsT,values,derivs,numerical_deriv);
388 59 : ofile_values.fmtField(output_fmt);
389 59 : ofile_derivs.fmtField(output_fmt);
390 56341 : for(unsigned int i=0; i<args.size(); i++) {
391 37482 : ofile_values.printField("arg",args[i]);
392 37482 : ofile_derivs.printField("arg",args[i]);
393 528451 : for(unsigned int k=0; k<getNumberOfBasisFunctions(); k++) {
394 1019420 : ofile_values.printField(getBasisFunctionLabel(k),values[i][k]);
395 1019420 : ofile_derivs.printField("d_"+getBasisFunctionLabel(k),derivs[i][k]);
396 : }
397 18741 : ofile_values.printField();
398 18741 : ofile_derivs.printField();
399 : }
400 59 : ofile_values.fmtField();
401 59 : ofile_derivs.fmtField();
402 :
403 59 : delete value_pntr[0]; value_pntr.clear();
404 :
405 59 : }
406 :
407 :
408 : }
409 : }
|
