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 357 : void TargetDistribution::registerKeywords( Keywords& keys ) {
41 357 : Action::registerKeywords(keys);
42 1428 : 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 1071 : 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 1071 : 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 357 : }
46 :
47 :
48 342 : 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 342 : bias_cutoff_active_(false)
72 : {
73 : //
74 684 : if(keywords.exists("WELLTEMPERED_FACTOR")) {
75 253 : double welltempered_factor=0.0;
76 506 : parse("WELLTEMPERED_FACTOR",welltempered_factor);
77 : //
78 253 : if(welltempered_factor>0.0) {
79 12 : TargetDistModifer* pntr = new WellTemperedModifer(welltempered_factor);
80 6 : targetdist_modifer_pntrs_.push_back(pntr);
81 : }
82 247 : else if(welltempered_factor<0.0) {
83 0 : plumed_merror(getName()+": negative value in WELLTEMPERED_FACTOR does not make sense");
84 : }
85 : }
86 : //
87 684 : if(keywords.exists("SHIFT_TO_ZERO")) {
88 482 : parseFlag("SHIFT_TO_ZERO",shift_targetdist_to_zero_);
89 241 : 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 684 : if(keywords.exists("NORMALIZE")) {
96 215 : bool force_normalization=false;
97 430 : parseFlag("NORMALIZE",force_normalization);
98 215 : 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 342 : }
105 :
106 :
107 684 : TargetDistribution::~TargetDistribution() {
108 342 : if(targetdist_grid_pntr_!=NULL) {
109 342 : delete targetdist_grid_pntr_;
110 : }
111 342 : if(log_targetdist_grid_pntr_!=NULL) {
112 342 : delete log_targetdist_grid_pntr_;
113 : }
114 702 : for(unsigned int i=0; i<targetdist_modifer_pntrs_.size(); i++) {
115 6 : delete targetdist_modifer_pntrs_[i];
116 : }
117 342 : }
118 :
119 :
120 352 : double TargetDistribution::getBeta() const {
121 352 : plumed_massert(vesbias_pntr_!=NULL,"The VesBias has to be linked to use TargetDistribution::getBeta()");
122 352 : return vesbias_pntr_->getBeta();
123 : }
124 :
125 :
126 355 : void TargetDistribution::setDimension(const unsigned int dimension) {
127 355 : plumed_massert(dimension_==0,"setDimension: the dimension of the target distribution has already been set");
128 355 : dimension_=dimension;
129 355 : }
130 :
131 :
132 44 : void TargetDistribution::linkVesBias(VesBias* vesbias_pntr_in) {
133 44 : vesbias_pntr_ = vesbias_pntr_in;
134 44 : action_pntr_ = static_cast<Action*>(vesbias_pntr_in);
135 44 : }
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 35 : void TargetDistribution::linkFesGrid(Grid* fes_grid_pntr_in) {
154 35 : fes_grid_pntr_ = fes_grid_pntr_in;
155 35 : }
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 342 : 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 342 : if(getDimension()==0) {
177 73 : setDimension(arguments.size());
178 : }
179 : unsigned int dimension = getDimension();
180 342 : plumed_massert(arguments.size()==dimension,"TargetDistribution::setupGrids: mismatch between number of values given for grid parameters");
181 342 : plumed_massert(min.size()==dimension,"TargetDistribution::setupGrids: mismatch between number of values given for grid parameters");
182 342 : plumed_massert(max.size()==dimension,"TargetDistribution::setupGrids: mismatch between number of values given for grid parameters");
183 342 : plumed_massert(nbins.size()==dimension,"TargetDistribution::setupGrids: mismatch between number of values given for grid parameters");
184 342 : grid_args_=arguments;
185 684 : targetdist_grid_pntr_ = new Grid("targetdist",arguments,min,max,nbins,false,false);
186 684 : log_targetdist_grid_pntr_ = new Grid("log_targetdist",arguments,min,max,nbins,false,false);
187 342 : setupAdditionalGrids(arguments,min,max,nbins);
188 342 : }
189 :
190 :
191 308 : void TargetDistribution::calculateStaticDistributionGrid() {
192 308 : if(static_grid_calculated && !bias_cutoff_active_) {return;}
193 : // plumed_massert(isStatic(),"this should only be used for static distributions");
194 288 : plumed_massert(targetdist_grid_pntr_!=NULL,"the grids have not been setup using setupGrids");
195 288 : plumed_massert(log_targetdist_grid_pntr_!=NULL,"the grids have not been setup using setupGrids");
196 899422 : for(Grid::index_t l=0; l<targetdist_grid_pntr_->getSize(); l++)
197 : {
198 449567 : std::vector<double> argument = targetdist_grid_pntr_->getPoint(l);
199 449567 : double value = getValue(argument);
200 449567 : targetdist_grid_pntr_->setValue(l,value);
201 449567 : log_targetdist_grid_pntr_->setValue(l,-std::log(value));
202 : }
203 288 : log_targetdist_grid_pntr_->setMinToZero();
204 288 : static_grid_calculated = true;
205 : }
206 :
207 :
208 812 : double TargetDistribution::integrateGrid(const Grid* grid_pntr) {
209 2436 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(grid_pntr);
210 : double sum = 0.0;
211 5053712 : for(Grid::index_t l=0; l<grid_pntr->getSize(); l++) {
212 2526450 : sum += integration_weights[l]*grid_pntr->getValue(l);
213 : }
214 812 : return sum;
215 : }
216 :
217 :
218 78 : double TargetDistribution::normalizeGrid(Grid* grid_pntr) {
219 78 : double normalization = TargetDistribution::integrateGrid(grid_pntr);
220 78 : grid_pntr->scaleAllValuesAndDerivatives(1.0/normalization);
221 78 : return normalization;
222 : }
223 :
224 :
225 28 : Grid TargetDistribution::getMarginalDistributionGrid(Grid* grid_pntr, const std::vector<std::string>& args) {
226 28 : plumed_massert(grid_pntr->getDimension()>1,"doesn't make sense calculating the marginal distribution for a one-dimensional distribution");
227 28 : 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 56 : std::vector<std::string> argnames = grid_pntr->getArgNames();
230 28 : std::vector<unsigned int> args_index(0);
231 168 : for(unsigned int i=0; i<argnames.size(); i++) {
232 280 : for(unsigned int l=0; l<args.size(); l++) {
233 112 : if(argnames[i]==args[l]) {args_index.push_back(i);}
234 : }
235 : }
236 28 : plumed_massert(args.size()==args_index.size(),"getMarginalDistributionGrid: problem with the arguments of the marginal");
237 : //
238 28 : MarginalWeight* Pw = new MarginalWeight();
239 28 : Grid proj_grid = grid_pntr->project(args,Pw);
240 28 : 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 28 : double intVol = grid_pntr->getBinVolume();
245 140 : for(unsigned int l=0; l<args_index.size(); l++) {
246 84 : intVol/=grid_pntr->getDx()[args_index[l]];
247 : }
248 28 : proj_grid.scaleAllValuesAndDerivatives(intVol);
249 : //
250 28 : 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 713 : void TargetDistribution::updateTargetDist() {
260 : //
261 713 : updateGrid();
262 : //
263 1444 : for(unsigned int i=0; i<targetdist_modifer_pntrs_.size(); i++) {
264 6 : applyTargetDistModiferToGrid(targetdist_modifer_pntrs_[i]);
265 : }
266 : //
267 713 : if(bias_cutoff_active_) {updateBiasCutoffForTargetDistGrid();}
268 : //
269 713 : if(shift_targetdist_to_zero_ && !(bias_cutoff_active_)) {setMinimumOfTargetDistGridToZero();}
270 713 : if(force_normalization_ && !(bias_cutoff_active_) ) {normalizeTargetDistGrid();}
271 : //
272 : // if(check_normalization_ && !force_normalization_ && !shift_targetdist_to_zero_){
273 713 : if(check_normalization_ && !(bias_cutoff_active_)) {
274 614 : double normalization = integrateGrid(targetdist_grid_pntr_);
275 : const double normalization_thrshold = 0.1;
276 614 : 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 713 : if(check_nonnegative_) {
284 : const double nonnegative_thrshold = -0.02;
285 710 : double grid_min_value = targetdist_grid_pntr_->getMinValue();
286 710 : 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 713 : if(check_nan_inf_) {checkNanAndInf();}
294 : //
295 713 : }
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 5224 : 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 42418 : 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 22817 : 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 16 : 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 : Grid* restart_grid = Grid::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 8 : VesTools::copyGridValues(restart_grid,targetdist_grid_pntr_);
374 8 : updateLogTargetDistGrid();
375 8 : delete restart_grid;
376 8 : }
377 :
378 1 : void TargetDistribution::clearLogTargetDistGrid() {
379 1 : log_targetdist_grid_pntr_->clear();
380 1 : }
381 :
382 :
383 0 : void TargetDistribution::checkNanAndInf() {
384 0 : for(Grid::index_t l=0; l<targetdist_grid_pntr_->getSize(); l++)
385 : {
386 0 : double value = targetdist_grid_pntr_->getValue(l);
387 0 : if(std::isnan(value) || std::isinf(value)) {
388 0 : std::string vs; Tools::convert(value,vs);
389 0 : std::vector<double> p = targetdist_grid_pntr_->getPoint(l);
390 0 : std::string ps; Tools::convert(p[0],ps);
391 0 : ps = "(" + ps;
392 0 : for(unsigned int k=1; k<p.size(); k++) {
393 0 : std::string t1; Tools::convert(p[k],t1);
394 0 : ps = ps + "," + t1;
395 : }
396 0 : ps = ps + ")";
397 0 : plumed_merror(getName()+": problem with target distribution, the value at " + ps + " is " + vs);
398 : }
399 : }
400 0 : }
401 :
402 : }
403 : }
|
