diff --git a/CHANGES/v2.4.md b/CHANGES/v2.4.md
index 8bb4723f726e8c4fcbe5e82a01d6d138bb22ff73..c890624f22051191e05fe68b8a62943840a3e165 100644
--- a/CHANGES/v2.4.md
+++ b/CHANGES/v2.4.md
@@ -35,10 +35,10 @@ Changes from version 2.3 which are relevant for users:
   - A new DRR module have been included, contributed by Haochuan Chen and Haohao Fu.
     This module implements the following methods:
     - \ref DRR
-    - \ref drr_tools
+    - \ref drr_tool
 - New collective variables:
   - \ref DIMER (thanks to Marco Nava).
-  - \ref EFFSolv : EEF1 implicit solvent solvation energy
+  - \ref EFFSOLV : EEF1 implicit solvent solvation energy
   - \ref ADAPTIVE_PATH : Adaptive path variables using the method from \cite BerndAdaptivePath
 - New actions:
   - \ref INENVELOPE
diff --git a/src/bias/MaxEnt.cpp b/src/bias/MaxEnt.cpp
index c1786defa7309f7b37165418ef3fdbb4359654e0..d8ef1e29a1a1a65d52a8437ff8439356325724be 100644
--- a/src/bias/MaxEnt.cpp
+++ b/src/bias/MaxEnt.cpp
@@ -48,7 +48,7 @@ namespace bias {
 
 //+PLUMEDOC BIAS MAXENT
 /*
-Add a linear biasing potential on one or more variables f_{i}\left(\boldsymbol{x}\right) satisfying the maximum entropy principle as proposed in Ref. \cite cesari2016maxent .
+Add a linear biasing potential on one or more variables \f$f_{i}\left(\boldsymbol{x}\right)\f$ satisfying the maximum entropy principle as proposed in Ref. \cite cesari2016maxent .
 
 \warning
     Notice that syntax is still under revision and might change
diff --git a/src/drr/DRR.cpp b/src/drr/DRR.cpp
index 324c53b2b52d664c6f207f976fe86d5cd3efbd8b..76be92568884af1e889c5190d16928edc2b1e351 100644
--- a/src/drr/DRR.cpp
+++ b/src/drr/DRR.cpp
@@ -102,7 +102,7 @@ void DRRForceGrid::fillTable(const std::vector<std::vector<double>> &in) {
 
 DRRForceGrid::DRRForceGrid()
   : suffix(""), ndims(0), dimensions(0), sampleSize(0), forceSize(0),
-    headers(""), table(0), forces(0), samples(0), endpoints(0), shifts(0) {}
+    headers(""), table(0), forces(0), samples(0), endpoints(0), shifts(0), outputunit(1.0) {}
 
 DRRForceGrid::DRRForceGrid(const std::vector<DRRAxis> &p_dimensions,
                            const std::string &p_suffix, bool initializeTable)
@@ -132,6 +132,7 @@ DRRForceGrid::DRRForceGrid(const std::vector<DRRAxis> &p_dimensions,
   forceSize = sampleSize * ndims;
   forces.resize(forceSize, 0.0);
   samples.resize(sampleSize, 0);
+  outputunit = 1.0;
   // For 1D pmf
   if (ndims == 1) {
     endpoints.resize(dimensions[0].nbins + 1, 0);
@@ -301,24 +302,20 @@ DRRForceGrid::getCountsLogDerivative(const std::vector<double> &pos) const {
 // }
 
 void DRRForceGrid::write1DPMF(std::string filename) const {
-  std::stringstream ssp;
-  std::fstream fspmf;
   filename += suffix + ".pmf";
+  FILE *ppmf;
+  ppmf = fopen(filename.c_str(), "w");
   const double w = dimensions[0].getWidth();
   double pmf = 0;
-  ssp << std::fixed << std::setprecision(OUTPUTPRECISION) << endpoints[0] << ' '
-      << pmf << '\n';
+  fprintf(ppmf, "%.9f %.9f\n", endpoints[0], pmf);
   for (size_t i = 0; i < dimensions[0].nbins; ++i) {
     std::vector<double> pos(1, 0);
     pos[0] = table[0][i];
     const std::vector<double> f = getGradient(pos, true);
-    ssp << endpoints[i + 1];
     pmf += f[0] * w;
-    ssp << ' ' << pmf << '\n';
+    fprintf(ppmf, "%.9f %.9f\n", endpoints[i + 1], pmf);
   }
-  fspmf.open(filename.c_str(), std::ios_base::out);
-  fspmf.write(ssp.str().c_str(), ssp.str().length());
-  fspmf.close();
+  fclose(ppmf);
 }
 
 // Write the gradients to a .count file.
@@ -342,34 +339,36 @@ void DRRForceGrid::write1DPMF(std::string filename) const {
 // }
 
 void DRRForceGrid::writeAll(const std::string &filename) const {
-  std::stringstream ssc, ssg;
-  std::fstream fscount, fsgrad;
   std::string countname = filename + suffix + ".count";
   std::string gradname = filename + suffix + ".grad";
-  ssc << headers << std::left << std::fixed
-      << std::setprecision(OUTPUTPRECISION);
-  ssg << headers << std::left << std::fixed
-      << std::setprecision(OUTPUTPRECISION);
   std::vector<double> pos(ndims, 0);
+  FILE *pGrad, *pCount;
+  pGrad = fopen(gradname.c_str(), "w");
+  pCount = fopen(countname.c_str(), "w");
+  char *buffer1, *buffer2;
+  buffer1 = (char*)malloc((sizeof(double))*sampleSize*ndims);
+  buffer2 = (char*)malloc((sizeof(double))*sampleSize*ndims);
+  setvbuf(pGrad, buffer1, _IOFBF, (sizeof(double))*sampleSize*ndims);
+  setvbuf(pCount, buffer2, _IOFBF, (sizeof(double))*sampleSize*ndims);
+  fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pGrad);
+  fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pCount);
   for (size_t i = 0; i < sampleSize; ++i) {
     for (size_t j = 0; j < ndims; ++j) {
       pos[j] = table[j][i];
-      ssc << ' ' << table[j][i];
-      ssg << ' ' << table[j][i];
+      fprintf(pGrad, " %.9f", table[j][i]);
+      fprintf(pCount, " %.9f", table[j][i]);
     }
-    ssc << ' ' << getCount(pos, true) << '\n';
-    //     ssc << ' ' << samples[i] << '\n';
+    fprintf(pCount, " %lu\n", getCount(pos, true));
     std::vector<double> f = getGradient(pos, true);
-    for (const auto &i : f)
-      ssg << ' ' << i;
-    ssg << '\n';
+    for (size_t j = 0; j < ndims; ++j) {
+      fprintf(pGrad, " %.9f", (f[j] / outputunit));
+    }
+    fprintf(pGrad, "\n");
   }
-  fscount.open(countname.c_str(), std::ios_base::out);
-  fsgrad.open(gradname.c_str(), std::ios_base::out);
-  fscount.write(ssc.str().c_str(), ssc.str().length());
-  fsgrad.write(ssg.str().c_str(), ssg.str().length());
-  fscount.close();
-  fsgrad.close();
+  fclose(pGrad);
+  fclose(pCount);
+  free(buffer1);
+  free(buffer2);
   if (ndims == 1) {
     write1DPMF(filename);
   }
diff --git a/src/drr/DRR.h b/src/drr/DRR.h
index 32c8d6528cb6b312f5126bcc5b430e9e84350b4e..e516678cd0041575071033dcddbda8fb93e78e91 100644
--- a/src/drr/DRR.h
+++ b/src/drr/DRR.h
@@ -195,7 +195,9 @@ public:
                                     const std::vector<DRRAxis> &dB);
   /// Get suffix
   std::string getSuffix() const { return suffix; }
-  // Destructor
+  /// Set unit for .grad output
+  void setOutputUnit(double unit) { outputunit = unit; }
+  /// Destructor
   virtual ~DRRForceGrid() {}
 
 protected:
@@ -227,6 +229,8 @@ protected:
   /// The abf_intergrate program has precision requirement.
   /// I test 9 and it just works.
   static const size_t OUTPUTPRECISION = 9;
+  /// For set different output units
+  double outputunit;
 
   /// Miscellaneous helper functions
   static size_t index1D(const DRRAxis &c, double x);
@@ -270,6 +274,7 @@ protected:
     }
     fillTable(mp);
     headers = ss.str();
+    outputunit = 1.0;
     // For 1D pmf
     if (ndims == 1) {
       endpoints.resize(dimensions[0].nbins + 1, 0);
diff --git a/src/drr/drrtool.cpp b/src/drr/drrtool.cpp
index 41f23fff731d5de37db372c919a312a1a946442b..08bdd0ec9bd5e8011400474517def8097e9c4dc9 100644
--- a/src/drr/drrtool.cpp
+++ b/src/drr/drrtool.cpp
@@ -46,13 +46,13 @@ namespace drr {
 
 The following command will extract .grad and .count files.
 \verbatim
-plumed drrtool --extract eabf.drrstate
+plumed drr_tool --extract eabf.drrstate
 \endverbatim
 
 The following command will merge windows of two .drrstate file, and output the
 .grad and .count files.
 \verbatim
-plumed drrtool --merge win1.drrstate,win2.drrstate
+plumed drr_tool --merge win1.drrstate,win2.drrstate
 \endverbatim
 
 After getting the .grad and .count file, you can do numerical integration by
@@ -61,6 +61,8 @@ https://github.com/Colvars/colvars/tree/master/colvartools
 \verbatim
 abf_integrate eabf.czar.grad
 \endverbatim
+\note
+The abf_integrate in colvartools is in kcal/mol, so it may be better to use --units kcal/mol when running drr_tool
 
 */
 //+ENDPLUMEDOC
@@ -76,6 +78,7 @@ public:
 
 private:
   bool verbosity;
+  Units units;
   const std::string suffix{".drrstate"};
 };
 
@@ -85,6 +88,7 @@ void drrtool::registerKeywords(Keywords &keys) {
   CLTool::registerKeywords(keys);
   keys.add("optional", "--extract", "Extract drrstate file(s)");
   keys.add("optional", "--merge", "Merge eABF windows");
+  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");
   keys.addFlag("-v", false, "Verbose output");
 }
 
@@ -96,6 +100,9 @@ drrtool::drrtool(const CLToolOptions &co) : CLTool(co) {
 int drrtool::main(FILE *in, FILE *out, Communicator &pc) {
   parseFlag("-v", verbosity);
   std::vector<std::string> stateFilesToExtract;
+  std::string unitname;
+  parse("--units",unitname);
+  units.setEnergy( unitname );
   bool doextract = parseVector("--extract", stateFilesToExtract);
   if (doextract) {
     extractdrr(stateFilesToExtract);
@@ -124,7 +131,10 @@ void drrtool::extractdrr(const std::vector<std::string> &filename) {
     ia >> step >> fict >> vfict >> vfict_laststep >> ffict >> abfgrid >>
        czarestimator;
     in.close();
+    abfgrid.setOutputUnit(units.getEnergy());
+    czarestimator.setOutputUnit(units.getEnergy());
     if (verbosity) {
+      std::cout << "Output units factor: " << units.getEnergy() << '\n';
       std::cout << "Dumping information of extended variables..." << '\n';
       std::cout << "Step: " << step << '\n';
       for (size_t i = 0; i < fict.size(); ++i) {
@@ -168,6 +178,8 @@ void drrtool::mergewindows(const std::vector<std::string> &filename) {
     CZAR czarestimator;
     ia >> step >> fict >> vfict >> vfict_laststep >> ffict >> abfgrid >>
        czarestimator;
+    abfgrid.setOutputUnit(units.getEnergy());
+    czarestimator.setOutputUnit(units.getEnergy());
     abfs.push_back(abfgrid);
     czars.push_back(czarestimator);
     in.close();