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