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