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