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         475 : void TargetDistribution::registerKeywords( Keywords& keys ) {
      41         475 :   Action::registerKeywords(keys);
      42         475 :   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         950 :   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         950 :   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         475 : }
      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         407 :   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         348 :   plumed_massert(targetdist_grid_pntr_,"the grids have not been setup using setupGrids");
     184         348 :   plumed_massert(log_targetdist_grid_pntr_,"the grids have not been setup using setupGrids");
     185      467955 :   for(Grid::index_t l=0; l<targetdist_grid_pntr_->getSize(); l++) {
     186      467607 :     std::vector<double> argument = targetdist_grid_pntr_->getPoint(l);
     187      467607 :     double value = getValue(argument);
     188      467607 :     targetdist_grid_pntr_->setValue(l,value);
     189      467607 :     log_targetdist_grid_pntr_->setValue(l,-std::log(value));
     190             :   }
     191         348 :   log_targetdist_grid_pntr_->setMinToZero();
     192         348 :   static_grid_calculated = true;
     193             : }
     194             : 
     195             : 
     196         901 : double TargetDistribution::integrateGrid(const Grid* grid_pntr) {
     197        1802 :   std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(grid_pntr);
     198             :   double sum = 0.0;
     199     2579140 :   for(Grid::index_t l=0; l<grid_pntr->getSize(); l++) {
     200     2578239 :     sum += integration_weights[l]*grid_pntr->getValue(l);
     201             :   }
     202         901 :   return sum;
     203             : }
     204             : 
     205             : 
     206          90 : double TargetDistribution::normalizeGrid(Grid* grid_pntr) {
     207          90 :   double normalization = TargetDistribution::integrateGrid(grid_pntr);
     208          90 :   grid_pntr->scaleAllValuesAndDerivatives(1.0/normalization);
     209          90 :   return normalization;
     210             : }
     211             : 
     212             : 
     213          28 : Grid TargetDistribution::getMarginalDistributionGrid(Grid* grid_pntr, const std::vector<std::string>& args) {
     214          28 :   plumed_massert(grid_pntr->getDimension()>1,"doesn't make sense calculating the marginal distribution for a one-dimensional distribution");
     215          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");
     216             :   //
     217          28 :   std::vector<std::string> argnames = grid_pntr->getArgNames();
     218          28 :   std::vector<unsigned int> args_index(0);
     219          84 :   for(unsigned int i=0; i<argnames.size(); i++) {
     220         112 :     for(unsigned int l=0; l<args.size(); l++) {
     221          56 :       if(argnames[i]==args[l]) {
     222          28 :         args_index.push_back(i);
     223             :       }
     224             :     }
     225             :   }
     226          28 :   plumed_massert(args.size()==args_index.size(),"getMarginalDistributionGrid: problem with the arguments of the marginal");
     227             :   //
     228             :   auto Pw = Tools::make_unique<MarginalWeight>();
     229          28 :   Grid proj_grid = grid_pntr->project(args,Pw.get());
     230             :   Pw.reset();
     231             :   //
     232             :   // scale with the bin volume used for the integral such that the
     233             :   // marginals are properly normalized to 1.0
     234          28 :   double intVol = grid_pntr->getBinVolume();
     235          56 :   for(unsigned int l=0; l<args_index.size(); l++) {
     236          28 :     intVol/=grid_pntr->getDx()[args_index[l]];
     237             :   }
     238          28 :   proj_grid.scaleAllValuesAndDerivatives(intVol);
     239             :   //
     240          28 :   return proj_grid;
     241          56 : }
     242             : 
     243             : 
     244           8 : Grid TargetDistribution::getMarginal(const std::vector<std::string>& args) {
     245           8 :   return TargetDistribution::getMarginalDistributionGrid(targetdist_grid_pntr_.get(),args);
     246             : }
     247             : 
     248             : 
     249         802 : void TargetDistribution::updateTargetDist() {
     250             :   //
     251         802 :   updateGrid();
     252             :   //
     253         808 :   for(unsigned int i=0; i<targetdist_modifer_pntrs_.size(); i++) {
     254           6 :     applyTargetDistModiferToGrid(targetdist_modifer_pntrs_[i].get());
     255             :   }
     256             :   //
     257         802 :   if(bias_cutoff_active_) {
     258          24 :     updateBiasCutoffForTargetDistGrid();
     259             :   }
     260             :   //
     261         802 :   if(shift_targetdist_to_zero_ && !(bias_cutoff_active_)) {
     262           3 :     setMinimumOfTargetDistGridToZero();
     263             :   }
     264         802 :   if(force_normalization_ && !(bias_cutoff_active_) ) {
     265          87 :     normalizeTargetDistGrid();
     266             :   }
     267         802 :   if(check_normalization_ && !(bias_cutoff_active_)) {
     268         691 :     double normalization = integrateGrid(targetdist_grid_pntr_.get());
     269             :     const double normalization_thrshold = 0.1;
     270         691 :     if(normalization < 1.0-normalization_thrshold || normalization > 1.0+normalization_thrshold) {
     271             :       std::string norm_str;
     272           9 :       Tools::convert(normalization,norm_str);
     273           9 :       std::string msg = "the target distribution grid is not properly normalized, integrating over the grid gives: " + norm_str + " - You can avoid this problem by using the NORMALIZE keyword";
     274           9 :       warning(msg);
     275             :     }
     276             :   }
     277             :   //
     278         802 :   if(check_nonnegative_) {
     279             :     const double nonnegative_thrshold = -0.02;
     280         799 :     double grid_min_value = targetdist_grid_pntr_->getMinValue();
     281         799 :     if(grid_min_value<nonnegative_thrshold) {
     282             :       std::string grid_min_value_str;
     283           0 :       Tools::convert(grid_min_value,grid_min_value_str);
     284           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";
     285           0 :       warning(msg);
     286             :     }
     287             :   }
     288             :   //
     289         802 :   if(check_nan_inf_) {
     290           0 :     checkNanAndInf();
     291             :   }
     292             :   //
     293         802 : }
     294             : 
     295             : 
     296          24 : void TargetDistribution::updateBiasCutoffForTargetDistGrid() {
     297          24 :   plumed_massert(vesbias_pntr_!=NULL,"The VesBias has to be linked to use updateBiasCutoffForTargetDistGrid()");
     298          24 :   plumed_massert(vesbias_pntr_->biasCutoffActive(),"updateBiasCutoffForTargetDistGrid() should only be used if the bias cutoff is active");
     299          24 :   plumed_massert(getBiasWithoutCutoffGridPntr()!=NULL,"the bias without cutoff grid has to be linked");
     300             :   //
     301          48 :   std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(targetdist_grid_pntr_.get());
     302             :   double norm = 0.0;
     303        2624 :   for(Grid::index_t l=0; l<targetdist_grid_pntr_->getSize(); l++) {
     304        2600 :     double value = targetdist_grid_pntr_->getValue(l);
     305        2600 :     double bias = getBiasWithoutCutoffGridPntr()->getValue(l);
     306        2600 :     double deriv_factor_swf = 0.0;
     307        2600 :     double swf = vesbias_pntr_->getBiasCutoffSwitchingFunction(bias,deriv_factor_swf);
     308             :     // this comes from the p(s)
     309        2600 :     value *= swf;
     310        2600 :     norm += integration_weights[l]*value;
     311             :     // this comes from the derivative of V(s)
     312        2600 :     value *= deriv_factor_swf;
     313        2600 :     targetdist_grid_pntr_->setValue(l,value);
     314             :   }
     315          24 :   targetdist_grid_pntr_->scaleAllValuesAndDerivatives(1.0/norm);
     316          24 : }
     317             : 
     318           6 : void TargetDistribution::applyTargetDistModiferToGrid(TargetDistModifer* modifer_pntr) {
     319          12 :   std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(targetdist_grid_pntr_.get());
     320             :   double norm = 0.0;
     321       21212 :   for(Grid::index_t l=0; l<targetdist_grid_pntr_->getSize(); l++) {
     322       21206 :     double value = targetdist_grid_pntr_->getValue(l);
     323       21206 :     std::vector<double> cv_values = targetdist_grid_pntr_->getPoint(l);
     324       21206 :     value = modifer_pntr->getModifedTargetDistValue(value,cv_values);
     325       21206 :     norm += integration_weights[l]*value;
     326       21206 :     targetdist_grid_pntr_->setValue(l,value);
     327       21206 :     log_targetdist_grid_pntr_->setValue(l,-std::log(value));
     328             :   }
     329           6 :   targetdist_grid_pntr_->scaleAllValuesAndDerivatives(1.0/norm);
     330           6 :   log_targetdist_grid_pntr_->setMinToZero();
     331           6 : }
     332             : 
     333             : 
     334          29 : void TargetDistribution::updateLogTargetDistGrid() {
     335       44070 :   for(Grid::index_t l=0; l<targetdist_grid_pntr_->getSize(); l++) {
     336       44041 :     log_targetdist_grid_pntr_->setValue(l,-std::log(targetdist_grid_pntr_->getValue(l)));
     337             :   }
     338          29 :   log_targetdist_grid_pntr_->setMinToZero();
     339          29 : }
     340             : 
     341             : 
     342           3 : void TargetDistribution::setMinimumOfTargetDistGridToZero() {
     343           3 :   targetDistGrid().setMinToZero();
     344           3 :   normalizeTargetDistGrid();
     345           3 :   updateLogTargetDistGrid();
     346           3 : }
     347             : 
     348             : 
     349           8 : void TargetDistribution::readInRestartTargetDistGrid(const std::string& grid_fname) {
     350           8 :   plumed_massert(isDynamic(),"this should only be used for dynamically updated target distributions!");
     351           8 :   IFile gridfile;
     352           8 :   if(!gridfile.FileExist(grid_fname)) {
     353           0 :     plumed_merror(getName()+": problem with reading previous target distribution when restarting, cannot find file " + grid_fname);
     354             :   }
     355           8 :   gridfile.open(grid_fname);
     356          16 :   std::unique_ptr<GridBase> restart_grid = GridBase::create("targetdist",grid_args_,gridfile,false,false,false);
     357           8 :   if(restart_grid->getSize()!=targetdist_grid_pntr_->getSize()) {
     358           0 :     plumed_merror(getName()+": problem with reading previous target distribution when restarting, the grid is not of the correct size!");
     359             :   }
     360           8 :   VesTools::copyGridValues(restart_grid.get(),targetdist_grid_pntr_.get());
     361           8 :   updateLogTargetDistGrid();
     362           8 : }
     363             : 
     364           1 : void TargetDistribution::clearLogTargetDistGrid() {
     365           1 :   log_targetdist_grid_pntr_->clear();
     366           1 : }
     367             : 
     368             : 
     369           0 : void TargetDistribution::checkNanAndInf() {
     370           0 :   for(Grid::index_t l=0; l<targetdist_grid_pntr_->getSize(); l++) {
     371           0 :     double value = targetdist_grid_pntr_->getValue(l);
     372           0 :     if(std::isnan(value) || std::isinf(value)) {
     373             :       std::string vs;
     374           0 :       Tools::convert(value,vs);
     375           0 :       std::vector<double> p = targetdist_grid_pntr_->getPoint(l);
     376             :       std::string ps;
     377           0 :       Tools::convert(p[0],ps);
     378           0 :       ps = "(" + ps;
     379           0 :       for(unsigned int k=1; k<p.size(); k++) {
     380             :         std::string t1;
     381           0 :         Tools::convert(p[k],t1);
     382           0 :         ps = ps + "," + t1;
     383             :       }
     384           0 :       ps = ps + ")";
     385           0 :       plumed_merror(getName()+": problem with target distribution, the value at " + ps + " is " + vs);
     386             :     }
     387             :   }
     388           0 : }
     389             : 
     390             : }
     391             : }