Sparse Sampling of Water Density Fluctuations near Liquid-Vapor Coexistence

TL;DR

This paper introduces a generalized sparse sampling method for efficiently characterizing water density fluctuation free energetics near liquid-vapor coexistence, significantly reducing computational costs compared to traditional umbrella sampling.

Contribution

The authors extend and improve sparse sampling techniques with heuristics and adaptive strategies, enabling accurate free energy estimation in complex water systems.

Findings

Sparse sampling is nearly 100 times faster than umbrella sampling.

Method accurately characterizes cavitation and capillary evaporation.

Approach reduces computational cost for complex free energy landscapes.

Abstract

The free energetics of water density fluctuations in bulk water, at interfaces, and in hydrophobic confinement inform the hydration of hydrophobic solutes as well as their interactions and assembly. The characterization of such free energetics is typically performed using enhanced sampling techniques such as umbrella sampling. In umbrella sampling, order parameter distributions obtained from adjacent biased simulations must overlap in order to estimate free energy differences between biased ensembles. Many biased simulations are typically required to ensure such overlap, which exacts a steep computational cost. We recently introduced a sparse sampling method, which circumvents the overlap requirement by using thermodynamic integration to estimate free energy differences between biased ensembles. Here we build upon and generalize sparse sampling for characterizing the free energetics of water density fluctuations in systems near liquid-vapor coexistence. We also introduce sensible heuristics for choosing the biasing potential parameters and strategies for adaptively refining them, which facilitate the estimation of such free energetics accurately and efficiently. We illustrate the method by characterizing the free energetics of cavitation in a large volume in bulk water. We also use sparse sampling to characterize the free energetics of capillary evaporation for water confined between two hydrophobic plates. In both cases, sparse sampling is nearly two orders of magnitude faster than umbrella sampling. Given its efficiency, the sparse sampling method is particularly well suited for characterizing free energy landscapes for systems wherein umbrella sampling is prohibitively expensive.

Figures

Equations

$$F_v(\tilde{N})=F_v^{\bar{\lambda }}(\tilde{N})-U_{\bar{\lambda }}(\tilde{N})+F_{\bar{\lambda }},$$

$$F_v(\tilde{N})=F_v^{\varphi }(\tilde{N})-U_{\varphi }(\tilde{N})+F_{\varphi }.$$

$$F_{\varphi }={\int }_0^{\varphi }⟨\widetilde{N_v}{⟩}_{{\varphi }^{\prime }}d{\varphi }^{\prime }.$$

$$F_v(\tilde{N})=F_v^{\kappa ,N^{\ast }}(\tilde{N})-U_{\kappa ,N^{\ast }}(\tilde{N})+F_{\kappa ,N^{\ast }}$$

$$F_{N^{\ast }}={\int }_{N_{\mathrm{ref}}^{\ast }}^{N^{\ast }}⟨d{U}_{\kappa ,{N^{\ast }}^{\prime }}/d{N^{\ast }}^{\prime }{⟩}_{\kappa ,{N^{\ast }}^{\prime }}\mathrm{d}{N^{\ast }}^{\prime }={\int }_{N_{\mathrm{ref}}^{\ast }}^{N^{\ast }}\kappa ({N^{\ast }}^{\prime }-⟨\widetilde{N_v}{⟩}_{\kappa ,{N^{\ast }}^{\prime }})\mathrm{d}{N^{\ast }}^{\prime }$$

$$F_{\kappa ,N_{\mathrm{ref}}^{\ast }}={\int }_0^{\kappa }⟨d{U}_{{\kappa }^{\prime },N_{\mathrm{ref}}^{\ast }}/d{\kappa }^{\prime }{⟩}_{{\kappa }^{\prime },N_{\mathrm{ref}}^{\ast }}\mathrm{d}{\kappa }^{\prime }=\frac{1}{2}{\int }_0^{\kappa }⟨(\widetilde{N_v}-N_{\mathrm{ref}}^{\ast }{)}^2{⟩}_{{\kappa }^{\prime },N_{\mathrm{ref}}^{\ast }}\mathrm{d}{\kappa }^{\prime }.$$

$$\frac{d^2{F}_v^{\kappa ,N^{\ast }}}{d{\tilde{N}}^2}=\frac{d^2{F}_v}{d{\tilde{N}}^2}+\frac{d^2{U}_{\kappa ,N^{\ast }}}{d{\tilde{N}}^2}=\frac{d^2{F}_v}{d{\tilde{N}}^2}+\kappa >0.$$

$$\frac{d{F}_v}{d\tilde{N}}(\tilde{N}=⟨\widetilde{N_v}{⟩}_{\kappa ,N^{\ast }})\approx -\kappa (⟨\widetilde{N_v}{⟩}_{\kappa ,N^{\ast }}-N^{\ast })$$

$$P_v(\tilde{N})=⟨\delta ({\tilde{N}}_v-\tilde{N}){⟩}_0=\frac{1}{Q_0}\int d\overline{\text{R}}\delta ({\tilde{N}}_v-\tilde{N})e^{-\beta {\mathcal{H}}_0(\overline{\text{R}})}$$

$$P_v^{\bar{\lambda }}(\tilde{N})=⟨\delta ({\tilde{N}}_v-\tilde{N}){⟩}_{\bar{\lambda }}=\frac{1}{Q_{\bar{\lambda }}}\int d\overline{\text{R}}\delta ({\tilde{N}}_v-\tilde{N})e^{-\beta {\mathcal{H}}_{\bar{\lambda }}(\overline{\text{R}})},$$

$$P_v(\tilde{N})=\left(\frac{Q_{\bar{\lambda }}}{Q_0}\right)\frac{1}{Q_{\bar{\lambda }}}\int d\overline{\text{R}}\delta ({\tilde{N}}_v-\tilde{N})e^{-\beta {\mathcal{H}}_{\bar{\lambda }}}{e}^{\beta U_{\bar{\lambda }}({\tilde{N}}_v)}$$

$$P_v(\tilde{N})=\left(\frac{Q_{\bar{\lambda }}}{Q_0}\right)e^{\beta U_{\bar{\lambda }}(\tilde{N})}\frac{1}{Q_{\bar{\lambda }}}\int d\overline{\text{R}}\delta ({\tilde{N}}_v-\tilde{N})e^{-\beta {\mathcal{H}}_{\bar{\lambda }}}$$

$$P_v(\tilde{N})=e^{-\beta F_{\bar{\lambda }}}{e}^{\beta U_{\bar{\lambda }}(\tilde{N})}{P}_v^{\bar{\lambda }}(\tilde{N})$$

$$F_v(\tilde{N})=F_v^{\bar{\lambda }}(\tilde{N})-U_{\bar{\lambda }}(\tilde{N})+F_{\bar{\lambda }}$$

$$\frac{\partial F_v}{\partial \tilde{N}}=\frac{\partial F_v^{\kappa ,N^{\ast }}}{\partial \tilde{N}}-\frac{\partial U_{\kappa ,N^{\ast }}}{\partial \tilde{N}}=\frac{\partial F_v^{\kappa ,N^{\ast }}}{\partial \tilde{N}}-\kappa (\tilde{N}-N^{\ast }).$$

Full text

University of Pennsylvania] Department of Chemical and Biomolecular Engineering, University of Pennsylvania, Philadelphia, PA 19104, USA

Keywords: free energy method, umbrella sampling, thermodynamic integration, hydrophobic confinement

1 Introduction

At ambient conditions, water and many other liquids are close to coexistence with their vapor phase 1, 2, 3, 4, 5, 6, 7, 8, 9, 10. Liquid water in the vicinity of hydrophobic surfaces is destabilized further, situating interfacial waters at the edge of a dewetting transition, and rendering them sensitive to unfavorable perturbations 11, 12, 13, 14, 15, 16, 17, 18, 19. Moreover, when confined between hydrophobic surfaces, liquid water can become metastable (or even unstable) with respect to its vapor 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37. Such interfacial or confined waters, which are situated near liquid-vapor coexistence, play an important role in diverse processes, ranging from colloidal assembly and protein interactions, to heterogeneous vapor nucleation and the loss or recovery of superhydrophobicity 38, 39, 40, 41, 42, 43, 36, 23, 44, 45, 46. In particular, displacing the interfacial or confined waters can incur a substantial free energetic cost, and influence the kinetics of these process 47, 48, 35, 49. Thus, characterizing the free energetics of water density fluctuations can shed light on the mechanistic pathways associated with wetting-dewetting transitions in such systems, and also inform the magnitude of the associated barriers 50, 51, 52, 53, 33, 45, 54, 55, 36.

To characterize the free energetics of rare water density fluctuations, it becomes necessary to employ enhanced sampling techniques, such as umbrella sampling, which are computationally expensive. However, umbrella sampling can become prohibitively expensive when the simulations themselves are very expensive or when a large number of biased simulations are needed 56, 39, 57, 58, 59, 60, 61. To address these challenges and facilitate the computationally efficient characterization of the free energetics of water density fluctuations, here we build upon a sparse sampling method that we previously introduced 62, and generalize it to study systems near liquid-vapor coexistence. In the following sections, we will first use model (analytical) free energy landscapes to illustrate the sparse sampling method, the reasons underlying its successes, and situations when it is challenged. We then extend the sparse sampling method to tackle such challenging scenarios. Both the accuracy and the efficiency of sparse sampling rely on the choice of biasing potential parameters; we introduce sensible heuristics for choosing such parameters and strategies for adaptively refining them. We then highlight the efficiency of the sparse sampling method by using it to characterize the free energetics of water density fluctuations in a large volume in bulk water. Finally, we apply sparse sampling to a particularly challenging system, wherein water is confined between two hydrophobic surfaces, and the free energetics of water density fluctuations display two basins that are separated by a barrier.

2 Illustrating Sparse Sampling Using Model Landscapes

2.1 Umbrella Sampling vs Sparse Sampling

Characterizing the free energetics of order parameter fluctuations, which are too rare to be observed in unbiased molecular simulations, requires the use of non-Boltzmann sampling techniques such as umbrella sampling. To facilitate the sampling of such rare order parameter fluctuations, umbrella sampling prescribes the use of biasing potentials, which enhance the likelihood with which otherwise improbable fluctuations are sampled. Although umbrella sampling provides a powerful way to characterize the free energetics of rare order parameter fluctuations, it exacts a steep computational cost, which can become prohibitive, either when the biased simulations themselves are expensive 56, 39, 57, 58, 63, 61, 60, or when a large number of biased simulations are needed to span the order parameter range of interest. To address this challenge, we recently introduced a sparse sampling method 62, which employs sparsely separated biased simulations, and can be orders of magnitude more efficient than conventional umbrella sampling.

Although the sparse sampling method is generally applicable to any order parameter, here we focus on an order parameter that represents the smoothed (or coarse-grained) number of water molecules, $\tilde{N_{v}}$, in an observation volume, $v$, of interest. We bias $\tilde{N_{v}}$ in lieu of the closely related (discrete) number of waters, $N_{v}$, in $v$, because performing biased molecular dynamics (MD) simulations with a discrete order parameter results in impulsive forces. The precise definition of $\tilde{N_{v}}$ as well as its dependence on the atomic positions, $\overline{\textbf{R}}$, can be found in ref. 64. We wish to estimate the free energetics of water density (or rather $\tilde{N_{v}}$) fluctuations, $F_{v}(\tilde{N})=-k_{\rm{B}}T\ln P_{v}(\tilde{N})$, where $k_{\rm B}$ is the Boltzmann constant, $T$ is the system temperature, and $P_{v}(\tilde{N})=\langle\delta(\tilde{N_{v}}-\tilde{N})\rangle_{0}$, is the probability of observing $\tilde{N}$ coarse-grained waters in $v$; here, $\langle\mathcal{O}(\overline{\textbf{R}})\rangle_{0}\equiv\int d\overline{\textbf{R}}~{}\mathcal{O}(\overline{\textbf{R}})\exp(-\beta\mathcal{H}_{0})/Q_{0}$ represents the ensemble average of $\mathcal{O}(\overline{\textbf{R}})$, given the generalized Hamiltonian, $\mathcal{H}_{0}(\overline{\textbf{R}})$, and $Q_{0}\equiv\int d\overline{\textbf{R}}~{}\exp(-\beta\mathcal{H}_{0})$ is the corresponding partition function. To facilitate the sampling of $\tilde{N_{v}}(\overline{\textbf{R}})$ over the entire $\tilde{N}$-range of interest, we use potentials, $U_{\bar{\lambda}}(\tilde{N_{v}})$, which bias $\tilde{N_{v}}$, and are parametrized by the vector, $\bar{\lambda}$; that is, we perform biased simulations using the Hamiltonians, $\mathcal{H}_{\bar{\lambda}}(\overline{\textbf{R}})=\mathcal{H}_{0}(\overline{\textbf{R}})+U_{\bar{\lambda}}(\tilde{N_{v}}(\overline{\textbf{R}}))$. Using such biased simulations, we can estimate averages in the biased ensembles, $\langle\mathcal{O}(\overline{\textbf{R}})\rangle_{\bar{\lambda}}$, and in particular, we estimate the biased distributions, $P_{v}^{\bar{\lambda}}(\tilde{N})=\langle\delta(\tilde{N_{v}}-\tilde{N})\rangle_{\bar{\lambda}}$, and the corresponding free energetics, $F_{v}^{\bar{\lambda}}(\tilde{N})=-k_{\rm{B}}T\ln P_{v}^{\bar{\lambda}}(\tilde{N})$. For the $\tilde{N}$-range sampled by a biased simulation, $F_{v}(\tilde{N})$ can then be related to $F_{v}^{\bar{\lambda}}(\tilde{N})$ using the exact result

[TABLE]

where $F_{\bar{\lambda}}\equiv-k_{\rm{B}}T\ln\left(\frac{Q_{\bar{\lambda}}}{Q_{0}}\right)$ is the free energy difference between the biased and the unbiased ensembles, and $Q_{\bar{\lambda}}$ is the partition function corresponding to $\mathcal{H}_{\bar{\lambda}}$. The derivation of Equation 1 is included as supplementary information. Using $F_{v}^{\bar{\lambda}}(\tilde{N})$ obtained from the biased simulation, and the known functional form of $U_{\bar{\lambda}}(\tilde{N})$, $F_{v}(\tilde{N})$ can then be obtained to within the unknown constant offset, $F_{\bar{\lambda}}$.

In umbrella sampling, estimates of $F_{\bar{\lambda}}$ for the different biased ensembles are obtained by requiring an overlap in the range of $\tilde{N}$-values sampled by adjacent biased simulations, and in particular, by matching $F_{v}(\tilde{N})$-values obtained from different biased simulations in the overlap regions; algorithms such as the Weighted Histogram Analysis Method (WHAM) or Multi-state Bennet Acceptance Ratio (MBAR) are typically used for this purpose 65, 66, 67, 68, 69, 70. Such an umbrella sampling strategy is powerful and has facilitated the characterization of the free energetics of water density fluctuations in numerous contexts 16, 17, 18, 5, 19, 71. However, to satisfy the above overlap requirement, it becomes necessary to run many long biased simulations. In practice, the associated computational cost limits the size of $v$ as well as the level of simulation detail (e.g., classical vs ab initio MD 61, 60) for which $F_{v}(\tilde{N})$ can be characterized. In contrast, the sparse sampling method does not require overlap between adjacent biased distributions 62. Instead, sparse sampling employs thermodynamic integration to obtain $F_{\bar{\lambda}}$, and in principle, it can facilitate estimation of $F_{v}(\tilde{N})$ at sparsely distributed $\tilde{N}$-values orders of magnitude faster than conventional umbrella sampling. Below, we first describe the sparse sampling method, as introduced in ref. 62, and discuss its advantages and shortcomings; then, we generalize the sparse sampling method to address those shortcomings.

2.2 Sparse Sampling with a Linear Biasing Potential

In ref. 62, we introduced sparse sampling using a linear biasing potential, $U_{\phi}(\tilde{N_{v}})=\phi\tilde{N_{v}}$, parameterized by $\phi$; that is, $\bar{\lambda}=[\phi]$, and the corresponding biased ensemble Hamiltonian, $\mathcal{H}_{\phi}=\mathcal{H}_{0}+U_{\phi}(\tilde{N_{v}})$. Equation 1 then becomes:

[TABLE]

Instead of relying on an overlap in the $\tilde{N}$-values sampled by adjacent biased ensembles, sparse sampling prescribes estimating $F_{\phi}$ by integrating $dF_{\phi}/d\phi\equiv\langle dU_{\phi}/d\phi\rangle_{\phi}=\langle\tilde{N_{v}}\rangle_{\phi}$, that is, by using the identity:

[TABLE]

In practice, only a small number of biased simulations that sample well-separated $\tilde{N_{v}}$-values are performed; hence the term sparse sampling. Following Equation 3, the corresponding $F_{\phi}$-estimates are then obtained by numerically integrating the averages, $\langle\tilde{N_{v}}\rangle_{\phi}$, obtained from those simulations. To minimize integration errors in estimating $F_{\phi}$, the sparse sampling method thus relies on capturing the variation of $\langle\tilde{N_{v}}\rangle_{\phi}$, the average thermodynamic force (and the integrand in thermodynamic integration), with $\phi$. In other words, to accurately estimate $F_{\phi}$, the biased simulations must capture the functional dependence of $\langle\tilde{N_{v}}\rangle_{\phi}$ on $\phi$. Indeed, our choice of a linear biasing potential, $U_{\phi}(\tilde{N_{v}})=\phi\tilde{N_{v}}$, was informed by this requirement. In particular, because $d\langle\tilde{N_{v}}\rangle_{\phi}/d\phi=-\beta\langle\delta\tilde{N}_{v}^{2}\rangle_{\phi}\leq 0$, this choice ensures that the average thermodynamic force, $\langle\tilde{N_{v}}\rangle_{\phi}$, decreases monotonically with $\phi$, and for certain systems, the decrease can even be linear in $\phi$ 72. Here, $\langle\delta\tilde{N}_{v}^{2}\rangle_{\phi}$ is the variance of $P_{v}^{\phi}(\tilde{N})$. In ref.62, we used Equations 2 and 3 to efficiently estimate $F_{v}(\tilde{N})$ in small volumes in bulk water and at interfaces, and in large volumes in heterogeneous interfacial environments such as protein hydration shells. However, as highlighted in ref. 62, and discussed in further detail below, the use of a linear potential is not suitable for sparse sampling systems near liquid-vapor coexistence.

2.3 Sparse Sampling Systems Near and Far from Coexistence

To illustrate the differences between systems near and far from liquid-vapor coexistence, we consider two model systems characterized by distinct analytical free energy profiles (Figure 1a). System 1 is representative of a system far from liquid-vapor coexistence, and is monostable: $F_{v}^{\rm(1)}(\tilde{N})=\frac{1}{2}\kappa_{\rm m}(\tilde{N}-n^{(1)}_{\rm liq})^{2}$; the corresponding distribution, $P_{v}(\tilde{N})$, is Gaussian. In contrast, system 2 resembles a system close to liquid-vapor coexistence, and is bistable: $F_{v}^{\rm(2)}(\tilde{N})=\kappa_{\rm b}(\tilde{N}-n^{(2)}_{\rm vap})^{2}(\tilde{N}-n^{(2)}_{\rm liq})^{2}-\phi_{\rm b}\tilde{N}$, with distinct basins at low and high $\tilde{N}$-values. We have designed $F_{v}^{\rm(1)}(\tilde{N})$ and $F_{v}^{\rm(2)}(\tilde{N})$ so that the high $\tilde{N}$ (liquid) basins for the two model systems are similar to one another, and to a $3~{}$nm spherical observation volume in bulk water, which we will consider in the following section. To do so, we choose $\beta\kappa_{\rm m}=0.00318$, $n^{(1)}_{\rm liq}=3720$, $\beta\kappa_{\rm b}=1.2\times 10^{-10}$, $n^{(2)}_{\rm vap}=500$, $n^{(2)}_{\rm liq}=3600$ and $\beta\phi_{\rm b}=0.3$, where $\beta\equiv 1/k_{\rm{B}}T$. Although $F_{v}^{\rm(1)}(\tilde{N})$ and $F_{v}^{\rm(2)}(\tilde{N})$ are representative of commonly encountered free energy profiles far from and near liquid-vapor coexistence respectively, we note that these model systems are purely illustrative. For typical systems of interest, such $F_{v}(\tilde{N})$-profiles will not be known a priori; rather, our goal will be to characterize them using sparse sampling.

The first step in characterizing the $F_{v}(\tilde{N})$-profiles of the above systems using sparse sampling is to obtain the averages, $\langle\tilde{N_{v}}\rangle_{\phi}$, for several $\phi$-values. In lieu of performing biased simulations, we infer the biased free energetics, $F_{v}^{\phi}(\tilde{N})$, by reweighting the $F_{v}(\tilde{N})$-profiles for the two systems using Equation 2. The locations of the minima in $F_{v}^{\phi}(\tilde{N})$ then provide us with the averages, $\langle\tilde{N_{v}}\rangle_{\phi}$, that would be obtained from the corresponding biased simulations. We assume that all biased simulations are initialized with either high $\tilde{N_{v}}$-values ($\tilde{N_{v}}^{0}\approx n^{(1)}_{\rm liq}$) or low $\tilde{N_{v}}$-values ($\tilde{N_{v}}^{0}\approx 0$). Thus, if $F_{v}^{\phi}(\tilde{N})$ displays multiple minima separated by a barrier, we assign to $\langle\tilde{N_{v}}\rangle_{\phi}$ the minimum that is closest to $\tilde{N_{v}}^{0}$, noting that in realistic simulations, a barrier of only a few $k_{\rm{B}}T$ is sufficient to localize the system in the basin it was initialized in. In this instance, $\langle\tilde{N_{v}}\rangle_{\phi}$ does not represent a true biased ensemble average, but an estimate of the average that would be obtained from a biased simulation, and one that may suffer from serious inaccuracies due to the loss of ergodicity, which accompanies a bistable $F_{v}^{\phi}(\tilde{N})$.

In Figure 1b (symbols), we show estimates of $\langle\tilde{N_{v}}\rangle_{\phi}$, computed accordingly for equally spaced $\phi$-values. Although $\langle\tilde{N_{v}}\rangle_{\phi}$ decreases monotonically with $\phi$ for both systems, there are significant differences in the $\langle\tilde{N_{v}}\rangle_{\phi}$ vs. $\phi$ response of the two systems. For system 1, the dependence of $\langle\tilde{N_{v}}\rangle_{\phi}$ on $\phi$ is particularly simple; $\langle\tilde{N_{v}}\rangle_{\phi}$ decreases linearly with $\phi$. Moreover, identical estimates of $\langle\tilde{N_{v}}\rangle_{\phi}$ are obtained regardless of how the biased simulations are initialized. The integral, $F_{\phi}=\int_{0}^{\phi}\langle\tilde{N_{v}}\rangle_{\phi^{\prime}}d\phi^{\prime}$ (Equation 3), can thus be estimated accurately for each of the 8 equally-spaced $\phi$-values. Then, Equation 2 can be used to obtain estimates of $F_{v}(\tilde{N})$ at $\tilde{N}=\langle\tilde{N_{v}}\rangle_{\phi}$. These estimates are shown in Figure 1c, and are in excellent agreement with $F_{v}^{\rm(1)}(\tilde{N})$. The linear dependence of $\langle\tilde{N_{v}}\rangle_{\phi}$ on $\phi$ also means that the $\tilde{N_{v}}$-values sampled in the biased simulations are distributed across the $\tilde{N}$-range of interest, enabling determination of $F_{v}(\tilde{N})$ over that entire range. Thus, the sparse sampling method in conjunction with a linear potential is particularly well suited for characterizing the free energetics of systems far from coexistence.

In contrast with the gradual decrease of $\langle\tilde{N_{v}}\rangle_{\phi}$ with increasing $\phi$ for system 1, $\langle\tilde{N_{v}}\rangle_{\phi}$ decreases sharply over a narrow range of $\phi$-values for system 2, as shown in Figure 1b. Thus, for uniformly spaced $\phi$-values, high and low $\tilde{N_{v}}$-values are sampled well in the biased simulations, but intermediate values are not. Consequently, there is a large range of $\tilde{N}$-values ($200\lesssim\tilde{N}\lesssim 3400$) for which estimates of $F_{v}(\tilde{N})$ cannot be obtained (Figure 1c). Moreover, the gap in our knowledge of the functional form of $\langle\tilde{N_{v}}\rangle_{\phi}$ versus $\phi$ results in significant errors in our estimation of $F_{\phi}$ for the higher $\phi$-values, and thereby in estimates of $F_{v}(\tilde{N})$ for the lower $\tilde{N}$-values (Figure 1c). To make matters worse, for the $\phi$-range over which $\langle\tilde{N_{v}}\rangle_{\phi}$ decreases from high to low values, biased simulations initialized with low and high $\tilde{N_{v}}$-values result in very different estimates of $\langle\tilde{N_{v}}\rangle_{\phi}$; that is, hysteresis is observed in $\langle\tilde{N_{v}}\rangle_{\phi}$ versus $\phi$ (Figure 1b).

To understand this contrast between systems 1 and 2, we recognize that the $\tilde{N_{v}}$-values sampled in a biased simulation will be in the vicinity of a minimum in $F^{\phi}_{v}(\tilde{N})$. Extrema of $F^{\phi}_{v}(\tilde{N})$ obey $\partial F^{\phi}_{v}/\partial\tilde{N}=0$, and because $F^{\phi}_{v}(\tilde{N})=F_{v}(\tilde{N})+U_{\phi}(\tilde{N})-F_{\phi}$ (Equation 2), the extrema will satisfy $\partial F_{v}/\partial\tilde{N}=-\partial U_{\phi}/\partial\tilde{N}=-\phi$. In Figure 1d, we plot $\partial F_{v}/\partial\tilde{N}$ for both systems; for system 1, it is linear and increases monotonically with $\tilde{N}$, whereas for system 2, it varies non-monotonically with $\tilde{N}$. For system 1, there can thus be only one solution to $\partial F_{v}/\partial\tilde{N}=-\phi$, and thereby only one minimum in $F^{\phi}_{v}(\tilde{N})$ for any value of $\phi$ (Figure 1d). This is indeed the case as shown in Figure 1e. The corresponding biased distributions, $P^{\phi}_{v}(\tilde{N})$, will thus be unimodal, facilitating the efficient estimation of the thermodynamic force, $\langle\tilde{N_{v}}\rangle_{\phi}$, for all values of $\phi$. Moreover, these arguments apply to any system for which $\partial F_{v}/\partial\tilde{N}$ increases monotonically with $\tilde{N}$, and correspondingly, $F_{v}(\tilde{N})$ is convex over the entire range of $\tilde{N}$-values. In ref. 62, we showed that heterogeneous surfaces, such as the protein hydration shells, which display a wide range of chemistries from hydrophobic to hydrophilic, indeed satisfy these criteria.

In contrast, the non-monotonic variation of $\partial F_{v}/\partial\tilde{N}$ with $\tilde{N}$ for system 2 allows for the possibility of three solutions to $\partial F_{v}/\partial\tilde{N}=-\phi$ in the range $\phi_{\rm min}\leq\phi\leq\phi_{\rm max}$ (with $\beta\phi_{\rm min}=-0.4$ and $\beta\phi_{\rm max}=1$). In this $\phi$-range, $F^{\phi}_{v}(\tilde{N})$ exhibits two minima and a maxima, that is, two basins separated by a barrier, as shown in Figure 1f. The presence of a sufficiently large barrier in $F^{\phi}_{v}(\tilde{N})$ will localize the system to the basin it was initialized in, leading to hysteresis and erroneous estimates of $\langle\tilde{N_{v}}\rangle_{\phi}$. Moreover, for the wide range of $\tilde{N}$-values for which $\partial F_{v}/\partial\tilde{N}$ decreases, $\partial^{2}F^{\phi}_{v}/\partial\tilde{N}^{2}=\partial^{2}F_{v}/\partial\tilde{N}^{2}<0$. As a result, $F^{\phi}_{v}(\tilde{N})$ will have negative curvature in this $\tilde{N}$-range, precluding its sampling in the biased simulations, and resulting in a sharp decrease in $\langle\tilde{N_{v}}\rangle_{\phi}$ as $\phi$ is increased from $\phi_{\rm min}$ to $\phi_{\rm max}$. These issues arise from the presence of locally concave regions of $F_{v}(\tilde{N})$, and hinder the sparse sampling method from obtaining accurate estimates of $F_{v}(\tilde{N})$ over the entire $\tilde{N}$-range of interest. Because concave regions in $F_{v}(\tilde{N})$ are a common feature of systems near coexistence, the use of a linear potential is not appropriate for sparse sampling such systems.

2.4 Extending Sparse Sampling to Systems near Coexistence

To generalize the sparse sampling method for systems near coexistence, here we explore the use of a biasing potential with a different (non-linear) functional form. In particular, we consider the frequently used parabolic potential, $U_{\kappa,N^{*}}(\tilde{N_{v}})=\frac{1}{2}\kappa(\tilde{N_{v}}-N^{*})^{2}$, which is parameterized by $\bar{\lambda}=[\kappa,N^{*}]$. For sufficiently large $\kappa$-values, such a biasing potential can regularize the corresponding biased free energy profile, $F_{v}^{\kappa,N^{*}}(\tilde{N})$, ensuring that it is monostable. To illustrate this point, we first rewrite Equation 1 for this potential:

[TABLE]

For $F_{v}^{\kappa,N^{*}}(\tilde{N})$ to be monostable, $dF^{\kappa,N^{*}}_{v}/d\tilde{N}=0$ ought to have only one solution, and correspondingly, so should $dF_{v}/d\tilde{N}=-dU_{\kappa,N^{*}}/d\tilde{N}=-\kappa(\tilde{N}-N^{*})$. In Figure 2a, we again plot $dF^{(2)}_{v}/d\tilde{N}$ for system 2 (as in Figure 1d); however, we now search for its intersection with $-\kappa(\tilde{N}-N^{*})$ with $\beta\kappa=0.003$. In particular, we see that lines corresponding to $-\kappa(\tilde{N}-N^{*})$ for several different $N^{*}$-values intersect $dF^{(2)}_{v}/d\tilde{N}$ only once when a sufficiently large $\kappa$-value is chosen. As a result, the corresponding biased free energy profiles all display a single basin (Figure 2b). Moreover, the biased simulations facilitate sampling of $\tilde{N_{v}}$ over the entire range of $\tilde{N}$-values, and the corresponding $\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$-values increase systematically with $N^{*}$ (Figure 2c). In fact, the monotonic increase of $\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$ with $N^{*}$ is guaranteed because $\partial\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}/\partial N^{*}=\beta\kappa\langle\delta\tilde{N}_{v}^{2}\rangle_{\kappa,N^{*}}\geq 0$.

Using Equation 4, we can then estimate $F_{v}(\tilde{N})$ at the $\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$-values shown in Figure 2c. We can obtain $F_{v}^{\kappa,N^{*}}(\tilde{N})$ directly from the biased simulations, and to a good approximation, $\beta F_{v}^{\kappa,N^{*}}(\tilde{N}=\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}})\approx\frac{1}{2}\log(2\pi\langle\delta\tilde{N}_{v}^{2}\rangle_{\kappa,N^{*}})$. Similarly, $U_{\kappa,N^{*}}(\tilde{N})$ is readily obtained as: $U_{\kappa,N^{*}}(\tilde{N}=\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}})=\frac{1}{2}\kappa(\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}-N^{*})^{2}$. The first two terms on the right hand side of Equation 4, thus obtained, are shown in Figure 2d. In the sparse sampling approach, the final term, $F_{\kappa,N^{*}}$, is obtained using thermodynamic integration. The estimation of $F_{\kappa,N^{*}}$ is simplified somewhat if a single $\kappa$-value is adopted, and $N^{*}$ is varied in the biased simulations. In particular, by using $dF_{\kappa,N^{*}}/dN^{*}=\langle dU_{\kappa,N^{*}}/dN^{*}\rangle_{\kappa,N^{*}}$, and performing thermodynamic integration from a reference value, $N^{*}_{\rm ref}$, to the $N^{*}$ of interest, $F_{N^{*}}\equiv F_{\kappa,N^{*}}-F_{\kappa,N^{*}_{\rm ref}}$ can be obtained as:

[TABLE]

In this way, $F_{\kappa,N^{*}}=F_{N^{*}}+F_{\kappa,N^{*}_{\rm ref}}$ is obtained to within an unknown constant, $F_{\kappa,N^{*}_{\rm ref}}$, by using estimates of $\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$ obtained from the biased simulations. The constant, $F_{\kappa,N^{*}_{\rm ref}}$, if desired, can be obtained by thermodynamic integration from 0 to $\kappa$ as:

[TABLE]

Alternatively, if $N^{*}_{\rm ref}$ is chosen to be a basin of $F_{v}(\tilde{N})$, there will be significant overlap between the unbiased and the (reference) biased ensembles, allowing $F_{\kappa,N^{*}_{\rm ref}}$ to be estimated using free energy perturbation 73.

For system 2, we first plot the thermodynamic force (the integrand in Equation 5) as a function of $N^{*}$ in Figure 2e. We note that the integrand varies non-monotonically with $N^{*}$ in contrast with the corresponding integrand for a linear biasing potential (Figure 1b). We then choose $N^{*}_{\rm ref}$ to be $n^{(1)}_{\rm liq}$ (which is a minimum of $F^{(2)}_{v}(\tilde{N}))$, and use numerical integration to obtain $F_{N^{*}}$ for all the simulated $N^{*}$-values (Figure 2f). Having obtained $F_{N^{*}}$, we can then use Equation 4 to obtain sparse sampling estimates of $F_{v}(\tilde{N})$ at $\tilde{N}=\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$ to within the constant, $F_{\kappa,N^{*}_{\rm ref}}$. We do not estimate $F_{\kappa,N^{*}_{\rm ref}}$, but instead shift $F_{v}(\tilde{N})$ vertically to set the zero of $F_{v}(\tilde{N})$ at $\tilde{N}=N^{*}_{\rm ref}$. As shown in Figure 2g, the free energy profile, $F_{v}(\tilde{N})$, thus obtained, is in excellent agreement with $F_{v}^{\rm(2)}(\tilde{N})$, highlighting that when the sparse sampling method is used in conjunction with a parabolic potential, it can be used to characterize the free energetics of systems near coexistence. A discussion of the similarities and differences between the sparse sampling method and related techniques, such as Umbrella Integration (UI) 74, is included as supplementary information.

2.5 Choosing the Biasing Potential Parameters Judiciously

The sparse sampling method can be used to efficiently sample $F_{v}(\tilde{N})$; however, its practical implementation relies on the judicious choice of the biasing potential parameters, $\kappa$ and $N^{*}$. As shown in Figure 2e, the average thermodynamic force, $dF_{N^{*}}/dN^{*}$, does not vary monotonically with $N^{*}$; thus, care must be exercised in choosing $N^{*}$-values that capture its functional form. To do so, we propose choosing $N^{*}$-values in an adaptive fashion. To begin with, we employ equally spaced $N^{*}$-values, which span the desired $\tilde{N}$-range (from 0 to $n_{\rm{liq}}$). The average thermodynamic force, $dF_{N^{*}}/dN^{*}=\langle dU_{\kappa,N^{*}}/dN^{*}\rangle_{\kappa,N^{*}}=\kappa(N^{*}-\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}})$, obtained from the biased simulations, is then inspected, and in regions where it changes relatively abruptly, $N^{*}$-values are chosen to perform additional simulations. In section 3, we will illustrate this strategy to estimate $F_{v}(\tilde{N})$ in a large spherical volume in bulk water.

Care must also be exercised in our choice of $\kappa$. In Figures 2a and 2b, we showed that when a parabolic potential with a sufficiently large $\kappa$ is used, $dF_{v}/d\tilde{N}=-\kappa(\tilde{N}-N^{*})$ has only one solution for all choices of $N^{*}$, which then implies that $dF_{v}^{\kappa,N^{*}}/d\tilde{N}=0$ has only one solution, and the biased free energetics, $F_{v}^{\kappa,N^{*}}(\tilde{N})$, are monostable for all $N^{*}$. However, if $\kappa$ is chosen to be too small, $dF_{v}/d\tilde{N}=-\kappa(\tilde{N}-N^{*})$ can have 3 roots for certain values of $N^{*}$, with the corresponding $F_{v}^{\kappa,N^{*}}(\tilde{N})$ being bistable, as shown in Figures 3a and 3b (green) for $\beta\kappa=0.0005$. To ensure that $F_{v}^{\kappa,N^{*}}(\tilde{N})$ is monostable for all $N^{*}$-values, $\kappa$, the slope of the lines in Figure 3a, must be larger in magnitude than the largest negative slope of $dF_{v}/d\tilde{N}$. In other words, a good choice of $\kappa$ is one that ensures the convexity of $F_{v}^{\kappa,N^{*}}(\tilde{N})$ for all $\tilde{N}$ by satisfying the criterion:

[TABLE]

Thus, if $F^{\prime\prime}_{\rm max}$ is the largest value that $-d^{2}F_{v}^{\kappa,N^{*}}/d\tilde{N}^{2}$ takes, we ought to choose $\kappa>F^{\prime\prime}_{\rm max}$. However, given that we wish to estimate $F_{v}(\tilde{N})$, the values taken by $d^{2}F_{v}/d\tilde{N}^{2}$ are not known to us a priori, and neither is $F^{\prime\prime}_{\rm max}$. Moreover, choosing excessively large values of $\kappa$ is also not recommended, because they can lead to significant errors in estimates of the average thermodynamic force, $\kappa(N^{*}-\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}})$, by amplifying small errors in our estimates of $\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$ 68. Such errors in the thermodynamic force are then propagated to $F_{N^{*}}$ through thermodynamic integration (Equation 5), and eventually onto $F_{v}(\tilde{N})$ (Equation 4). Thus, although it is important to choose a sufficiently large $\kappa$-value that exceeds $F^{\prime\prime}_{\rm max}$, it is not advisable to choose an excessively large $\kappa$-value because the length of time for which biased simulations must be run to obtain comparable errors in $F_{v}(\tilde{N})$ grows with $\kappa^{2}$.

Given these competing requirements, how do we optimally choose $\kappa$? As a rule of thumb, we propose choosing $\beta\kappa$ to be a multiple of $\langle\delta\tilde{N}_{v}^{2}\rangle_{0}^{-1}$, that is, $\beta\kappa=\alpha\langle\delta\tilde{N}_{v}^{2}\rangle_{0}^{-1}$ with $\alpha$ in the range, $3\leq\alpha\leq 5$. The choice is underpinned by the observation that $\beta F^{\prime\prime}_{\rm max}$, the maximum negative value assumed by $\beta d^{2}F_{v}/d\tilde{N}^{2}$, is often comparable to its (positive) value in the liquid basin, $\langle\delta\tilde{N}_{v}^{2}\rangle_{0}^{-1}$, which is readily accessible from an unbiased simulation. We have found that $\kappa$ chosen in accordance with this rule of thumb tends to satisfy the criterion in Equation 7 (that is, $\kappa>F^{\prime\prime}_{\rm max}$). However, in the unlikely event that our initial choice of $\kappa$ is too small, our biased simulations will exhibit a number of symptoms. In particular, if $\kappa<F^{\prime\prime}_{\rm max}$, for certain $N^{*}$-values, $d^{2}F_{v}^{\kappa,N^{*}}/d\tilde{N}^{2}$ will be negative over a range of $\tilde{N}$-values, which will not be sampled in the biased simulations. Consequently, there will be a sharp change in $\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$ as a function of $N^{*}$, as shown in Figure 3c. A second symptom of a small $\kappa$ is the appearance of hysteresis in $\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$ (Figure 3c), making it prudent to perform two sets biased simulations that are initialized in the liquid and vapor basins, respectively. Thus, when $\kappa$ is chosen to be smaller than $F^{\prime\prime}_{\rm max}$, we encounter all of the same challenges that were encountered when using a linear biasing potential.

If our initial choice of $\kappa$ turns out to be too small, how do we choose a revised $\kappa$-value? Can our initial set of biased simulations inform this revised choice? To address these questions, we recognize that the response curves can also provide estimates of $dF_{v}/d\tilde{N}$ at $\tilde{N}=\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$ through:

[TABLE]

To obtain this equation, we recognize that $dF_{v}/d\tilde{N}=dF_{v}^{\kappa,N^{*}}/d\tilde{N}-\kappa(\tilde{N}-N^{*})$ (Equation 4), and that $dF_{v}^{\kappa,N^{*}}/d\tilde{N}$ ought to be 0 for $\tilde{N}=\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$. Estimates of $dF_{v}/d\tilde{N}$ thus obtained are shown in Figure 3d, and display two clusters at low and high $\tilde{N}$-values. Although $dF_{v}/d\tilde{N}$ increases with increasing $\tilde{N}$ within the clusters, it decreases as we move from the low-$\tilde{N}$ cluster to the high-$\tilde{N}$ cluster, implying that $d^{2}F_{v}/d\tilde{N}^{2}<0$ for intermediate $\tilde{N}$-values. By identifying points that display the sharpest decrease in $dF_{v}/d\tilde{N}$, and estimating the magnitude of the slope of the line connecting those points, $F^{\prime\prime}_{\rm est}$ (Figure 3d), we can obtain a rough estimate of $F^{\prime\prime}_{\rm max}$; in fact, $F^{\prime\prime}_{\rm est}$ provides a lower bound on $F^{\prime\prime}_{\rm max}$, that is, $F^{\prime\prime}_{\rm est}\lesssim F^{\prime\prime}_{\rm max}$. Thus, if our choice of $\kappa$ is smaller than $F^{\prime\prime}_{\rm est}$ — that is, if $\kappa<F^{\prime\prime}_{\rm est}\lesssim F^{\prime\prime}_{\rm max}$ — it was clearly too small. To be safe, we propose the heuristic that $\kappa$ be larger than some multiple of $F^{\prime\prime}_{\rm est}$; that is, $\kappa>\alpha F^{\prime\prime}_{\rm est}$, with $3\leq\alpha\leq 5$. If this criterion is violated and either of the two low-$\kappa$ symptoms, discussed above, are observed, we recommend that the revised $\kappa$ be chosen according to: $\kappa_{\rm new}=\alpha~{}\rm{max}(\kappa_{\rm old},F^{\prime\prime}_{\rm est})$; that is, the revised (larger) value of $\kappa$ ought to be chosen to be $\alpha$ times either the previous (smaller) value of $\kappa$ or $F^{\prime\prime}_{\rm est}$, whichever is greater. Thus, inspecting $dF_{v}/d\tilde{N}$ not only provides us with a third symptom of a small $\kappa$, but it also provides a way to choose the revised (higher) value of $\kappa$. Simulations with this higher value of $\kappa$ can then be repeated at all the $N^{*}$-values. However, using a single $\kappa$-value for all biased simulations is not necessary. We can also combine our initial low-$\kappa$ simulations (that are outside the hysteresis region) with the newer high-$\kappa$ simulations to obtain $F_{v}(\tilde{N})$. In section 4, we will illustrate how to integrate simulations with different $\kappa$-values.

3 Fluctuations in a Large Volume in Bulk Water

Here we use sparse sampling to characterize the free energetics, $F_{v}(\tilde{N})$, of observing $\tilde{N}$ waters in a spherical volume of radius, $R_{v}=3$ nm, in bulk water. Because the volume contains a large number, $n_{\rm liq}\approx 3700$, water molecules on average (Figure 4a), a large number of simulations are needed to obtain $F_{v}(\tilde{N})$ using umbrella sampling. Indeed, 147 biased simulations, each run for 1 ns, were needed to ensure overlap between adjacent biased ensembles, and obtain the free energy profile, $F_{v}(\tilde{N})$, shown in Figure 4d (red curve). Sparse sampling provides an efficient alternative to obtain $F_{v}(\tilde{N})$ at sparsely sampled $\tilde{N}$-values; given the proximity of bulk water at ambient conditions to liquid-vapor coexistence, we make use of a parabolic potential, $U_{\kappa,N^{*}}$. We first perform an unbiased simulation to estimate, $\langle\delta\tilde{N}_{v}^{2}\rangle_{0}$, and choose $\beta\kappa=\alpha\langle\delta\tilde{N}_{v}^{2}\rangle_{0}^{-1}=0.016$ using $\alpha=5$. We then perform biased simulations with this value of $\kappa$, and 7 sparsely spaced $N^{*}$-values that span from 0 to $n^{(1)}_{\rm liq}$.

The average number of water molecules, $\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$, obtained from the 7 biased simulations increase as $N^{*}$ is increased (Figure 4b, green squares). However, the thermodynamic force, $dF_{N^{*}}/dN^{*}=\kappa({N^{*}}-\langle\tilde{N_{v}}\rangle_{\kappa,{N^{*}}})$, which must be integrated to obtain $F_{N^{*}}$ (Equation 5), varies non-monotonically with $N^{*}$, first increasing at low $N^{*}$, then decreasing gradually, only to increase abruptly at the highest $N^{*}$-values (Figure 4c). This abrupt increase suggests that the functional form of the integrand may not be adequately captured in the high $N^{*}$-region. Nevertheless, the sparse sampled $F_{v}(\tilde{N})$-values (Figure 4d, green squares), obtained using the $\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$-values from these 7 simulations alone, not only capture $F_{v}(\tilde{N})$ qualitatively, but also display semi-quantitative agreement with the exact $F_{v}(\tilde{N})$ obtained from umbrella sampling (red curve), with the error in the free energies being roughly 15%.

To better capture the functional form of $dF_{N^{*}}/dN^{*}$, and to improve the accuracy with which the $F_{N^{*}}$-values, and thereby the sparse sampled $F_{v}(\tilde{N})$ profile are estimated, we include four additional simulations; two each between the three highest $N^{*}$-values. This increases the resolution with which we are able to characterize $dF_{N^{*}}/dN^{*}$ in a targeted fashion, providing estimates in the $N^{*}$-range where they are needed the most. The corresponding results, shown by the yellow circles in Figures 4b-d, highlight the improvement in our knowledge of the functional form of the $dF_{N^{*}}/dN^{*}$ (Figure 4c), and the corresponding improvement in our estimation of $F_{v}(\tilde{N})$ (Figure 4d), which result from the four additional simulations. This procedure of inspecting the variation of the integrand with $N^{*}$, and employing additional simulations in regions of substantive variation to adaptively augment our characterization of the functional form of the integrand can be repeated further to obtain even more accurate estimates of $F_{v}(\tilde{N})$. Here, we perform two more simulations at such $N^{*}$-values, and obtain estimates of $F_{v}(\tilde{N})$ (blue triangles), which are not only in excellent quantitative agreement with the umbrella sampling results, but also differ only marginally from estimates in the previous iteration (yellow circles), suggesting that convergence has been achieved. However, in contrast with the umbrella sampling, which incurred a computational overhead of 147 ns, obtaining the sparse sampled $F_{v}(\tilde{N})$ (blue triangles) in Figure 4d required only 13 biased simulations, each run for 0.2 ns, for a total simulation time of 2.6 ns. Thus, we were able to estimate $F_{v}(\tilde{N})$ roughly 2 orders of magnitude faster using sparse sampling.

Although no sharp jumps in $\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$ were observed upon increasing $N^{*}$ in Figure 4b, we nevertheless estimate $\partial F_{v}/\partial\tilde{N}$ at those $\tilde{N}=\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$-values (Figure 4e, symbols) to ensure that our choice of $\kappa$ was sufficiently large. By connecting the points that display the sharpest decrease in $\partial F_{v}/\partial\tilde{N}$ with a line (Figure 4e, black line), we estimate $\beta F^{\prime\prime}_{\rm est}\approx 0.004$, and correspondingly, $\kappa/F^{\prime\prime}_{\rm est}\approx 4.2$, suggesting that our initial choice of $\kappa$ was sufficiently large.

4 Water in hydrophobic confinement

Here we use sparse sampling to characterize the free energetics of water density fluctuations in confinement between two square hydrophobic plates, roughly 4 nm by 4 nm in size, that are separated by 1.6 nm (Figure 5a). For this separation, $F_{v}(\tilde{N})$ is known to feature two basins 35: a liquid (high $\tilde{N}$) basin and a vapor (low $\tilde{N}$) basin (Figure 5c, red curve). Moreover, in the vicinity of the barrier that separates the two basins, that is, near the maximum in $F_{v}(\tilde{N})$, there appears to be a “kink” or a sharp change in the slope of $F_{v}(\tilde{N})$. In other words, over a small range of $\tilde{N}$-values, $F_{v}(\tilde{N})$ displays high negative curvature. Recent work has attributed such a kink in $F_{v}(\tilde{N})$ to the presence of different dewetted morphologies on either side of the kink 35, 36. In particular, Remsing et al. showed that the low-$\tilde{N}$ (vapor) side of the kink features vapor tubes spanning the confined region between the two plates, whereas the high-$\tilde{N}$ (liquid) side of the kink features isolated cavities that appear adjacent to one plate or the other 35. Similar behavior has also been reported for dewetting on nanotextured hydrophobic surfaces, wherein multiple kinks were observed in $F_{v}(\tilde{N})$, and were associated with transitions between different dewetted morphologies. 36

To sparse sample the system shown in Figure 5a, we first choose $\beta\kappa=0.016$, following our rule of thumb (section 2.5). Given that this choice makes use of the curvature of $F_{v}(\tilde{N})$ in the liquid basin, it is unlikely to yield a $\kappa$-value that exceeds the large negative curvature in $F_{v}(\tilde{N})$ in the vicinity of a kink. We thus recognize that our initial choice of $\kappa$ is likely to result in bistable biased free energy profiles, $F_{v}^{\kappa,N^{*}}(\tilde{N})$; however, because the large negative curvature in $F_{v}(\tilde{N})$ is localized to a small $\tilde{N}$-region near the kink, $F_{v}^{\kappa,N^{*}}(\tilde{N})$ ought to be bistable for only a small range of $N^{*}$-values. Using our initial choice of $\kappa$, we run biased simulations for 12 sparsely distributed $N^{*}$-values; the resulting values of $\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$ are shown in Figure 5b as a function of $N^{*}$. The increase in $\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$ with $N^{*}$ appears to be gradual, apart from a relatively sharp increase near $N^{*}=500$. By integrating $dF_{N^{*}}/dN^{*}=\kappa(N^{*}-\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}})$ with respect to $N^{*}$ to obtain estimates of $F_{N^{*}}$, and using Equation 4, we then obtain estimates of $F_{v}(\tilde{N})$ at $\tilde{N}=\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$ (Figure 5c, blue circles). The results display qualitative, but not quantitative agreement with the exact umbrella sampling results (Figure 5c, red line).

To obtain more accurate estimates of $F_{v}(\tilde{N})$, and in particular, to ascertain whether the sharp increase in $\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$ near $N^{*}=500$ is symptomatic of bistability in $F_{v}^{\kappa,N^{*}}(\tilde{N})$, we plot estimates of $dF_{v}/d\tilde{N}$ at $\tilde{N}=\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$ in Figure 5d. We find that $dF_{v}/d\tilde{N}$ decreases sharply in the vicinity of $\tilde{N}=500$, and that the magnitude of the corresponding slope, $\beta F^{\prime\prime}_{\rm est}\approx 0.0086$, and correspondingly, $\kappa/F^{\prime\prime}_{\rm est}\approx 1.86$, which is less than our chosen threshold of $\alpha=3$. Taken together with the sharp increase in $\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$, this suggests that our initial choice of $\kappa$ was not large enough to prevent bistability in $F_{v}^{\kappa,N^{*}}(\tilde{N})$ near $N^{*}=500$. Thus, either the biased simulations must be extended until a converged $F_{v}^{\kappa,N^{*}}(\tilde{N})$ is obtained with sufficient sampling of both its basins, or a larger value of $\kappa$ ought to be chosen to ensure that the biased free energy profiles are monostable. This is a non-trivial choice, which must take into consideration the possibility of hidden barriers, that is, barriers in collective variables that are orthogonal to $\tilde{N_{v}}$, and the anticipated height of those barriers. Here we follow the prescription outlined in section 2.5, and choose a revised (higher) value of $\beta\kappa=0.048$. However, given that the sharp decrease in $dF_{v}/d\tilde{N}$ occurs over only a small range of $\tilde{N}$-values (Figure 5d), we do not repeat all the biased simulations using this higher value of $\kappa$. Instead, we repeat our biased simulations with the higher $\kappa$-value only for $N^{*}$-values around 500; the $\kappa$ and $N^{*}$-values that we adopt for the old and new biased simulations are shown in Figure 5e. The higher $\kappa$ simulations allow us to augment our knowledge of $dF_{v}/d\tilde{N}$, and obtain a revised estimate for $F^{\prime\prime}_{\rm est}$ (Figure 5f). Comparing the revised values, $\beta F^{\prime\prime}_{\rm est}\approx 0.0096$, and $\beta\kappa=0.048$, we get a revised estimate of $\kappa/F^{\prime\prime}_{\rm est}\approx 5$, which is above our threshold of $\alpha=3$. Thus, we expect our revised choice of $\kappa$ to be large enough to facilitate sampling of $\tilde{N_{v}}$ in the vicinity of the kink, without the appearance of bistability in $F_{v}^{\kappa,N^{*}}(\tilde{N})$ or hysteresis in $\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$.

Given that we no longer employ a single value of $\kappa$, in addition to estimating how $F_{\kappa,N^{*}}$ varies with $N^{*}$ (for a given $\kappa$), we must also estimate the change in $F_{\kappa,N^{*}}$ as $\kappa$ is changed (for the two $N^{*}$ values shown in Figure 5e). We find that the latter can be readily estimated using the Bennett Acceptance Ratio (BAR) method 75 because changing $\kappa$ (while keeping $N^{*}$ the same) leads to good overlap in $\tilde{N_{v}}$-values sampled by the biased ensembles with different $\kappa$-values. In Figure 5g, we plot estimates of $F_{\kappa,N^{*}}$ for all the biased ensembles considered. With all the $F_{\kappa,N^{*}}$-values in hand, Equation 4 can then be used to obtain $F_{v}(\tilde{N})$ at $\tilde{N}=\langle\tilde{N_{v}}\rangle_{\kappa,N^{*}}$; the resulting sparse sampled $F_{v}(\tilde{N})$-estimates are shown in Figure 5h, and display good agreement with the umbrella sampling result. Finally, we note that the total computational cost for performing sparse sampling is roughly 3 ns (16 simulations, 200 ps each), whereas the cost for performing umbrella sampling is roughly 230 ns (38 simulations, 6 ns each), once again resulting in a nearly two orders of magnitude improvement in computational efficiency.

5 Conclusions

A characterization of the free energetics of water density fluctuations in bulk water, at interfaces, and in hydrophobic confinement has substantially advanced our understanding of hydrophobic hydration, interactions, and assembly 76, 77, 78, 4, 14, 79, 16, 17, 18, 80, 71, 81. Such a characterization typically requires enhanced sampling techniques, such as umbrella sampling, which can be expensive 16, 64. In particular, the requirement that order parameter distributions obtained from adjacent biased simulations display overlap, exacts a steep computational cost, and in certain contexts, can render umbrella sampling prohibitively expensive. We recently introduced a sparse sampling method, which circumvents the overlap requirement central to umbrella sampling, and instead uses thermodynamic integration to estimate free energy differences between the biased and unbiased ensembles 62. By employing only a few sparsely separated biased simulations, the sparse sampling method was able to estimate free energy profiles at a small fraction of the computational cost exacted by umbrella sampling.

Although the sparse sampling method, as introduced in ref. 62, is efficient in its estimation of $F_{v}(\tilde{N})$, it is only suitable for systems with convex free energy landscapes, thereby excluding the important class of systems which are at or close to coexistence. To estimate $F_{v}(\tilde{N})$ for systems near liquid-vapor coexistence, we generalized the sparse sampling method to sample free energy landscapes with concave regions. Using model systems with analytical landscapes, we first highlighted the challenges associated with sparse sampling systems near coexistence. We then illustrated the use of a harmonic potential, $U_{\kappa,N^{*}}=\frac{1}{2}\kappa(\tilde{N_{v}}-N^{*})^{2}$, to regularize the free energy landscape. Importantly, both the accuracy and efficiency of sparse sampling rely critically on the choice of the force constant, $\kappa$, as well as the $N^{*}$-values for which biased simulations are performed. In order to minimize computational expense while ensuring reasonable accuracy, we proposed the use of sensible heuristics which guide the initial choice of the biasing potential parameters, $\kappa$ and $N^{*}$. We also introduced criteria for determining whether these initial choices needed to be refined, and proposed strategies for adaptively choosing new values of $\kappa$ and $N^{*}$, when needed.

To illustrate the efficiency of the sparse sampling method, we studied two realistic systems close to liquid-vapor coexistence. First, we characterized $F_{v}(\tilde{N})$ in a large volume containing roughly 3700 waters in bulk water at ambient conditions Although liquid water is stable relative to its vapor at ambient conditions, it can undergo a cavitation transition undergo sufficient tension, which gives rise to concave features in $F_{v}(\tilde{N})$ 82, 83, 55, 84. By choosing $N^{*}$-values adaptively, we obtained an accurate characterization of the free energetics of this process. We also computed $F_{v}(\tilde{N})$ for water confined between two hydrophobic plates — a system with distinct liquid and vapor basins. Due to the presence of a small region of high negative curvature in $F_{v}(\tilde{N})$, this system provides a more stringent test to the sparse sampling method. By choosing $\kappa$-values adaptively, sparse sampling again facilitates an accurate characterization $F_{v}(\tilde{N})$. In both cases, the computational expense of estimating $F_{v}(\tilde{N})$ is roughly two orders of magnitude smaller than umbrella sampling.

Although we focus on the free energetics of water density fluctuations, the sparse sampling method itself is general, and can be readily used with any order parameter. Moreover, it can also be extended to higher dimensions, and used to sample multiple order parameters. The efficiency of the sparse sampling method makes it particularly attractive for characterizing the free energy landscapes of systems for which umbrella sampling is prohibitively expensive. Such instances may arise when the (biased) simulations themselves are inherently expensive; for example, it can be quite expensive to simulate crystalline systems with slow dynamics, or to explicitly account for polarizability or electronic effects 59, 63, 85, 86. Alternatively, when a large number of biased simulations are needed to satisfy the overlap requirement, and span the order parameter range of interest, umbrella sampling can again become very expensive; examples include the order parameter range of interest itself being large, inherent system fluctuations being small, and a high dimensional order parameter phase space 87, 88, 89. In each of these contexts, wherein performing extensive umbrella sampling simulations may not be feasible, we hope that the estimation of free energy profiles will be facilitated by the sparse sampling method presented here.

{acknowledgement}

The authors gratefully acknowledge financial support from the University of Pennsylvania Materials Research Science and Engineering Center (NSF UPENN MRSEC DMR 1720530), National Science Foundation grants CBET 1511437 and CBET 1652646, as well as the Charles E. Kaufman Foundation grant KA-2015-79204. A.J.P. is thankful to Talid Sinno for numerous helpful discussions.

Simulation Details We perform all-atom molecular dynamics (MD) simulations using the GROMACS package 90, suitably modified to incorporate the biasing potentials of interest. The leap frog algorithm 91 was used to integrate the equations of motion with a 2 fs time-step, and periodic boundary conditions were employed in all dimensions. The SPC/E model of water 92 was employed in all simulations, with the hydrogen bonds of water molecules being constrained using the SHAKE algorithm 93. The Particle Mesh Ewald algorithm 94 was used to compute the long range electrostatic interactions, and short range interactions Lennard-Jones and electrostatic interactions were truncated at 1 nm. The canonical velocity-rescaling thermostat 95 was used to maintain a system temperature of 300 K.

Spherical Volume in Bulk Water

The spherical observation volume with a radius, $R_{v}=3$ nm, was placed at the center of a cubic simulation box with a side of 10 nm. Simulations were performed in the NPT-ensemble at a pressure of 1 bar, which was maintained using the Parrinello-Rahman barostat 96. To allow for equilibration, the first 100 ps were discarded from each biased simulation.

Water in Hydrophobic Confinement

The simulation setup for the hydrophobic plates is the same as that used in ref. 35. The plates are 4 nm by 4 nm in size and are separated by 1.6 nm, so that the confinement region is spanned by a $4\times 4\times 1.6$ nm3 cuboid observation volume that is positioned between the plates. Simulations were performed in NVT-ensemble with a buffering water-vapor interface, which keeps the system at its coexistence pressure 51, 16, 64. To allow for equilibration, the first 500 ps were discarded from each biased simulation.

6 Relating Free Energetics in the Biased and Unbiased Ensembles

Here we include a derivation of Equation 1 in the main text, which often serves as the starting point for enhanced sampling methods such as umbrella sampling. The probability, $P_{v}(\tilde{N})$, of observing $\tilde{N}$ coarse-grained waters in the unbiased ensemble with the generalized Hamiltonian, $\mathcal{H}_{0}$, is given by:

[TABLE]

where $Q_{0}\equiv\int d\overline{\textbf{R}}~{}e^{-\beta\mathcal{H}_{0}(\overline{\textbf{R}})}$ is the partition function of the unbiased ensemble, and $\overline{\textbf{R}}$ represents the positions of all atoms in the system. Upon addition of a biasing potential, the statistics of the biased ensemble are dictated by the Hamiltonian, $\mathcal{H}_{\bar{\lambda}}(\overline{\textbf{R}})=\mathcal{H}_{0}(\overline{\textbf{R}})+U_{\bar{\lambda}}(\tilde{N}_{v}(\overline{\textbf{R}}))$. In particular, the biased distribution,

[TABLE]

where $Q_{\bar{\lambda}}\equiv\int d\overline{\textbf{R}}~{}e^{-\beta\mathcal{H}_{\bar{\lambda}}(\overline{\textbf{R}})}$ is the partition function of the biased ensemble. To relate $P_{v}(\tilde{N})$ and $P_{v}^{\bar{\lambda}}(\tilde{N})$, we can rewrite Equation 9 as:

[TABLE]

Recognizing that the delta function allows us to rewrite $e^{\beta U_{\bar{\lambda}}({\tilde{N}}_{v})}$ as $e^{\beta U_{\bar{\lambda}}({\tilde{N}})}$, and pull it outside the integral (because it no longer depends on $\overline{\textbf{R}}$), we get

[TABLE]

Using Equation 10 and defining the free energy difference between the biased and the unbiased ensembles, $F_{\bar{\lambda}}\equiv(-1/\beta)\ln(Q_{\bar{\lambda}}/Q_{0})$, we then get:

[TABLE]

Finally, taking the logarithm of both sides and multiplying by $-(1/\beta)$, we get the desired result:

[TABLE]

7 Comparing Sparse Sampling to Umbrella Integration

Here, we compare sparse sampling with the umbrella integration (UI) method of Kastner and Thiel 74, Kastner:JCP:2006. Although both methods employ biased simulations in conjunction with thermodynamic integration to obtain free energy profiles, there are important differences both in the spirit of the two methods, and in their respective implementations. In lieu of estimating the offset, $F_{\kappa,N^{*}}$, UI seeks to eliminate it from the equation by considering the derivative of Equation 4 of the main text. Thus, the central equation for UI becomes:

[TABLE]

In UI, the biased free energetics, $F_{v}^{\kappa,N^{*}}(\tilde{N})$ are further assumed to be parabolic, so that they can be estimated using the mean and the variance of the biased distributions. Estimates of $\partial F_{v}/\partial\tilde{N}$ obtained from different biased distributions are then combined using a sensible weighted averaging scheme, and integrated numerically to obtain estimates of $F_{v}(\tilde{N})$. Formulated in this way, UI was proposed to circumvent many of the early challenges associated with WHAM, such as the choice of bin width or the estimation of errors in $F_{\kappa,N^{*}}$ (these issues have since been addressed by more recent reformulations of WHAM, such as MBAR 67 or UWHAM 69.) By using only the biased means and variances, whose convergences could be readily ascertained, UI provided a useful alternative for analyzing umbrella sampling simulations.

In contrast with UI, sparse sampling seeks to estimate $F_{v}(\tilde{N})$ using a small number of short simulations that are sparsely distributed across the order parameter space. In this way, sparse sampling seeks to circumvent the overlap requirement intrinsic to umbrella sampling, and provide a computationally efficient alternative for estimating free energy landscapes. In sparse sampling, $F_{\kappa,N^{*}}$ is explicitly estimated by integrating its derivative, $dF_{\kappa,N^{*}}/dN^{*}$ with respect to $N^{*}$ (in contrast, $\partial F_{v}/\partial\tilde{N}$ is integrated with respect to $\tilde{N}$ in UI). Thus, we rely on a knowledge of the functional form of $dF_{\kappa,N^{*}}/dN^{*}$ on $N^{*}$, which is determined approximately at first, and is refined further using subsequent simulations. Finally, we note that as $\kappa\to\infty$, the two method converge; however, as discussed in the main text, errors in $dF_{\kappa,N^{*}}/dN^{*}$ increase with increasing $\kappa$, making it expedient to choose a small $\kappa$-value that is nevertheless large enough to result in monostable biased free energy profiles. To facilitate such an optimal choice of $\kappa$, sparse sampling makes use of estimates of $\partial F_{v}/\partial\tilde{N}$ at $\tilde{N}=\langle N_{v}\rangle_{\kappa,N^{*}}$.