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 356 : void TargetDistribution::registerKeywords( Keywords& keys ) {
41 356 : Action::registerKeywords(keys);
42 356 : 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 356 : 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 356 : 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 356 : }
46 :
47 :
48 341 : 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 341 : bias_cutoff_active_(false)
72 : {
73 : //
74 341 : if(keywords.exists("WELLTEMPERED_FACTOR")) {
75 253 : double welltempered_factor=0.0;
76 253 : parse("WELLTEMPERED_FACTOR",welltempered_factor);
77 : //
78 253 : if(welltempered_factor>0.0) {
79 6 : 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 341 : if(keywords.exists("SHIFT_TO_ZERO")) {
88 241 : 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 341 : if(keywords.exists("NORMALIZE")) {
96 215 : bool force_normalization=false;
97 215 : 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 3 : setForcedNormalization();
101 : }
102 : }
103 :
104 341 : }
105 :
106 :
107 682 : TargetDistribution::~TargetDistribution() {
108 341 : if(targetdist_grid_pntr_!=NULL) {
109 341 : delete targetdist_grid_pntr_;
110 : }
111 341 : if(log_targetdist_grid_pntr_!=NULL) {
112 341 : delete log_targetdist_grid_pntr_;
113 : }
114 347 : for(unsigned int i=0; i<targetdist_modifer_pntrs_.size(); i++) {
115 6 : delete targetdist_modifer_pntrs_[i];
116 : }
117 341 : }
118 :
119 :
120 341 : double TargetDistribution::getBeta() const {
121 341 : plumed_massert(vesbias_pntr_!=NULL,"The VesBias has to be linked to use TargetDistribution::getBeta()");
122 341 : return vesbias_pntr_->getBeta();
123 : }
124 :
125 :
126 354 : void TargetDistribution::setDimension(const unsigned int dimension) {
127 354 : plumed_massert(dimension_==0,"setDimension: the dimension of the target distribution has already been set");
128 354 : dimension_=dimension;
129 354 : }
130 :
131 :
132 43 : void TargetDistribution::linkVesBias(VesBias* vesbias_pntr_in) {
133 43 : vesbias_pntr_ = vesbias_pntr_in;
134 43 : action_pntr_ = static_cast<Action*>(vesbias_pntr_in);
135 43 : }
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 34 : void TargetDistribution::linkFesGrid(Grid* fes_grid_pntr_in) {
154 34 : fes_grid_pntr_ = fes_grid_pntr_in;
155 34 : }
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 3 : setBiasWithoutCutoffGridNeeded();
167 3 : 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 341 : 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 341 : if(getDimension()==0) {
177 72 : setDimension(arguments.size());
178 : }
179 341 : unsigned int dimension = getDimension();
180 341 : plumed_massert(arguments.size()==dimension,"TargetDistribution::setupGrids: mismatch between number of values given for grid parameters");
181 341 : plumed_massert(min.size()==dimension,"TargetDistribution::setupGrids: mismatch between number of values given for grid parameters");
182 341 : plumed_massert(max.size()==dimension,"TargetDistribution::setupGrids: mismatch between number of values given for grid parameters");
183 341 : plumed_massert(nbins.size()==dimension,"TargetDistribution::setupGrids: mismatch between number of values given for grid parameters");
184 341 : grid_args_=arguments;
185 341 : targetdist_grid_pntr_ = new Grid("targetdist",arguments,min,max,nbins,false,false);
186 341 : log_targetdist_grid_pntr_ = new Grid("log_targetdist",arguments,min,max,nbins,false,false);
187 341 : setupAdditionalGrids(arguments,min,max,nbins);
188 341 : }
189 :
190 :
191 308 : void TargetDistribution::calculateStaticDistributionGrid() {
192 616 : 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 449855 : 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 449567 : }
203 288 : log_targetdist_grid_pntr_->setMinToZero();
204 288 : static_grid_calculated = true;
205 : }
206 :
207 :
208 801 : double TargetDistribution::integrateGrid(const Grid* grid_pntr) {
209 801 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(grid_pntr);
210 801 : double sum = 0.0;
211 2417251 : for(Grid::index_t l=0; l<grid_pntr->getSize(); l++) {
212 2416450 : sum += integration_weights[l]*grid_pntr->getValue(l);
213 : }
214 801 : 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 22 : Grid TargetDistribution::getMarginalDistributionGrid(Grid* grid_pntr, const std::vector<std::string>& args) {
226 22 : plumed_massert(grid_pntr->getDimension()>1,"doesn't make sense calculating the marginal distribution for a one-dimensional distribution");
227 22 : 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 22 : std::vector<std::string> argnames = grid_pntr->getArgNames();
230 44 : std::vector<unsigned int> args_index(0);
231 66 : for(unsigned int i=0; i<argnames.size(); i++) {
232 88 : for(unsigned int l=0; l<args.size(); l++) {
233 44 : if(argnames[i]==args[l]) {args_index.push_back(i);}
234 : }
235 : }
236 22 : plumed_massert(args.size()==args_index.size(),"getMarginalDistributionGrid: problem with the arguments of the marginal");
237 : //
238 22 : MarginalWeight* Pw = new MarginalWeight();
239 22 : Grid proj_grid = grid_pntr->project(args,Pw);
240 22 : 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 22 : double intVol = grid_pntr->getBinVolume();
245 44 : for(unsigned int l=0; l<args_index.size(); l++) {
246 22 : intVol/=grid_pntr->getDx()[args_index[l]];
247 : }
248 22 : proj_grid.scaleAllValuesAndDerivatives(intVol);
249 : //
250 44 : 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 702 : void TargetDistribution::updateTargetDist() {
260 : //
261 702 : updateGrid();
262 : //
263 708 : for(unsigned int i=0; i<targetdist_modifer_pntrs_.size(); i++) {
264 6 : applyTargetDistModiferToGrid(targetdist_modifer_pntrs_[i]);
265 : }
266 : //
267 702 : if(bias_cutoff_active_) {updateBiasCutoffForTargetDistGrid();}
268 : //
269 702 : if(shift_targetdist_to_zero_ && !(bias_cutoff_active_)) {setMinimumOfTargetDistGridToZero();}
270 702 : if(force_normalization_ && !(bias_cutoff_active_) ) {normalizeTargetDistGrid();}
271 : //
272 : // if(check_normalization_ && !force_normalization_ && !shift_targetdist_to_zero_){
273 702 : if(check_normalization_ && !(bias_cutoff_active_)) {
274 603 : double normalization = integrateGrid(targetdist_grid_pntr_);
275 603 : const double normalization_thrshold = 0.1;
276 603 : 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 6 : warning(msg);
280 : }
281 : }
282 : //
283 702 : if(check_nonnegative_) {
284 699 : const double nonnegative_thrshold = -0.02;
285 699 : double grid_min_value = targetdist_grid_pntr_->getMinValue();
286 699 : 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 702 : if(check_nan_inf_) {checkNanAndInf();}
294 : //
295 702 : }
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 24 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(targetdist_grid_pntr_);
306 24 : 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 6 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(targetdist_grid_pntr_);
331 6 : 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 21206 : }
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 8 : 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 0 : }
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 : }
|
