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