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