{'fraction_non_alphanumeric': 0.030529229220833683, 'fraction_numerical': 0.013671055942296402, 'mean_word_length': 4.735449735449736, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Population annealing molecular dynamics with adaptive temperature steps

Population annealing is a novel generalized-ensemble simulation scheme used in large-scale parallel Monte Carlo simulations of disordered spin systems and similar problems. In a recent publication we proposed a generalization of this method to molecular dynamics simulations of biopolymers. In the present article we review this work and introduce a scheme for automatically choosing the temperature steps based on the observed distribution of potential energies in the population of simulated replica.


Introduction
Due to a rugged free-energy landscape, studies of the folding behavior of peptides are among the computationally most difficult problems in computer simulations. In such problems, the systems get trapped in local minima that are separated by high barriers in free energy that are increasingly difficult to overcome by activated dynamics as the temperature is gradually lowered. Several numerical approaches have been developed to overcome such trapping problems and the ensuing loss of ergodicity in simulations of biopolymers. Parallel tempering (PT), sometimes also known as replica exchange, has previously been shown to have the potential of successfully sampling the full configuration space when combined with molecular dynamics (MD) simulations of peptides [1,2]. In PT, one considers a number of replicas (≤ 10 2 ) of the same system, which are, a priori, simulated independently at different temperatures, ranging from unphysiologically high temperatures where configurations equilibrate quickly, forgetting their previous history, down to the physiologically relevant temperatures, where the peptide is in a folded state. The simulations running at neighboring temperatures are subsequently permitted to exchange configurations according to an acceptance probability taking into account the Boltzmann weight at the two temperatures. This setup results in a diffusive motion of replicas in temperature space, thus allowing the replicated system to escape from local free-energy minima by heating copies up to the highest simulated temperatures before they cool down again, typically arriving in a different local minimum from the one occupied previously. Although PT can easily be implemented in parallel, the number of processors which can effectively be used is limited, as a setup with too many closely spaced temperatures does not improve the equilibration and even slows down the random walk in temperature space. In contrast, population annealing (PA) [3,4,5,6] is a simulation scheme introduced in the context of Monte Carlo (MC) simulations that can make use of very large numbers of parallel workers and hence physical cores and massively parallel accelerators such as GPUs [7] without ensuing problems of bad parallel scaling or decreasing simulational efficiency. In this method one considers a multitude (typically at least 10 4 ) of independent replicas which are initially equilibrated at a high (or potentially even infinite) temperature to get an ensemble of configurations that is perfectly or nearly perfectly uncorrelated. The temperature is then lowered in a sequence of steps, similar to what is done in simulated annealing for a single replica [8]. However, in contrast to simulated annealing, the population is additionally resampled to the new temperature after a number of simulation steps. This resampling leads to a significant speed-up in equilibration and eliminates energetically stuck configurations (reminiscent of the "go with the winners" strategy of Ref. [9]). Observables in PA are calculated using an average over the ensemble, whereas in PT averages are taken over the time series at a fixed temperature. Population annealing has been shown to perform well, e.g., in the simulation of spin glasses [4,6]. We recently showed how PA can be combined with MD (PAMD) to efficiently simulate macromolecular systems such as the folding of the peptide met-enkephalin [10]. In the present work we demonstrate how the schedule of temperature steps in this scheme does not need to be pre-determined, but it is much preferable to decide about the next temperature to include as one goes along by utilizing the information about the distribution of potential energies resulting from the parallel simulation of a population of replicas.

Model and method
The algorithm for PAMD as described in Ref. [10] comprises the following steps: (i) Set up an equilibrium ensemble of R independent copies of the system at some high temperature T 0 . (ii) Choose the next temperature T i from a pre-defined sequence. (iii) Resample the ensemble of systems to the new temperature T i < T i−1 by replicating each copy a number of times proportional to the relative Boltzmann weight In our previous work [10], studying the folding of the penta-peptide met-enkephalin with the amino-acid sequence Tyr-Gly-Gly-Phe-Met, we adapted a temperature set (T i = 700, 585, 489, 409, 342, 286, 239, 200 K) from a PT simulation of the same peptide [2]. However, finding an appropriate temperature set in PT simulations is a non-trivial task, for which a number of solutions are being discussed [11,12,13]. It is clear, for instance, that too large temperature steps in PT lead to low acceptance rates of swap moves, whereas too many small steps lead to a slowing down resulting from the scaling properties of the temperature random walk. Similar difficulties arise in PAMD with a fixed temperature protocol. In population annealing, however, it is possible to exploit the fact that one has access to the energy histogram at temperature T i−1 from the population of replicas at the point of having to decide about the next step, T i , in the temperature scheme. To benefit from this additional information, we apply single-histogram reweighing [14] to get an estimate of the energy histogram at T i . We then require that the histograms of potential energy at T i and T i−1 have a certain prescribed overlap α at each step. The actual overlap can be estimated before taking the real temperature step from the following expression [7], Here E j is the energy of the jth replica and serves as a normalization. Before the simulation one then only has to choose a fixed value of α = α , and any given value of α results in a certain temperature sequence. In the above algorithm, we hence propose to replace step (ii) by the following adaptive variant: (ii') Find the next temperature T i in such a way, that α(T i−1 , T i ) = α .
The question of what are reasonable values of α is investigated in the following section.
In the simulations reported below we used OpenMM [15] to perform the MD, relying on the Langevin thermostat and a velocity-Verlet like integrator. The interactions of the peptide are modeled via the AMBER force-field ff94 [16]. It is worthwhile to point out that for PAMD simulations to work properly the thermostat has to be stochastic as otherwise all of the resampled copies of a population member will follow exactly the same trajectory in further evolving the system, thus defying the decorrelating purpose of the resampling step.

Results
We fixed the relevant PAMD parameters at R = 10 000 replicas, a time step of dt = 0.5 fs, and a friction coefficient γ = 1/ps. For the constant energy overlap simulations, we chose T 0 = 700 K as the initial high temperature and set T N = 200 K as the target temperature. The start configurations at T 0 were taken from a canonical simulation of length 25 ns, followed by short simulations for each replica to further decorrelate the population. To find a suitable set of overlaps α and corresponding number θ of MD simulation steps per temperature, we tested α = 0.3, 0.5, 0.8 and θ = 2000, 4000, 6000, 8000. In Fig. 1(a)  As is evident from Fig. 1(a), α = 0.3 leads to a temperature set which is very close to the one used in Refs. [2,10]. Thus the simulation with α = 0.3 and θ = 4000 is close to equivalent to our previous study in [10], and it is hence used as a benchmark. In the regime of sufficiently large θ to ensure equilibration in between temperature steps, the adaptive stepping algorithm is insensitive against the chosen θ and mainly depends on α . This is also reflected in Fig. 1(b), where we plot the measured temperature T m against the externally imposed (thermostat) temperature T . Also plotted is a solid line, corresponding to the expected behavior in properly equilibrated simulations, namely T m = T . Note that this type of plot is mostly used in MD simulations, where the intrinsic temperature is readily defined in terms of the kinetic energy, while in MC simulations proper equilibration is mostly ascertained by checking the stationarity of time series of suitably chosen observables.
In Fig. 2 the histogram of potential energy at T = 200 K is presented for the different values of α and θ used in Figs. 1(a) and (b). The histograms were obtained by means of the weighted histogram analysis method (WHAM) [17]. One advantage of PAMD simulations is, that one can easily seed the WHAM with good estimates of the partition functions Z, because [5,7] This allows one to only perform very few iteration steps in the WHAM while still getting very good estimates of the partition functions. From Fig. 2 it is clear that all considered simulation protocols produce an energy histogram that is well equilibrated.
The quality of sampling of configuration space can be more readily judged from Ramachandran plots, showing the dihedral angles φ and ψ along the backbone of the peptide. Such plots are shown in Fig. 3 for the center amino acid GLY3. Panel (a) shows the result of a PT simulation with the same parameters as used in Ref. [10] for reference, panels (b) and (c) correspond to PAMD simulations with adaptive temperature steps using α = 0.5, θ = 4000  (b) and α = 0.8, θ = 2000 (c), respectively. All three simulations show broadly comparable sampling of the configuration space. It appears from a comparison of (b) and (c), that it is favorable to rather choose a larger number of steps θ than increasing the overlap α . A similar improvement as in (b), however, is achieved by increasing the population size R, see Fig. 3(d). However, this comparison does not take into account the computational resources used. Relative to the PT simulation shown in (a), the adaptive PAMD simulations in (b) and (c) used 1.525 and 1.875 times more computing cycles, respectively, whereas the non-adaptive but larger PAMD simulation shown in (d) also used 1.875 times the cycles of (a). The sampling of (b) and (d) is in both cases slightly superior to (c). One should keep in mind that the wall-clock time in PAMD with increasing number of replicas R can be kept very low when sufficient parallel resources are available. The minimum wallclock time in PA is predominantly limited only by the short simulation time t s needed for a single replica where t is the total simulation time.

Conclusion and Outlook
We have performed population annealing molecular dynamics (PAMD) simulations with adaptive temperature steps for the folding of the penta-peptide met-enkephalin. It was shown that using the criterion of a constant energy overlap significantly simplifies the task of choosing a proper temperature set, a notorious problem also encountered in parallel tempering simulations. We find reliable results in a wide range of possible overlap values α with similar computational efforts if the number θ of MD steps is correspondingly adapted. In the regime where θ is large enough, the precise value of θ has little influence on the actual schedule proposed. If additional parallel computational resources are available, it seems to be preferable to use those to simulate more replicas in the population rather than increasing θ or α , but the precise optimum in this multi-dimensional parameter space will depend on the system under consideration.