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 119 : 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 119 : 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 119 : 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 119 : grid_min_(nargs_),
67 119 : grid_max_(nargs_),
68 119 : 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 952 : targetdist_pntr_(NULL)
79 : {
80 119 : plumed_massert(args_pntrs_.size()==basisf_pntrs_.size(),"number of arguments and basis functions do not match");
81 563 : for(unsigned int k=0; k<nargs_; k++) {nbasisf_[k]=basisf_pntrs_[k]->getNumberOfBasisFunctions();}
82 : //
83 119 : if(bias_coeffs_pntr_==NULL) {
84 0 : bias_coeffs_pntr_ = new CoeffsVector(label_+".coeffs",args_pntrs_,basisf_pntrs_,mycomm_,true);
85 : }
86 357 : plumed_massert(bias_coeffs_pntr_->numberOfDimensions()==basisf_pntrs_.size(),"dimension of coeffs does not match with number of basis functions ");
87 : //
88 119 : ncoeffs_ = bias_coeffs_pntr_->numberOfCoeffs();
89 119 : targetdist_averages_pntr_ = new CoeffsVector(*bias_coeffs_pntr_);
90 :
91 119 : std::string targetdist_averages_label = bias_coeffs_pntr_->getLabel();
92 119 : if(targetdist_averages_label.find("coeffs")!=std::string::npos) {
93 357 : 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 119 : targetdist_averages_pntr_->setLabels(targetdist_averages_label);
99 : //
100 415 : for(unsigned int k=0; k<nargs_; k++) {
101 592 : grid_min_[k] = basisf_pntrs_[k]->intervalMinStr();
102 444 : grid_max_[k] = basisf_pntrs_[k]->intervalMaxStr();
103 : }
104 119 : }
105 :
106 357 : LinearBasisSetExpansion::~LinearBasisSetExpansion() {
107 119 : if(targetdist_averages_pntr_!=NULL) {
108 119 : delete targetdist_averages_pntr_;
109 : }
110 119 : if(bias_grid_pntr_!=NULL) {
111 119 : delete bias_grid_pntr_;
112 : }
113 119 : if(bias_withoutcutoff_grid_pntr_!=NULL) {
114 3 : delete bias_withoutcutoff_grid_pntr_;
115 : }
116 119 : if(fes_grid_pntr_!=NULL) {
117 80 : delete fes_grid_pntr_;
118 : }
119 119 : }
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 82 : void LinearBasisSetExpansion::linkVesBias(VesBias* vesbias_pntr_in) {
132 82 : vesbias_pntr_ = vesbias_pntr_in;
133 82 : action_pntr_ = static_cast<Action*>(vesbias_pntr_in);
134 :
135 82 : }
136 :
137 :
138 0 : void LinearBasisSetExpansion::linkAction(Action* action_pntr_in) {
139 0 : action_pntr_ = action_pntr_in;
140 0 : }
141 :
142 :
143 119 : void LinearBasisSetExpansion::setGridBins(const std::vector<unsigned int>& grid_bins_in) {
144 119 : plumed_massert(grid_bins_in.size()==nargs_,"the number of grid bins given doesn't match the number of arguments");
145 119 : plumed_massert(bias_grid_pntr_==NULL,"setGridBins should be used before setting up the grids, otherwise it doesn't work");
146 119 : plumed_massert(fes_grid_pntr_==NULL,"setGridBins should be used before setting up the grids, otherwise it doesn't work");
147 119 : grid_bins_=grid_bins_in;
148 119 : }
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 202 : Grid* LinearBasisSetExpansion::setupGeneralGrid(const std::string& label_suffix, const bool usederiv) {
158 : bool use_spline = false;
159 606 : Grid* grid_pntr = new Grid(label_+"."+label_suffix,args_pntrs_,grid_min_,grid_max_,grid_bins_,use_spline,usederiv);
160 202 : return grid_pntr;
161 : }
162 :
163 :
164 155 : void LinearBasisSetExpansion::setupBiasGrid(const bool usederiv) {
165 155 : if(bias_grid_pntr_!=NULL) {return;}
166 238 : bias_grid_pntr_ = setupGeneralGrid("bias",usederiv);
167 119 : if(biasCutoffActive()) {
168 6 : bias_withoutcutoff_grid_pntr_ = setupGeneralGrid("bias_withoutcutoff",usederiv);
169 : }
170 : }
171 :
172 :
173 111 : void LinearBasisSetExpansion::setupFesGrid() {
174 111 : if(fes_grid_pntr_!=NULL) {return;}
175 80 : if(bias_grid_pntr_==NULL) {
176 32 : setupBiasGrid(true);
177 : }
178 160 : 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 783 : void LinearBasisSetExpansion::updateBiasGrid() {
190 783 : plumed_massert(bias_grid_pntr_!=NULL,"the bias grid is not defined");
191 1455 : if(action_pntr_!=NULL && getStepOfLastBiasGridUpdate()==action_pntr_->getStep()) {
192 : return;
193 : }
194 3785769 : for(Grid::index_t l=0; l<bias_grid_pntr_->getSize(); l++) {
195 1892607 : std::vector<double> forces(nargs_);
196 1892607 : std::vector<double> args = bias_grid_pntr_->getPoint(l);
197 1892607 : bool all_inside=true;
198 1892607 : double bias=getBiasAndForces(args,all_inside,forces);
199 : //
200 1892607 : if(biasCutoffActive()) {
201 600 : vesbias_pntr_->applyBiasCutoff(bias,forces);
202 : }
203 : //
204 3785214 : if(bias_grid_pntr_->hasDerivatives()) {
205 1454196 : bias_grid_pntr_->setValueAndDerivatives(l,bias,forces);
206 : }
207 : else {
208 438411 : bias_grid_pntr_->setValue(l,bias);
209 : }
210 : //
211 : }
212 555 : if(vesbias_pntr_!=NULL) {
213 444 : vesbias_pntr_->setCurrentBiasMaxValue(bias_grid_pntr_->getMaxValue());
214 : }
215 555 : if(action_pntr_!=NULL) {
216 444 : 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 : return;
227 : }
228 : //
229 4823 : 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 4622 : 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 507 : void LinearBasisSetExpansion::updateFesGrid() {
281 507 : plumed_massert(fes_grid_pntr_!=NULL,"the FES grid is not defined");
282 507 : updateBiasGrid();
283 1014 : if(action_pntr_!=NULL && getStepOfLastFesGridUpdate() == action_pntr_->getStep()) {
284 : return;
285 : }
286 : //
287 : double bias2fes_scalingf = -1.0;
288 2908234 : for(Grid::index_t l=0; l<fes_grid_pntr_->getSize(); l++) {
289 1453896 : double fes_value = bias2fes_scalingf*bias_grid_pntr_->getValue(l);
290 1453896 : if(log_targetdist_grid_pntr_!=NULL) {
291 1204866 : fes_value += kBT()*log_targetdist_grid_pntr_->getValue(l);
292 : }
293 1453896 : fes_grid_pntr_->setValue(l,fes_value);
294 : }
295 442 : fes_grid_pntr_->setMinToZero();
296 442 : if(action_pntr_!=NULL) {
297 442 : setStepOfLastFesGridUpdate(action_pntr_->getStep());
298 : }
299 : }
300 :
301 :
302 202 : void LinearBasisSetExpansion::writeBiasGridToFile(OFile& ofile, const bool append_file) const {
303 202 : plumed_massert(bias_grid_pntr_!=NULL,"the bias grid is not defined");
304 202 : if(append_file) {ofile.enforceRestart();}
305 202 : bias_grid_pntr_->writeToFile(ofile);
306 202 : }
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 161 : void LinearBasisSetExpansion::writeFesGridToFile(OFile& ofile, const bool append_file) const {
317 161 : plumed_massert(fes_grid_pntr_!=NULL,"the FES grid is not defined");
318 161 : if(append_file) {ofile.enforceRestart();}
319 161 : fes_grid_pntr_->writeToFile(ofile);
320 161 : }
321 :
322 :
323 36 : void LinearBasisSetExpansion::writeFesProjGridToFile(const std::vector<std::string>& proj_arg, OFile& ofile, const bool append_file) const {
324 36 : plumed_massert(fes_grid_pntr_!=NULL,"the FES grid is not defined");
325 36 : FesWeight* Fw = new FesWeight(beta_);
326 72 : Grid proj_grid = fes_grid_pntr_->project(proj_arg,Fw);
327 36 : proj_grid.setMinToZero();
328 36 : if(append_file) {ofile.enforceRestart();}
329 36 : proj_grid.writeToFile(ofile);
330 36 : delete Fw;
331 36 : }
332 :
333 :
334 72 : void LinearBasisSetExpansion::writeTargetDistGridToFile(OFile& ofile, const bool append_file) const {
335 72 : if(targetdist_grid_pntr_==NULL) {return;}
336 72 : if(append_file) {ofile.enforceRestart();}
337 72 : targetdist_grid_pntr_->writeToFile(ofile);
338 : }
339 :
340 :
341 72 : void LinearBasisSetExpansion::writeLogTargetDistGridToFile(OFile& ofile, const bool append_file) const {
342 72 : if(log_targetdist_grid_pntr_==NULL) {return;}
343 72 : if(append_file) {ofile.enforceRestart();}
344 72 : log_targetdist_grid_pntr_->writeToFile(ofile);
345 : }
346 :
347 :
348 20 : void LinearBasisSetExpansion::writeTargetDistProjGridToFile(const std::vector<std::string>& proj_arg, OFile& ofile, const bool append_file) const {
349 20 : if(targetdist_grid_pntr_==NULL) {return;}
350 20 : if(append_file) {ofile.enforceRestart();}
351 40 : Grid proj_grid = TargetDistribution::getMarginalDistributionGrid(targetdist_grid_pntr_,proj_arg);
352 20 : 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 1900406 : 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 1900406 : unsigned int nargs = args_values.size();
372 1900406 : plumed_assert(coeffs_pntr_in->numberOfDimensions()==nargs);
373 1900406 : plumed_assert(basisf_pntrs_in.size()==nargs);
374 1900406 : plumed_assert(forces.size()==nargs);
375 1900406 : plumed_assert(coeffsderivs_values.size()==coeffs_pntr_in->numberOfCoeffs());
376 :
377 1900406 : std::vector<double> args_values_trsfrm(nargs);
378 : // std::vector<bool> inside_interval(nargs,true);
379 1900406 : all_inside = true;
380 : //
381 3800812 : std::vector< std::vector <double> > bf_values(nargs);
382 3800812 : std::vector< std::vector <double> > bf_derivs(nargs);
383 : //
384 9370914 : for(unsigned int k=0; k<nargs; k++) {
385 18676270 : bf_values[k].assign(basisf_pntrs_in[k]->getNumberOfBasisFunctions(),0.0);
386 14941016 : bf_derivs[k].assign(basisf_pntrs_in[k]->getNumberOfBasisFunctions(),0.0);
387 3735254 : bool curr_inside=true;
388 7470508 : 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 3735254 : if(!curr_inside) {all_inside=false;}
391 3735254 : forces[k]=0.0;
392 : }
393 : //
394 : size_t stride=1;
395 : size_t rank=0;
396 1900406 : if(comm_in!=NULL)
397 : {
398 1900406 : stride=comm_in->Get_size();
399 1900406 : rank=comm_in->Get_rank();
400 : }
401 : // loop over coeffs
402 1900406 : double bias=0.0;
403 265330016 : for(size_t i=rank; i<coeffs_pntr_in->numberOfCoeffs(); i+=stride) {
404 131714805 : 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 657517621 : for(unsigned int k=0; k<nargs; k++) {
408 788704224 : bf_curr*=bf_values[k][indices[k]];
409 : }
410 131714805 : bias+=coeff*bf_curr;
411 131714805 : coeffsderivs_values[i] = bf_curr;
412 657517621 : for(unsigned int k=0; k<nargs; k++) {
413 : double der = 1.0;
414 1313450636 : for(unsigned int l=0; l<nargs; l++) {
415 1050021026 : if(l!=k) {der*=bf_values[l][indices[l]];}
416 788704224 : else {der*=bf_derivs[l][indices[l]];}
417 : }
418 525802816 : 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 1900406 : if(comm_in!=NULL) {
425 : // coeffsderivs_values is not summed as the mpi Sum is done later on for the averages
426 1900406 : comm_in->Sum(bias);
427 1900406 : comm_in->Sum(forces);
428 : }
429 3800812 : return bias;
430 : }
431 :
432 :
433 828566 : 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 828566 : unsigned int nargs = args_values.size();
435 828566 : plumed_assert(coeffs_pntr_in->numberOfDimensions()==nargs);
436 828566 : plumed_assert(basisf_pntrs_in.size()==nargs);
437 :
438 828566 : std::vector<double> args_values_trsfrm(nargs);
439 828566 : std::vector< std::vector <double> > bf_values;
440 : //
441 4102964 : for(unsigned int k=0; k<nargs; k++) {
442 4911597 : std::vector<double> tmp_val(basisf_pntrs_in[k]->getNumberOfBasisFunctions());
443 1637199 : std::vector<double> tmp_der(tmp_val.size());
444 1637199 : bool inside=true;
445 3274398 : basisf_pntrs_in[k]->getAllValues(args_values[k],args_values_trsfrm[k],inside,tmp_val,tmp_der);
446 1637199 : bf_values.push_back(tmp_val);
447 : }
448 : //
449 : size_t stride=1;
450 : size_t rank=0;
451 828566 : 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 161378026 : for(size_t i=rank; i<coeffs_pntr_in->numberOfCoeffs(); i+=stride) {
458 80274730 : std::vector<unsigned int> indices=coeffs_pntr_in->getIndices(i);
459 : double bf_curr=1.0;
460 400933524 : for(unsigned int k=0; k<nargs; k++) {
461 480988191 : bf_curr*=bf_values[k][indices[k]];
462 : }
463 80274730 : basisset_values[i] = bf_curr;
464 : }
465 : //
466 828566 : if(comm_in!=NULL) {
467 0 : comm_in->Sum(basisset_values);
468 : }
469 828566 : }
470 :
471 :
472 379 : double LinearBasisSetExpansion::getBasisSetValue(const std::vector<double>& args_values, const size_t index, std::vector<BasisFunctions*>& basisf_pntrs_in, CoeffsVector* coeffs_pntr_in) {
473 379 : unsigned int nargs = args_values.size();
474 379 : plumed_assert(coeffs_pntr_in->numberOfDimensions()==nargs);
475 379 : plumed_assert(basisf_pntrs_in.size()==nargs);
476 :
477 379 : std::vector<double> args_values_trsfrm(nargs);
478 379 : std::vector< std::vector <double> > bf_values;
479 : //
480 1363 : for(unsigned int k=0; k<nargs; k++) {
481 1476 : std::vector<double> tmp_val(basisf_pntrs_in[k]->getNumberOfBasisFunctions());
482 492 : std::vector<double> tmp_der(tmp_val.size());
483 492 : bool inside=true;
484 984 : basisf_pntrs_in[k]->getAllValues(args_values[k],args_values_trsfrm[k],inside,tmp_val,tmp_der);
485 492 : bf_values.push_back(tmp_val);
486 : }
487 : //
488 379 : std::vector<unsigned int> indices=coeffs_pntr_in->getIndices(index);
489 : double bf_value=1.0;
490 1363 : for(unsigned int k=0; k<nargs; k++) {
491 1476 : bf_value*=bf_values[k][indices[k]];
492 : }
493 379 : return bf_value;
494 : }
495 :
496 :
497 42 : void LinearBasisSetExpansion::setupUniformTargetDistribution() {
498 84 : std::vector< std::vector <double> > bf_integrals(0);
499 42 : std::vector<double> targetdist_averages(ncoeffs_,0.0);
500 : //
501 146 : for(unsigned int k=0; k<nargs_; k++) {
502 156 : bf_integrals.push_back(basisf_pntrs_[k]->getUniformIntegrals());
503 : }
504 : //
505 3230 : for(size_t i=0; i<ncoeffs_; i++) {
506 1594 : std::vector<unsigned int> indices=bias_coeffs_pntr_->getIndices(i);
507 : double value = 1.0;
508 7258 : for(unsigned int k=0; k<nargs_; k++) {
509 8496 : value*=bf_integrals[k][indices[k]];
510 : }
511 1594 : targetdist_averages[i]=value;
512 : }
513 42 : TargetDistAverages() = targetdist_averages;
514 42 : }
515 :
516 :
517 40 : void LinearBasisSetExpansion::setupTargetDistribution(TargetDistribution* targetdist_pntr_in) {
518 40 : targetdist_pntr_ = targetdist_pntr_in;
519 : //
520 40 : targetdist_pntr_->setupGrids(args_pntrs_,grid_min_,grid_max_,grid_bins_);
521 80 : targetdist_grid_pntr_ = targetdist_pntr_->getTargetDistGridPntr();
522 40 : log_targetdist_grid_pntr_ = targetdist_pntr_->getLogTargetDistGridPntr();
523 : //
524 40 : if(targetdist_pntr_->isDynamic()) {
525 34 : vesbias_pntr_->enableDynamicTargetDistribution();
526 : }
527 : //
528 80 : if(targetdist_pntr_->biasGridNeeded()) {
529 0 : setupBiasGrid(true);
530 0 : targetdist_pntr_->linkBiasGrid(bias_grid_pntr_);
531 : }
532 80 : if(targetdist_pntr_->biasWithoutCutoffGridNeeded()) {
533 3 : setupBiasGrid(true);
534 3 : targetdist_pntr_->linkBiasWithoutCutoffGrid(bias_withoutcutoff_grid_pntr_);
535 : }
536 80 : if(targetdist_pntr_->fesGridNeeded()) {
537 31 : setupFesGrid();
538 31 : targetdist_pntr_->linkFesGrid(fes_grid_pntr_);
539 : }
540 : //
541 40 : targetdist_pntr_->updateTargetDist();
542 40 : calculateTargetDistAveragesFromGrid(targetdist_grid_pntr_);
543 40 : }
544 :
545 :
546 331 : void LinearBasisSetExpansion::updateTargetDistribution() {
547 331 : plumed_massert(targetdist_pntr_!=NULL,"the target distribution hasn't been setup!");
548 331 : plumed_massert(targetdist_pntr_->isDynamic(),"this should only be used for dynamically updated target distributions!");
549 331 : if(targetdist_pntr_->biasGridNeeded()) {updateBiasGrid();}
550 331 : if(biasCutoffActive()) {updateBiasWithoutCutoffGrid();}
551 662 : if(targetdist_pntr_->fesGridNeeded()) {updateFesGrid();}
552 331 : targetdist_pntr_->updateTargetDist();
553 331 : calculateTargetDistAveragesFromGrid(targetdist_grid_pntr_);
554 331 : }
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 379 : void LinearBasisSetExpansion::calculateTargetDistAveragesFromGrid(const Grid* targetdist_grid_pntr) {
572 379 : plumed_assert(targetdist_grid_pntr!=NULL);
573 379 : std::vector<double> targetdist_averages(ncoeffs_,0.0);
574 1137 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(targetdist_grid_pntr);
575 379 : Grid::index_t stride=mycomm_.Get_size();
576 379 : Grid::index_t rank=mycomm_.Get_rank();
577 1657511 : for(Grid::index_t l=rank; l<targetdist_grid_pntr->getSize(); l+=stride) {
578 828566 : std::vector<double> args_values = targetdist_grid_pntr->getPoint(l);
579 828566 : std::vector<double> basisset_values(ncoeffs_);
580 : // parallelization done over the grid -> should NOT use parallel in getBasisSetValues!!
581 828566 : getBasisSetValues(args_values,basisset_values,false);
582 828566 : double weight = integration_weights[l]*targetdist_grid_pntr->getValue(l);
583 161378026 : for(unsigned int i=0; i<ncoeffs_; i++) {
584 160549460 : targetdist_averages[i] += weight*basisset_values[i];
585 : }
586 : }
587 379 : mycomm_.Sum(targetdist_averages);
588 : // the overall constant;
589 379 : targetdist_averages[0] = getBasisSetConstant();
590 379 : TargetDistAverages() = targetdist_averages;
591 379 : }
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 1893100 : bool LinearBasisSetExpansion::biasCutoffActive() const {
609 3347752 : 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 4839 : }
|
