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