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