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