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 118 : 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 118 : kbt_(1.0/beta_),
59 : args_pntrs_(args_pntrs_in),
60 118 : 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 354 : targetdist_pntr_(NULL)
79 : {
80 118 : plumed_massert(args_pntrs_.size()==basisf_pntrs_.size(),"number of arguments and basis functions do not match");
81 118 : for(unsigned int k=0; k<nargs_; k++) {nbasisf_[k]=basisf_pntrs_[k]->getNumberOfBasisFunctions();}
82 : //
83 118 : if(bias_coeffs_pntr_==NULL) {
84 0 : bias_coeffs_pntr_ = new CoeffsVector(label_+".coeffs",args_pntrs_,basisf_pntrs_,mycomm_,true);
85 : }
86 118 : plumed_massert(bias_coeffs_pntr_->numberOfDimensions()==basisf_pntrs_.size(),"dimension of coeffs does not match with number of basis functions ");
87 : //
88 118 : ncoeffs_ = bias_coeffs_pntr_->numberOfCoeffs();
89 118 : targetdist_averages_pntr_ = new CoeffsVector(*bias_coeffs_pntr_);
90 :
91 118 : std::string targetdist_averages_label = bias_coeffs_pntr_->getLabel();
92 118 : if(targetdist_averages_label.find("coeffs")!=std::string::npos) {
93 118 : targetdist_averages_label.replace(targetdist_averages_label.find("coeffs"), std::string("coeffs").length(), "targetdist_averages");
94 : }
95 : else {
96 0 : targetdist_averages_label += "_targetdist_averages";
97 : }
98 118 : targetdist_averages_pntr_->setLabels(targetdist_averages_label);
99 : //
100 264 : for(unsigned int k=0; k<nargs_; k++) {
101 146 : grid_min_[k] = basisf_pntrs_[k]->intervalMinStr();
102 146 : grid_max_[k] = basisf_pntrs_[k]->intervalMaxStr();
103 118 : }
104 118 : }
105 :
106 236 : LinearBasisSetExpansion::~LinearBasisSetExpansion() {
107 118 : if(targetdist_averages_pntr_!=NULL) {
108 118 : delete targetdist_averages_pntr_;
109 : }
110 118 : if(bias_grid_pntr_!=NULL) {
111 118 : delete bias_grid_pntr_;
112 : }
113 118 : if(bias_withoutcutoff_grid_pntr_!=NULL) {
114 3 : delete bias_withoutcutoff_grid_pntr_;
115 : }
116 118 : if(fes_grid_pntr_!=NULL) {
117 79 : delete fes_grid_pntr_;
118 : }
119 118 : }
120 :
121 :
122 0 : bool LinearBasisSetExpansion::isStaticTargetDistFileOutputActive() const {
123 0 : bool output_static_targetdist_files=true;
124 0 : if(vesbias_pntr_!=NULL) {
125 0 : output_static_targetdist_files = vesbias_pntr_->isStaticTargetDistFileOutputActive();
126 : }
127 0 : return output_static_targetdist_files;
128 : }
129 :
130 :
131 81 : void LinearBasisSetExpansion::linkVesBias(VesBias* vesbias_pntr_in) {
132 81 : vesbias_pntr_ = vesbias_pntr_in;
133 81 : action_pntr_ = static_cast<Action*>(vesbias_pntr_in);
134 :
135 81 : }
136 :
137 :
138 0 : void LinearBasisSetExpansion::linkAction(Action* action_pntr_in) {
139 0 : action_pntr_ = action_pntr_in;
140 0 : }
141 :
142 :
143 118 : void LinearBasisSetExpansion::setGridBins(const std::vector<unsigned int>& grid_bins_in) {
144 118 : plumed_massert(grid_bins_in.size()==nargs_,"the number of grid bins given doesn't match the number of arguments");
145 118 : plumed_massert(bias_grid_pntr_==NULL,"setGridBins should be used before setting up the grids, otherwise it doesn't work");
146 118 : plumed_massert(fes_grid_pntr_==NULL,"setGridBins should be used before setting up the grids, otherwise it doesn't work");
147 118 : grid_bins_=grid_bins_in;
148 118 : }
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 200 : Grid* LinearBasisSetExpansion::setupGeneralGrid(const std::string& label_suffix, const bool usederiv) {
158 200 : bool use_spline = false;
159 200 : Grid* grid_pntr = new Grid(label_+"."+label_suffix,args_pntrs_,grid_min_,grid_max_,grid_bins_,use_spline,usederiv);
160 200 : return grid_pntr;
161 : }
162 :
163 :
164 153 : void LinearBasisSetExpansion::setupBiasGrid(const bool usederiv) {
165 306 : if(bias_grid_pntr_!=NULL) {return;}
166 118 : bias_grid_pntr_ = setupGeneralGrid("bias",usederiv);
167 118 : if(biasCutoffActive()) {
168 3 : bias_withoutcutoff_grid_pntr_ = setupGeneralGrid("bias_withoutcutoff",usederiv);
169 : }
170 : }
171 :
172 :
173 109 : void LinearBasisSetExpansion::setupFesGrid() {
174 218 : if(fes_grid_pntr_!=NULL) {return;}
175 79 : if(bias_grid_pntr_==NULL) {
176 31 : setupBiasGrid(true);
177 : }
178 79 : fes_grid_pntr_ = setupGeneralGrid("fes",false);
179 : }
180 :
181 :
182 7 : void LinearBasisSetExpansion::setupFesProjGrid() {
183 7 : if(fes_grid_pntr_==NULL) {
184 1 : setupFesGrid();
185 : }
186 7 : }
187 :
188 :
189 764 : void LinearBasisSetExpansion::updateBiasGrid() {
190 764 : plumed_massert(bias_grid_pntr_!=NULL,"the bias grid is not defined");
191 764 : if(action_pntr_!=NULL && getStepOfLastBiasGridUpdate()==action_pntr_->getStep()) {
192 984 : return;
193 : }
194 1783151 : for(Grid::index_t l=0; l<bias_grid_pntr_->getSize(); l++) {
195 1782607 : std::vector<double> forces(nargs_);
196 3565214 : std::vector<double> args = bias_grid_pntr_->getPoint(l);
197 1782607 : bool all_inside=true;
198 1782607 : double bias=getBiasAndForces(args,all_inside,forces);
199 : //
200 1782607 : if(biasCutoffActive()) {
201 600 : vesbias_pntr_->applyBiasCutoff(bias,forces);
202 : }
203 : //
204 1782607 : if(bias_grid_pntr_->hasDerivatives()) {
205 1344196 : bias_grid_pntr_->setValueAndDerivatives(l,bias,forces);
206 : }
207 : else {
208 438411 : bias_grid_pntr_->setValue(l,bias);
209 : }
210 : //
211 1782607 : }
212 544 : if(vesbias_pntr_!=NULL) {
213 433 : vesbias_pntr_->setCurrentBiasMaxValue(bias_grid_pntr_->getMaxValue());
214 : }
215 544 : if(action_pntr_!=NULL) {
216 433 : 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 27 : if(action_pntr_!=NULL && getStepOfLastBiasWithoutCutoffGridUpdate()==action_pntr_->getStep()) {
226 31 : return;
227 : }
228 : //
229 2423 : for(Grid::index_t l=0; l<bias_withoutcutoff_grid_pntr_->getSize(); l++) {
230 2400 : std::vector<double> forces(nargs_);
231 4800 : 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 2400 : 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 2400 : }
241 : //
242 23 : double bias_max = bias_withoutcutoff_grid_pntr_->getMaxValue();
243 23 : double bias_min = bias_withoutcutoff_grid_pntr_->getMinValue();
244 23 : double shift = 0.0;
245 23 : bool bias_shifted=false;
246 23 : if(bias_min < 0.0) {
247 22 : shift += -bias_min;
248 22 : bias_shifted=true;
249 22 : BiasCoeffs()[0] -= bias_min;
250 22 : bias_max -= bias_min;
251 : }
252 23 : if(bias_max > vesbias_pntr_->getBiasCutoffValue()) {
253 22 : shift += -(bias_max-vesbias_pntr_->getBiasCutoffValue());
254 22 : bias_shifted=true;
255 22 : BiasCoeffs()[0] -= (bias_max-vesbias_pntr_->getBiasCutoffValue());
256 22 : 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 2322 : for(Grid::index_t l=0; l<bias_withoutcutoff_grid_pntr_->getSize(); l++) {
262 2300 : 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 23 : vesbias_pntr_->setCurrentBiasMaxValue(bias_max);
273 : }
274 23 : if(action_pntr_!=NULL) {
275 23 : setStepOfLastBiasWithoutCutoffGridUpdate(action_pntr_->getStep());
276 : }
277 : }
278 :
279 :
280 491 : void LinearBasisSetExpansion::updateFesGrid() {
281 491 : plumed_massert(fes_grid_pntr_!=NULL,"the FES grid is not defined");
282 491 : updateBiasGrid();
283 491 : if(action_pntr_!=NULL && getStepOfLastFesGridUpdate() == action_pntr_->getStep()) {
284 551 : return;
285 : }
286 : //
287 431 : double bias2fes_scalingf = -1.0;
288 1344327 : for(Grid::index_t l=0; l<fes_grid_pntr_->getSize(); l++) {
289 1343896 : double fes_value = bias2fes_scalingf*bias_grid_pntr_->getValue(l);
290 1343896 : if(log_targetdist_grid_pntr_!=NULL) {
291 1094866 : fes_value += kBT()*log_targetdist_grid_pntr_->getValue(l);
292 : }
293 1343896 : fes_grid_pntr_->setValue(l,fes_value);
294 : }
295 431 : fes_grid_pntr_->setMinToZero();
296 431 : if(action_pntr_!=NULL) {
297 431 : setStepOfLastFesGridUpdate(action_pntr_->getStep());
298 : }
299 : }
300 :
301 :
302 199 : void LinearBasisSetExpansion::writeBiasGridToFile(OFile& ofile, const bool append_file) const {
303 199 : plumed_massert(bias_grid_pntr_!=NULL,"the bias grid is not defined");
304 199 : if(append_file) {ofile.enforceRestart();}
305 199 : bias_grid_pntr_->writeToFile(ofile);
306 199 : }
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 158 : void LinearBasisSetExpansion::writeFesGridToFile(OFile& ofile, const bool append_file) const {
317 158 : plumed_massert(fes_grid_pntr_!=NULL,"the FES grid is not defined");
318 158 : if(append_file) {ofile.enforceRestart();}
319 158 : fes_grid_pntr_->writeToFile(ofile);
320 158 : }
321 :
322 :
323 30 : void LinearBasisSetExpansion::writeFesProjGridToFile(const std::vector<std::string>& proj_arg, OFile& ofile, const bool append_file) const {
324 30 : plumed_massert(fes_grid_pntr_!=NULL,"the FES grid is not defined");
325 30 : FesWeight* Fw = new FesWeight(beta_);
326 30 : Grid proj_grid = fes_grid_pntr_->project(proj_arg,Fw);
327 30 : proj_grid.setMinToZero();
328 30 : if(append_file) {ofile.enforceRestart();}
329 30 : proj_grid.writeToFile(ofile);
330 30 : delete Fw;
331 30 : }
332 :
333 :
334 69 : void LinearBasisSetExpansion::writeTargetDistGridToFile(OFile& ofile, const bool append_file) const {
335 138 : if(targetdist_grid_pntr_==NULL) {return;}
336 69 : if(append_file) {ofile.enforceRestart();}
337 69 : targetdist_grid_pntr_->writeToFile(ofile);
338 : }
339 :
340 :
341 69 : void LinearBasisSetExpansion::writeLogTargetDistGridToFile(OFile& ofile, const bool append_file) const {
342 138 : if(log_targetdist_grid_pntr_==NULL) {return;}
343 69 : if(append_file) {ofile.enforceRestart();}
344 69 : log_targetdist_grid_pntr_->writeToFile(ofile);
345 : }
346 :
347 :
348 14 : void LinearBasisSetExpansion::writeTargetDistProjGridToFile(const std::vector<std::string>& proj_arg, OFile& ofile, const bool append_file) const {
349 28 : if(targetdist_grid_pntr_==NULL) {return;}
350 14 : if(append_file) {ofile.enforceRestart();}
351 14 : Grid proj_grid = TargetDistribution::getMarginalDistributionGrid(targetdist_grid_pntr_,proj_arg);
352 14 : 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 1790385 : 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 1790385 : unsigned int nargs = args_values.size();
372 1790385 : plumed_assert(coeffs_pntr_in->numberOfDimensions()==nargs);
373 1790385 : plumed_assert(basisf_pntrs_in.size()==nargs);
374 1790385 : plumed_assert(forces.size()==nargs);
375 1790385 : plumed_assert(coeffsderivs_values.size()==coeffs_pntr_in->numberOfCoeffs());
376 :
377 1790385 : std::vector<double> args_values_trsfrm(nargs);
378 : // std::vector<bool> inside_interval(nargs,true);
379 1790385 : all_inside = true;
380 : //
381 3580770 : std::vector< std::vector <double> > bf_values(nargs);
382 3580770 : std::vector< std::vector <double> > bf_derivs(nargs);
383 : //
384 5305597 : for(unsigned int k=0; k<nargs; k++) {
385 3515212 : bf_values[k].assign(basisf_pntrs_in[k]->getNumberOfBasisFunctions(),0.0);
386 3515212 : bf_derivs[k].assign(basisf_pntrs_in[k]->getNumberOfBasisFunctions(),0.0);
387 3515212 : bool curr_inside=true;
388 3515212 : 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 3515212 : if(!curr_inside) {all_inside=false;}
391 3515212 : forces[k]=0.0;
392 : }
393 : //
394 1790385 : size_t stride=1;
395 1790385 : size_t rank=0;
396 1790385 : if(comm_in!=NULL)
397 : {
398 1790385 : stride=comm_in->Get_size();
399 1790385 : rank=comm_in->Get_rank();
400 : }
401 : // loop over coeffs
402 1790385 : double bias=0.0;
403 122613111 : for(size_t i=rank; i<coeffs_pntr_in->numberOfCoeffs(); i+=stride) {
404 120822726 : std::vector<unsigned int> indices=coeffs_pntr_in->getIndices(i);
405 120822726 : double coeff = coeffs_pntr_in->getValue(i);
406 120822726 : double bf_curr=1.0;
407 361939976 : for(unsigned int k=0; k<nargs; k++) {
408 241117250 : bf_curr*=bf_values[k][indices[k]];
409 : }
410 120822726 : bias+=coeff*bf_curr;
411 120822726 : coeffsderivs_values[i] = bf_curr;
412 361939976 : for(unsigned int k=0; k<nargs; k++) {
413 241117250 : double der = 1.0;
414 722823548 : for(unsigned int l=0; l<nargs; l++) {
415 481706298 : if(l!=k) {der*=bf_values[l][indices[l]];}
416 241117250 : else {der*=bf_derivs[l][indices[l]];}
417 : }
418 241117250 : 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 120822726 : }
423 : //
424 1790385 : if(comm_in!=NULL) {
425 : // coeffsderivs_values is not summed as the mpi Sum is done later on for the averages
426 1790385 : comm_in->Sum(bias);
427 1790385 : comm_in->Sum(forces);
428 : }
429 3580770 : return bias;
430 : }
431 :
432 :
433 718566 : 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 718566 : unsigned int nargs = args_values.size();
435 718566 : plumed_assert(coeffs_pntr_in->numberOfDimensions()==nargs);
436 718566 : plumed_assert(basisf_pntrs_in.size()==nargs);
437 :
438 718566 : std::vector<double> args_values_trsfrm(nargs);
439 1437132 : std::vector< std::vector <double> > bf_values;
440 : //
441 2135765 : for(unsigned int k=0; k<nargs; k++) {
442 1417199 : std::vector<double> tmp_val(basisf_pntrs_in[k]->getNumberOfBasisFunctions());
443 2834398 : std::vector<double> tmp_der(tmp_val.size());
444 1417199 : bool inside=true;
445 1417199 : basisf_pntrs_in[k]->getAllValues(args_values[k],args_values_trsfrm[k],inside,tmp_val,tmp_der);
446 1417199 : bf_values.push_back(tmp_val);
447 1417199 : }
448 : //
449 718566 : size_t stride=1;
450 718566 : size_t rank=0;
451 718566 : 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 70103296 : for(size_t i=rank; i<coeffs_pntr_in->numberOfCoeffs(); i+=stride) {
458 69384730 : std::vector<unsigned int> indices=coeffs_pntr_in->getIndices(i);
459 69384730 : double bf_curr=1.0;
460 207934127 : for(unsigned int k=0; k<nargs; k++) {
461 138549397 : bf_curr*=bf_values[k][indices[k]];
462 : }
463 69384730 : basisset_values[i] = bf_curr;
464 69384730 : }
465 : //
466 718566 : if(comm_in!=NULL) {
467 0 : comm_in->Sum(basisset_values);
468 718566 : }
469 718566 : }
470 :
471 :
472 368 : double LinearBasisSetExpansion::getBasisSetValue(const std::vector<double>& args_values, const size_t index, std::vector<BasisFunctions*>& basisf_pntrs_in, CoeffsVector* coeffs_pntr_in) {
473 368 : unsigned int nargs = args_values.size();
474 368 : plumed_assert(coeffs_pntr_in->numberOfDimensions()==nargs);
475 368 : plumed_assert(basisf_pntrs_in.size()==nargs);
476 :
477 368 : std::vector<double> args_values_trsfrm(nargs);
478 736 : std::vector< std::vector <double> > bf_values;
479 : //
480 838 : for(unsigned int k=0; k<nargs; k++) {
481 470 : std::vector<double> tmp_val(basisf_pntrs_in[k]->getNumberOfBasisFunctions());
482 940 : std::vector<double> tmp_der(tmp_val.size());
483 470 : bool inside=true;
484 470 : basisf_pntrs_in[k]->getAllValues(args_values[k],args_values_trsfrm[k],inside,tmp_val,tmp_der);
485 470 : bf_values.push_back(tmp_val);
486 470 : }
487 : //
488 736 : std::vector<unsigned int> indices=coeffs_pntr_in->getIndices(index);
489 368 : double bf_value=1.0;
490 838 : for(unsigned int k=0; k<nargs; k++) {
491 470 : bf_value*=bf_values[k][indices[k]];
492 : }
493 736 : return bf_value;
494 : }
495 :
496 :
497 42 : void LinearBasisSetExpansion::setupUniformTargetDistribution() {
498 42 : std::vector< std::vector <double> > bf_integrals(0);
499 84 : std::vector<double> targetdist_averages(ncoeffs_,0.0);
500 : //
501 94 : for(unsigned int k=0; k<nargs_; k++) {
502 52 : bf_integrals.push_back(basisf_pntrs_[k]->getUniformIntegrals());
503 : }
504 : //
505 1636 : for(size_t i=0; i<ncoeffs_; i++) {
506 1594 : std::vector<unsigned int> indices=bias_coeffs_pntr_->getIndices(i);
507 1594 : double value = 1.0;
508 4426 : for(unsigned int k=0; k<nargs_; k++) {
509 2832 : value*=bf_integrals[k][indices[k]];
510 : }
511 1594 : targetdist_averages[i]=value;
512 1594 : }
513 84 : TargetDistAverages() = targetdist_averages;
514 42 : }
515 :
516 :
517 39 : void LinearBasisSetExpansion::setupTargetDistribution(TargetDistribution* targetdist_pntr_in) {
518 39 : targetdist_pntr_ = targetdist_pntr_in;
519 : //
520 39 : targetdist_pntr_->setupGrids(args_pntrs_,grid_min_,grid_max_,grid_bins_);
521 39 : targetdist_grid_pntr_ = targetdist_pntr_->getTargetDistGridPntr();
522 39 : log_targetdist_grid_pntr_ = targetdist_pntr_->getLogTargetDistGridPntr();
523 : //
524 39 : if(targetdist_pntr_->isDynamic()) {
525 33 : vesbias_pntr_->enableDynamicTargetDistribution();
526 : }
527 : //
528 39 : if(targetdist_pntr_->biasGridNeeded()) {
529 0 : setupBiasGrid(true);
530 0 : targetdist_pntr_->linkBiasGrid(bias_grid_pntr_);
531 : }
532 39 : if(targetdist_pntr_->biasWithoutCutoffGridNeeded()) {
533 3 : setupBiasGrid(true);
534 3 : targetdist_pntr_->linkBiasWithoutCutoffGrid(bias_withoutcutoff_grid_pntr_);
535 : }
536 39 : if(targetdist_pntr_->fesGridNeeded()) {
537 30 : setupFesGrid();
538 30 : targetdist_pntr_->linkFesGrid(fes_grid_pntr_);
539 : }
540 : //
541 39 : targetdist_pntr_->updateTargetDist();
542 39 : calculateTargetDistAveragesFromGrid(targetdist_grid_pntr_);
543 39 : }
544 :
545 :
546 321 : void LinearBasisSetExpansion::updateTargetDistribution() {
547 321 : plumed_massert(targetdist_pntr_!=NULL,"the target distribution hasn't been setup!");
548 321 : plumed_massert(targetdist_pntr_->isDynamic(),"this should only be used for dynamically updated target distributions!");
549 321 : if(targetdist_pntr_->biasGridNeeded()) {updateBiasGrid();}
550 321 : if(biasCutoffActive()) {updateBiasWithoutCutoffGrid();}
551 321 : if(targetdist_pntr_->fesGridNeeded()) {updateFesGrid();}
552 321 : targetdist_pntr_->updateTargetDist();
553 321 : calculateTargetDistAveragesFromGrid(targetdist_grid_pntr_);
554 321 : }
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 368 : void LinearBasisSetExpansion::calculateTargetDistAveragesFromGrid(const Grid* targetdist_grid_pntr) {
572 368 : plumed_assert(targetdist_grid_pntr!=NULL);
573 368 : std::vector<double> targetdist_averages(ncoeffs_,0.0);
574 736 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(targetdist_grid_pntr);
575 368 : Grid::index_t stride=mycomm_.Get_size();
576 368 : Grid::index_t rank=mycomm_.Get_rank();
577 718934 : for(Grid::index_t l=rank; l<targetdist_grid_pntr->getSize(); l+=stride) {
578 718566 : std::vector<double> args_values = targetdist_grid_pntr->getPoint(l);
579 1437132 : std::vector<double> basisset_values(ncoeffs_);
580 : // parallelization done over the grid -> should NOT use parallel in getBasisSetValues!!
581 718566 : getBasisSetValues(args_values,basisset_values,false);
582 718566 : double weight = integration_weights[l]*targetdist_grid_pntr->getValue(l);
583 70103296 : for(unsigned int i=0; i<ncoeffs_; i++) {
584 69384730 : targetdist_averages[i] += weight*basisset_values[i];
585 : }
586 718566 : }
587 368 : mycomm_.Sum(targetdist_averages);
588 : // the overall constant;
589 368 : targetdist_averages[0] = getBasisSetConstant();
590 736 : TargetDistAverages() = targetdist_averages;
591 368 : }
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 37 : 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 1783089 : bool LinearBasisSetExpansion::biasCutoffActive() const {
609 1783089 : if(vesbias_pntr_!=NULL) {return vesbias_pntr_->biasCutoffActive();}
610 438448 : 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 0 : 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 4821 : }
|
