ABF using kernelized gradient grid

colvartools/poisson_integrator.cpp

@@ -1,42 +1,104 @@
 #include <iostream>
 #include <fstream>
+#include <sys/stat.h>
 
 #include "colvargrid.h"
 #include "colvargrid_integrate.h"
 #include "colvarproxy.h"
+#include "CLI11.hpp"
 
+int main(int argc, char *argv[])
+{
+    bool weighted = false;
+    bool debug = false;
+    int itmax = 1000;
+    cvm::real err;
+    cvm::real tol = 1e-3;
 
-int main (int argc, char *argv[]) {
+    std::string gradfile;
+    std::string countfile;
 
-  if (argc < 2) {
-    std::cerr << "\n\nOne argument needed: gradient multicol file name.\n";
-    return 1;
-  }
+    CLI::App app("Colvars gradient grid integrator");
 
-  colvarproxy *proxy = new colvarproxy();
-  proxy->colvars = new colvarmodule(proxy); // This could be omitted if we used the colvarproxy_stub class
+    app.add_option("gradfile", gradfile, "Gradient multicol file name")
+        ->required()
+        ->check(CLI::ExistingFile);
 
-  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; }
+    app.add_option("countfile", countfile, "Count file name");
+    app.add_flag("--weighted,!--unweighted", weighted, "Enable or disable weighting")
+        ->default_val(weighted);
 
-  int itmax = 10000;
-  cvm::real err;
-  cvm::real tol = 1e-8;
+    app.add_option("--tol", tol, "Tolerance value");
+    app.add_option("--itmax", itmax, "Maximum iterations");
+    app.add_flag("--debug", debug, "Write debugging output and files");
 
-  colvargrid_integrate fes(grad_ptr);
-  fes.set_div();
-  fes.integrate(itmax, tol, err);
-  fes.set_zero_minimum();
+    CLI11_PARSE(app, argc, argv);
 
-  if (fes.num_variables() < 3) {
-    std::cout << "\nWriting integrated free energy in multicol format to " + gradfile + ".int\n";
-    fes.write_multicol(std::string(gradfile + ".int"), "integrated free energy");
-  } else { // Write 3D grids to more convenient DX format
-    std::cout << "\nWriting integrated free energy in OpenDX format to " + gradfile + ".int.dx\n";
-    fes.write_opendx(std::string(gradfile + ".int.dx"), "integrated free energy");
-  }
+    // (Moved after parsing so we don't allocate memory if help/error is printed)
+    colvarproxy *proxy = new colvarproxy();
+    proxy->colvars = new colvarmodule(proxy);
+    std::shared_ptr<colvar_grid_scalar> count_ptr;
 
-  delete proxy;
-  return 0;
+    // Deduce countfile if not provided by user
+    if (countfile.empty())
+    {
+        size_t pos = gradfile.rfind(std::string(".czar.grad"));
+        if (pos != std::string::npos)
+        {
+            countfile = gradfile.substr(0, pos) + ".zcount";
+        }
+        else
+        {
+            pos = gradfile.rfind(std::string(".grad"));
+            if (pos != std::string::npos)
+            {
+                countfile = gradfile.substr(0, pos) + ".count";
+            }
+        }
+    }
+
+    if (countfile.size())
+    {
+        struct stat buffer;
+        if (stat(countfile.c_str(), &buffer) == 0)
+        {
+            std::cout << "Found associated count file " << countfile << ", reading...\n";
+            count_ptr.reset(new colvar_grid_scalar(countfile));
+            if (!count_ptr || count_ptr->nd == 0)
+            {
+                // catch constructor failure
+                cvm::error("Error reading count grid.");
+                return cvm::get_error();
+            }
+            if (debug) count_ptr->write_multicol("counts.dat");
+        }
+    }
+
+    std::cout << "Reading gradient file " << gradfile << std::endl;
+    std::shared_ptr<colvar_grid_gradient> grad_ptr = std::make_shared<colvar_grid_gradient>(gradfile, count_ptr);
+
+    if (!grad_ptr || grad_ptr->nd == 0)
+    {
+        // catch constructor failure
+        cvm::error("Error reading gradient grid.");
+        return cvm::get_error();
+    }
+
+    if (debug) grad_ptr->write_multicol("gradient_in.dat");
+
+    colvargrid_integrate fes(grad_ptr, weighted);
+    fes.integrate(itmax, tol, err, true);
+    fes.set_zero_minimum();
+    if (fes.num_variables() < 3)
+    {
+        std::cout << "\nWriting integrated fes in multicol format to " + gradfile + ".int\n";
+        fes.write_multicol(std::string(gradfile + ".int"), "integrated fes");
+    } else {
+        // Write 3D grids to more convenient DX format
+        std::cout << "\nWriting integrated free energy in OpenDX format to " + gradfile + ".int.dx\n";
+        fes.write_opendx(std::string(gradfile + ".int.dx"), "integrated free energy");
+    }
+
+    delete proxy;
+    return 0;
 }

colvartools/poisson_integrator_conv.cpp

@@ -17,14 +17,15 @@ int main (int argc, char *argv[]) {
   proxy->colvars = new colvarmodule(proxy); // This could be omitted if we used the colvarproxy_stub class
 
   std::string gradfile (argv[1]);
-  std::shared_ptr<colvar_grid_gradient> grad_ptr = std::make_shared<colvar_grid_gradient>(gradfile);
+  std::shared_ptr<colvar_grid_count> const nullpointer = nullptr;
+  std::shared_ptr<colvar_grid_gradient> grad_ptr = std::make_shared<colvar_grid_gradient>(gradfile, nullpointer);
   if (cvm::get_error()) { return -1; }
 
   cvm::real err = 1.;
   cvm::real tol = 1e-10;
 
   colvargrid_integrate fes(grad_ptr);
-  fes.set_div();
+  fes.set_unweighted_div();
 
   // Load reference
   colvar_grid_scalar ref(gradfile + ".ref");

doc/colvars-code-refs.bib

@@ -242,6 +242,13 @@ @Article{Henin2021
   url     = {https://doi.org/10.1021/acs.jctc.1c00593},
 }
 
+% Poisson equation resolution using discrete difference scheme
+@misc{ChenFDM101,
+  author = {Long Chen},
+  title = {Finite difference methods for Poisson equation},
+  url = {https://www.math.uci.edu/~chenlong/226/FDM.pdf},
+}
+
 % VMD engine
 @article{Humphrey1996,
   title = {{VMD}: visual molecular dynamics},

doc/colvars-refman-main.tex

@@ -5862,6 +5862,27 @@
   }
 \end{itemize}
 
+\cvsubsubsec{Kernel grid ABF}
+ABF stores energy gradient sums and counts on two separate grids in order to estimate the free energy gradient. If the
+grids are big, whether for accuracy or because of the dimensionality curse filling them will be long. In order to reduce
+this time we use kernel grid to leverage the free energy function regularity. At each time step $t$ : the system is in
+the state $(x(t), \xi(t), \nabla U(x(t))$, with $x(t)$ the atomic coordinate, $\xi(t)$ the value of the collective
+variables, $\nabla U(x(t))$ the value of the energy of the system. With kernel grid, instead of updating only the cell of
+each grid corresponding to $\xi(t)$ we'll also update the surrounding cells with a weight based on their distance (kernel)
+to $\xi(t)$. Now in order to estimate the free energy gradient at a bin, the bin averaging will actually corresponds to a
+weighted average of the energy gradients of nearby datapoints seen during simulation.
+
+The weights are calculated using a kernel of the form
+\begin{equation}
+  K(\tilde{\xi}, \tilde{\xi_{cell}}) = \frac{1}{Z}e^{\frac{||\tilde{\xi} - \tilde{\xi_{cell}}||^2}{\sigma}}
+\end{equation}
+Where \tilde{\xi} is the coordinate of the collective variable in the grid. i.e.
+\begin{equation}
+  \tilde{\xi^i} = \frac{\xi^i - boundary_{lower}^i}{boundary_{upper}^i - boundary_{lower}^i}
+\end{equation}
+In order to activate kernel grid \refkey{smoothing}{abf|smoothing} to set the value of
+the kernel's bandwith $\sigma$ to a value, if the value is > 0 kernel grid is in use.
+
 \cvsubsubsec{Multiple-walker ABF}{sec:colvarbias_abf_shared}
 \label{sec:mw-ABF}
 
@@ -5947,19 +5968,19 @@
 {\texttt{yes}}
 {
   This option, active by default when the dimension of the colvar space is 3 or less, enables the calculation of an integrated free energy surface every time ABF output files are written.
-  In dimension 2 or 3, integration is performed by solving a Poisson equation:~\cite{Henin2021}
+  In dimension 2 or 3, by default integration is performed by solving an unweighted Poisson equation:~\cite{Henin2021}
   \begin{equation}
    \nabla^2 A_t = \nabla \cdot G_t
   \end{equation}
-  wehere $G_t$ is the estimated gradient at time $t$, and $A_t$ is corresponding free energy surface.
+  where $G_t$ is the estimated gradient at time $t$, and $A_t$ is corresponding free energy surface.
   The free energy surface is written under the file name \texttt{<outputName>.pmf},
   in a plain text format (see \ref{sec:colvar_multicolumn_grid}) that can be read by most data plotting
   and analysis programs (e.g.{} Gnuplot).
   Periodic boundary conditions are applied to periodic coordinates, and Neumann boundary
   conditions otherwise (imposed free energy gradient at the boundary of the domain).
   The grid used for free energy discretization is extended by one point along
-  non-periodic coordinates, but not along periodic coordinates. See ref.~\cite{Henin2021}
-  for details.
+  non-periodic coordinates, but not along periodic coordinates. See ref.~\cite{Henin2021}, ~\cite{ChenFDM101}
+  for implementation details.
 }
 
 \item \keydef{integrateTol}{\texttt{abf}}{%
@@ -5981,6 +6002,22 @@
   is stopped when the number of iterations reaches \texttt{integrateMaxIterations},
   unless the RMS error has reached \texttt{integrateTol} before.
 }
+\item \keydef{integrateWeighted}{\texttt{abf}}{%
+  Use weighted method to integrate free energy surface from ABF gradients}
+{boolean}
+{\texttt{no}}{
+  This option is disabled by default. If activated, the calculation of the free energy will be performed using weighted
+  integration. This means solving the following weighted Poisson equation:
+    \begin{equation}
+      \nabla\cdot \theta \nabla A_t = \nabla \cdot G_t.
+    \end{equation}
+  With $\theta$, the weights, corresponding to how much an area has been visited during the simulation. The weights values
+  are all strictly superior to 0 so the problem is well-posed. For any cell with zero visits, its weight is assigned to
+  the value corresponding to the 10th percentile of the visit-count distribution. Moreover, visit counts are capped at
+  the 50th percentile measured from the upper tail (i.e., the 50th percentile from the lower tail). The choice of the
+  percentile is hard code. The goal is that the most visited cells which corresponds to the cells where the gradient is
+  the best known have more importance in the free energy calculation.
+  }
 \end{itemize}
 
 

misc_interfaces/stubs/colvarproxy_stub.cpp

@@ -74,7 +74,7 @@ bool colvarproxy_stub::total_forces_enabled() const
 
 bool colvarproxy_stub::total_forces_same_step() const
 {
-  return total_force_requested;
+  return true;
 }