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 130 : 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 130 : CoeffsVector* bias_coeffs_pntr_in):
54 130 : label_(label),
55 130 : action_pntr_(NULL),
56 130 : vesbias_pntr_(NULL),
57 130 : mycomm_(cc),
58 130 : serial_(false),
59 130 : beta_(beta_in),
60 130 : kbt_(1.0/beta_),
61 130 : args_pntrs_(args_pntrs_in),
62 130 : nargs_(args_pntrs_.size()),
63 130 : basisf_pntrs_(basisf_pntrs_in),
64 130 : nbasisf_(basisf_pntrs_.size()),
65 130 : bias_coeffs_pntr_(bias_coeffs_pntr_in),
66 130 : ncoeffs_(0),
67 130 : grid_min_(nargs_),
68 130 : grid_max_(nargs_),
69 130 : grid_bins_(nargs_,100),
70 130 : targetdist_grid_label_("targetdist"),
71 130 : step_of_last_biasgrid_update(-1000),
72 130 : step_of_last_biaswithoutcutoffgrid_update(-1000),
73 130 : step_of_last_fesgrid_update(-1000),
74 130 : log_targetdist_grid_pntr_(NULL),
75 130 : targetdist_grid_pntr_(NULL),
76 130 : targetdist_pntr_(NULL),
77 260 : gfmt_("%14.9f") {
78 130 : plumed_massert(args_pntrs_.size()==basisf_pntrs_.size(),"number of arguments and basis functions do not match");
79 295 : for(unsigned int k=0; k<nargs_; k++) {
80 165 : nbasisf_[k]=basisf_pntrs_[k]->getNumberOfBasisFunctions();
81 : }
82 : //
83 130 : if(bias_coeffs_pntr_==NULL) {
84 0 : bias_coeffs_pntr_ = new CoeffsVector(label_+".coeffs",args_pntrs_,basisf_pntrs_,mycomm_,true);
85 : }
86 130 : plumed_massert(bias_coeffs_pntr_->numberOfDimensions()==basisf_pntrs_.size(),"dimension of coeffs does not match with number of basis functions ");
87 : //
88 130 : ncoeffs_ = bias_coeffs_pntr_->numberOfCoeffs();
89 130 : targetdist_averages_pntr_ = Tools::make_unique<CoeffsVector>(*bias_coeffs_pntr_);
90 :
91 130 : std::string targetdist_averages_label = bias_coeffs_pntr_->getLabel();
92 130 : if(targetdist_averages_label.find("coeffs")!=std::string::npos) {
93 260 : 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 130 : targetdist_averages_pntr_->setLabels(targetdist_averages_label);
98 : //
99 295 : for(unsigned int k=0; k<nargs_; k++) {
100 330 : grid_min_[k] = basisf_pntrs_[k]->intervalMinStr();
101 330 : grid_max_[k] = basisf_pntrs_[k]->intervalMaxStr();
102 : }
103 130 : }
104 :
105 130 : LinearBasisSetExpansion::~LinearBasisSetExpansion() {
106 390 : }
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 90 : void LinearBasisSetExpansion::linkVesBias(VesBias* vesbias_pntr_in) {
119 90 : vesbias_pntr_ = vesbias_pntr_in;
120 90 : action_pntr_ = static_cast<Action*>(vesbias_pntr_in);
121 :
122 90 : }
123 :
124 :
125 0 : void LinearBasisSetExpansion::linkAction(Action* action_pntr_in) {
126 0 : action_pntr_ = action_pntr_in;
127 0 : }
128 :
129 :
130 129 : void LinearBasisSetExpansion::setGridBins(const std::vector<unsigned int>& grid_bins_in) {
131 129 : plumed_massert(grid_bins_in.size()==nargs_,"the number of grid bins given doesn't match the number of arguments");
132 129 : plumed_massert(!bias_grid_pntr_,"setGridBins should be used before setting up the grids, otherwise it doesn't work");
133 129 : plumed_massert(!fes_grid_pntr_,"setGridBins should be used before setting up the grids, otherwise it doesn't work");
134 129 : grid_bins_=grid_bins_in;
135 129 : }
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 220 : std::unique_ptr<Grid> LinearBasisSetExpansion::setupGeneralGrid(const std::string& label_suffix, const bool usederiv) {
145 220 : bool use_spline = false;
146 440 : auto grid_pntr = Tools::make_unique<Grid>(label_+"."+label_suffix,args_pntrs_,grid_min_,grid_max_,grid_bins_,use_spline,usederiv);
147 220 : return grid_pntr;
148 : }
149 :
150 :
151 166 : void LinearBasisSetExpansion::setupBiasGrid(const bool usederiv) {
152 166 : if(bias_grid_pntr_) {
153 : return;
154 : }
155 258 : bias_grid_pntr_ = setupGeneralGrid("bias",usederiv);
156 129 : if(biasCutoffActive()) {
157 6 : 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 176 : fes_grid_pntr_ = setupGeneralGrid("fes",false);
170 : }
171 :
172 :
173 16 : void LinearBasisSetExpansion::setupFesProjGrid() {
174 8 : if(!fes_grid_pntr_) {
175 1 : setupFesGrid();
176 : }
177 8 : }
178 :
179 :
180 829 : void LinearBasisSetExpansion::updateBiasGrid() {
181 829 : plumed_massert(bias_grid_pntr_,"the bias grid is not defined");
182 829 : if(action_pntr_!=NULL && getStepOfLastBiasGridUpdate()==action_pntr_->getStep()) {
183 : return;
184 : }
185 1982290 : for(Grid::index_t l=0; l<bias_grid_pntr_->getSize(); l++) {
186 1981698 : std::vector<double> forces(nargs_);
187 1981698 : std::vector<double> args = bias_grid_pntr_->getPoint(l);
188 1981698 : bool all_inside=true;
189 1981698 : double bias=getBiasAndForces(args,all_inside,forces);
190 : //
191 1981698 : if(biasCutoffActive()) {
192 600 : vesbias_pntr_->applyBiasCutoff(bias,forces);
193 : }
194 : //
195 1981698 : if(bias_grid_pntr_->hasDerivatives()) {
196 1473981 : bias_grid_pntr_->setValueAndDerivatives(l,bias,forces);
197 : } else {
198 507717 : bias_grid_pntr_->setValue(l,bias);
199 : }
200 : //
201 : }
202 592 : if(vesbias_pntr_!=NULL) {
203 475 : vesbias_pntr_->setCurrentBiasMaxValue(bias_grid_pntr_->getMaxValue());
204 : }
205 592 : if(action_pntr_!=NULL) {
206 475 : 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 212 : void LinearBasisSetExpansion::writeBiasGridToFile(OFile& ofile, const bool append_file) const {
291 212 : plumed_massert(bias_grid_pntr_,"the bias grid is not defined");
292 212 : if(append_file) {
293 0 : ofile.enforceRestart();
294 : }
295 212 : bias_grid_pntr_->setOutputFmt(gfmt_);
296 212 : bias_grid_pntr_->writeToFile(ofile);
297 212 : }
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(gfmt_);
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(gfmt_);
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(gfmt_);
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(gfmt_);
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(gfmt_);
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(gfmt_);
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 2011787 : 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 2011787 : unsigned int nargs = args_values.size();
390 2011787 : plumed_assert(coeffs_pntr_in->numberOfDimensions()==nargs);
391 2011787 : plumed_assert(basisf_pntrs_in.size()==nargs);
392 2011787 : plumed_assert(forces.size()==nargs);
393 2011787 : plumed_assert(coeffsderivs_values.size()==coeffs_pntr_in->numberOfCoeffs());
394 :
395 2011787 : std::vector<double> args_values_trsfrm(nargs);
396 : // std::vector<bool> inside_interval(nargs,true);
397 2011787 : all_inside = true;
398 : //
399 2011787 : std::vector< std::vector <double> > bf_values(nargs);
400 2011787 : std::vector< std::vector <double> > bf_derivs(nargs);
401 : //
402 5966255 : for(unsigned int k=0; k<nargs; k++) {
403 3954468 : bf_values[k].assign(basisf_pntrs_in[k]->getNumberOfBasisFunctions(),0.0);
404 3954468 : bf_derivs[k].assign(basisf_pntrs_in[k]->getNumberOfBasisFunctions(),0.0);
405 3954468 : bool curr_inside=true;
406 3954468 : basisf_pntrs_in[k]->getAllValues(args_values[k],args_values_trsfrm[k],curr_inside,bf_values[k],bf_derivs[k]);
407 : // inside_interval[k]=curr_inside;
408 3954468 : if(!curr_inside) {
409 6908 : all_inside=false;
410 : }
411 3954468 : forces[k]=0.0;
412 : }
413 : //
414 : size_t stride=1;
415 : size_t rank=0;
416 2011787 : if(comm_in!=NULL) {
417 2011787 : stride=comm_in->Get_size();
418 2011787 : rank=comm_in->Get_rank();
419 : }
420 : // loop over coeffs
421 2011787 : double bias=0.0;
422 144769949 : for(size_t i=rank; i<coeffs_pntr_in->numberOfCoeffs(); i+=stride) {
423 142758162 : std::vector<unsigned int> indices=coeffs_pntr_in->getIndices(i);
424 : double coeff = coeffs_pntr_in->getValue(i);
425 : double bf_curr=1.0;
426 435044808 : for(unsigned int k=0; k<nargs; k++) {
427 292286646 : bf_curr*=bf_values[k][indices[k]];
428 : }
429 142758162 : bias+=coeff*bf_curr;
430 142758162 : coeffsderivs_values[i] = bf_curr;
431 435044808 : for(unsigned int k=0; k<nargs; k++) {
432 : double der = 1.0;
433 898592256 : for(unsigned int l=0; l<nargs; l++) {
434 606305610 : if(l!=k) {
435 314018964 : der*=bf_values[l][indices[l]];
436 : } else {
437 292286646 : der*=bf_derivs[l][indices[l]];
438 : }
439 : }
440 292286646 : forces[k]-=coeff*der;
441 : // maybe faster but dangerous
442 : // forces[k]-=coeff*bf_curr*(bf_derivs[k][indices[k]]/bf_values[k][indices[k]]);
443 : }
444 : }
445 : //
446 2011787 : if(comm_in!=NULL) {
447 : // coeffsderivs_values is not summed as the mpi Sum is done later on for the averages
448 2011787 : comm_in->Sum(bias);
449 2011787 : comm_in->Sum(forces);
450 : }
451 2011787 : return bias;
452 2011787 : }
453 :
454 :
455 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) {
456 862315 : unsigned int nargs = args_values.size();
457 862315 : plumed_assert(coeffs_pntr_in->numberOfDimensions()==nargs);
458 862315 : plumed_assert(basisf_pntrs_in.size()==nargs);
459 :
460 862315 : std::vector<double> args_values_trsfrm(nargs);
461 : std::vector< std::vector <double> > bf_values;
462 : //
463 2583312 : for(unsigned int k=0; k<nargs; k++) {
464 1720997 : std::vector<double> tmp_val(basisf_pntrs_in[k]->getNumberOfBasisFunctions());
465 1720997 : std::vector<double> tmp_der(tmp_val.size());
466 1720997 : bool inside=true;
467 1720997 : basisf_pntrs_in[k]->getAllValues(args_values[k],args_values_trsfrm[k],inside,tmp_val,tmp_der);
468 1720997 : bf_values.push_back(tmp_val);
469 : }
470 : //
471 : size_t stride=1;
472 : size_t rank=0;
473 862315 : if(comm_in!=NULL) {
474 0 : stride=comm_in->Get_size();
475 0 : rank=comm_in->Get_rank();
476 : }
477 : // loop over basis set
478 95955507 : for(size_t i=rank; i<coeffs_pntr_in->numberOfCoeffs(); i+=stride) {
479 95093192 : std::vector<unsigned int> indices=coeffs_pntr_in->getIndices(i);
480 : double bf_curr=1.0;
481 298507612 : for(unsigned int k=0; k<nargs; k++) {
482 203414420 : bf_curr*=bf_values[k][indices[k]];
483 : }
484 95093192 : basisset_values[i] = bf_curr;
485 : }
486 : //
487 862315 : if(comm_in!=NULL) {
488 0 : comm_in->Sum(basisset_values);
489 : }
490 1724630 : }
491 :
492 :
493 408 : double LinearBasisSetExpansion::getBasisSetValue(const std::vector<double>& args_values, const size_t index, std::vector<BasisFunctions*>& basisf_pntrs_in, CoeffsVector* coeffs_pntr_in) {
494 408 : unsigned int nargs = args_values.size();
495 408 : plumed_assert(coeffs_pntr_in->numberOfDimensions()==nargs);
496 408 : plumed_assert(basisf_pntrs_in.size()==nargs);
497 :
498 408 : std::vector<double> args_values_trsfrm(nargs);
499 : std::vector< std::vector <double> > bf_values;
500 : //
501 938 : for(unsigned int k=0; k<nargs; k++) {
502 530 : std::vector<double> tmp_val(basisf_pntrs_in[k]->getNumberOfBasisFunctions());
503 530 : std::vector<double> tmp_der(tmp_val.size());
504 530 : bool inside=true;
505 530 : basisf_pntrs_in[k]->getAllValues(args_values[k],args_values_trsfrm[k],inside,tmp_val,tmp_der);
506 530 : bf_values.push_back(tmp_val);
507 : }
508 : //
509 408 : std::vector<unsigned int> indices=coeffs_pntr_in->getIndices(index);
510 : double bf_value=1.0;
511 938 : for(unsigned int k=0; k<nargs; k++) {
512 530 : bf_value*=bf_values[k][indices[k]];
513 : }
514 408 : return bf_value;
515 408 : }
516 :
517 :
518 45 : void LinearBasisSetExpansion::setupUniformTargetDistribution() {
519 45 : std::vector< std::vector <double> > bf_integrals(0);
520 45 : std::vector<double> targetdist_averages(ncoeffs_,0.0);
521 : //
522 100 : for(unsigned int k=0; k<nargs_; k++) {
523 110 : bf_integrals.push_back(basisf_pntrs_[k]->getUniformIntegrals());
524 : }
525 : //
526 1672 : for(size_t i=0; i<ncoeffs_; i++) {
527 1627 : std::vector<unsigned int> indices=bias_coeffs_pntr_->getIndices(i);
528 : double value = 1.0;
529 4492 : for(unsigned int k=0; k<nargs_; k++) {
530 2865 : value*=bf_integrals[k][indices[k]];
531 : }
532 1627 : targetdist_averages[i]=value;
533 : }
534 45 : TargetDistAverages() = targetdist_averages;
535 45 : }
536 :
537 :
538 45 : void LinearBasisSetExpansion::setupTargetDistribution(TargetDistribution* targetdist_pntr_in) {
539 45 : targetdist_pntr_ = targetdist_pntr_in;
540 : //
541 45 : targetdist_pntr_->setupGrids(args_pntrs_,grid_min_,grid_max_,grid_bins_);
542 45 : targetdist_grid_pntr_ = targetdist_pntr_->getTargetDistGridPntr();
543 45 : log_targetdist_grid_pntr_ = targetdist_pntr_->getLogTargetDistGridPntr();
544 : //
545 45 : if(targetdist_pntr_->isDynamic()) {
546 39 : vesbias_pntr_->enableDynamicTargetDistribution();
547 : }
548 : //
549 45 : if(targetdist_pntr_->biasGridNeeded()) {
550 0 : setupBiasGrid(true);
551 0 : targetdist_pntr_->linkBiasGrid(bias_grid_pntr_.get());
552 : }
553 45 : if(targetdist_pntr_->biasWithoutCutoffGridNeeded()) {
554 3 : setupBiasGrid(true);
555 3 : targetdist_pntr_->linkBiasWithoutCutoffGrid(bias_withoutcutoff_grid_pntr_.get());
556 : }
557 45 : if(targetdist_pntr_->fesGridNeeded()) {
558 36 : setupFesGrid();
559 36 : targetdist_pntr_->linkFesGrid(fes_grid_pntr_.get());
560 : }
561 : //
562 45 : targetdist_pntr_->updateTargetDist();
563 45 : calculateTargetDistAveragesFromGrid(targetdist_grid_pntr_);
564 45 : }
565 :
566 :
567 355 : void LinearBasisSetExpansion::updateTargetDistribution() {
568 355 : plumed_massert(targetdist_pntr_!=NULL,"the target distribution hasn't been setup!");
569 355 : plumed_massert(targetdist_pntr_->isDynamic(),"this should only be used for dynamically updated target distributions!");
570 355 : if(targetdist_pntr_->biasGridNeeded()) {
571 0 : updateBiasGrid();
572 : }
573 355 : if(biasCutoffActive()) {
574 21 : updateBiasWithoutCutoffGrid();
575 : }
576 355 : if(targetdist_pntr_->fesGridNeeded()) {
577 334 : updateFesGrid();
578 : }
579 355 : targetdist_pntr_->updateTargetDist();
580 355 : calculateTargetDistAveragesFromGrid(targetdist_grid_pntr_);
581 355 : }
582 :
583 :
584 8 : void LinearBasisSetExpansion::readInRestartTargetDistribution(const std::string& grid_fname) {
585 8 : targetdist_pntr_->readInRestartTargetDistGrid(grid_fname);
586 8 : if(biasCutoffActive()) {
587 1 : targetdist_pntr_->clearLogTargetDistGrid();
588 : }
589 8 : }
590 :
591 :
592 8 : void LinearBasisSetExpansion::restartTargetDistribution() {
593 8 : plumed_massert(targetdist_pntr_!=NULL,"the target distribution hasn't been setup!");
594 8 : plumed_massert(targetdist_pntr_->isDynamic(),"this should only be used for dynamically updated target distributions!");
595 8 : if(biasCutoffActive()) {
596 1 : updateBiasWithoutCutoffGrid();
597 : }
598 8 : calculateTargetDistAveragesFromGrid(targetdist_grid_pntr_);
599 8 : }
600 :
601 :
602 408 : void LinearBasisSetExpansion::calculateTargetDistAveragesFromGrid(const Grid* targetdist_grid_pntr) {
603 408 : plumed_assert(targetdist_grid_pntr!=NULL);
604 408 : std::vector<double> targetdist_averages(ncoeffs_,0.0);
605 816 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(targetdist_grid_pntr);
606 408 : Grid::index_t stride=mycomm_.Get_size();
607 408 : Grid::index_t rank=mycomm_.Get_rank();
608 862723 : for(Grid::index_t l=rank; l<targetdist_grid_pntr->getSize(); l+=stride) {
609 862315 : std::vector<double> args_values = targetdist_grid_pntr->getPoint(l);
610 862315 : std::vector<double> basisset_values(ncoeffs_);
611 : // parallelization done over the grid -> should NOT use parallel in getBasisSetValues!!
612 862315 : getBasisSetValues(args_values,basisset_values,false);
613 862315 : double weight = integration_weights[l]*targetdist_grid_pntr->getValue(l);
614 95955507 : for(unsigned int i=0; i<ncoeffs_; i++) {
615 95093192 : targetdist_averages[i] += weight*basisset_values[i];
616 : }
617 : }
618 408 : mycomm_.Sum(targetdist_averages);
619 : // the overall constant;
620 408 : targetdist_averages[0] = getBasisSetConstant();
621 408 : TargetDistAverages() = targetdist_averages;
622 408 : }
623 :
624 :
625 39 : void LinearBasisSetExpansion::setBiasMinimumToZero() {
626 39 : plumed_massert(bias_grid_pntr_,"setBiasMinimumToZero can only be used if the bias grid is defined");
627 39 : updateBiasGrid();
628 39 : BiasCoeffs()[0]-=bias_grid_pntr_->getMinValue();
629 39 : }
630 :
631 :
632 0 : void LinearBasisSetExpansion::setBiasMaximumToZero() {
633 0 : plumed_massert(bias_grid_pntr_,"setBiasMaximumToZero can only be used if the bias grid is defined");
634 0 : updateBiasGrid();
635 0 : BiasCoeffs()[0]-=bias_grid_pntr_->getMaxValue();
636 0 : }
637 :
638 :
639 1982225 : bool LinearBasisSetExpansion::biasCutoffActive() const {
640 1982225 : if(vesbias_pntr_!=NULL) {
641 1474469 : return vesbias_pntr_->biasCutoffActive();
642 : } else {
643 : return false;
644 : }
645 : }
646 :
647 :
648 0 : double LinearBasisSetExpansion::calculateReweightFactor() const {
649 0 : plumed_massert(targetdist_grid_pntr_!=NULL,"calculateReweightFactor only be used if the target distribution grid is defined");
650 0 : plumed_massert(bias_grid_pntr_,"calculateReweightFactor only be used if the bias grid is defined");
651 : double sum = 0.0;
652 0 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(targetdist_grid_pntr_);
653 : //
654 0 : for(Grid::index_t l=0; l<targetdist_grid_pntr_->getSize(); l++) {
655 0 : sum += integration_weights[l] * targetdist_grid_pntr_->getValue(l) * exp(+beta_*bias_grid_pntr_->getValue(l));
656 : }
657 0 : if(sum==0.0) {
658 : sum=std::numeric_limits<double>::min();
659 : }
660 0 : return (1.0/beta_)*std::log(sum);
661 : }
662 :
663 :
664 :
665 :
666 : }
667 :
668 : }
|
