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