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