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 542 : keys.add("compulsory","ORDER","The order of the basis function expansion.");
39 542 : keys.add("compulsory","MINIMUM","The minimum of the interval on which the basis functions are defined.");
40 542 : keys.add("compulsory","MAXIMUM","The maximum of the interval on which the basis functions are defined.");
41 542 : keys.add("hidden","NGRID_POINTS","The number of grid points used for numerical integrals");
42 542 : keys.addFlag("DEBUG_INFO",false,"Print out more detailed information about the basis set. Useful for debugging.");
43 542 : 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 : {
78 249 : bf_keywords_.push_back(getName());
79 498 : if(keywords.exists("ORDER")) {
80 711 : parse("ORDER",norder_); addKeywordToList("ORDER",norder_);
81 : }
82 249 : nbasis_=norder_+1;
83 : //
84 : std::string str_imin; std::string str_imax;
85 496 : if(keywords.exists("MINIMUM") && keywords.exists("MAXIMUM")) {
86 741 : parse("MINIMUM",str_imin); addKeywordToList("MINIMUM",str_imin);
87 741 : 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 249 : if(!Tools::convertNoexcept(str_imin,interval_min_)) {
96 0 : plumed_merror(getName()+": cannot convert the value given in MINIMUM to a double");
97 : }
98 249 : if(!Tools::convertNoexcept(str_imax,interval_max_)) {
99 0 : plumed_merror(getName()+": cannot convert the value given in MAXIMUM to a double");
100 : }
101 249 : if(interval_min_>interval_max_) {plumed_merror(getName()+": MINIMUM and MAXIMUM are not correctly defined");}
102 : //
103 249 : parseFlag("DEBUG_INFO",print_debug_info_);
104 498 : if(keywords.exists("NUMERICAL_INTEGRALS")) {
105 336 : parseFlag("NUMERICAL_INTEGRALS",numerical_uniform_integrals_);
106 : }
107 498 : if(keywords.exists("NGRID_POINTS")) {
108 494 : parse("NGRID_POINTS",nbins_);
109 : }
110 : // log.printf(" %s \n",getKeywordString().c_str());
111 :
112 249 : }
113 :
114 :
115 79 : void BasisFunctions::setIntrinsicInterval(const double interval_intrinsic_min_in, const double interval_intrinsic_max_in) {
116 79 : interval_intrinsic_min_ = interval_intrinsic_min_in;
117 79 : interval_intrinsic_max_ = interval_intrinsic_max_in;
118 79 : VesTools::convertDbl2Str(interval_intrinsic_min_,interval_intrinsic_min_str_);
119 79 : VesTools::convertDbl2Str(interval_intrinsic_max_,interval_intrinsic_max_str_);
120 79 : plumed_massert(interval_intrinsic_min_<interval_intrinsic_max_,"setIntrinsicInterval: intrinsic intervals are not defined correctly");
121 79 : }
122 :
123 :
124 170 : void BasisFunctions::setIntrinsicInterval(const std::string& interval_intrinsic_min_str_in, const std::string& interval_intrinsic_max_str_in) {
125 170 : interval_intrinsic_min_str_ = interval_intrinsic_min_str_in;
126 170 : interval_intrinsic_max_str_ = interval_intrinsic_max_str_in;
127 170 : if(!Tools::convertNoexcept(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 170 : if(!Tools::convertNoexcept(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 170 : plumed_massert(interval_intrinsic_min_<interval_intrinsic_max_,"setIntrinsicInterval: intrinsic intervals are not defined correctly");
134 170 : }
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::convertNoexcept(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::convertNoexcept(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 249 : void BasisFunctions::setupInterval() {
160 : // if(!intervalBounded()){plumed_merror("setupInterval() only works for bounded interval");}
161 249 : interval_intrinsic_range_ = interval_intrinsic_max_-interval_intrinsic_min_;
162 249 : interval_intrinsic_mean_ = 0.5*(interval_intrinsic_max_+interval_intrinsic_min_);
163 249 : interval_range_ = interval_max_-interval_min_;
164 249 : interval_mean_ = 0.5*(interval_max_+interval_min_);
165 249 : argT_derivf_ = interval_intrinsic_range_/interval_range_;
166 249 : }
167 :
168 :
169 78 : void BasisFunctions::setupLabels() {
170 884 : for(unsigned int i=0; i < nbasis_; i++) {
171 806 : std::string is; Tools::convert(i,is);
172 1612 : bf_labels_[i]=bf_label_prefix_+is+"(s)";
173 : }
174 78 : }
175 :
176 :
177 79 : void BasisFunctions::setupUniformIntegrals() {
178 79 : numerical_uniform_integrals_=true;
179 79 : numericalUniformIntegrals();
180 79 : }
181 :
182 :
183 249 : void BasisFunctions::setupBF() {
184 249 : if(interval_intrinsic_min_>interval_intrinsic_max_) {plumed_merror("setupBF: default intervals are not correctly set");}
185 249 : setupInterval();
186 249 : setupLabels();
187 249 : if(bf_labels_.size()==1) {plumed_merror("setupBF: the labels of the basis functions are not correct.");}
188 249 : if(!numerical_uniform_integrals_) {setupUniformIntegrals();}
189 6 : else {numericalUniformIntegrals();}
190 249 : if(uniform_integrals_.size()==1) {plumed_merror("setupBF: the integrals of the basis functions is not correct.");}
191 249 : if(type_=="Undefined") {plumed_merror("setupBF: the type of the basis function is not defined.");}
192 249 : if(description_=="Undefined") {plumed_merror("setupBF: the description of the basis function is not defined.");}
193 249 : has_been_set=true;
194 249 : printInfo();
195 249 : }
196 :
197 :
198 249 : void BasisFunctions::printInfo() const {
199 249 : if(!has_been_set) {plumed_merror("the basis set has not be setup correctly");}
200 249 : log.printf(" One-dimensional basis set\n");
201 249 : log.printf(" Description: %s\n",description_.c_str());
202 249 : log.printf(" Type: %s\n",type_.c_str());
203 249 : if(periodic_) {log.printf(" The basis functions are periodic\n");}
204 249 : log.printf(" Order of basis set: %u\n",norder_);
205 249 : 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 249 : log.printf(" Interval of basis set: %s to %s\n",interval_min_str_.c_str(),interval_max_str_.c_str());
208 249 : log.printf(" Description of basis functions:\n");
209 4293 : for(unsigned int i=0; i < nbasis_; i++) {log.printf(" %2u %10s\n",i,bf_labels_[i].c_str());}
210 : //
211 249 : if(print_debug_info_) {
212 38 : log.printf(" Debug information:\n");
213 : // log.printf(" Default interval of basis set: [%f,%f]\n",interval_intrinsic_min_,interval_intrinsic_max_);
214 38 : log.printf(" Intrinsic interval of basis set: [%s,%s]\n",interval_intrinsic_min_str_.c_str(),interval_intrinsic_max_str_.c_str());
215 38 : 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 38 : log.printf(" Defined interval of basis set: [%s,%s]\n",interval_min_str_.c_str(),interval_max_str_.c_str());
218 38 : log.printf(" Defined interval of basis set: range=%f, mean=%f\n",interval_range_,interval_mean_);
219 38 : log.printf(" Derivative factor due to interval translation: %f\n",argT_derivf_);
220 38 : log.printf(" Integral of basis functions over the interval:\n");
221 38 : if(numerical_uniform_integrals_) {log.printf(" Note: calculated numerically\n");}
222 558 : for(unsigned int i=0; i < nbasis_; i++) {log.printf(" %2u %16.10f\n",i,uniform_integrals_[i]);}
223 38 : log.printf(" --------------------------\n");
224 : }
225 249 : }
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 85 : void BasisFunctions::numericalUniformIntegrals() {
240 85 : std::vector<std::string> grid_min(1); grid_min[0]=intervalMinStr();
241 85 : std::vector<std::string> grid_max(1); grid_max[0]=intervalMaxStr();
242 85 : std::vector<unsigned int> grid_bins(1); grid_bins[0]=nbins_;
243 170 : std::vector<std::unique_ptr<Value>> arguments(1); arguments[0]= Tools::make_unique<Value>(nullptr,"arg",false);
244 105 : if(arePeriodic()) {arguments[0]->setDomain(intervalMinStr(),intervalMaxStr());}
245 75 : else {arguments[0]->setNotPeriodic();}
246 170 : auto uniform_grid = Tools::make_unique<Grid>("uniform",Tools::unique2raw(arguments),grid_min,grid_max,grid_bins,false,false);
247 : //
248 85 : double inverse_normalization = 1.0/(intervalMax()-intervalMin());
249 85245 : for(Grid::index_t l=0; l<uniform_grid->getSize(); l++) {
250 85160 : uniform_grid->setValue(l,inverse_normalization);
251 : }
252 85 : uniform_integrals_ = numericalTargetDistributionIntegralsFromGrid(uniform_grid.get());
253 170 : }
254 :
255 :
256 261 : std::vector<double> BasisFunctions::numericalTargetDistributionIntegralsFromGrid(const Grid* grid_pntr) const {
257 261 : plumed_massert(grid_pntr!=NULL,"the grid is not defined");
258 261 : plumed_massert(grid_pntr->getDimension()==1,"the target distribution grid should be one dimensional");
259 : //
260 261 : std::vector<double> targetdist_integrals(nbasis_,0.0);
261 522 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(grid_pntr);
262 :
263 138336 : for(Grid::index_t k=0; k < grid_pntr->getSize(); k++) {
264 138075 : double arg = grid_pntr->getPoint(k)[0];
265 138075 : std::vector<double> bf_values(nbasis_);
266 138075 : std::vector<double> bf_derivs(nbasis_);
267 138075 : bool inside=true;
268 138075 : double argT=0.0;
269 138075 : getAllValues(arg,argT,inside,bf_values,bf_derivs);
270 3336092 : for(unsigned int i=0; i < nbasis_; i++) {
271 3198017 : targetdist_integrals[i] += (integration_weights[k] * grid_pntr->getValue(k)) * bf_values[i];
272 : }
273 : }
274 : // assume that the first function is the constant
275 261 : bool inside=true;
276 261 : double argT=0.0;
277 261 : targetdist_integrals[0] = getValue(0.0,0,argT,inside);
278 261 : return targetdist_integrals;
279 : }
280 :
281 :
282 247 : std::vector<double> BasisFunctions::getTargetDistributionIntegrals(const TargetDistribution* targetdist_pntr) const {
283 247 : if(targetdist_pntr==NULL) {
284 : return getUniformIntegrals();
285 : }
286 : else {
287 : Grid* targetdist_grid = targetdist_pntr->getTargetDistGridPntr();
288 176 : return numericalTargetDistributionIntegralsFromGrid(targetdist_grid);
289 : }
290 : }
291 :
292 :
293 142 : std::string BasisFunctions::getKeywordString() const {
294 142 : std::string str_keywords=bf_keywords_[0];
295 762 : for(unsigned int i=1; i<bf_keywords_.size(); i++) {str_keywords+=" "+bf_keywords_[i];}
296 142 : return str_keywords;
297 : }
298 :
299 :
300 261 : double BasisFunctions::getValue(const double arg, const unsigned int n, double& argT, bool& inside_range) const {
301 261 : plumed_massert(n<numberOfBasisFunctions(),"getValue: n is outside range of the defined order of the basis set");
302 261 : inside_range=true;
303 261 : std::vector<double> tmp_values(numberOfBasisFunctions());
304 261 : std::vector<double> tmp_derivs(numberOfBasisFunctions());
305 261 : getAllValues(arg, argT, inside_range, tmp_values, tmp_derivs);
306 522 : return tmp_values[n];
307 : }
308 :
309 :
310 0 : void BasisFunctions::getAllValuesNumericalDerivs(const double arg, double& argT, bool& inside_range, std::vector<double>& values, std::vector<double>& derivs) const {
311 : // use forward difference, unless very close to the boundary
312 : double delta = sqrt(epsilon);
313 0 : if((arg+delta)>intervalMax()) {
314 : delta *= -1.0;
315 : }
316 0 : inside_range=true;
317 0 : std::vector<double> values_delta(numberOfBasisFunctions());
318 0 : std::vector<double> derivs_dummy(numberOfBasisFunctions());
319 0 : getAllValues(arg+delta, argT, inside_range, values_delta, derivs_dummy);
320 0 : getAllValues(arg, argT, inside_range, values, derivs_dummy);
321 0 : for(unsigned int i=0; i<numberOfBasisFunctions(); i++) {
322 0 : derivs[i] = (values_delta[i]-values[i])/delta;
323 : }
324 0 : }
325 :
326 :
327 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 {
328 71 : argsT.resize(args.size());
329 : values.clear();
330 : derivs.clear();
331 23019 : for(unsigned int i=0; i<args.size(); i++) {
332 22948 : std::vector<double> tmp_values(getNumberOfBasisFunctions());
333 22948 : std::vector<double> tmp_derivs(getNumberOfBasisFunctions());
334 22948 : bool inside_interval=true;
335 22948 : if(!numerical_deriv) {
336 22948 : getAllValues(args[i],argsT[i],inside_interval,tmp_values,tmp_derivs);
337 : } else {
338 0 : getAllValuesNumericalDerivs(args[i],argsT[i],inside_interval,tmp_values,tmp_derivs);
339 : }
340 22948 : values.push_back(tmp_values);
341 22948 : derivs.push_back(tmp_derivs);
342 : }
343 71 : }
344 :
345 :
346 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 {
347 71 : std::vector<std::string> min(1); min[0]=min_in;
348 71 : std::vector<std::string> max(1); max[0]=max_in;
349 71 : std::vector<unsigned int> nbins(1); nbins[0]=nbins_in;
350 71 : std::vector<std::unique_ptr<Value>> value_pntr(1);
351 142 : value_pntr[0]= Tools::make_unique<Value>(nullptr,"arg",false);
352 117 : if(arePeriodic() && !ignore_periodicity) {value_pntr[0]->setDomain(intervalMinStr(),intervalMaxStr());}
353 48 : else {value_pntr[0]->setNotPeriodic();}
354 142 : Grid args_grid = Grid("grid",Tools::unique2raw(value_pntr),min,max,nbins,false,false);
355 :
356 71 : std::vector<double> args(args_grid.getSize(),0.0);
357 23019 : for(unsigned int i=0; i<args.size(); i++) {
358 22948 : args[i] = args_grid.getPoint(i)[0];
359 : }
360 : std::vector<double> argsT;
361 : std::vector<std::vector<double> > values;
362 : std::vector<std::vector<double> > derivs;
363 :
364 142 : ofile_values.addConstantField("bf_keywords").printField("bf_keywords","{"+getKeywordString()+"}");
365 142 : ofile_derivs.addConstantField("bf_keywords").printField("bf_keywords","{"+getKeywordString()+"}");
366 :
367 213 : ofile_values.addConstantField("min").printField("min",intervalMinStr());
368 213 : ofile_values.addConstantField("max").printField("max",intervalMaxStr());
369 :
370 213 : ofile_derivs.addConstantField("min").printField("min",intervalMinStr());
371 213 : ofile_derivs.addConstantField("max").printField("max",intervalMaxStr());
372 :
373 142 : ofile_values.addConstantField("nbins").printField("nbins",static_cast<int>(args_grid.getNbin()[0]));
374 142 : ofile_derivs.addConstantField("nbins").printField("nbins",static_cast<int>(args_grid.getNbin()[0]));
375 :
376 71 : if(arePeriodic()) {
377 52 : ofile_values.addConstantField("periodic").printField("periodic","true");
378 52 : ofile_derivs.addConstantField("periodic").printField("periodic","true");
379 : }
380 : else {
381 90 : ofile_values.addConstantField("periodic").printField("periodic","false");
382 90 : ofile_derivs.addConstantField("periodic").printField("periodic","false");
383 : }
384 :
385 71 : getMultipleValue(args,argsT,values,derivs,numerical_deriv);
386 71 : ofile_values.fmtField(output_fmt_values);
387 71 : ofile_derivs.fmtField(output_fmt_derivs);
388 23019 : for(unsigned int i=0; i<args.size(); i++) {
389 45896 : ofile_values.printField("arg",args[i]);
390 22948 : ofile_derivs.printField("arg",args[i]);
391 410060 : for(unsigned int k=0; k<getNumberOfBasisFunctions(); k++) {
392 774224 : ofile_values.printField(getBasisFunctionLabel(k),values[i][k]);
393 774224 : ofile_derivs.printField("d_"+getBasisFunctionLabel(k),derivs[i][k]);
394 : }
395 22948 : ofile_values.printField();
396 22948 : ofile_derivs.printField();
397 : }
398 71 : ofile_values.fmtField();
399 71 : ofile_derivs.fmtField();
400 :
401 213 : }
402 :
403 :
404 : }
405 : }
|
