Comparing generalized ensemble methods for sampling of systems with many degrees of freedom

Enhanced sampling refers to simulation methods that generate configurations that are not easily accessed from standard molecular dynamics (MD) trajectories, but which are needed in order to generate a meaningful statistical average.¹ The popular temperature replica exchange method (TREx) utilizes multiple replicas to generate configurations at higher temperatures to facilitate jumping between minima on the potential energy surface of a lower target temperature replica of interest.^2–4 Hamiltonian replica exchange (HREx) operates under similar principles and protocols to TREx, but using judicious modification of the Hamiltonian to increase access of the target replica to important local minima.^5,6 Overall, the REx approaches have been shown to improve sampling efficiency by increasing the rate of convergence of properties compared to averages accumulated over a standard MD trajectory at the target temperature or Hamiltonian.⁷ In addition, they have the attractive feature that tens to hundreds of replicas can be run in parallel to improve computational efficiency on CPUs, although REx is somewhat more cumbersome on GPUs.

Although the REx approaches are in practice, an improvement over straight MD, we have shown previously that in fact their sampling efficiency can be improved even further.⁸ Since the REx methods rely on multiple intermediate replicas to more efficiently pass trial moves to the target replica, that in turn can lead to diffusiveness in the exploration of configuration space; i.e., in the limit of an infinite number of replicas the swap acceptance probability would be 1, but exploration through parameter (temperature or forces) space would become a random walk. To address this problem, we have previously introduced a non-equilibrium alternative to TREx known as Temperature Cool Walking (TCW).⁸ TCW uses only one high temperature replica to generate non-local trial moves to ensure better ergodic sampling, while also being formulated to satisfy detailed balance to reach the correct limiting distribution appropriate to the lower target temperature. Furthermore, TCW correctly weights the configuration exchange from high temperature at any stage of cooling, thereby finding the sweet spot for optimal overlap with the target temperature distribution to increase the acceptance rates. The TCW algorithm was shown to have superior sampling capabilities compared to TREx on a model problem of a rough 1D potential energy surface for which ergodicity metrics can be analytically quantified.⁸ 

In this work we apply TCW to a more challenging set of applications of all atom simulations of peptides of increasing size in explicit water. We demonstrate its efficiency and ability to more rapidly converge to the proper limiting distribution when compared to standard TREx as well as compared to two altered Hamiltonian methods,⁹ Coulomb REx (CREx) introduced by Itoh and Okumura,¹⁰ and Dielectric Walking (DW) which is formulated in the same spirit as the original TCW approach but using the dielectric constant instead of temperature as the annealing variable. In the section titled “Cool Walking Theory,” we describe the non-equilibrium TCW and DW methods in some detail. In the section titled “Computational Details” we present the results of the various methods applied to a simple 1D potential, as well as the all-atom systems alanine dipeptide,^10,11 met-enkephalin,¹² and amyloid-β peptides,^13–15 all simulated with explicit water, where TCW is shown to be superior to the other methods. We conclude in the section titled “Discussion and Conclusion” with a discussion and summary of results.

Here we review the TCW method⁸ that in turn helps describe the new DW method introduced for the first time here, in which variations in temperature are replaced by variations in the dielectric constant. For TCW, only two explicit replicas are defined, corresponding to simultaneously propagating molecular dynamics (or Monte Carlo) at a high temperature T_H in which better ergodic sampling is ensured, in order to benefit the sampling convergence needed at the target state of low temperature T_L. The TCW algorithm uses a simulated annealing protocol to bring a trial configuration from the high temperature replica close enough to the distribution of the low temperature space to more readily accept the trial move. To formulate an acceptance criterion that satisfies detailed balance, we must perform an equal and opposite simulated annealing heating protocol on a low temperature conformation.

Algorithmically, we define a jump probability p_J that defines the frequency at which high temperature configurations will launch a statistical cooling run to serve as trial moves for the low temperature replica. When a random deviate ξ between 0 and 1 satisfies the condition that ξ < p_J, the current configuration for the high temperature replica, x_H, and low temperature replica, x_L, is stored. We continue to propagate the high temperature replica during the cooling of the configuration x_H to ensure that, in between successive statistical cooling attempts, the high temperature replica has sufficient time to fully decorrelate. We then perform simulated annealing on the configuration x_H, propagating it over a schedule of intermediate temperatures T_H−j, with j = 1 through H − 1.

We begin by capturing the probability for transitioning x_H from T_H to T_H−1, i.e., P_cw (x_H → x_H−1), by propagating x_H at this new temperature T_H−1 for a number of MD steps to generate a configuration x_H−1. If we desired to provide a trial move for the low temperature replica from the T_H−1 replica, the transition process for the first step of cooling, P_cw (x_H → x_H−1), is seen to simply be the ratio of the Boltzmann factors e^−βU(x) for x = x_H evaluated at T_H−1 and T_H. We then capture the next transition probability P_cw (x_H−1 → x_H−2) for cooling between T_H−1 and T_H−2, by evolving x_H−1 at the lower temperature T_H−2 to generate x_H−2, etc. Thus at any point in the cooling process, we can define a total transition probability, P_cw (x_H → x_H−j), which in general has the form

where each $Pqxr$ is the canonical Boltzmann probability of the configuration x_r in the ensemble at temperature T_q.

After the evaluation of every transition probability to the temperature T_H−j, we generate another random deviate ξ and compare it to the exchange probability p_E at the temperature T_H−j. If the random deviate ξ is less than p_E, we must similarly generate a transition probability for performing the reverse process, i.e., heating from T_H−j to T_H in order to satisfy detailed balance. To do so, we take the initial low replica configuration x_L and evolve it using MD with a temperature of T_H−j to generate a configuration $xH−j′$⁠. We use the notation x′ to describe the configurations whose temperature is increased from their original replica. The temperature of the $xH−j′$ configuration is then increased to T_H using the same schedule and process as used for cooling, ultimately yielding the configuration $xH′$⁠. This allows us to define a general transition probability for heating $Pcw (xH−j′→xH′)$⁠,

that allows us to an acceptance criterion that satisfies detailed balance.

Before we define the acceptance criteria, we note that we can further increase the acceptance rate by defining a window of states protocol that evaluates transition probabilities between volumes of phase space as opposed to points in phase space to minimize “bad” fluctuations that diminish acceptance rates. To describe the window of states component, we consider an intermediate temperature T_H−j from which we attempt to impose the configuration x_H−j onto the low temperature replica. From the cooling trajectory we accumulate the following window of weights:

and simultaneously for the heating trajectory

where the indices k through m in each summation refer to the individual conformations sampled. The two transition probabilities in Eqs. (1) and (2), as well as the two weighting factors in Eqs. (3) and (4), are combined in a standard Metropolis acceptance rule for imposing the cooled configuration x_H−j onto the low temperature replica, thus replacing the original configuration x_L,

In summary, the acceptance rule has three main parts: first, the simple Metropolis criterion for exchanging x_L with x_H−j; second, the ratio of transition probabilities for moving from temperature T_H to T_H−j and the reverse process; and finally the ratio of the window of states data.

If the move is accepted based on Eq. (5), we exit the cooling cycle and continue propagation of the two replicas, after updating the configuration of the low temperature replica to x_H−j. If the move is rejected, we continue by cooling x_H−j to the next temperature T_H−j−1 and continue through the cooling cycle until either an exchange is accepted or the cooling schedule completes without an exchange, at which point we return to propagation of the two replicas until the next CW cycle is attempted. We find it optimal to set p_E such that at least one exchange attempt is performed per CW cycle, on the order of 3%–10%. Setting p_E too high increases computational expense, as a new heating cycle must be performed for every exchange attempt during a cooling cycle, whereas setting p_E too low can result in a waste of a cooling trajectory if insufficient exchange attempts are made. Note that the window of states weighting is not necessary for the maintenance of detailed balance and does not have to be used, although we have found it to improve acceptance rates; addition of the window of states protocol was found to decrease the amount of annealing required by 25% in the original TCW paper on a simple 1D potential.⁸ 

Acceptance ratios for TCW moves tend to increase with more stringent schedules of the simulated annealing protocol used, and can readily achieve acceptance rates of 20%–30% or greater. Depending on the size of the system, we find a range of 25–50 fs of annealing at each temperature in the schedule to be sufficient for adequately equilibrating configurations to the lower temperatures as they are annealed. Given this, and the fact that every exchange results in imposition of a configuration that originated from the high temperature replica on the low temperature replica (as compared to REx, which requires many successive exchanges between intermediate replicas to “pass” a configuration all the way from the highest to the lowest temperature), p_J can be much lower in TCW than TREx, typically around 0.1–0.5 ps⁻¹. We have found that a ratio of 8 fs of additional high temperature replica propagation for decorrelation per 1 fs of cooling gives us the best increase in rate of convergence, though lower ratios are also workable.

We emphasize that TCW does not perform a complete configuration “exchange” as in standard REx protocols; instead we impose the annealed configuration from the high temperature replica onto the target replica, without imposing the configuration of the low temperature replica onto the high temperature replica. This assumes that the high temperature replica is an infinite reservoir of configurational states that decorrelate rapidly on the time scale between swap attempts. As we will see this is a good approximation for TCW, provided the maximum temperature is adequately high, but needs to be re-evaluated when the CW procedure is applied to other system parameters, such as modifying the Hamiltonian in our DW method.

In particular, what we have outlined in Eqs. (1)-(5) for TCW applies to a modification where temperature is fixed, but instead the Hamiltonian changes to define a series of dielectric constants that scale the permanent electrostatics (DW), allowing the protein to exit minima and pass through maxima caused primarily by electrostatic forces. In this case, the “ergodic sampling” replica could correspond to a higher dielectric than the target replica, which is decreased as we move x_H through annealing to x_H−j. Finally we can, in principle, also combine the TCW and DW approach (TCW-DW). In this method, the non-physical replica has both increased the temperature and altered protein dielectric, in which both decrease when moving the configuration x_H to x_H−j. If modification of the protein dielectric is indeed a viable way to increase sampling, when it is combined with temperature we might anticipate an increase in rate of convergence beyond that of standard TCW.

We have implemented TCW, DW, and TCW-DW into the OpenMM software package so that others can access the methodology presented here; we have also implemented the TREx and CREx methods in OpenMM in order to facilitate direct comparison. The methods have been tested on a number of systems presented here including a 1D potential,^8,16 as well as atomistic simulations of alanine dipeptide,^10,11 the 5-residue met-enkephalin peptide,¹⁷ and the 42-residue Aβ42 peptide,^13–15 all in explicit water.

Our first model system is a 1-dimensional potential energy function defined as

where the coefficients are chosen on the interval [ − 1, 1], and the length of the simulation box in reduced units is L = 10. The coefficients are the same as used in the original papers,^8,16 but are given in Table S1 with a corresponding plot of the potential energy (Fig. S1) provided for the benefit of the reader. For this simple system we have used a hybrid Monte Carlo (HMC) method to propagate sampling on its energy surface, which is described in the original TCW paper. For DW we use only two replicas which are propagated by HMC at values of ε = 30 and ε = 1, and then trial moves are generated for the low dielectric replica using a (roughly) geometric schedule for a statistical unscreening process (ε = 30, ε = 29, ε = 24, and ε = 1) and using a similar geometric schedule for the screening process (ε = 1, ε = 3, ε = 10, and ε = 30) needed for detailed balance. For CRex, we found 4 replicas to be optimal, with a geometric schedule of dielectric replicas (ε = 1, ε = 3, ε = 10, and ε = 30). All HMC parameters and swap attempt rates that we test for the DW and CREx methods are the same in the comparison, and the DW and CREx methods differ only in the number and the way in which replicas interact.

We also consider three successively larger peptides in explicit water to evaluate the TCW, DW, and CREx methods, to see whether they actually perform better than the standard TREx protocol. For alanine dipeptide, we use the Amber ff99sb force field¹⁸ for the protein and TIP3P for the water model¹⁹ to compare to previous studies^10,11 that use this combination of force fields. For met-enkephalin, we use the same combination of force field and water model, again to compare to previous studies.¹⁷ For Aβ42 we used the Amber ff99sb force field¹⁸ and TIP4P-Ew water model²⁰ since we have published results^13–15 using TREx with this force field combination. All simulations used a 1 fs time step with SHAKE constraints to freeze out hydrogen vibrations, while an Andersen thermostat maintained the temperature. Ewald was used for calculating long-range electrostatic forces, with a cutoff of 9.5 Å for the real space electrostatics and Lennard-Jones forces. The replica exchange attempt frequencies were once per 1 ps for alanine dipeptide and once per 500 fs for met-enkephalin and Aβ42. For alanine dipeptide and met-enkephalin, we used the LEaP module to prepare the initial extended structure of the peptide, and then OpenMM to solvate it within a periodic cube of water appropriate to the system size (233 and 499 water molecules for alanine dipeptide and met-enkephalin, respectively). Using these starting states we equilibrated the temperature or dielectric replicas for 100 ps, and generated multiple independent production runs of 50-100 ns. For Aβ42, we used a pdb file generated from TREx simulations previously run in our lab, and it includes three sodium ions to neutralize charges on the peptide within the box. The initial structure used for Aβ42 is shown in Fig. S2.

We first consider the performance of CREx and DW on a rough one-dimensional energy surface, which we have done previously for TCW and TREx.⁸ This simple model system has proved useful to confirm detailed balance of the implemented code and to measure the sampling efficiency, i.e., the time required to reach convergence as measured by the ergodicity factor, χ(t),

where

and ε is the relevant dielectric constant, and Z is the configuration integral, and in which the exact probability distribution is analytically solvable. Figure 1 plots χ²_1D(t) vs. time for the 1D case, in which it is evident that the DW algorithm converges more quickly to the correct limiting distribution compared to the CREx method, and once again illustrates that any replica exchange approach (temperature or Hamiltonian) can be plagued by problems of intermediate replicas that retard the movement of configurations generated by the most ergodic replica from reaching the target replica.

We next consider the alanine dipeptide in explicit solvent as the first fully atomistic test of our TCW and DW approaches, and compare their performance to TREx and CREx. Given the relatively small conformational space defined by the ϕ and ψ dihedral angles of the dipeptide, we make the assumption that two independent microsecond long standard MD trajectories, using the same force field, will serve as an exact benchmark for testing the rate of convergence of the four sampling methods, since it should be a sufficient amount of time for the peptide to fully cover its conformational space. A good first test of the different enhanced sampling methods is to access the relatively favorable left-handed α − helix (ϕ = 60°, ψ = 50°) that must overcome barriers of ∼10-15 k_BT from the right-handed α − helix (ϕ = − 60°, ψ = − 50°), extended–β (ϕ = − 180°, ψ = 180°), or the polyproline II (ϕ = − 150°, ψ = 150°) conformations. This is an interesting test case for the dielectric methods in particular since the barrier is electrostatic in origin, i.e., in which the two peptide oxygens need to come into close contact to execute the transition.

For TREx and CREx, the initial spacing of temperature and dielectric replicas, respectively, was the same as that published in the CREx paper.¹⁰ However, we found that the temperature replica spacing was non-optimal and puts TREx at a significant disadvantage compared to CREx. Therefore, we adjusted the maximum temperature and replica spacings to obtain more standard acceptance rates of ∼20% for configurational exchanges between replicas using TREx (Table S2). In both cases, replicas were run for 50 ns each. For TCW, the maximum temperature was set to be the same as that for TREx and followed an annealing schedule that was largely the same as the TREx replica spacings. For DW, the maximum peptide dielectric was not the same as that for CREx, which was infinity, but instead it was set to 16; repeated testing of DW did not demonstrate any significant difference once the maximum dielectric was set at or above this value. The “descreening” schedule used for DW closely followed the same replica spacings used in CREx. The two replicas for TCW and DW were also run for 50 ns each. Replica spacings and all exchange probabilities for TREx and CREx, as well as cooling and heating schedules for TCW and DW, are presented in Table S2.

Figure 2 shows the free energy plots for the alanine dipeptide in which it is evident that all 4 methods qualitatively reproduce the result generated from the exhaustive MD benchmark. To measure the rate of convergence, the Ramachandran plots are discretized into a 73 × 73 grid, and the square root of the total sum of squared differences between the Ramachandran free energy values, A, of the enhanced sampling techniques over time relative to our reference MD Ramachandran plot is calculated,

and normalized such that it varies between 0 and 1. Five independent trajectories of each kind of enhanced sampling method were run to obtain a measure of the average performance for the $χADt$ metric for alanine dipeptide. Figure 3 shows the convergence profile over 50 ns for each method. The TCW method clearly converges the fastest, reproducibly dropping below a $χADt$ value of 0.1 within the first 10 ns of simulation, while TREx is noticeably more slowly convergent and variable between trajectories. The DW and CREx methods are roughly competitive with each other, and exhibit significantly greater variance in their rate of convergence depending on the trajectory compared to the temperature methods. We attempted multiple variations of the DW protocol, including shaping the schedule differently, or having both replicas run unperturbed but have the peptide dielectric increase and then decrease back to one through the course of the annealing schedule. Neither resulted in improvement over the published data here. We also combined DW with TCW to see if the additional alteration of the temperature would improve the rate of convergence relative to standard DW, but observed no significant improvement relative to TCW. We determine at this stage that the CREx and DW approaches are not competitive to the temperature enhanced sampling methods, even for problems dominated by electrostatic barriers, and we move forward on only testing TCW and TREx on successively larger all-atom systems. We return to this point in the discussion. Figures S3(a) and S3(b) are trajectory movies for representative tests of TCW and TREx, respectively. All molecular renderings and movies were generated using the visual molecular dynamics package.²¹ 

We next compare TREx and TCW on the met-enkephalin peptide. The temperature schedule for TCW and temperatures and acceptance rates for exchanges between replicas for TREx are given in Table S2, and we ran 5 independent 50 ns trajectories for each method. Met-enkephalin, having many more degrees of freedom, would be more difficult to fully sample from a standard MD trajectory to serve as a gold standard for convergence. Instead, we measure the rate of convergence of the independent trajectories to each other; if the trajectories are initiated from very different initial conditions, the expectation is that any significant difference would indicate that a given method suffers from inferior sampling efficiency over the 50 ns period. Trajectory movies for TCW and TREx are included in Figures S4(a) and S4(b), respectively.

We begin with the measured variance between different trajectories for free energy values generated across the Ramachandran plots for each of the 5 residues in met-enkephalin. However, since higher deviations are likely to occur in the high energy regions where sampling is poorer, we have plotted a normalized version where the standard deviation at a point is divided by the average free energy at that point. This will more heavily weight the lower free energy regions of the Ramachandran plot, and thus the degree to which the 5 independent trajectories converge to a stable free energy value for the most important minima. Figures 4 and 5 show the two normalized standard deviation Ramachandran plots for both methods for representative residues tyrosine-1 and methionine-5, respectively; all additional Ramachandran plots for the other amino acids are available in the supplementary material (Figures S5-S7). Across all residues for met-enkephalin, the TREx methods exhibit higher standard deviations than for TCW. This indicates that the free energy of the low energy basins is not as well converged after 50 ns for TREx, whereas for TCW there is clearly convergence to a stable free energy value in each of the basins.

We can also get a better sense of global sampling of the met-enkephalin peptide conformations using a principal component analysis. In this case we grouped all five trajectories together to generate an average conformation, center all conformations to this average, redefine all conformations as deviations from the average, and then calculate the eigenvectors that describe the major modes of the collective ensemble. Each individual conformation of a given trajectory is then projected onto the principal modes, and in Figure 6 we plot these projections onto the first principal component against the projections onto the second for two trajectories each of TCW and TREx. All additional plots for the remaining trajectories are available in the supplementary material (Figures S8 and S9). While the TREx trajectories exhibit very limited, localized sampling of a small and closely related group of conformations, conformational sampling is evidently more uniform in the TCW tests. Though TREx is likely identifying the local minimum states, as evidenced by the high probabilities for certain regions, the strong disagreement over which the states are highly probable between independent trajectories indicates that each trajectory spends a disproportionate amount of time trapped in a few minima.

The PCA results can also be used to develop a time-dependent metric for the rate of convergence of two independent trajectories for each method. The Kullback-Leibler divergence (KLD) quantifies the difference between two probability distributions, and we use it here to calculate the degree of difference in the exploration of principal component space between pairs of independent trajectories over time. If two trajectories reproducibly cover the same regions in phase space, the KLD will decrease toward 0. We calculated the KLD for each pair of trajectories, monitoring their overlap in sampling of the first two principal components.

The KLD values for the two principal components are plotted in Figure 7. For the first principal component, PC-1, the greatest change occurs over the first 30 ns, indicating that much of phase space is covered by then, with steady convergence to the global metric for the remaining 20 ns (Figure 7(a)). The KLD metric converges faster for the TCW method, decreasing to a final value of about 0.015 over the 50 ns sample, indicating that the independent trajectories are thus very close to converging to each other with respect to PC-1. In contrast, the TREx result converges somewhat more slowly in the first 30 ns, and it remains at a higher value after 50 ns relative to TCW, indicating that the independent trajectories for TREx are farther from convergence for PC-1 than that found for TCW. Similar behavior is evident for PC-2 in which the KLD metric for TCW much more rapidly attains lower values than is observed for TREx (Figure 7(b)).

To illustrate that the TCW method is computationally tractable for much larger and more difficult systems, we have evaluated the structural ensemble of the Aβ42 IDP in explicit water. Because Aβ42 is an intrinsically disordered peptide, its structural ensemble requires enhanced sampling techniques. We have previously used TREx using 58 replicas, running each replica out to 50 ns across a large number of CPUs to take advantage of fine-grained parallelization of energy and forces, to create two independent structural ensembles of the Aβ42 peptide.^13–15 For the Aβ42 peptide sampling reported here, we extended these two independent TREx trajectories to 100 ns, to which we compare to two independent 100 ns trajectories with TCW on a single GPU. Details of the replicas and temperature schedule are reported in Table S3. Trajectory movies from TCW and TREx are linked in Figures 8(a) and 8(b) (Multimedia view), respectively. Figure 9 compares the radius of gyration (R_g) distribution in which we see that the Aβ42 monomer R_g is more extended for the TCW ensemble than that observed for the TREx ensemble. Secondary structure propensities for each residue are presented in Figure 10 in which we consider helices (α, 3-10, and π in Figure 10(a)), all β content (bridges, hairpins, and sheets in Figure 10(b)), as well as all turn types (Figure 10(c)). Overall there is a qualitative change in the distribution of secondary structure propensities, in which there are more localized regions of helical and turn propensities and much less β content in the TCW ensembles, consistent with the R_g trends. We will report on the comparison to recent NMR data obtained on Aβ42 in a future publication.

We have compared four different methods for their ability to enhance sampling and thereby increase the rate of convergence to the limiting structural ensemble, two of which use temperature to overcome conformational barriers^2,8 and two that alter the potential energy surface by scaling the electrostatic charges.¹⁰ We first tested the methods on two simple systems, an analytic rough 1D energy surface and an all-atom alanine dipeptide in water simulation, of which the relatively small phase space allows for the ability to define quantitative convergence metrics. While the DW method is shown to be better than CREx for the 1D potential, the DW and CREx methods perform poorly with respect to the temperature methods for the alanine dipeptide, even though this should be a good case for both dielectric-based methods where the dominate barrier is electrostatic in origin. In essence, screening of electrostatics merely gives rise to new potential energy barriers that inhibit ergodic sampling, and we can only conclude that the dielectric CREx and DW methods are far inferior to standard TREx.

In this quantitative comparison for alanine dipeptide we also show that TCW is superior to TREx, and this outcome was shown to hold for met-enkephalin¹⁷ based on the KLD metric that measures the convergence of independent trajectories to converge onto the two primary principal components of the peptide. Lastly, we compared the results of TCW and TREx on Aβ42, the Alzheimer’s peptide, known for its role in Alzheimer’s disease,²² in which we find non-negligible differences in the disordered ensemble compared to standard TREx.

Finally, we consider the computational efficiency of the two temperature methods. The TREx approach requires tens of replicas and a CPU cluster with a real communication backbone to perform swaps and can be made highly efficient on CPUs. In contrast, TCW always requires only two replicas and is better suited to GPU systems. In order to disseminate our new enhanced sampling approaches applied to complex systems such as peptides and proteins in water, TCW and DW have been made available through the OpenMM Omnia software consortium;²³ for comparison, the TREx and CREx methods have also been made available through Omnia as well. The current GPU implementation, however, does not have true parallelization for the propagation of the high and low temperature walkers simultaneously, and we are currently working on a parallel implementation in OpenMM so that computational efficiency does not inhibit the excellent sampling efficiency of the TCW method.

See the supplementary material for 1D potential parameters (Table S1) and probability (Figure S1), replicas, exchange probabilities, schedules for all methods applied to alanine dipeptide (Table S2), replicas, exchange probabilities, schedules for TREx, and TCW applied to Aβ42 (Table S3). Initial structure used for Aβ42 production (Figure S2). Trajectory movies for alanine dipeptide (Figure S3) and met-enkephalin (Figure S4). Normalized standard deviation between trajectories quantified for the Ramachandran phi, psi angles for Gly-2, Gly-3, Phe-4 (Figures S5-S7). Principal component analysis tests for three additional TCW runs (Figure S8) and TREx runs (Figure S9).

T.H.G. acknowledges support from the NSF under Grant No. CHE-1363320. J.L. acknowledges partial support under the NIH Molecular Biophysics TG, T32-GM008295. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.