Implement reflecting boundary conditions for metadynamics bias

doc/colvars-code-refs.bib

@@ -156,6 +156,7 @@ @article{Fiorin2020
 % Scripted functions (Tcl)
 % ABMD bias
 % Updated multiple-walker ABF implementation
+% Metadynamics Hills reflection at the boundaries
 @article{Fiorin2024,
   author = {Fiorin, Giacomo and Marinelli, Fabrizio and Forrest, Lucy R. and Chen, Haochuan and Chipot, Christophe and Kohlmeyer, Axel and Santuz, Hubert and H{\'e}nin, J{\'e}rôme},
   title = {Expanded Functionality and Portability for the Colvars Library},

doc/colvars-refman-main.tex

@@ -6249,66 +6249,6 @@
 where $t_{e}$ is the time after which the bias potential grows (approximately) evenly during the simulation and $t_{tot}$ is the total simulation time. The free energy calculated according to eq.~\ref{eq:colvars_meta_fes_av} can thus be obtained averaging on time multiple time-dependent free energy estimates, that can be printed out through the keyword \texttt{keepFreeEnergyFiles}. An alternative is to obtain the free energy profiles by summing the hills added during the simulation; the hills trajectory can be printed out by enabling the option \texttt{writeHillsTrajectory}.
 
 
-\cvsubsubsec{Treatment of the PMF boundaries}{sec:colvarbias_meta_boundaries}
-
-In typical scenarios the Gaussian hills of a metadynamics potential are interpolated and summed together onto a grid, which is much more efficient than computing each hill independently at every step (the keyword \refkey{useGrids}{metadynamics|useGrids} is \texttt{on} by default).
-This numerical approximation typically yields negligible errors in the resulting PMF \cite{Fiorin2013}.
-However, due to the finite thickness of the Gaussian function, the metadynamics potential would suddenly vanish each time a variable exceeds its grid boundaries.
-
-To avoid such discontinuity the Colvars metadynamics code will keep an explicit copy of each hill that straddles a grid's boundary, and will use it to compute metadynamics forces outside the grid.
-This measure is taken to protect the accuracy and stability of a metadynamics simulation, except in cases of ``natural'' boundaries (for example, the $[0:180]$ interval of an \texttt{angle} colvar) or when the flags \refkey{hardLowerBoundary}{colvar|hardLowerBoundary} and \refkey{hardUpperBoundary}{colvar|hardUpperBoundary} are explicitly set by the user.
-Unfortunately, processing explicit hills alongside the potential and force grids could easily become inefficient, slowing down the simulation and increasing the state file's size.
-
-In general, it is a good idea to \emph{define a repulsive potential to avoid hills from coming too close to the grid's boundaries}, for example as a \texttt{harmonicWalls} restraint (see \ref{sec:colvarbias_harmonic_walls}).\\
-
-\noindent\textbf{Example:} Using harmonic walls to protect the grid's boundaries.
-\begin{cvexampleinput}
-\-colvar~\{\\
-\-\-~~name~r\\
-\-\-~~distance~\{~...~\}\\
-\-\-~~upperBoundary~15.0\\
-\-\-~~width~0.2\\
-\-\}\\
-\\
-\-metadynamics~\{\\
-\-\-~~name~meta\_r\\
-\-\-~~colvars~r\\
-\-\-~~hillWeight~0.001\\
-\-\-~~hillWidth~2.0\\
-\-\}\\
-\\
-\-harmonicWalls~\{\\
-\-\-~~name~wall\_r\\
-\-\-~~colvars~r\\
-\-\-~~upperWalls~13.0\\
-\-\-~~upperWallConstant~2.0\\
-\-\}
-\end{cvexampleinput}
-
-In the colvar \texttt{r}, the \texttt{distance} function used has a \texttt{lowerBoundary} automatically set to 0 by default, thus the keyword \texttt{lowerBoundary} itself is not mandatory and \texttt{hardLowerBoundary} is set to \texttt{yes} internally.
-However, \texttt{upperBoundary} does not have such a ``natural'' choice of value.
-The metadynamics potential \texttt{meta\_r} will individually process any hill whose center is too close to the \texttt{upperBoundary}, more precisely within fewer grid points than 6 times the Gaussian $\sigma$ parameter plus one.
-It goes without saying that if the colvar \texttt{r} represents a distance between two freely-moving molecules, it will cross this ``threshold'' rather frequently.
-
-In this example, where the value of \texttt{hillWidth} ($2\sigma$) amounts to 2 grid points, the threshold is 6+1 = 7 grid points away from \texttt{upperBoundary}.
-In explicit units, the \texttt{width} of $r$ is $w_r =$ 0.2~\AA, and the threshold is 15.0 - 7$\times$0.2 = 13.6~\AA.
-
-The \texttt{wall\_r} restraint included in the example prevents this: the position of its \texttt{upperWall} is 13~\AA{}, i.e.{} 3 grid points below the buffer's threshold (13.6~\AA).
-For the chosen value of \texttt{upperWallConstant}, the energy of the \texttt{wall\_r} bias at \texttt{r} = $r_{\mathrm{upper}}$ = 13.6~\AA{} is:
-\begin{equation*}
-  E^* = \frac{1}{2} \- k \left(\frac{r - r_{\mathrm{upper}}}{w_r}\right)^2 = \frac{1}{2} \- 2.0 \left(-3\right)^2 = 9~\mathrm{kcal/mol}
-\end{equation*}
-which results in a relative probability $\exp(-E^*/\kappa_{\mathrm{B}}T) \simeq$ $3\times{}10^{-7}$ that \texttt{r} crosses the threshold.
-The probability that \texttt{r} exceeds \texttt{upperBoundary}, which is further away, has also become vanishingly small.
-At that point, you may want to set \texttt{hardUpperBoundary} to \texttt{yes} for \texttt{r}, and let \texttt{meta\_r} know that no special treatment near the grid's boundaries will be needed.
-
-\emph{What is the impact of the wall restraint onto the PMF?} Not a very complicated one: the PMF reconstructed by metadynamics will simply show a sharp increase in free-energy where the wall potential kicks in (\texttt{r}~$>$ 13~\AA{}).
-You may then choose between using the PMF only up until that point and discard the rest, or subtracting the energy of the \texttt{harmonicWalls} restraint from the PMF itself.
-Keep in mind, however, that the statistical convergence of metadynamics may be less accurate where the wall potential is strong.
-
-In summary, although it would be simpler to set the wall's position \texttt{upperWall} and the grid's boundary \texttt{upperBoundary} to the same number, the finite width of the Gaussian hills calls for setting the former strictly within the latter.
-
-
 \cvsubsubsec{Required metadynamics keywords}{sec:colvarbias_meta_basics}
 
 To enable a metadynamics-based calculation, a \texttt{metadynamics \{...\}} block must be included in the Colvars configuration file.
@@ -6379,6 +6319,123 @@
 \end{itemize}
 
 
+\cvsubsubsec{Treatment of the PMF boundaries}{sec:colvarbias_meta_boundaries}
+
+Similar to other PMF calculation methods, it is common for sampling to be restricted to a predefined region of the CV space. This restriction can arise either from the ``natural'' boundaries of the CVs (e.g., the $[0:180]$ interval of an \texttt{angle} collective variable) or from user-imposed limits to focus on a specific region of interest, such as defining a repulsive potential between boundaries using \texttt{harmonicWalls}. In both cases, the finite size of the Gaussian hills causes significant systematic errors near the boundaries \cite{Crespo2010,Fiorin2024}, resulting in artificial spikes and oscillations in the bias potential. Over time, this leads to large forces (also by driving the system towards high-energy regions of \texttt{harmonicWalls} restraints), ultimately generating simulation instabilities.
+To prevent such errors, colvars implements a specific methodology\cite{Fiorin2024} whereby the hills are reflected beyond the boundaries and the biasing forces are set to zero in the direction orthogonal to the boundaries once crossed. This results in a continuous biasing potential that remains constant in the direction orthogonal to the boundaries once crossed \cite{Fiorin2024}. This feature can be enabled by setting the keyword \refkey{useHillsReflection}{metadynamics|useHillsReflection} to \texttt{on} and is applied by default on all non-periodic active CVs, except for non-scalar variables, which, if present, will lead to deactivation of the boundary correction.
+Normally, the Gaussian hills in a metadynamics potential are interpolated and summed on a grid, which yields negligible errors in the PMF \cite{Fiorin2013} and is more efficient than computing each hill independently at every step (the \texttt{useGrids} keyword is \texttt{on} by default). When both \texttt{useGrids} and \refkey{useHillsReflection}{metadynamics|useHillsReflection} are active, the reflective boudary correction is automatically applied at the grid edge`s, preventing discontinuities outside the boundaries.
+It is recommended to activate \refkey{useHillsReflection}{metadynamics|useHillsReflection} both in presence of ``natural'' boundaries or user-defined \texttt{harmonicWalls}. The latter walls should be placed at the same boundaries where the reflective boudary correction is applied (e.g. at the CVs or grid boundaries).
+When \refkey{useHillsReflection}{metadynamics|useHillsReflection} is inactive, Colvars keeps an explicit copy of each hill crossing a grid boundary to compute metadynamics forces outside the grid. This approach slows the simulation and increases the state file size over time, but is necessary to avoid critical issues, like boundary trapping. This feature is disabled when 'natural' boundaries or explicit \texttt{hardLowerBoundary} and \texttt{hardUpperBoundary} flags are set. If \refkey{useHillsReflection}{metadynamics|useHillsReflection} is off, it’s recommended to use \texttt{harmonicWalls} restraints (see \ref{sec:colvarbias_harmonic_walls}) placed at least 8 Gaussian sigmas within the boundaries to minimize inefficiencies. While these measures prevent boundary trapping, they do not address systematic errors, which are fully resolved only by enabling \refkey{useHillsReflection}{metadynamics|useHillsReflection}.
+
+\noindent\textbf{Example:} Using hills reflection to avoid systematic erros at the grid's boundaries and harmonic walls to limit the sampling within those.
+\begin{cvexampleinput}
+\-colvar~\{\\
+\-\-~~name~r\\
+\-\-~~distance~\{~...~\}\\
+\-\-~~upperBoundary~13.0\\
+\-\-~~width~0.2\\
+\-\}\\
+\\
+\-metadynamics~\{\\
+\-\-~~name~meta\_r\\
+\-\-~~colvars~r\\
+\-\-~~hillWeight~0.001\\
+\-\-~~hillWidth~2.0\\
+\-\-~~useHillsReflection~on\\
+\-\}\\
+\\
+\-harmonicWalls~\{\\
+\-\-~~name~wall\_r\\
+\-\-~~colvars~r\\
+\-\-~~upperWalls~13.0\\
+\-\-~~upperWallConstant~2.0\\
+\-\}
+\end{cvexampleinput}
+
+In the colvar \texttt{r}, the \texttt{distance} function used has a \texttt{lowerBoundary} automatically set to 0 by default, thus the keyword \texttt{lowerBoundary} itself is not mandatory and \refkey{hardLowerBoundary}{colvar|hardLowerBoundary} is set to \texttt{yes} internally.
+However, \texttt{upperBoundary} does not have such a ``natural'' choice of value.
+The activation of the reflective boudary correction via \refkey{useHillsReflection}{metadynamics|useHillsReflection} set to \texttt{on}, will lead to hills that are reflected at distances of 0 and 13.
+
+\emph{What is the impact of the reflective boudary correction onto the PMF?} Not a very complicated one: the PMF reconstructed by metadynamics will simply show a slight flattening very close to the boundaries (\texttt{r}~$=$ 0~\AA{} and \texttt{r}~$=$ 13~\AA{}), with the biasing forces becoming negligible. Those forces are completely removed after crossing the boundaries (\texttt{r}~$>$ 13~\AA{}) and the metadynamics bias potential remains constant. Thus, the only biasing forces applied beyond the upper boundary are those from the \texttt{harmonicWalls} restraints.
+
+Although the reflective boudary correction typically leads to a reasonable estimate of the PMF, a more rigorous evaluation can be obtained by rewighting or mean force estimates, such as the Force Correction Analisys Method\cite{Marinelli2021}.
+
+In summary, enabling the \refkey{useHillsReflection}{metadynamics|useHillsReflection} keyword corrects boundary artifacts, ensuring a more stable metadynamics simulation and improved PMF accuracy.
+
+While \refkey{useHillsReflection}{metadynamics|useHillsReflection} is in most cases the only required keyword to activate the reflective boundary correction, a complete list of associated keywords is provided below.
+
+\begin{itemize}
+
+\item %
+  \labelkey{metadynamics|useHillsReflection}
+  \keydef
+    {useHillsReflection}{%
+    \texttt{metadynamics}}{%
+    use reflective boundary correction}{%
+    boolean}{%
+    \texttt{off}}{%
+    This keyword instructs Colvars to use reflective boundary correction.
+    When active, this is applied by default on all non-periodic variables.
+    When \texttt{useGrids} is \texttt{on}, the boundary correction is applied
+    by default at the grid boundaries. Optionally, the CVs and boundaries on
+    which this correction is applied can be specified via the keywords
+    \texttt{reflectionLowLimitUseCVs}, \texttt{reflectionUpLimitUseCVs},
+    \texttt{reflectionLowLimit} and \texttt{reflectionUpLimit}. Currently,
+    this option is implemented for all types of variables except the non-scalar
+    types (\texttt{distanceDir} or \texttt{orientation}), and it becomes inactive
+    if those are active. To enhance efficiency, the boundary correction is applied
+    only to Gaussians that lie within a specified cutoff distance from the boundaries
+    (defaulting to 6 times the Gaussian's sigma, though this can be adjusted using
+    the keyword \texttt{reflectionRange})}
+
+\item %
+  \labelkey{metadynamics|reflectionLowLimitUseCVs}
+  \keydef
+    {reflectionLowLimitUseCVs}{%
+    \texttt{metadynamics}}{% 
+    Collective variables involved}{%
+    space-separated list of colvar names}{%
+    All non-periodic variables}{%
+    This option selects by name all the variables to which a reflective boundary
+    correction is applied at their lower boundaries.}
+
+
+\item %
+  \labelkey{metadynamics|reflectionUpLimitUseCVs}
+  \keydef
+    {reflectionUpLimitUseCVs}{%
+    \texttt{metadynamics}}{% 
+    Collective variables involved}{%
+    space-separated list of colvar names}{%
+    All non-periodic variables}{%
+    This option selects by name all the variables to which a reflective boundary
+    correction is applied at their upper boundaries.}
+
+\item %
+  \labelkey{metadynamics|reflectionLowLimit}
+  \keydef
+    {reflectionLowLimit}{%
+    \texttt{metadynamics}}{%
+    Lower boundaries for the reflective boundary correction (one for each colvar)}{%
+    space-separated list of decimals}{%
+    Grid lower boundaries if \refkey{useGrids}{metadynamics|useGrids} is \texttt{on}}{%
+    This option sets the values of the lower boundaries for the CVs specified in
+    \refkey{reflectionLowLimitUseCVs}{metadynamics|reflectionLowLimitUseCVs}.}
+
+\item %
+  \labelkey{metadynamics|reflectionUpLimit}
+  \keydef
+    {reflectionUpLimit}{%
+    \texttt{metadynamics}}{%
+    Upper boundaries for the reflective boundary correction (one for each colvar)}{%
+    space-separated list of decimals}{%
+    Grid upper boundaries if \refkey{useGrids}{metadynamics|useGrids} is \texttt{on}}{%
+    This option sets the values of the upper boundaries for the CVs specified in
+    \refkey{reflectionUpLimitUseCVs}{metadynamics|reflectionUpLimitUseCVs}.}
+
+\end{itemize}
+
+
 \cvsubsubsec{Output files}{sec:colvarbias_meta_output}
 
 When interpolating grids are enabled (default behavior), the PMF is written by default every \texttt{colvarsRestartFrequency} steps to the file \outputName\texttt{.pmf} in multicolumn text format (\ref{sec:colvar_multicolumn_grid}).

doc/colvars-refman.bib

@@ -84,6 +84,16 @@ @ARTICLE{Crespo2010
   url = {https://link.aps.org/doi/10.1103/PhysRevE.81.055701},
 }
 
+@article{Marinelli2021,
+  author = {Marinelli, Fabrizio and Faraldo-Gómez, José D.},
+  title = {Force-Correction Analysis Method for Derivation of Multidimensional Free-Energy Landscapes from Adaptively Biased Replica Simulations},
+  journal = {J. Chem. Theory Comput.},
+  volume = {17},
+  pages = {6775-6788},
+  year = {2021},
+  url = {https://doi.org/10.1021/acs.jctc.1c00586},
+}
+
 @ARTICLE{Darve2008,
   author = {Eric Darve and David Rodr{\'i}guez-G{\'o}mez and Andrew Pohorille},
   title = {Adaptive biasing force method for scalar and vector free energy calculations},