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 "BasisFunctions.h"
      24             : #include "TargetDistribution.h"
      25             : #include "VesBias.h"
      26             : #include "VesTools.h"
      27             : #include "GridIntegrationWeights.h"
      28             : 
      29             : #include "tools/Grid.h"
      30             : #include "tools/Tools.h"
      31             : 
      32             : 
      33             : namespace PLMD {
      34             : namespace ves {
      35             : 
      36         216 : void BasisFunctions::registerKeywords(Keywords& keys) {
      37         216 :   Action::registerKeywords(keys);
      38         864 :   keys.add("compulsory","ORDER","The order of the basis function expansion.");
      39         864 :   keys.add("compulsory","MINIMUM","The minimum of the interval on which the basis functions are defined.");
      40         864 :   keys.add("compulsory","MAXIMUM","The maximum of the interval on which the basis functions are defined.");
      41         864 :   keys.add("hidden","NGRID_POINTS","The number of grid points used for numerical integrals");
      42         648 :   keys.addFlag("DEBUG_INFO",false,"Print out more detailed information about the basis set. Useful for debugging.");
      43         648 :   keys.addFlag("NUMERICAL_INTEGRALS",false,"Calculate basis function integral for the uniform distribution numerically. Useful for debugging.");
      44         216 : }
      45             : 
      46             : 
      47         203 : BasisFunctions::BasisFunctions(const ActionOptions&ao):
      48             :   Action(ao),
      49             :   print_debug_info_(false),
      50             :   has_been_set(false),
      51             :   description_("Undefined"),
      52             :   type_("Undefined"),
      53             :   norder_(0),
      54             :   nbasis_(1),
      55             :   bf_label_prefix_("f"),
      56             :   bf_labels_(nbasis_,"f0"),
      57             :   periodic_(false),
      58             :   orthogonal_(false),
      59             :   interval_bounded_(true),
      60             :   interval_intrinsic_min_str_("1.0"),
      61             :   interval_intrinsic_max_str_("-1.0"),
      62             :   interval_intrinsic_min_(1.0),
      63             :   interval_intrinsic_max_(-1.0),
      64             :   interval_intrinsic_range_(0.0),
      65             :   interval_intrinsic_mean_(0.0),
      66             :   interval_min_str_(""),
      67             :   interval_max_str_(""),
      68             :   interval_min_(0.0),
      69             :   interval_max_(0.0),
      70             :   interval_range_(0.0),
      71             :   interval_mean_(0.0),
      72             :   argT_derivf_(1.0),
      73             :   numerical_uniform_integrals_(false),
      74             :   nbins_(1001),
      75             :   uniform_integrals_(nbasis_,0.0),
      76             :   vesbias_pntr_(NULL),
      77         609 :   action_pntr_(NULL)
      78             : {
      79         203 :   bf_keywords_.push_back(getName());
      80         406 :   if(keywords.exists("ORDER")) {
      81         573 :     parse("ORDER",norder_); addKeywordToList("ORDER",norder_);
      82             :   }
      83         203 :   nbasis_=norder_+1;
      84             :   //
      85             :   std::string str_imin; std::string str_imax;
      86         607 :   if(keywords.exists("MINIMUM") && keywords.exists("MAXIMUM")) {
      87         804 :     parse("MINIMUM",str_imin); addKeywordToList("MINIMUM",str_imin);
      88         804 :     parse("MAXIMUM",str_imax); addKeywordToList("MAXIMUM",str_imax);
      89             :   }
      90             :   else {
      91             :     str_imin = "-1.0";
      92             :     str_imax = "1.0";
      93             :   }
      94             :   interval_min_str_ = str_imin;
      95             :   interval_max_str_ = str_imax;
      96         203 :   if(!Tools::convert(str_imin,interval_min_)) {
      97           0 :     plumed_merror(getName()+": cannot convert the value given in MINIMUM to a double");
      98             :   }
      99         203 :   if(!Tools::convert(str_imax,interval_max_)) {
     100           0 :     plumed_merror(getName()+": cannot convert the value given in MAXIMUM to a double");
     101             :   }
     102         203 :   if(interval_min_>interval_max_) {plumed_merror(getName()+": MINIMUM and MAXIMUM are not correctly defined");}
     103             :   //
     104         406 :   parseFlag("DEBUG_INFO",print_debug_info_);
     105         406 :   if(keywords.exists("NUMERICAL_INTEGRALS")) {
     106         346 :     parseFlag("NUMERICAL_INTEGRALS",numerical_uniform_integrals_);
     107             :   }
     108         406 :   if(keywords.exists("NGRID_POINTS")) {
     109         402 :     parse("NGRID_POINTS",nbins_);
     110             :   }
     111             :   // log.printf(" %s \n",getKeywordString().c_str());
     112             : 
     113         203 : }
     114             : 
     115             : 
     116          32 : void BasisFunctions::setIntrinsicInterval(const double interval_intrinsic_min_in, const double interval_intrinsic_max_in) {
     117          32 :   interval_intrinsic_min_ = interval_intrinsic_min_in;
     118          32 :   interval_intrinsic_max_ = interval_intrinsic_max_in;
     119          32 :   VesTools::convertDbl2Str(interval_intrinsic_min_,interval_intrinsic_min_str_);
     120          32 :   VesTools::convertDbl2Str(interval_intrinsic_max_,interval_intrinsic_max_str_);
     121          32 :   plumed_massert(interval_intrinsic_min_<interval_intrinsic_max_,"setIntrinsicInterval: intrinsic intervals are not defined correctly");
     122          32 : }
     123             : 
     124             : 
     125         171 : void BasisFunctions::setIntrinsicInterval(const std::string& interval_intrinsic_min_str_in, const std::string& interval_intrinsic_max_str_in) {
     126         171 :   interval_intrinsic_min_str_ = interval_intrinsic_min_str_in;
     127         171 :   interval_intrinsic_max_str_ = interval_intrinsic_max_str_in;
     128         171 :   if(!Tools::convert(interval_intrinsic_min_str_,interval_intrinsic_min_)) {
     129           0 :     plumed_merror("setIntrinsicInterval: cannot convert string value given for the minimum of the intrinsic interval to a double");
     130             :   }
     131         171 :   if(!Tools::convert(interval_intrinsic_max_str_,interval_intrinsic_max_)) {
     132           0 :     plumed_merror("setIntrinsicInterval: cannot convert string value given for the maximum of the intrinsic interval to a double");
     133             :   }
     134         171 :   plumed_massert(interval_intrinsic_min_<interval_intrinsic_max_,"setIntrinsicInterval: intrinsic intervals are not defined correctly");
     135         171 : }
     136             : 
     137             : 
     138           0 : void BasisFunctions::setInterval(const double interval_min_in, const double interval_max_in) {
     139           0 :   interval_min_ = interval_min_in;
     140           0 :   interval_max_ = interval_max_in;
     141           0 :   VesTools::convertDbl2Str(interval_min_,interval_min_str_);
     142           0 :   VesTools::convertDbl2Str(interval_max_,interval_max_str_);
     143           0 :   plumed_massert(interval_min_<interval_max_,"setInterval: intervals are not defined correctly");
     144           0 : }
     145             : 
     146             : 
     147           2 : void BasisFunctions::setInterval(const std::string& interval_min_str_in, const std::string& interval_max_str_in) {
     148           2 :   interval_min_str_ = interval_min_str_in;
     149           2 :   interval_max_str_ = interval_max_str_in;
     150           2 :   if(!Tools::convert(interval_min_str_,interval_min_)) {
     151           0 :     plumed_merror("setInterval: cannot convert string value given for the minimum of the interval to a double");
     152             :   }
     153           2 :   if(!Tools::convert(interval_max_str_,interval_max_)) {
     154           0 :     plumed_merror("setInterval: cannot convert string value given for the maximum of the interval to a double");
     155             :   }
     156           2 :   plumed_massert(interval_min_<interval_max_,"setInterval: intervals are not defined correctly");
     157           2 : }
     158             : 
     159             : 
     160         203 : void BasisFunctions::setupInterval() {
     161             :   // if(!intervalBounded()){plumed_merror("setupInterval() only works for bounded interval");}
     162         203 :   interval_intrinsic_range_ = interval_intrinsic_max_-interval_intrinsic_min_;
     163         203 :   interval_intrinsic_mean_  = 0.5*(interval_intrinsic_max_+interval_intrinsic_min_);
     164         203 :   interval_range_ = interval_max_-interval_min_;
     165         203 :   interval_mean_  = 0.5*(interval_max_+interval_min_);
     166         203 :   argT_derivf_ = interval_intrinsic_range_/interval_range_;
     167         203 : }
     168             : 
     169             : 
     170          81 : void BasisFunctions::setupLabels() {
     171         950 :   for(unsigned int i=0; i < nbasis_; i++) {
     172         869 :     std::string is; Tools::convert(i,is);
     173        3476 :     bf_labels_[i]=bf_label_prefix_+is+"(s)";
     174             :   }
     175          81 : }
     176             : 
     177             : 
     178          32 : void BasisFunctions::setupUniformIntegrals() {
     179          32 :   numerical_uniform_integrals_=true;
     180          32 :   numericalUniformIntegrals();
     181          32 : }
     182             : 
     183             : 
     184         203 : void BasisFunctions::setupBF() {
     185         203 :   if(interval_intrinsic_min_>interval_intrinsic_max_) {plumed_merror("setupBF: default intervals are not correctly set");}
     186         203 :   setupInterval();
     187         203 :   setupLabels();
     188         203 :   if(bf_labels_.size()==1) {plumed_merror("setupBF: the labels of the basis functions are not correct.");}
     189         203 :   if(!numerical_uniform_integrals_) {setupUniformIntegrals();}
     190           6 :   else {numericalUniformIntegrals();}
     191         203 :   if(uniform_integrals_.size()==1) {plumed_merror("setupBF: the integrals of the basis functions is not correct.");}
     192         406 :   if(type_=="Undefined") {plumed_merror("setupBF: the type of the basis function is not defined.");}
     193         406 :   if(description_=="Undefined") {plumed_merror("setupBF: the description of the basis function is not defined.");}
     194         203 :   has_been_set=true;
     195         203 :   printInfo();
     196         203 : }
     197             : 
     198             : 
     199         203 : void BasisFunctions::printInfo() const {
     200         203 :   if(!has_been_set) {plumed_merror("the basis set has not be setup correctly");}
     201         203 :   log.printf("  One-dimensional basis set\n");
     202         203 :   log.printf("   Description: %s\n",description_.c_str());
     203         203 :   log.printf("   Type: %s\n",type_.c_str());
     204         203 :   if(periodic_) {log.printf("   The basis functions are periodic\n");}
     205         203 :   if(orthogonal_) {
     206         354 :     if(getInnerProductWeightStr()=="1") {
     207         173 :       log.printf("   The basis functions are orthogonal\n");
     208             :     }
     209             :     else {
     210          12 :       log.printf("   The basis functions are orthogonal with respect to the weight %s\n",getInnerProductWeightStr().c_str());
     211             :     }
     212             :   }
     213             :   else {
     214          26 :     log.printf("   The basis functions are not orthogonal\n");
     215             :   }
     216         203 :   log.printf("   Order of basis set: %u\n",norder_);
     217         203 :   log.printf("   Number of basis functions: %u\n",nbasis_);
     218             :   // log.printf("   Interval of basis set: %f to %f\n",interval_min_,interval_max_);
     219         203 :   log.printf("   Interval of basis set: %s to %s\n",interval_min_str_.c_str(),interval_max_str_.c_str());
     220         203 :   log.printf("   Description of basis functions:\n");
     221        2412 :   for(unsigned int i=0; i < nbasis_; i++) {log.printf("    %2u       %10s\n",i,bf_labels_[i].c_str());}
     222             :   //
     223         203 :   if(print_debug_info_) {
     224          44 :     log.printf("  Debug information:\n");
     225             :     // log.printf("   Default interval of basis set: [%f,%f]\n",interval_intrinsic_min_,interval_intrinsic_max_);
     226          44 :     log.printf("   Intrinsic interval of basis set: [%s,%s]\n",interval_intrinsic_min_str_.c_str(),interval_intrinsic_max_str_.c_str());
     227          44 :     log.printf("   Intrinsic interval of basis set: range=%f,  mean=%f\n",interval_intrinsic_range_,interval_intrinsic_mean_);
     228             :     // log.printf("   Defined interval of basis set: [%f,%f]\n",interval_min_,interval_max_);
     229          44 :     log.printf("   Defined interval of basis set: [%s,%s]\n",interval_min_str_.c_str(),interval_max_str_.c_str());
     230          44 :     log.printf("   Defined interval of basis set: range=%f,  mean=%f\n",interval_range_,interval_mean_);
     231          44 :     log.printf("   Derivative factor due to interval translation: %f\n",argT_derivf_);
     232          44 :     log.printf("   Integral of basis functions over the interval:\n");
     233          44 :     if(numerical_uniform_integrals_) {log.printf("   Note: calculated numerically\n");}
     234        1288 :     for(unsigned int i=0; i < nbasis_; i++) {log.printf("    %2u       %16.10f\n",i,uniform_integrals_[i]);}
     235          44 :     log.printf("   --------------------------\n");
     236             :   }
     237         203 : }
     238             : 
     239             : 
     240           0 : void BasisFunctions::linkVesBias(VesBias* vesbias_pntr_in) {
     241           0 :   vesbias_pntr_ = vesbias_pntr_in;
     242           0 :   action_pntr_ = static_cast<Action*>(vesbias_pntr_in);
     243           0 : }
     244             : 
     245             : 
     246           0 : void BasisFunctions::linkAction(Action* action_pntr_in) {
     247           0 :   action_pntr_ = action_pntr_in;
     248           0 : }
     249             : 
     250             : 
     251          38 : void BasisFunctions::numericalUniformIntegrals() {
     252          76 :   std::vector<std::string> grid_min(1); grid_min[0]=intervalMinStr();
     253         114 :   std::vector<std::string> grid_max(1); grid_max[0]=intervalMaxStr();
     254          76 :   std::vector<unsigned int> grid_bins(1); grid_bins[0]=nbins_;
     255         114 :   std::vector<Value*> arguments(1); arguments[0]= new Value(NULL,"arg",false);
     256          50 :   if(arePeriodic()) {arguments[0]->setDomain(intervalMinStr(),intervalMaxStr());}
     257          34 :   else {arguments[0]->setNotPeriodic();}
     258          76 :   Grid* uniform_grid = new Grid("uniform",arguments,grid_min,grid_max,grid_bins,false,false);
     259             :   //
     260          38 :   double inverse_normalization = 1.0/(intervalMax()-intervalMin());
     261       38110 :   for(Grid::index_t l=0; l<uniform_grid->getSize(); l++) {
     262       38072 :     uniform_grid->setValue(l,inverse_normalization);
     263             :   }
     264          76 :   uniform_integrals_ = numericalTargetDistributionIntegralsFromGrid(uniform_grid);
     265          38 :   delete arguments[0]; arguments.clear();
     266          76 :   delete uniform_grid;
     267          38 : }
     268             : 
     269             : 
     270         188 : std::vector<double> BasisFunctions::numericalTargetDistributionIntegralsFromGrid(const Grid* grid_pntr) const {
     271         188 :   plumed_massert(grid_pntr!=NULL,"the grid is not defined");
     272         188 :   plumed_massert(grid_pntr->getDimension()==1,"the target distribution grid should be one dimensional");
     273             :   //
     274         188 :   std::vector<double> targetdist_integrals(nbasis_,0.0);
     275         564 :   std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(grid_pntr);
     276             : 
     277       83369 :   for(Grid::index_t k=0; k < grid_pntr->getSize(); k++) {
     278      249543 :     double arg = grid_pntr->getPoint(k)[0];
     279       83181 :     std::vector<double> bf_values(nbasis_);
     280       83181 :     std::vector<double> bf_derivs(nbasis_);
     281       83181 :     bool inside=true;
     282       83181 :     double argT=0.0;
     283       83181 :     getAllValues(arg,argT,inside,bf_values,bf_derivs);
     284     1145995 :     for(unsigned int i=0; i < nbasis_; i++) {
     285     3437985 :       targetdist_integrals[i] += (integration_weights[k] * grid_pntr->getValue(k)) * bf_values[i];
     286             :     }
     287             :   }
     288             :   // assume that the first function is the constant
     289         188 :   bool inside=true;
     290         188 :   double argT=0.0;
     291         188 :   targetdist_integrals[0] = getValue(0.0,0,argT,inside);
     292         188 :   return targetdist_integrals;
     293             : }
     294             : 
     295             : 
     296         233 : std::vector<double> BasisFunctions::getTargetDistributionIntegrals(const TargetDistribution* targetdist_pntr) const {
     297         233 :   if(targetdist_pntr==NULL) {
     298             :     return getUniformIntegrals();
     299             :   }
     300             :   else {
     301             :     Grid* targetdist_grid = targetdist_pntr->getTargetDistGridPntr();
     302         150 :     return numericalTargetDistributionIntegralsFromGrid(targetdist_grid);
     303             :   }
     304             : }
     305             : 
     306             : 
     307         166 : std::string BasisFunctions::getKeywordString() const {
     308         166 :   std::string str_keywords=bf_keywords_[0];
     309        2420 :   for(unsigned int i=1; i<bf_keywords_.size(); i++) {str_keywords+=" "+bf_keywords_[i];}
     310         166 :   return str_keywords;
     311             : }
     312             : 
     313             : 
     314         188 : double BasisFunctions::getValue(const double arg, const unsigned int n, double& argT, bool& inside_range) const {
     315         188 :   plumed_massert(n<numberOfBasisFunctions(),"getValue: n is outside range of the defined order of the basis set");
     316         188 :   inside_range=true;
     317         188 :   std::vector<double> tmp_values(numberOfBasisFunctions());
     318         188 :   std::vector<double> tmp_derivs(numberOfBasisFunctions());
     319         188 :   getAllValues(arg, argT, inside_range, tmp_values, tmp_derivs);
     320         564 :   return tmp_values[n];
     321             : }
     322             : 
     323             : 
     324        5716 : void BasisFunctions::getAllValuesNumericalDerivs(const double arg, double& argT, bool& inside_range, std::vector<double>& values, std::vector<double>& derivs) const {
     325             :   // use forward difference, unless very close to the boundary
     326             :   double delta = sqrt(epsilon);
     327       11432 :   if((arg+delta)>intervalMax()) {
     328             :     delta *= -1.0;
     329             :   }
     330        5716 :   inside_range=true;
     331        5716 :   std::vector<double> values_delta(numberOfBasisFunctions());
     332        5716 :   std::vector<double> derivs_dummy(numberOfBasisFunctions());
     333        5716 :   getAllValues(arg+delta, argT, inside_range, values_delta, derivs_dummy);
     334        5716 :   getAllValues(arg, argT, inside_range, values, derivs_dummy);
     335      200678 :   for(unsigned int i=0; i<numberOfBasisFunctions(); i++) {
     336      292443 :     derivs[i] = (values_delta[i]-values[i])/delta;
     337             :   }
     338        5716 : }
     339             : 
     340             : 
     341          83 : void BasisFunctions::getMultipleValue(const std::vector<double>& args, std::vector<double>& argsT, std::vector<std::vector<double> >& values, std::vector<std::vector<double> >& derivs, const bool numerical_deriv) const {
     342          83 :   argsT.resize(args.size());
     343             :   values.clear();
     344             :   derivs.clear();
     345       53096 :   for(unsigned int i=0; i<args.size(); i++) {
     346       26465 :     std::vector<double> tmp_values(getNumberOfBasisFunctions());
     347       26465 :     std::vector<double> tmp_derivs(getNumberOfBasisFunctions());
     348       26465 :     bool inside_interval=true;
     349       26465 :     if(!numerical_deriv) {
     350       41498 :       getAllValues(args[i],argsT[i],inside_interval,tmp_values,tmp_derivs);
     351             :     } else {
     352        5716 :       getAllValuesNumericalDerivs(args[i],argsT[i],inside_interval,tmp_values,tmp_derivs);
     353             :     }
     354       26465 :     values.push_back(tmp_values);
     355       26465 :     derivs.push_back(tmp_derivs);
     356             :   }
     357          83 : }
     358             : 
     359             : 
     360          83 : void BasisFunctions::writeBasisFunctionsToFile(OFile& ofile_values, OFile& ofile_derivs, const std::string& min_in, const std::string& max_in, unsigned int nbins_in, const bool ignore_periodicity, const std::string& output_fmt, const bool numerical_deriv) const {
     361          83 :   std::vector<std::string> min(1); min[0]=min_in;
     362         166 :   std::vector<std::string> max(1); max[0]=max_in;
     363         166 :   std::vector<unsigned int> nbins(1); nbins[0]=nbins_in;
     364          83 :   std::vector<Value*> value_pntr(1);
     365         249 :   value_pntr[0]= new Value(NULL,"arg",false);
     366         137 :   if(arePeriodic() && !ignore_periodicity) {value_pntr[0]->setDomain(intervalMinStr(),intervalMaxStr());}
     367          65 :   else {value_pntr[0]->setNotPeriodic();}
     368         249 :   Grid args_grid = Grid("grid",value_pntr,min,max,nbins,false,false);
     369             : 
     370          83 :   std::vector<double> args(args_grid.getSize(),0.0);
     371       53096 :   for(unsigned int i=0; i<args.size(); i++) {
     372       79395 :     args[i] = args_grid.getPoint(i)[0];
     373             :   }
     374             :   std::vector<double> argsT;
     375          83 :   std::vector<std::vector<double> > values;
     376          83 :   std::vector<std::vector<double> > derivs;
     377             : 
     378         581 :   ofile_values.addConstantField("bf_keywords").printField("bf_keywords","{"+getKeywordString()+"}");
     379         581 :   ofile_derivs.addConstantField("bf_keywords").printField("bf_keywords","{"+getKeywordString()+"}");
     380             : 
     381         332 :   ofile_values.addConstantField("min").printField("min",intervalMinStr());
     382         332 :   ofile_values.addConstantField("max").printField("max",intervalMaxStr());
     383             : 
     384         332 :   ofile_derivs.addConstantField("min").printField("min",intervalMinStr());
     385         332 :   ofile_derivs.addConstantField("max").printField("max",intervalMaxStr());
     386             : 
     387         415 :   ofile_values.addConstantField("nbins").printField("nbins",static_cast<int>(args_grid.getNbin()[0]));
     388         415 :   ofile_derivs.addConstantField("nbins").printField("nbins",static_cast<int>(args_grid.getNbin()[0]));
     389             : 
     390          83 :   if(arePeriodic()) {
     391          84 :     ofile_values.addConstantField("periodic").printField("periodic","true");
     392          84 :     ofile_derivs.addConstantField("periodic").printField("periodic","true");
     393             :   }
     394             :   else {
     395         248 :     ofile_values.addConstantField("periodic").printField("periodic","false");
     396         248 :     ofile_derivs.addConstantField("periodic").printField("periodic","false");
     397             :   }
     398             : 
     399          83 :   getMultipleValue(args,argsT,values,derivs,numerical_deriv);
     400          83 :   ofile_values.fmtField(output_fmt);
     401          83 :   ofile_derivs.fmtField(output_fmt);
     402       53013 :   for(unsigned int i=0; i<args.size(); i++) {
     403       52930 :     ofile_values.printField("arg",args[i]);
     404       52930 :     ofile_derivs.printField("arg",args[i]);
     405      901292 :     for(unsigned int k=0; k<getNumberOfBasisFunctions(); k++) {
     406     1696724 :       ofile_values.printField(getBasisFunctionLabel(k),values[i][k]);
     407     1696724 :       ofile_derivs.printField("d_"+getBasisFunctionLabel(k),derivs[i][k]);
     408             :     }
     409       26465 :     ofile_values.printField();
     410       26465 :     ofile_derivs.printField();
     411             :   }
     412          83 :   ofile_values.fmtField();
     413          83 :   ofile_derivs.fmtField();
     414             : 
     415         166 :   delete value_pntr[0]; value_pntr.clear();
     416             : 
     417          83 : }
     418             : 
     419             : 
     420          83 : std::vector<std::vector<double> > BasisFunctions::getAllInnerProducts() const {
     421             :   // setup temp grid
     422         166 :   std::vector<std::string> grid_min(1); grid_min[0]=intervalMinStr();
     423         249 :   std::vector<std::string> grid_max(1); grid_max[0]=intervalMaxStr();
     424         166 :   std::vector<unsigned int> grid_bins(1); grid_bins[0]=nbins_;
     425         249 :   std::vector<Value*> arguments(1); arguments[0]= new Value(NULL,"arg",false);
     426         146 :   if(arePeriodic()) {arguments[0]->setDomain(intervalMinStr(),intervalMaxStr());}
     427          62 :   else {arguments[0]->setNotPeriodic();}
     428         166 :   Grid* grid = new Grid("uniform",arguments,grid_min,grid_max,grid_bins,false,false);
     429       83228 :   for(Grid::index_t l=0; l<grid->getSize(); l++) {
     430       83145 :     grid->setValue(l,1.0);
     431             :   }
     432             :   //
     433          83 :   std::vector<std::vector<double> > inner_products = getAllInnerProducts(grid);
     434          83 :   delete grid;
     435          83 :   return inner_products;
     436             :   //
     437             : }
     438             : 
     439             : 
     440          83 : std::vector<std::vector<double> > BasisFunctions::getAllInnerProducts(const Grid* grid_pntr) const {
     441             : 
     442         249 :   std::vector<std::vector<double> > inner_products(numberOfBasisFunctions(), std::vector<double>(numberOfBasisFunctions()));
     443        2828 :   for(unsigned int i=0; i<numberOfBasisFunctions(); i++) {
     444       28323 :     for(unsigned int j=i; j<numberOfBasisFunctions(); j++) {
     445       26992 :       inner_products[i][j] = inner_products[j][i] = getInnerProduct(i,j,grid_pntr);
     446             :     }
     447             :   }
     448          83 :   return inner_products;
     449             : }
     450             : 
     451             : 
     452       13496 : double BasisFunctions::getInnerProduct(const unsigned int n, const unsigned int m, const Grid* grid_pntr) const {
     453       13496 :   plumed_massert(grid_pntr->getDimension()==1,"the grid must be one-dimensional");
     454       40488 :   std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(grid_pntr);
     455             :   double sum = 0.0;
     456    13534023 :   for(Grid::index_t k=0; k < grid_pntr->getSize(); k++) {
     457    40561581 :     double arg = grid_pntr->getPoint(k)[0];
     458             :     double argT;
     459    13520527 :     bool inside_range=true;
     460    13520527 :     std::vector<double> values(numberOfBasisFunctions());
     461    13520527 :     std::vector<double> derivs(numberOfBasisFunctions());
     462    13520527 :     getAllValues(arg, argT, inside_range, values, derivs);
     463    13520527 :     plumed_massert(inside_range,"the basis functions values must be inside the range of the defined interval!");
     464             :     //
     465    40561581 :     sum += integration_weights[k]*values[n]*values[m];
     466             :     // disable weight for the moment
     467             :     // sum += integration_weights[k]*values[n]*values[m]*getInnerProductWeight(arg);
     468             :     //
     469             :   }
     470       13496 :   return sum;
     471             : }
     472             : 
     473             : 
     474          83 : void BasisFunctions::writeInnerProductsToFiles(OFile& ofile, const std::string& output_fmt) const {
     475             : 
     476          83 :   std::string int_fmt = "%8d";
     477          83 :   ofile.fmtField(output_fmt);
     478             : 
     479          83 :   std::string weight = getInnerProductWeightStr();
     480             : 
     481         166 :   std::vector<std::vector<double> > inner_products = getAllInnerProducts();
     482          83 :   char* s1 = new char[20];
     483             : 
     484             :   // diagonal part
     485          83 :   ofile.printf("#! Diagonal Part\n");
     486         249 :   ofile.addConstantField("weight").printField("weight",weight);
     487        2828 :   for(unsigned int i=0; i<numberOfBasisFunctions(); i++) {
     488        3993 :     sprintf(s1,int_fmt.c_str(),i); ofile.printField("i",s1);
     489        3993 :     sprintf(s1,int_fmt.c_str(),i); ofile.printField("j",s1);
     490        3993 :     ofile.printField("inner_product",inner_products[i][i]);
     491        1331 :     ofile.printField();
     492             :   }
     493             : 
     494          83 :   ofile.clearFields();
     495             :   // off-diagonal part
     496          83 :   ofile.printf("#! Off-Diagonal Part\n");
     497         249 :   ofile.addConstantField("weight").printField("weight",weight);
     498        1497 :   for(unsigned int i=0; i<numberOfBasisFunctions(); i++) {
     499       26992 :     for(unsigned int j=i+1; j<numberOfBasisFunctions(); j++) {
     500       36495 :       sprintf(s1,int_fmt.c_str(),i); ofile.printField("i",s1);
     501       36495 :       sprintf(s1,int_fmt.c_str(),j); ofile.printField("j",s1);
     502       36495 :       ofile.printField("inner_product",inner_products[i][j]);
     503       12165 :       ofile.printField();
     504             :     }
     505             :   }
     506          83 : }
     507             : 
     508             : 
     509             : 
     510             : }
     511             : }