diff --git a/doc/colvars-code-refs.bib b/doc/colvars-code-refs.bib
index 8e2570f0e..6b96b1fc2 100644
--- a/doc/colvars-code-refs.bib
+++ b/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},
diff --git a/doc/colvars-refman-main.tex b/doc/colvars-refman-main.tex
index 256da757a..1e200a9fe 100644
--- a/doc/colvars-refman-main.tex
+++ b/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}).
diff --git a/doc/colvars-refman.bib b/doc/colvars-refman.bib
index 8ae6dbc64..2bbb3532d 100644
--- a/doc/colvars-refman.bib
+++ b/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},
diff --git a/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.colvars.out b/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.colvars.out
new file mode 100644
index 000000000..3b313108e
--- /dev/null
+++ b/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.colvars.out
@@ -0,0 +1,207 @@
+colvars: ----------------------------------------------------------------------
+colvars: Please cite Fiorin et al, Mol Phys 2013:
+colvars:   https://doi.org/10.1080/00268976.2013.813594
+colvars: as well as all other papers listed below for individual features used.
+colvars: SMP parallelism is enabled (num threads = 4).
+colvars: This version was built with the C++11 standard or higher.
+colvars: Redefining the Tcl "cv" command to the new script interface.
+colvars: The restart output state file will be "test.tmp.colvars.state".
+colvars: The final output state file will be "test.colvars.state".
+colvars: ----------------------------------------------------------------------
+colvars: Reading new configuration from file "test.in":
+colvars: # units = "" [default]
+colvars: # indexFile = "index.ndx"
+colvars: The following index groups are currently defined:
+colvars:   Protein (104 atoms)
+colvars:   Protein_noH (51 atoms)
+colvars:   Protein_Backbone (40 atoms)
+colvars:   Protein_C-alpha (10 atoms)
+colvars:   RMSD_atoms (10 atoms)
+colvars:   Protein_C-alpha_1_2 (2 atoms)
+colvars:   Protein_C-alpha_9_10 (2 atoms)
+colvars:   Protein_C-alpha_1 (1 atoms)
+colvars:   group1 (4 atoms)
+colvars:   Protein_C-alpha_2 (1 atoms)
+colvars:   group2 (4 atoms)
+colvars:   Protein_C-alpha_3 (1 atoms)
+colvars:   group3 (4 atoms)
+colvars:   Protein_C-alpha_4 (1 atoms)
+colvars:   group4 (4 atoms)
+colvars:   Protein_C-alpha_5 (1 atoms)
+colvars:   group5 (4 atoms)
+colvars:   Protein_C-alpha_6 (1 atoms)
+colvars:   group6 (4 atoms)
+colvars:   Protein_C-alpha_7 (1 atoms)
+colvars:   group7 (4 atoms)
+colvars:   Protein_C-alpha_8 (1 atoms)
+colvars:   group8 (4 atoms)
+colvars:   Protein_C-alpha_9 (1 atoms)
+colvars:   group9 (4 atoms)
+colvars:   Protein_C-alpha_10 (1 atoms)
+colvars:   group10 (4 atoms)
+colvars:   heavy_atoms (51 atoms)
+colvars: # smp = on [default]
+colvars: # colvarsTrajFrequency = 1
+colvars: # colvarsRestartFrequency = 10
+colvars: # scriptedColvarForces = off [default]
+colvars: # scriptingAfterBiases = off [default]
+colvars: # sourceTclFile = "" [default]
+colvars: ----------------------------------------------------------------------
+colvars:   Initializing a new collective variable.
+colvars:   # name = "one"
+colvars:   Initializing a new "distance" component.
+colvars:   # name = "" [default]
+colvars:   # componentCoeff = 1 [default]
+colvars:   # componentExp = 1 [default]
+colvars:   # period = 0 [default]
+colvars:   # wrapAround = 0 [default]
+colvars:   # forceNoPBC = off [default]
+colvars:   # scalable = on [default]
+colvars:     Initializing atom group "group1".
+colvars:     # name = "" [default]
+colvars:     # centerToOrigin = off [default]
+colvars:     # centerToReference = off [default]
+colvars:     # rotateToReference = off [default]
+colvars:     # atomsOfGroup = "" [default]
+colvars:     # indexGroup = "group1"
+colvars:     # psfSegID =  [default]
+colvars:     # atomsFile = "" [default]
+colvars:     # dummyAtom = ( 0 , 0 , 0 ) [default]
+colvars:     # enableFitGradients = on [default]
+colvars:     Enabling scalable calculation for group "group1".
+colvars:     # printAtomIDs = off [default]
+colvars:     Atom group "group1" defined with 4 atoms requested: total mass = 54.028, total charge = -0.72.
+colvars:     Initializing atom group "group2".
+colvars:     # name = "" [default]
+colvars:     # centerToOrigin = off [default]
+colvars:     # centerToReference = off [default]
+colvars:     # rotateToReference = off [default]
+colvars:     # atomsOfGroup = "" [default]
+colvars:     # indexGroup = "group2"
+colvars:     # psfSegID =  [default]
+colvars:     # atomsFile = "" [default]
+colvars:     # dummyAtom = ( 0 , 0 , 0 ) [default]
+colvars:     # enableFitGradients = on [default]
+colvars:     Enabling scalable calculation for group "group2".
+colvars:     # printAtomIDs = off [default]
+colvars:     Atom group "group2" defined with 4 atoms requested: total mass = 54.028, total charge = -0.4.
+colvars:   # oneSiteSystemForce = off [default]
+colvars:   # oneSiteTotalForce = off [default]
+colvars:   All components initialized.
+colvars:   # timeStepFactor = 1 [default]
+colvars:   # width = 0.5
+colvars:   # lowerBoundary = 3.2
+colvars:   # upperBoundary = 10.2
+colvars:   # hardLowerBoundary = on
+colvars:   # hardUpperBoundary = on
+colvars:   # expandBoundaries = off [default]
+colvars:   # extendedLagrangian = off [default]
+colvars:   # outputValue = on [default]
+colvars:   # outputVelocity = off [default]
+colvars:   # outputTotalForce = off [default]
+colvars:   # outputAppliedForce = on
+colvars:   # subtractAppliedForce = off [default]
+colvars:   # runAve = off [default]
+colvars:   # corrFunc = off [default]
+colvars: ----------------------------------------------------------------------
+colvars: Collective variables initialized, 1 in total.
+colvars: ----------------------------------------------------------------------
+colvars:   Initializing a new "harmonicwalls" instance.
+colvars:   # name = "wall_one"
+colvars:   # colvars = { one }
+colvars:   # stepZeroData = off [default]
+colvars:   # outputEnergy = off [default]
+colvars:   # outputFreq = 10 [default]
+colvars:   # timeStepFactor = 1 [default]
+colvars:   # writeTISamples = off [default]
+colvars:   # writeTIPMF = off [default]
+colvars:   # forceConstant = 1
+colvars:   # decoupling = off [default]
+colvars:   # targetForceConstant = -1 [default]
+colvars:   # lowerWalls = { 3.2 }
+colvars:   # upperWalls = { 10.2 }
+colvars:   # lowerWallConstant = 1 [default]
+colvars:   # upperWallConstant = 1 [default]
+colvars:   The lower wall force constant for colvar "one" will be rescaled to 4 according to the specified width (0.5).
+colvars:   The upper wall force constant for colvar "one" will be rescaled to 4 according to the specified width (0.5).
+colvars: ----------------------------------------------------------------------
+colvars:   Initializing a new "metadynamics" instance.
+colvars:   # name = "metadynamics1" [default]
+colvars:   # colvars = { one }
+colvars:   # stepZeroData = off [default]
+colvars:   # outputEnergy = off [default]
+colvars:   # outputFreq = 10 [default]
+colvars:   # timeStepFactor = 1 [default]
+colvars:   # writeTISamples = off [default]
+colvars:   # writeTIPMF = off [default]
+colvars:   # hillWeight = 0.001
+colvars:   # newHillFrequency = 10
+colvars:   # gaussianSigmas =  [default]
+colvars:   # hillWidth = 1.25331
+colvars:   Half-widths of the Gaussian hills (sigma's):
+colvars:   one: 0.313329
+colvars:   # multipleReplicas = off [default]
+colvars:   # useGrids = on [default]
+colvars:   # gridsUpdateFrequency = 10 [default]
+colvars:   # rebinGrids = off [default]
+colvars:   # writeFreeEnergyFile = on [default]
+colvars:   # keepHills = off [default]
+colvars:   # keepFreeEnergyFiles = off [default]
+colvars:   # writeHillsTrajectory = off [default]
+colvars:   # wellTempered = off [default]
+colvars:   # biasTemperature = -1 [default]
+colvars:   # useHillsReflection = on
+colvars:   # reflectionRange = 6 [default]
+colvars:   Reflection range is 6.
+colvars:   # reflectionLowLimitUseCVs =  [default]
+colvars:   Using all non-periodic variables for lower limits of reflection 
+colvars:   # reflectionUpLimitUseCVs =  [default]
+colvars:   Using all non-periodic variables for upper limits of reflection 
+colvars:   # reflectionLowLimit = { 3.2 } [default]
+colvars:   Reflection condition is applied on a lower limit for CV 0.
+colvars:   Reflection condition lower limit for this CV is 3.2.
+colvars:   # reflectionUpLimit = { 10.2 } [default]
+colvars:   Reflection condition is applied on an upper limit for CV 0.
+colvars:   Reflection condition upper limit for this CV is 10.2.
+colvars:   # ebMeta = off [default]
+colvars: ----------------------------------------------------------------------
+colvars: Collective variables biases initialized, 2 in total.
+colvars: ----------------------------------------------------------------------
+colvars: Collective variables module (re)initialized.
+colvars: ----------------------------------------------------------------------
+colvars: Updating NAMD interface:
+colvars: updating atomic data (0 atoms).
+colvars: updating group data (2 scalable groups, 8 atoms in total).
+colvars: updating grid object data (0 grid objects in total).
+colvars: 
+colvars: SUMMARY OF COLVARS FEATURES USED SO FAR AND THEIR CITATIONS:
+colvars: 
+colvars: - Colvars module:
+colvars: - Colvars-NAMD interface:
+colvars: - Metadynamics colvar bias implementation:
+colvars: - Optimal rotation via flexible fitting:
+colvars: - distance colvar component:
+colvars: - harmonicWalls colvar bias implementation:
+colvars:   Fiorin2013 https://doi.org/10.1080/00268976.2013.813594
+colvars: 
+colvars: - NAMD engine:
+colvars: - Scalable center-of-mass computation (NAMD):
+colvars:   Phillips2020 https://doi.org/10.1063/5.0014475
+colvars: 
+colvars: updating target temperature (T = 0 K).
+colvars: Updating NAMD interface:
+colvars: updating atomic data (0 atoms).
+colvars: updating group data (2 scalable groups, 8 atoms in total).
+colvars: updating grid object data (0 grid objects in total).
+colvars: updating target temperature (T = 0 K).
+colvars: Current simulation parameters: initial step = 0, integration timestep = 1
+colvars: Updating atomic parameters (masses, charges, etc).
+colvars: Re-initialized atom group for variable "one":0/0. 4 atoms: total mass = 54.028, total charge = -0.72.
+colvars: Re-initialized atom group for variable "one":0/1. 4 atoms: total mass = 54.028, total charge = -0.4.
+colvars: The restart output state file will be "test.tmp.colvars.state".
+colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj".
+colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj".
+colvars: Saving collective variables state to "test.tmp.colvars.state".
+colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj".
+colvars: Saving collective variables state to "test.tmp.colvars.state".
+colvars: Saving collective variables state to "test.colvars.state".
diff --git a/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.colvars.state.stripped b/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.colvars.state.stripped
new file mode 100644
index 000000000..ca7eeb6cc
--- /dev/null
+++ b/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.colvars.state.stripped
@@ -0,0 +1,52 @@
+configuration {
+  step           20
+  dt 1.000000e+00
+}
+
+colvar {
+  name one
+  x       3.2168853380692
+}
+
+restraint {
+  configuration {
+step 20
+name wall_one
+  }
+}
+
+metadynamics {
+  configuration {
+step 20
+name metadynamics1
+  }
+
+hills_energy
+grid_parameters {
+  n_colvars 1
+  lower_boundaries  3.2
+  upper_boundaries  10.2
+  widths  0.5
+  sizes  14
+}
+ 2.90793197471122e-03 2.29589040998682e-04 1.43098718458920e-06
+ 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00
+ 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00
+ 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00
+ 0.00000000000000e+00 0.00000000000000e+00
+
+hills_energy_gradients
+grid_parameters {
+  n_colvars 1
+  lower_boundaries  3.2
+  upper_boundaries  10.2
+  widths  0.5
+  sizes  14
+}
+ -7.38291182480299e-03 -1.74872845262777e-03 -1.81663860935017e-05
+ 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00
+ 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00
+ 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00
+ 0.00000000000000e+00 0.00000000000000e+00
+}
+
diff --git a/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.colvars.traj b/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.colvars.traj
new file mode 100644
index 000000000..84d326415
--- /dev/null
+++ b/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.colvars.traj
@@ -0,0 +1,22 @@
+# step         one                   fa_one                 
+           0    3.20554673468334e+00  0.00000000000000e+00  
+           1    3.20437148316484e+00  0.00000000000000e+00  
+           2    3.20384028588952e+00  0.00000000000000e+00  
+           3    3.20396721186862e+00  0.00000000000000e+00  
+           4    3.20472817286016e+00  0.00000000000000e+00  
+           5    3.20606308065050e+00  0.00000000000000e+00  
+           6    3.20787989000757e+00  0.00000000000000e+00  
+           7    3.21005997940810e+00  0.00000000000000e+00  
+           8    3.21246460651160e+00  0.00000000000000e+00  
+           9    3.21494247383722e+00  0.00000000000000e+00  
+          10    3.21733851230199e+00  3.69111079117202e-03  
+          11    3.21950373044117e+00  3.69111079117202e-03  
+          12    3.22130516783262e+00  3.69111079117202e-03  
+          13    3.22263513141286e+00  3.69111079117202e-03  
+          14    3.22341825276891e+00  3.69111079117202e-03  
+          15    3.22361556235479e+00  3.69111079117202e-03  
+          16    3.22322531690825e+00  3.69111079117202e-03  
+          17    3.22228092687886e+00  3.69111079117202e-03  
+          18    3.22084669842434e+00  3.69111079117202e-03  
+          19    3.21901213410055e+00  3.69111079117202e-03  
+          20    3.21688533806922e+00  7.38291182480299e-03  
diff --git a/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.pmf b/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.pmf
new file mode 100644
index 000000000..7aa384a05
--- /dev/null
+++ b/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.pmf
@@ -0,0 +1,17 @@
+# 1
+#  3.20000000000000e+00  5.00000000000000e-01         14  0
+
+  3.45000000000000e+00  -0.00000000000000e+00
+  3.95000000000000e+00   2.67834293371254e-03
+  4.45000000000000e+00   2.90650098752663e-03
+  4.95000000000000e+00   2.90793197471122e-03
+  5.45000000000000e+00   2.90793197471122e-03
+  5.95000000000000e+00   2.90793197471122e-03
+  6.45000000000000e+00   2.90793197471122e-03
+  6.95000000000000e+00   2.90793197471122e-03
+  7.45000000000000e+00   2.90793197471122e-03
+  7.95000000000000e+00   2.90793197471122e-03
+  8.45000000000000e+00   2.90793197471122e-03
+  8.95000000000000e+00   2.90793197471122e-03
+  9.45000000000000e+00   2.90793197471122e-03
+  9.95000000000000e+00   2.90793197471122e-03
diff --git a/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.restart.colvars.out b/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.restart.colvars.out
new file mode 100644
index 000000000..056ab6f7e
--- /dev/null
+++ b/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.restart.colvars.out
@@ -0,0 +1,215 @@
+colvars: ----------------------------------------------------------------------
+colvars: Please cite Fiorin et al, Mol Phys 2013:
+colvars:   https://doi.org/10.1080/00268976.2013.813594
+colvars: as well as all other papers listed below for individual features used.
+colvars: SMP parallelism is enabled (num threads = 4).
+colvars: This version was built with the C++11 standard or higher.
+colvars: Redefining the Tcl "cv" command to the new script interface.
+colvars: The restart output state file will be "test.restart.tmp.colvars.state".
+colvars: The final output state file will be "test.restart.colvars.state".
+colvars: ----------------------------------------------------------------------
+colvars: Reading new configuration from file "test.in":
+colvars: # units = "" [default]
+colvars: # indexFile = "index.ndx"
+colvars: The following index groups are currently defined:
+colvars:   Protein (104 atoms)
+colvars:   Protein_noH (51 atoms)
+colvars:   Protein_Backbone (40 atoms)
+colvars:   Protein_C-alpha (10 atoms)
+colvars:   RMSD_atoms (10 atoms)
+colvars:   Protein_C-alpha_1_2 (2 atoms)
+colvars:   Protein_C-alpha_9_10 (2 atoms)
+colvars:   Protein_C-alpha_1 (1 atoms)
+colvars:   group1 (4 atoms)
+colvars:   Protein_C-alpha_2 (1 atoms)
+colvars:   group2 (4 atoms)
+colvars:   Protein_C-alpha_3 (1 atoms)
+colvars:   group3 (4 atoms)
+colvars:   Protein_C-alpha_4 (1 atoms)
+colvars:   group4 (4 atoms)
+colvars:   Protein_C-alpha_5 (1 atoms)
+colvars:   group5 (4 atoms)
+colvars:   Protein_C-alpha_6 (1 atoms)
+colvars:   group6 (4 atoms)
+colvars:   Protein_C-alpha_7 (1 atoms)
+colvars:   group7 (4 atoms)
+colvars:   Protein_C-alpha_8 (1 atoms)
+colvars:   group8 (4 atoms)
+colvars:   Protein_C-alpha_9 (1 atoms)
+colvars:   group9 (4 atoms)
+colvars:   Protein_C-alpha_10 (1 atoms)
+colvars:   group10 (4 atoms)
+colvars:   heavy_atoms (51 atoms)
+colvars: # smp = on [default]
+colvars: # colvarsTrajFrequency = 1
+colvars: # colvarsRestartFrequency = 10
+colvars: # scriptedColvarForces = off [default]
+colvars: # scriptingAfterBiases = off [default]
+colvars: # sourceTclFile = "" [default]
+colvars: ----------------------------------------------------------------------
+colvars:   Initializing a new collective variable.
+colvars:   # name = "one"
+colvars:   Initializing a new "distance" component.
+colvars:   # name = "" [default]
+colvars:   # componentCoeff = 1 [default]
+colvars:   # componentExp = 1 [default]
+colvars:   # period = 0 [default]
+colvars:   # wrapAround = 0 [default]
+colvars:   # forceNoPBC = off [default]
+colvars:   # scalable = on [default]
+colvars:     Initializing atom group "group1".
+colvars:     # name = "" [default]
+colvars:     # centerToOrigin = off [default]
+colvars:     # centerToReference = off [default]
+colvars:     # rotateToReference = off [default]
+colvars:     # atomsOfGroup = "" [default]
+colvars:     # indexGroup = "group1"
+colvars:     # psfSegID =  [default]
+colvars:     # atomsFile = "" [default]
+colvars:     # dummyAtom = ( 0 , 0 , 0 ) [default]
+colvars:     # enableFitGradients = on [default]
+colvars:     Enabling scalable calculation for group "group1".
+colvars:     # printAtomIDs = off [default]
+colvars:     Atom group "group1" defined with 4 atoms requested: total mass = 54.028, total charge = -0.72.
+colvars:     Initializing atom group "group2".
+colvars:     # name = "" [default]
+colvars:     # centerToOrigin = off [default]
+colvars:     # centerToReference = off [default]
+colvars:     # rotateToReference = off [default]
+colvars:     # atomsOfGroup = "" [default]
+colvars:     # indexGroup = "group2"
+colvars:     # psfSegID =  [default]
+colvars:     # atomsFile = "" [default]
+colvars:     # dummyAtom = ( 0 , 0 , 0 ) [default]
+colvars:     # enableFitGradients = on [default]
+colvars:     Enabling scalable calculation for group "group2".
+colvars:     # printAtomIDs = off [default]
+colvars:     Atom group "group2" defined with 4 atoms requested: total mass = 54.028, total charge = -0.4.
+colvars:   # oneSiteSystemForce = off [default]
+colvars:   # oneSiteTotalForce = off [default]
+colvars:   All components initialized.
+colvars:   # timeStepFactor = 1 [default]
+colvars:   # width = 0.5
+colvars:   # lowerBoundary = 3.2
+colvars:   # upperBoundary = 10.2
+colvars:   # hardLowerBoundary = on
+colvars:   # hardUpperBoundary = on
+colvars:   # expandBoundaries = off [default]
+colvars:   # extendedLagrangian = off [default]
+colvars:   # outputValue = on [default]
+colvars:   # outputVelocity = off [default]
+colvars:   # outputTotalForce = off [default]
+colvars:   # outputAppliedForce = on
+colvars:   # subtractAppliedForce = off [default]
+colvars:   # runAve = off [default]
+colvars:   # corrFunc = off [default]
+colvars: ----------------------------------------------------------------------
+colvars: Collective variables initialized, 1 in total.
+colvars: ----------------------------------------------------------------------
+colvars:   Initializing a new "harmonicwalls" instance.
+colvars:   # name = "wall_one"
+colvars:   # colvars = { one }
+colvars:   # stepZeroData = off [default]
+colvars:   # outputEnergy = off [default]
+colvars:   # outputFreq = 10 [default]
+colvars:   # timeStepFactor = 1 [default]
+colvars:   # writeTISamples = off [default]
+colvars:   # writeTIPMF = off [default]
+colvars:   # forceConstant = 1
+colvars:   # decoupling = off [default]
+colvars:   # targetForceConstant = -1 [default]
+colvars:   # lowerWalls = { 3.2 }
+colvars:   # upperWalls = { 10.2 }
+colvars:   # lowerWallConstant = 1 [default]
+colvars:   # upperWallConstant = 1 [default]
+colvars:   The lower wall force constant for colvar "one" will be rescaled to 4 according to the specified width (0.5).
+colvars:   The upper wall force constant for colvar "one" will be rescaled to 4 according to the specified width (0.5).
+colvars: ----------------------------------------------------------------------
+colvars:   Initializing a new "metadynamics" instance.
+colvars:   # name = "metadynamics1" [default]
+colvars:   # colvars = { one }
+colvars:   # stepZeroData = off [default]
+colvars:   # outputEnergy = off [default]
+colvars:   # outputFreq = 10 [default]
+colvars:   # timeStepFactor = 1 [default]
+colvars:   # writeTISamples = off [default]
+colvars:   # writeTIPMF = off [default]
+colvars:   # hillWeight = 0.001
+colvars:   # newHillFrequency = 10
+colvars:   # gaussianSigmas =  [default]
+colvars:   # hillWidth = 1.25331
+colvars:   Half-widths of the Gaussian hills (sigma's):
+colvars:   one: 0.313329
+colvars:   # multipleReplicas = off [default]
+colvars:   # useGrids = on [default]
+colvars:   # gridsUpdateFrequency = 10 [default]
+colvars:   # rebinGrids = off [default]
+colvars:   # writeFreeEnergyFile = on [default]
+colvars:   # keepHills = off [default]
+colvars:   # keepFreeEnergyFiles = off [default]
+colvars:   # writeHillsTrajectory = off [default]
+colvars:   # wellTempered = off [default]
+colvars:   # biasTemperature = -1 [default]
+colvars:   # useHillsReflection = on
+colvars:   # reflectionRange = 6 [default]
+colvars:   Reflection range is 6.
+colvars:   # reflectionLowLimitUseCVs =  [default]
+colvars:   Using all non-periodic variables for lower limits of reflection 
+colvars:   # reflectionUpLimitUseCVs =  [default]
+colvars:   Using all non-periodic variables for upper limits of reflection 
+colvars:   # reflectionLowLimit = { 3.2 } [default]
+colvars:   Reflection condition is applied on a lower limit for CV 0.
+colvars:   Reflection condition lower limit for this CV is 3.2.
+colvars:   # reflectionUpLimit = { 10.2 } [default]
+colvars:   Reflection condition is applied on an upper limit for CV 0.
+colvars:   Reflection condition upper limit for this CV is 10.2.
+colvars:   # ebMeta = off [default]
+colvars: ----------------------------------------------------------------------
+colvars: Collective variables biases initialized, 2 in total.
+colvars: ----------------------------------------------------------------------
+colvars: Collective variables module (re)initialized.
+colvars: ----------------------------------------------------------------------
+colvars: Updating NAMD interface:
+colvars: updating atomic data (0 atoms).
+colvars: updating group data (2 scalable groups, 8 atoms in total).
+colvars: updating grid object data (0 grid objects in total).
+colvars: 
+colvars: SUMMARY OF COLVARS FEATURES USED SO FAR AND THEIR CITATIONS:
+colvars: 
+colvars: - Colvars module:
+colvars: - Colvars-NAMD interface:
+colvars: - Metadynamics colvar bias implementation:
+colvars: - Optimal rotation via flexible fitting:
+colvars: - distance colvar component:
+colvars: - harmonicWalls colvar bias implementation:
+colvars:   Fiorin2013 https://doi.org/10.1080/00268976.2013.813594
+colvars: 
+colvars: - NAMD engine:
+colvars: - Scalable center-of-mass computation (NAMD):
+colvars:   Phillips2020 https://doi.org/10.1063/5.0014475
+colvars: 
+colvars: updating target temperature (T = 0 K).
+colvars: ----------------------------------------------------------------------
+colvars: Loading state from text file "test.colvars.state".
+colvars:   Restarting collective variable "one" from value: 3.21689
+colvars:   Restarted harmonicwalls bias "wall_one" with step number 20.
+colvars:     successfully read the biasing potential and its gradients from grids.
+colvars:     successfully read 0 explicit hills from state.
+colvars:   Restarted metadynamics bias "metadynamics1" with step number 20.
+colvars: ----------------------------------------------------------------------
+colvars: Updating NAMD interface:
+colvars: updating atomic data (0 atoms).
+colvars: updating group data (2 scalable groups, 8 atoms in total).
+colvars: updating grid object data (0 grid objects in total).
+colvars: updating target temperature (T = 0 K).
+colvars: Current simulation parameters: initial step = 20, integration timestep = 1
+colvars: Updating atomic parameters (masses, charges, etc).
+colvars: Re-initialized atom group for variable "one":0/0. 4 atoms: total mass = 54.028, total charge = -0.72.
+colvars: Re-initialized atom group for variable "one":0/1. 4 atoms: total mass = 54.028, total charge = -0.4.
+colvars: The restart output state file will be "test.restart.tmp.colvars.state".
+colvars: Synchronizing (emptying the buffer of) trajectory file "test.restart.colvars.traj".
+colvars: Synchronizing (emptying the buffer of) trajectory file "test.restart.colvars.traj".
+colvars: Saving collective variables state to "test.restart.tmp.colvars.state".
+colvars: Synchronizing (emptying the buffer of) trajectory file "test.restart.colvars.traj".
+colvars: Saving collective variables state to "test.restart.tmp.colvars.state".
+colvars: Saving collective variables state to "test.restart.colvars.state".
diff --git a/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.restart.colvars.state.stripped b/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.restart.colvars.state.stripped
new file mode 100644
index 000000000..54fade905
--- /dev/null
+++ b/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.restart.colvars.state.stripped
@@ -0,0 +1,52 @@
+configuration {
+  step           40
+  dt 1.000000e+00
+}
+
+colvar {
+  name one
+  x         3.21407617475
+}
+
+restraint {
+  configuration {
+step 40
+name wall_one
+  }
+}
+
+metadynamics {
+  configuration {
+step 40
+name metadynamics1
+  }
+
+hills_energy
+grid_parameters {
+  n_colvars 1
+  lower_boundaries  3.2
+  upper_boundaries  10.2
+  widths  0.5
+  sizes  14
+}
+ 5.81689037532263e-03 4.58131531298534e-04 2.84164743658542e-06
+ 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00
+ 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00
+ 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00
+ 0.00000000000000e+00 0.00000000000000e+00
+
+hills_energy_gradients
+grid_parameters {
+  n_colvars 1
+  lower_boundaries  3.2
+  upper_boundaries  10.2
+  widths  0.5
+  sizes  14
+}
+ -1.47827984492445e-02 -3.49284785612392e-03 -3.61088189994010e-05
+ 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00
+ 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00
+ 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00
+ 0.00000000000000e+00 0.00000000000000e+00
+}
+
diff --git a/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.restart.colvars.traj b/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.restart.colvars.traj
new file mode 100644
index 000000000..c0034b7fa
--- /dev/null
+++ b/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.restart.colvars.traj
@@ -0,0 +1,22 @@
+# step         one                   fa_one                 
+          20    3.21688533806922e+00  7.38291182480299e-03  
+          21    3.21458588882441e+00  7.38291182480299e-03  
+          22    3.21223714587457e+00  7.38291182480299e-03  
+          23    3.20995862787481e+00  7.38291182480299e-03  
+          24    3.20785876515602e+00  7.38291182480299e-03  
+          25    3.20602880105322e+00  7.38291182480299e-03  
+          26    3.20453847572721e+00  7.38291182480299e-03  
+          27    3.20343382065677e+00  7.38291182480299e-03  
+          28    3.20273699133106e+00  7.38291182480299e-03  
+          29    3.20244768843368e+00  7.38291182480299e-03  
+          30    3.20254574820329e+00  1.10871252733508e-02  
+          31    3.20299458251580e+00  1.10871252733508e-02  
+          32    3.20374508637636e+00  1.10871252733508e-02  
+          33    3.20474018224177e+00  1.10871252733508e-02  
+          34    3.20591955787466e+00  1.10871252733508e-02  
+          35    3.20722423546519e+00  1.10871252733508e-02  
+          36    3.20860052667966e+00  1.10871252733508e-02  
+          37    3.21000305295022e+00  1.10871252733508e-02  
+          38    3.21139678731726e+00  1.10871252733508e-02  
+          39    3.21275826748767e+00  1.10871252733508e-02  
+          40    3.21407617475000e+00  1.47827984492445e-02  
diff --git a/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.restart.pmf b/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.restart.pmf
new file mode 100644
index 000000000..a360859b5
--- /dev/null
+++ b/namd/tests/library/000_distance-grid_metadynamics-reflection/AutoDiff/test.restart.pmf
@@ -0,0 +1,17 @@
+# 1
+#  3.20000000000000e+00  5.00000000000000e-01         14  0
+
+  3.45000000000000e+00  -0.00000000000000e+00
+  3.95000000000000e+00   5.35875884402410e-03
+  4.45000000000000e+00   5.81404872788604e-03
+  4.95000000000000e+00   5.81689037532263e-03
+  5.45000000000000e+00   5.81689037532263e-03
+  5.95000000000000e+00   5.81689037532263e-03
+  6.45000000000000e+00   5.81689037532263e-03
+  6.95000000000000e+00   5.81689037532263e-03
+  7.45000000000000e+00   5.81689037532263e-03
+  7.95000000000000e+00   5.81689037532263e-03
+  8.45000000000000e+00   5.81689037532263e-03
+  8.95000000000000e+00   5.81689037532263e-03
+  9.45000000000000e+00   5.81689037532263e-03
+  9.95000000000000e+00   5.81689037532263e-03
diff --git a/namd/tests/library/000_distance-grid_metadynamics-reflection/namd-version.txt b/namd/tests/library/000_distance-grid_metadynamics-reflection/namd-version.txt
new file mode 100644
index 000000000..fd296a0f7
--- /dev/null
+++ b/namd/tests/library/000_distance-grid_metadynamics-reflection/namd-version.txt
@@ -0,0 +1,3 @@
+Info: NAMD 3.0.1 for Linux-x86_64-multicore
+colvars: Initializing the collective variables module, version 2023-12-05.
+colvars: Using NAMD interface, version "2023-12-05".
diff --git a/namd/tests/library/000_distance-grid_metadynamics-reflection/test.in b/namd/tests/library/000_distance-grid_metadynamics-reflection/test.in
new file mode 100644
index 000000000..97b3e5b99
--- /dev/null
+++ b/namd/tests/library/000_distance-grid_metadynamics-reflection/test.in
@@ -0,0 +1,46 @@
+colvarsTrajFrequency 1
+colvarsRestartFrequency 10
+indexFile index.ndx
+
+colvar {
+
+    name one
+
+    outputAppliedForce on
+
+    # use a non-trivial width to test bias behavior
+    width 0.5
+    # The lower boundary is already defined at 0 for a distance function
+    lowerBoundary 3.2
+    upperBoundary 10.2
+
+    # not required to define hard limits but recommended to avoid unnecessary grid expansion
+    # likely to be activated by default with refletion in the future
+    hardUpperBoundary on
+    hardLowerBoundary on
+
+    distance {
+        group1 {
+            indexGroup group1
+        }
+        group2 {
+            indexGroup group2
+        }
+    }
+} 
+
+metadynamics {
+    colvars        one
+    hillWeight     0.001
+    hillWidth      1.2533141373155001  # Old default
+    newHillFrequency 10
+    useHillsReflection on
+}
+
+harmonicWalls {
+    name wall_one
+    colvars one
+    forceConstant 1.0
+    lowerWalls 3.2
+    upperWalls 10.2
+}
diff --git a/src/colvarbias_meta.cpp b/src/colvarbias_meta.cpp
index 1131c88ec..a0c92b2ce 100644
--- a/src/colvarbias_meta.cpp
+++ b/src/colvarbias_meta.cpp
@@ -156,6 +156,7 @@ int colvarbias_meta::init(std::string const &conf)
 
   error_code |= init_replicas_params(conf);
   error_code |= init_well_tempered_params(conf);
+  error_code |= init_reflection_params(conf);
   error_code |= init_ebmeta_params(conf);
 
   if (cvm::debug())
@@ -311,6 +312,304 @@ int colvarbias_meta::init_ebmeta_params(std::string const &conf)
   return error_code;
 }
 
+int colvarbias_meta::init_reflection_params(std::string const &conf)
+{
+  nrefvarsl=0;
+  nrefvarsu=0;
+  size_t nonpvars;
+  size_t i;
+  size_t icount=0;
+  size_t j;
+  size_t jcount;
+
+  // in future remove the line below and uncomment the three following ones as reflection should be default with grids
+
+  use_reflection=false;
+
+  //if (use_grids) {
+  //  use_reflection=true;
+  //}
+
+  for ( i = 0; i < num_variables(); i++ ) {
+     if (!variables(i)->is_enabled(f_cv_periodic)) {
+       icount++;
+     }
+  }
+  nonpvars=icount;
+
+  get_keyval(conf, "useHillsReflection", use_reflection, use_reflection);
+
+  for ( i = 0; i < num_variables(); i++ ) {
+       if (variables(i)->value().type()!=colvarvalue::type_scalar) {
+	 use_reflection=false;
+         cvm::log("Note: CV "+cvm::to_str(variables(i)->name)+" is not of scalar type. Hills reflection has been disabled as it can be used only with scalar variables.\n");
+       }
+  }
+
+  if (use_reflection) {
+
+    get_keyval(conf, "reflectionRange", reflection_int, 6.0);
+    cvm::log("Reflection range is "+cvm::to_str(reflection_int)+".\n");
+
+    std::vector<std::string> colvar_names_l;
+    std::vector<int> check_dup(num_variables());
+    if (get_keyval(conf, "reflectionLowLimitUseCVs", colvar_names_l)) {
+      nrefvarsl=colvar_names_l.size();
+      reflection_llimit_cv.resize(nrefvarsl);
+
+      for (i = 0; i < num_variables(); i++) check_dup[i]=0;
+      for (j = 0; j < nrefvarsl; j++) {
+         icount=0;
+         for (i = 0; i < num_variables(); i++) {
+            if (colvar_names_l[j].compare(variables(i)->name)==0) {
+              reflection_llimit_cv[j]=i;
+              icount++;
+              check_dup[i]++;
+            }
+         }
+         if (icount!=1) cvm::error("Error: CV name "+colvar_names_l[j]+" mismatch or duplication for lower reflection limit \n", COLVARS_INPUT_ERROR);
+      }
+      for (i = 0; i < num_variables(); i++) {
+         if (check_dup[i]>1) cvm::error("Error: CV name duplication for lower reflection limit \n", COLVARS_INPUT_ERROR);
+      }
+      if(nrefvarsl>num_variables()) cvm::error("Error: number of CVs with active lower reflection limit is > num_variables  \n", COLVARS_INPUT_ERROR);
+      cvm::log("Using lower limits reflection on "+cvm::to_str(nrefvarsl)+" variables.\n");
+    } else {
+      nrefvarsl=nonpvars;
+      reflection_llimit_cv.resize(nrefvarsl);
+      icount=0;
+      for (i = 0; i < num_variables(); i++) {
+         if (!variables(i)->is_enabled(f_cv_periodic)) {
+           reflection_llimit_cv[icount]=i;
+           icount++;
+         }
+      }
+      cvm::log("Using all non-periodic variables for lower limits of reflection \n");
+    }
+
+    if (reflection_llimit.size()==0) {
+      reflection_llimit.resize(nrefvarsl);
+    }
+
+    std::vector<std::string> colvar_names_u;
+    if (get_keyval(conf, "reflectionUpLimitUseCVs", colvar_names_u)) {
+      nrefvarsu=colvar_names_u.size();
+      reflection_ulimit_cv.resize(nrefvarsu);
+      for (i = 0; i < num_variables(); i++) check_dup[i]=0;
+      for (j = 0; j < nrefvarsu; j++) {
+         icount=0;
+         for (i = 0; i < num_variables(); i++) {
+            if (colvar_names_u[j].compare(variables(i)->name)==0) {
+              reflection_ulimit_cv[j]=i;
+              icount++;
+              check_dup[i]++;
+            }
+         }
+         if (icount!=1) cvm::error("Error: CV name "+cvm::to_str(colvar_names_u[j])+" mismatch or duplication for upper reflection limit \n", COLVARS_INPUT_ERROR);
+      }
+      for (i = 0; i < num_variables(); i++) {
+         if (check_dup[i]>1) cvm::error("Error: CV name duplication for upper reflection limit \n", COLVARS_INPUT_ERROR);
+      }
+      if(nrefvarsu>num_variables()) cvm::error("Error: number of CVs with active upper reflection limit is > num_variables  \n", COLVARS_INPUT_ERROR);
+      cvm::log("Using upper limits reflection on "+cvm::to_str(nrefvarsu)+" variables.\n");
+    } else {
+      nrefvarsu=nonpvars;
+      reflection_ulimit_cv.resize(nrefvarsu);
+      icount=0;
+      for (i = 0; i < num_variables(); i++) {
+         if (!variables(i)->is_enabled(f_cv_periodic)) {
+           reflection_ulimit_cv[icount]=i;
+           icount++;
+         }
+      }
+      cvm::log("Using all non-periodic variables for upper limits of reflection \n");
+    }
+
+    if (reflection_ulimit.size()==0) {
+      reflection_ulimit.resize(nrefvarsu);
+    }
+
+    // use reflection only with scalar variables
+
+    for (i = 0; i < nrefvarsl; i++) {
+       if (reflection_llimit_cv[i]>=static_cast<int>(num_variables()) || reflection_llimit_cv[i]<0) {
+         cvm::error("Error: CV number is negative or >= num_variables  \n", COLVARS_INPUT_ERROR);
+       }
+       j=reflection_llimit_cv[i];
+       if (variables(j)->is_enabled(f_cv_periodic)) {
+         cvm::log("Warning: you are using hills reflection with a periodic variable, make sure you are using it far from periodic boundaries \n");
+       }
+    }
+
+    for (i = 0; i < nrefvarsu; i++) {
+       if (reflection_ulimit_cv[i]>=static_cast<int>(num_variables()) || reflection_ulimit_cv[i]<0) {
+         cvm::error("Error: CV number is negative or >= num_variables  \n", COLVARS_INPUT_ERROR);
+       }
+       j=reflection_ulimit_cv[i];
+       if (variables(j)->is_enabled(f_cv_periodic)) {
+         cvm::log("Warning: you are using hills reflection with a periodic variable, make sure you are using it far from periodic boundaries \n");
+       }
+    }
+
+    // if grids are defined, set by default reflection boundaries as grid boundaries
+    if (use_grids) {
+      for (i = 0; i < nrefvarsl; i++) {
+         icount=reflection_llimit_cv[i];
+         reflection_llimit[i]=hills_energy->lower_boundaries[icount].real_value;
+      }
+      for (i = 0; i < nrefvarsu; i++) {
+         icount=reflection_ulimit_cv[i];
+         reflection_ulimit[i]=hills_energy->upper_boundaries[icount].real_value;
+      }
+    }
+
+    if (nrefvarsl>0) {
+      if (get_keyval(conf, "reflectionLowLimit", reflection_llimit, reflection_llimit)) {
+        for (i = 0; i < nrefvarsl; i++) {
+           if (use_grids) {
+             icount=reflection_llimit_cv[i];
+             cvm:: real bound=hills_energy->lower_boundaries[icount].real_value;
+             if (reflection_llimit[i] != bound && reflection_llimit[i] < bound+variables(icount)->width) {
+               cvm::error("Error: please set lower reflection limit for CV "+cvm::to_str(icount)+" either at grid lower boundary ("+cvm::to_str(bound)+") or well above it (above "+cvm::to_str(bound+variables(icount)->width)+").\n", COLVARS_INPUT_ERROR);
+             }
+           }
+        }
+      } else {
+        if (!use_grids) {
+          cvm::error("Error: Lower limits for reflection not provided.\n", COLVARS_INPUT_ERROR);
+        }
+      }
+    }
+
+    for (i = 0; i < nrefvarsl; i++) {
+       cvm::log("Reflection condition is applied on a lower limit for CV "+cvm::to_str(variables(reflection_llimit_cv[i])->name)+".\n");
+       cvm::log("Reflection condition lower limit for this CV is "+cvm::to_str(reflection_llimit[i])+".\n");
+    }
+
+    if (nrefvarsu>0) {
+      if (get_keyval(conf, "reflectionUpLimit", reflection_ulimit, reflection_ulimit)) {
+        for (i = 0; i < nrefvarsu; i++) {
+           if (use_grids) {
+             icount=reflection_ulimit_cv[i];
+             cvm:: real bound=hills_energy->upper_boundaries[icount].real_value;
+             if (reflection_ulimit[i] != bound && reflection_ulimit[i] > bound-variables(icount)->width) {
+               cvm::error("Error: please set upper reflection limit for CV "+cvm::to_str(icount)+" either at grid upper boundary ("+cvm::to_str(bound)+") or well below it (below "+cvm::to_str(bound-variables(icount)->width)+").\n", COLVARS_INPUT_ERROR);
+             }
+           }
+        }
+      } else {
+        if (!use_grids) {
+          cvm::error("Error: Upper limits for reflection not provided.\n", COLVARS_INPUT_ERROR);
+        }
+      }
+    }
+
+    for (i = 0; i < nrefvarsu; i++) {
+       cvm::log("Reflection condition is applied on an upper limit for CV "+cvm::to_str(variables(reflection_ulimit_cv[i])->name)+".\n");
+       cvm::log("Reflection condition upper limit for this CV is "+cvm::to_str(reflection_ulimit[i])+".\n");
+    }
+
+    // multimensional reflection
+
+    // generate reflection states
+    size_t sum;
+    sum=1;
+    size_t nstates;
+    size_t count;
+
+    if (reflection_usel.size()==0) {
+      reflection_usel.resize(num_variables(),std::vector<bool>(2));
+    }
+
+    if (reflection_l.size()==0) {
+      reflection_l.resize(num_variables(),std::vector<cvm::real>(2));
+    }
+
+    for (j = 1; j < num_variables(); j++) {
+       reflection_usel[j][0]=false;
+       reflection_l[j][0]=0.0;
+       reflection_usel[j][1]=false;
+       reflection_l[j][1]=0.0;
+    }
+
+    for (i = 0; i < nrefvarsl; i++) {
+       j=reflection_llimit_cv[i];
+       reflection_usel[j][0]=true;
+       reflection_l[j][0]=reflection_llimit[i];
+    }
+
+    for (i = 0; i < nrefvarsu; i++) {
+       j=reflection_ulimit_cv[i];
+       reflection_usel[j][1]=true;
+       reflection_l[j][1]=reflection_ulimit[i];
+    }
+
+//  Generate all possible reflection states (e.g. through faces, edges and vertex).
+//  Consider for example a cube, the states are:
+//  [0,0,1]
+//  [0,1,0] [0,1,1]
+//  [1,0,0] [1,0,1] [1,1,0] [1,1,1]
+//  where 1 means reflect on that coordinate and 0 do not reflect.
+//  These states can be generated as:
+//  ref_state[0][0]=1
+//  ref_state[1][0]=10  ref_state[1][1]=11
+//  ref_state[2][0]=100 ref_state[2][1]=101 ref_state[2][2]=110 ref_state[2][3]=111
+//  going down along the rows the size ref_state[j].size() is the number of previous states
+//  (j-1) plus one.
+//  A specific state instead can be generated starting from a power of 10 and then summing
+//  the states of the previous rows:
+//  ref_state[1][1]=ref_state[1][0]+ref_state[0][0]
+//  ref_state[2][1]=ref_state[2][0]+ref_state[0][0]
+//  ref_state[2][2]=ref_state[2][0]+ref_state[1][0]
+//  ref_state[2][3]=ref_state[2][0]+ref_state[1][1]
+
+    if (ref_state.size()==0) {
+      ref_state.resize(num_variables(),std::vector<size_t>(1));
+    }
+    ref_state[0][0]=1;
+    for (j = 1; j < num_variables(); j++) {
+      sum*=10;
+      nstates=0;
+      for (jcount = 0; jcount < j; jcount++) {
+            nstates+=ref_state[jcount].size();
+      }
+      nstates++;
+      ref_state[j].resize(nstates);
+      ref_state[j][0]=sum;
+      count=0;
+      for (jcount = 0; jcount < j; jcount++) {
+         for ( icount = 0; icount < ref_state[jcount].size(); icount++) {
+            count++;
+            ref_state[j][count]=ref_state[j][0]+ref_state[jcount][icount];
+         }
+      }
+    }
+  }
+
+  if (which_int_llimit_cv.size()==0) {
+    which_int_llimit_cv.resize(num_variables());
+  }
+  for ( i = 0; i < num_variables(); i++) {
+     which_int_llimit_cv[i]=-1;
+  }
+  for ( i = 0; i < nrefvarsl; i++) {
+     j=reflection_llimit_cv[i];
+     which_int_llimit_cv[j]=i;
+  }
+
+  if (which_int_ulimit_cv.size()==0) {
+    which_int_ulimit_cv.resize(num_variables());
+  }
+  for ( i = 0; i < num_variables(); i++) {
+     which_int_ulimit_cv[i]=-1;
+  }
+  for ( i = 0; i < nrefvarsu; i++) {
+     j=reflection_ulimit_cv[i];
+     which_int_ulimit_cv[j]=i;
+  }
+
+  return COLVARS_OK;
+}
 
 colvarbias_meta::~colvarbias_meta()
 {
@@ -365,6 +664,167 @@ colvarbias_meta::add_hill(colvarbias_meta::hill const &h)
   return hills.end();
 }
 
+bool colvarbias_meta::check_reflection_limits(bool &ah)
+{
+  size_t i;
+  size_t icount;
+  for ( i = 0; i < nrefvarsl; i++) {
+     icount=reflection_llimit_cv[i];
+     if (colvar_values[icount]<reflection_llimit[i]) {
+       ah=false;
+     }
+  }
+  for ( i = 0; i < nrefvarsu; i++) {
+     icount=reflection_ulimit_cv[i];
+     if (colvar_values[icount]>reflection_ulimit[i]) {
+       ah=false;
+     }
+  }
+  return ah;
+}
+
+int colvarbias_meta::reflect_hill_multid(cvm::real const &h_scale)
+{
+  size_t i = 0;
+  size_t j;
+  size_t jcount;
+  size_t startsum;
+  int getsum;
+  int check_val;
+  size_t numberref;
+  size_t startsumk;
+  int upordown;
+  size_t nkstates;
+  size_t kstate;
+  size_t k;
+  size_t kcount;
+  size_t countstate;
+  bool hill_add;
+  int getsumk;
+  int checkstate;
+  size_t state;
+  int valk;
+  cvm:: real tmps;
+  colvarvalue tmp;
+  colvarvalue unitary;
+  cvm:: real reflection_limit;
+  cvm:: real tmpd;
+
+  std::vector<colvarvalue> curr_cv_values(num_variables());
+  for (i = 0; i < num_variables(); i++) {
+    curr_cv_values[i].type(variables(i)->value());
+  }
+  for (i = 0; i < num_variables(); i++) {
+      curr_cv_values[i] = colvar_values[i];
+  }
+
+  // sum over all possible reflection states previously generated,
+  // see init
+
+  for ( j = 0; j < num_variables(); j++) {
+     startsum=1;
+     for (i = 0; i < j; i++) {
+        startsum*=10;
+     }
+     for (jcount = 0; jcount < ref_state[j].size(); jcount++) {
+           getsum=startsum;
+           check_val=ref_state[j][jcount];
+           numberref=0;
+           startsumk=1;
+           for (i = 0; i <= j; i++) {
+              upordown=std::floor(check_val/getsum);
+              if(check_val-getsum>=0) check_val=check_val-getsum;
+              getsum=getsum/10;
+              if (upordown==1) {
+                numberref++;
+                if(numberref>1) startsumk*=10;
+              }
+           }
+
+           // sum over all possible lower and upper boudary combinations
+           // exploiting kstate=ref_state[k][kcount]:
+           // for just one reflection these are 0(lower boundary) and 1(upper boundary)
+           // for two reflections these are 0 1 10 11 (0,0 0,1 1,0 1,1)
+           // where 0 is reflect on the two lower boudaries of the two coordinates etc.
+
+           nkstates=2;
+           kstate=0;
+           for ( k = 0; k <numberref; k++ ) {
+              if (k>0)  nkstates=ref_state[k].size();
+              for ( kcount = 0; kcount < nkstates; kcount++) {
+                 if (k==0 && kcount==1) {
+                   kstate=1;
+                 } else if (k>0) {
+                   kstate=ref_state[k][kcount];
+                 }
+
+                 getsum=startsum;
+                 countstate=0;
+                 check_val=ref_state[j][jcount];
+                 hill_add=true;
+                 getsumk=startsumk;
+                 checkstate=kstate;
+                 for (i = 0; i <= j; i++) {
+                    upordown=std::floor(check_val/getsum);
+                    state=num_variables()-1-j+countstate;
+                    countstate++;
+                    if(check_val-getsum>=0) check_val=check_val-getsum;
+                    getsum=getsum/10;
+                    if (upordown==1) {
+                      tmps=colvar_sigmas[state];
+                      tmp=curr_cv_values[state]; // store original current cv value
+                      unitary=curr_cv_values[state];
+                      unitary.set_ones();
+                      valk=std::floor(checkstate/getsumk);
+                      if(checkstate-getsumk>=0) checkstate=checkstate-getsumk;
+                      getsumk=getsumk/10;
+                      reflection_limit=reflection_l[state][valk];
+                      tmpd=reflection_limit-cvm::real(curr_cv_values[state]);
+                      tmpd=std::sqrt(tmpd*tmpd);
+                      if (tmpd<reflection_int*tmps && reflection_usel[state][valk] ) { // do mirror within selected range in case upordown=1
+                        curr_cv_values[state]=2.0*reflection_limit*unitary-tmp; // reflected cv value
+                      } else {
+                        hill_add=false;
+                      }
+                    }
+                 }
+                 if (hill_add) {
+                   switch (comm) {
+
+                   case single_replica:
+
+                     add_hill(hill(cvm::step_absolute(), hill_weight*h_scale, curr_cv_values, colvar_sigmas));
+
+                     break;
+
+                   case multiple_replicas:
+                     add_hill(hill(cvm::step_absolute(), hill_weight*h_scale, curr_cv_values, colvar_sigmas, replica_id));
+                     std::ostream &replica_hills_os =
+                       cvm::proxy->output_stream(replica_hills_file, "replica hills file");
+                     if (replica_hills_os) {
+                       write_hill(replica_hills_os, hills.back());
+                     } else {
+                       return cvm::error("Error: in metadynamics bias \""+this->name+"\""+
+                                         ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
+                                         " while writing hills for the other replicas.\n", COLVARS_FILE_ERROR);
+                     }
+                     break;
+                   }
+
+                   for (i = 0; i < num_variables(); i++) {
+                      curr_cv_values[i] = colvar_values[i]; // go back to previous values
+                   }
+                 } else {
+                   for (i = 0; i < num_variables(); i++) {
+                      curr_cv_values[i] = colvar_values[i]; // go back to previous values
+                   }
+                 }
+              }
+           }
+     }
+  }
+  return COLVARS_OK;
+}
 
 std::list<colvarbias_meta::hill>::const_iterator
 colvarbias_meta::delete_hill(hill_iter &h)
@@ -551,7 +1011,7 @@ int colvarbias_meta::update_bias()
     cvm::real hills_scale=1.0;
 
     if (ebmeta) {
-      hills_scale *= 1.0/target_dist->value(target_dist->get_colvars_index());
+      hills_scale *= 1.0/target_dist->value(target_dist->get_colvars_index_bound());
       if(cvm::step_absolute() <= ebmeta_equil_steps) {
         cvm::real const hills_lambda =
           (cvm::real(ebmeta_equil_steps - cvm::step_absolute())) /
@@ -571,29 +1031,50 @@ int colvarbias_meta::update_bias()
       hills_scale *= cvm::exp(-1.0*hills_energy_sum_here/(bias_temperature*proxy->boltzmann()));
     }
 
-    switch (comm) {
+    // Whether add a hill
+    bool add_a_hill=true;
 
-    case single_replica:
+    // Do not add hills beyond reflection borders
+    // as just reflected hills must be present
+    // beyond those boundaries
 
-      add_hill(hill(cvm::step_absolute(), hill_weight*hills_scale,
-                    colvar_values, colvar_sigmas));
+    // Check reflection borders: if beyond borders do not add hill
 
-      break;
+    add_a_hill=check_reflection_limits(add_a_hill);
 
-    case multiple_replicas:
-      add_hill(hill(cvm::step_absolute(), hill_weight*hills_scale,
-                    colvar_values, colvar_sigmas, replica_id));
-      std::ostream &replica_hills_os =
-        cvm::proxy->output_stream(replica_hills_file, "replica hills file");
-      if (replica_hills_os) {
-        write_hill(replica_hills_os, hills.back());
-      } else {
-        return cvm::error("Error: in metadynamics bias \""+this->name+"\""+
-                          ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
-                          " while writing hills for the other replicas.\n", COLVARS_FILE_ERROR);
+    if (add_a_hill) {
+
+      switch (comm) {
+
+      case single_replica:
+
+        add_hill(hill(cvm::step_absolute(), hill_weight*hills_scale,
+                      colvar_values, colvar_sigmas));
+
+        break;
+
+      case multiple_replicas:
+        add_hill(hill(cvm::step_absolute(), hill_weight*hills_scale,
+                      colvar_values, colvar_sigmas, replica_id));
+        std::ostream &replica_hills_os =
+          cvm::proxy->output_stream(replica_hills_file, "replica hills file");
+        if (replica_hills_os) {
+          write_hill(replica_hills_os, hills.back());
+        } else {
+          return cvm::error("Error: in metadynamics bias \""+this->name+"\""+
+                            ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
+                            " while writing hills for the other replicas.\n", COLVARS_FILE_ERROR);
+        }
+        break;
       }
-      break;
+
+      // add reflected hills if required
+      if (use_reflection) {
+        reflect_hill_multid(hills_scale);
+      }
+
     }
+
   }
 
   return COLVARS_OK;
@@ -604,7 +1085,9 @@ int colvarbias_meta::update_grid_data()
 {
   if ((cvm::step_absolute() % grids_freq) == 0) {
     // map the most recent gaussians to the grids
-    project_hills(new_hills_begin, hills.end(), hills_energy.get(), hills_energy_gradients.get());
+    project_hills(new_hills_begin, hills.end(), hills_energy.get(), hills_energy_gradients.get(),
+                  which_int_llimit_cv, which_int_ulimit_cv,
+                  reflection_llimit, reflection_ulimit);
     new_hills_begin = hills.end();
 
     // TODO: we may want to condense all into one replicas array,
@@ -614,7 +1097,9 @@ int colvarbias_meta::update_grid_data()
         replicas[ir]->project_hills(replicas[ir]->new_hills_begin,
                                     replicas[ir]->hills.end(),
                                     replicas[ir]->hills_energy.get(),
-                                    replicas[ir]->hills_energy_gradients.get());
+                                    replicas[ir]->hills_energy_gradients.get(),
+                                    which_int_llimit_cv, which_int_ulimit_cv,
+                                    reflection_llimit, reflection_ulimit);
         replicas[ir]->new_hills_begin = replicas[ir]->hills.end();
       }
     }
@@ -628,24 +1113,33 @@ int colvarbias_meta::calc_energy(std::vector<colvarvalue> const *values)
 {
   size_t ir = 0;
 
-  for (ir = 0; ir < replicas.size(); ir++) {
-    replicas[ir]->bias_energy = 0.0;
-  }
-
-  bool index_ok = false;
+  size_t i;
+  int icount;
   std::vector<int> curr_bin;
 
-  if (use_grids) {
+  curr_values = values ? *values : colvar_values;
 
-    curr_bin = values ?
-      hills_energy->get_colvars_index(*values) :
-      hills_energy->get_colvars_index();
-
-    index_ok = hills_energy->index_ok(curr_bin);
+  if (use_reflection) {
+    for (i = 0; i < num_variables(); i++) {
+       icount=which_int_llimit_cv[i];
+       if (icount>-1 && curr_values[i]<reflection_llimit[icount] ) {
+         curr_values[i]=reflection_llimit[icount];
+       }
+       icount=which_int_ulimit_cv[i];
+       if (icount>-1 && curr_values[i]>reflection_ulimit[icount] ) {
+         curr_values[i]=reflection_ulimit[icount];
+       }
+    }
+    curr_bin = hills_energy->get_colvars_index_bound(curr_values);
+  } else {
+    curr_bin = hills_energy->get_colvars_index(curr_values);
+  }
 
+  for (ir = 0; ir < replicas.size(); ir++) {
+    replicas[ir]->bias_energy = 0.0;
   }
 
-  if ( index_ok ) {
+  if (hills_energy->index_ok(curr_bin)) {
     // index is within the grid: get the energy from there
     for (ir = 0; ir < replicas.size(); ir++) {
 
@@ -660,11 +1154,13 @@ int colvarbias_meta::calc_energy(std::vector<colvarvalue> const *values)
     }
   } else {
     // off the grid: compute analytically only the hills at the grid's edges
-    for (ir = 0; ir < replicas.size(); ir++) {
-      calc_hills(replicas[ir]->hills_off_grid.begin(),
-                 replicas[ir]->hills_off_grid.end(),
-                 bias_energy,
-                 values);
+    if (!use_reflection) {
+      for (ir = 0; ir < replicas.size(); ir++) {
+        calc_hills(replicas[ir]->hills_off_grid.begin(),
+                   replicas[ir]->hills_off_grid.end(),
+                   bias_energy,
+                   &curr_values);
+      }
     }
   }
 
@@ -675,7 +1171,7 @@ int colvarbias_meta::calc_energy(std::vector<colvarvalue> const *values)
     calc_hills(replicas[ir]->new_hills_begin,
                replicas[ir]->hills.end(),
                bias_energy,
-               values);
+               &curr_values);
     if (cvm::debug()) {
       cvm::log("Hills energy = "+cvm::to_str(bias_energy)+".\n");
     }
@@ -688,42 +1184,61 @@ int colvarbias_meta::calc_energy(std::vector<colvarvalue> const *values)
 int colvarbias_meta::calc_forces(std::vector<colvarvalue> const *values)
 {
   size_t ir = 0, ic = 0;
+  int icount;
+  std::vector<int> curr_bin;
+  curr_values = values ? *values : colvar_values;
+  std::vector<bool> add_force(num_variables());
   for (ir = 0; ir < replicas.size(); ir++) {
     for (ic = 0; ic < num_variables(); ic++) {
       replicas[ir]->colvar_forces[ic].reset();
     }
   }
 
-  bool index_ok = false;
-  std::vector<int> curr_bin;
-
-  if (use_grids) {
-
-    curr_bin = values ?
-      hills_energy->get_colvars_index(*values) :
-      hills_energy->get_colvars_index();
-
-    index_ok = hills_energy->index_ok(curr_bin);
+  for (ic = 0; ic < num_variables(); ic++) {
+    add_force[ic]=true;
+    icount=which_int_llimit_cv[ic];
+    if (icount>-1) {
+      if ( curr_values[ic]<reflection_llimit[icount] ) {
+        curr_values[ic]=reflection_llimit[icount];
+        add_force[ic]=false;
+      }
+    }
+    icount=which_int_ulimit_cv[ic];
+    if (icount>-1) {
+      if ( curr_values[ic]>reflection_ulimit[icount] ) {
+        add_force[ic]=false;
+        curr_values[ic]=reflection_ulimit[icount];
+      }
+    }
+  }
 
+  if (use_reflection) {
+    curr_bin = hills_energy->get_colvars_index_bound(curr_values);
+  } else {
+    curr_bin = hills_energy->get_colvars_index(curr_values);
   }
 
-  if ( index_ok ) {
+  if (hills_energy->index_ok(curr_bin)) {
     for (ir = 0; ir < replicas.size(); ir++) {
       cvm::real const *f = &(replicas[ir]->hills_energy_gradients->value(curr_bin));
       for (ic = 0; ic < num_variables(); ic++) {
         // the gradients are stored, not the forces
-        colvar_forces[ic].real_value += -1.0 * f[ic];
+        if (add_force[ic]) {
+          colvar_forces[ic].real_value += -1.0 * f[ic];
+        }
       }
     }
   } else {
     // off the grid: compute analytically only the hills at the grid's edges
-    for (ir = 0; ir < replicas.size(); ir++) {
-      for (ic = 0; ic < num_variables(); ic++) {
-        calc_hills_force(ic,
-                         replicas[ir]->hills_off_grid.begin(),
-                         replicas[ir]->hills_off_grid.end(),
-                         colvar_forces,
-                         values);
+    if (!use_reflection) {
+      for (ir = 0; ir < replicas.size(); ir++) {
+        for (ic = 0; ic < num_variables(); ic++) {
+           calc_hills_force(ic,
+                            replicas[ir]->hills_off_grid.begin(),
+                            replicas[ir]->hills_off_grid.end(),
+                            colvar_forces,
+                            &curr_values);
+        }
       }
     }
   }
@@ -739,11 +1254,13 @@ int colvarbias_meta::calc_forces(std::vector<colvarvalue> const *values)
 
   for (ir = 0; ir < replicas.size(); ir++) {
     for (ic = 0; ic < num_variables(); ic++) {
-      calc_hills_force(ic,
-                       replicas[ir]->new_hills_begin,
-                       replicas[ir]->hills.end(),
-                       colvar_forces,
-                       values);
+      if (add_force[ic]) {
+        calc_hills_force(ic,
+                         replicas[ir]->new_hills_begin,
+                         replicas[ir]->hills.end(),
+                         colvar_forces,
+                         &curr_values);
+      }
       if (cvm::debug()) {
         cvm::log("Hills forces = "+cvm::to_str(colvar_forces)+".\n");
       }
@@ -864,6 +1381,10 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter  h_first,
                                     colvarbias_meta::hill_iter  h_last,
                                     colvar_grid_scalar         *he,
                                     colvar_grid_gradient       *hg,
+                                    std::vector<int> const &w_int_llimit_cv,
+                                    std::vector<int> const &w_int_ulimit_cv,
+                                    std::vector<cvm::real> const &ref_llimit,
+                                    std::vector<cvm::real> const &ref_ulimit,
                                     bool print_progress)
 {
   if (cvm::debug())
@@ -874,6 +1395,7 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter  h_first,
   // TODO: improve it by looping over a small subgrid instead of the whole grid
 
   std::vector<colvarvalue> new_colvar_values(num_variables());
+
   std::vector<cvm::real> colvar_forces_scalar(num_variables());
 
   std::vector<int> he_ix = he->new_index();
@@ -887,12 +1409,29 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter  h_first,
   if (hg != NULL) {
 
     // loop over the points of the grid
+    size_t i;
+    int icount;
+    std::vector<bool> add_force(num_variables());
     for ( ;
           (he->index_ok(he_ix)) && (hg->index_ok(hg_ix));
           count++) {
-      size_t i;
       for (i = 0; i < num_variables(); i++) {
+        add_force[i]=true;
         new_colvar_values[i] = he->bin_to_value_scalar(he_ix[i], i);
+        icount=w_int_llimit_cv[i];
+        if (icount>-1 ){
+          if ( new_colvar_values[i]<ref_llimit[icount] ) {
+            new_colvar_values[i]=ref_llimit[icount];
+            add_force[i]=false;
+          }
+        }
+        icount=w_int_ulimit_cv[i];
+        if (icount>-1){
+          if( new_colvar_values[i]>ref_ulimit[icount] ) {
+            new_colvar_values[i]=ref_ulimit[icount];
+            add_force[i]=false;
+          }
+        }
       }
 
       // loop over the hills and increment the energy grid locally
@@ -902,7 +1441,9 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter  h_first,
 
       for (i = 0; i < num_variables(); i++) {
         hills_forces_here[i].reset();
-        calc_hills_force(i, h_first, h_last, hills_forces_here, &new_colvar_values);
+        if (add_force[i]){
+          calc_hills_force(i, h_first, h_last, hills_forces_here, &new_colvar_values);
+        }
         colvar_forces_scalar[i] = hills_forces_here[i].real_value;
       }
       hg->acc_force(hg_ix, &(colvar_forces_scalar.front()));
@@ -1460,8 +2001,10 @@ void colvarbias_meta::rebin_grids_after_restart()
       // if there are hills, recompute the new grids from them
       cvm::log("Rebinning the energy and forces grids from "+
                cvm::to_str(hills.size())+" hills (this may take a bit)...\n");
-      project_hills(hills.begin(), hills.end(), new_hills_energy.get(),
-                    new_hills_energy_gradients.get(), true);
+      project_hills(hills.begin(), hills.end(),
+                    new_hills_energy.get(), new_hills_energy_gradients.get(),
+                    which_int_llimit_cv, which_int_ulimit_cv,
+                    reflection_llimit, reflection_ulimit, true);
       cvm::log("rebinning done.\n");
 
     } else {
@@ -1819,7 +2362,10 @@ template <typename OST> OST &colvarbias_meta::write_state_data_template_(OST &os
 
     // this is a very good time to project hills, if you haven't done
     // it already!
-    project_hills(new_hills_begin, hills.end(), hills_energy.get(), hills_energy_gradients.get());
+    project_hills(new_hills_begin, hills.end(),
+                  hills_energy.get(), hills_energy_gradients.get(),
+                  which_int_llimit_cv, which_int_ulimit_cv,
+                  reflection_llimit, reflection_ulimit);
     new_hills_begin = hills.end();
 
     // write down the grids to the restart file
diff --git a/src/colvarbias_meta.h b/src/colvarbias_meta.h
index 57aa21ed6..8e745c187 100644
--- a/src/colvarbias_meta.h
+++ b/src/colvarbias_meta.h
@@ -48,6 +48,7 @@ class colvarbias_meta
   virtual int init_replicas_params(std::string const &conf);
   virtual int init_well_tempered_params(std::string const &conf);
   virtual int init_ebmeta_params(std::string const &conf);
+  virtual int init_reflection_params(std::string const &conf);
 
   virtual int clear_state_data();
 
@@ -154,6 +155,14 @@ class colvarbias_meta
   /// the next hill in the list)
   std::list<hill>::const_iterator delete_hill(hill_iter &h);
 
+
+  /// \brief Check is current colvar value is within inversion or
+  /// reflection limits to assess whether to add a hill
+  bool check_reflection_limits(bool &ah);
+
+  /// \brief Multidimensional routine to reflect hills
+  int reflect_hill_multid(cvm::real const &h_scale);
+
   /// \brief Calculate the values of the hills, incrementing
   /// bias_energy
   virtual void calc_hills(hill_iter  h_first,
@@ -210,6 +219,38 @@ class colvarbias_meta
   /// \brief Biasing temperature in well-tempered metadynamics
   cvm::real  bias_temperature;
 
+  // hills reflection
+
+  /// \brief Whether using hills reflection
+  bool use_reflection;
+
+  /// \brief For which variables reflection limits are on
+
+  std::vector<int> reflection_llimit_cv;
+  std::vector<int> reflection_ulimit_cv;
+  cvm::real reflection_int;
+  size_t nrefvarsl;
+  size_t nrefvarsu;
+
+  /// \brief Limits for reflection
+  std::vector<cvm::real> reflection_llimit;
+  std::vector<cvm::real> reflection_ulimit;
+
+  /// \brief Multidimensional reflection : store cvs to use and pointers to the limits
+  std::vector<std::vector<bool> > reflection_usel;
+  std::vector<std::vector<cvm::real> > reflection_l;
+
+  /// \brief Multidimensional reflection states
+  std::vector<std::vector<size_t> > ref_state;
+
+  /// \brief For which variables hills forces beyond the boundaries(interval) must be removed
+
+  std::vector<int> which_int_llimit_cv;
+  std::vector<int> which_int_ulimit_cv;
+
+  /// \brief Current value of colvars to be modifed for calculation of energy and forces with interval
+  std::vector<colvarvalue> curr_values;
+
   /// Ensemble-biased metadynamics (EBmeta) flag
   bool       ebmeta;
 
@@ -232,8 +273,11 @@ class colvarbias_meta
   std::shared_ptr<colvar_grid_gradient> hills_energy_gradients;
 
   /// Project the selected hills onto grids
-  void project_hills(hill_iter h_first, hill_iter h_last, colvar_grid_scalar *ge,
-                     colvar_grid_gradient *gf, bool print_progress = false);
+  void project_hills(hill_iter h_first, hill_iter h_last,
+                     colvar_grid_scalar *ge, colvar_grid_gradient *gf,
+                     std::vector<int> const &w_int_llimit_cv, std::vector<int> const &w_int_ulimit_cv,
+                     std::vector<cvm::real> const &ref_llimit, std::vector<cvm::real> const &ref_ulimit,              
+                     bool print_progress = false);
 
 
   // Multiple Replicas variables and functions
diff --git a/src/colvargrid.h b/src/colvargrid.h
index e31c424a5..a0ac538b5 100644
--- a/src/colvargrid.h
+++ b/src/colvargrid.h
@@ -714,6 +714,17 @@ template <class T> class colvar_grid : public colvar_grid_params, public colvarp
 
   /// \brief Get the bin indices corresponding to the provided values of
   /// the colvars and assign first or last bin if out of boundaries
+  inline std::vector<int> const get_colvars_index_bound(std::vector<colvarvalue> const &values) const
+  {
+    std::vector<int> index = new_index();
+    for (size_t i = 0; i < nd; i++) {
+      index[i] = value_to_bin_scalar_bound(values[i], i);
+    }
+    return index;
+  }
+
+  /// \brief Get the bin indices corresponding to the current values of
+  /// the colvars and assign first or last bin if out of boundaries
   inline std::vector<int> const get_colvars_index_bound() const
   {
     std::vector<int> index = new_index();