{'fraction_non_alphanumeric': 0.03365925413169508, 'fraction_numerical': 0.011306279022814456, 'mean_word_length': 4.450715296337054, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Activated sampling in complex materials at finite temperature: the properly-obeying-probability activation-relaxation technique

While the dynamics of many complex systems is dominated by activated events, there are very few simulation methods that take advantage of this fact. Most of these procedures are restricted to relatively simple systems or, as with the activation-relaxation technique (ART), sample the conformation space efficiently at the cost of a correct thermodynamical description. We present here an extension of ART, the properly-obeying-probability ART (POP-ART), that obeys detailed balance and samples correctly the thermodynamic ensemble. Testing POP-ART on two model systems, a vacancy and an interstitial in crystalline silicon, we show that this method recovers the proper thermodynamical weights associated with the various accessible states and is significantly faster than MD in the diffusion of a vacancy below 700 K.


I. INTRODUCTION
Throughout physics, chemistry and biology, a large proportion of atomistic processes take place on time scales many orders of magnitude longer than the typical period of atomic vibrations. These processes are out of reach of straightforward molecular dynamics (MD), which cannot exceed simulation times equivalent to the microseconds at best. It is not surprizing, therefore, that considerable effort has been devoted in the last few years, with some degree of success, to develop algorithms that overcome this limitation [1,2,3,4,5,6,7,8,9,10].
These methods can be separated into two classes. Activated methods, such as the activation-relaxation technique (ART) [1,9] and related approaches [2,7,8] sample the energy landscape of complex systems by identifying minima connected by minimum energy pathways. This family of methods is very efficient for sampling conformations. Recently, ART was found to be the most efficient method for high-dimensional problems [11]. However, because of an unknown bias in the selection of events for these methods, it is not possible to ensure a proper thermodynamic sampling. While this is not a major limitation for sampling states or even identifying pathways, such as protein folding trajectories, for example, it is sufficiently severe to prevent the use of these methods to sample equilibrium or dynamical quantities.
The second class of methods is based on molecular dynamics and corresponds to methods such as hyper-MD [3], temperature-assisted dynamics (TAD) [4], selfguided dynamics [6] and biased thermodynamics [10]. Until now, the application of these methods has been mostly restricted to simple systems, with a limited num- * Permanent address ber of relatively well-characterized barriers or with a welldefined reaction coordinate. Recently, Choudhary and Clancy have proposed modifications of hyper-MD that could allow its application to disordered materials [12]. It appears, however, that applying a significant boost in hyper-MD could break the thermodynamical character of the algorithm, placing this method in the first category of activated methods.
In this paper, we present an algorithm that offers correct thermodynamical sampling without suffering from the usual exponential slowing down with decreasing temperature. The properly-obeying-probability ART (POP-ART) samples the thermodynamically relevant parts of phase space, hopping efficiently over high barriers separating low-energy basins. We apply this algorithm to two test cases, the diffusion of a self-interstitial and the self-diffusion of a vacancy in Stillinger-Weber silicon [13], to verify the correctness of the method, and to assess its efficiency.
The paper is organized as follows. We start with a brief discussion of limitations of standard activated methods, such as ART. We then introduce POP-ART and show how it can overcome these limitations and ensure a proper thermodynamical sampling. The justification for the various steps needed to construct activated pathways with detailed balance is then discussed in detail and a summary of the algorithm is given. The algorithm is tested in a study of the diffusion of an interstitial and a vacancy in Stillinger-Weber silicon. In Appendix A, we discuss a physical interpretation of the Jacobian used in POP-ART; in Appendix B, we present an analytical calculation for a simple model potential that provides further insights into the method;

II. SAMPLING THE ENERGY LANDSCAPE USING ACTIVATED METHODS
The energy landscape of a system of M atoms can be represented as an N ≡ 3M -dimensional hypersurface, with the "height" indicating the value of the potential of the configuration at a given set of atomic coordinates. In a dynamical regime dominated by rare events, a system spends most of its time oscillating near a local energy minimum, hopping over an energy barrier only when a thermal fluctuation transfers sufficient energy onto the corresponding mode. Since the probability of energy transfer decreases exponentially with its size, the activated trajectory will tend to cross near the lowest-energy point on the ridge, corresponding to a first-order saddle point.
It is possible to reconstruct these trajectories, as a sequence of local minima separated by transition points, using the activation-relaxation technique [1] or related methods [2,7,8]. In its latest form, called ART nouveau [9], this method works in three steps: (1) Starting from a local minimum, a deformation is applied until the lowest curvature of the Hessian matrix, given by becomes negative. (2) The configuration is then pushed along the corresponding direction while the energy in the perpendicular hyperplane is minimized until it reaches a first-order saddle point. (3) The configuration is then pushed slightly further, away from this saddle point, and its energy is minimized using a standard minimization technique.
ART and similar methods have been applied with success to study the topology of the energy landscape and activated mechanisms in a wide range of materials including amorphous and crystalline semiconductors [14,15,16], glassy materials [17], atomic clusters [9,18], and proteins [19,20,21,22,23]. A recent study has shown that ART compares favorably with other related techniques in terms of efficiency and completeness of finding activated mechanisms [11].
It is tempting to use ART approaches to study dynamical trajectories, generating an on-the-fly catalog of events that can be selected using kinetic Monte Carlo. This approach was used by Henkelman and Jónsson in the simulation of growth on an Al (100) surface [24], and by Middleton and Wales [25] for a binary Lennard Jones glass. Comparing with molecular dynamics, Middleton and Wales observed that the kMC results with an on-the-fly catalog are qualitatively incorrect. It is likely that this discrepancy is at least partly due to an uncontrolled selection bias for these types of methods, making it impossible to ensure a proper sampling of the barriers generating the catalog of kMC events.
These limitations are fundamental and cannot be lifted by simply doubling or tripling the sampling. It is essential to incorporate thermodynamics at the core of the algorithm. One method for achieving this is discussed in the next section.

III. THE POP-ART APPROACH
We start by separating the configuration space into two regions: the basin and the saddle regions. Basins are defined in a way that ensures that they contain most of the thermodynamical weight at a given temperature, while the saddle regions are visited only rarely and in passing. The dividing (hyper-)surface between these two regions is chosen to be the surface where the lowest curvature of the potential energy surface equals a threshold value λ 0 . The basins represent the parts of the configuration space where the lowest curvature has a value above λ 0 ; they form a series of disconnected regions surrounding local minima. The saddle region is on the other side of the threshold and includes most other local extrema, such as first-and higher-order saddle points (see Fig. 1). The use of this criterion for separating the configuration space is convenient as the status of any point in the configuration space can be decided locally, without relaxing to a nearby stationary point. For a given threshold, it is always possible that the negative eigenvalue associated with a particular saddle is higher (lower by absolute value) than the chosen threshold, such that it belongs to a basin. As will be seen from the treatment of the method below, this does not invalidate the algorithm but may even be used to one's advantage.
Having separated the energy landscape, we define motion in each of these regions. All the motion within the basin is performed with conventional MD at the desired temperature. Once the configuration hits the dividing surface, the MD is halted, the configuration is brought through the saddle region to a new basin at the same energy, according to the activation rules described below, and the MD is resumed at the same temperature. All steps respect detailed balance and the overall trajectory samples the basins according to the proper thermodynamical ensemble.
The activated part of the algorithm is composed of two steps: (1) the activation trajectory is first generated, from one basin to the other, and then (2) the free energy difference between the beginning and the end of this trajectory is calculated. The latter information is used for the accept-reject step.
In the next subsections, we discuss these two steps in detail before presenting the algorithm as it is currently implemented.

A. The activated trajectory
As in ART, an activation trajectory is created by moving along the eigenvector corresponding to the lowest eigenvalue of the Hessian. Unlike ART, however, there is no relaxation in the perpendicular hyperplane. Instead, all atoms are moved in such a way as to keep the total potential energy constant. This is easily achieved as the configuration is thermalized, with roughly k B T /2 of available potential energy per degree of freedom above that in the local energy minimum. More specifically, the activated trajectory is generated by iterating the following equation: where, h i is the normalized eigenvector at x i corresponding to the lowest Hessian eigenvalue, F i is the total Ndimensional force at x i , ∆τ is a constant that determines the size of the increment, and c is a multiplicative constant, chosen to project the trajectory onto the hypersurface of constant potential energy. The orientation of h 0 is chosen initially so as to point in the direction of more negative curvature, i.e., away from the initial basin; it is updated at each step by requiring that the inner product of the local eigenvector h i with that at the previous step, h i−1 , be always positive. Values of h and F at point x i+1 are obtained iteratively, i.e., initially h i+1 = h i and F i+1 = F i are used in Eq. (2) to obtain a value of x i+1 , then values of h and F are calculated at the new point and inserted into Eq. (2) to get the next iteration, etc.
Unlike in ART, there is no separate relaxation stage, and Eq. (2) is iterated across the saddle region until the new basin is reached, i.e., until the lowest eigenvalue passes the threshold (from below, this time). At this point, the activation-relaxation phase is stopped and MD is resumed starting with the new configuration (see Fig. 1). This ensures that the path generated from x 0 to x p is fully reversible: a configuration in basin p reaching x p would trigger the activation, bringing it to the other end of this path, at x 0 . Reversibility, in a weak sense, is ensured by the symmetric criterion for entering and leaving the saddle region as well as by keeping the path on a hypersurface of constant energy: for each transition, its inverse is also possible. In addition to reversibility, we have to ensure that the relative probabilities of these transitions are correctly weighted; this is discussed below.
The conservation of energy requires that the length of the velocity vector as MD is restarted is equal to that at the beginning of the activation (| v p | = | v 0 |); the direction should be chosen so as to point inside the new basin, but otherwise is arbitrary. After entering the new basin, MD is continued for a very small number of steps -10 or so -to prevent the system from a quick recrossing back to the original basin. This is implemented by letting the system bounce back against the constant eigenvalue surface. After these few steps, the MD stage continues until the system crosses the basin boundary, and another activation is begun, etc.
We note that the activation path does not always lead to a new basin. In our simulations, as is illustrated in Fig. 1, it is not rare to see the trajectory form a circular path, coming back very close to the initial point x 0 , after a long excursion in the saddle region. It is also possible that the trajectory returns to the same basin but not at the starting point. This does not invalidate the algorithm but makes it less efficient.

B. Calculating the event free energy
Once we have a trajectory, it is necessary to compute the free energy difference between its beginning and its end.
Consider the diagram in Fig. 2. This shows schematically a few nearby activation paths between different basin boundaries (which are (N − 1)-dimensional hypersurfaces). Points within area dS 1 in the figure move to points within area dS 2 . If an ensemble with the microcanonical distribution is considered, the density of flux of the trajectories through a hypersurface is the same for any hypersurface at all points having the same potential energy. Then the ratio of the rate of the direct transition (from dS 1 to dS 2 ) to the rate of the inverse transition (from dS 2 to dS 1 ) is equal to the ratio of the areas, dS 1 /dS 2 . If we want the microcanonical distribution to be preserved, this ratio should equal 1. In general, however, it is not unity and we need to add an additional acceptance/rejection step, for instance, accepting a particular activation transition with a Metropolis-like probability P acc = min(dS 2 /dS 1 , 1) for the transition from dS 1 to dS 2 in the figure.
The activation transition can be considered as a transformation between points on different basin boundaries (for instance, point x 0 is transformed into point x p , area dS 1 in Fig. 2 is transformed into area dS 2 ). The ratio J = dS 2 /dS 1 corresponds therefore to the Jacobian of this transformation and it is all we need to ensure detailed balance.

The Jacobian of the activation transformation -the boundary factor
Imagine again a tube formed by nearby activation trajectories (Fig. 2). If the trajectories forming the tube start within area dS 1 on the first basin boundary, then the cross-section at the beginning of the tube is dS ′ 1 = dS 1 cos α 1 , where α 1 is the angle between the normal to the initial basin boundary and the activation trajectory at its start. Likewise, the cross-section at the end of the tube is dS ′ 2 = dS 2 cos α 2 , where α 2 is the angle between the normal to the final basin boundary and the activation trajectory at its end. Including all contributions, the Jacobian can then be written as where We call J b the boundary factor and J xs the cross-section factor.
We start by calculating the boundary part J b . First, we note that Eq. (2) for the activation trajectory is the discretized version of This allows us to get an estimate of c, the factor preserving the total energy during the activation. We write the change in potential energy as where F = ( F · h). Since we want to keep U constant, dU dτ = 0 and F 2 is a sum over N components of the force and thus scales as O(N ), the system size, since all modes are roughly equally excited. For its part, F represents just one component along the activated direction and does not grow with system size. These observations imply therefore that c scales as O(1/N ). Next, we will show that the eigenvector h is nearly tangent to the trajectory. Using Eq. (6), the angle β between the activation trajectory and the eigenvector h is given by (9) Since, as discussed above, F 2 ≪ F 2 for big systems, cF , c 2 F 2 ⊥ ≪ 1, and cos β is nearly 1. Thus we can replace the direction of the trajectory with the direction of h when calculating angles α 1 and α 2 .
Note that since the basin boundary is by definition the constant-eigenvalue surface, the normal to it is parallel to ∇λ, where λ is the lowest eigenvalue. Then where all quantities are evaluated at the beginning of the activation trajectory for α 1 and at its end for α 2 . Thus in order to calculate α 1 and α 2 , we need a way to find ∇λ numerically. The most efficient method is as follows.
By definition, at point x in the configuration space, whereĤ( x) is the Hessian operator at point x. Considering x = (x 1 , x 2 , . . . , x N ) as a set of N parameters {x i }, we can apply the Hellmann-Feynman theorem and find, in vector form and with Einstein's summation convention: where e i are unit vectors along the coordinate axes.
This expression is simply the second derivative of F along the direction of the eigenvector h with the negative sign, i.e., It can be used directly for numerical evaluation. We simply need to compute the force F at three nearby points along the direction of h for each boundary in order to obtain an accurate evaluation of the boundary factor J b .

Analysis of the cross-section Jacobian
The second factor in the total Jacobian is the crosssection factor J xs . To evaluate it, we need to see how the cross-section of an infinitesimally narrow tube formed by activation trajectories changes between the two basin boundaries. Describing the evolution in the configuration space as a function of τ by the equation then, as τ is incremented by dτ , point x transforms into x + f ( x)dτ . The Jacobian of that transformation is given by the determinant of the matrix For an infinitesimal volume δV ( x) around point x, the rate of change simply becomes Which can be solved formally: Going back to the continuous version of our evolution equation, Eq. (6), we note that |c F | = F /| F | ≪ 1. Therefore, the speed along the activation trajectory is nearly constant (equal to one). Thus the infinitesimal volume δV will only change its size and shape in the transverse directions, but will not shrink or expand in the longitudinal direction. Then the tube cross-section ratio between any two points on the trajectory is the same as the volume ratio between the same two points. The logarithm of the cross-section contribution to the Jacobian is then Equation (17), together with Eq. (13) for the boundary factor, can be used in principle to calculate the activation Jacobian. However, straightforward evaluation of the divergences entering Eq. (17) by calculating numerically the derivatives of h and F along N orthogonal directions for many points on the trajectory is computationally very costly and any usable method will require further careful analysis and making certain reasonable approximations, as discussed below. The logarithm of the cross-section factor in the Jacobian is an integral along the activation trajectory: where Compare now the second and the third terms in Eq. (19) to show that the third term can be neglected. In the second term, div F is the trace of the Hessian H taken with the negative sign and is therefore O(N ). Since c is O(1/N ), the second term in Eq. (19) is O(1). Now consider the third term. Using Eq. (8), If we use the coordinate system in which axes are parallel to the eigenvectors of the Hessian at point x (in particular, the zeroth axis is parallel to h), then ∂F i /∂x j = −λ i δ ij , where λ i is the ith eigenvalue of the Hessian. Then the first term in Eq. (20) is and is thus negligible compared to the second term of Eq. (19). In the second term of Eq. (20), In the last expression, the first term is clearly O(1/N ) and thus negligible; the second term is also negligible (this will be so even under a completely unrealistic assumption that all of ∂h i /∂x j are O(1), since F i are of random signs). Thus the third term in Eq. (19) can always be neglected for big enough N and we end up with In Appendix A, we will discuss the physical meaning of the second term in Eq. (23), using the harmonic approximation.

C. Implementation of the POP-ART algorithm
The actual implementation of POP-ART, as used to obtain the results presented in the next section, incorporates the following steps: 1. We start with MD at finite temperature and first equilibrate by rescaling the velocities. We use a 1 fs step and compute the lowest eigenvalue every 10 steps. After we have crossed the basin boundary defined by the threshold, we retrace our MD path and identify the crossing time with an accuracy of 1 fs.
2. We then apply Eq. (2) and generate the event from one basin to another, stopping at the threshold and saving configurations along the way. We take ∆τ = 0.01Å.
3. From the first and final configurations of the activation path, we compute the boundary factor, J b , using Eqs. (10) and (13). 4. We then evaluate the cross-section Jacobian J xs by integrating j, as defined by Eq. (23), over the pathway. We now have the full free energy difference between the entry and exit points.
Straightforward evaluation of the divergences entering Eq. (23) by calculating numerically the derivatives of h and F along N orthogonal directions for many points on the trajectory is computationally very demanding. It is possible, however, to lower this cost significantly while keeping a reasonable accuracy. First, we use a 15iteration Lanczos scheme which allows us to obtain the eigenvector h within O(N ). The divergence of the eigenvector can also be obtained with O(N ), admittedly with a much larger prefactor, provided the potential is sufficiently short-ranged. For two atomic coordinates i and j belonging to atoms which are well outside of each other's interaction range, ∂h i /∂x j ≈ 0. We can exploit this property, to obtain div h = N −1 i=0 ∂h i /∂x i with less than O(N ) force evaluations. For example, two terms in this summation can be obtained with one eigenvector computation: ∂h i /∂x i +∂h j /∂x j ≈ m={i,j} (h m ( x + ∆ e i + ∆ e j ) − h m ( x))/∆. This trick can easily be extended as long as the added coordinates are sufficiently far apart. For the two systems studied here, the unit cell is divided into 25 groups of 40 non-interacting atoms each. The total cost of evaluating div h adds up to 75 Lanczos recursions with only 5 iterations each. 5. The previous two steps provide J b and J xs and thus the free energy difference between the first and the last state on the activation trajectory, which is then used in a Metropolis accept/reject move. If the event is rejected, the component of the velocity normal to the basin boundary is reversed, and MD is continued. If the event is accepted, we continue MD in the new basin, using the initial velocity (reversing the component of the velocity normal to the basin boundary, if necessary). As mentioned above, we run 10 steps to bring the configuration away from the border.

IV. SIMULATION RESULTS
To verify the thermodynamical correctness of the POP-ART method and to investigate its computational efficiency, we have studied the diffusion by interstitials and vacancies in a silicon crystal, described by the Stillinger-Weber potential [13] in systems of respectively 1001 and 999 atoms, and a cubic simulation cell of 27.136Å 3 .

A. Interstitial diffusion
First, we look at the interstitial in a 1001-atom cell of Si. In principle, many interstitial sites are possible in this system, but only three of them are significantly populated: the lowest-energy configuration and two configurations with almost the same energy (within 0.01 eV), about 0.75 eV above the first one [26]. We are interested in computing the probability of being in one of the higher-energy states.
If the population of the higher-energy states were determined by the energy difference alone, this would amount to a population of the higher-energy states of only 0.07% at 1200 K. However, there are degeneracies and the potential wells of the higher-energy states are much flatter than that for the low-energy state, leading to a noticeable entropy difference. Because of this difficulty, we extract the thermodynamical equilibrium between these two states with MD. This forces us to perform the tests at a relatively high temperature. Here, we report results for 1000 K and 1200 K.
These conditions are not ideal for POP-ART since at such high temperatures the jumps between the minima are rather frequent and straightforward MD is quite efficient. But given the significance of both energetic and entropic contributions, as well as some anharmonic effects present at such high temperatures, it is a very good test for the accuracy (rather than efficiency) of POP-ART.
At 1000 K and 1200 K, the system can spend a nonnegligible amount of time outside the basins, i.e., in the saddle region (which would not be the case in systems more appropriate for POP-ART), leading to a different value of the ratio of time spent in the upper vs. the lowerenergy states. We therefore need to distinguish carefully between the probability of being in the attraction region of a given minimum and the time spent in a particular basin (understood as defined in this paper) as a fraction of the total time spent within all basins. It is the latter quantity that we use to compare to the POP-ART result.
At 1000 and 1200 K, the MD result for this quantity is, respectively, 1.6 ± 0.1% and 3.6 ± 0.1%, determined as an average over 25 runs, each of which lasted 10 ns and an assignment to a basin, with threshold λ = −2 and −1 eV/Å 2 , done every 100 fs.
For POP-ART the fraction is obtained as an average over 5 runs, of 10000 iterations of the POP-ART algorithm each. The assignment to a basin is done every 100 fs, excluding the 10 steps for moving away from the border after each generated event. At 1000 and 1200 K, we obtain the respective ratios of 1.3 ± 0.3% and 3.5 ± 0.3%. Clearly POP-ART samples accurately the thermodynamical weight of the high and low-energy local mimina.

B. Vacancy diffusion
Having established the accuracy of POP-ART, we now characterize its efficiency. For this, we consider a vacancy in a 999-atom cell of Si. Vacancy diffusion is associated with a single activation barrier of 0.43 eV [27]. Assessing the efficiency of POP-ART relative to MD can be cleanly done in this system, since the speed of phase space exploration is simply given by the vacancy's diffusion coefficient, which follows the Arrhenius law. The comparison with MD is done on the basis of the number of calls to the force routine, since that takes ∼99% of the computer time. Figure 3 shows an Arrhenius plot of the diffusion rate per million force operations obtained by MD and by POP-ART, as a function of temperature. The diffusion rate is defined as the vacancy hopping rate, which is linearly related to the square displacement per unit time. The value of λ used in the POP-ART simulation is selected to optimize the diffusion rate. The threshold values used in the simulation are λ = −2 eV/Å 2 for T > 750 K, λ = −1 eV/Å 2 for 600 K < T < 750 K and λ = 0 eV/Å 2 for T < 600 K. Clearly, the diffusion rate per force evaluation obtained with POP-ART does not show activated behavior and provides a significant boost at temperatures below 700 K, reaching a factor of more than 4 orders of magnitude at room temperature. Interestingly, the diffusion rate per force evaluation is not constant with POP-ART but can be fitted by a 1/T curve. The slowing down with decreasing temperature is related to the fact that before POP-ART attempts to find its way to a new basin, it needs to reach a given curvature threshold. At a very low temperature, even this will be a very rare event.

V. DISCUSSION AND CONCLUSION
Efficient sampling of slow systems is one of the main challenges in computational physics today. For disordered systems such as glasses, for example, standard thermodynamical methods fail because of the very small phase space occupied by the relevant configurations. Activated methods, such as ART, overcome these difficulties by generating physically-possible trajectories through the conformation space but they do not offer a proper thermodynamical sampling.
The properly-obeying-probability activationrelaxation technique (POP-ART) lifts these limitations by generating activated trajectories with proper thermodynamical weighting. Mixing molecular dynamics with activation over barriers, this algorithm respects detailed balance and samples in a well-defined thermodynamical ensemble.
To verify the correctness of POP-ART, we sampled the states visited by an interstitial in c-Si. In this system, the higher-energy states of the interstitial are energetically suppressed but entropically favored. Comparing with MD, we found that POP-ART samples the highenergy states with the proper probability, demonstrating its correctness. In order to assess the efficiency of this method, we looked at the vacancy diffusion also in c-Si. In this case, POP-ART is found to outperform MD below 700 K, and it is about four orders of magnitude faster at room temperature.
One of the advantages of POP-ART is that all the information it needs is local. This makes it possible to apply a number of approximations to increase further its efficiency without sacrificing the sampling. POP-ART can also be extended to reproduce the correct activated dynamics; this is currently examined and will be reported in a further publication.
Using the harmonic approximation, we can reveal better the physical meaning of the second term in Eq. (23). First of all, using the expression for c, Eq. (8), we get At relatively low temperatures, the system is well described by the harmonic approximation. The potential around the minimum can then be rewritten as where x i represent the normal modes, and the force becomes Using the same approximation for F 2 and neglecting anharmonicity, we get: Since the spectrum should be dense, x 2 i can be replaced with their thermal averages x 2 i = k B T /k i , giving and After integrating over the whole activation trajectory, we obtain where Then the contribution of the c div F term to the Jacobian J is a Boltzmann factor. ∆E is essentially the energy change that would have occurred along the activation trajectory, if it were parallel to h everywhere and the energycorrecting c F term in Eq. (2) was not there. The transition probability between two minima should contain both energetic and entropic contributions. Given the above result, it is tempting to associate the div h term with the entropic and the c div F with the energetic contribution; however, an example considered in Appendix B shows that the reality is more complex: in fact, ∆E , defined as above, is temperature-dependent and the c div F term therefore contains both energetic and the entropic contributions.
Appendix B: the Jacobian in a model potential To get some insight into the physical meaning of the activation Jacobian and its different components, consider the following model, defined by the potential: Here U 0 (x 0 ) is a function with two minima and a maximum between them, so that coordinate x 0 describes the activated mode (i.e., is the "reaction coordinate"), and the other N − 1 degrees of freedom serve as the "heat bath". "Force constants" k i (x 0 ) are assumed to remain positive for all x 0 of interest (e.g., between the minima). Our model does not represent the most general situation. In particular, the eigenmode-following transition path (such as, e.g., ART would find) is a straight line (coinciding with the 0th axis); also, on that line, an eigenvector for a particular mode has the same direction (parallel to a coordinate axis) everywhere. However, the model is interesting enough: as the frequencies of the bath modes (determined by k i (x 0 )) can change along the activation path, there are both energetic and entropic contributions to the probability of being in a particular place along the reaction coordinate. Indeed, the probability density of having the zeroth coordinate equal to x 0 is or where the free energy and the entropy We will assume in what follows that k i (x 0 ) are linear functions, i.e., The matrix elements of the Hessian for the potential given by Eq. (33) are and the remaining elements are zero. The force components are For a Hessian with only H ii and H 0i = H i0 non-zero, where C is the normalization constant, and the eigenvalue λ is the lowest solution of the following equation: In a real physical system, the entropy change along an activation trajectory will always remain finite and of the same order of magnitude as the thermal energy, as the system size increases. In our model, this will mean that most k (1) i are small enough. Note that this is essentially the same as the assumption of most modes being nearly harmonic that we have used when approximating the c div F term in the Jacobian. Specifically, For simplicity, we will assume in addition that This will be the case, in particular, at low enough T , when the magnitudes of most x i are small. Under these conditions, The constant-eigenvalue surfaces are then orthogonal to the 0th axis and the cosine of the angle between the normal to a constant-eigenvalue surface and h is nearly 1, so the boundary contribution to the activation Jacobian can be neglected. Consider now the cross-section factor. We need to calculate j, as given by Eq. (23). First of all, , . . .
and so div h = Further, using Eq. (40), We can now use Eq. (29) to calculate c div F ; it is, however, instructive to repeat the considerations, to see how exactly the appropriate approximations are made in this particular case. We obtain (where an approximation similar to Eq. (26) is obtained by neglecting the term associated with the anharmonic activated mode) and In the last expression, bearing in mind the condition (44), the second term in the parentheses is of the same order of magnitude as the first one, and then both of them are negligible compared to the second sum (this is an approximation similar to the one used in Eq. (27)). We therefore obtain As explained in Appendix A, we can replace x 2 i with their thermal averages and then c div F ≈ − U ′ 0 (x 0 ) (55) Then finally and the Jacobian J is where ∆U 0 is the change in U 0 , ∆S is the change in entropy S [as given by Eq. (37)] between the two basin boundaries, and ∆F = ∆U 0 − T ∆S -exactly as expected.