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