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 "VesBias.h"
24 : #include "BasisFunctions.h"
25 : #include "CoeffsVector.h"
26 : #include "CoeffsMatrix.h"
27 : #include "Optimizer.h"
28 : #include "FermiSwitchingFunction.h"
29 : #include "VesTools.h"
30 : #include "TargetDistribution.h"
31 :
32 : #include "tools/Communicator.h"
33 : #include "core/ActionSet.h"
34 : #include "core/PlumedMain.h"
35 : #include "tools/File.h"
36 :
37 :
38 : namespace PLMD {
39 : namespace ves {
40 :
41 90 : VesBias::VesBias(const ActionOptions&ao):
42 : Action(ao),
43 : Bias(ao),
44 90 : ncoeffssets_(0),
45 90 : sampled_averages(0),
46 90 : sampled_cross_averages(0),
47 90 : use_multiple_coeffssets_(false),
48 90 : coeffs_fnames(0),
49 90 : ncoeffs_total_(0),
50 90 : optimizer_pntr_(NULL),
51 90 : optimize_coeffs_(false),
52 90 : compute_hessian_(false),
53 90 : diagonal_hessian_(true),
54 90 : aver_counters(0),
55 90 : kbt_(0.0),
56 90 : targetdist_pntrs_(0),
57 90 : dynamic_targetdist_(false),
58 90 : grid_bins_(0),
59 90 : grid_min_(0),
60 90 : grid_max_(0),
61 90 : bias_filename_(""),
62 90 : fes_filename_(""),
63 90 : targetdist_filename_(""),
64 90 : coeffs_id_prefix_("c-"),
65 90 : bias_file_fmt_("14.9f"),
66 90 : fes_file_fmt_("14.9f"),
67 90 : targetdist_file_fmt_("14.9f"),
68 90 : targetdist_restart_file_fmt_("30.16e"),
69 90 : filenames_have_iteration_number_(false),
70 90 : bias_fileoutput_active_(false),
71 90 : fes_fileoutput_active_(false),
72 90 : fesproj_fileoutput_active_(false),
73 90 : dynamic_targetdist_fileoutput_active_(false),
74 90 : static_targetdist_fileoutput_active_(true),
75 90 : bias_cutoff_active_(false),
76 90 : bias_cutoff_value_(0.0),
77 90 : bias_current_max_value(0.0),
78 90 : calc_reweightfactor_(false),
79 270 : optimization_threshold_(0.0)
80 : {
81 90 : log.printf(" VES bias, please read and cite ");
82 180 : log << plumed.cite("Valsson and Parrinello, Phys. Rev. Lett. 113, 090601 (2014)");
83 90 : log.printf("\n");
84 :
85 90 : kbt_=getkBT();
86 90 : if(kbt_>0.0) {
87 90 : log.printf(" KbT: %f\n",kbt_);
88 : }
89 : // NOTE: the check for that the temperature is given is done when linking the optimizer later on.
90 :
91 180 : if(keywords.exists("COEFFS")) {
92 180 : parseVector("COEFFS",coeffs_fnames);
93 : }
94 :
95 180 : if(keywords.exists("GRID_BINS")) {
96 180 : parseMultipleValues<unsigned int>("GRID_BINS",grid_bins_,getNumberOfArguments(),100);
97 : }
98 :
99 90 : if(keywords.exists("GRID_MIN") && keywords.exists("GRID_MAX")) {
100 0 : parseMultipleValues("GRID_MIN",grid_min_,getNumberOfArguments());
101 0 : parseMultipleValues("GRID_MAX",grid_max_,getNumberOfArguments());
102 : }
103 :
104 : std::vector<std::string> targetdist_labels;
105 180 : if(keywords.exists("TARGET_DISTRIBUTION")) {
106 180 : parseVector("TARGET_DISTRIBUTION",targetdist_labels);
107 90 : if(targetdist_labels.size()>1) {
108 0 : plumed_merror(getName()+" with label "+getLabel()+": multiple target distribution labels not allowed");
109 : }
110 : }
111 0 : else if(keywords.exists("TARGET_DISTRIBUTIONS")) {
112 0 : parseVector("TARGET_DISTRIBUTIONS",targetdist_labels);
113 : }
114 :
115 90 : std::string error_msg = "";
116 180 : targetdist_pntrs_ = VesTools::getPointersFromLabels<TargetDistribution*>(targetdist_labels,plumed.getActionSet(),error_msg);
117 90 : if(error_msg.size()>0) {plumed_merror("Problem with target distribution in "+getName()+": "+error_msg);}
118 :
119 135 : for(unsigned int i=0; i<targetdist_pntrs_.size(); i++) {
120 45 : targetdist_pntrs_[i]->linkVesBias(this);
121 : }
122 :
123 :
124 90 : if(getNumberOfArguments()>2) {
125 : disableStaticTargetDistFileOutput();
126 : }
127 :
128 :
129 180 : if(keywords.exists("BIAS_FILE")) {
130 180 : parse("BIAS_FILE",bias_filename_);
131 90 : if(bias_filename_.size()==0) {
132 180 : bias_filename_ = "bias." + getLabel() + ".data";
133 : }
134 : }
135 180 : if(keywords.exists("FES_FILE")) {
136 180 : parse("FES_FILE",fes_filename_);
137 90 : if(fes_filename_.size()==0) {
138 180 : fes_filename_ = "fes." + getLabel() + ".data";
139 : }
140 : }
141 180 : if(keywords.exists("TARGETDIST_FILE")) {
142 180 : parse("TARGETDIST_FILE",targetdist_filename_);
143 90 : if(targetdist_filename_.size()==0) {
144 180 : targetdist_filename_ = "targetdist." + getLabel() + ".data";
145 : }
146 : }
147 : //
148 180 : if(keywords.exists("BIAS_FILE_FMT")) {
149 0 : parse("BIAS_FILE_FMT",bias_file_fmt_);
150 : }
151 180 : if(keywords.exists("FES_FILE_FMT")) {
152 0 : parse("FES_FILE_FMT",fes_file_fmt_);
153 : }
154 180 : if(keywords.exists("TARGETDIST_FILE_FMT")) {
155 0 : parse("TARGETDIST_FILE_FMT",targetdist_file_fmt_);
156 : }
157 180 : if(keywords.exists("TARGETDIST_RESTART_FILE_FMT")) {
158 0 : parse("TARGETDIST_RESTART_FILE_FMT",targetdist_restart_file_fmt_);
159 : }
160 :
161 : //
162 180 : if(keywords.exists("BIAS_CUTOFF")) {
163 90 : double cutoff_value=0.0;
164 90 : parse("BIAS_CUTOFF",cutoff_value);
165 90 : if(cutoff_value<0.0) {
166 0 : plumed_merror("the value given in BIAS_CUTOFF doesn't make sense, it should be larger than 0.0");
167 : }
168 : //
169 90 : if(cutoff_value>0.0) {
170 3 : double fermi_lambda=10.0;
171 3 : parse("BIAS_CUTOFF_FERMI_LAMBDA",fermi_lambda);
172 3 : setupBiasCutoff(cutoff_value,fermi_lambda);
173 3 : log.printf(" Employing a bias cutoff of %f (the lambda value for the Fermi switching function is %f), see and cite ",cutoff_value,fermi_lambda);
174 6 : log << plumed.cite("McCarty, Valsson, Tiwary, and Parrinello, Phys. Rev. Lett. 115, 070601 (2015)");
175 3 : log.printf("\n");
176 : }
177 : }
178 :
179 :
180 180 : if(keywords.exists("PROJ_ARG")) {
181 : std::vector<std::string> proj_arg;
182 90 : for(int i=1;; i++) {
183 212 : if(!parseNumberedVector("PROJ_ARG",i,proj_arg)) {break;}
184 : // checks
185 16 : if(proj_arg.size() > (getNumberOfArguments()-1) ) {
186 0 : plumed_merror("PROJ_ARG must be a subset of ARG");
187 : }
188 : //
189 32 : for(unsigned int k=0; k<proj_arg.size(); k++) {
190 : bool found = false;
191 24 : for(unsigned int l=0; l<getNumberOfArguments(); l++) {
192 24 : if(proj_arg[k]==getPntrToArgument(l)->getName()) {
193 : found = true;
194 : break;
195 : }
196 : }
197 16 : if(!found) {
198 0 : std::string s1; Tools::convert(i,s1);
199 0 : std::string error = "PROJ_ARG" + s1 + ": label " + proj_arg[k] + " is not among the arguments given in ARG";
200 0 : plumed_merror(error);
201 : }
202 : }
203 : //
204 16 : projection_args_.push_back(proj_arg);
205 16 : }
206 90 : }
207 :
208 180 : if(keywords.exists("CALC_REWEIGHT_FACTOR")) {
209 0 : parseFlag("CALC_REWEIGHT_FACTOR",calc_reweightfactor_);
210 0 : if(calc_reweightfactor_) {
211 0 : addComponent("rct"); componentIsNotPeriodic("rct");
212 0 : updateReweightFactor();
213 : }
214 : }
215 :
216 180 : if(keywords.exists("OPTIMIZATION_THRESHOLD")) {
217 90 : parse("OPTIMIZATION_THRESHOLD",optimization_threshold_);
218 90 : if(optimization_threshold_ < 0.0) {
219 0 : plumed_merror("OPTIMIZATION_THRESHOLD should be a postive value");
220 : }
221 90 : if(optimization_threshold_!=0.0) {
222 1 : log.printf(" Employing a threshold value of %e for optimization of localized basis functions.\n",optimization_threshold_);
223 : }
224 : }
225 :
226 90 : }
227 :
228 :
229 90 : VesBias::~VesBias() {
230 180 : }
231 :
232 :
233 92 : void VesBias::registerKeywords( Keywords& keys ) {
234 92 : Bias::registerKeywords(keys);
235 184 : keys.add("optional","TEMP","the system temperature - this is needed if the MD code does not pass the temperature to PLUMED.");
236 : //
237 184 : keys.reserve("optional","COEFFS","read in the coefficients from files.");
238 : //
239 184 : keys.reserve("optional","TARGET_DISTRIBUTION","the label of the target distribution to be used.");
240 184 : keys.reserve("optional","TARGET_DISTRIBUTIONS","the label of the target distribution to be used. Here you are allows to use multiple labels.");
241 : //
242 184 : keys.reserve("optional","GRID_BINS","the number of bins used for the grid. The default value is 100 bins per dimension.");
243 184 : keys.reserve("optional","GRID_MIN","the lower bounds used for the grid.");
244 184 : keys.reserve("optional","GRID_MAX","the upper bounds used for the grid.");
245 : //
246 184 : keys.add("optional","BIAS_FILE","filename of the file on which the bias should be written out. By default it is bias.LABEL.data. Note that suffixes indicating the iteration number (iter-#) are added to the filename when optimizing coefficients.");
247 184 : keys.add("optional","FES_FILE","filename of the file on which the FES should be written out. By default it is fes.LABEL.data. Note that suffixes indicating the iteration number (iter-#) are added to the filename when optimizing coefficients.");
248 184 : keys.add("optional","TARGETDIST_FILE","filename of the file on which the target distribution should be written out. By default it is targetdist.LABEL.data. Note that suffixes indicating the iteration number (iter-#) are added to the filename when optimizing coefficients and the target distribution is dynamic.");
249 : //
250 : // keys.add("optional","BIAS_FILE_FMT","the format of the bias files, by default it is %14.9f.");
251 : // keys.add("optional","FES_FILE_FMT","the format of the FES files, by default it is %14.9f.");
252 : // keys.add("optional","TARGETDIST_FILE_FMT","the format of the target distribution files, by default it is %14.9f.");
253 : // keys.add("hidden","TARGETDIST_RESTART_FILE_FMT","the format of the target distribution files that are used for restarting, by default it is %30.16e.");
254 : //
255 184 : keys.reserve("optional","BIAS_CUTOFF","cutoff the bias such that it only fills the free energy surface up to certain level F_cutoff, here you should give the value of the F_cutoff.");
256 184 : keys.reserve("optional","BIAS_CUTOFF_FERMI_LAMBDA","the lambda value used in the Fermi switching function for the bias cutoff (BIAS_CUTOFF), the default value is 10.0.");
257 : //
258 184 : keys.reserve("numbered","PROJ_ARG","arguments for doing projections of the FES or the target distribution.");
259 : //
260 184 : keys.reserveFlag("CALC_REWEIGHT_FACTOR",false,"enable the calculation of the reweight factor c(t). You should also give a stride for updating the reweight factor in the optimizer by using the REWEIGHT_FACTOR_STRIDE keyword if the coefficients are updated.");
261 184 : keys.add("optional","OPTIMIZATION_THRESHOLD","Threshold value to turn off optimization of localized basis functions.");
262 :
263 92 : }
264 :
265 :
266 92 : void VesBias::useInitialCoeffsKeywords(Keywords& keys) {
267 92 : keys.use("COEFFS");
268 92 : }
269 :
270 :
271 92 : void VesBias::useTargetDistributionKeywords(Keywords& keys) {
272 184 : plumed_massert(!keys.exists("TARGET_DISTRIBUTIONS"),"you cannot use both useTargetDistributionKeywords and useMultipleTargetDistributionKeywords");
273 92 : keys.use("TARGET_DISTRIBUTION");
274 92 : }
275 :
276 :
277 0 : void VesBias::useMultipleTargetDistributionKeywords(Keywords& keys) {
278 0 : plumed_massert(!keys.exists("TARGET_DISTRIBUTION"),"you cannot use both useTargetDistributionKeywords and useMultipleTargetDistributionKeywords");
279 0 : keys.use("TARGET_DISTRIBUTIONS");
280 0 : }
281 :
282 :
283 92 : void VesBias::useGridBinKeywords(Keywords& keys) {
284 92 : keys.use("GRID_BINS");
285 92 : }
286 :
287 :
288 0 : void VesBias::useGridLimitsKeywords(Keywords& keys) {
289 0 : keys.use("GRID_MIN");
290 0 : keys.use("GRID_MAX");
291 0 : }
292 :
293 :
294 92 : void VesBias::useBiasCutoffKeywords(Keywords& keys) {
295 92 : keys.use("BIAS_CUTOFF");
296 92 : keys.use("BIAS_CUTOFF_FERMI_LAMBDA");
297 92 : }
298 :
299 :
300 92 : void VesBias::useProjectionArgKeywords(Keywords& keys) {
301 92 : keys.use("PROJ_ARG");
302 92 : }
303 :
304 :
305 0 : void VesBias::useReweightFactorKeywords(Keywords& keys) {
306 0 : keys.use("CALC_REWEIGHT_FACTOR");
307 0 : keys.addOutputComponent("rct","CALC_REWEIGHT_FACTOR","the reweight factor c(t).");
308 0 : }
309 :
310 :
311 0 : void VesBias::addCoeffsSet(const std::vector<std::string>& dimension_labels,const std::vector<unsigned int>& indices_shape) {
312 0 : auto coeffs_pntr_tmp = Tools::make_unique<CoeffsVector>("coeffs",dimension_labels,indices_shape,comm,true);
313 0 : initializeCoeffs(std::move(coeffs_pntr_tmp));
314 0 : }
315 :
316 :
317 90 : void VesBias::addCoeffsSet(std::vector<Value*>& args,std::vector<BasisFunctions*>& basisf) {
318 90 : auto coeffs_pntr_tmp = Tools::make_unique<CoeffsVector>("coeffs",args,basisf,comm,true);
319 90 : initializeCoeffs(std::move(coeffs_pntr_tmp));
320 90 : }
321 :
322 :
323 0 : void VesBias::addCoeffsSet(std::unique_ptr<CoeffsVector> coeffs_pntr_in) {
324 0 : initializeCoeffs(std::move(coeffs_pntr_in));
325 0 : }
326 :
327 :
328 90 : void VesBias::initializeCoeffs(std::unique_ptr<CoeffsVector> coeffs_pntr_in) {
329 : //
330 90 : coeffs_pntr_in->linkVesBias(this);
331 : //
332 : std::string label;
333 90 : if(!use_multiple_coeffssets_ && ncoeffssets_==1) {
334 0 : plumed_merror("you are not allowed to use multiple coefficient sets");
335 : }
336 : //
337 180 : label = getCoeffsSetLabelString("coeffs",ncoeffssets_);
338 90 : coeffs_pntr_in->setLabels(label);
339 :
340 90 : coeffs_pntrs_.emplace_back(std::move(coeffs_pntr_in));
341 : auto aver_ps_tmp = Tools::make_unique<CoeffsVector>(*coeffs_pntrs_.back());
342 180 : label = getCoeffsSetLabelString("targetdist_averages",ncoeffssets_);
343 90 : aver_ps_tmp->setLabels(label);
344 90 : aver_ps_tmp->setValues(0.0);
345 90 : targetdist_averages_pntrs_.emplace_back(std::move(aver_ps_tmp));
346 : //
347 : auto gradient_tmp = Tools::make_unique<CoeffsVector>(*coeffs_pntrs_.back());
348 180 : label = getCoeffsSetLabelString("gradient",ncoeffssets_);
349 90 : gradient_tmp->setLabels(label);
350 90 : gradient_pntrs_.emplace_back(std::move(gradient_tmp));
351 : //
352 180 : label = getCoeffsSetLabelString("hessian",ncoeffssets_);
353 90 : auto hessian_tmp = Tools::make_unique<CoeffsMatrix>(label,coeffs_pntrs_.back().get(),comm,diagonal_hessian_);
354 :
355 90 : hessian_pntrs_.emplace_back(std::move(hessian_tmp));
356 : //
357 : std::vector<double> aver_sampled_tmp;
358 90 : aver_sampled_tmp.assign(coeffs_pntrs_.back()->numberOfCoeffs(),0.0);
359 90 : sampled_averages.push_back(aver_sampled_tmp);
360 : //
361 : std::vector<double> cross_aver_sampled_tmp;
362 90 : cross_aver_sampled_tmp.assign(hessian_pntrs_.back()->getSize(),0.0);
363 90 : sampled_cross_averages.push_back(cross_aver_sampled_tmp);
364 : //
365 90 : aver_counters.push_back(0);
366 : //
367 90 : ncoeffssets_++;
368 180 : }
369 :
370 :
371 90 : bool VesBias::readCoeffsFromFiles() {
372 90 : plumed_assert(ncoeffssets_>0);
373 180 : plumed_massert(keywords.exists("COEFFS"),"you are not allowed to use this function as the COEFFS keyword is not enabled");
374 : bool read_coeffs = false;
375 90 : if(coeffs_fnames.size()>0) {
376 4 : plumed_massert(coeffs_fnames.size()==ncoeffssets_,"COEFFS keyword is of the wrong size");
377 4 : if(ncoeffssets_==1) {
378 4 : log.printf(" Read in coefficients from file ");
379 : }
380 : else {
381 0 : log.printf(" Read in coefficients from files:\n");
382 : }
383 8 : for(unsigned int i=0; i<ncoeffssets_; i++) {
384 4 : IFile ifile;
385 4 : ifile.link(*this);
386 4 : ifile.open(coeffs_fnames[i]);
387 8 : if(!ifile.FieldExist(coeffs_pntrs_[i]->getDataLabel())) {
388 0 : std::string error_msg = "Problem with reading coefficients from file " + ifile.getPath() + ": no field with name " + coeffs_pntrs_[i]->getDataLabel() + "\n";
389 0 : plumed_merror(error_msg);
390 : }
391 4 : size_t ncoeffs_read = coeffs_pntrs_[i]->readFromFile(ifile,false,false);
392 4 : coeffs_pntrs_[i]->setIterationCounterAndTime(0,getTime());
393 4 : if(ncoeffssets_==1) {
394 8 : log.printf("%s (read %zu of %zu values)\n", ifile.getPath().c_str(),ncoeffs_read,coeffs_pntrs_[i]->numberOfCoeffs());
395 : }
396 : else {
397 0 : log.printf(" coefficient %u: %s (read %zu of %zu values)\n",i,ifile.getPath().c_str(),ncoeffs_read,coeffs_pntrs_[i]->numberOfCoeffs());
398 : }
399 4 : ifile.close();
400 4 : }
401 : read_coeffs = true;
402 : }
403 90 : return read_coeffs;
404 : }
405 :
406 :
407 22810 : void VesBias::updateGradientAndHessian(const bool use_mwalkers_mpi) {
408 45620 : for(unsigned int k=0; k<ncoeffssets_; k++) {
409 : //
410 22810 : comm.Sum(sampled_averages[k]);
411 22810 : comm.Sum(sampled_cross_averages[k]);
412 22810 : if(use_mwalkers_mpi) {
413 : double walker_weight=1.0;
414 120 : if(aver_counters[k]==0) {walker_weight=0.0;}
415 120 : multiSimSumAverages(k,walker_weight);
416 : }
417 : // NOTE: this assumes that all walkers have the same TargetDist, might change later on!!
418 22810 : Gradient(k).setValues( TargetDistAverages(k) - sampled_averages[k] );
419 22810 : Hessian(k) = computeCovarianceFromAverages(k);
420 22810 : Hessian(k) *= getBeta();
421 :
422 22810 : if(optimization_threshold_ != 0.0) {
423 390 : for(size_t c_id=0; c_id < sampled_averages[k].size(); ++c_id) {
424 380 : if(fabs(sampled_averages[k][c_id]) < optimization_threshold_) {
425 219 : Gradient(k).setValue(c_id, 0.0);
426 219 : Hessian(k).setValue(c_id, c_id, 0.0);
427 : }
428 : }
429 : }
430 : //
431 : Gradient(k).activate();
432 : Hessian(k).activate();
433 : //
434 : // Check the total number of samples (from all walkers) and deactivate the Gradient and Hessian if it
435 : // is zero
436 22810 : unsigned int total_samples = aver_counters[k];
437 22810 : if(use_mwalkers_mpi) {
438 120 : if(comm.Get_rank()==0) {multi_sim_comm.Sum(total_samples);}
439 120 : comm.Bcast(total_samples,0);
440 : }
441 22810 : if(total_samples==0) {
442 : Gradient(k).deactivate();
443 95 : Gradient(k).clear();
444 : Hessian(k).deactivate();
445 95 : Hessian(k).clear();
446 : }
447 : //
448 : std::fill(sampled_averages[k].begin(), sampled_averages[k].end(), 0.0);
449 : std::fill(sampled_cross_averages[k].begin(), sampled_cross_averages[k].end(), 0.0);
450 22810 : aver_counters[k]=0;
451 : }
452 22810 : }
453 :
454 :
455 120 : void VesBias::multiSimSumAverages(const unsigned int c_id, const double walker_weight) {
456 120 : plumed_massert(walker_weight>=0.0,"the weight of the walker cannot be negative!");
457 120 : if(walker_weight!=1.0) {
458 7860 : for(size_t i=0; i<sampled_averages[c_id].size(); i++) {
459 7800 : sampled_averages[c_id][i] *= walker_weight;
460 : }
461 7860 : for(size_t i=0; i<sampled_cross_averages[c_id].size(); i++) {
462 7800 : sampled_cross_averages[c_id][i] *= walker_weight;
463 : }
464 : }
465 : //
466 120 : if(comm.Get_rank()==0) {
467 120 : multi_sim_comm.Sum(sampled_averages[c_id]);
468 120 : multi_sim_comm.Sum(sampled_cross_averages[c_id]);
469 120 : double norm_weights = walker_weight;
470 120 : multi_sim_comm.Sum(norm_weights);
471 120 : if(norm_weights>0.0) {norm_weights=1.0/norm_weights;}
472 8580 : for(size_t i=0; i<sampled_averages[c_id].size(); i++) {
473 8460 : sampled_averages[c_id][i] *= norm_weights;
474 : }
475 8580 : for(size_t i=0; i<sampled_cross_averages[c_id].size(); i++) {
476 8460 : sampled_cross_averages[c_id][i] *= norm_weights;
477 : }
478 : }
479 120 : comm.Bcast(sampled_averages[c_id],0);
480 120 : comm.Bcast(sampled_cross_averages[c_id],0);
481 120 : }
482 :
483 :
484 23556 : void VesBias::addToSampledAverages(const std::vector<double>& values, const unsigned int c_id) {
485 : /*
486 : use the following online equation to calculate the average and covariance
487 : (see https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Covariance)
488 : xm[n+1] = xm[n] + (x[n+1]-xm[n])/(n+1)
489 : */
490 23556 : double counter_dbl = static_cast<double>(aver_counters[c_id]);
491 : size_t ncoeffs = numberOfCoeffs(c_id);
492 23556 : std::vector<double> deltas(ncoeffs,0.0);
493 23556 : size_t stride = comm.Get_size();
494 23556 : size_t rank = comm.Get_rank();
495 : // update average and diagonal part of Hessian
496 1830228 : for(size_t i=rank; i<ncoeffs; i+=stride) {
497 : size_t midx = getHessianIndex(i,i,c_id);
498 1806672 : deltas[i] = (values[i]-sampled_averages[c_id][i])/(counter_dbl+1); // (x[n+1]-xm[n])/(n+1)
499 1806672 : sampled_averages[c_id][i] += deltas[i];
500 1806672 : sampled_cross_averages[c_id][midx] += (values[i]*values[i]-sampled_cross_averages[c_id][midx])/(counter_dbl+1);
501 : }
502 23556 : comm.Sum(deltas);
503 : // update off-diagonal part of the Hessian
504 23556 : if(!diagonal_hessian_) {
505 0 : for(size_t i=rank; i<ncoeffs; i+=stride) {
506 0 : for(size_t j=(i+1); j<ncoeffs; j++) {
507 : size_t midx = getHessianIndex(i,j,c_id);
508 0 : sampled_cross_averages[c_id][midx] += (values[i]*values[j]-sampled_cross_averages[c_id][midx])/(counter_dbl+1);
509 : }
510 : }
511 : }
512 : // NOTE: the MPI sum for sampled_averages and sampled_cross_averages is done later
513 23556 : aver_counters[c_id] += 1;
514 23556 : }
515 :
516 :
517 0 : void VesBias::setTargetDistAverages(const std::vector<double>& coeffderivs_aver_ps, const unsigned int coeffs_id) {
518 0 : TargetDistAverages(coeffs_id) = coeffderivs_aver_ps;
519 0 : TargetDistAverages(coeffs_id).setIterationCounterAndTime(this->getIterationCounter(),this->getTime());
520 0 : }
521 :
522 :
523 453 : void VesBias::setTargetDistAverages(const CoeffsVector& coeffderivs_aver_ps, const unsigned int coeffs_id) {
524 453 : TargetDistAverages(coeffs_id).setValues( coeffderivs_aver_ps );
525 453 : TargetDistAverages(coeffs_id).setIterationCounterAndTime(this->getIterationCounter(),this->getTime());
526 453 : }
527 :
528 :
529 0 : void VesBias::setTargetDistAveragesToZero(const unsigned int coeffs_id) {
530 0 : TargetDistAverages(coeffs_id).setAllValuesToZero();
531 0 : TargetDistAverages(coeffs_id).setIterationCounterAndTime(this->getIterationCounter(),this->getTime());
532 0 : }
533 :
534 :
535 175 : void VesBias::checkThatTemperatureIsGiven() {
536 175 : if(kbt_==0.0) {
537 0 : std::string err_msg = "VES bias " + getLabel() + " of type " + getName() + ": the temperature is needed so you need to give it using the TEMP keyword as the MD engine does not pass it to PLUMED.";
538 0 : plumed_merror(err_msg);
539 : }
540 175 : }
541 :
542 :
543 1005 : unsigned int VesBias::getIterationCounter() const {
544 : unsigned int iteration = 0;
545 1005 : if(optimizeCoeffs()) {
546 : iteration = getOptimizerPntr()->getIterationCounter();
547 : }
548 : else {
549 113 : iteration = getCoeffsPntrs()[0]->getIterationCounter();
550 : }
551 1005 : return iteration;
552 : }
553 :
554 :
555 85 : void VesBias::linkOptimizer(Optimizer* optimizer_pntr_in) {
556 : //
557 85 : if(optimizer_pntr_==NULL) {
558 85 : optimizer_pntr_ = optimizer_pntr_in;
559 : }
560 : else {
561 0 : std::string err_msg = "VES bias " + getLabel() + " of type " + getName() + " has already been linked with optimizer " + optimizer_pntr_->getLabel() + " of type " + optimizer_pntr_->getName() + ". You cannot link two optimizer to the same VES bias.";
562 0 : plumed_merror(err_msg);
563 : }
564 85 : checkThatTemperatureIsGiven();
565 85 : optimize_coeffs_ = true;
566 85 : filenames_have_iteration_number_ = true;
567 85 : }
568 :
569 :
570 82 : void VesBias::enableHessian(const bool diagonal_hessian) {
571 82 : compute_hessian_=true;
572 82 : diagonal_hessian_=diagonal_hessian;
573 82 : sampled_cross_averages.clear();
574 164 : for (unsigned int i=0; i<ncoeffssets_; i++) {
575 82 : std::string label = getCoeffsSetLabelString("hessian",i);
576 164 : hessian_pntrs_[i] = Tools::make_unique<CoeffsMatrix>(label,coeffs_pntrs_[i].get(),comm,diagonal_hessian_);
577 : //
578 : std::vector<double> cross_aver_sampled_tmp;
579 82 : cross_aver_sampled_tmp.assign(hessian_pntrs_[i]->getSize(),0.0);
580 82 : sampled_cross_averages.push_back(cross_aver_sampled_tmp);
581 : }
582 82 : }
583 :
584 :
585 0 : void VesBias::disableHessian() {
586 0 : compute_hessian_=false;
587 0 : diagonal_hessian_=true;
588 0 : sampled_cross_averages.clear();
589 0 : for (unsigned int i=0; i<ncoeffssets_; i++) {
590 0 : std::string label = getCoeffsSetLabelString("hessian",i);
591 0 : hessian_pntrs_[i] = Tools::make_unique<CoeffsMatrix>(label,coeffs_pntrs_[i].get(),comm,diagonal_hessian_);
592 : //
593 : std::vector<double> cross_aver_sampled_tmp;
594 0 : cross_aver_sampled_tmp.assign(hessian_pntrs_[i]->getSize(),0.0);
595 0 : sampled_cross_averages.push_back(cross_aver_sampled_tmp);
596 : }
597 0 : }
598 :
599 :
600 442 : std::string VesBias::getCoeffsSetLabelString(const std::string& type, const unsigned int coeffs_id) const {
601 442 : std::string label_prefix = getLabel() + ".";
602 442 : std::string label_postfix = "";
603 442 : if(use_multiple_coeffssets_) {
604 0 : Tools::convert(coeffs_id,label_postfix);
605 0 : label_postfix = "-" + label_postfix;
606 : }
607 1326 : return label_prefix+type+label_postfix;
608 : }
609 :
610 :
611 567 : std::unique_ptr<OFile> VesBias::getOFile(const std::string& filepath, const bool multi_sim_single_file, const bool enforce_backup) {
612 : auto ofile_pntr = Tools::make_unique<OFile>();
613 567 : std::string fp = filepath;
614 567 : ofile_pntr->link(*static_cast<Action*>(this));
615 567 : if(enforce_backup) {ofile_pntr->enforceBackup();}
616 567 : if(multi_sim_single_file) {
617 56 : unsigned int r=0;
618 56 : if(comm.Get_rank()==0) {r=multi_sim_comm.Get_rank();}
619 56 : comm.Bcast(r,0);
620 56 : if(r>0) {fp="/dev/null";}
621 112 : ofile_pntr->enforceSuffix("");
622 : }
623 567 : ofile_pntr->open(fp);
624 567 : return ofile_pntr;
625 0 : }
626 :
627 :
628 0 : void VesBias::setGridBins(const std::vector<unsigned int>& grid_bins_in) {
629 0 : plumed_massert(grid_bins_in.size()==getNumberOfArguments(),"the number of grid bins given doesn't match the number of arguments");
630 0 : grid_bins_=grid_bins_in;
631 0 : }
632 :
633 :
634 0 : void VesBias::setGridBins(const unsigned int nbins) {
635 0 : std::vector<unsigned int> grid_bins_in(getNumberOfArguments(),nbins);
636 0 : grid_bins_=grid_bins_in;
637 0 : }
638 :
639 :
640 0 : void VesBias::setGridMin(const std::vector<double>& grid_min_in) {
641 0 : plumed_massert(grid_min_in.size()==getNumberOfArguments(),"the number of lower bounds given for the grid doesn't match the number of arguments");
642 0 : grid_min_=grid_min_in;
643 0 : }
644 :
645 :
646 0 : void VesBias::setGridMax(const std::vector<double>& grid_max_in) {
647 0 : plumed_massert(grid_max_in.size()==getNumberOfArguments(),"the number of upper bounds given for the grid doesn't match the number of arguments");
648 0 : grid_max_=grid_max_in;
649 0 : }
650 :
651 :
652 383 : std::string VesBias::getCurrentOutputFilename(const std::string& base_filename, const std::string& suffix) const {
653 383 : std::string filename = base_filename;
654 383 : if(suffix.size()>0) {
655 82 : filename = FileBase::appendSuffix(filename,"."+suffix);
656 : }
657 383 : if(filenamesIncludeIterationNumber()) {
658 756 : filename = FileBase::appendSuffix(filename,"."+getIterationFilenameSuffix());
659 : }
660 383 : return filename;
661 : }
662 :
663 :
664 192 : std::string VesBias::getCurrentTargetDistOutputFilename(const std::string& suffix) const {
665 192 : std::string filename = targetdist_filename_;
666 192 : if(suffix.size()>0) {
667 204 : filename = FileBase::appendSuffix(filename,"."+suffix);
668 : }
669 192 : if(filenamesIncludeIterationNumber() && dynamicTargetDistribution()) {
670 348 : filename = FileBase::appendSuffix(filename,"."+getIterationFilenameSuffix());
671 : }
672 192 : return filename;
673 : }
674 :
675 :
676 552 : std::string VesBias::getIterationFilenameSuffix() const {
677 : std::string iter_str;
678 552 : Tools::convert(getIterationCounter(),iter_str);
679 552 : iter_str = "iter-" + iter_str;
680 552 : return iter_str;
681 : }
682 :
683 :
684 0 : std::string VesBias::getCoeffsSetFilenameSuffix(const unsigned int coeffs_id) const {
685 0 : std::string suffix = "";
686 0 : if(use_multiple_coeffssets_) {
687 0 : Tools::convert(coeffs_id,suffix);
688 0 : suffix = coeffs_id_prefix_ + suffix;
689 : }
690 0 : return suffix;
691 : }
692 :
693 :
694 3 : void VesBias::setupBiasCutoff(const double bias_cutoff_value, const double fermi_lambda) {
695 : //
696 : double fermi_exp_max = 100.0;
697 : //
698 : std::string str_bias_cutoff_value; VesTools::convertDbl2Str(bias_cutoff_value,str_bias_cutoff_value);
699 : std::string str_fermi_lambda; VesTools::convertDbl2Str(fermi_lambda,str_fermi_lambda);
700 : std::string str_fermi_exp_max; VesTools::convertDbl2Str(fermi_exp_max,str_fermi_exp_max);
701 6 : std::string swfunc_keywords = "FERMI R_0=" + str_bias_cutoff_value + " FERMI_LAMBDA=" + str_fermi_lambda + " FERMI_EXP_MAX=" + str_fermi_exp_max;
702 : //
703 3 : std::string swfunc_errors="";
704 3 : bias_cutoff_swfunc_pntr_ = Tools::make_unique<FermiSwitchingFunction>();
705 3 : bias_cutoff_swfunc_pntr_->set(swfunc_keywords,swfunc_errors);
706 3 : if(swfunc_errors.size()>0) {
707 0 : plumed_merror("problem with setting up Fermi switching function: " + swfunc_errors);
708 : }
709 : //
710 3 : bias_cutoff_value_=bias_cutoff_value;
711 3 : bias_cutoff_active_=true;
712 : enableDynamicTargetDistribution();
713 3 : }
714 :
715 :
716 3263 : double VesBias::getBiasCutoffSwitchingFunction(const double bias, double& deriv_factor) const {
717 3263 : plumed_massert(bias_cutoff_active_,"The bias cutoff is not active so you cannot call this function");
718 3263 : double arg = -(bias-bias_current_max_value);
719 3263 : double deriv=0.0;
720 3263 : double value = bias_cutoff_swfunc_pntr_->calculate(arg,deriv);
721 : // as FermiSwitchingFunction class has different behavior from normal SwitchingFunction class
722 : // I was having problems with NaN as it was dividing with zero
723 : // deriv *= arg;
724 3263 : deriv_factor = value-bias*deriv;
725 3263 : return value;
726 : }
727 :
728 :
729 567 : bool VesBias::useMultipleWalkers() const {
730 : bool use_mwalkers_mpi=false;
731 567 : if(optimizeCoeffs() && getOptimizerPntr()->useMultipleWalkers()) {
732 : use_mwalkers_mpi=true;
733 : }
734 567 : return use_mwalkers_mpi;
735 : }
736 :
737 :
738 0 : void VesBias::updateReweightFactor() {
739 0 : if(calc_reweightfactor_) {
740 0 : double value = calculateReweightFactor();
741 0 : getPntrToComponent("rct")->set(value);
742 : }
743 0 : }
744 :
745 :
746 0 : double VesBias::calculateReweightFactor() const {
747 0 : plumed_merror(getName()+" with label "+getLabel()+": calculation of the reweight factor c(t) has not been implemented for this type of VES bias");
748 : return 0.0;
749 : }
750 :
751 :
752 : }
753 : }
|
