{'fraction_non_alphanumeric': 0.03130381303813038, 'fraction_numerical': 0.007011070110701107, 'mean_word_length': 4.469559367642113, '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': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Enhanced sampling of transition states

The free energy landscapes of several fundamental processes are characterized by high barriers separating long-lived metastable states. In order to explore these type of landscapes enhanced sampling methods are used. While many such methods are able to obtain sufficient sampling in order to draw the free energy, the transition states are often sparsely sampled. We propose an approach based on the Variationally Enhanced Sampling Method to enhance sampling in the transition region. To this effect, we introduce a dynamic target distribution which uses the derivative of the instantaneous free energy surface to locate the transition regions on the fly and modulate the probability of sampling different regions. Finally, we exemplify the effectiveness of this approach in enriching the number of configurations in the transition state region in the cases of a chemical reaction and of a nucleation process.


Introduction
One of the crucial problems in the atomistic simulation of systems of chemical-physical relevance is the study of rare events. Typical rare events are chemical reactions, phase transitions, and biomolecule conformational changes. These events take place on a time scale that cannot be reached by straightforward molecular dynamics (MD) and require enhanced sampling methods to be studied. A vast class of such methods 1-7 relies on the identification of an appropriate set of collective variables (CVs) whose fluctuations are then suitably enhanced. The CVs, here indicated as s = s(R) , are functions of the atomic coordinates R, and are meant to encode the slow degrees of freedom of the system. In our group we have focused on two such approaches, namely Metadynamics 2 and Variationally Enhanced Sampling (VES). 6 One of the outputs of these calculations is the free energy F (s) as a function of the chosen CVs. With adroitly chosen CVs, the F (s) separates the metastable states and the apparent TS does coincide with the real one. Luckily much progress has been made in finding 8 and ameliorating CVs 9-11 and we shall assume here that we are in a situation such that the CVs pass close to the real transition state (TS).
However, even with good CVs, in a standard Metadynamics or VES calculation the transition state is not as well sampled as one would like it to be. In fact in these methods based on a dynamical approach in which the F (s) is progressively filled by a bias, most of the sampling effort is spent in exploring the F (s) minima and the transition states (TS) are sampled only when the system moves from one basin to another. This is an event that, while accelerated by the bias, is comparatively rare. Given the relevance that the transition state plays in understanding the process under study and eventually guiding it towards a desired end, it would be of great help to have a method that partitions sampling more evenly between the different stationary points of F (s) and harvests a much larger set of transition state configurations.
In order to achieve this result we shall use one of the most attractive features of VES, namely the possibility of searching for a bias such that the distribution in the biased ensemble is equal to a pre-assigned one. This leads to an enhanced sampling of the transition state region whose properties can then be studied in great detail. We apply the method thus developed to the study of a bimolecular nucleophilic substitution (S N 2) reaction and to a liquid-vapor transition, which is a prototypical nucleation process.

Variationally Enhanced Sampling
In VES, the bias potential V (s) is obtained by optimizing the functional where T is the temperature, k B is the Boltzmann constant, β = 1/k B T and p(s) is a target distribution which is to be normalized. It has been shown 6 that Ω[V (s)] is a convex functional and at its minimum the value for the bias V (s) is, within an irrelevant constant, From equation (2) it follows that at the minimum the CV distribution in the biased ensemble, is equal to the pre-assigned target distribution p(s).
Optimization In order to use the variational principle and minimize the functional, the bias V (s) is expressed in terms of a set of parameters α = {α 1 , α 2 , α 3 , ...}. This turns the functional Ω[V (α)] into a function Ω(α). This function Ω(α) is then optimized using the gradient ∂Ω(α) where the terms · · · V (α) and · · · p are expectation values calculated in the ensemble biased by V (α) and over the distribution p(s) respectively. However, the gradient evaluated from such a statistical calculation is inherently noisy and thus, a stochastic optimization algorithm 12 is necessary to obtain a faster and smoother convergence.
Here, we take the bias potential V (s; α) to be a linear expansion in a finite number of basis functions, and use the Bach and Moulines averaged stochastic gradient descent algorithm 13 to update iteratively the parameters of the bias expansion.
Target distribution A straightforward choice of p(s) is the uniform distribution, that at convergence gives V (s) = −F (s). However, this distribution forces the system to visit with equal probability all points in the CV space, including uninteresting high energy regions, often resulting in poor convergence. Although it is possible to construct a p(s) that focuses on important regions of the CV space, this requires apriori information on F (s), a quantity that is in general unknown at the beginning of the calculation. However, Valsson and Parrinello 7 introduced a self-consistent approach to include F (s) in the definition of p(s). They suggested a p(s) that at convergence results in, where γ ≥ 1 is a parameter called bias factor. Unlike VES with a fixed p(s), this approach requires two coupled iterative schemes. In addition to updating the parameters for V (s; α) every ∆t V MD steps, the p(s) is updated every ∆t TD number of iterations of the optimization of Ω(α). At first, p (0) (s) is typically taken to be a uniform distribution and the first guess for F (s) is zero. The F (s) at any given step (k + 1) is then obtained as, In the asymptotic limit, one has V (s) = −(1 − 1 γ ) F (s) modulo a constant as in the welltempered ensemble.
It has been shown that, even if the instantaneous estimate of F (s) is initially inaccurate, VES quickly converges to the correct F (s). 6,7 Additionally, in many practical cases, the well-tempered variant was found to converge much faster than a VES simulation that uses a uniform p(s). In figure 1, we illustrate with the case of a one-dimensional double-well potential of barrier height 40 k B T , the probability density of the CV obtained with uniform and well-tempered p(s) as compared to the unbiased Boltzmann probability distribution. It can be seen that in a well-tempered simulation, the system visits the TS less frequently than in a uniform target simulation. However, as discussed above one has to pay the price of a slower convergence.

Targeting the transition state
For simplicity sake, we deal here only with cases in which the CV is unidimensional. However, an extension to a multidimensional set of CVs is possible. For a well chosen CV, the TS is located at the apparent F (s) TS. Hence, sampling of the TS region can be enhanced by designing a p(s) that gives a higher probability of sampling in a region around the apparent TS.
To do so, we utilize the instantaneous F (s) similarly to what is done in well-tempered VES. We take into account the fact that the apparent TS is a F (s) stationary point, where dF (s) ds = 0. A p(s) constructed to enhance the probability of sampling the stationary points should, in principle, be able to boost TS sampling. The most obvious choice of such a distribution would be a Gaussian-like distribution similar to the well-tempered p(s) that has as argument dF (s) ds . However, we found it difficult to describe the rapidly varying tail of such a function with a finite basis expansion unless γ was taken so large that the resulting distribution was practically uniform. Hence, we ameliorated the behaviour at the tail by using a Lorentzian-like function that decays polynomially where Z p is a normalization factor, and ζ is a parameter of units s −1 that regulates the broadening of the distribution. This p L (s) will weight uniformly all the stationary states of F (s): minima, TS, and if present, inflection points. However for the purpose of our calculation, it is desirable to give more weight to the TS. In order to do this we take advantage of the fact that at a local maximum not only dF (s) ds vanishes but also the second derivative satisfies the condition d 2 F (s) ds 2 < 0. For this reason, we multiply p L (s) in equation (8) where 0 ≤ a ≤ 1. This function goes from the value S a,λ (s) = 1, when d 2 F (s) In figure 2, we show the effect of changing the parameters a and λ in the symmetric double well case.
In order to obtain these results and the others reported below we had to solve the problem of computing the first and second derivatives of F (s). If not done with care this can introduce unwanted noise. To this effect we expand F (s) in the same basis set {f k (s)} used to express V (s) (see equation (5)). In this way, the s dependence appears only in the functions f k (s) whose derivatives can be computed analytically.
We now illustrate how this approach works in practice on the paradigmatic cases of a bimolecular nucleophilic substitution reaction and the condensation of a Lennard-Jones gas.

Bimolecular nucleophilic substitution reaction
Nucleophilic substitution reactions constitute a popular class of chemical reactions. The asymmetric substitution reaction of fluoromethane and chloromethane CH 3 F + Cl -CH 3 Cl + Fis a simple example of such reactions and has been well studied over the years. Following a previous enhanced sampling study of the same reaction, 14 we chose the difference between the C-Cl (d Cl ) and C-F (d F ) atomic distances as CV.
Here, we contrast four different simulations that use VES with different target distributions. The first is a well-tempered one with γ = 50, the second uses a p(s) that is uniform in the interval −2.0 ≤ s ≤ 2.9. The third simulation uses p L (s) as defined in equation (8)  We first converge the bias coefficients (α) and use the resulting static potential to perform the simulation. We then obtain the F (s) using an umbrella-like reweighting. As expected

Liquid-vapour phase transition
Nucleation is a phenomena of interest in a broad range of fields. A commonly used model in homogeneous nucleation studies is that of the Classical Nucleation Theory (CNT). It states that the energetics of growth of a stable nucleus (or cluster) in a metastable bulk is shaped by an interplay of two factors. While the creation of the interface between this cluster and the bulk phase hinders the cluster growth, the formation of these more stable regions favours it. These two terms are proportional to the surface area and the volume of the cluster respectively. As a result of these competing contributions, cluster growth is hindered until a critical size is reached. Beyond this critical size, the stability gained from the volume term surpasses the cost of interface formation and the cluster grows spontaneously. Due to the existence of this barrier, it is often difficult to determine the critical cluster size 15 16 and obtain sufficient sampling in the region even using enhanced sampling techniques.
Here, we compare the results of our TS-target VES simulation with uniform VES simulation in the case of condensation of Argon. The details of the simulation can be found in the SI. Following Piaggi et al., 17 we use a CV defined as s = (n l ) 1/3 where n l is the number of liquid-like atoms. In the case in which there is only one cluster present s is roughly proportional to the radius of the cluster. We further restrict the growth of clusters beyond 120 atoms using a harmonic restraining potential having κ = 0.1. The parameters for the uniform VES simulation are similar to those of Piaggi et al. 17 while that of the TS-target are a = 0.02, λ = 0.1 and ζ = 30. Since we want to focus predominantly in the critical region we use a very small a. Here also, we first converge the bias coefficients (α) for the two cases using 10 walker simulations and then perform a single walker simulation with a static bias potential. The F (s) is obtained by an umbrella-like reweighting and is shown in figure 4.
Further, we analyze the shapes of the largest clusters obtained during the simulations.
To do so, first we identify the largest clusters of liquid-like atoms obtained from both the simulations using the depth first search clustering algorithm described in Ref. 18. Clusters spanning more than half the box size and clusters with less than 6 atoms are rejected. We  We diagonalize the matrix to obtain the three eigenvalues (λ 1 ≥ λ 2 ≥ λ 3 ) and subsequently calculate the shape anisotropy value defined as k = 3 The value k = 0 indicates a perfectly spherical cluster whereas the value k = 1 corresponds to a cluster where all the particles are linearly arranged. Figure 5 shows how even clusters with the same number of particles could have very different shapes and hence, a range of anisotropy values.
In figure 6, the probability density of the CV space is shown. The probability density has been estimated using multivariate gaussian kernels. It can be concluded from figure 6 that the proposed target distribution has been able to sample more effectively the critical cluster size. With the uniform target distribution, the system tends to sample more the uninteresting region with very small clusters due to the low energy of the region while with the TS-target distribution, the high-energy region around the critical cluster has been sampled more in the same computational time. Also, With the TS-target, the shape of the critical cluster is thus estimated more accurately. Only clusters having at least 10 atoms have been considered in the plots it is evident from both the figures that as the cluster size increases the clusters tend to be more spherical. Hence, while the results obtained conform to previously observed trends, we are now able to quantify more confidently about the critical cluster properties.
Last but not least in this case an improvement in convergence rate is observed. In order to quantify the convergence of F (s), we use an error metric similar to that used by Valsson and Parrinello. 7 whereF Multiple F (s) are obtained by reweighting the Uniform and TS-target 10 walkers VES simulations at regular intervals using the reweighting scheme proposed in Ref. 19. The reference F (s) is obtained by an umbrella-like reweighting of the simulation that used the converged static bias obtained from the 10 walker Uniform VES simulation. It can be seen from figure   7 that the TS-target does converges faster.  Figure 7: The error (ε) calculated using Eq. (11). Multiple F (s) were obtained by reweighting 10 walker simulations at intervals of 5 × 10 4 MD steps. The reference F (n 1/3 ) was obtained by reweighting a simulation using the converged bias from the 10 Walker Uniform VES simulation.

Conclusion
We have proposed a very general approach based on the Variationally Enhanced Sampling method to enhance sampling in the transition region, and exemplified with the cases of a S N 2 reaction and condensation of Lennard-Jones gas how our method can be successfully applied to very diverse scenarios. Furthermore, our approach is capable of locating the transition states on the fly and acquiring an increased sampling of the transition region despite the process being a rare-event. This of course under the condition of dealing with a good quality CV. Taking into account the ability of our method to increase the frequency with which the transition region is visited, we believe that this could be especially useful in studying reactions that take place in condensed phases where the effect of the environment on the transition state is relevant and yet difficult to describe with a single conformation.
Nucleation is one example of such class of reactions. We believe our approach in these cases may obtain a much faster F (s) exploration.