Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :     Copyright (c) 2017 of Haochuan Chen (excluding colvar_UIestimator.h)
       3             :     Copyright (c) 2017 of Haohao Fu (colvar_UIestimator.h)
       4             : 
       5             :     This program is free software: you can redistribute it and/or modify
       6             :     it under the terms of the GNU Lesser General Public License as published
       7             :     by the Free Software Foundation, either version 3 of the License, or
       8             :     (at your option) any later version.
       9             : 
      10             :     This program is distributed in the hope that it will be useful,
      11             :     but WITHOUT ANY WARRANTY; without even the implied warranty of
      12             :     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      13             :     GNU Lesser General Public License for more details.
      14             : 
      15             :     You should have received a copy of the GNU Lesser General Public License
      16             :     along with this program.  If not, see <http://www.gnu.org/licenses/>.
      17             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      18             : #ifdef __PLUMED_HAS_BOOST_SERIALIZATION
      19             : #include "cltools/CLTool.h"
      20             : #include "cltools/CLToolRegister.h"
      21             : #include "config/Config.h"
      22             : #include "core/ActionRegister.h"
      23             : #include "DRR.h"
      24             : #include "tools/Tools.h"
      25             : #include "tools/Units.h"
      26             : #include <boost/archive/binary_iarchive.hpp>
      27             : #include <boost/archive/binary_oarchive.hpp>
      28             : #include <boost/serialization/vector.hpp>
      29             : #include <cstdlib>
      30             : #include <fstream>
      31             : #include <iostream>
      32             : #include <string>
      33             : 
      34             : using namespace PLMD;
      35             : using namespace cltools;
      36             : 
      37             : namespace PLMD {
      38             : namespace drr {
      39             : 
      40             : //+PLUMEDOC EABFMOD_TOOLS drr_tool
      41             : /*
      42             :  - Extract .grad and .count files from the binary output .drrstate
      43             :  - Merge windows
      44             : 
      45             : \par Examples
      46             : 
      47             : The following command will extract .grad and .count files.
      48             : \verbatim
      49             : plumed drr_tool --extract eabf.drrstate
      50             : \endverbatim
      51             : 
      52             : The following command will merge windows of two .drrstate file, and output the
      53             : .grad and .count files.
      54             : \verbatim
      55             : plumed drr_tool --merge win1.drrstate,win2.drrstate
      56             : \endverbatim
      57             : 
      58             : After getting the .grad and .count file, you can do numerical integration by
      59             : using abf_integrate tool from
      60             : https://github.com/Colvars/colvars/tree/master/colvartools
      61             : \verbatim
      62             : abf_integrate eabf.czar.grad
      63             : \endverbatim
      64             : \note
      65             : The abf_integrate in colvartools is in kcal/mol, so it may be better to use --units kcal/mol when running drr_tool
      66             : 
      67             : */
      68             : //+ENDPLUMEDOC
      69             : 
      70             : using std::vector;
      71             : using std::string;
      72             : 
      73          12 : class drrtool : public CLTool {
      74             : public:
      75             :   static void registerKeywords(Keywords &keys);
      76             :   explicit drrtool(const CLToolOptions &co);
      77             :   int main(FILE *in, FILE *out, Communicator &pc);
      78             :   void extractdrr(const vector<string> &filename);
      79             :   void mergewindows(const vector<string> &filename);
      80             :   void calcDivergence(const vector<string> &filename);
      81           0 :   string description() const { return "Extract or merge the drrstate files."; }
      82             : 
      83             : private:
      84             :   bool verbosity;
      85             :   Units units;
      86             :   const string suffix{".drrstate"};
      87             : };
      88             : 
      89        7364 : PLUMED_REGISTER_CLTOOL(drrtool, "drr_tool")
      90             : 
      91        1839 : void drrtool::registerKeywords(Keywords &keys) {
      92        1839 :   CLTool::registerKeywords(keys);
      93        7356 :   keys.add("optional", "--extract", "Extract drrstate file(s)");
      94        7356 :   keys.add("optional", "--merge", "Merge eABF windows");
      95        7356 :   keys.add("optional", "--divergence", "Calculate divergence of gradient field (experimental)");
      96        9195 :   keys.add("compulsory","--units","kj/mol","the units of energy can be kj/mol, kcal/mol, j/mol, eV or the conversion factor from kj/mol");
      97        5517 :   keys.addFlag("-v", false, "Verbose output");
      98        1839 : }
      99             : 
     100           4 : drrtool::drrtool(const CLToolOptions &co) : CLTool(co) {
     101           4 :   inputdata = commandline;
     102           4 :   verbosity = false;
     103           4 : }
     104             : 
     105           4 : int drrtool::main(FILE *in, FILE *out, Communicator &pc) {
     106           8 :   parseFlag("-v", verbosity);
     107           4 :   vector<string> stateFilesToExtract;
     108             :   string unitname;
     109           8 :   parse("--units",unitname);
     110           4 :   units.setEnergy( unitname );
     111           8 :   bool doextract = parseVector("--extract", stateFilesToExtract);
     112           4 :   if (doextract) {
     113           2 :     extractdrr(stateFilesToExtract);
     114             :   }
     115           4 :   vector<string> stateFilesToMerge;
     116           8 :   bool domerge = parseVector("--merge", stateFilesToMerge);
     117           4 :   if (domerge) {
     118           1 :     mergewindows(stateFilesToMerge);
     119             :   }
     120           4 :   vector<string> stateFilesToDivergence;
     121           8 :   bool dodivergence = parseVector("--divergence", stateFilesToDivergence);
     122           4 :   if (dodivergence) {
     123           1 :     calcDivergence(stateFilesToDivergence);
     124             :   }
     125           4 :   return 0;
     126             : }
     127             : 
     128           2 : void drrtool::extractdrr(const vector<string> &filename) {
     129           6 :   #pragma omp parallel for
     130           4 :   for (size_t j = 0; j < filename.size(); ++j) {
     131           4 :     std::ifstream in;
     132           4 :     in.open(filename[j]);
     133             :     boost::archive::binary_iarchive ia(in);
     134             :     long long int step;
     135             :     vector<double> fict;
     136             :     vector<double> vfict;
     137             :     vector<double> vfict_laststep;
     138             :     vector<double> ffict;
     139             :     ABF abfgrid;
     140             :     CZAR czarestimator;
     141             :     ia >> step >> fict >> vfict >> vfict_laststep >> ffict >> abfgrid >>
     142             :        czarestimator;
     143           2 :     in.close();
     144           2 :     abfgrid.setOutputUnit(units.getEnergy());
     145             :     czarestimator.setOutputUnit(units.getEnergy());
     146           2 :     if (verbosity) {
     147           1 :       std::cout << "Output units factor: " << units.getEnergy() << '\n';
     148             :       std::cout << "Dumping information of extended variables..." << '\n';
     149           1 :       std::cout << "Step: " << step << '\n';
     150           3 :       for (size_t i = 0; i < fict.size(); ++i) {
     151           1 :         std::cout << "Dimension[" << i + 1 << "]:\n"
     152           1 :                   << "  Coordinate: " << fict[i] << '\n'
     153           1 :                   << "  Velocity: " << vfict[i] << '\n'
     154           1 :                   << "  Velocity(laststep): " << vfict_laststep[i] << '\n'
     155           1 :                   << "  Force: " << ffict[i] << '\n';
     156             :       }
     157             :       std::cout << "Dumping counts and gradients from grids..." << '\n';
     158             :     }
     159           2 :     string outputname(filename[j]);
     160           4 :     outputname = outputname.substr(0, outputname.length() - suffix.length());
     161           2 :     if (verbosity)
     162             :       std::cout << "Writing ABF(naive) estimator files..." << '\n';
     163           2 :     abfgrid.writeAll(outputname);
     164           2 :     if (verbosity)
     165             :       std::cout << "Writing CZAR estimator files..." << '\n';
     166           2 :     czarestimator.writeAll(outputname);
     167             :   }
     168           2 : }
     169             : 
     170           1 : void drrtool::mergewindows(const vector<string> &filename) {
     171           1 :   if (filename.size() < 2) {
     172             :     std::cerr << "ERROR! You need at least two .drrstate file to merge windows!" << std::endl;
     173           0 :     std::abort();
     174             :   }
     175             :   // Read grid into abfs and czars;
     176           1 :   vector<ABF> abfs;
     177           1 :   vector<CZAR> czars;
     178           3 :   for (auto it_fn = filename.begin(); it_fn != filename.end(); ++it_fn) {
     179           4 :     std::ifstream in;
     180           2 :     in.open((*it_fn));
     181             :     boost::archive::binary_iarchive ia(in);
     182             :     long long int step;
     183             :     vector<double> fict;
     184             :     vector<double> vfict;
     185             :     vector<double> vfict_laststep;
     186             :     vector<double> ffict;
     187             :     ABF abfgrid;
     188             :     CZAR czarestimator;
     189             :     ia >> step >> fict >> vfict >> vfict_laststep >> ffict >> abfgrid >>
     190             :        czarestimator;
     191           2 :     abfgrid.setOutputUnit(units.getEnergy());
     192             :     czarestimator.setOutputUnit(units.getEnergy());
     193           2 :     abfs.push_back(abfgrid);
     194           2 :     czars.push_back(czarestimator);
     195           2 :     in.close();
     196             :   }
     197           1 :   CZAR cmerged = CZAR::mergewindow(czars[0], czars[1]);
     198           1 :   ABF amerged = ABF::mergewindow(abfs[0], abfs[1]);
     199           1 :   for (size_t i = 2; i < czars.size(); ++i) {
     200           0 :     cmerged = CZAR::mergewindow(cmerged, czars[i]);
     201           0 :     amerged = ABF::mergewindow(amerged, abfs[i]);
     202             :   }
     203             :   // Generate new file name for merged grad and count
     204           2 :   vector<string> tmp_name = filename;
     205           6 :   std::transform(std::begin(tmp_name), std::end(tmp_name), std::begin(tmp_name), [&](string s) {return s.substr(0, s.find(suffix));});
     206           7 :   string mergename = std::accumulate(std::begin(tmp_name), std::end(tmp_name), string(""), [](string a, string b) {return a + b + "+";});
     207           2 :   mergename = mergename.substr(0, mergename.size() - 1);
     208           1 :   cmerged.writeAll(mergename);
     209           1 :   amerged.writeAll(mergename);
     210           1 : }
     211             : 
     212           1 : void drrtool::calcDivergence(const vector<string> &filename) {
     213           3 :   #pragma omp parallel for
     214           2 :   for (size_t j = 0; j < filename.size(); ++j) {
     215           2 :     std::ifstream in;
     216           2 :     in.open(filename[j]);
     217             :     boost::archive::binary_iarchive ia(in);
     218             :     long long int step;
     219             :     vector<double> fict;
     220             :     vector<double> vfict;
     221             :     vector<double> vfict_laststep;
     222             :     vector<double> ffict;
     223             :     ABF abfgrid;
     224             :     CZAR czarestimator;
     225             :     ia >> step >> fict >> vfict >> vfict_laststep >> ffict >> abfgrid >>
     226             :        czarestimator;
     227           1 :     in.close();
     228           1 :     abfgrid.setOutputUnit(units.getEnergy());
     229             :     czarestimator.setOutputUnit(units.getEnergy());
     230           1 :     if (verbosity) {
     231           0 :       std::cout << "Output units factor: " << units.getEnergy() << '\n';
     232             :       std::cout << "Dumping information of extended variables..." << '\n';
     233           0 :       std::cout << "Step: " << step << '\n';
     234           0 :       for (size_t i = 0; i < fict.size(); ++i) {
     235           0 :         std::cout << "Dimension[" << i + 1 << "]:\n"
     236           0 :                   << "  Coordinate: " << fict[i] << '\n'
     237           0 :                   << "  Velocity: " << vfict[i] << '\n'
     238           0 :                   << "  Velocity(laststep): " << vfict_laststep[i] << '\n'
     239           0 :                   << "  Force: " << ffict[i] << '\n';
     240             :       }
     241             :       std::cout << "Dumping counts and gradients from grids..." << '\n';
     242             :     }
     243           1 :     string outputname(filename[j]);
     244           2 :     outputname = outputname.substr(0, outputname.length() - suffix.length());
     245           1 :     abfgrid.writeDivergence(outputname);
     246           1 :     czarestimator.writeDivergence(outputname);
     247             :   }
     248           1 : }
     249             : 
     250             : } // End of namespace
     251        5517 : }
     252             : 
     253             : #endif