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 "TargetDistribution.h"
24 : #include "TargetDistModifer.h"
25 :
26 : #include "VesBias.h"
27 : #include "GridIntegrationWeights.h"
28 : #include "VesTools.h"
29 :
30 : #include "core/Value.h"
31 : #include "tools/Grid.h"
32 : #include "tools/File.h"
33 : #include "tools/Keywords.h"
34 :
35 : #include "GridProjWeights.h"
36 :
37 : namespace PLMD {
38 : namespace ves {
39 :
40 400 : void TargetDistribution::registerKeywords( Keywords& keys ) {
41 400 : Action::registerKeywords(keys);
42 1600 : keys.reserve("optional","WELLTEMPERED_FACTOR","Broaden the target distribution such that it is taken as [p(s)]^(1/\\f$\\gamma\\f$) where \\f$\\gamma\\f$ is the well tempered factor given here. If this option is active the distribution will be automatically normalized.");
43 1200 : keys.reserveFlag("SHIFT_TO_ZERO",false,"Shift the minimum value of the target distribution to zero. This can for example be used to avoid negative values in the target distribution. If this option is active the distribution will be automatically normalized.");
44 1200 : keys.reserveFlag("NORMALIZE",false,"Renormalized the target distribution over the intervals on which it is defined to make sure that it is properly normalized to 1. In most cases this should not be needed as the target distributions should be normalized. The code will issue a warning (but still run) if this is needed for some reason.");
45 400 : }
46 :
47 :
48 382 : TargetDistribution::TargetDistribution(const ActionOptions&ao):
49 : Action(ao),
50 : type_(static_targetdist),
51 : force_normalization_(false),
52 : check_normalization_(true),
53 : check_nonnegative_(true),
54 : check_nan_inf_(false),
55 : shift_targetdist_to_zero_(false),
56 : dimension_(0),
57 : grid_args_(0),
58 : targetdist_grid_pntr_(NULL),
59 : log_targetdist_grid_pntr_(NULL),
60 : targetdist_modifer_pntrs_(0),
61 : action_pntr_(NULL),
62 : vesbias_pntr_(NULL),
63 : needs_bias_grid_(false),
64 : needs_bias_withoutcutoff_grid_(false),
65 : needs_fes_grid_(false),
66 : bias_grid_pntr_(NULL),
67 : bias_withoutcutoff_grid_pntr_(NULL),
68 : fes_grid_pntr_(NULL),
69 : static_grid_calculated(false),
70 : allow_bias_cutoff_(true),
71 382 : bias_cutoff_active_(false)
72 : {
73 : //
74 764 : if(keywords.exists("WELLTEMPERED_FACTOR")) {
75 278 : double welltempered_factor=0.0;
76 556 : parse("WELLTEMPERED_FACTOR",welltempered_factor);
77 : //
78 278 : if(welltempered_factor>0.0) {
79 12 : TargetDistModifer* pntr = new WellTemperedModifer(welltempered_factor);
80 6 : targetdist_modifer_pntrs_.push_back(pntr);
81 : }
82 272 : else if(welltempered_factor<0.0) {
83 0 : plumed_merror(getName()+": negative value in WELLTEMPERED_FACTOR does not make sense");
84 : }
85 : }
86 : //
87 764 : if(keywords.exists("SHIFT_TO_ZERO")) {
88 532 : parseFlag("SHIFT_TO_ZERO",shift_targetdist_to_zero_);
89 266 : if(shift_targetdist_to_zero_) {
90 3 : if(bias_cutoff_active_) {plumed_merror(getName()+": using SHIFT_TO_ZERO with bias cutoff is not allowed.");}
91 3 : check_nonnegative_=false;
92 : }
93 : }
94 : //
95 764 : if(keywords.exists("NORMALIZE")) {
96 240 : bool force_normalization=false;
97 480 : parseFlag("NORMALIZE",force_normalization);
98 240 : if(force_normalization) {
99 3 : if(shift_targetdist_to_zero_) {plumed_merror(getName()+" with label "+getLabel()+": using NORMALIZE with SHIFT_TO_ZERO is not needed, the target distribution will be automatically normalized.");}
100 : setForcedNormalization();
101 : }
102 : }
103 :
104 382 : }
105 :
106 :
107 764 : TargetDistribution::~TargetDistribution() {
108 382 : if(targetdist_grid_pntr_!=NULL) {
109 382 : delete targetdist_grid_pntr_;
110 : }
111 382 : if(log_targetdist_grid_pntr_!=NULL) {
112 382 : delete log_targetdist_grid_pntr_;
113 : }
114 394 : for(unsigned int i=0; i<targetdist_modifer_pntrs_.size(); i++) {
115 6 : delete targetdist_modifer_pntrs_[i];
116 : }
117 382 : }
118 :
119 :
120 401 : double TargetDistribution::getBeta() const {
121 401 : plumed_massert(vesbias_pntr_!=NULL,"The VesBias has to be linked to use TargetDistribution::getBeta()");
122 401 : return vesbias_pntr_->getBeta();
123 : }
124 :
125 :
126 395 : void TargetDistribution::setDimension(const unsigned int dimension) {
127 395 : plumed_massert(dimension_==0,"setDimension: the dimension of the target distribution has already been set");
128 395 : dimension_=dimension;
129 395 : }
130 :
131 :
132 50 : void TargetDistribution::linkVesBias(VesBias* vesbias_pntr_in) {
133 50 : vesbias_pntr_ = vesbias_pntr_in;
134 50 : action_pntr_ = static_cast<Action*>(vesbias_pntr_in);
135 50 : }
136 :
137 :
138 0 : void TargetDistribution::linkAction(Action* action_pntr_in) {
139 0 : action_pntr_ = action_pntr_in;
140 0 : }
141 :
142 :
143 0 : void TargetDistribution::linkBiasGrid(Grid* bias_grid_pntr_in) {
144 0 : bias_grid_pntr_ = bias_grid_pntr_in;
145 0 : }
146 :
147 :
148 3 : void TargetDistribution::linkBiasWithoutCutoffGrid(Grid* bias_withoutcutoff_grid_pntr_in) {
149 3 : bias_withoutcutoff_grid_pntr_ = bias_withoutcutoff_grid_pntr_in;
150 3 : }
151 :
152 :
153 41 : void TargetDistribution::linkFesGrid(Grid* fes_grid_pntr_in) {
154 41 : fes_grid_pntr_ = fes_grid_pntr_in;
155 41 : }
156 :
157 :
158 3 : void TargetDistribution::setupBiasCutoff() {
159 3 : if(!allow_bias_cutoff_) {
160 0 : plumed_merror(getName()+" with label "+getLabel()+": this target distribution does not support a bias cutoff");
161 : }
162 3 : if(targetdist_modifer_pntrs_.size()>0) {
163 0 : plumed_merror(getName()+" with label "+getLabel()+": using a bias cutoff with a target distribution modifer like WELLTEMPERED_FACTOR is not allowed");
164 : }
165 3 : bias_cutoff_active_=true;
166 : setBiasWithoutCutoffGridNeeded();
167 : setDynamic();
168 : // as the p(s) includes the derivative factor so normalization
169 : // check can be misleading
170 3 : check_normalization_=false;
171 3 : force_normalization_=false;
172 3 : }
173 :
174 :
175 382 : void TargetDistribution::setupGrids(const std::vector<Value*>& arguments, const std::vector<std::string>& min, const std::vector<std::string>& max, const std::vector<unsigned int>& nbins) {
176 382 : if(getDimension()==0) {
177 78 : setDimension(arguments.size());
178 : }
179 : unsigned int dimension = getDimension();
180 382 : plumed_massert(arguments.size()==dimension,"TargetDistribution::setupGrids: mismatch between number of values given for grid parameters");
181 382 : plumed_massert(min.size()==dimension,"TargetDistribution::setupGrids: mismatch between number of values given for grid parameters");
182 382 : plumed_massert(max.size()==dimension,"TargetDistribution::setupGrids: mismatch between number of values given for grid parameters");
183 382 : plumed_massert(nbins.size()==dimension,"TargetDistribution::setupGrids: mismatch between number of values given for grid parameters");
184 382 : grid_args_=arguments;
185 764 : targetdist_grid_pntr_ = new Grid("targetdist",arguments,min,max,nbins,false,false);
186 764 : log_targetdist_grid_pntr_ = new Grid("log_targetdist",arguments,min,max,nbins,false,false);
187 382 : setupAdditionalGrids(arguments,min,max,nbins);
188 382 : }
189 :
190 :
191 342 : void TargetDistribution::calculateStaticDistributionGrid() {
192 684 : if(static_grid_calculated && !bias_cutoff_active_) {return;}
193 : // plumed_massert(isStatic(),"this should only be used for static distributions");
194 322 : plumed_massert(targetdist_grid_pntr_!=NULL,"the grids have not been setup using setupGrids");
195 322 : plumed_massert(log_targetdist_grid_pntr_!=NULL,"the grids have not been setup using setupGrids");
196 459801 : for(Grid::index_t l=0; l<targetdist_grid_pntr_->getSize(); l++)
197 : {
198 459801 : std::vector<double> argument = targetdist_grid_pntr_->getPoint(l);
199 459801 : double value = getValue(argument);
200 459801 : targetdist_grid_pntr_->setValue(l,value);
201 459801 : log_targetdist_grid_pntr_->setValue(l,-std::log(value));
202 : }
203 322 : log_targetdist_grid_pntr_->setMinToZero();
204 322 : static_grid_calculated = true;
205 : }
206 :
207 :
208 888 : double TargetDistribution::integrateGrid(const Grid* grid_pntr) {
209 2664 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(grid_pntr);
210 : double sum = 0.0;
211 2584272 : for(Grid::index_t l=0; l<grid_pntr->getSize(); l++) {
212 2583384 : sum += integration_weights[l]*grid_pntr->getValue(l);
213 : }
214 888 : return sum;
215 : }
216 :
217 :
218 88 : double TargetDistribution::normalizeGrid(Grid* grid_pntr) {
219 88 : double normalization = TargetDistribution::integrateGrid(grid_pntr);
220 88 : grid_pntr->scaleAllValuesAndDerivatives(1.0/normalization);
221 88 : return normalization;
222 : }
223 :
224 :
225 26 : Grid TargetDistribution::getMarginalDistributionGrid(Grid* grid_pntr, const std::vector<std::string>& args) {
226 26 : plumed_massert(grid_pntr->getDimension()>1,"doesn't make sense calculating the marginal distribution for a one-dimensional distribution");
227 26 : plumed_massert(args.size()<grid_pntr->getDimension(),"the number of arguments for the marginal distribution should be less than the dimension of the full distribution");
228 : //
229 26 : std::vector<std::string> argnames = grid_pntr->getArgNames();
230 26 : std::vector<unsigned int> args_index(0);
231 156 : for(unsigned int i=0; i<argnames.size(); i++) {
232 156 : for(unsigned int l=0; l<args.size(); l++) {
233 104 : if(argnames[i]==args[l]) {args_index.push_back(i);}
234 : }
235 : }
236 26 : plumed_massert(args.size()==args_index.size(),"getMarginalDistributionGrid: problem with the arguments of the marginal");
237 : //
238 26 : MarginalWeight* Pw = new MarginalWeight();
239 26 : Grid proj_grid = grid_pntr->project(args,Pw);
240 26 : delete Pw;
241 : //
242 : // scale with the bin volume used for the integral such that the
243 : // marginals are proberly normalized to 1.0
244 26 : double intVol = grid_pntr->getBinVolume();
245 104 : for(unsigned int l=0; l<args_index.size(); l++) {
246 78 : intVol/=grid_pntr->getDx()[args_index[l]];
247 : }
248 26 : proj_grid.scaleAllValuesAndDerivatives(intVol);
249 : //
250 26 : return proj_grid;
251 : }
252 :
253 :
254 8 : Grid TargetDistribution::getMarginal(const std::vector<std::string>& args) {
255 8 : return TargetDistribution::getMarginalDistributionGrid(targetdist_grid_pntr_,args);
256 : }
257 :
258 :
259 789 : void TargetDistribution::updateTargetDist() {
260 : //
261 789 : updateGrid();
262 : //
263 1590 : for(unsigned int i=0; i<targetdist_modifer_pntrs_.size(); i++) {
264 6 : applyTargetDistModiferToGrid(targetdist_modifer_pntrs_[i]);
265 : }
266 : //
267 789 : if(bias_cutoff_active_) {updateBiasCutoffForTargetDistGrid();}
268 : //
269 789 : if(shift_targetdist_to_zero_ && !(bias_cutoff_active_)) {setMinimumOfTargetDistGridToZero();}
270 789 : if(force_normalization_ && !(bias_cutoff_active_) ) {normalizeTargetDistGrid();}
271 : //
272 : // if(check_normalization_ && !force_normalization_ && !shift_targetdist_to_zero_){
273 789 : if(check_normalization_ && !(bias_cutoff_active_)) {
274 680 : double normalization = integrateGrid(targetdist_grid_pntr_);
275 : const double normalization_thrshold = 0.1;
276 680 : if(normalization < 1.0-normalization_thrshold || normalization > 1.0+normalization_thrshold) {
277 3 : std::string norm_str; Tools::convert(normalization,norm_str);
278 6 : std::string msg = "the target distribution grid is not proberly normalized, integrating over the grid gives: " + norm_str + " - You can avoid this problem by using the NORMALIZE keyword";
279 3 : warning(msg);
280 : }
281 : }
282 : //
283 789 : if(check_nonnegative_) {
284 : const double nonnegative_thrshold = -0.02;
285 786 : double grid_min_value = targetdist_grid_pntr_->getMinValue();
286 786 : if(grid_min_value<nonnegative_thrshold) {
287 0 : std::string grid_min_value_str; Tools::convert(grid_min_value,grid_min_value_str);
288 0 : std::string msg = "the target distribution grid has negative values, the lowest value is: " + grid_min_value_str + " - You can avoid this problem by using the SHIFT_TO_ZERO keyword";
289 0 : warning(msg);
290 : }
291 : }
292 : //
293 789 : if(check_nan_inf_) {checkNanAndInf();}
294 : //
295 789 : }
296 :
297 :
298 24 : void TargetDistribution::updateBiasCutoffForTargetDistGrid() {
299 24 : plumed_massert(vesbias_pntr_!=NULL,"The VesBias has to be linked to use updateBiasCutoffForTargetDistGrid()");
300 24 : plumed_massert(vesbias_pntr_->biasCutoffActive(),"updateBiasCutoffForTargetDistGrid() should only be used if the bias cutoff is active");
301 : // plumed_massert(targetdist_grid_pntr_!=NULL,"the grids have not been setup using setupGrids");
302 : // plumed_massert(log_targetdist_grid_pntr_!=NULL,"the grids have not been setup using setupGrids");
303 24 : plumed_massert(getBiasWithoutCutoffGridPntr()!=NULL,"the bias without cutoff grid has to be linked");
304 : //
305 72 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(targetdist_grid_pntr_);
306 : double norm = 0.0;
307 2624 : for(Grid::index_t l=0; l<targetdist_grid_pntr_->getSize(); l++)
308 : {
309 2600 : double value = targetdist_grid_pntr_->getValue(l);
310 2600 : double bias = getBiasWithoutCutoffGridPntr()->getValue(l);
311 2600 : double deriv_factor_swf = 0.0;
312 2600 : double swf = vesbias_pntr_->getBiasCutoffSwitchingFunction(bias,deriv_factor_swf);
313 : // this comes from the p(s)
314 2600 : value *= swf;
315 2600 : norm += integration_weights[l]*value;
316 : // this comes from the derivative of V(s)
317 2600 : value *= deriv_factor_swf;
318 2600 : targetdist_grid_pntr_->setValue(l,value);
319 : // double log_value = log_targetdist_grid_pntr_->getValue(l) - std::log(swf);
320 : // log_targetdist_grid_pntr_->setValue(l,log_value);
321 : }
322 24 : targetdist_grid_pntr_->scaleAllValuesAndDerivatives(1.0/norm);
323 : // log_targetdist_grid_pntr_->setMinToZero();
324 24 : }
325 :
326 6 : void TargetDistribution::applyTargetDistModiferToGrid(TargetDistModifer* modifer_pntr) {
327 : // plumed_massert(targetdist_grid_pntr_!=NULL,"the grids have not been setup using setupGrids");
328 : // plumed_massert(log_targetdist_grid_pntr_!=NULL,"the grids have not been setup using setupGrids");
329 : //
330 18 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(targetdist_grid_pntr_);
331 : double norm = 0.0;
332 21212 : for(Grid::index_t l=0; l<targetdist_grid_pntr_->getSize(); l++)
333 : {
334 21206 : double value = targetdist_grid_pntr_->getValue(l);
335 21206 : std::vector<double> cv_values = targetdist_grid_pntr_->getPoint(l);
336 21206 : value = modifer_pntr->getModifedTargetDistValue(value,cv_values);
337 21206 : norm += integration_weights[l]*value;
338 21206 : targetdist_grid_pntr_->setValue(l,value);
339 21206 : log_targetdist_grid_pntr_->setValue(l,-std::log(value));
340 : }
341 6 : targetdist_grid_pntr_->scaleAllValuesAndDerivatives(1.0/norm);
342 6 : log_targetdist_grid_pntr_->setMinToZero();
343 6 : }
344 :
345 :
346 11 : void TargetDistribution::updateLogTargetDistGrid() {
347 11414 : for(Grid::index_t l=0; l<targetdist_grid_pntr_->getSize(); l++)
348 : {
349 11403 : log_targetdist_grid_pntr_->setValue(l,-std::log(targetdist_grid_pntr_->getValue(l)));
350 : }
351 11 : log_targetdist_grid_pntr_->setMinToZero();
352 11 : }
353 :
354 :
355 3 : void TargetDistribution::setMinimumOfTargetDistGridToZero() {
356 3 : targetDistGrid().setMinToZero();
357 3 : normalizeTargetDistGrid();
358 3 : updateLogTargetDistGrid();
359 3 : }
360 :
361 :
362 8 : void TargetDistribution::readInRestartTargetDistGrid(const std::string& grid_fname) {
363 8 : plumed_massert(isDynamic(),"this should only be used for dynamically updated target distributions!");
364 8 : IFile gridfile;
365 8 : if(!gridfile.FileExist(grid_fname)) {
366 0 : plumed_merror(getName()+": problem with reading previous target distribution when restarting, cannot find file " + grid_fname);
367 : }
368 8 : gridfile.open(grid_fname);
369 16 : std::unique_ptr<GridBase> restart_grid = GridBase::create("targetdist",grid_args_,gridfile,false,false,false);
370 8 : if(restart_grid->getSize()!=targetdist_grid_pntr_->getSize()) {
371 0 : plumed_merror(getName()+": problem with reading previous target distribution when restarting, the grid is not of the correct size!");
372 : }
373 16 : VesTools::copyGridValues(restart_grid.get(),targetdist_grid_pntr_);
374 16 : updateLogTargetDistGrid();
375 8 : }
376 :
377 1 : void TargetDistribution::clearLogTargetDistGrid() {
378 1 : log_targetdist_grid_pntr_->clear();
379 1 : }
380 :
381 :
382 0 : void TargetDistribution::checkNanAndInf() {
383 0 : for(Grid::index_t l=0; l<targetdist_grid_pntr_->getSize(); l++)
384 : {
385 0 : double value = targetdist_grid_pntr_->getValue(l);
386 0 : if(std::isnan(value) || std::isinf(value)) {
387 0 : std::string vs; Tools::convert(value,vs);
388 0 : std::vector<double> p = targetdist_grid_pntr_->getPoint(l);
389 0 : std::string ps; Tools::convert(p[0],ps);
390 0 : ps = "(" + ps;
391 0 : for(unsigned int k=1; k<p.size(); k++) {
392 0 : std::string t1; Tools::convert(p[k],t1);
393 0 : ps = ps + "," + t1;
394 : }
395 0 : ps = ps + ")";
396 0 : plumed_merror(getName()+": problem with target distribution, the value at " + ps + " is " + vs);
397 : }
398 : }
399 0 : }
400 :
401 : }
402 : }
|
