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