New tool poisson_integrator_conv

colvartools/Makefile

@@ -2,26 +2,36 @@ SHELL=/bin/sh
 #########################
 # adjust as needed
 # settings for linux with GCC
+COLVARS_SRC=../src
 CXX=g++
-CXXFLAGS=-Wall -O2
+CXXFLAGS=-Wall -O2 -I$(COLVARS_SRC)
 EXT=
 # settings for mingw cross-compiler to windows (32-bit)
 #CXX=i688-w64-mingw64-g++
 #CXXFLAGS=-Wall -O2
 #EXT=.exe
 #########################
 
-all: abf_integrate$(EXT)
+all: poisson_integrator$(EXT) abf_integrate$(EXT)
 
 clean:
-	-rm *~ *.o abf_integrate$(EXT) *.exe
+	-rm *~ *.o abf_integrate$(EXT) poisson_integrator$(EXT) *.exe
 
 abf_integrate$(EXT): abf_integrate.o abf_data.o
 	$(CXX) -o $@ $(CXXFLAGS) $^
 
+poisson_integrator$(EXT): poisson_integrator.o libcolvars
+	$(CXX) -o $@ $(CXXFLAGS) poisson_integrator.o $(COLVARS_SRC)/libcolvars.a
+
+poisson_integrator_conv: poisson_integrator_conv.o libcolvars
+	$(CXX) -o $@ $(CXXFLAGS) poisson_integrator_conv.o $(COLVARS_SRC)/libcolvars.a
+
 %.o: %.cpp
 	$(CXX) -o $@ -c $(CXXFLAGS) $<
 
+libcolvars:
+	$(MAKE) -C $(COLVARS_SRC) libcolvars.a
+
 # dependencies
 abf_integrate.o: abf_integrate.cpp abf_data.h
 abf_data.o: abf_data.cpp abf_data.h

colvartools/README.md

@@ -5,7 +5,9 @@ This directory contains both standalone tools and Colvars scripts.
 ## Standalone tools
 | File name | Summary |
 | ------------- | ------------- |
-| **abf_integrate** | Post-process gradient files produced by ABF and related methods, to generate a PMF. Superseded by builtin integration for dimensions 2 and 3, still needed for higher-dimension PMFs. Build using the provided **Makefile**.|
+| **abf_integrate** | Process free energy gradient from ABF/TI to generate a free energy surface using a legacy MCMC procedure. Superseded by Poisson integration (builtin or standalone) for dimensions 2 and 3, still needed for higher-dimension PMFs. Build using the provided **Makefile**.|
+| **poisson_integrator** | Process free energy gradient from ABF/TI to generate a free energy surface as the solution of a Poisson equation. Build using the provided **Makefile**.|
+| **poisson_integrator_conv** | Process free energy gradient from ABF/TI to generate a free energy surface as the solution of a Poisson equation. Build using the provided **Makefile**.|
 | **noe_to_colvars.py** | Parse an X-PLOR style list of assign commands for NOE restraints.|
 | **plot_colvars_traj.py** | Select variables from a Colvars trajectory file and optionally plot them as a 1D graph as a function of time or of one of the variables.|
 | **quaternion2rmatrix.tcl** | As the name says.|
@@ -16,4 +18,4 @@ This directory contains both standalone tools and Colvars scripts.
 | File name | Summary |
 | ------------- | ------------- |
 | **abmd.tcl** | Adiabatic Biased MD after Marchi & Ballone JCP 1999 / Paci & Karplus JMB 1999; implemented in 20 lines of Tcl.|
-| **pathCV.tcl** | Path-based collective variables after Branduardi et al. (2007). Optimized implementation that calculates only the necessary distances.|
+| **pathCV.tcl** | Path-based collective variables after Branduardi et al. (2007). Optimized implementation that calculates only the necessary distances.|

colvartools/poisson_integrator.cpp

@@ -8,7 +8,7 @@
 int main (int argc, char *argv[]) {
 
   if (argc < 2) {
-    std::cerr << "One argument needed: file name.\n";
+    std::cerr << "\n\nOne argument needed: gradient multicol file name.\n";
     return 1;
   }
 
@@ -17,18 +17,24 @@ int main (int argc, char *argv[]) {
 
   std::string gradfile (argv[1]);
   std::shared_ptr<colvar_grid_gradient> grad_ptr = std::make_shared<colvar_grid_gradient>(gradfile);
+  if (cvm::get_error()) { return -1; }
 
-  int itmax = 1000;
+  int itmax = 10000;
   cvm::real err;
-  cvm::real tol = 1e-6;
+  cvm::real tol = 1e-8;
 
   integrate_potential potential(grad_ptr);
   potential.set_div();
   potential.integrate(itmax, tol, err);
   potential.set_zero_minimum();
 
-  potential.write_multicol(std::string(gradfile + ".int"),
-                           "integrated potential");
+  if (potential.num_variables() < 3) {
+    std::cout << "\nWriting integrated potential in multicol format to " + gradfile + ".int\n";
+    potential.write_multicol(std::string(gradfile + ".int"), "integrated potential");
+  } else { // Write 3D grids to more convenient DX format
+    std::cout << "\nWriting integrated potential in OpenDX format to " + gradfile + ".int.dx\n";
+    potential.write_opendx(std::string(gradfile + ".int.dx"), "integrated potential");
+  }
 
   delete colvars;
   return 0;

colvartools/poisson_integrator_conv.cpp

@@ -0,0 +1,64 @@
+#include <iostream>
+
+#include "colvargrid.h"
+#include "colvarproxy.h"
+
+// Integrate provided gradients while monitoring convergence towards a provided scalar grid
+// (typically the result of a previous integration)
+
+int main (int argc, char *argv[])
+{
+  colvarproxy *proxy = new colvarproxy();
+  colvarmodule *colvars = new colvarmodule(proxy);
+
+  if (argc < 2) {
+    std::cerr << "\n\nOne argument needed: gradient multicol file name.\n";
+    return 1;
+  }
+
+  std::string gradfile (argv[1]);
+  std::shared_ptr<colvar_grid_gradient> grad_ptr = std::make_shared<colvar_grid_gradient>(gradfile);
+  if (cvm::get_error()) { return -1; }
+
+  cvm::real err = 1.;
+  cvm::real tol = 1e-10;
+
+  integrate_potential potential(grad_ptr);
+  potential.set_div();
+
+  // Load reference
+  colvar_grid_scalar ref(gradfile + ".ref");
+  if (cvm::get_error()) { return -1; }
+
+  if (ref.number_of_points() != potential.number_of_points()) {
+    cvm::error("Reference grid has wrong number of points: " + cvm::to_str(ref.number_of_points()) + "\n");
+    return -1;
+  }
+
+  // Monitor convergence
+  int rounds = 100;
+  int steps_per_round = 10;
+
+  std::cout << 0 << " " << 0 << " " << ref.grid_rmsd(potential) << std::endl;
+
+  for (int i = 0; i < rounds && err > tol; i++) {
+    potential.reset(0.);
+    potential.integrate(steps_per_round * (i+1), tol, err);
+    potential.set_zero_minimum();
+
+    std::cout << (i+1)*steps_per_round << " " <<  err << " " << ref.grid_rmsd(potential) << std::endl;
+
+    char buff[100];
+    snprintf(buff, sizeof(buff), "%04i", steps_per_round * (i+1));
+  }
+
+  if (potential.num_variables() < 3) {
+    std::cout << "\nWriting integrated potential in multicol format to " + gradfile + ".int\n";
+    potential.write_multicol(std::string(gradfile + ".int"), "integrated potential");
+  } else { // Write 3D grids to more convenient DX format
+    std::cout << "\nWriting integrated potential in OpenDX format to " + gradfile + ".int.dx\n";
+    potential.write_opendx(std::string(gradfile + ".int.dx"), "integrated potential");
+  }
+  delete colvars;
+  return 0;
+}

src/colvargrid.cpp

@@ -142,6 +142,62 @@ colvar_grid_scalar::colvar_grid_scalar(std::vector<colvar *> &colvars, bool marg
 {
 }
 
+
+colvar_grid_scalar::colvar_grid_scalar(std::string const &filename)
+  : colvar_grid<cvm::real>(),
+    samples(NULL)
+{
+  std::istream &is = cvm::main()->proxy->input_stream(filename, "multicol scalar file");
+  if (!is) {
+    return;
+  }
+
+  // Data in the header: nColvars, then for each
+  // xiMin, dXi, nPoints, periodic flag
+
+  std::string  hash;
+  size_t i;
+
+  if ( !(is >> hash) || (hash != "#") ) {
+    cvm::error("Error reading grid at position "+
+                cvm::to_str(static_cast<size_t>(is.tellg()))+
+                " in stream(read \"" + hash + "\")\n");
+    return;
+  }
+
+  is >> nd;
+  mult = 1;
+  std::vector<cvm::real> lower_in(nd), widths_in(nd);
+  std::vector<int>       nx_in(nd);
+  std::vector<int>       periodic_in(nd);
+
+  for (i = 0; i < nd; i++ ) {
+    if ( !(is >> hash) || (hash != "#") ) {
+      cvm::error("Error reading grid at position "+
+                  cvm::to_str(static_cast<size_t>(is.tellg()))+
+                  " in stream(read \"" + hash + "\")\n");
+      return;
+    }
+
+    is >> lower_in[i] >> widths_in[i] >> nx_in[i] >> periodic_in[i];
+  }
+
+  this->setup(nx_in, 0., mult);
+
+  widths = widths_in;
+
+  for (i = 0; i < nd; i++ ) {
+    lower_boundaries.push_back(colvarvalue(lower_in[i]));
+    periodic.push_back(static_cast<bool>(periodic_in[i]));
+  }
+
+  // Reset the istream for read_multicol, which expects the whole file
+  is.clear();
+  is.seekg(0);
+  read_multicol(is);
+  cvm::main()->proxy->close_input_stream(filename);
+}
+
 colvar_grid_scalar::~colvar_grid_scalar()
 {
 }
@@ -351,12 +407,11 @@ colvar_grid_gradient::colvar_grid_gradient(std::vector<colvar *> &colvars, std::
 }
 
 
-colvar_grid_gradient::colvar_grid_gradient(std::string &filename)
+colvar_grid_gradient::colvar_grid_gradient(std::string const &filename)
   : colvar_grid<cvm::real>(),
     samples(NULL)
 {
-  std::istream &is = cvm::main()->proxy->input_stream(filename,
-                                                      "gradient file");
+  std::istream &is = cvm::main()->proxy->input_stream(filename, "gradient file");
   if (!is) {
     return;
   }

src/colvargrid.h

@@ -1262,6 +1262,9 @@ class colvar_grid_scalar : public colvar_grid<cvm::real>
   colvar_grid_scalar(std::vector<colvar *> &colvars,
                      bool add_extra_bin = false);
 
+  /// Constructor from a multicol file
+  colvar_grid_scalar(std::string const &filename);
+
   /// Accumulate the value
   inline void acc_value(std::vector<int> const &ix,
                         cvm::real const &new_value,
@@ -1573,7 +1576,7 @@ class colvar_grid_gradient : public colvar_grid<cvm::real>
   colvar_grid_gradient(std::vector<colvar *>  &colvars);
 
   /// Constructor from a multicol file
-  colvar_grid_gradient(std::string &filename);
+  colvar_grid_gradient(std::string const &filename);
 
   /// Constructor from a vector of colvars and a pointer to the count grid
   colvar_grid_gradient(std::vector<colvar *> &colvars, std::shared_ptr<colvar_grid_count> samples_in);