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 : {
63 : private:
64 : unsigned int ncoeffssets_;
65 : std::vector<std::unique_ptr<CoeffsVector>> coeffs_pntrs_;
66 : std::vector<std::unique_ptr<CoeffsVector>> targetdist_averages_pntrs_;
67 : std::vector<std::unique_ptr<CoeffsVector>> gradient_pntrs_;
68 : std::vector<std::unique_ptr<CoeffsMatrix>> hessian_pntrs_;
69 : std::vector<std::vector<double> > sampled_averages;
70 : std::vector<std::vector<double> > sampled_cross_averages;
71 : bool use_multiple_coeffssets_;
72 : //
73 : std::vector<std::string> coeffs_fnames;
74 : //
75 : size_t ncoeffs_total_;
76 : //
77 : Optimizer* optimizer_pntr_;
78 : bool optimize_coeffs_;
79 : //
80 : bool compute_hessian_;
81 : bool diagonal_hessian_;
82 : //
83 : std::vector<unsigned int> aver_counters;
84 : //
85 : double kbt_;
86 : //
87 : std::vector<TargetDistribution*> targetdist_pntrs_;
88 : bool dynamic_targetdist_;
89 : //
90 : std::vector<unsigned int> grid_bins_;
91 : std::vector<double> grid_min_;
92 : std::vector<double> grid_max_;
93 : //
94 : std::string bias_filename_;
95 : std::string fes_filename_;
96 : std::string targetdist_filename_;
97 : std::string targetdist_averages_filename_;
98 : std::string coeffs_id_prefix_;
99 : //
100 : std::string bias_file_fmt_;
101 : std::string fes_file_fmt_;
102 : std::string targetdist_file_fmt_;
103 : std::string targetdist_restart_file_fmt_;
104 : //
105 : bool filenames_have_iteration_number_;
106 : //
107 : bool bias_fileoutput_active_;
108 : bool fes_fileoutput_active_;
109 : bool fesproj_fileoutput_active_;
110 : bool dynamic_targetdist_fileoutput_active_;
111 : bool static_targetdist_fileoutput_active_;
112 : //
113 :
114 : bool bias_cutoff_active_;
115 : double bias_cutoff_value_;
116 : double bias_current_max_value;
117 : std::unique_ptr<FermiSwitchingFunction> bias_cutoff_swfunc_pntr_;
118 : //
119 : std::vector< std::vector<std::string> > projection_args_;
120 : //
121 : bool calc_reweightfactor_;
122 : //
123 : double optimization_threshold_;
124 : private:
125 : void initializeCoeffs(std::unique_ptr<CoeffsVector>);
126 : std::vector<double> computeCovarianceFromAverages(const unsigned int) const;
127 : void multiSimSumAverages(const unsigned int, const double walker_weight=1.0);
128 : protected:
129 : //
130 : void checkThatTemperatureIsGiven();
131 : //
132 : void addCoeffsSet(const std::vector<std::string>&,const std::vector<unsigned int>&);
133 : void addCoeffsSet(std::vector<Value*>&,std::vector<BasisFunctions*>&);
134 : void addCoeffsSet(std::unique_ptr<CoeffsVector>);
135 : //
136 : std::string getCoeffsSetLabelString(const std::string&, const unsigned int coeffs_id = 0) const;
137 : void clearCoeffsPntrsVector() {coeffs_pntrs_.clear();}
138 : void addToSampledAverages(const std::vector<double>&, const unsigned int c_id = 0);
139 : void setTargetDistAverages(const std::vector<double>&, const unsigned int coeffs_id = 0);
140 : void setTargetDistAverages(const CoeffsVector&, const unsigned int coeffs_id= 0);
141 : void setTargetDistAveragesToZero(const unsigned int coeffs_id= 0);
142 : //
143 : bool readCoeffsFromFiles();
144 : //
145 : template<class T>
146 : bool parseMultipleValues(const std::string&, std::vector<T>&, unsigned int);
147 : template<class T>
148 : bool parseMultipleValues(const std::string&, std::vector<T>&, unsigned int, const T&);
149 : //
150 : public:
151 : static void registerKeywords(Keywords&);
152 : explicit VesBias(const ActionOptions&ao);
153 : ~VesBias();
154 : //
155 : static void useInitialCoeffsKeywords(Keywords&);
156 : static void useTargetDistributionKeywords(Keywords&);
157 : static void useMultipleTargetDistributionKeywords(Keywords&);
158 : static void useGridBinKeywords(Keywords&);
159 : static void useGridLimitsKeywords(Keywords&);
160 : static void useBiasCutoffKeywords(Keywords&);
161 : static void useProjectionArgKeywords(Keywords&);
162 : static void useReweightFactorKeywords(Keywords&);
163 : //
164 267 : std::vector<CoeffsVector*> getCoeffsPntrs() const {return Tools::unique2raw(coeffs_pntrs_);}
165 85 : std::vector<CoeffsVector*> getTargetDistAveragesPntrs() const {return Tools::unique2raw(targetdist_averages_pntrs_);}
166 85 : std::vector<CoeffsVector*> getGradientPntrs()const {return Tools::unique2raw(gradient_pntrs_);}
167 82 : std::vector<CoeffsMatrix*> getHessianPntrs() const {return Tools::unique2raw(hessian_pntrs_);}
168 93 : std::vector<TargetDistribution*> getTargetDistributionPntrs() const {return targetdist_pntrs_;}
169 : //
170 : CoeffsVector* getCoeffsPntr(const unsigned int coeffs_id = 0) const {return coeffs_pntrs_[coeffs_id].get();}
171 : CoeffsVector* getTargetDistAveragesPntr(const unsigned int coeffs_id = 0) const {return targetdist_averages_pntrs_[coeffs_id].get();}
172 : CoeffsVector* getGradientPntr(const unsigned int coeffs_id = 0)const {return gradient_pntrs_[coeffs_id].get();}
173 : CoeffsMatrix* getHessianPntr(const unsigned int coeffs_id = 0) const {return hessian_pntrs_[coeffs_id].get();}
174 : //
175 90 : unsigned int getNumberOfTargetDistributionPntrs() const {return targetdist_pntrs_.size();}
176 : //
177 22810 : size_t numberOfCoeffs(const unsigned int coeffs_id = 0) const {return coeffs_pntrs_[coeffs_id]->numberOfCoeffs();}
178 : size_t totalNumberOfCoeffs() const {return ncoeffs_total_;}
179 3 : unsigned int numberOfCoeffsSets() const {return ncoeffssets_;}
180 85 : double getKbT() const {return kbt_;}
181 : double getBeta() const;
182 : //
183 : CoeffsVector& Coeffs(const unsigned int coeffs_id = 0) const {return *coeffs_pntrs_[coeffs_id];}
184 453 : CoeffsVector& TargetDistAverages(const unsigned int coeffs_id = 0) const {return *targetdist_averages_pntrs_[coeffs_id];}
185 : CoeffsVector& Gradient(const unsigned int coeffs_id = 0) const {return *gradient_pntrs_[coeffs_id];}
186 : CoeffsMatrix& Hessian(const unsigned int coeffs_id = 0) const {return *hessian_pntrs_[coeffs_id];}
187 : //
188 : size_t getCoeffsIndex(const std::vector<unsigned int>& indices, const unsigned int coeffs_id = 0) const;
189 : std::vector<unsigned int> getCoeffsIndices(const size_t index, const unsigned int coeffs_id = 0) const;
190 : size_t getHessianIndex(const size_t index1, const size_t index2, const unsigned int coeffs_id = 0) const;
191 : //
192 : bool computeHessian() const {return compute_hessian_;}
193 : bool diagonalHessian() const {return diagonal_hessian_;}
194 : //
195 1572 : bool optimizeCoeffs() const {return optimize_coeffs_;}
196 1426 : Optimizer* getOptimizerPntr() const {return optimizer_pntr_;}
197 : bool useMultipleWalkers() const;
198 : //
199 : unsigned int getIterationCounter() const;
200 : //
201 : void updateGradientAndHessian(const bool);
202 : void clearGradientAndHessian() {};
203 : //
204 0 : virtual void updateTargetDistributions() {};
205 0 : virtual void restartTargetDistributions() {};
206 : //
207 : void linkOptimizer(Optimizer*);
208 : void enableHessian(const bool diagonal_hessian=true);
209 : void disableHessian();
210 : //
211 : void enableMultipleCoeffsSets() {use_multiple_coeffssets_=true;}
212 : //
213 42 : void enableDynamicTargetDistribution() {dynamic_targetdist_=true;}
214 : void disableDynamicTargetDistribution() {dynamic_targetdist_=false;}
215 273 : bool dynamicTargetDistribution() const {return dynamic_targetdist_;}
216 : //
217 90 : std::vector<unsigned int> getGridBins() const {return grid_bins_;}
218 : void setGridBins(const std::vector<unsigned int>&);
219 : void setGridBins(const unsigned int);
220 : std::vector<double> getGridMax() const {return grid_max_;}
221 : void setGridMax(const std::vector<double>&);
222 : std::vector<double> getGridMin() const {return grid_min_;}
223 : void setGridMin(const std::vector<double>&);
224 : //
225 575 : bool filenamesIncludeIterationNumber() const {return filenames_have_iteration_number_;}
226 3 : void enableIterationNumberInFilenames() {filenames_have_iteration_number_=true;}
227 : //
228 : std::string getIterationFilenameSuffix() const;
229 : std::string getCoeffsSetFilenameSuffix(const unsigned int coeffs_id) const;
230 : std::string getCurrentOutputFilename(const std::string&, const std::string& suffix="") const;
231 : std::string getBiasOutputFilename() const {return bias_filename_;}
232 : std::string getCurrentBiasOutputFilename(const std::string& suffix="") const;
233 : std::string getFesOutputFilename() const {return fes_filename_;}
234 : std::string getCurrentFesOutputFilename(const std::string& suffix="") const;
235 : std::string getTargetDistOutputFilename() const {return targetdist_filename_;}
236 : std::string getCurrentTargetDistOutputFilename(const std::string& suffix="") const;
237 : //
238 83 : void enableBiasFileOutput() {bias_fileoutput_active_=true;}
239 : void disableBiasFileOutput() {bias_fileoutput_active_=false;}
240 : bool isBiasFileOutputActive() const {return bias_fileoutput_active_;}
241 : std::string getBiasFileFmt() const {return bias_file_fmt_;}
242 : //
243 83 : void enableFesFileOutput() {fes_fileoutput_active_=true;}
244 : void disableFesFileOutput() {fes_fileoutput_active_=false;}
245 : bool isFesFileOutputActive() const {return fes_fileoutput_active_;}
246 : std::string getFesFileFmt() const {return fes_file_fmt_;}
247 : //
248 17 : void enableFesProjFileOutput() {fesproj_fileoutput_active_=true;}
249 : void disableFesFileProjOutput() {fesproj_fileoutput_active_=false;}
250 : bool isFesProjFileOutputActive() const {return fesproj_fileoutput_active_;}
251 : //
252 41 : void enableDynamicTargetDistFileOutput() {dynamic_targetdist_fileoutput_active_=true;}
253 : void disableDynamicTargetDistFileOutput() {dynamic_targetdist_fileoutput_active_=false;}
254 : bool isDynamicTargetDistFileOutputActive() const {return dynamic_targetdist_fileoutput_active_;}
255 : std::string getTargetDistFileFmt() const {return targetdist_file_fmt_;}
256 : std::string getTargetDistRestartFileFmt() const {return targetdist_restart_file_fmt_;}
257 : //
258 : void enableStaticTargetDistFileOutput() {static_targetdist_fileoutput_active_=true;}
259 46 : void disableStaticTargetDistFileOutput() {static_targetdist_fileoutput_active_=false;}
260 50 : bool isStaticTargetDistFileOutputActive() const {return static_targetdist_fileoutput_active_;}
261 : //
262 : std::vector< std::vector<std::string> > getProjectionArguments() const {return projection_args_;}
263 56 : std::vector<std::string> getProjectionArgument(unsigned int i) const {return projection_args_[i];}
264 122 : unsigned int getNumberOfProjectionArguments() const {return projection_args_.size();}
265 : //
266 : void setupBiasCutoff(const double, const double);
267 1498465 : bool biasCutoffActive() const {return bias_cutoff_active_;}
268 45 : double getBiasCutoffValue() const {return bias_cutoff_value_;}
269 498 : void setCurrentBiasMaxValue(const double max_value) {bias_current_max_value=max_value;}
270 : double getCurrentBiasMaxValue() const {return bias_current_max_value;}
271 : double getBiasCutoffSwitchingFunction(const double, double&) const;
272 : double getBiasCutoffSwitchingFunction(const double) const;
273 : void applyBiasCutoff(double&, std::vector<double>&) const;
274 : void applyBiasCutoff(double&, std::vector<double>&, std::vector<double>&) const;
275 : //
276 : std::unique_ptr<OFile> getOFile(const std::string& filename, const bool multi_sim_single_file=false, const bool enforce_backup=true);
277 : //
278 0 : virtual void setupBiasFileOutput() {};
279 0 : virtual void writeBiasToFile() {};
280 0 : virtual void resetBiasFileOutput() {};
281 : //
282 0 : virtual void setupFesFileOutput() {};
283 0 : virtual void writeFesToFile() {};
284 0 : virtual void resetFesFileOutput() {};
285 : //
286 0 : virtual void setupFesProjFileOutput() {};
287 0 : virtual void writeFesProjToFile() {};
288 0 : virtual void resetFesProjFileOutput() {};
289 : //
290 44 : virtual void setupTargetDistFileOutput() {};
291 0 : virtual void writeTargetDistToFile() {};
292 0 : virtual void resetTargetDistFileOutput() {};
293 : //
294 9 : virtual void setupTargetDistProjFileOutput() {};
295 0 : virtual void writeTargetDistProjToFile() {};
296 0 : virtual void resetTargetDistProjFileOutput() {};
297 : //
298 : void updateReweightFactor();
299 : virtual double calculateReweightFactor() const;
300 0 : bool isReweightFactorCalculated() const {return calc_reweightfactor_;}
301 : };
302 :
303 :
304 : inline
305 : size_t VesBias::getCoeffsIndex(const std::vector<unsigned int>& indices, const unsigned int coeffs_id) const {return coeffs_pntrs_[coeffs_id]->getIndex(indices);}
306 :
307 : inline
308 : std::vector<unsigned int> VesBias::getCoeffsIndices(const size_t index, const unsigned int coeffs_id) const {return coeffs_pntrs_[coeffs_id]->getIndices(index);}
309 :
310 : inline
311 3607712 : size_t VesBias::getHessianIndex(const size_t index1, const size_t index2, const unsigned int coeffs_id) const {return hessian_pntrs_[coeffs_id]->getMatrixIndex(index1,index2);}
312 :
313 :
314 : inline
315 23277 : double VesBias::getBeta() const {
316 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.");
317 23277 : return 1.0/kbt_;
318 : }
319 :
320 :
321 : inline
322 : std::string VesBias::getCurrentBiasOutputFilename(const std::string& suffix) const {
323 178 : return getCurrentOutputFilename(bias_filename_,suffix);
324 : }
325 :
326 :
327 : inline
328 : std::string VesBias::getCurrentFesOutputFilename(const std::string& suffix) const {
329 205 : return getCurrentOutputFilename(fes_filename_,suffix);
330 : }
331 :
332 :
333 : inline
334 : double VesBias::getBiasCutoffSwitchingFunction(const double bias) const {
335 : double dummy=0.0;
336 : return getBiasCutoffSwitchingFunction(bias,dummy);
337 : }
338 :
339 :
340 : inline
341 600 : void VesBias::applyBiasCutoff(double& bias, std::vector<double>& forces) const {
342 600 : std::vector<double> dummy(0);
343 600 : applyBiasCutoff(bias,forces,dummy);
344 600 : }
345 :
346 :
347 : inline
348 663 : void VesBias::applyBiasCutoff(double& bias, std::vector<double>& forces, std::vector<double>& coeffsderivs_values) const {
349 663 : double deriv_factor_sf=0.0;
350 663 : double value_sf = getBiasCutoffSwitchingFunction(bias,deriv_factor_sf);
351 663 : bias *= value_sf;
352 1326 : for(unsigned int i=0; i<forces.size(); i++) {
353 663 : forces[i] *= deriv_factor_sf;
354 : }
355 : //
356 1398 : for(unsigned int i=0; i<coeffsderivs_values.size(); i++) {
357 735 : coeffsderivs_values[i] *= deriv_factor_sf;
358 : }
359 663 : }
360 :
361 :
362 : inline
363 22810 : std::vector<double> VesBias::computeCovarianceFromAverages(const unsigned int c_id) const {
364 : size_t ncoeffs = numberOfCoeffs(c_id);
365 22810 : std::vector<double> covariance(sampled_cross_averages[c_id].size(),0.0);
366 : // diagonal part
367 1823850 : for(size_t i=0; i<ncoeffs; i++) {
368 : size_t midx = getHessianIndex(i,i,c_id);
369 1801040 : covariance[midx] = sampled_cross_averages[c_id][midx] - sampled_averages[c_id][i]*sampled_averages[c_id][i];
370 : }
371 22810 : if(!diagonal_hessian_) {
372 0 : for(size_t i=0; i<ncoeffs; i++) {
373 0 : for(size_t j=(i+1); j<ncoeffs; j++) {
374 : size_t midx = getHessianIndex(i,j,c_id);
375 0 : covariance[midx] = sampled_cross_averages[c_id][midx] - sampled_averages[c_id][i]*sampled_averages[c_id][j];
376 : }
377 : }
378 : }
379 22810 : return covariance;
380 : }
381 :
382 :
383 : template<class T>
384 180 : bool VesBias::parseMultipleValues(const std::string& keyword, std::vector<T>& values, unsigned int nvalues) {
385 180 : plumed_assert(nvalues>0);
386 180 : plumed_assert(values.size()==0);
387 : bool identical_values=false;
388 : //
389 180 : parseVector(keyword,values);
390 180 : if(values.size()==1 && nvalues>1) {
391 0 : values.resize(nvalues,values[0]);
392 : identical_values=true;
393 : }
394 180 : if(values.size()>0 && values.size()!=nvalues) {
395 0 : std::string s1; Tools::convert(nvalues,s1);
396 0 : plumed_merror("Error in " + keyword + " keyword: either give 1 common parameter value or " + s1 + " separate parameter values");
397 : }
398 180 : return identical_values;
399 : }
400 :
401 : template<class T>
402 90 : bool VesBias::parseMultipleValues(const std::string& keyword, std::vector<T>& values, unsigned int nvalues, const T& default_value) {
403 90 : bool identical_values = parseMultipleValues(keyword,values,nvalues);
404 90 : if(values.size()==0) {
405 0 : values.resize(nvalues,default_value);
406 : identical_values=true;
407 : }
408 90 : return identical_values;
409 : }
410 :
411 :
412 : }
413 : }
414 :
415 : #endif
|
