The present paper proposes an adaptive biasing potential for the computation of free energy landscapes. It is motivated by statistical learning arguments and unifies the tasks of biasing the molecular dynamics to escape free energy wells and estimating the free energy function, under the same objective. It offers rigorous convergence diagnostics even though history dependent, non-Markovian dynamics are employed. It makes use of a greedy optimization scheme in order to obtain sparse representations of the free energy function which can be particularly useful in multidimensional cases. It employs embarrassingly parallelizable sampling schemes that are based on adaptive Sequential Monte Carlo and can be readily coupled with legacy molecular dynamics simulators. The sequential nature of the learning and sampling scheme enables the efficient calculation of free energy functions parametrized by the temperature. The characteristics and capabilities of the proposed method are demonstrated in three numerical examples.
Free energy is a central concept in thermodynamics and in the study of several systems in biology, chemistry and physics [8]. It represents a rigorous way to coarse-grain systems consisting of very large numbers of atomistic degrees of freedom, to probe states not accessible experimentally, to characterize global changes as well as investigate relative stabilities. In most applications, a bruteforce computation based on sampling the atomistic positions is impractical or infeasible as the free energy barriers to overcome are so large that the system remains trapped in metastable free energy sets [40,42,8,56].
Equilibrium techniques for computing free energy surfaces such as Thermodynamic Integration [29] or Adaptive Integration [49,21] require the simulation of very long atomistic trajectories in order to achieve equilibrium and lack convergence diagnostics. Techniques based on non-equilibrium path sampling [26,27,23,25] lack adaptivity and require the user to specify a particular path on the reaction coordinate space connecting two energetically important free energy regions, which can be non-trivial a task [20]. Furthermore, sampling along these paths correctly might necessitate advanced and quite involved techniques [10]. More recently proposed adaptive biasing potential [2,55,32,1,18] and adaptive biasing force [11,12,41,53,24] techniques are capable of dynamically utilizing information obtained from the atomistic trajectories to bias the current dynamics in order to facilitate the escape from metastable sets [35]. They are able to automatically discover important regions of the reaction coordinate space. Since they rely on history-dependent, non-Markovian dynamics, it is not a priori clear, and in which sense, the system reaches a stationary state, although some work has been done along theses lines in [4] for Langevin-type systems and [39,35].
We propose an adaptive biasing potential technique where the two tasks of biasing the dynamics and estimating the free energy landscape are unified under the same objective of minimizing the Kullback-Leibler divergence between appropriately selected distributions on the extended space that includes atomic coordinates and the collective variables [37,38]. This framework provides a natural way for selecting the basis functions used in the approximation of the free energy and obtaining sparse representations which is critical when multi-dimensional collective variables are used. It allows the analyst to utilize and correct any prior information on the free energy landscape and provides an efficient manner of obtaining good estimates at various temperatures. The scheme proposed is embarrassingly parallelizable and relies on adaptive Sequential Monte Carlo procedures which enable efficient sampling from the high-dimensional and potentially multi-modal distributions of interest.
2 Methodology -A statistical learning approach for adaptively calculating free energies For clarity of the presentation, we will first introduce our method for the so-called alchemical case and generalize it later for the reaction coordinate case. Consider a molecular system with generalized coordinates q ∈ M ⊂ R N following a Boltzmann-like distribution which in turn depends on some parameters z ∈ D ⊂ R d p(q|z) ∝ exp (-βV (q; z))
where V (q; z) is the potential energy of the system and β plays the role of inverse temperature. The free energy A(z) is defined, up to an additive constant, by: A(z) = -β -1 exp (-βV (q; z)) dq
Our goal is to compute the function A(z) over the whole domain D.
Let Â(z; θ) be an estimate of A(z) parametrized by θ ∈ Θ ⊂ R K . We adopt a statistical perspective of learning A(z) from simulation data. A popular approach to carrying out regression tasks and functional approximations relies on kernel models [28]. Kernel regression models have proven successful in highdimensional scenaria where d is in the order of 10 or 100 [51,52]. The unknown function is selected from a Reproducing Kernel Hilbert Space (RKHS) H K induced by a semi-positive definite kernel K(•, •). We adopt representations with respect to a kernel function K(•, •) [28]:
where z j are points in D which are selected as described in the sequence,
In order to fix the additive constant, we select a point z 0 ∈ D such that:
In relevant literature different types of kernel functions have been used such as thin plate splines, multiquadrics, or Gaussians. While all these functions can be employed in the framework presented, we focus our discussion here on Gaussian kernels which also have an intuitive parametrization with regards to the scale of variability of  as quantified by the bandwidth parameters
Gaussian kernels in the context of free energy approximations have also been used in [32,37,18].
We define a joint probability distribution on the generalized coordinates q and the parameters z as follows:
where 1 D (z) is the indicator function on D and Z(θ) is the normalization constant, i.e.:
It is n
This content is AI-processed based on open access ArXiv data.
