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