LCOV - code coverage report
Current view: top level - function - FuncSumHills.cpp (source / functions) Hit Total Coverage
Test: plumed test coverage Lines: 278 318 87.4 %
Date: 2020-11-18 11:20:57 Functions: 18 21 85.7 %

          Line data    Source code
       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2012-2019 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed 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             :    plumed 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 plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "ActionRegister.h"
      23             : #include "Function.h"
      24             : #include "tools/Exception.h"
      25             : #include "tools/Communicator.h"
      26             : #include "tools/BiasRepresentation.h"
      27             : #include "tools/File.h"
      28             : #include "tools/Tools.h"
      29             : #include "tools/Stopwatch.h"
      30             : #include <iostream>
      31             : 
      32             : using namespace std;
      33             : 
      34             : namespace PLMD {
      35             : namespace function {
      36             : 
      37             : 
      38             : //+PLUMEDOC FUNCTION FUNCSUMHILLS
      39             : /*
      40             : This function is intended to be called by the command line tool sum_hills
      41             : and it is meant to integrate a HILLS file or an HILLS file interpreted as
      42             : a histogram i a variety of ways. Therefore it is not expected that you use this
      43             : during your dynamics (it will crash!)
      44             : 
      45             : In the future one could implement periodic integration during the metadynamics
      46             : or straightforward MD as a tool to check convergence
      47             : 
      48             : \par Examples
      49             : 
      50             : There are currently no examples for this keyword.
      51             : 
      52             : */
      53             : //+ENDPLUMEDOC
      54             : 
      55             : class FilesHandler {
      56             :   vector <string> filenames;
      57             :   vector <IFile*>  ifiles;
      58             :   Action *action;
      59             :   Log *log;
      60             :   bool parallelread;
      61             :   unsigned beingread;
      62             :   bool isopen;
      63             : public:
      64             :   FilesHandler(const vector<string> &filenames, const bool &parallelread,  Action &myaction, Log &mylog);
      65             :   bool readBunch(BiasRepresentation *br, int stride);
      66             :   bool scanOneHill(BiasRepresentation *br, IFile *ifile );
      67             :   void getMinMaxBin(vector<Value*> vals, Communicator &cc, vector<double> &vmin, vector<double> &vmax, vector<unsigned> &vbin);
      68             :   void getMinMaxBin(vector<Value*> vals, Communicator &cc, vector<double> &vmin, vector<double> &vmax, vector<unsigned> &vbin, vector<double> &histosigma);
      69             :   ~FilesHandler();
      70             : };
      71          22 : FilesHandler::FilesHandler(const vector<string> &filenames, const bool &parallelread, Action &action, Log &mylog ):filenames(filenames),log(&mylog),parallelread(parallelread),beingread(0),isopen(false) {
      72          11 :   this->action=&action;
      73          58 :   for(unsigned i=0; i<filenames.size(); i++) {
      74          12 :     IFile *ifile = new IFile();
      75          12 :     ifile->link(action);
      76          12 :     ifiles.push_back(ifile);
      77          24 :     plumed_massert((ifile->FileExist(filenames[i])), "the file "+filenames[i]+" does not exist " );
      78             :   }
      79             : 
      80          11 : }
      81             : 
      82          22 : FilesHandler::~FilesHandler() {
      83          58 :   for(unsigned i=0; i<ifiles.size(); i++) delete ifiles[i];
      84          11 : }
      85             : 
      86             : // note that the FileHandler is completely transparent respect to the biasrepresentation
      87             : // no check are made at this level
      88          12 : bool FilesHandler::readBunch(BiasRepresentation *br, int stride = -1) {
      89             :   bool morefiles; morefiles=true;
      90          12 :   if(parallelread) {
      91           0 :     (*log)<<"  doing parallelread \n";
      92           0 :     plumed_merror("parallelread is not yet implemented !!!");
      93             :   } else {
      94          12 :     (*log)<<"  doing serialread \n";
      95             :     // read one by one hills
      96             :     // is the type defined? if not, assume it is a gaussian
      97             :     IFile *ff;
      98          24 :     ff=ifiles[beingread];
      99          12 :     if(!isopen) {
     100          22 :       (*log)<<"  opening file "<<filenames[beingread]<<"\n";
     101          22 :       ff->open(filenames[beingread]); isopen=true;
     102             :     }
     103             :     int n;
     104             :     while(true) {
     105             :       bool fileisover=true;
     106        5566 :       while(scanOneHill(br,ff)) {
     107             :         // here do the dump if needed
     108        5554 :         n=br->getNumberOfKernels();
     109        5554 :         if(stride>0 && n%stride==0 && n!=0  ) {
     110           1 :           (*log)<<"  done with this chunk: now with "<<n<<" kernels  \n";
     111             :           fileisover=false;
     112             :           break;
     113             :         }
     114             :       }
     115             :       if(fileisover) {
     116          24 :         (*log)<<"  closing file "<<filenames[beingread]<<"\n";
     117          12 :         ff->close();
     118          12 :         isopen=false;
     119          12 :         (*log)<<"  now total "<<br->getNumberOfKernels()<<" kernels \n";
     120          12 :         beingread++;
     121          24 :         if(beingread<ifiles.size()) {
     122           2 :           ff=ifiles[beingread]; ff->open(filenames[beingread]);
     123           2 :           (*log)<<"  opening file "<<filenames[beingread]<<"\n";
     124           1 :           isopen=true;
     125             :         } else {
     126             :           morefiles=false;
     127          11 :           (*log)<<"  final chunk: now with "<<n<<" kernels  \n";
     128             :           break;
     129             :         }
     130             :       }
     131             :       // if there are no more files to read and this file is over then quit
     132             :       if(fileisover && !morefiles) {break;}
     133             :       // if you are in the middle of a file and you are here
     134             :       // then means that you read what you need to read
     135           2 :       if(!fileisover ) {break;}
     136             :     }
     137             :   }
     138          12 :   return morefiles;
     139             : }
     140           3 : void FilesHandler::getMinMaxBin(vector<Value*> vals, Communicator &cc, vector<double> &vmin, vector<double> &vmax, vector<unsigned> &vbin) {
     141             :   // create the representation (no grid)
     142           6 :   BiasRepresentation br(vals,cc);
     143             :   // read all the kernels
     144           3 :   readBunch(&br);
     145             :   // loop over the kernels and get the support
     146           3 :   br.getMinMaxBin(vmin,vmax,vbin);
     147           3 : }
     148           1 : void FilesHandler::getMinMaxBin(vector<Value*> vals, Communicator &cc, vector<double> &vmin, vector<double> &vmax, vector<unsigned> &vbin, vector<double> &histosigma) {
     149           2 :   BiasRepresentation br(vals,cc,histosigma);
     150             :   // read all the kernels
     151           1 :   readBunch(&br);
     152             :   // loop over the kernels and get the support
     153           1 :   br.getMinMaxBin(vmin,vmax,vbin);
     154             :   //for(unsigned i=0;i<vals.size();i++){cerr<<"XXX "<<vmin[i]<<" "<<vmax[i]<<" "<<vbin[i]<<"\n";}
     155           1 : }
     156        5566 : bool FilesHandler::scanOneHill(BiasRepresentation *br, IFile *ifile ) {
     157             :   double dummy;
     158       11132 :   if(ifile->scanField("time",dummy)) {
     159             :     //(*log)<<"   scanning one hill: "<<dummy<<" \n";
     160       16662 :     if(ifile->FieldExist("biasf")) ifile->scanField("biasf",dummy);
     161       11108 :     if(ifile->FieldExist("clock")) ifile->scanField("clock",dummy);
     162             :     // keep this intermediate function in case you need to parse more data in the future
     163        5554 :     br->pushKernel(ifile);
     164             :     //(*log)<<"  read hill\n";
     165        5554 :     if(br->hasSigmaInInput())ifile->allowIgnoredFields();
     166        5554 :     ifile->scanField();
     167             :     return true;
     168             :   } else {
     169             :     return false;
     170             :   }
     171             : }
     172             : 
     173             : 
     174         900 : double  mylog( double v1 ) {
     175         900 :   return log(v1);
     176             : }
     177             : 
     178        1800 : double  mylogder( double v1 ) {
     179        1800 :   return 1./v1;
     180             : }
     181             : 
     182             : 
     183             : 
     184             : class FuncSumHills :
     185             :   public Function
     186             : {
     187             :   vector<string> hillsFiles,histoFiles;
     188             :   vector<string> proj;
     189             :   int initstride;
     190             :   bool iscltool,integratehills,integratehisto,parallelread;
     191             :   bool negativebias;
     192             :   bool nohistory;
     193             :   bool minTOzero;
     194             :   bool doInt;
     195             :   double lowI_;
     196             :   double uppI_;
     197             :   double beta;
     198             :   string outhills,outhisto,fmt;
     199             :   BiasRepresentation *biasrep;
     200             :   BiasRepresentation *historep;
     201             : public:
     202             :   explicit FuncSumHills(const ActionOptions&);
     203             :   ~FuncSumHills();
     204             :   void calculate(); // this probably is not needed
     205             :   bool checkFilesAreExisting(const vector<string> & hills );
     206             :   static void registerKeywords(Keywords& keys);
     207             : };
     208             : 
     209        6459 : PLUMED_REGISTER_ACTION(FuncSumHills,"FUNCSUMHILLS")
     210             : 
     211           8 : void FuncSumHills::registerKeywords(Keywords& keys) {
     212           8 :   Function::registerKeywords(keys);
     213          16 :   keys.use("ARG");
     214          32 :   keys.add("optional","HILLSFILES"," source file for hills creation(may be the same as HILLS)"); // this can be a vector!
     215          32 :   keys.add("optional","HISTOFILES"," source file for histogram creation(may be the same as HILLS)"); // also this can be a vector!
     216          32 :   keys.add("optional","HISTOSIGMA"," sigmas for binning when the histogram correction is needed    ");
     217          32 :   keys.add("optional","PROJ"," only with sumhills: the projection on the cvs");
     218          32 :   keys.add("optional","KT"," only with sumhills: the kt factor when projection on cvs");
     219          32 :   keys.add("optional","GRID_MIN","the lower bounds for the grid");
     220          32 :   keys.add("optional","GRID_MAX","the upper bounds for the grid");
     221          32 :   keys.add("optional","GRID_BIN","the number of bins for the grid");
     222          32 :   keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
     223          32 :   keys.add("optional","INTERVAL","set monodimensional INTERVAL");
     224          32 :   keys.add("optional","OUTHILLS"," output file for hills ");
     225          32 :   keys.add("optional","OUTHISTO"," output file for histogram ");
     226          32 :   keys.add("optional","INITSTRIDE"," stride if you want an initial dump ");
     227          32 :   keys.add("optional","STRIDE"," stride when you do it on the fly ");
     228          24 :   keys.addFlag("ISCLTOOL",true,"use via plumed commandline: calculate at read phase and then go");
     229          24 :   keys.addFlag("PARALLELREAD",false,"read parallel HILLS file");
     230          24 :   keys.addFlag("NEGBIAS",false,"dump  negative bias ( -bias )   instead of the free energy: needed in welltempered with flexible hills ");
     231          24 :   keys.addFlag("NOHISTORY",false,"to be used with INITSTRIDE:  it splits the bias/histogram in pieces without previous history  ");
     232          24 :   keys.addFlag("MINTOZERO",false,"translate the resulting bias/histogram to have the minimum to zero  ");
     233          32 :   keys.add("optional","FMT","the format that should be used to output real numbers");
     234           8 : }
     235             : 
     236           7 : FuncSumHills::FuncSumHills(const ActionOptions&ao):
     237             :   Action(ao),
     238             :   Function(ao),
     239             :   initstride(-1),
     240             :   iscltool(false),
     241             :   integratehills(false),
     242             :   integratehisto(false),
     243             :   parallelread(false),
     244             :   negativebias(false),
     245             :   nohistory(false),
     246             :   minTOzero(false),
     247             :   doInt(false),
     248             :   lowI_(-1.),
     249             :   uppI_(-1.),
     250             :   beta(-1.),
     251             :   fmt("%14.9f"),
     252             :   biasrep(NULL),
     253          42 :   historep(NULL)
     254             : {
     255             : 
     256             :   // format
     257          14 :   parse("FMT",fmt);
     258           7 :   log<<"  Output format is "<<fmt<<"\n";
     259             :   // here read
     260             :   // Grid Stuff
     261           0 :   vector<std::string> gmin;
     262          14 :   parseVector("GRID_MIN",gmin);
     263           7 :   if(gmin.size()!=getNumberOfArguments() && gmin.size()!=0) error("not enough values for GRID_MIN");
     264           7 :   plumed_massert(gmin.size()==getNumberOfArguments() || gmin.size()==0,"need GRID_MIN argument for this") ;
     265           0 :   vector<std::string> gmax;
     266          21 :   parseVector("GRID_MAX",gmax);
     267           7 :   if(gmax.size()!=getNumberOfArguments() && gmax.size()!=0) error("not enough values for GRID_MAX");
     268           7 :   plumed_massert(gmax.size()==getNumberOfArguments() || gmax.size()==0,"need GRID_MAX argument for this") ;
     269             :   vector<unsigned> gbin;
     270             :   vector<double>   gspacing;
     271          21 :   parseVector("GRID_BIN",gbin);
     272           7 :   plumed_massert(gbin.size()==getNumberOfArguments() || gbin.size()==0,"need GRID_BIN argument for this") ;
     273           9 :   if(gbin.size()!=getNumberOfArguments() && gbin.size()!=0) error("not enough values for GRID_BIN");
     274          14 :   parseVector("GRID_SPACING",gspacing);
     275           7 :   plumed_massert(gspacing.size()==getNumberOfArguments() || gspacing.size()==0,"need GRID_SPACING argument for this") ;
     276          14 :   if(gspacing.size()!=getNumberOfArguments() && gspacing.size()!=0) error("not enough values for GRID_SPACING");
     277           7 :   if(gspacing.size()!=0 && gbin.size()==0) {
     278           0 :     log<<"  The number of bins will be estimated from GRID_SPACING\n";
     279           7 :   } else if(gspacing.size()!=0 && gbin.size()!=0) {
     280           0 :     log<<"  You specified both GRID_BIN and GRID_SPACING\n";
     281           0 :     log<<"  The more conservative (highest) number of bins will be used for each variable\n";
     282             :   }
     283           7 :   if(gspacing.size()!=0) for(unsigned i=0; i<getNumberOfArguments(); i++) {
     284           0 :       if(gbin.size()==0) gbin.assign(getNumberOfArguments(),1);
     285             :       double a,b;
     286           0 :       Tools::convert(gmin[i],a);
     287           0 :       Tools::convert(gmax[i],b);
     288           0 :       unsigned n=((b-a)/gspacing[i])+1;
     289           0 :       if(gbin[i]<n) gbin[i]=n;
     290             :     }
     291             : 
     292             :   // Inteval keyword
     293           7 :   vector<double> tmpI(2);
     294          14 :   parseVector("INTERVAL",tmpI);
     295           7 :   if(tmpI.size()!=2&&tmpI.size()!=0) error("both a lower and an upper limits must be provided with INTERVAL");
     296           7 :   else if(tmpI.size()==2) {
     297           0 :     lowI_=tmpI.at(0);
     298           0 :     uppI_=tmpI.at(1);
     299           0 :     if(getNumberOfArguments()!=1) error("INTERVAL limits correction works only for monodimensional metadynamics!");
     300           0 :     if(uppI_<lowI_) error("The Upper limit must be greater than the Lower limit!");
     301           0 :     doInt=true;
     302             :   }
     303           7 :   if(doInt) {
     304           0 :     log << "  Upper and Lower limits boundaries for the bias are activated at " << lowI_ << " - " << uppI_<<"\n";
     305           0 :     log << "  Using the same values as boundaries for the grid if not other value was defined (default: 200 bins)\n";
     306           0 :     std::ostringstream strsmin, strsmax;
     307           0 :     strsmin << lowI_;
     308           0 :     strsmax << uppI_;
     309           0 :     if(gmin.size()==0) gmin.push_back(strsmin.str());
     310           0 :     if(gmax.size()==0) gmax.push_back(strsmax.str());
     311           0 :     if(gbin.size()==0) gbin.push_back(200);
     312             :   }
     313             : 
     314             : 
     315             :   // hills file:
     316          14 :   parseVector("HILLSFILES",hillsFiles);
     317           7 :   if(hillsFiles.size()==0) {
     318           1 :     integratehills=false; // default behaviour
     319             :   } else {
     320           6 :     integratehills=true;
     321          33 :     for(unsigned i=0; i<hillsFiles.size(); i++) log<<"  hillsfile  : "<<hillsFiles[i]<<"\n";
     322             :   }
     323             :   // histo file:
     324          14 :   parseVector("HISTOFILES",histoFiles);
     325           7 :   if(histoFiles.size()==0) {
     326           6 :     integratehisto=false;
     327             :   } else {
     328           1 :     integratehisto=true;
     329           5 :     for(unsigned i=0; i<histoFiles.size(); i++) log<<"  histofile  : "<<histoFiles[i]<<"\n";
     330             :   }
     331             :   vector<double> histoSigma;
     332           7 :   if(integratehisto) {
     333           2 :     parseVector("HISTOSIGMA",histoSigma);
     334           8 :     for(unsigned i=0; i<histoSigma.size(); i++) log<<"  histosigma  : "<<histoSigma[i]<<"\n";
     335             :   }
     336             : 
     337             :   // needs a projection?
     338           7 :   proj.clear();
     339          14 :   parseVector("PROJ",proj);
     340           7 :   if(integratehills) {
     341           6 :     plumed_massert(proj.size()<getNumberOfArguments()," The number of projection must be less than the full list of arguments ");
     342             :   }
     343           7 :   if(integratehisto) {
     344           1 :     plumed_massert(proj.size()<=getNumberOfArguments()," The number of projection must be less or equal to the full list of arguments ");
     345             :   }
     346           8 :   if(integratehisto&&proj.size()==0) {
     347           5 :     for(unsigned i=0; i<getNumberOfArguments(); i++) proj.push_back(getPntrToArgument(i)->getName());
     348             :   }
     349             : 
     350             :   // add some automatic hills width: not in case stride is defined
     351             :   // since when you start from zero the automatic size will be zero!
     352          10 :   if(gmin.size()==0 || gmax.size()==0) {
     353           4 :     log<<"   \n";
     354           4 :     log<<"  No boundaries defined: need to do a prescreening of hills \n";
     355             :     std::vector<Value*> tmphillsvalues, tmphistovalues;
     356           4 :     if(integratehills) {
     357          21 :       for(unsigned i=0; i<getNumberOfArguments(); i++)tmphillsvalues.push_back( getPntrToArgument(i) );
     358             :     }
     359           4 :     if(integratehisto) {
     360           5 :       for(unsigned i=0; i<getNumberOfArguments(); i++) {
     361             :         std::string ss = getPntrToArgument(i)->getName();
     362          16 :         for(unsigned j=0; j<proj.size(); j++) {
     363           8 :           if(proj[j]==ss) tmphistovalues.push_back( getPntrToArgument(i) );
     364             :         }
     365             :       }
     366             :     }
     367             : 
     368           4 :     if(integratehills) {
     369           6 :       FilesHandler hillsHandler(hillsFiles,parallelread,*this, log);
     370             :       vector<double> vmin,vmax;
     371             :       vector<unsigned> vbin;
     372           6 :       hillsHandler.getMinMaxBin(tmphillsvalues,comm,vmin,vmax,vbin);
     373           3 :       log<<"  found boundaries from hillsfile: \n";
     374           3 :       gmin.resize(vmin.size());
     375           3 :       gmax.resize(vmax.size());
     376           3 :       if(gbin.size()==0) {
     377           2 :         gbin=vbin;
     378             :       } else {
     379           1 :         log<<"  found nbins in input, this overrides the automatic choice \n";
     380             :       }
     381          15 :       for(unsigned i=0; i<getNumberOfArguments(); i++) {
     382          12 :         Tools::convert(vmin[i],gmin[i]);
     383           6 :         Tools::convert(vmax[i],gmax[i]);
     384          30 :         log<<"  variable "<< getPntrToArgument(i)->getName()<<" min: "<<gmin[i]<<" max: "<<gmax[i]<<" nbin: "<<gbin[i]<<"\n";
     385             :       }
     386             :     }
     387             :     // if at this stage bins are not there then do it with histo
     388           4 :     if(gmin.size()==0) {
     389           2 :       FilesHandler histoHandler(histoFiles,parallelread,*this, log);
     390             :       vector<double> vmin,vmax;
     391             :       vector<unsigned> vbin;
     392           2 :       histoHandler.getMinMaxBin(tmphistovalues,comm,vmin,vmax,vbin,histoSigma);
     393           1 :       log<<"  found boundaries from histofile: \n";
     394           1 :       gmin.resize(vmin.size());
     395           1 :       gmax.resize(vmax.size());
     396           1 :       if(gbin.size()==0) {
     397           0 :         gbin=vbin;
     398             :       } else {
     399           1 :         log<<"  found nbins in input, this overrides the automatic choice \n";
     400             :       }
     401           8 :       for(unsigned i=0; i<proj.size(); i++) {
     402           2 :         Tools::convert(vmin[i],gmin[i]);
     403           2 :         Tools::convert(vmax[i],gmax[i]);
     404          10 :         log<<"  variable "<< proj[i] <<" min: "<<gmin[i]<<" max: "<<gmax[i]<<" nbin: "<<gbin[i]<<"\n";
     405             :       }
     406             :     }
     407           4 :     log<<"  done!\n";
     408           4 :     log<<"   \n";
     409             :   }
     410             : 
     411             : 
     412           7 :   if( proj.size() != 0 || integratehisto==true  ) {
     413           4 :     parse("KT",beta);
     414          10 :     for(unsigned i=0; i<proj.size(); i++) log<<"  projection "<<i<<" : "<<proj[i]<<"\n";
     415             :     // this should be only for projection or free energy from histograms
     416           2 :     plumed_massert(beta>0.,"if you make a projection or a histogram correction then you need KT flag!");
     417           2 :     beta=1./beta;
     418           2 :     log<<"  beta is "<<beta<<"\n";
     419             :   }
     420             :   // is a cltool: then you start and then die
     421          14 :   parseFlag("ISCLTOOL",iscltool);
     422             :   //
     423          14 :   parseFlag("NEGBIAS",negativebias);
     424             :   //
     425          14 :   parseFlag("PARALLELREAD",parallelread);
     426             :   // stride
     427          14 :   parse("INITSTRIDE",initstride);
     428             :   // output suffix or names
     429           7 :   if(initstride<0) {
     430           6 :     log<<"  Doing only one integration: no stride \n";
     431             :     outhills="fes.dat"; outhisto="histo.dat";
     432             :   }
     433             :   else {
     434             :     outhills="fes_"; outhisto="histo_";
     435           1 :     log<<"  Doing integration slices every "<<initstride<<" kernels\n";
     436           2 :     parseFlag("NOHISTORY",nohistory);
     437           1 :     if(nohistory)log<<"  nohistory: each stride block has no memory of the previous block\n";
     438             :   }
     439          14 :   parseFlag("MINTOZERO",minTOzero);
     440           7 :   if(minTOzero)log<<"  mintozero: bias/histogram will be translated to have the minimum value equal to zero\n";
     441             :   //what might it be this?
     442             :   // here start
     443             :   // want something right now?? do it and return
     444             :   // your argument is a set of cvs
     445             :   // then you need: a hills / a colvar-like file (to do a histogram)
     446             :   // create a bias representation for this
     447           7 :   if(iscltool) {
     448             : 
     449             :     std::vector<Value*> tmphillsvalues, tmphistovalues;
     450           7 :     if(integratehills) {
     451          30 :       for(unsigned i=0; i<getNumberOfArguments(); i++) {
     452             :         // allocate a new value from the old one: no deriv here
     453             :         // if we are summing hills then all the arguments are needed
     454          24 :         tmphillsvalues.push_back( getPntrToArgument(i) );
     455             :       }
     456             :     }
     457           7 :     if(integratehisto) {
     458           5 :       for(unsigned i=0; i<getNumberOfArguments(); i++) {
     459             :         std::string ss = getPntrToArgument(i)->getName();
     460          16 :         for(unsigned j=0; j<proj.size(); j++) {
     461           8 :           if(proj[j]==ss) tmphistovalues.push_back( getPntrToArgument(i) );
     462             :         }
     463             :       }
     464             :     }
     465             : 
     466             :     // check if the files exists
     467           7 :     if(integratehills) {
     468           6 :       checkFilesAreExisting(hillsFiles);
     469           6 :       biasrep=new BiasRepresentation(tmphillsvalues,comm, gmin, gmax, gbin, doInt, lowI_, uppI_);
     470           6 :       if(negativebias) {
     471           1 :         biasrep->setRescaledToBias(true);
     472           1 :         log<<"  required the -bias instead of the free energy \n";
     473           1 :         if(initstride<0) {outhills="negativebias.dat";}
     474             :         else {outhills="negativebias_";}
     475             :       }
     476             :     }
     477             : 
     478          14 :     parse("OUTHILLS",outhills);
     479          14 :     parse("OUTHISTO",outhisto);
     480           7 :     if(integratehills)log<<"  output file for fes/bias  is :  "<<outhills<<"\n";
     481           7 :     if(integratehisto)log<<"  output file for histogram is :  "<<outhisto<<"\n";
     482           7 :     checkRead();
     483             : 
     484           7 :     log<<"\n";
     485           7 :     log<<"  Now calculating...\n";
     486           7 :     log<<"\n";
     487             : 
     488             :     // here it defines the column to be histogrammed, tmpvalues should be only
     489             :     // the list of the collective variable one want to consider
     490           7 :     if(integratehisto) {
     491           1 :       checkFilesAreExisting(histoFiles);
     492           1 :       historep=new BiasRepresentation(tmphistovalues,comm,gmin,gmax,gbin,histoSigma);
     493             :     }
     494             : 
     495             :     // decide how to source hills ( serial/parallel )
     496             :     // here below the input control
     497             :     // say how many hills and it will read them from the
     498             :     // bunch of files provided, will update the representation
     499             :     // of hills (i.e. a list of hills and the associated grid)
     500             : 
     501             :     // decide how to source colvars ( serial parallel )
     502             :     FilesHandler *hillsHandler;
     503             :     FilesHandler *histoHandler;
     504             : 
     505             :     hillsHandler=NULL;
     506             :     histoHandler=NULL;
     507             : 
     508           7 :     if(integratehills)  hillsHandler=new FilesHandler(hillsFiles,parallelread,*this, log);
     509           7 :     if(integratehisto)  histoHandler=new FilesHandler(histoFiles,parallelread,*this, log);
     510             : 
     511             :     Stopwatch sw;
     512             : 
     513          14 :     sw.start("0 Summing hills");
     514             : 
     515             :     // read a number of hills and put in the bias representation
     516             :     int nfiles=0;
     517           7 :     bool ibias=integratehills; bool ihisto=integratehisto;
     518             :     while(true) {
     519           8 :       if(  integratehills  && ibias  ) {
     520           7 :         if(nohistory) {biasrep->clear(); log<<"  clearing history before reading a new block\n";};
     521           7 :         log<<"  reading hills: \n";
     522           7 :         ibias=hillsHandler->readBunch(biasrep,initstride) ; log<<"\n";
     523             :       }
     524             : 
     525           8 :       if(  integratehisto  && ihisto ) {
     526           1 :         if(nohistory) {historep->clear(); log<<"  clearing history before reading a new block\n";};
     527           1 :         log<<"  reading histogram: \n";
     528           1 :         ihisto=histoHandler->readBunch(historep,initstride) ;  log<<"\n";
     529             :       }
     530             : 
     531             :       // dump: need to project?
     532           8 :       if(proj.size()!=0) {
     533             : 
     534           3 :         if(integratehills) {
     535             : 
     536           2 :           log<<"  Bias: Projecting on subgrid... \n";
     537           2 :           BiasWeight Bw(beta);
     538           4 :           Grid biasGrid=*(biasrep->getGridPtr());
     539           4 :           Grid smallGrid=biasGrid.project(proj,&Bw);
     540           4 :           OFile gridfile; gridfile.link(*this);
     541           4 :           std::ostringstream ostr; ostr<<nfiles;
     542             :           string myout;
     543          10 :           if(initstride>0) { myout=outhills+ostr.str()+".dat" ;} else {myout=outhills;}
     544           2 :           log<<"  Bias: Writing subgrid on file "<<myout<<" \n";
     545           2 :           gridfile.open(myout);
     546           2 :           if(minTOzero) smallGrid.setMinToZero();
     547             :           smallGrid.setOutputFmt(fmt);
     548           2 :           smallGrid.writeToFile(gridfile);
     549           2 :           gridfile.close();
     550           2 :           if(!ibias)integratehills=false;// once you get to the final bunch just give up
     551             :         }
     552             :         // this should be removed
     553           3 :         if(integratehisto) {
     554             : 
     555           1 :           log<<"  Histo: Projecting on subgrid... \n";
     556           2 :           Grid histoGrid=*(historep->getGridPtr());
     557             : 
     558           2 :           OFile gridfile; gridfile.link(*this);
     559           2 :           std::ostringstream ostr; ostr<<nfiles;
     560             :           string myout;
     561           1 :           if(initstride>0) { myout=outhisto+ostr.str()+".dat" ;} else {myout=outhisto;}
     562           1 :           log<<"  Histo: Writing subgrid on file "<<myout<<" \n";
     563           1 :           gridfile.open(myout);
     564             : 
     565           1 :           histoGrid.applyFunctionAllValuesAndDerivatives(&mylog,&mylogder);
     566           1 :           histoGrid.scaleAllValuesAndDerivatives(-1./beta);
     567           1 :           if(minTOzero) histoGrid.setMinToZero();
     568             :           histoGrid.setOutputFmt(fmt);
     569           1 :           histoGrid.writeToFile(gridfile);
     570           1 :           gridfile.close();
     571             : 
     572           1 :           if(!ihisto)integratehisto=false;// once you get to the final bunch just give up
     573             :         }
     574             : 
     575             :       } else {
     576             : 
     577           5 :         if(integratehills) {
     578             : 
     579          10 :           Grid biasGrid=*(biasrep->getGridPtr());
     580           5 :           biasGrid.scaleAllValuesAndDerivatives(-1.);
     581             : 
     582          10 :           OFile gridfile; gridfile.link(*this);
     583          10 :           std::ostringstream ostr; ostr<<nfiles;
     584             :           string myout;
     585           5 :           if(initstride>0) { myout=outhills+ostr.str()+".dat" ;} else {myout=outhills;}
     586           5 :           log<<"  Writing full grid on file "<<myout<<" \n";
     587           5 :           gridfile.open(myout);
     588             : 
     589           5 :           if(minTOzero) biasGrid.setMinToZero();
     590             :           biasGrid.setOutputFmt(fmt);
     591           5 :           biasGrid.writeToFile(gridfile);
     592           5 :           gridfile.close();
     593             :           // rescale back prior to accumulate
     594           5 :           if(!ibias)integratehills=false;// once you get to the final bunch just give up
     595             :         }
     596           5 :         if(integratehisto) {
     597             : 
     598           0 :           Grid histoGrid=*(historep->getGridPtr());
     599             :           // do this if you want a free energy from a grid, otherwise do not
     600           0 :           histoGrid.applyFunctionAllValuesAndDerivatives(&mylog,&mylogder);
     601           0 :           histoGrid.scaleAllValuesAndDerivatives(-1./beta);
     602             : 
     603           0 :           OFile gridfile; gridfile.link(*this);
     604           0 :           std::ostringstream ostr; ostr<<nfiles;
     605             :           string myout;
     606           0 :           if(initstride>0) { myout=outhisto+ostr.str()+".dat" ;} else {myout=outhisto;}
     607           0 :           log<<"  Writing full grid on file "<<myout<<" \n";
     608           0 :           gridfile.open(myout);
     609             : 
     610             :           // also this is usefull only for free energy
     611           0 :           if(minTOzero) histoGrid.setMinToZero();
     612             :           histoGrid.setOutputFmt(fmt);
     613           0 :           histoGrid.writeToFile(gridfile);
     614           0 :           gridfile.close();
     615             : 
     616           0 :           if(!ihisto)integratehisto=false; // once you get to the final bunch just give up
     617             :         }
     618             :       }
     619           8 :       if ( !ibias && !ihisto) break; //when both are over then just quit
     620             : 
     621           1 :       nfiles++;
     622           1 :     }
     623             : 
     624           7 :     if( hillsHandler ) delete hillsHandler;
     625           7 :     if( histoHandler ) delete histoHandler;
     626             : 
     627          14 :     sw.stop("0 Summing hills");
     628             : 
     629           7 :     log<<sw;
     630             : 
     631             :     return;
     632             :   }
     633             :   // just an initialization but you need to do something on the fly?: need to connect with a metad run and its grid representation
     634             :   // your argument is a metad run
     635             :   // if the grid does not exist crash and say that you need some data
     636             :   // otherwise just link with it
     637             : 
     638             : }
     639             : 
     640           0 : void FuncSumHills::calculate() {
     641             :   // this should be connected only with a grid representation to metadynamics
     642             :   // at regular time just dump it
     643           0 :   plumed_merror("You should have never got here: this stuff is not yet implemented!");
     644             : }
     645             : 
     646          21 : FuncSumHills::~FuncSumHills() {
     647           7 :   if(historep) delete historep;
     648           7 :   if(biasrep) delete biasrep;
     649          14 : }
     650             : 
     651           7 : bool FuncSumHills::checkFilesAreExisting(const vector<string> & hills ) {
     652           7 :   plumed_massert(hills.size()!=0,"the number of  files provided should be at least one" );
     653           7 :   IFile *ifile = new IFile();
     654           7 :   ifile->link(*this);
     655          38 :   for(unsigned i=0; i< hills.size(); i++) {
     656          16 :     plumed_massert(ifile->FileExist(hills[i]),"missing file "+hills[i]);
     657             :   }
     658           7 :   delete ifile;
     659           7 :   return true;
     660             : 
     661             : }
     662             : 
     663             : }
     664             : 
     665        4839 : }
     666             : 
     667             : 

Generated by: LCOV version 1.13