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 "LinearBasisSetExpansion.h"
24 : #include "VesBias.h"
25 : #include "CoeffsVector.h"
26 : #include "VesTools.h"
27 : #include "GridIntegrationWeights.h"
28 : #include "BasisFunctions.h"
29 : #include "TargetDistribution.h"
30 :
31 :
32 : #include "tools/Keywords.h"
33 : #include "tools/Grid.h"
34 : #include "tools/Communicator.h"
35 :
36 : #include "GridProjWeights.h"
37 :
38 : #include <limits>
39 :
40 : namespace PLMD {
41 : namespace ves {
42 :
43 0 : void LinearBasisSetExpansion::registerKeywords(Keywords& keys) {
44 0 : }
45 :
46 :
47 131 : LinearBasisSetExpansion::LinearBasisSetExpansion(
48 : const std::string& label,
49 : const double beta_in,
50 : Communicator& cc,
51 : const std::vector<Value*>& args_pntrs_in,
52 : std::vector<BasisFunctions*>& basisf_pntrs_in,
53 131 : CoeffsVector* bias_coeffs_pntr_in):
54 131 : label_(label),
55 131 : action_pntr_(NULL),
56 131 : vesbias_pntr_(NULL),
57 131 : mycomm_(cc),
58 131 : serial_(false),
59 131 : beta_(beta_in),
60 131 : kbt_(1.0/beta_),
61 131 : args_pntrs_(args_pntrs_in),
62 131 : nargs_(args_pntrs_.size()),
63 131 : basisf_pntrs_(basisf_pntrs_in),
64 131 : nbasisf_(basisf_pntrs_.size()),
65 131 : bias_coeffs_pntr_(bias_coeffs_pntr_in),
66 131 : ncoeffs_(0),
67 131 : grid_min_(nargs_),
68 131 : grid_max_(nargs_),
69 131 : grid_bins_(nargs_,100),
70 131 : targetdist_grid_label_("targetdist"),
71 131 : step_of_last_biasgrid_update(-1000),
72 131 : step_of_last_biaswithoutcutoffgrid_update(-1000),
73 131 : step_of_last_fesgrid_update(-1000),
74 131 : log_targetdist_grid_pntr_(NULL),
75 131 : targetdist_grid_pntr_(NULL),
76 131 : targetdist_pntr_(NULL),
77 131 : grid_fmt_("%14.9f") {
78 131 : plumed_massert(args_pntrs_.size()==basisf_pntrs_.size(),"number of arguments and basis functions do not match");
79 297 : for(unsigned int k=0; k<nargs_; k++) {
80 166 : nbasisf_[k]=basisf_pntrs_[k]->getNumberOfBasisFunctions();
81 : }
82 : //
83 131 : if(bias_coeffs_pntr_==NULL) {
84 0 : bias_coeffs_pntr_ = new CoeffsVector(label_+".coeffs",args_pntrs_,basisf_pntrs_,mycomm_,true);
85 : }
86 131 : plumed_massert(bias_coeffs_pntr_->numberOfDimensions()==basisf_pntrs_.size(),"dimension of coeffs does not match with number of basis functions ");
87 : //
88 131 : ncoeffs_ = bias_coeffs_pntr_->numberOfCoeffs();
89 131 : targetdist_averages_pntr_ = Tools::make_unique<CoeffsVector>(*bias_coeffs_pntr_);
90 :
91 131 : std::string targetdist_averages_label = bias_coeffs_pntr_->getLabel();
92 131 : if(targetdist_averages_label.find("coeffs")!=std::string::npos) {
93 131 : targetdist_averages_label.replace(targetdist_averages_label.find("coeffs"), std::string("coeffs").length(), "targetdist_averages");
94 : } else {
95 : targetdist_averages_label += "_targetdist_averages";
96 : }
97 131 : targetdist_averages_pntr_->setLabels(targetdist_averages_label);
98 : //
99 297 : for(unsigned int k=0; k<nargs_; k++) {
100 332 : grid_min_[k] = basisf_pntrs_[k]->intervalMinStr();
101 332 : grid_max_[k] = basisf_pntrs_[k]->intervalMaxStr();
102 : }
103 393 : }
104 :
105 131 : LinearBasisSetExpansion::~LinearBasisSetExpansion() {
106 393 : }
107 :
108 :
109 0 : bool LinearBasisSetExpansion::isStaticTargetDistFileOutputActive() const {
110 : bool output_static_targetdist_files=true;
111 0 : if(vesbias_pntr_!=NULL) {
112 : output_static_targetdist_files = vesbias_pntr_->isStaticTargetDistFileOutputActive();
113 : }
114 0 : return output_static_targetdist_files;
115 : }
116 :
117 :
118 91 : void LinearBasisSetExpansion::linkVesBias(VesBias* vesbias_pntr_in) {
119 91 : vesbias_pntr_ = vesbias_pntr_in;
120 91 : action_pntr_ = static_cast<Action*>(vesbias_pntr_in);
121 :
122 91 : }
123 :
124 :
125 0 : void LinearBasisSetExpansion::linkAction(Action* action_pntr_in) {
126 0 : action_pntr_ = action_pntr_in;
127 0 : }
128 :
129 :
130 130 : void LinearBasisSetExpansion::setGridBins(const std::vector<unsigned int>& grid_bins_in) {
131 130 : plumed_massert(grid_bins_in.size()==nargs_,"the number of grid bins given doesn't match the number of arguments");
132 130 : plumed_massert(!bias_grid_pntr_,"setGridBins should be used before setting up the grids, otherwise it doesn't work");
133 130 : plumed_massert(!fes_grid_pntr_,"setGridBins should be used before setting up the grids, otherwise it doesn't work");
134 130 : grid_bins_=grid_bins_in;
135 130 : }
136 :
137 :
138 39 : void LinearBasisSetExpansion::setGridBins(const unsigned int nbins) {
139 39 : std::vector<unsigned int> grid_bins_in(nargs_,nbins);
140 39 : setGridBins(grid_bins_in);
141 39 : }
142 :
143 :
144 221 : std::unique_ptr<Grid> LinearBasisSetExpansion::setupGeneralGrid(const std::string& label_suffix, const bool usederiv) {
145 221 : bool use_spline = false;
146 221 : auto grid_pntr = Tools::make_unique<Grid>(label_+"."+label_suffix,args_pntrs_,grid_min_,grid_max_,grid_bins_,use_spline,usederiv);
147 221 : return grid_pntr;
148 : }
149 :
150 :
151 167 : void LinearBasisSetExpansion::setupBiasGrid(const bool usederiv) {
152 167 : if(bias_grid_pntr_) {
153 : return;
154 : }
155 260 : bias_grid_pntr_ = setupGeneralGrid("bias",usederiv);
156 130 : if(biasCutoffActive()) {
157 9 : bias_withoutcutoff_grid_pntr_ = setupGeneralGrid("bias_withoutcutoff",usederiv);
158 : }
159 : }
160 :
161 :
162 120 : void LinearBasisSetExpansion::setupFesGrid() {
163 120 : if(fes_grid_pntr_) {
164 : return;
165 : }
166 88 : if(!bias_grid_pntr_) {
167 37 : setupBiasGrid(true);
168 : }
169 264 : fes_grid_pntr_ = setupGeneralGrid("fes",false);
170 : }
171 :
172 :
173 8 : void LinearBasisSetExpansion::setupFesProjGrid() {
174 8 : if(!fes_grid_pntr_) {
175 1 : setupFesGrid();
176 : }
177 8 : }
178 :
179 :
180 830 : void LinearBasisSetExpansion::updateBiasGrid() {
181 830 : plumed_massert(bias_grid_pntr_,"the bias grid is not defined");
182 830 : if(action_pntr_!=NULL && getStepOfLastBiasGridUpdate()==action_pntr_->getStep()) {
183 : return;
184 : }
185 1982342 : for(Grid::index_t l=0; l<bias_grid_pntr_->getSize(); l++) {
186 1981749 : std::vector<double> forces(nargs_);
187 1981749 : std::vector<double> args = bias_grid_pntr_->getPoint(l);
188 1981749 : bool all_inside=true;
189 1981749 : double bias=getBiasAndForces(args,all_inside,forces);
190 : //
191 1981749 : if(biasCutoffActive()) {
192 600 : vesbias_pntr_->applyBiasCutoff(bias,forces);
193 : }
194 : //
195 1981749 : if(bias_grid_pntr_->hasDerivatives()) {
196 1474032 : bias_grid_pntr_->setValueAndDerivatives(l,bias,forces);
197 : } else {
198 507717 : bias_grid_pntr_->setValue(l,bias);
199 : }
200 : //
201 : }
202 593 : if(vesbias_pntr_!=NULL) {
203 476 : vesbias_pntr_->setCurrentBiasMaxValue(bias_grid_pntr_->getMaxValue());
204 : }
205 593 : if(action_pntr_!=NULL) {
206 476 : setStepOfLastBiasGridUpdate(action_pntr_->getStep());
207 : }
208 : }
209 :
210 :
211 27 : void LinearBasisSetExpansion::updateBiasWithoutCutoffGrid() {
212 27 : plumed_massert(bias_withoutcutoff_grid_pntr_,"the bias without cutoff grid is not defined");
213 27 : plumed_massert(biasCutoffActive(),"the bias cutoff has to be active");
214 27 : plumed_massert(vesbias_pntr_!=NULL,"has to be linked to a VesBias to work");
215 27 : if(action_pntr_!=NULL && getStepOfLastBiasWithoutCutoffGridUpdate()==action_pntr_->getStep()) {
216 : return;
217 : }
218 : //
219 2423 : for(Grid::index_t l=0; l<bias_withoutcutoff_grid_pntr_->getSize(); l++) {
220 2400 : std::vector<double> forces(nargs_);
221 2400 : std::vector<double> args = bias_withoutcutoff_grid_pntr_->getPoint(l);
222 2400 : bool all_inside=true;
223 2400 : double bias=getBiasAndForces(args,all_inside,forces);
224 2400 : if(bias_withoutcutoff_grid_pntr_->hasDerivatives()) {
225 2400 : bias_withoutcutoff_grid_pntr_->setValueAndDerivatives(l,bias,forces);
226 : } else {
227 0 : bias_withoutcutoff_grid_pntr_->setValue(l,bias);
228 : }
229 : }
230 : //
231 23 : double bias_max = bias_withoutcutoff_grid_pntr_->getMaxValue();
232 23 : double bias_min = bias_withoutcutoff_grid_pntr_->getMinValue();
233 : double shift = 0.0;
234 : bool bias_shifted=false;
235 23 : if(bias_min < 0.0) {
236 22 : shift += -bias_min;
237 : bias_shifted=true;
238 22 : BiasCoeffs()[0] -= bias_min;
239 22 : bias_max -= bias_min;
240 : }
241 23 : if(bias_max > vesbias_pntr_->getBiasCutoffValue()) {
242 22 : shift += -(bias_max-vesbias_pntr_->getBiasCutoffValue());
243 : bias_shifted=true;
244 22 : BiasCoeffs()[0] -= (bias_max-vesbias_pntr_->getBiasCutoffValue());
245 22 : bias_max -= (bias_max-vesbias_pntr_->getBiasCutoffValue());
246 : }
247 23 : if(bias_shifted) {
248 : // this should be done inside a grid function really,
249 : // need to define my grid class for that
250 2322 : for(Grid::index_t l=0; l<bias_withoutcutoff_grid_pntr_->getSize(); l++) {
251 2300 : if(bias_withoutcutoff_grid_pntr_->hasDerivatives()) {
252 2300 : std::vector<double> zeros(nargs_,0.0);
253 2300 : bias_withoutcutoff_grid_pntr_->addValueAndDerivatives(l,shift,zeros);
254 : } else {
255 0 : bias_withoutcutoff_grid_pntr_->addValue(l,shift);
256 : }
257 : }
258 : }
259 23 : if(vesbias_pntr_!=NULL) {
260 : vesbias_pntr_->setCurrentBiasMaxValue(bias_max);
261 : }
262 23 : if(action_pntr_!=NULL) {
263 23 : setStepOfLastBiasWithoutCutoffGridUpdate(action_pntr_->getStep());
264 : }
265 : }
266 :
267 :
268 539 : void LinearBasisSetExpansion::updateFesGrid() {
269 539 : plumed_massert(fes_grid_pntr_,"the FES grid is not defined");
270 539 : updateBiasGrid();
271 539 : if(action_pntr_!=NULL && getStepOfLastFesGridUpdate() == action_pntr_->getStep()) {
272 : return;
273 : }
274 : //
275 : double bias2fes_scalingf = -1.0;
276 1474154 : for(Grid::index_t l=0; l<fes_grid_pntr_->getSize(); l++) {
277 1473681 : double fes_value = bias2fes_scalingf*bias_grid_pntr_->getValue(l);
278 1473681 : if(log_targetdist_grid_pntr_!=NULL) {
279 1224051 : fes_value += kBT()*log_targetdist_grid_pntr_->getValue(l);
280 : }
281 1473681 : fes_grid_pntr_->setValue(l,fes_value);
282 : }
283 473 : fes_grid_pntr_->setMinToZero();
284 473 : if(action_pntr_!=NULL) {
285 473 : setStepOfLastFesGridUpdate(action_pntr_->getStep());
286 : }
287 : }
288 :
289 :
290 213 : void LinearBasisSetExpansion::writeBiasGridToFile(OFile& ofile, const bool append_file) const {
291 213 : plumed_massert(bias_grid_pntr_,"the bias grid is not defined");
292 213 : if(append_file) {
293 0 : ofile.enforceRestart();
294 : }
295 213 : bias_grid_pntr_->setOutputFmt(grid_fmt_);
296 213 : bias_grid_pntr_->writeToFile(ofile);
297 213 : }
298 :
299 :
300 5 : void LinearBasisSetExpansion::writeBiasWithoutCutoffGridToFile(OFile& ofile, const bool append_file) const {
301 5 : plumed_massert(bias_withoutcutoff_grid_pntr_,"the bias without cutoff grid is not defined");
302 5 : if(append_file) {
303 0 : ofile.enforceRestart();
304 : }
305 5 : bias_withoutcutoff_grid_pntr_->setOutputFmt(grid_fmt_);
306 5 : bias_withoutcutoff_grid_pntr_->writeToFile(ofile);
307 5 : }
308 :
309 :
310 169 : void LinearBasisSetExpansion::writeFesGridToFile(OFile& ofile, const bool append_file) const {
311 169 : plumed_massert(fes_grid_pntr_!=NULL,"the FES grid is not defined");
312 169 : if(append_file) {
313 0 : ofile.enforceRestart();
314 : }
315 169 : fes_grid_pntr_->setOutputFmt(grid_fmt_);
316 169 : fes_grid_pntr_->writeToFile(ofile);
317 169 : }
318 :
319 :
320 36 : void LinearBasisSetExpansion::writeFesProjGridToFile(const std::vector<std::string>& proj_arg, OFile& ofile, const bool append_file) const {
321 36 : plumed_massert(fes_grid_pntr_,"the FES grid is not defined");
322 36 : auto Fw = Tools::make_unique<FesWeight>(beta_);
323 36 : Grid proj_grid = fes_grid_pntr_->project(proj_arg,Fw.get());
324 36 : proj_grid.setMinToZero();
325 36 : if(append_file) {
326 0 : ofile.enforceRestart();
327 : }
328 36 : proj_grid.setOutputFmt(grid_fmt_);
329 36 : proj_grid.writeToFile(ofile);
330 36 : }
331 :
332 :
333 82 : void LinearBasisSetExpansion::writeTargetDistGridToFile(OFile& ofile, const bool append_file) const {
334 82 : if(targetdist_grid_pntr_==NULL) {
335 : return;
336 : }
337 82 : if(append_file) {
338 0 : ofile.enforceRestart();
339 : }
340 82 : targetdist_grid_pntr_->setOutputFmt(grid_fmt_);
341 82 : targetdist_grid_pntr_->writeToFile(ofile);
342 : }
343 :
344 :
345 82 : void LinearBasisSetExpansion::writeLogTargetDistGridToFile(OFile& ofile, const bool append_file) const {
346 82 : if(log_targetdist_grid_pntr_==NULL) {
347 : return;
348 : }
349 82 : if(append_file) {
350 0 : ofile.enforceRestart();
351 : }
352 82 : log_targetdist_grid_pntr_->setOutputFmt(grid_fmt_);
353 82 : log_targetdist_grid_pntr_->writeToFile(ofile);
354 : }
355 :
356 :
357 20 : void LinearBasisSetExpansion::writeTargetDistProjGridToFile(const std::vector<std::string>& proj_arg, OFile& ofile, const bool append_file) const {
358 20 : if(targetdist_grid_pntr_==NULL) {
359 0 : return;
360 : }
361 20 : if(append_file) {
362 0 : ofile.enforceRestart();
363 : }
364 20 : Grid proj_grid = TargetDistribution::getMarginalDistributionGrid(targetdist_grid_pntr_,proj_arg);
365 20 : proj_grid.setOutputFmt(grid_fmt_);
366 20 : proj_grid.writeToFile(ofile);
367 20 : }
368 :
369 :
370 0 : void LinearBasisSetExpansion::writeTargetDistributionToFile(const std::string& filename) const {
371 0 : OFile of1;
372 0 : OFile of2;
373 0 : if(action_pntr_!=NULL) {
374 0 : of1.link(*action_pntr_);
375 0 : of2.link(*action_pntr_);
376 : }
377 0 : of1.enforceBackup();
378 0 : of2.enforceBackup();
379 0 : of1.open(filename);
380 0 : of2.open(FileBase::appendSuffix(filename,".log"));
381 0 : writeTargetDistGridToFile(of1);
382 0 : writeLogTargetDistGridToFile(of2);
383 0 : of1.close();
384 0 : of2.close();
385 0 : }
386 :
387 :
388 2011898 : double LinearBasisSetExpansion::getBiasAndForces(const std::vector<double>& args_values, bool& all_inside, std::vector<double>& forces, std::vector<double>& coeffsderivs_values, std::vector<BasisFunctions*>& basisf_pntrs_in, CoeffsVector* coeffs_pntr_in, Communicator* comm_in) {
389 2011898 : unsigned int nargs = args_values.size();
390 2011898 : plumed_assert(coeffs_pntr_in->numberOfDimensions()==nargs);
391 2011898 : plumed_assert(basisf_pntrs_in.size()==nargs);
392 2011898 : plumed_assert(forces.size()==nargs);
393 2011898 : plumed_assert(coeffsderivs_values.size()==coeffs_pntr_in->numberOfCoeffs());
394 :
395 2011898 : std::vector<double> args_values_trsfrm(nargs);
396 2011898 : all_inside = true;
397 : //
398 2011898 : std::vector< std::vector <double> > bf_values(nargs);
399 2011898 : std::vector< std::vector <double> > bf_derivs(nargs);
400 : //
401 5966477 : for(unsigned int k=0; k<nargs; k++) {
402 3954579 : bf_values[k].assign(basisf_pntrs_in[k]->getNumberOfBasisFunctions(),0.0);
403 3954579 : bf_derivs[k].assign(basisf_pntrs_in[k]->getNumberOfBasisFunctions(),0.0);
404 3954579 : bool curr_inside=true;
405 3954579 : basisf_pntrs_in[k]->getAllValues(args_values[k],args_values_trsfrm[k],curr_inside,bf_values[k],bf_derivs[k]);
406 3954579 : if(!curr_inside) {
407 6932 : all_inside=false;
408 : }
409 3954579 : forces[k]=0.0;
410 : }
411 : //
412 : size_t stride=1;
413 : size_t rank=0;
414 2011898 : if(comm_in!=NULL) {
415 2011898 : stride=comm_in->Get_size();
416 2011898 : rank=comm_in->Get_rank();
417 : }
418 : // loop over coeffs
419 2011898 : double bias=0.0;
420 144770504 : for(size_t i=rank; i<coeffs_pntr_in->numberOfCoeffs(); i+=stride) {
421 142758606 : std::vector<unsigned int> indices=coeffs_pntr_in->getIndices(i);
422 : double coeff = coeffs_pntr_in->getValue(i);
423 : double bf_curr=1.0;
424 435045696 : for(unsigned int k=0; k<nargs; k++) {
425 292287090 : bf_curr*=bf_values[k][indices[k]];
426 : }
427 142758606 : bias+=coeff*bf_curr;
428 142758606 : coeffsderivs_values[i] = bf_curr;
429 435045696 : for(unsigned int k=0; k<nargs; k++) {
430 : double der = 1.0;
431 898593144 : for(unsigned int l=0; l<nargs; l++) {
432 606306054 : if(l!=k) {
433 314018964 : der*=bf_values[l][indices[l]];
434 : } else {
435 292287090 : der*=bf_derivs[l][indices[l]];
436 : }
437 : }
438 292287090 : forces[k]-=coeff*der;
439 : // maybe faster but dangerous
440 : // forces[k]-=coeff*bf_curr*(bf_derivs[k][indices[k]]/bf_values[k][indices[k]]);
441 : }
442 : }
443 : //
444 2011898 : if(comm_in!=NULL) {
445 : // coeffsderivs_values is not summed as the mpi Sum is done later on for the averages
446 2011898 : comm_in->Sum(bias);
447 2011898 : comm_in->Sum(forces);
448 : }
449 2011898 : return bias;
450 2011898 : }
451 :
452 :
453 862315 : void LinearBasisSetExpansion::getBasisSetValues(const std::vector<double>& args_values, std::vector<double>& basisset_values, std::vector<BasisFunctions*>& basisf_pntrs_in, CoeffsVector* coeffs_pntr_in, Communicator* comm_in) {
454 862315 : unsigned int nargs = args_values.size();
455 862315 : plumed_assert(coeffs_pntr_in->numberOfDimensions()==nargs);
456 862315 : plumed_assert(basisf_pntrs_in.size()==nargs);
457 :
458 862315 : std::vector<double> args_values_trsfrm(nargs);
459 862315 : std::vector< std::vector <double> > bf_values;
460 : //
461 2583312 : for(unsigned int k=0; k<nargs; k++) {
462 1720997 : std::vector<double> tmp_val(basisf_pntrs_in[k]->getNumberOfBasisFunctions());
463 1720997 : std::vector<double> tmp_der(tmp_val.size());
464 1720997 : bool inside=true;
465 1720997 : basisf_pntrs_in[k]->getAllValues(args_values[k],args_values_trsfrm[k],inside,tmp_val,tmp_der);
466 1720997 : bf_values.push_back(tmp_val);
467 : }
468 : //
469 : size_t stride=1;
470 : size_t rank=0;
471 862315 : if(comm_in!=NULL) {
472 0 : stride=comm_in->Get_size();
473 0 : rank=comm_in->Get_rank();
474 : }
475 : // loop over basis set
476 95955507 : for(size_t i=rank; i<coeffs_pntr_in->numberOfCoeffs(); i+=stride) {
477 95093192 : std::vector<unsigned int> indices=coeffs_pntr_in->getIndices(i);
478 : double bf_curr=1.0;
479 298507612 : for(unsigned int k=0; k<nargs; k++) {
480 203414420 : bf_curr*=bf_values[k][indices[k]];
481 : }
482 95093192 : basisset_values[i] = bf_curr;
483 : }
484 : //
485 862315 : if(comm_in!=NULL) {
486 0 : comm_in->Sum(basisset_values);
487 : }
488 1724630 : }
489 :
490 :
491 408 : double LinearBasisSetExpansion::getBasisSetValue(const std::vector<double>& args_values, const size_t index, std::vector<BasisFunctions*>& basisf_pntrs_in, CoeffsVector* coeffs_pntr_in) {
492 408 : unsigned int nargs = args_values.size();
493 408 : plumed_assert(coeffs_pntr_in->numberOfDimensions()==nargs);
494 408 : plumed_assert(basisf_pntrs_in.size()==nargs);
495 :
496 408 : std::vector<double> args_values_trsfrm(nargs);
497 408 : std::vector< std::vector <double> > bf_values;
498 : //
499 938 : for(unsigned int k=0; k<nargs; k++) {
500 530 : std::vector<double> tmp_val(basisf_pntrs_in[k]->getNumberOfBasisFunctions());
501 530 : std::vector<double> tmp_der(tmp_val.size());
502 530 : bool inside=true;
503 530 : basisf_pntrs_in[k]->getAllValues(args_values[k],args_values_trsfrm[k],inside,tmp_val,tmp_der);
504 530 : bf_values.push_back(tmp_val);
505 : }
506 : //
507 408 : std::vector<unsigned int> indices=coeffs_pntr_in->getIndices(index);
508 : double bf_value=1.0;
509 938 : for(unsigned int k=0; k<nargs; k++) {
510 530 : bf_value*=bf_values[k][indices[k]];
511 : }
512 408 : return bf_value;
513 408 : }
514 :
515 :
516 46 : void LinearBasisSetExpansion::setupUniformTargetDistribution() {
517 46 : std::vector< std::vector <double> > bf_integrals(0);
518 46 : std::vector<double> targetdist_averages(ncoeffs_,0.0);
519 : //
520 102 : for(unsigned int k=0; k<nargs_; k++) {
521 112 : bf_integrals.push_back(basisf_pntrs_[k]->getUniformIntegrals());
522 : }
523 : //
524 1677 : for(size_t i=0; i<ncoeffs_; i++) {
525 1631 : std::vector<unsigned int> indices=bias_coeffs_pntr_->getIndices(i);
526 : double value = 1.0;
527 4500 : for(unsigned int k=0; k<nargs_; k++) {
528 2869 : value*=bf_integrals[k][indices[k]];
529 : }
530 1631 : targetdist_averages[i]=value;
531 : }
532 46 : TargetDistAverages() = targetdist_averages;
533 46 : }
534 :
535 :
536 45 : void LinearBasisSetExpansion::setupTargetDistribution(TargetDistribution* targetdist_pntr_in) {
537 45 : targetdist_pntr_ = targetdist_pntr_in;
538 : //
539 45 : targetdist_pntr_->setupGrids(args_pntrs_,grid_min_,grid_max_,grid_bins_);
540 45 : targetdist_grid_pntr_ = targetdist_pntr_->getTargetDistGridPntr();
541 45 : log_targetdist_grid_pntr_ = targetdist_pntr_->getLogTargetDistGridPntr();
542 : //
543 45 : if(targetdist_pntr_->isDynamic()) {
544 39 : vesbias_pntr_->enableDynamicTargetDistribution();
545 : }
546 : //
547 45 : if(targetdist_pntr_->biasGridNeeded()) {
548 0 : setupBiasGrid(true);
549 0 : targetdist_pntr_->linkBiasGrid(bias_grid_pntr_.get());
550 : }
551 45 : if(targetdist_pntr_->biasWithoutCutoffGridNeeded()) {
552 3 : setupBiasGrid(true);
553 3 : targetdist_pntr_->linkBiasWithoutCutoffGrid(bias_withoutcutoff_grid_pntr_.get());
554 : }
555 45 : if(targetdist_pntr_->fesGridNeeded()) {
556 36 : setupFesGrid();
557 36 : targetdist_pntr_->linkFesGrid(fes_grid_pntr_.get());
558 : }
559 : //
560 45 : targetdist_pntr_->updateTargetDist();
561 45 : calculateTargetDistAveragesFromGrid(targetdist_grid_pntr_);
562 45 : }
563 :
564 :
565 355 : void LinearBasisSetExpansion::updateTargetDistribution() {
566 355 : plumed_massert(targetdist_pntr_!=NULL,"the target distribution hasn't been setup!");
567 355 : plumed_massert(targetdist_pntr_->isDynamic(),"this should only be used for dynamically updated target distributions!");
568 355 : if(targetdist_pntr_->biasGridNeeded()) {
569 0 : updateBiasGrid();
570 : }
571 355 : if(biasCutoffActive()) {
572 21 : updateBiasWithoutCutoffGrid();
573 : }
574 355 : if(targetdist_pntr_->fesGridNeeded()) {
575 334 : updateFesGrid();
576 : }
577 355 : targetdist_pntr_->updateTargetDist();
578 355 : calculateTargetDistAveragesFromGrid(targetdist_grid_pntr_);
579 355 : }
580 :
581 :
582 8 : void LinearBasisSetExpansion::readInRestartTargetDistribution(const std::string& grid_fname) {
583 8 : targetdist_pntr_->readInRestartTargetDistGrid(grid_fname);
584 8 : if(biasCutoffActive()) {
585 1 : targetdist_pntr_->clearLogTargetDistGrid();
586 : }
587 8 : }
588 :
589 :
590 8 : void LinearBasisSetExpansion::restartTargetDistribution() {
591 8 : plumed_massert(targetdist_pntr_!=NULL,"the target distribution hasn't been setup!");
592 8 : plumed_massert(targetdist_pntr_->isDynamic(),"this should only be used for dynamically updated target distributions!");
593 8 : if(biasCutoffActive()) {
594 1 : updateBiasWithoutCutoffGrid();
595 : }
596 8 : calculateTargetDistAveragesFromGrid(targetdist_grid_pntr_);
597 8 : }
598 :
599 :
600 408 : void LinearBasisSetExpansion::calculateTargetDistAveragesFromGrid(const Grid* targetdist_grid_pntr) {
601 408 : plumed_assert(targetdist_grid_pntr!=NULL);
602 408 : std::vector<double> targetdist_averages(ncoeffs_,0.0);
603 816 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(targetdist_grid_pntr);
604 408 : Grid::index_t stride=mycomm_.Get_size();
605 408 : Grid::index_t rank=mycomm_.Get_rank();
606 862723 : for(Grid::index_t l=rank; l<targetdist_grid_pntr->getSize(); l+=stride) {
607 862315 : std::vector<double> args_values = targetdist_grid_pntr->getPoint(l);
608 862315 : std::vector<double> basisset_values(ncoeffs_);
609 : // parallelization done over the grid -> should NOT use parallel in getBasisSetValues!!
610 862315 : getBasisSetValues(args_values,basisset_values,false);
611 862315 : double weight = integration_weights[l]*targetdist_grid_pntr->getValue(l);
612 95955507 : for(unsigned int i=0; i<ncoeffs_; i++) {
613 95093192 : targetdist_averages[i] += weight*basisset_values[i];
614 : }
615 : }
616 408 : mycomm_.Sum(targetdist_averages);
617 : // the overall constant;
618 408 : targetdist_averages[0] = getBasisSetConstant();
619 408 : TargetDistAverages() = targetdist_averages;
620 408 : }
621 :
622 :
623 39 : void LinearBasisSetExpansion::setBiasMinimumToZero() {
624 39 : plumed_massert(bias_grid_pntr_,"setBiasMinimumToZero can only be used if the bias grid is defined");
625 39 : updateBiasGrid();
626 39 : BiasCoeffs()[0]-=bias_grid_pntr_->getMinValue();
627 39 : }
628 :
629 :
630 0 : void LinearBasisSetExpansion::setBiasMaximumToZero() {
631 0 : plumed_massert(bias_grid_pntr_,"setBiasMaximumToZero can only be used if the bias grid is defined");
632 0 : updateBiasGrid();
633 0 : BiasCoeffs()[0]-=bias_grid_pntr_->getMaxValue();
634 0 : }
635 :
636 :
637 1982277 : bool LinearBasisSetExpansion::biasCutoffActive() const {
638 1982277 : if(vesbias_pntr_!=NULL) {
639 1474521 : return vesbias_pntr_->biasCutoffActive();
640 : } else {
641 : return false;
642 : }
643 : }
644 :
645 :
646 0 : double LinearBasisSetExpansion::calculateReweightFactor() const {
647 0 : plumed_massert(targetdist_grid_pntr_!=NULL,"calculateReweightFactor only be used if the target distribution grid is defined");
648 0 : plumed_massert(bias_grid_pntr_,"calculateReweightFactor only be used if the bias grid is defined");
649 : double sum = 0.0;
650 0 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(targetdist_grid_pntr_);
651 : //
652 0 : for(Grid::index_t l=0; l<targetdist_grid_pntr_->getSize(); l++) {
653 0 : sum += integration_weights[l] * targetdist_grid_pntr_->getValue(l) * exp(+beta_*bias_grid_pntr_->getValue(l));
654 : }
655 0 : if(sum==0.0) {
656 : sum=std::numeric_limits<double>::min();
657 : }
658 0 : return (1.0/beta_)*std::log(sum);
659 : }
660 :
661 :
662 :
663 :
664 : }
665 :
666 : }
|
