# Consensus-based construction of high-dimensional free energy surface

Liyao Lyu
Department of Computational Mathematics, Science & Engineering, Michigan State
University, MI 48824, USA &Huan Lei
Department of Computational Mathematics, Science & Engineering, Michigan State
University, MI 48824, USA
Department of Statistics & Probability, Michigan State University, MI 48824,
USA

One essential problem in quantifying the collective behaviors of molecular
systems lies in the accurate construction of free energy surfaces (FESs). The
main challenges arise from the prevalence of energy barriers and the high
dimensionality. Existing approaches are often based on sophisticated enhanced
sampling methods to establish efficient exploration of the full phase space.
On the other hand, the collection of optimal sample points for the numerical
approximation of FESs remains largely under-explored, where the discretization
error could become dominant for systems with a large number of collective
variables (CVs). We propose a consensus sampling based approach by
reformulating the construction as a minimax problem which simultaneously
optimizes the function representation and the training set. In particular, the
maximization step establishes a stochastic interacting particle system to
achieve the adaptive sampling of the max-residue regime by modulating the
exploitation of the Laplace approximation of the current loss function and the
exploration of the uncharted phase space; the minimization step updates the
FES approximation with the new training set. By iteratively solving the
minimax problem, the present method essentially achieves an adversarial
learning of the FESs with unified tasks for both phase space exploration and
posterior error enhanced sampling. We demonstrate the method by constructing
the FESs of molecular systems with a number of CVs up to $30$.

_Keywords_ Free energy $|$ Coarse-graining $|$ Machine learning $|$ Molecular
simulation $|$ Adversarial learning

Molecular dynamics (MD) provides an essential tool [1] to access the micro-
scale insights of complex processes in many scientific applications, including
chemical engineering, material synthesis, and drug design. To further predict
the collective behaviors, it is often desirable to construct the free energy
surfaces (FESs) with respect to a set of functional-relevant collective
variables (CVs). However, the accurate construction of FES remains a practical
challenge due to the high dimensionality and the prevalence of energy
barriers. Several ingenious methods based on importance sampling have been
developed, such as umbrella sampling [2], histogram reweighting [3],
metadynamics [4, 5], variational enhanced sampling [6, 7, 8], weighted
ensemble and its variants [9, 10, 11], hyperdynamics [12], orthogonal space
random walk [13] and adaptive biasing force [14, 15, 16]. Despite the efficacy
in enhancing the phase space exploration, these methods are based on iterative
estimation of the biased probability density functions (PDFs) (or biased
force), where the efficiency generally diminishes as the number of CVs
increases. Strategies such as bias-exchange [17] and parallel bias
metadynamics [18, 19] have been developed to alleviate this issue. Both
methods factorize the biasing potential as a product of several low-
dimensional functions and can work effectively when the underlying FESs admit
such structures which, however, may not always be the case for high-
dimensional problems. Several other methods [20, 21, 22, 23, 24, 25, 26, 27,
28] achieve enhanced sampling irrespective of pre-defined CVs (see review
[29]). However, the construction of explicit high-dimensional FESs remains a
challenging problem.

Alternative to the histogram methods, the temperature-accelerated molecular
dynamics [30] and the adiabatic free energy dynamics [31, 32] are developed by
introducing an extended dynamics of the CVs with an artificially high
temperature to overcome the energy barriers. The methods avoid dealing with
the biased potentials and have shown the promise of efficient exploration of
the high-dimensional CV space [33]. On the other hand, to further construct
the FES, the sampling points are deposited in a greedy way [34]; the
effectiveness for high-dimensional problems remains under-explored.

With the recent advancements in machine learning, kernel methods and deep
neural networks (DNNs) have demonstrated promising results in representing
multi-dimensional FES [35, 36, 37, 8, 38]. However, as also noted in Ref.
[39], one essential problem for constructing high-dimensional FES lies in how
to optimize the training set. Ideally, to efficiently approximate a high-
dimensional function, the sampling should enable certain adaptivity based on
the _posterior_ error so that the training set can be updated according to the
particular underlying FES. However, such _posterior_ error-based sampling
remains largely open for FES construction since the approximation residual can
neither be easily queried nor sampled within the full phase space. In
practice, the training samples are often collected either as pre-defined grid
points or in a heuristic manner. As the number of CVs increases, such choices
may induce pronounced approximation error and suffer from performance
degradation, since most phase space regimes essentially become
thermodynamically inaccessible. Remarkably, the DeepVes [8] enables the
efficient sampling of the CV space by jointly constructing the bias potential
and the target distribution. On the other hand, the accurate reconstruction of
explicit FES further requires the estimation of high-dimensional PDFs from the
samples. Recent efforts [40, 41] seek the direct approximation of the
transition operator instead of estimating PDFs using generative models, when
representative configurations are known _a priori_. In addition, the
reinforced dynamics [42, 43] proposes an efficient approach to impose
adaptivity by using an uncertainty indicator as an indirect measure of the
construction error to bias MD simulations, which relies on calculating the
standard deviation of the predictions from multiple DNNs trained on the same
dataset.

In this work, we present a consensus-based enhanced sampling (CES) method to
efficiently construct high-dimensional FES for complex MD systems. A unique
feature is that the method enables adaptive posterior error-based sampling
such that the FES approximation and the sample points can be simultaneously
optimized. Also, unlike the reinforced dynamics, the method does not rely on
the training of multiple DNN representations. The main idea is to formulate
the construction as a min-max problem, where the max-problem seeks a residue-
based distribution to establish adaptive sampling in the vicinity of the
explored phase space regime, and the min-problem optimizes the DNN parameters
for the FES representation. Iteratively optimizing both the training set and
the DNN representation achieves an adversarial construction of the FES
pertaining to the kinetic transition regime. For the maximization step, it is
worth mentioning that the query of the thermodynamically accessible phase
space is non-trivial. In particular, existing approaches [44, 45] in the
spirit of adversarial generative models can not be directly applied, since the
global residual error is unknown _a prior_. Instead, we establish a consensus-
based sampling [46] (see also Ref. [47, 48, 49, 50] ) in the form of a
stochastic interacting particle system governed by a McKean stochastic
differential equation. A quadratic potential is adaptively constructed to
probe the local maximum error regime by exploiting the Laplace approximation
under a low-temperature limit. Meanwhile, a coherent noise term is introduced
to efficiently explore the full CV space under a high-temperature limit and
yield the target sampling points for the maximization problem, by which the
minimization of the DNN residue can be solved jointly. We demonstrate the
effectiveness of the proposed method by constructing the FES of several
biomolecule systems with the number of CVs up to 30. Figure 1 sketches the
workflow of the proposed method.

Figure 1: The workflow of the present CES-based method. In the minimization
step, given a collection of sampling points, the reference force (i.e., the
gradient of the underlying FES) can be calculated using the restrained
dynamics; a comparison with the force inferred from the DNN approximation
yields a posterior error. In the maximization step, the posterior error
determines a residue-based distribution with entropy regularization. An
interacting particle system following McKean stochastic dynamics is used to
achieve an adaptive sampling of the max-residue regime by modulating the
exploitation of the Laplace approximation of the current residue-based
distribution and the exploration of the uncharted phase space. The FES can be
accurately reconstructed after several iterations of the minimization and
maximization step.

## Consensus-based construction method

### Free energy and mean forces

We consider a full model with micro-scale coordinates
$\mathbf{r}\in\mathbb{R}^{N}$ whose dynamics is governed by potential
$U(\mathbf{r}):\mathbb{R}^{N}\rightarrow\mathbb{R}$ under temperature $T$.
Suppose we are interested in CVs
$\mathbf{s}(\mathbf{r}):\mathbb{R}^{N}\to\Gamma$ with
$\Gamma\subset\mathbb{R}^{M}$, the FES $A(\mathbf{z})$ of the CVs is defined
by

$\displaystyle A(\mathbf{z})=-\frac{1}{\beta}\ln\rho(\mathbf{z}),$ (1)

where $\beta=1/k_{B}T$ is the inverse of the thermal temperature,
$\rho(\mathbf{z})=\frac{1}{Z}\int\exp{(-\beta
U(\mathbf{r}))}\delta(\mathbf{s}(\mathbf{r})-\mathbf{z})\mathrm{d}\mathbf{r}$
is the marginal density for $\mathbf{s}(\mathbf{r})=\mathbf{z}$,
$Z=\int\exp{(-\beta U(\mathbf{r}))}\mathrm{d}\mathbf{r}$ is the partition
function, and $\delta(\cdot)$ represents the Dirac delta function. For high-
dimensional CVs, the direct estimation of $\rho(\mathbf{z})$ often becomes
numerically challenging. An alternative approach is to fit the mean force
$\mathbf{F}(\mathbf{z}):=-\nabla A(\mathbf{z})$ at various sample points,
which can be estimated via the restrained dynamics [51] by introducing a
harmonic term into the full potential, i.e.,

$U_{k}(\mathbf{r},\mathbf{z})=U(\mathbf{r})+\frac{k}{2}(\mathbf{s}(\mathbf{r})-\mathbf{z})^{T}(\mathbf{s}(\mathbf{r})-\mathbf{z}),$
(2)

where $k$ represents the magnitude of the restrained potential. As shown in
Ref [30], the mean force can be computed by
$\mathbf{F}(\mathbf{z})=\lim_{k\to\infty}\mathbf{F}^{k}(\mathbf{z})$, where
$\mathbf{F}^{k}(\mathbf{z})$ is defined by

$\mathbf{F}^{k}(\mathbf{z})=\frac{1}{Z_{k}(\mathbf{z})}\int
k(\mathbf{s}(\mathbf{r})-\mathbf{z})\exp{(-\beta
U_{k}(\mathbf{r},\mathbf{z}))}\mathrm{d}\mathbf{r},$ (3)

and can be sampled as the local first-moment estimation.

In principle, given a collection of sample points $\mathbf{z}$,
$A(\mathbf{z})$ can be re-constructed (up to a constant) by matching the mean
force $-\nabla A(\mathbf{z})$ at the individual points. However, for complex
MD systems, the sampling over the phase space could be highly non-trivial due
to the prevalence of local energy minima; the training set is often determined
_a prior_ as pre-selected points or in a greedy manner. As the number of CVs
increases, the empirical random samples may introduce pronounced
discretization error of the loss function over the full phase space. To
efficiently construct $A(\mathbf{z})$ in the thermodynamically accessible
regime, it is desirable to simultaneously optimize the training set and the
FES approximation through certain adaptive sampling based on the _posterior_
residual error. This motivates the present method illustrated below.

### Min-Max formulation

Let $A_{\mathcal{N}}(\mathbf{z})$ denote the DNN representation of the FES
$A(\mathbf{z})$, which is parameterized by minimizing the loss function

$\mathcal{L}_{\mathcal{N}}(\mathbf{z})=\|\nabla_{\mathbf{z}}A_{\mathcal{N}}(\mathbf{z})+\mathbf{F}(\mathbf{z}))\|^{2}$
(4)

for $\mathbf{z}\in\Gamma$. To solve the problem, we introduce a sampling
distribution $q(\mathbf{z})\in W$, where $W$ is a function space that defines
a proper constraint of the PDF $q(\mathbf{z})$ on $\Gamma$, and define the
weighted loss

$(\mathcal{L}_{\mathcal{N}},q)=\int_{\Gamma}\mathcal{L}_{\mathcal{N}}(\mathbf{z})q(\mathbf{z})\mathrm{d}\mathbf{z}.$
(5)

A desired distribution intends to maximize the discrepancy in the dataset for
a given network $A_{\mathcal{N}}(\mathbf{z})$. Accordingly, we define the
maximum problem as

$G[\mathcal{L}_{\mathcal{N}}]=\max_{q}(\mathcal{L}_{\mathcal{N}},q).$ (6)

Since $(\mathcal{L}_{\mathcal{N}},q)$ is always non-negative, a good
approximation of the original free energy surface
$A_{\mathcal{N}}(\mathbf{z})$ (up to a constant) can be obtained by solving
the following problems

$\min_{A_{\mathcal{N}}}\max_{q}(\mathcal{L}_{\mathcal{N}},q).$ (7)

######  Proposition 1.

Assuming that there exists a solution of
$\mathcal{L}_{\mathcal{N}}(\mathbf{z})=0$ for $\mathbf{z}\in\Gamma$, then
$A^{\ast}$ is a solution if and only if it solves (7).

###### Proof.

Suppose $A^{\ast}$ is the solution for
$\mathcal{L}_{\mathcal{N}}(\mathbf{z})=0$, it satisfies
$(\mathcal{L}_{\mathcal{N}^{\ast}},q)=0$ for any $q\in W$, i.e.,
$G[\mathcal{L}_{\mathcal{N}^{\ast}}]=0$. Therefore, $A^{\ast}$ is a solution
for the minimax problem (7). On the other hand, if $\hat{A}$ is the minimizer
for problem (7) but not the solution for
$\mathcal{L}_{\mathcal{N}}(\mathbf{z})=0$, then there exists $\hat{q}\in V$
such that $(\mathcal{L}_{\hat{\mathcal{N}}},\hat{q})>0$. However,
$(\mathcal{L}_{\mathcal{N}^{\ast}},q)=0$ for all $q$, which contradicts the
assumption that $\hat{A}$ is the minimizer. ∎

Proposition 1 shows that the direct construction of FES $A(\mathbf{z})$ can be
reformulated as an adversarial learning of the optimal solution
$A_{\mathcal{N}}(\mathbf{z})$ for the min-max problem (7). Accordingly, the
training consists of two components: the minimization part optimizes the DNN
representation with the current training set; the maximization part explores
the regime of the largest residue for the current DNN representation and
essentially establishes an adaptive sampling of the training set based on the
_posterior_ error.

To numerically solve the max-problem, certain regularization needs to be
introduced. Otherwise, the max-problem will simply yield a delta measure,
i.e., $\delta(\mathbf{z}-\mathbf{z}^{*})$, where
$\mathbf{z}^{\ast}=\arg\max\mathcal{L}_{\mathcal{N}}(\mathbf{z})$. Since the
sampling needs to simultaneously identify the max-residual regime and explore
the uncharted phase space, we introduce the entropy-based regularization [52,
53] (see also Refs. [54, 55, 44] for gradient-based regularization), and the
max-problem is reformulated by

$\min_{q}\int(-\mathcal{L}_{\mathcal{N}}(\mathbf{z})+\kappa_{h}^{-1}\ln
q(\mathbf{z}))q(\mathbf{z})\mathrm{d}\mathbf{z}.$ (8)

The problem is convex for a PDF $q$ with a unique global minimum at
$q^{*}(\mathbf{z})=\exp(-\kappa_{h}\mathcal{L}_{\mathcal{N}}^{-}(\mathbf{z}))/Z^{*}$,
where
$Z^{*}=\int\exp(-\kappa_{h}\mathcal{L}_{\mathcal{N}}^{-}(\mathbf{z}))\mathrm{d}\mathbf{z}$
and
$\mathcal{L}_{\mathcal{N}}^{-}(\mathbf{z})=-\mathcal{L}_{\mathcal{N}}(\mathbf{z})$.
The parameter $\kappa_{h}$ is a Lagrangian multiplier, balancing the focus
between the peak concentration and the scope of exploration, and is somewhat
similar to the inverse temperature in statistical mechanics. An elevated
$\kappa_{h}^{-1}$ induces a distribution closer to a uniform distribution.
Conversely, a diminished $\kappa_{h}^{-1}$ induces a concentrated distribution
near the max-residue point.

Given $q^{\ast}(\mathbf{z})$, the direct sampling remains challenging for
high-dimensional CVs. Inspired by the consensus-based sampling method [46], we
establish a stochastic particle system. By properly constructing the
conservative potential and noise terms, the particle distribution provides a
good approximation of $q^{\ast}$ as illustrated below.

### Exploitation and exploration in the max-problem

To approximate the target distribution $q^{\ast}(\mathbf{z})$, particularly in
the vicinity of the max-residual point $\mathbf{z}^{*}$, we exploit Laplace’s
principle in the large deviations theory, i.e.,

$\lim_{\kappa\to\infty}\left(-\frac{1}{\kappa}\log\left(\int\exp{(-\kappa
f(\mathbf{z}))}\mathrm{d}\rho^{*}(\mathbf{z})\right)\right)=f(\mathbf{z}^{*})$
(9)

holds true for any compactly supported probability measure $\rho^{\ast}$,
where $\mathbf{z}^{\ast}\in\operatorname{supp}(\rho^{*})$ uniquely minimizes
the function $f$. This enables us to identify the max-residual point from a
collection of samples $\left\\{\mathbf{z}^{i}\right\\}_{i=1}^{N_{w}}$ by the
first-order momentum $\mathbf{m}$ under the weighted density function
$p(\mathbf{z})$, i.e.,

$\mathbf{m}=\int\mathbf{z}p(\mathbf{z})\mathrm{d}\mathbf{z}\approx\sum_{i=1}^{N_{w}}\mathbf{z}^{i}\hat{p}(\mathbf{z}^{i})\quad\hat{p}(\mathbf{z})=\frac{\exp{(-\kappa_{l}\mathcal{L}_{\mathcal{N}}^{-}(\mathbf{z}))}}{\sum_{i=1}^{N_{w}}\exp{(-\kappa_{l}\mathcal{L}_{\mathcal{N}}^{-}(\mathbf{z}^{i})})},$
(10)

where $\kappa_{l}^{-1}$ represents a low temperature limit. However, the
integration is subject to the so-called curse of dimensionality as the number
of CVs increases. Instead, we treat sampler $\mathbf{z}^{i}$ as a random
walker $\mathbf{z}^{i}_{t}$ governed by the following McKean stochastic
differential equation

$\dot{\mathbf{z}}^{i}_{t}=-\frac{1}{\gamma}\nabla_{\mathbf{z}}G(\mathbf{z}^{i}_{t})+\sqrt{\frac{2}{\kappa_{h}\gamma}}\mathbf{\xi}_{i}(t),$
(11)

where
$G(\mathbf{z}_{t})=\frac{1}{2}(\mathbf{z}_{t}-\mathbf{m}_{t})^{T}V_{t}^{-1}(\mathbf{z}_{t}-\mathbf{m}_{t})$
denotes an adaptively constructed conservative potential function and $\gamma$
is the friction coefficient. The formulations of $\mathbf{m}_{t},V_{t}$ are
specified in (12) with the rationale discussed in the next section.
Consequently, $G(\mathbf{z}_{t})$ navigates the random walkers (i.e.,
individual particles) towards $\mathbf{m}_{t}$, which represents the region of
large residual error. The second term in Eq. (11) represents a stochastic term
where $\gamma$ represents the friction coefficient and $\mathbf{\xi}(t)$
represents the standard Gaussian white noise characterized by zero mean and
covariance
$\mathbb{E}[\mathbf{\xi}_{i}(t)\mathbf{\xi}_{j}(t^{\prime})]=\delta_{ij}\delta(t-t^{\prime})$.

The coupling of the conservative and stochastic terms maintains a relatively
high temperature $\kappa_{h}^{-1}$, and a large friction coefficient $\gamma$
is applied such that the distribution is almost always Gaussian during the
evolution. As shown in the following section, the distribution of walkers
$q_{t}(\mathbf{z})$ converges to $\propto\exp(-\kappa_{h}G(\mathbf{z}))$
characterized by $\mathbf{m}_{\infty}$ and $V_{\infty}$. Accordingly, the
balance between exploitation and exploration is controlled using two
temperatures $\kappa_{l}^{-1}$ and $\kappa_{h}^{-1}$. As $\kappa_{l}^{-1}$
decreases, the random walker distribution concentrates near the max-residual
points of the current model, reflecting the role of exploitation. Conversely,
as $\kappa_{h}^{-1}$ increases, the distribution smoothens progressively,
enhancing the exploration of the uncharted regions.

### Adaptive parameter estimation and FES construction

In this section, we show that the sampling distribution governed by Eq. (11)
converges to a steady state as an approximation of $q^{\ast}$ with the proper
choices of $\mathbf{m}_{t}$ and $V_{t}$ given by

$\displaystyle\mathbf{m}_{t}$
$\displaystyle=\sum_{i=1}^{N_{w}}\mathbf{z}^{i}_{t}\hat{p}(\mathbf{z}^{i}_{t}),$
(12) $\displaystyle V_{t}$
$\displaystyle=\kappa_{t}\sum_{i=1}^{N_{w}}(\mathbf{z}^{i}_{t}-\mathbf{m}_{t})(\mathbf{z}^{i}_{t}-\mathbf{m}_{t})^{T}\hat{p}(\mathbf{z}^{i}_{t}),$

where $\hat{p}(\mathbf{z})$ is defined by (10). In particular, with
$\kappa_{t}=\kappa_{l}+\kappa_{h}$, (11) converges the target residue-based
distribution.

######  Proposition 2.

Suppose $\mathcal{L}_{\mathcal{N}}^{-}(\mathbf{z})$ takes a local quadratic
approximation in form of
$\frac{1}{2}(\mathbf{z}-\boldsymbol{\mu})^{T}\Sigma^{-1}(\mathbf{z}-\boldsymbol{\mu})$,
$q_{t}\to\frac{\exp{(-\kappa_{h}\mathcal{L}_{\mathcal{N}}^{-}(\mathbf{z}))}}{\int\exp{(-\kappa_{h}\mathcal{L}_{\mathcal{N}}^{-}(\mathbf{z}))}\mathrm{d}\mathbf{z}}$
as $t\to\infty$, given $\kappa_{t}=\kappa_{l}+\kappa_{h}$.

We refer to Appendix for the proof. The stochastic dynamics (11) and (12)
differs from the regular Langevin dynamics; it is nonlocal and similar to the
one in Ref. [46] except that a parameter $\kappa_{h}$ appears in the target
distribution that modulates the exploration of the sampling dynamics. The
quadratic assumption of the loss function
$\mathcal{L}_{\mathcal{N}}^{-}(\mathbf{z})$ is due to the fact that we are
mainly interested in the regime near the max-residue point. Under the low-
temperature limit, the local regime can be well characterized by the first and
second moments following Laplace’s principle. While $\mathbf{m}$ identifies
the extremal point, $V$ recognizes the anisotropic nature among the different
CVs. In this study, for the sake of computational efficiency, we further
simplify $V$ by only considering the diagonal entries denoted as $\mathbf{v}$.
Moreover, we utilize a moving average to ensure stable estimation

$\displaystyle\mathbf{m}_{t+1}$
$\displaystyle=\beta_{1}\mathbf{m}_{t}+(1-\beta_{1})\sum_{i=1}^{N_{w}}\mathbf{z}^{i}_{t}\hat{p}(\mathbf{z}^{i}_{t}),$
(13) $\displaystyle\mathbf{v}_{t+1}$
$\displaystyle=\beta_{2}\mathbf{v}_{t}+(1-\beta_{2})\kappa_{t}\sum_{i=1}^{N_{w}}(\mathbf{z}^{i}_{t}-\mathbf{m})\odot(\mathbf{z}^{i}_{t}-\mathbf{m})\hat{p}(\mathbf{z}^{i}_{t}),$

where $\beta_{1}$ and $\beta_{2}$ are hyper-parameters. To maintain unbiased
estimation, normalization is implemented as follows:
$\mathbf{m}=\frac{\mathbf{m}_{t}}{1-\beta_{1}^{t}}$ and
$\mathbf{v}=\frac{\mathbf{v}_{t}}{1-\beta_{2}^{t}}$.

With the training samples obtained from the aforementioned maximization step,
the DNN representation of the FES is optimized using the Adam stochastic
gradient descent method [56] for the minimization step. The loss function of
the updated DNN representation $A_{\mathcal{N}}(\mathbf{z})$, in turn,
navigates the consensus-based adaptive sampling for the updated maximization
step. The min-max problem is solved iteratively to achieve comprehensive
sampling of the full phase space.

So far, the construction process is based on the function approximation of
$A(\mathbf{z})$ over the full regime. However, we note that the kinetic
processes of a molecular system are generally characterized by the local
minima and saddles points on the FES. On the other hand, the regimes of high
free energy are less relevant. To accurately construct these thermodynamically
accessible regimes, we modify the loss function $\mathcal{L}$ as

$\mathcal{L}_{\mathcal{N}}(\mathbf{z})=\frac{\left\|\nabla_{\mathbf{z}}A_{\mathcal{N}}(\mathbf{z})+\mathbf{F}(\mathbf{z})\right\|^{2}}{\left\|\mathbf{F}(\mathbf{z})\right\|^{2}+e},$
(14)

for all the biomolecule systems except the toy example of the 1D Rastrigin
function. $e$ is a small value to regularize the denominator. We note that a
similar formulation is used in Ref. [34] to quantify the model accuracy after
the FES is constructed. Alternatively, in this study,
$\mathcal{L}_{\mathcal{N}}(\mathbf{z})$ directly involves the construction
process for the adaptive sampling of the training set. A detailed algorithm is
presented in Algorithm 1.

Algorithm 1 Consensus-based enhance sampling.

Initial sampling point $\mathbf{z}^{i}_{0}$, for $i=1,\ldots,N_{w}$

Initial DNN parameter $\theta_{0}$

The number of training iterations $N_{train}$

The number of data collected $N_{data}$ in each training iteration

$j\leftarrow 0,t\leftarrow 0$

$T\leftarrow\lceil\frac{N_{data}}{N_{w}}\rceil$

while $j<N_{train}$ do

while $t\leq T$ do

Calculate the mean force $\mathbf{F}^{i}_{t}$ at $\mathbf{z}^{i}_{t}$

Calculate the predicted force
$\mathcal{F}_{\theta}(\mathbf{z}^{i}_{t})=\nabla_{\mathbf{z}}A_{\mathcal{N}}(\mathbf{z}^{i}_{t};\theta_{j})$

$L^{i}\leftarrow\mathcal{L}_{\mathcal{N}}(\mathbf{z}^{i}_{t})$

$w^{i}\leftarrow\frac{\exp{(\kappa_{l}L^{i}})}{\sum_{i}\exp{(\kappa_{l}L^{i}})}$

$\mathbf{m}_{t+1}\leftarrow\beta_{1}\mathbf{m}_{t}+(1-\beta_{1})\sum_{i}\mathbf{z}^{i}_{t}w^{i}$

$\mathbf{v}_{t+1}\leftarrow\beta_{2}\mathbf{v}_{t}+(\kappa_{l}+\kappa_{h})(1-\beta_{2})\sum_{i}(\mathbf{z}^{i}_{t}-\mathbf{m}_{t})^{2}w^{i}$

$\mathbf{m}\leftarrow\frac{\mathbf{m}_{t+1}}{1-\beta_{1}^{t}}$

$\mathbf{v}\leftarrow\frac{\mathbf{v}_{t+1}}{1-\beta_{2}^{t}}$

$\mathbf{z}^{i}_{t+1}\leftarrow\mathbf{z}^{i}_{t}-\frac{\alpha}{\gamma}(\mathbf{z}_{t}^{i}-\mathbf{m})\odiv\mathbf{v}+\sqrt{\frac{2\alpha}{\gamma\kappa_{h}}}\mathbf{\eta}_{t},\eta_{t}\sim\mathcal{N}(0,1)$

$t\leftarrow t+1$

end while

Save the training dataset
$\mathcal{D}_{j}=\\{\mathbf{z}^{i}_{t},\mathbf{F}^{i}_{t}\\}_{t=0}^{T}$

Optimize $\theta_{j+1}$ using the generated training set $\mathcal{D}_{l}$ for
$l=0,\ldots,j$.

$j\leftarrow j+1$

end while

## Numerical examples

### One-dimensional Rastrigin function

To illustrate the essential idea of the present method, we start with the
1-dimensional Rastrigin function:

$f(x)=x^{2}-\cos(2\pi x),\quad x\in[-3,3].$ (15)

Figure 2: Adaptive sampling and construction of the 1D Rastrigin function.
Left: The reference function $f(x)$ and the constructed approximations
$f_{\theta_{i}}(x)$ obtained at different iteration steps. The relative
$l_{2}$ error is less than $6\times 10^{-3}$ after $12$ iterations. Right: The
residual function $|f(x)-f_{\theta_{i}}(x)|$. The symbols represent the
locations identified by each adaptive sampling (i.e., the maximization) step
where new points will be added for the next training (i.e., the minimization)
step.

Instead of a neural network, we simply use a piecewise polynomial function
$f_{\theta}(x)$ for this 1D problem, where $\theta$ represents the fitting
parameters. Accordingly, the loss function is directly defined as
$\|f-f_{\theta_{i}}\|$. Initially, we set $f_{\theta_{0}}(x)\equiv 8$ with
consistent boundary condition. We use the proposed sampling method with 10
walkers to estimate the first and second moments, yielding $m_{1}=1\times
10^{-4}$ and $V_{1}=0.022$ and therefore the first and second derivative
$f^{\prime}(m_{1})=0$ and $f^{\prime\prime}(m_{1})=\frac{1}{V_{1}}=40.18$.
This value is extremely close to the actual minimizer at $0$ (i.e., the max-
residue point) and the second derivative $f^{\prime\prime}(m_{1})=41.08$.
Accordingly, we add new data points near $x=1\times 10^{-4}$ into the training
set consisting of the boundary points $f(-3)=f(3)=8$, which enables us to
construct an updated approximation $f_{\theta_{1}}(x)$. Similarly, with the
approximation $f_{\theta_{i-1}}(x)$, we conduct another sampling process
(i.e., the maximization step) and include the obtained samples near $m_{i}$
into the training set, yielding an updated approximation $f_{\theta_{i}}(x)$.
As shown in Figure 2, for each iteration, the sampling step can pinpoint the
regime where the error is most pronounced. The underlying function $f(x)$ can
be fully recovered after $12$ iterations.

### Two-dimensional FES

Figure 3: The 2D FES for the Ala2 molecule on the $\phi-\psi$ plane
constructed by (a) Metadynamics (reference) (b) VES (c) RiD (d) CES. The
accuracy and the computational cost are shown in Table 1.

We use the alanine dipeptide (Ace-Ala-Nme), referred to as Ala2, as a
benchmark problem. The molecule is solvated in $383$ TIP3P water molecules
similar to Ref. [42]. The full MD system is simulated in a canonical ensemble
under temperature $300K$ using the Amber99-SB force field [57] with a time
step is 2 fs. We refer to the Appendix for the simulation details.

We choose the CVs as the torsion angles: $\phi$ (C, N, Cα, C) and $\psi$ (N,
Cα, C, N). For comparison, we also construct the FESs using the VES [8] and
RiD [42] method with the same setup parameters presented therein. For the CES
method, we use $10$ walkers to explore the configuration space. The initial
points of the walkers at the $k$th iteration are chosen to be the final $10$
points of a biased dynamics with the bias potential
$-A_{\mathcal{N}}(\mathbf{z};\theta_{k-1})$. In the sampling stage, the
inverse of the low and high temperatures are set to be $\kappa_{l}=10$ and
$\kappa_{h}=1$, respectively. For each sample point, restrained dynamics is
conducted with $\kappa=500$ for $5000$ steps to compute the average force. The
timestep $\alpha$ for the dynamics of the random walkers (11) is set to be
$0.1$. Figure 3 shows the FESs constructed by the three different methods and
the reference obtained by the metadynamics [4] using a long simulation time.
As presented in Table 1, the CES method yields smaller approximation error and
meanwhile requires lower computational cost. The better performance is not
unexpected since the sample points can be adaptively optimized based on the
construction error.

It is worth mentioning that the VES method was initially developed to achieve
efficient enhanced sampling in the high-dimensional phase space; it is
somewhat unfair to use it for direct FES construction since the estimation of
the probability density function from the obtained samples could be
numerically challenging for high-dimensional cases. On the other hand, the
present CES method enables the explicit construction of FES if the analytical
form is needed. Also, we note that the accuracy of the RiD method can be
further improved using more iterations at the cost of larger simulation and
training overhead. Below, we test the methods in more complex systems.

Method | Accuracy | Time
---|---|---
$l_{2}$ error |  $l_{\infty}$ error | Simulation | Train
VES | $5.39$ | $21.03$ | $47.5$
RiD | $3.15$ | $11.04$ | $17.98$ |  $0.22$ (GPU)
CES | $1.88$ | $10.68$ | $0.23\times 10$ |  $0.18$ (CPU)
$0.13$ (GPU)
Table 1: The accuracy of the constructed 2D FES (the Ala2 molecule) and
computational time (in hours, the same below) for the VES, RiD, and CES
methods. The $l_{2}$ and $l_{\infty}$ error are computed up to $40$ KJ/mol.
The reference solution is constructed by the metadynamics. The simulation time
of the CES method is multiplied by $10$ since 10 walkers are used.

### Three-dimensional FES

Next, we consider a s-(1)-phenylethyl (s1pe) peptoid in an aqueous environment
similar to Refs. [58, 52]. The full system consists of one biomolecule and 546
water molecules in a $(2.9{\rm nm})^{3}$ dodecahedron box. The CHARMM general
force field (CGenFF) [58] is used for the biomolecule and the TIP3P model [59]
is used for the water molecules. The system temperature is set to be 298 K
with a time step is $2$ fs. We refer to the Appendix for details.

The CVs are the three torsion angles $\omega$ (Cα, C, N, Cα), $\phi$, and
$\psi$, where the latter two are the same as the Ala2 molecule. The FES is
constructed by both the RiD and CES methods. The setup of the RiD method is
the same as [52]. For the CES method, we use 20 walkers and set the inverse of
the low and high temperatures to be $\kappa_{l}=10$ and $\kappa_{h}=2$,
respectively. The initial points of these walkers at each iteration are chosen
in the same method as the two-dimensional problem. The timestep $\alpha$ for
the dynamics of the random walkers (11) is set to be $0.1$. For each sample
point, restrained dynamics with $\kappa=500$ is conducted for $10000$ steps to
compute the mean force.

For visualization, we project the constructed FES onto a two-dimensional plane
and fix the third variable. Figure 4 shows the projected FES on the
$\omega-\phi$ and $\phi-\psi$ plane obtained from the CES and the RiD methods.
For each projection, the reference is constructed as a 2D FES using the
metadynamics [4]. Similar to the previous 2D case, the present CES method
yields higher accuracy with lower computational cost.

Figure 4: The 3D FES for molecule s1pe projected on the $\omega-\phi$ with $\psi=1.5$ (left column) and $\phi-\psi$ planes with $\omega=1.5$ (right-column). (a-b) The 2D FES constructed using metadynamics with the third variable restrained (reference); (c-d) projection of the 3D FES constructed by the RiD method; (e-f) projection of the 3D FES constructed by the CES method.  Method | Restraint | Accuracy | Time
---|---|---|---
$l_{2}$ error |  $l_{\infty}$ error | Simulation | Train
RiD | $\psi=1.5$ | $5.76$ | $25.72$ | $423.33$ | $8$ (GPU)
$\omega=1.5$ | $12.04$ | $49.13$
CES | $\psi=1.5$ | $2.44$ | $11.21$ | $4.81\times 20$ | $0.84$ (GPU)
$\omega=1.5$ | $3.89$ | $28.80$
Table 2: The accuracy of the constructed 3D FES (the s1pe molecule) and
computational time for RiD and CES methods. The $l_{2}$ and $l_{\infty}$ error
are calculated up to $40$ KJ/mol. For each case, the FES is projected onto a
2D plane with the third variable fixed; the reference solution is constructed
as a 2D FES by the metadynamics. The simulation time of the CES method is
multiplied by $20$ since 20 walkers are used.

### Nine-dimensional FES

Furthermore, we consider a more complex molecule, the peptoid trimer (s1pe)3,
solvated in a (4.2 nm)3 dodecahedron box with 1622 TIP3P water molecules. The
force field and other simulation setups are similar to the s1pe molecule.

The chosen CVs consist of the $9$ torsion angles $\omega,\phi,\psi$ that are
defined in Ala2 and s1pe case, associated with the different Cα atoms and
denoted as
$\omega_{1},\phi_{1},\psi_{1},\omega_{2},\phi_{2},\psi_{2},\omega_{3},\phi_{3},\psi_{3}$.
We use 64 walkers for this case and the initial conditions of the walkers at
each iteration are chosen with the same method as the previous cases. The
inverse of the low and the high temperature are set to be $\kappa_{l}=100$ and
$\kappa_{h}=2$, respectively. The timestep $\alpha$ of the dynamics of the
random walkers is set to be $0.1$. The FES is constructed by the CSE method
using 28 iterations of sampling and training, which requires $225.53\times
64=14434.35$ CPU hours for simulation and $6.06$ GPU hours for training. For
comparison, the RiD method uses $17900$ CPU hours for simulation and $15.44$
GPU hours for training. Similar to the above 3D problem, the constructed 9D
FES is projected onto various 2D planes with the remaining variables fixed.
For each 2D projection, the reference is constructed as a 2D FES using
metadynamics. Fig. 5 shows the projection on the $\omega_{1}-\phi_{1}$ and
$\omega_{1}-\psi_{1}$ plane (see Appendix for visualization of other 2D
projections). Compared with the Rid method, the present CES method yields
higher accuracy with lower computational cost.

The different performances could be possibly due to the distinct ways of
imposing the sampling adaptivity for the two methods. Specifically, the Rid
method trains a replica of DNNs on the same sample set and uses the standard
deviation of multiple DNNs’ predictions as an indirect measure of the
construction error. The significant sampling error of the mean force (i.e.,
$-\nabla A(\mathbf{z})$) could be overlooked by using the standard deviation
as the uncertainty indicator, and further propagate into the construction
error. For instance, a biased sampling of mean force may lead to a consistent
biased prediction among multiple DNNs. As a result, even if the standard
deviation of multiple DNNs’ predictions for a sample point is small, the
construction error could be still pronounced in that region. Instead, the CES
method directly uses the construction error to impose the adaptivity and
therefore achieve enhanced sampling of the regions lacking accuracy.

Figure 5: The 9D FES for molecule (s1pe)3 projected on the
$\omega_{1}-\phi_{1}$ (left column) and $\omega_{1}-\psi_{1}$ plane (right
column). (a-b) The 2D FES constructed using metadynamics with the remaining
variable restrained (reference); (c-d) projection of the 9D FES constructed by
the RiD method using $17900$ CPU hours for simulation and $15.44$ GPU hours
for training; (e-f) projection of the 9D FES constructed by the CES method
using $14434.35$ CPU hours for simulation and $6.06$ GPU hours for training.

### Thirty dimensional FES

Finally, we consider polyalanine-15 (Ace-(Ala)15-Nme), denoted as Ala16,
solvated in $2258$ water molecules. The full system is simulated temperature
300K in a ($4.62$ nm)3 dodecahedron box using the Amber99-SB force field with
a time step is 2 fs. We refer to the Appendix for the simulation details.

The chosen CVs consist of torsion angles $\phi$ and $\psi$, defined in Ala2
case, associated with the different Cαs, denoted as
$\\{\phi_{i},\psi_{i}\\}_{i=1}^{15}$. We use 64 walkers for this case and the
initial value of the walkers at each iteration is given by the biased
simulation before. The inverse of the low and high temperatures are set to be
$\kappa_{l}=20$ and $\kappa_{h}=5$, respectively. The timestep $\alpha$ for
dynamics of the random walkers 11 is set to be $0.1$. The FES is constructed
using $100$ iterations of sampling and training. Similar to the previous case,
the obtained FES is plotted on a two-dimensional plane while the remaining
variables are fixed and the reference is constructed as a 2D FES using the
metadynamics. Fig. 6 shows the projection on the $\phi_{1}-\psi_{1}$,
$\phi_{2}-\phi_{3}$ and $\psi_{2}-\psi_{3}$ plane. For all the cases, the
projection of the 30-dimensional FES shows good agreement with the
2-dimensional reference solution. We have also examined the projection on
other planes; the prediction shows good agreement with the reference solution
as well. We refer to the Appendix for details.

Figure 6: The 30D FES for molecule Ala16 projected on various 2D planes. (a-b)
$\phi_{1}-\psi_{1}$ (c-d) $\phi_{2}-\phi_{3}$ (e-f) $\psi_{2}-\psi_{3}$.
(a-c-e) 2D FES constructed by metadynamics (reference) (b-d-f) The 30D FES
constructed by the present CES method.

## Conclusion

We have presented a consensus-based approach to construct high-dimensional
FESs by reformulating the construction task as a minimax optimization problem.
Rather than seeking the direct fitting, this method essentially establishes an
adversarial learning of FESs by simultaneously optimizing the target function
approximation and the training set. While the common approaches mainly focus
on the efficient exploration of the phase space in the presence of local
minima, the present method further accounts for the discretization error that
has been broadly overlooked in FES construction. Adaptive sampling of the max-
residue regime is achieved through the consensus-based sampling of a
posteriori residue-induced distribution in the form of a stochastic particle
system in the CV space. Given the fact the sampling only relies on the first
and second-moment estimation, the method could be particularly efficient for
high-dimensional problems. While the numerical results of biomolecular systems
have demonstrated the effectiveness for FES construction, the present
framework of unifying the residual minimization and max-residue enhanced
sampling is quite general for model reduction of complex systems, e.g., the
stochastic reduced dynamics (e.g., see Ref. [60]) with a state-dependent
memory for multiscale physical systems.

##  Appendix A Proof of Proposition 2

The sampling distribution is governed by the following McKean stochastic
differential equation

$\dot{\mathbf{z}}^{i}_{t}=-\frac{1}{\gamma}\nabla_{\mathbf{z}}G(\mathbf{z}^{i}_{t})+\sqrt{\frac{2}{\kappa_{h}\gamma}}\mathbf{\xi}_{i}(t),$
(16)

where
$G(\mathbf{z}_{t})=\frac{1}{2}(\mathbf{z}_{t}-\mathbf{m}_{t})^{T}V_{t}^{-1}(\mathbf{z}_{t}-\mathbf{m}_{t})$
and $\gamma$ is the friction coefficient. We want to show it converges to a
steady state as an approximation of $q^{\ast}$ with the proper choices of
$\mathbf{m}_{t}$ and $V_{t}$ given by

$\displaystyle\mathbf{m}_{t}$
$\displaystyle=\sum_{i=1}^{N_{w}}\mathbf{z}^{i}_{t}\hat{p}(\mathbf{z}^{i}_{t}),$
(17) $\displaystyle V_{t}$
$\displaystyle=\kappa_{t}\sum_{i=1}^{N_{w}}(\mathbf{z}^{i}_{t}-\mathbf{m}_{t})(\mathbf{z}^{i}_{t}-\mathbf{m}_{t})^{T}\hat{p}(\mathbf{z}^{i}_{t}),$

where
$\hat{p}(\mathbf{z})=\frac{\exp{(-\kappa_{l}\mathcal{L}_{\mathcal{N}}^{-}(\mathbf{z}))}}{\sum_{i=1}^{N_{w}}\exp{(-\kappa_{l}\mathcal{L}_{\mathcal{N}}^{-}(\mathbf{z}^{i})})}$
and $\kappa_{t}=\kappa_{l}+\kappa_{h}$.

Proposition 2. Suppose $\mathcal{L}_{\mathcal{N}}^{-}(\mathbf{z})$ takes a
local quadratic approximation in form of
$\frac{1}{2}(\mathbf{z}-\boldsymbol{\mu})^{T}\Sigma^{-1}(\mathbf{z}-\boldsymbol{\mu})$,
$q_{t}\to\frac{\exp{(-\kappa_{h}\mathcal{L}_{\mathcal{N}}^{-}(\mathbf{z}))}}{\int\exp{(-\kappa_{h}\mathcal{L}_{\mathcal{N}}^{-}(\mathbf{z}))}\mathrm{d}\mathbf{z}}$
as $t\to\infty$, given $\kappa_{t}=\kappa_{l}+\kappa_{h}$.

###### Proof.

Let $q_{\infty}(\mathbf{z})$ denote the invariant distribution of the McKean
SDE (16). Then $q_{\infty}(\mathbf{z})$ must be the invariant distribution of
the following McKean stochastic dynamics

$\dot{\mathbf{z}}=-\frac{1}{\gamma}V^{-1}_{\kappa_{l},\infty}(\mathbf{z}-\mathbf{m}_{\kappa_{l},\infty})+\sqrt{\frac{2}{\gamma\kappa_{h}}}\xi(t),$
(18)

where $\mathbf{m}_{\kappa_{l},\infty}$ and
$\kappa_{t}^{-1}V_{\kappa_{l},\infty}$ are the mean and the covariance matrix
of the re-weighted density
$\propto~{}q_{\infty}(\mathbf{z})e^{-\kappa_{l}\mathcal{L}_{\mathcal{N}}^{-}(\mathbf{z})}$.
With the fluctuation-dissipation relation for Eq. (18), we can show
$q_{\infty}(\mathbf{z})$ follows the Gaussian distribution with mean
$\mathbf{m}_{\kappa_{l},\infty}$ and covariance matrix
$\kappa_{h}^{-1}V_{\kappa_{l},\infty}$.

Since
$\mathcal{L}_{\mathcal{N}}(\mathbf{z})=\frac{1}{2}(\mathbf{z}-\boldsymbol{\mu})^{T}\Sigma^{-1}(\mathbf{z}-\boldsymbol{\mu})$
is quadratic, the re-weighted density of a Gaussian distribution
$q(\mathbf{z})\sim\mathcal{N}(\mathbf{m},V)$ remains Gaussian, i.e.,
$q(\mathbf{z})e^{-\kappa_{l}\mathcal{L}_{\mathcal{N}}^{-}(\mathbf{z})}\propto~{}\mathcal{N}({\mathbf{m}}_{\kappa_{l}},{V}_{\kappa_{l}})$,
where ${\mathbf{m}}_{\kappa_{l}}$ and ${V}_{\kappa_{l}}$ are defined by

$\begin{split}\mathbf{m}_{\kappa_{l}}&=(V^{-1}+\kappa_{l}\Sigma^{-1})^{-1}(\kappa_{l}\Sigma^{-1}\boldsymbol{\mu}+V^{-1}\mathbf{m}),\\\
V_{\kappa_{l}}&=(V^{-1}+\kappa_{l}\Sigma^{-1})^{-1}.\end{split}$ (19)

In particular, we choose $\mathbf{m}=\mathbf{m}_{\kappa_{l},\infty}$ and
$V=\kappa_{h}^{-1}V_{\kappa_{l},\infty}$, then Eq. (19) yields

$\begin{split}\mathbf{m}_{\kappa_{l},\infty}&=(\kappa_{h}V_{\kappa_{l},\infty}^{-1}+\kappa_{l}\Sigma^{-1})^{-1}(\kappa_{l}\Sigma^{-1}\boldsymbol{\mu}+\kappa_{h}V_{\kappa_{l},\infty}^{-1}\mathbf{m}_{\kappa_{l},\infty}),\\\
\kappa_{t}^{-1}V_{\kappa_{l},\infty}&=(\kappa_{h}V_{\kappa_{l},\infty}^{-1}+\kappa_{l}\Sigma^{-1})^{-1}.\end{split}$

It is easy to show that by choosing $\kappa_{t}=\kappa_{l}+\kappa_{h}$,
$\mathbf{m}_{\kappa_{l},\infty}$ and $V_{\kappa_{l},\infty}$ recovers
$\boldsymbol{\mu}$ and $\Sigma$, respectively, and the invariant density takes
the form

$q_{\infty}(\mathbf{z})\sim\mathcal{N}\left(\boldsymbol{\mu},\kappa_{h}^{-1}\Sigma\right).$

∎

##  Appendix B Simulation setup

All the MD simulations are performed using the package GROMACS 2019.2 [61].
The simulation is carried out on Intel(R) Xeon(R) Platinum 8260 CPU.

###  B.1 ala2 (the two-dimensional problem)

The Ace-Ala-Nme (ala2) molecule is modeled by the Amber99SB force field [57].
The molecules are dissolved in 383 water molecules in a periodic simulation
cell. The cut-off radius of the van der Waals interaction is 0.9 nm. The
Coulomb interaction is treated with the smooth particle mesh Ewald method with
a real space cutoff of 0.9 nm and a reciprocal space grid spacing of 0.12 nm.
The system is integrated with the leap-frog scheme at time step 2 fs. The
temperature of the system is set to 300 K by a velocity-rescale thermostat
[62] with a relaxation time of 0.2 ps. The Parrinello-Rahman barostat [63]
with a relaxation time scale of 1.5 ps and a compressibility of $4.5\times
10^{-5}$ bar-1 is coupled to the system to control the pressure to 1 bar. The
hydrogen atom is constrained by the LINCS algorithm[64] and the H–O bond and
H–O–H angle of water molecules are constrained by the SETTLE algorithm[65].

###  B.2 s1pe (three-dimensional problem)

The s-(1)-phenylethyl (s1pe) molecule is modeled by the CHARMM general force
field (CGenFF) [58]. The molecule is dissolved in 546 TIP3P water molecules in
a (2.69 nm)3 dodecahedron box. The cut-off radius of the van der Waals
interaction is 1 nm. The Coulomb interaction is treated with the smooth
particle mesh Ewald method with a real space cutoff of 1 nm and a reciprocal
space grid spacing of 0.12 nm. The system is integrated with the leap-frog
scheme at time step 2 fs. The temperature of the system is set to 298 K by a
velocity-rescale thermostat [62] with a relaxation time of 0.2 ps. The
Parrinello-Rahman barostat pressure couple, the hydrogen atom constraint, the
H–O bond constraint, and the H–O–H angle of water molecules constraint are the
same as the previous one.

###  B.3 (s1pe)3 (the nine-dimensional problem)

The (s1pe)3 molecule is modeled by the CHARMM general force field (CGenFF)
[58]. The molecule is solvated in a (4.2 nm)3 dodecahedron box with 1622 TIP3P
water molecules. Other simulation setups are similar to the s1pe molecule.

###  B.4 ala16 (the thirty-dimensional problem)

The Ace-(Ala)15-Nme (ala16) molecule is modeled by the Amber99SB force field
[57]. The molecule is solvated in 2258 water molecules in a (4.62 nm)3
dodecahedron box. Other simulation setups are similar to the ala2 molecule.

##  Appendix C Traning

The training data is collected during the sampling process by restricted
dynamics. The initial $10\%$ steps of the restricted dynamics are used as
equilibrium and the rest $90\%$ steps are used to calculate the mean force.
The FES are parameterized as a fully connected neural network. The depth and
width of the NN are shown in Table 3. The NNs are trained by Adam for $100000$
steps with a learning rate $1\times 10^{-3}$. For each training step, $5000$
sampling points are randomly selected from the data set. All the training
process is carried out on v100 with 32768 MB memory.

molecule | depth | width
---|---|---
ala2 | 3 | 48
s1pe | 4 | 64
(s1pe)3 | 4 | 512
ala16 | 4 | 640
Table 3: The depth and width of the NN used to parameterize the FES of
different molecules.

##  Appendix D Additional result

###  D.1 Additional result of the three-dimensional FES (s1pe)

Figure 7-9 shows the additional 2D projections of the 3D FES (molecule s1pe)
constructed by the Rid and presented CES method. For each case, the reference
is constructed as a 2D FES using the metadynamics. The CES method yields a
better agreement with the reference.

Figure 7: The 3D FES for the s1pe molecule projected on the $\phi-\psi$ plane
with the third variable $\omega=-1$ (first row) and $\omega=-2$ (second row).
(a) 2D FES by metadynamics (reference) (b) RiD (c) CES.

Figure 8: The 3D FES of the s1pe molecule on the $\omega-\psi$ plane with the
third variable $\phi=0$ (first row) and $\phi=-1$ (second row). (a) 2D FES by
metadynamics (reference) (b) RiD (c) CES.

Figure 9: The 3D FES of the s1pe molecule on the $\omega-\phi$ plane with the
third variable $\psi=0$ (first row) and $\psi=-1$ (second row). (a) 2D FES by
metadynamics (reference) (b) RiD (c) CES.

###  D.2 Additional results of the nine-dimensional FES (s1pe)3

The addition 2D projections of the 9D FES for molecule (s1pe)3 are presented
in Figure 10. Similar to the 3D case, the FES constructed by the present CES
method shows a better agreement with the reference constructed as a 2D FES
using metadynamics.

Figure 10: The 9D FES for the (s1pe)3 molecule projected on the
$\omega_{1}-\phi_{1},\omega_{2}-\phi_{2},\omega_{2}-\psi_{2},\phi_{2}-\psi_{2},\omega_{3}-\phi_{3},\omega_{3}-\psi_{3},\phi_{3}-\psi_{3}$
plane from top to bottom. (a) 2D FES by metadynamics (reference) (b) RiD (c)
CES.

###  D.3 Additional results of the thirty-dimensional FES (ala16)

The addition 2D projections of the 30D FES for molecule ala16 are presented in
Figure 11. We note that the projection on the $\phi_{5}-\psi_{5}$ plane is
significantly different from other projections such as the $\phi_{1}-\psi_{1}$
plane in the main context. While the ala16 molecule consists of 15 sequential
alanine residues, the FES for individual $\phi-\psi$ projections shows
different features. The numerical results of the present CES method show good
agreement with the reference.

Figure 11: The 30D FES of the ala16 projected on the $\phi_{5}-\psi_{5}$
plane. (a) 2D FES constructed by metadynamics (reference) (b) 30D FES
constructed by the CES method.

##  Appendix E Code availablity

Python implementation of our codes are available at github
https://github.com/Lyuliyao/consensus-sampling-method-for-expolering-high-
dimensional-energy-surface.

## References

  * [1] Daan Frenkel and Berend Smit.  Understanding molecular simulation: from algorithms to applications.  1996\.
  * [2] Glenn M Torrie and John P Valleau.  Nonphysical sampling distributions in monte carlo free-energy estimation: Umbrella sampling.  Journal of Computational Physics, 23(2):187–199, 1977.
  * [3] Shankar Kumar, John M. Rosenberg, Djamal Bouzida, Robert H. Swendsen, and Peter A. Kollman.  The weighted histogram analysis method for free-energy calculations on biomolecules. i. the method.  Journal of Computational Chemistry, 13(8):1011–1021, 1992.
  * [4] Alessandro Laio and Michele Parrinello.  Escaping free-energy minima.  Proceedings of the National Academy of Sciences, 99(20):12562–12566, 2002.
  * [5] Alessandro Barducci, Giovanni Bussi, and Michele Parrinello.  Well-tempered metadynamics: a smoothly converging and tunable free-energy method.  Physical Review Letters, 100(2):020603, 2008.
  * [6] Omar Valsson and Michele Parrinello.  Variational approach to enhanced sampling and free energy calculations.  Phys. Rev. Lett., 113:090601, Aug 2014.
  * [7] Patrick Shaffer, Omar Valsson, and Michele Parrinello.  Enhanced, targeted sampling of high-dimensional free-energy landscapes using variationally enhanced sampling, with an application to chignolin.  Proceedings of the National Academy of Sciences, 113(5):1150–1155, 2016.
  * [8] Luigi Bonati, Yue-Yu Zhang, and Michele Parrinello.  Neural networks-based variationally enhanced sampling.  Proceedings of the National Academy of Sciences, 116(36):17641–17647, 2019.
  * [9] Bin W Zhang, David Jasnow, and Daniel M Zuckerman.  Efficient and verified simulation of a path ensemble for conformational change in a united-residue model of calmodulin.  Proceedings of the National Academy of Sciences, 104(46):18043–18048, 2007.
  * [10] Surl-Hee Ahn, Jay W Grate, and Eric F Darve.  Efficiently sampling conformations and pathways using the concurrent adaptive sampling (cas) algorithm.  The Journal of Chemical Physics, 147(7), 2017.
  * [11] Surl-Hee Ahn, Anupam A Ojha, Rommie E Amaro, and J Andrew McCammon.  Gaussian-accelerated molecular dynamics with the weighted ensemble method: A hybrid method improves thermodynamic and kinetic sampling.  Journal of Chemical Theory and Computation, 17(12):7938–7951, 2021\.
  * [12] Arthur F Voter.  Hyperdynamics: Accelerated molecular dynamics of infrequent events.  Physical Review Letters, 78(20):3908, 1997.
  * [13] Lianqing Zheng, Mengen Chen, and Wei Yang.  Random walk in orthogonal space to achieve efficient free-energy simulation of complex systems.  Proceedings of the National Academy of Sciences, 105(51):20227–20232, 2008.
  * [14] Eric Darve and Andrew Pohorille.  Calculating free energies using average force.  The Journal of Chemical Physics, 115(20):9169–9183, 2001.
  * [15] Eric Darve, David Rodríguez-Gómez, and Andrew Pohorille.  Adaptive biasing force method for scalar and vector free energy calculations.  The Journal of Chemical Physics, 128(14), 2008.
  * [16] Tony Lelièvre, Mathias Rousset, and Gabriel Stoltz.  Long-time convergence of an adaptive biasing force method.  Nonlinearity, 21(6):1155, 2008.
  * [17] Stefano Piana and Alessandro Laio.  A bias-exchange approach to protein folding.  The journal of Physical Chemistry B, 111(17):4553–4559, 2007.
  * [18] Jim Pfaendtner and Massimiliano Bonomi.  Efficient sampling of high-dimensional free-energy landscapes with parallel bias metadynamics.  Journal of Chemical Theory and Computation, 11(11):5062–5067, 2015\.
  * [19] Arushi Prakash, Christopher D Fu, Massimiliano Bonomi, and Jim Pfaendtner.  Biasing smarter, not harder, by partitioning collective variables into families in parallel bias metadynamics.  Journal of Chemical Theory and Computation, 14(10):4985–4990, 2018\.
  * [20] Yuji Sugita and Yuko Okamoto.  Replica-exchange molecular dynamics method for protein folding.  Chemical Physics Letters, 314(1-2):141–151, 1999.
  * [21] Enzo Marinari and Giorgio Parisi.  Simulated tempering: a new monte carlo scheme.  Europhysics letters, 19(6):451, 1992.
  * [22] Yi Qin Gao.  An integrate-over-temperature approach for enhanced sampling.  The Journal of Chemical Physics, 128(6), 2008.
  * [23] Lijiang Yang and Yi Qin Gao.  A selective integrated tempering method.  The Journal of Chemical Physics, 131(21), 2009.
  * [24] Jaegil Kim, John E Straub, and Thomas Keyes.  Statistical-temperature monte carlo and molecular dynamics algorithms.  Physical review letters, 97(5):050601, 2006.
  * [25] Donald Hamelberg, John Mongan, and J Andrew McCammon.  Accelerated molecular dynamics: a promising and efficient simulation method for biomolecules.  The Journal of Chemical Physics, 120(24):11919–11929, 2004.
  * [26] Yinglong Miao, Victoria A Feher, and J Andrew McCammon.  Gaussian accelerated molecular dynamics: Unconstrained enhanced sampling and free energy calculation.  Journal of Chemical Theory and Computation, 11(8):3584–3595, 2015\.
  * [27] Jianfeng Lu and Eric Vanden-Eijnden.  Infinite swapping replica exchange molecular dynamics leads to a simple simulation patch using mixture potentials.  The Journal of Chemical Physics, 138(8), 2013.
  * [28] Tang-Qing Yu, Jianfeng Lu, Cameron F Abrams, and Eric Vanden-Eijnden.  Multiscale implementation of infinite-swap replica exchange molecular dynamics.  Proceedings of the National Academy of Sciences, 113(42):11744–11749, 2016.
  * [29] Yi Isaac Yang, Qiang Shao, Jun Zhang, Lijiang Yang, and Yi Qin Gao.  Enhanced sampling in molecular dynamics.  The Journal of Chemical Physics, 151(7), 2019.
  * [30] Luca Maragliano and Eric Vanden-Eijnden.  A temperature accelerated method for sampling free energy and determining reaction pathways in rare events simulations.  Chemical physics letters, 426(1-3):168–175, 2006.
  * [31] Lula Rosso, Peter Mináry, Zhongwei Zhu, and Mark E Tuckerman.  On the use of the adiabatic molecular dynamics technique in the calculation of free energy profiles.  The Journal of Chemical Physics, 116(11):4389–4402, 2002.
  * [32] Jerry B Abrams and Mark E Tuckerman.  Efficient and direct generation of multidimensional free energy surfaces via adiabatic dynamics without coordinate transformations.  The Journal of Physical Chemistry B, 112(49):15742–15757, 2008\.
  * [33] Cameron F Abrams and Eric Vanden-Eijnden.  Large-scale conformational sampling of proteins using temperature-accelerated molecular dynamics.  Proceedings of the National Academy of Sciences, 107(11):4961–4966, 2010.
  * [34] Luca Maragliano and Eric Vanden-Eijnden.  Single-sweep methods for free energy calculations.  The Journal of Chemical Physics, 128(18):184110, 05 2008.
  * [35] Thomas Stecher, Noam Bernstein, and Gábor Csányi.  Free energy surface reconstruction from umbrella samples using gaussian process regression.  Journal of Chemical Theory and Computation, 10(9):4079–4097, 2014\.
  * [36] Letif Mones, Noam Bernstein, and Gábor Csányi.  Exploration, sampling, and reconstruction of free energy surfaces with gaussian process regression.  Journal of Chemical Theory and Computation, 12(10):5100–5110, 2016\.
  * [37] Elia Schneider, Luke Dai, Robert Q Topper, Christof Drechsel-Grau, and Mark E Tuckerman.  Stochastic neural network approach for learning high-dimensional free energy surfaces.  Physical Review Letters, 119(15):150601, 2017.
  * [38] Hythem Sidky and Jonathan K Whitmer.  Learning free energy landscapes using artificial neural networks.  The Journal of Chemical Physics, 148(10):104111, 2018.
  * [39] Joseph R. Cendagorta, Jocelyn Tolpin, Elia Schneider, Robert Q. Topper, and Mark E. Tuckerman.  Comparison of the performance of machine learning models in representing high-dimensional free energy surfaces and generating observables.  The Journal of Physical Chemistry B, 124(18):3647–3660, 2020.
  * [40] Frank Noé, Simon Olsson, Jonas Köhler, and Hao Wu.  Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning.  Science, 365(6457):eaaw1147, 2019.
  * [41] Marylou Gabrié, Grant M Rotskoff, and Eric Vanden-Eijnden.  Adaptive monte carlo augmented with normalizing flows.  Proceedings of the National Academy of Sciences, 119(10):e2109420119, 2022.
  * [42] Linfeng Zhang, Han Wang, and Weinan E.  Reinforced dynamics for enhanced sampling in large atomic and molecular systems.  The Journal of Chemical Physics, 148(12):124113, 2018.
  * [43] Dongdong Wang, Yanze Wang, Junhan Chang, Linfeng Zhang, Han Wang, and Weinan E.  Efficient sampling of high-dimensional free energy landscapes using adaptive reinforced dynamics.  Nature Computational Science, 2(1):20–29, 2022.
  * [44] Kejun Tang, Jiayu Zhai, Xiaoliang Wan, and Chao Yang.  Adversarial adaptive sampling: Unify pinn and optimal transport for the approximation of pdes.  arXiv preprint arXiv:2305.18702, 2023.
  * [45] Gang Bao, Xiaojing Ye, Yaohua Zang, and Haomin Zhou.  Numerical solution of inverse problems by weak adversarial networks.  Inverse Problems, 36(11):115003, 2020.
  * [46] Jose Antonio Carrillo, Franca Hoffmann, Andrew M Stuart, and Urbain Vaes.  Consensus-based sampling.  Studies in Applied Mathematics, 148(3):1069–1140, 2022.
  * [47] René Pinnau, Claudia Totzeck, Oliver Tse, and Stephan Martin.  A consensus-based model for global optimization and its mean-field limit.  Mathematical Models and Methods in Applied Sciences, 27(01):183–204, 2017.
  * [48] José A Carrillo, Young-Pil Choi, Claudia Totzeck, and Oliver Tse.  An analytical framework for consensus-based global optimization method.  Mathematical Models and Methods in Applied Sciences, 28(06):1037–1066, 2018.
  * [49] José A Carrillo, Shi Jin, Lei Li, and Yuhua Zhu.  A consensus-based global optimization method for high dimensional machine learning problems.  ESAIM: Control, Optimisation and Calculus of Variations, 27:S5, 2021\.
  * [50] Jingrun Chen, Shi Jin, and Liyao Lyu.  A consensus-based global optimization method with adaptive momentum estimation.  Communications in Computational Physics, 31(4):1296–1316, 2022\.
  * [51] Michael P Allen and Dominic J Tildesley.  Computer simulation of liquids.  Oxford University Press, 2017.
  * [52] Haoran Wang, Thaleia Zariphopoulou, and Xun Yu Zhou.  Reinforcement learning in continuous time and space: A stochastic control approach.  The Journal of Machine Learning Research, 21(1):8145–8178, 2020\.
  * [53] Xuefeng Gao, Zuo Quan Xu, and Xun Yu Zhou.  State-dependent temperature control for langevin diffusions.  SIAM Journal on Control and Optimization, 60(3):1250–1268, 2022\.
  * [54] Ishaan Gulrajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron C Courville.  Improved training of wasserstein gans.  Advances in neural information processing systems, 30, 2017.
  * [55] Takeru Miyato, Toshiki Kataoka, Masanori Koyama, and Yuichi Yoshida.  Spectral normalization for generative adversarial networks.  arXiv preprint arXiv:1802.05957, 2018.
  * [56] Diederik Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  International Conference on Learning Representations (ICLR), 12 2015\.
  * [57] Viktor Hornak, Robert Abel, Asim Okur, Bentley Strockbine, Adrian Roitberg, and Carlos Simmerling.  Comparison of multiple amber force fields and development of improved protein backbone parameters.  Proteins: Structure, Function, and Bioinformatics, 65(3):712–725, 2006.
  * [58] Laura J Weiser and Erik E Santiso.  A cgenff-based force field for simulations of peptoids with both cis and trans peptide bonds.  Journal of computational chemistry, 40(22):1946–1956, 2019.
  * [59] William L Jorgensen, Jayaraman Chandrasekhar, Jeffry D Madura, Roger W Impey, and Michael L Klein.  Comparison of simple potential functions for simulating liquid water.  The Journal of Chemical Physics, 79(2):926–935, 1983.
  * [60] Zhiyuan She, Pei Ge, and Huan Lei.  Data-driven construction of stochastic reduced dynamics encoded with non-markovian features.  The Journal of Chemical Physics, 158(3):034102, 2023.
  * [61] Lindahl, Abraham, Hess, and van der Spoel.  Gromacs 2019.2 source code, April 2019.
  * [62] Giovanni Bussi, Davide Donadio, and Michele Parrinello.  Canonical sampling through velocity rescaling.  The Journal of Chemical Physics, 126(1), 2007.
  * [63] Michele Parrinello and Aneesur Rahman.  Polymorphic transitions in single crystals: A new molecular dynamics method.  Journal of Applied Physics, 52(12):7182–7190, 1981.
  * [64] Berk Hess, Henk Bekker, Herman JC Berendsen, and Johannes GEM Fraaije.  Lincs: A linear constraint solver for molecular simulations.  Journal of Computational Chemistry, 18(12):1463–1472, 1997.
  * [65] Shuichi Miyamoto and Peter A Kollman.  Settle: An analytical version of the shake and rattle algorithm for rigid water models.  Journal of Computational Chemistry, 13(8):952–962, 1992.