diff --git a/colvartools/Makefile b/colvartools/Makefile index de5f6113c..6a55678d9 100644 --- a/colvartools/Makefile +++ b/colvartools/Makefile @@ -2,8 +2,9 @@ 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++ @@ -11,17 +12,26 @@ EXT= #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 diff --git a/colvartools/README.md b/colvartools/README.md index 8344961b3..fff5ffc95 100644 --- a/colvartools/README.md +++ b/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.| \ No newline at end of file +| **pathCV.tcl** | Path-based collective variables after Branduardi et al. (2007). Optimized implementation that calculates only the necessary distances.| diff --git a/colvartools/poisson_integrator.cpp b/colvartools/poisson_integrator.cpp index 05d9b9b61..4536e1f11 100644 --- a/colvartools/poisson_integrator.cpp +++ b/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 grad_ptr = std::make_shared(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; diff --git a/colvartools/poisson_integrator_conv.cpp b/colvartools/poisson_integrator_conv.cpp new file mode 100644 index 000000000..147c1f36a --- /dev/null +++ b/colvartools/poisson_integrator_conv.cpp @@ -0,0 +1,64 @@ +#include + +#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 grad_ptr = std::make_shared(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; +} diff --git a/src/colvargrid.cpp b/src/colvargrid.cpp index 11693a758..24077f60a 100644 --- a/src/colvargrid.cpp +++ b/src/colvargrid.cpp @@ -142,6 +142,62 @@ colvar_grid_scalar::colvar_grid_scalar(std::vector &colvars, bool marg { } + +colvar_grid_scalar::colvar_grid_scalar(std::string const &filename) + : colvar_grid(), + 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(is.tellg()))+ + " in stream(read \"" + hash + "\")\n"); + return; + } + + is >> nd; + mult = 1; + std::vector lower_in(nd), widths_in(nd); + std::vector nx_in(nd); + std::vector 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(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(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 &colvars, std:: } -colvar_grid_gradient::colvar_grid_gradient(std::string &filename) +colvar_grid_gradient::colvar_grid_gradient(std::string const &filename) : colvar_grid(), 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; } diff --git a/src/colvargrid.h b/src/colvargrid.h index 4cbbb1096..d452d56f3 100644 --- a/src/colvargrid.h +++ b/src/colvargrid.h @@ -1262,6 +1262,9 @@ class colvar_grid_scalar : public colvar_grid colvar_grid_scalar(std::vector &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 const &ix, cvm::real const &new_value, @@ -1573,7 +1576,7 @@ class colvar_grid_gradient : public colvar_grid colvar_grid_gradient(std::vector &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 &colvars, std::shared_ptr samples_in);