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 "GridIntegrationWeights.h"
      24             : 
      25             : #include "tools/Grid.h"
      26             : #include "tools/File.h"
      27             : #include "tools/Exception.h"
      28             : 
      29             : 
      30             : namespace PLMD {
      31             : namespace ves {
      32             : 
      33        2036 : std::vector<double> GridIntegrationWeights::getIntegrationWeights(const Grid* grid_pntr, const std::string& fname_weights_grid, const std::string& weights_type) {
      34        2036 :   std::vector<double> dx = grid_pntr->getDx();
      35        2036 :   std::vector<bool> isPeriodic = grid_pntr->getIsPeriodic();
      36        2036 :   std::vector<unsigned int> nbins = grid_pntr->getNbin();
      37             :   std::vector<std::vector<double> > weights_perdim;
      38        4580 :   for(unsigned int k=0; k<grid_pntr->getDimension(); k++) {
      39             :     std::vector<double> weights_tmp;
      40        2544 :     if(weights_type=="trapezoidal") {
      41        2544 :       weights_tmp = getOneDimensionalTrapezoidalWeights(nbins[k],dx[k],isPeriodic[k]);
      42             :     }
      43             :     else {
      44           0 :       plumed_merror("getIntegrationWeights: unknown weight type, the available type is trapezoidal");
      45             :     }
      46        2544 :     weights_perdim.push_back(weights_tmp);
      47             :   }
      48             : 
      49        2036 :   std::vector<double> weights_vector(grid_pntr->getSize(),0.0);
      50     5323391 :   for(Grid::index_t l=0; l<grid_pntr->getSize(); l++) {
      51     5321355 :     std::vector<unsigned int> ind = grid_pntr->getIndices(l);
      52             :     double value = 1.0;
      53    15690443 :     for(unsigned int k=0; k<grid_pntr->getDimension(); k++) {
      54    10369088 :       value *= weights_perdim[k][ind[k]];
      55             :     }
      56     5321355 :     weights_vector[l] = value;
      57             :   }
      58             : 
      59        2036 :   if(fname_weights_grid.size()>0) {
      60           0 :     Grid weights_grid = Grid(*grid_pntr);
      61           0 :     for(Grid::index_t l=0; l<weights_grid.getSize(); l++) {
      62           0 :       weights_grid.setValue(l,weights_vector[l]);
      63             :     }
      64           0 :     OFile ofile;
      65           0 :     ofile.enforceBackup();
      66           0 :     ofile.open(fname_weights_grid);
      67           0 :     weights_grid.writeToFile(ofile);
      68           0 :     ofile.close();
      69           0 :   }
      70             :   //
      71        2036 :   return weights_vector;
      72        2036 : }
      73             : 
      74             : 
      75          19 : void GridIntegrationWeights::getOneDimensionalIntegrationPointsAndWeights(std::vector<double>& points, std::vector<double>& weights, const unsigned int nbins, const double min, const double max, const std::string& weights_type) {
      76          19 :   double dx = (max-min)/(static_cast<double>(nbins)-1.0);
      77          19 :   points.resize(nbins);
      78       19038 :   for(unsigned int i=0; i<nbins; i++) {points[i] = min + i*dx;}
      79          19 :   if(weights_type=="trapezoidal") {
      80          19 :     weights = getOneDimensionalTrapezoidalWeights(nbins,dx,false);
      81             :   }
      82             :   else {
      83           0 :     plumed_merror("getOneDimensionalIntegrationWeights: unknown weight type, the available type is trapezoidal");
      84             :   }
      85          19 : }
      86             : 
      87             : 
      88        2563 : std::vector<double> GridIntegrationWeights::getOneDimensionalTrapezoidalWeights(const unsigned int nbins, const double dx, const bool periodic) {
      89        2563 :   std::vector<double> weights_1d(nbins);
      90      444294 :   for(unsigned int i=1; i<(nbins-1); i++) {
      91      441731 :     weights_1d[i] = dx;
      92             :   }
      93        2563 :   if(!periodic) {
      94        1202 :     weights_1d[0]= 0.5*dx;
      95        1202 :     weights_1d[(nbins-1)]= 0.5*dx;
      96             :   }
      97             :   else {
      98             :     // as for periodic arguments the first point should be counted twice as the
      99             :     // grid doesn't include its periodic copy
     100        1361 :     weights_1d[0]= dx;
     101        1361 :     weights_1d[(nbins-1)]= dx;
     102             :   }
     103        2563 :   return weights_1d;
     104             : }
     105             : 
     106             : 
     107             : }
     108             : }