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 : #ifndef __PLUMED_ves_VesBias_h
23 : #define __PLUMED_ves_VesBias_h
24 :
25 : #include "CoeffsVector.h"
26 : #include "CoeffsMatrix.h"
27 :
28 : #include "core/ActionPilot.h"
29 : #include "core/ActionWithValue.h"
30 : #include "core/ActionWithArguments.h"
31 : #include "bias/Bias.h"
32 :
33 : #include <vector>
34 : #include <string>
35 : #include <cmath>
36 : #include <memory>
37 :
38 :
39 : #define PLUMED_VES_VESBIAS_INIT(ao) Action(ao),VesBias(ao)
40 :
41 : namespace PLMD {
42 :
43 : class Value;
44 :
45 : namespace ves {
46 :
47 : class CoeffsVector;
48 : class CoeffsMatrix;
49 : class BasisFunctions;
50 : class Optimizer;
51 : class TargetDistribution;
52 : class FermiSwitchingFunction;
53 :
54 : /**
55 : \ingroup INHERIT
56 : Abstract base class for implementing biases the extents the normal Bias.h class
57 : to include functions related to the variational approach.
58 : */
59 :
60 : class VesBias:
61 : public bias::Bias {
62 : private:
63 : unsigned int ncoeffssets_;
64 : std::vector<std::unique_ptr<CoeffsVector>> coeffs_pntrs_;
65 : std::vector<std::unique_ptr<CoeffsVector>> targetdist_averages_pntrs_;
66 : std::vector<std::unique_ptr<CoeffsVector>> gradient_pntrs_;
67 : std::vector<std::unique_ptr<CoeffsMatrix>> hessian_pntrs_;
68 : std::vector<std::vector<double> > sampled_averages;
69 : std::vector<std::vector<double> > sampled_cross_averages;
70 : bool use_multiple_coeffssets_;
71 : //
72 : std::vector<std::string> coeffs_fnames;
73 : //
74 : size_t ncoeffs_total_;
75 : //
76 : Optimizer* optimizer_pntr_;
77 : bool optimize_coeffs_;
78 : //
79 : bool compute_hessian_;
80 : bool diagonal_hessian_;
81 : //
82 : std::vector<unsigned int> aver_counters;
83 : //
84 : double kbt_;
85 : //
86 : std::vector<TargetDistribution*> targetdist_pntrs_;
87 : bool dynamic_targetdist_;
88 : //
89 : std::vector<unsigned int> grid_bins_;
90 : std::vector<double> grid_min_;
91 : std::vector<double> grid_max_;
92 : //
93 : std::string bias_filename_;
94 : std::string fes_filename_;
95 : std::string targetdist_filename_;
96 : std::string targetdist_averages_filename_;
97 : std::string coeffs_id_prefix_;
98 : //
99 : std::string bias_file_fmt_;
100 : std::string fes_file_fmt_;
101 : std::string targetdist_file_fmt_;
102 : std::string targetdist_restart_file_fmt_;
103 : //
104 : bool filenames_have_iteration_number_;
105 : //
106 : bool bias_fileoutput_active_;
107 : bool fes_fileoutput_active_;
108 : bool fesproj_fileoutput_active_;
109 : bool dynamic_targetdist_fileoutput_active_;
110 : bool static_targetdist_fileoutput_active_;
111 : //
112 :
113 : bool bias_cutoff_active_;
114 : double bias_cutoff_value_;
115 : double bias_current_max_value;
116 : std::unique_ptr<FermiSwitchingFunction> bias_cutoff_swfunc_pntr_;
117 : //
118 : std::vector< std::vector<std::string> > projection_args_;
119 : //
120 : bool calc_reweightfactor_;
121 : //
122 : double optimization_threshold_;
123 : private:
124 : void initializeCoeffs(std::unique_ptr<CoeffsVector>);
125 : std::vector<double> computeCovarianceFromAverages(const unsigned int) const;
126 : void multiSimSumAverages(const unsigned int, const double walker_weight=1.0);
127 : protected:
128 : //
129 : void checkThatTemperatureIsGiven();
130 : //
131 : void addCoeffsSet(const std::vector<std::string>&,const std::vector<unsigned int>&);
132 : void addCoeffsSet(std::vector<Value*>&,std::vector<BasisFunctions*>&);
133 : void addCoeffsSet(std::unique_ptr<CoeffsVector>);
134 : //
135 : std::string getCoeffsSetLabelString(const std::string&, const unsigned int coeffs_id = 0) const;
136 : void clearCoeffsPntrsVector() {
137 : coeffs_pntrs_.clear();
138 : }
139 : void addToSampledAverages(const std::vector<double>&, const unsigned int c_id = 0);
140 : void setTargetDistAverages(const std::vector<double>&, const unsigned int coeffs_id = 0);
141 : void setTargetDistAverages(const CoeffsVector&, const unsigned int coeffs_id= 0);
142 : void setTargetDistAveragesToZero(const unsigned int coeffs_id= 0);
143 : //
144 : bool readCoeffsFromFiles();
145 : //
146 : template<class T>
147 : bool parseMultipleValues(const std::string&, std::vector<T>&, unsigned int);
148 : template<class T>
149 : bool parseMultipleValues(const std::string&, std::vector<T>&, unsigned int, const T&);
150 : //
151 : public:
152 : static void registerKeywords(Keywords&);
153 : explicit VesBias(const ActionOptions&ao);
154 : ~VesBias();
155 : //
156 : static void useInitialCoeffsKeywords(Keywords&);
157 : static void useTargetDistributionKeywords(Keywords&);
158 : static void useMultipleTargetDistributionKeywords(Keywords&);
159 : static void useGridBinKeywords(Keywords&);
160 : static void useGridLimitsKeywords(Keywords&);
161 : static void useBiasCutoffKeywords(Keywords&);
162 : static void useProjectionArgKeywords(Keywords&);
163 : static void useReweightFactorKeywords(Keywords&);
164 : //
165 : std::vector<CoeffsVector*> getCoeffsPntrs() const {
166 267 : return Tools::unique2raw(coeffs_pntrs_);
167 : }
168 : std::vector<CoeffsVector*> getTargetDistAveragesPntrs() const {
169 85 : return Tools::unique2raw(targetdist_averages_pntrs_);
170 : }
171 : std::vector<CoeffsVector*> getGradientPntrs()const {
172 85 : return Tools::unique2raw(gradient_pntrs_);
173 : }
174 : std::vector<CoeffsMatrix*> getHessianPntrs() const {
175 82 : return Tools::unique2raw(hessian_pntrs_);
176 : }
177 : std::vector<TargetDistribution*> getTargetDistributionPntrs() const {
178 93 : return targetdist_pntrs_;
179 : }
180 : //
181 : CoeffsVector* getCoeffsPntr(const unsigned int coeffs_id = 0) const {
182 : return coeffs_pntrs_[coeffs_id].get();
183 : }
184 : CoeffsVector* getTargetDistAveragesPntr(const unsigned int coeffs_id = 0) const {
185 : return targetdist_averages_pntrs_[coeffs_id].get();
186 : }
187 : CoeffsVector* getGradientPntr(const unsigned int coeffs_id = 0)const {
188 : return gradient_pntrs_[coeffs_id].get();
189 : }
190 : CoeffsMatrix* getHessianPntr(const unsigned int coeffs_id = 0) const {
191 : return hessian_pntrs_[coeffs_id].get();
192 : }
193 : //
194 : unsigned int getNumberOfTargetDistributionPntrs() const {
195 90 : return targetdist_pntrs_.size();
196 : }
197 : //
198 : size_t numberOfCoeffs(const unsigned int coeffs_id = 0) const {
199 22810 : return coeffs_pntrs_[coeffs_id]->numberOfCoeffs();
200 : }
201 : size_t totalNumberOfCoeffs() const {
202 : return ncoeffs_total_;
203 : }
204 : unsigned int numberOfCoeffsSets() const {
205 3 : return ncoeffssets_;
206 : }
207 : double getKbT() const {
208 85 : return kbt_;
209 : }
210 : double getBeta() const;
211 : //
212 : CoeffsVector& Coeffs(const unsigned int coeffs_id = 0) const {
213 : return *coeffs_pntrs_[coeffs_id];
214 : }
215 : CoeffsVector& TargetDistAverages(const unsigned int coeffs_id = 0) const {
216 453 : return *targetdist_averages_pntrs_[coeffs_id];
217 : }
218 : CoeffsVector& Gradient(const unsigned int coeffs_id = 0) const {
219 : return *gradient_pntrs_[coeffs_id];
220 : }
221 : CoeffsMatrix& Hessian(const unsigned int coeffs_id = 0) const {
222 : return *hessian_pntrs_[coeffs_id];
223 : }
224 : //
225 : size_t getCoeffsIndex(const std::vector<unsigned int>& indices, const unsigned int coeffs_id = 0) const;
226 : std::vector<unsigned int> getCoeffsIndices(const size_t index, const unsigned int coeffs_id = 0) const;
227 : size_t getHessianIndex(const size_t index1, const size_t index2, const unsigned int coeffs_id = 0) const;
228 : //
229 : bool computeHessian() const {
230 : return compute_hessian_;
231 : }
232 : bool diagonalHessian() const {
233 : return diagonal_hessian_;
234 : }
235 : //
236 : bool optimizeCoeffs() const {
237 1572 : return optimize_coeffs_;
238 : }
239 : Optimizer* getOptimizerPntr() const {
240 1426 : return optimizer_pntr_;
241 : }
242 : bool useMultipleWalkers() const;
243 : //
244 : unsigned int getIterationCounter() const;
245 : //
246 : void updateGradientAndHessian(const bool);
247 : void clearGradientAndHessian() {};
248 : //
249 0 : virtual void updateTargetDistributions() {};
250 0 : virtual void restartTargetDistributions() {};
251 : //
252 : void linkOptimizer(Optimizer*);
253 : void enableHessian(const bool diagonal_hessian=true);
254 : void disableHessian();
255 : //
256 : void enableMultipleCoeffsSets() {
257 : use_multiple_coeffssets_=true;
258 : }
259 : //
260 : void enableDynamicTargetDistribution() {
261 42 : dynamic_targetdist_=true;
262 39 : }
263 : void disableDynamicTargetDistribution() {
264 : dynamic_targetdist_=false;
265 : }
266 : bool dynamicTargetDistribution() const {
267 273 : return dynamic_targetdist_;
268 : }
269 : //
270 : std::vector<unsigned int> getGridBins() const {
271 90 : return grid_bins_;
272 : }
273 : void setGridBins(const std::vector<unsigned int>&);
274 : void setGridBins(const unsigned int);
275 : std::vector<double> getGridMax() const {
276 : return grid_max_;
277 : }
278 : void setGridMax(const std::vector<double>&);
279 : std::vector<double> getGridMin() const {
280 : return grid_min_;
281 : }
282 : void setGridMin(const std::vector<double>&);
283 : //
284 : bool filenamesIncludeIterationNumber() const {
285 575 : return filenames_have_iteration_number_;
286 : }
287 : void enableIterationNumberInFilenames() {
288 3 : filenames_have_iteration_number_=true;
289 : }
290 : //
291 : std::string getIterationFilenameSuffix() const;
292 : std::string getCoeffsSetFilenameSuffix(const unsigned int coeffs_id) const;
293 : std::string getCurrentOutputFilename(const std::string&, const std::string& suffix="") const;
294 : std::string getBiasOutputFilename() const {
295 : return bias_filename_;
296 : }
297 : std::string getCurrentBiasOutputFilename(const std::string& suffix="") const;
298 : std::string getFesOutputFilename() const {
299 : return fes_filename_;
300 : }
301 : std::string getCurrentFesOutputFilename(const std::string& suffix="") const;
302 : std::string getTargetDistOutputFilename() const {
303 : return targetdist_filename_;
304 : }
305 : std::string getCurrentTargetDistOutputFilename(const std::string& suffix="") const;
306 : //
307 : void enableBiasFileOutput() {
308 83 : bias_fileoutput_active_=true;
309 : }
310 : void disableBiasFileOutput() {
311 : bias_fileoutput_active_=false;
312 : }
313 : bool isBiasFileOutputActive() const {
314 : return bias_fileoutput_active_;
315 : }
316 : std::string getBiasFileFmt() const {
317 : return bias_file_fmt_;
318 : }
319 : //
320 : void enableFesFileOutput() {
321 83 : fes_fileoutput_active_=true;
322 : }
323 : void disableFesFileOutput() {
324 : fes_fileoutput_active_=false;
325 : }
326 : bool isFesFileOutputActive() const {
327 : return fes_fileoutput_active_;
328 : }
329 : std::string getFesFileFmt() const {
330 : return fes_file_fmt_;
331 : }
332 : //
333 : void enableFesProjFileOutput() {
334 17 : fesproj_fileoutput_active_=true;
335 : }
336 : void disableFesFileProjOutput() {
337 : fesproj_fileoutput_active_=false;
338 : }
339 : bool isFesProjFileOutputActive() const {
340 : return fesproj_fileoutput_active_;
341 : }
342 : //
343 : void enableDynamicTargetDistFileOutput() {
344 41 : dynamic_targetdist_fileoutput_active_=true;
345 : }
346 : void disableDynamicTargetDistFileOutput() {
347 : dynamic_targetdist_fileoutput_active_=false;
348 : }
349 : bool isDynamicTargetDistFileOutputActive() const {
350 : return dynamic_targetdist_fileoutput_active_;
351 : }
352 : std::string getTargetDistFileFmt() const {
353 : return targetdist_file_fmt_;
354 : }
355 : std::string getTargetDistRestartFileFmt() const {
356 : return targetdist_restart_file_fmt_;
357 : }
358 : //
359 : void enableStaticTargetDistFileOutput() {
360 : static_targetdist_fileoutput_active_=true;
361 : }
362 : void disableStaticTargetDistFileOutput() {
363 46 : static_targetdist_fileoutput_active_=false;
364 46 : }
365 : bool isStaticTargetDistFileOutputActive() const {
366 50 : return static_targetdist_fileoutput_active_;
367 : }
368 : //
369 : std::vector< std::vector<std::string> > getProjectionArguments() const {
370 : return projection_args_;
371 : }
372 : std::vector<std::string> getProjectionArgument(unsigned int i) const {
373 56 : return projection_args_[i];
374 : }
375 : unsigned int getNumberOfProjectionArguments() const {
376 122 : return projection_args_.size();
377 : }
378 : //
379 : void setupBiasCutoff(const double, const double);
380 : bool biasCutoffActive() const {
381 1498465 : return bias_cutoff_active_;
382 : }
383 : double getBiasCutoffValue() const {
384 45 : return bias_cutoff_value_;
385 : }
386 : void setCurrentBiasMaxValue(const double max_value) {
387 498 : bias_current_max_value=max_value;
388 498 : }
389 : double getCurrentBiasMaxValue() const {
390 : return bias_current_max_value;
391 : }
392 : double getBiasCutoffSwitchingFunction(const double, double&) const;
393 : double getBiasCutoffSwitchingFunction(const double) const;
394 : void applyBiasCutoff(double&, std::vector<double>&) const;
395 : void applyBiasCutoff(double&, std::vector<double>&, std::vector<double>&) const;
396 : //
397 : std::unique_ptr<OFile> getOFile(const std::string& filename, const bool multi_sim_single_file=false, const bool enforce_backup=true);
398 : //
399 0 : virtual void setupBiasFileOutput() {};
400 0 : virtual void writeBiasToFile() {};
401 0 : virtual void resetBiasFileOutput() {};
402 : //
403 0 : virtual void setupFesFileOutput() {};
404 0 : virtual void writeFesToFile() {};
405 0 : virtual void resetFesFileOutput() {};
406 : //
407 0 : virtual void setupFesProjFileOutput() {};
408 0 : virtual void writeFesProjToFile() {};
409 0 : virtual void resetFesProjFileOutput() {};
410 : //
411 44 : virtual void setupTargetDistFileOutput() {};
412 0 : virtual void writeTargetDistToFile() {};
413 0 : virtual void resetTargetDistFileOutput() {};
414 : //
415 9 : virtual void setupTargetDistProjFileOutput() {};
416 0 : virtual void writeTargetDistProjToFile() {};
417 0 : virtual void resetTargetDistProjFileOutput() {};
418 : //
419 : void updateReweightFactor();
420 : virtual double calculateReweightFactor() const;
421 : bool isReweightFactorCalculated() const {
422 0 : return calc_reweightfactor_;
423 : }
424 : };
425 :
426 :
427 : inline
428 : size_t VesBias::getCoeffsIndex(const std::vector<unsigned int>& indices, const unsigned int coeffs_id) const {
429 : return coeffs_pntrs_[coeffs_id]->getIndex(indices);
430 : }
431 :
432 : inline
433 : std::vector<unsigned int> VesBias::getCoeffsIndices(const size_t index, const unsigned int coeffs_id) const {
434 : return coeffs_pntrs_[coeffs_id]->getIndices(index);
435 : }
436 :
437 : inline
438 : size_t VesBias::getHessianIndex(const size_t index1, const size_t index2, const unsigned int coeffs_id) const {
439 3607712 : return hessian_pntrs_[coeffs_id]->getMatrixIndex(index1,index2);
440 : }
441 :
442 :
443 : inline
444 23277 : double VesBias::getBeta() const {
445 23277 : plumed_massert(kbt_!=0.0,"you are requesting beta=1/(kB*T) when kB*T has not been defined. You need to give the temperature using the TEMP keyword as the MD engine does not pass it to PLUMED.");
446 23277 : return 1.0/kbt_;
447 : }
448 :
449 :
450 : inline
451 : std::string VesBias::getCurrentBiasOutputFilename(const std::string& suffix) const {
452 178 : return getCurrentOutputFilename(bias_filename_,suffix);
453 : }
454 :
455 :
456 : inline
457 : std::string VesBias::getCurrentFesOutputFilename(const std::string& suffix) const {
458 205 : return getCurrentOutputFilename(fes_filename_,suffix);
459 : }
460 :
461 :
462 : inline
463 : double VesBias::getBiasCutoffSwitchingFunction(const double bias) const {
464 : double dummy=0.0;
465 : return getBiasCutoffSwitchingFunction(bias,dummy);
466 : }
467 :
468 :
469 : inline
470 600 : void VesBias::applyBiasCutoff(double& bias, std::vector<double>& forces) const {
471 600 : std::vector<double> dummy(0);
472 600 : applyBiasCutoff(bias,forces,dummy);
473 600 : }
474 :
475 :
476 : inline
477 663 : void VesBias::applyBiasCutoff(double& bias, std::vector<double>& forces, std::vector<double>& coeffsderivs_values) const {
478 663 : double deriv_factor_sf=0.0;
479 663 : double value_sf = getBiasCutoffSwitchingFunction(bias,deriv_factor_sf);
480 663 : bias *= value_sf;
481 1326 : for(unsigned int i=0; i<forces.size(); i++) {
482 663 : forces[i] *= deriv_factor_sf;
483 : }
484 : //
485 1398 : for(unsigned int i=0; i<coeffsderivs_values.size(); i++) {
486 735 : coeffsderivs_values[i] *= deriv_factor_sf;
487 : }
488 663 : }
489 :
490 :
491 : inline
492 22810 : std::vector<double> VesBias::computeCovarianceFromAverages(const unsigned int c_id) const {
493 : size_t ncoeffs = numberOfCoeffs(c_id);
494 22810 : std::vector<double> covariance(sampled_cross_averages[c_id].size(),0.0);
495 : // diagonal part
496 1823850 : for(size_t i=0; i<ncoeffs; i++) {
497 : size_t midx = getHessianIndex(i,i,c_id);
498 1801040 : covariance[midx] = sampled_cross_averages[c_id][midx] - sampled_averages[c_id][i]*sampled_averages[c_id][i];
499 : }
500 22810 : if(!diagonal_hessian_) {
501 0 : for(size_t i=0; i<ncoeffs; i++) {
502 0 : for(size_t j=(i+1); j<ncoeffs; j++) {
503 : size_t midx = getHessianIndex(i,j,c_id);
504 0 : covariance[midx] = sampled_cross_averages[c_id][midx] - sampled_averages[c_id][i]*sampled_averages[c_id][j];
505 : }
506 : }
507 : }
508 22810 : return covariance;
509 : }
510 :
511 :
512 : template<class T>
513 180 : bool VesBias::parseMultipleValues(const std::string& keyword, std::vector<T>& values, unsigned int nvalues) {
514 180 : plumed_assert(nvalues>0);
515 180 : plumed_assert(values.size()==0);
516 : bool identical_values=false;
517 : //
518 180 : parseVector(keyword,values);
519 180 : if(values.size()==1 && nvalues>1) {
520 0 : values.resize(nvalues,values[0]);
521 : identical_values=true;
522 : }
523 180 : if(values.size()>0 && values.size()!=nvalues) {
524 : std::string s1;
525 0 : Tools::convert(nvalues,s1);
526 0 : plumed_merror("Error in " + keyword + " keyword: either give 1 common parameter value or " + s1 + " separate parameter values");
527 : }
528 180 : return identical_values;
529 : }
530 :
531 : template<class T>
532 90 : bool VesBias::parseMultipleValues(const std::string& keyword, std::vector<T>& values, unsigned int nvalues, const T& default_value) {
533 90 : bool identical_values = parseMultipleValues(keyword,values,nvalues);
534 90 : if(values.size()==0) {
535 0 : values.resize(nvalues,default_value);
536 : identical_values=true;
537 : }
538 90 : return identical_values;
539 : }
540 :
541 :
542 : }
543 : }
544 :
545 : #endif
|
