Chapter 3
Properties of Rearrangement Pathways

  2.5.6 A Revised Connection Algorithm
  2.5.7 Applications to Isomerisation of LJ38 and LJ75
  2.5.8 A Dijkstra-based Selector
  2.5.9 Applications to Tryptophan Zippers
 2.6 Summary

There is a difference between knowing
the path and walking the path.
Andy and Larry Wachowski, The Matrix

3.1 Introduction

The number of elementary rearrangements increases exponentially with system size as for the number of transition states. For instance, there are approximately 30,000 such pathways on the PES of the 13-atom cluster bound by the Lennard-Jones potential. When permutation-inversion isomers are included, this number increases by a factor of order 2 × N! [9]. For PES’s that support such a large number of stationary points a whole range of properties is spanned by the corresponding rearrangement pathways. Understanding these properties can be helpful in answering questions such as why some rearrangement pathways are harder to find than others and whether there is any correlation between these properties that we could potentially utilise in studying PES’s.

Two activation barriers can be defined for each pathway in terms of the energy difference between the transition state and each of the minima. For non-degenerate rearrangements [9, 148] the two sides of the path are termed uphill and downhill, where the uphill barrier is the larger one, which leads to the higher minimum. The barriers and the normal modes of the minima and transition states can be used to calculate rate constants using harmonic transition state theory [45, 46, 48]. More sophisticated treatments based on anharmonic densities of states are possible but can be hard to reconcile with the coarse-grained view of the landscape adopted here, as a detailed knowledge of the basins of attraction is required [168, 169].

For each local minimum a catchment basin can be defined in terms of all the configurations from which steepest-descent paths lead to that minimum [170]. Some of these paths originate from transition states on the boundary of the catchment basin, which connect a given minimum to adjacent minima. The integrated path length for such rearrangements provides a measure of the separation between local minima, and may be related to the density of stationary points in configuration space. The integrated path length is usually approximated as the sum of Euclidean distances between configurations sampled along appropriate steepest-descent paths [9]. It provides a convenient coordinate for monitoring the progress of the reaction.

Calculated pathways can always be further classified mechanistically. For example, some rearrangements preserve the nearest-neighbour coordination shell for all the atoms. In previous studies of bulk models these cage-preserving pathways were generally found to outnumber the more localised cage-breaking processes, which are necessary for atomic transport [171]. It was found that the barriers for cage-breaking and cage-preserving processes were similar for bulk LJ systems, while the cage-breaking mechanisms have significantly higher barriers for bulk silicon modelled by the Stillinger-Weber potential [171].

For minima separated by increasing distances in configuration space, the pathways that connect them are likely to involve more and more elementary steps, and are not unique. Finding such paths in high-dimensional systems can become a challenging task [7, 8]. Some difficulties have been attributed to instabilities and inefficiencies in transition state searching algorithms [7, 107], as well as the existence of very different barrier height and path length scales [7]. A new algorithm for locating multi-step pathways in such cases was presented in the previous chapter.

In the present work we have used the doubly nudged elastic band (DNEB) method [7] in conjunction with eigenvector-following (EF) algorithms [1224] to locate rearrangement pathways in various systems. The LJ potential was used to describe the 13- and 75-atom Lennard-Jones clusters, LJ13 and LJ75. We have also considered a binary LJ (BLJ) system with parameters σAA = 1.0, σBB = 0.88, σAB = 0.8, 𝜖AA = 1.0, 𝜖BB = 0.5, 𝜖AB = 1.5, where A and B are atom types. The mixture with A : B ratio 80 : 20 provides a popular model bulk glass-former [171, 172]. We employed a periodically repeated cubic cell containing 205 A atoms and 51 B atoms. The density was fixed at 1.2σAA3 and the Stoddard-Ford scheme was used to prevent discontinuities [173].

The motivation for this work was our observation that construction of some multi-step pathways using the connection algorithm described in Chapter 2 is particularly difficult. We were unable to relate these difficulties to simple properties of the pathways such as the integrated path length, the uphill and downhill barriers, or the barrier and path length asymmetries. Instead, more precise measures of localisation and cooperativity are required, as shown in the following sections. It also seems likely that such tools may prove useful in analysing the dynamics of supercooled liquids, where processes such as intrabasin oscillations and interbasin hopping have been associated with different rearrangement mechanisms [174]. In particular, cooperativity is believed to play an important role at low temperatures in glass-forming systems [175], and dynamic heterogeneity may result in decoupling between structural relaxation and transport properties for supercooled liquids [176].

3.2 Localisation

The outcome of a pathway calculation for an atomic system will generally be a set of intermediate geometries, and the corresponding energies, for points along the two unique steepest-descent paths that link a transition state to two local minima. This discrete representation is a convenient starting point for our analysis of localisation and cooperativity. We number the structures along the path j = 1,2,,Nf starting from one of the two minima and reversing the other steepest-descent path, so that structure Nf corresponds to the other minimum. The transition state then lies somewhere between frames 1 and Nf. We define the three-dimensional vector ri(j)  to contain the Cartesian coordinates of atom i for structure j, i.e. ri(j) =(Xi(j),Y i(j),Zi(j)), where Xi(j) is the X coordinate of atom i in structure j, etc.

For each atom i we also define the displacement between structures j 1 and j as 

di(j) =|ri(j) ri(j 1)|. (3.1)

Hence the sum of displacements 

di = j=2Nf di(j), (3.2)

is an approximation to the integrated path length for atom i, which becomes increasingly accurate for smaller step sizes. The total integrated path length, s, is approximated as

s = j=2Nf i=1Ndi j2, (3.3)

where N is the total number of atoms. s is a characteristic property of the complete path, and is expected to correlate with parameters such as the curvature and barrier height for short paths [9, 177, 178].

The set {d1,d2,,dN} containing all N values of di will be denoted {d}, and analogous notation will be used for other sets below. We will also refer to the frequency distribution function, which can provide an alternative representation of such data [179]. For example, the frequency distribution function for a given continuous variable, x, tells us that x occurs in a certain interval (x) times.

Our objective in the present analysis is to provide a more detailed description of the degree oflocalisationandcooperativity corresponding to a given pathway. The first index we consider is Np, which is designed to provide an estimate of how many atoms participate in the rearrangement. We will refer to a rearrangement as localised if a small fraction of the atoms participate in the rearrangement, and as delocalised in the opposite limit. The second index we define, Nc, is intended to characterise the number of atoms that move simultaneously, i.e. cooperatively. We will refer to a rearrangement as cooperative if most of the atoms that participate in the rearrangement move simultaneously, and as uncooperative otherwise.

The nth moment about the mean for a data set {x1,x2,...,xM} is the expectation value of (xi x)n, where x = i=1MxiM and M is the number of elements in the set. Hence for the set {d} defined above we define the moments, mn, as

mn = 1 N i=1N d i dn. (3.4)

The kurtosis of the set {d} is then defined as the dimensionless ratio 

γ(d) = m4 (m2)2, (3.5)

and provides a measure of the shape of the frequency distribution function corresponding to {d}. If only one of the atoms moves, or all atoms except one move by the same amount, then γ(d) = N. Alternatively, if half the atoms move by the same amount whilst the others are stationary, then γ(d) = 1. Hence, a distribution with a broad peak and rapidly decaying tails will have a small kurtosis, γ 𝒪(1), while a distribution with a sharp peak and slowly decaying tails will have a larger value (Figure 3.1). The kurtosis can therefore identify distributions that contain large deviations from the average value [179]. For comparison, a Gaussian distribution has γ = 3 and a uniform distribution has γ = 1.8.

Figure 3.1:Two frequency distribution functions 1 and 2 of a continuous variable x are contrasted. Both functions have the same average m1 and standard deviation m2. However, due to the long tails, 1 has a significantly larger fourth moment m4 and, hence, a larger kurtosis, γ. If γ > 3 the distribution is said to be peaked or leptocurtic. Distributions with γ > 3 are known as platycurtic or heavy-tailed.

The above results show that the kurtosis γ(d) can be used to quantify the degree of localisation or delocalisation of a given rearrangement. However, it has the serious disadvantage that highly localised and delocalised mechanisms both have large values of γ(d). Since we are interested in estimating the number of atoms that move relative to the number with small or zero displacements, a better approach is to use moments taken about the origin, rather than about the mean, i.e.

mn = 1 N i=1Nd in, (3.6)

following Stillinger and Weber [152]. Note that while m1 = 0, the first moment m1 is the mean value. We therefore estimate the number of atoms that participate in the rearrangement, Np, as

Np = N γ(d), (3.7)


γ(d) = m4 m22. (3.8)

For the system with N atoms, if only one atom moves Np = 1, while if K atoms move by the same amount, Np = K.

A similar index to Np has been employed in previous work [9, 19, 180] using only the displacements between the two local minima, which corresponds to taking Nf = 2 in Equation 3.2. Using di values based upon a sum of displacements that approximates the integrated path length for atom i, rather than the overall displacement between the two minima, better reflects the character of the rearrangement, as it can account for the nonlinearity of the pathway. To describe this property more precisely we introduce a pathway nonlinearity index defined by 

α = s D s , (3.9)

where D is the Euclidean distance between the endpoints,

D = i=1N ri(Nf) ri(1)2. (3.10)

We calculated the α values for a database of 31,342 single transition state pathways of LJ75 (hereafter referred to as the LJ75 database). The average value of α was 0.4 with a standard deviation of 0.2, and, hence, there is a significant loss of information if γ is calculated only from the endpoints using Nf = 2. Comparison of the two indices for the LJ75 database revealed many examples where neglect of intermediate structures produces a misleading impression of the number of atoms that move. The definition in Equation 3.7 is therefore suggested as an improvement on previous indices of localisation [9, 19, 152, 180].

3.3 Cooperativity

Np atoms can participate in a rearrangement according to a continuous range of cooperativity. At one end of the scale there are rearrangements where Np atoms all move simultaneously [see Figure 3.2 (a)]. Although these paths exhibit the highest degree of correlated atomic motion they do not usually pose a problem for double-ended transition state search algorithms [7, 30]. Linear interpolation between the two minima tends to generate initial guesses that lie close to the true pathway, particularly if α 0. At the opposite extreme, atoms can move almost one at a time, following a dominopattern [see Figure 3.2 (b)]. Locating a transition state for such rearrangements may require a better initial guess, since linear interpolation effectively assumes that all the coordinates change at the same rate.

Figure 3.2:Comparison of cooperative (a) and uncooperative (b) rearrangements of the LJ75 cluster, for mechanisms that are localised mainly on two atoms. The displacement d as a function of the frame number j is shown for the two atoms that move the most. Panels (c) and (d) illustrate the potential energy V as a function of frame number j for the rearrangements in (a) and (b), respectively. Dashed circles indicate flatter parts of the energy profile, which correspond to the most cooperative regions of the pathway.

The degree of correlation in the atomic displacements can be quantified  by considering the displacementoverlap

Ok = j=2Nf Ok j = j=2Nf min[dc(1)(j),dc(2)(j),...,dc(k)(j)], (3.11)

where the index k indicates that O was calculated for k atoms numbered c(1), c(2),...,c(k). c is a k-dimensional vector that represents a particular choice of k atoms from N, and hence there are CNk = N!k!(N k)! possible values of Ok. The index O can be thought of as a measure of how the displacements of the atoms c(1), c(2), etc. overlap along the pathway. For example, if two atoms move at different times then O2 is small for this pair because the minimum displacement in Equation 3.11 is always small. However, if both atoms move in the same region of the path then O2 is larger.

We now explain how the statistics of the overlaps, Ok, can be used to extract a measure of cooperativity (Figure 3.3). Suppose that m atoms move simultaneously in a hypothetical rearrangement. Then all the overlaps Ok for k > m will be relatively small, because one or more atoms are included in the calculation whose motion is uncorrelated with the others. For overlaps Ok with k m the set of Ok for all possible choices of k atoms from N will exhibit some large values and some small. The large values occur when all the chosen atoms are members of the set that move cooperatively, while other choices give small values of Ok. Hence the kurtosis of the set {Ok}, γ(Ok), calculated from moments taken about the origin, will be large for k m, and small for k > m.

To obtain a measure of how many atoms move cooperatively we could therefore calculate γ(O2), γ(O3), etc. and look for the value of k where γ(Ok) falls in magnitude. However, to avoid an arbitrary cut-off, it is better to calculate the kurtosis of the set {γ(O2), γ(O3),...,γ(Ok)}, or γ[γ(O)] for short. There are N 2 members of this set, and by analogy with the definition of Np = Nγ(d), we could define a cooperativity index Nc = (N 2)γ[γ(O)] + 1. Then, if γ(O2) is large, and all the other γ(Ok) are small, we obtain γ(γ(O)) N 2 and Nc 2, correctly reflecting the number of atoms that move together.

Figure 3.3:γ(O2) plotted against γ(d) for the LJ75 pathway database. The figure shows how γ(O2) can discriminate between rearrangements that have similar values of γ(d) but different cooperativity. The data point for the most cooperative rearrangement localised on two atoms has Np = Nc = 2.

In practice, there are several problems with the above definition of Nc. Calculating Nc in this way quickly becomes costly as the number of atoms and/or number of frames in the pathway increases, because the number of elements in the set {Ok} varies combinatorially with k. Secondly, as k approaches N the distribution of all the possible values for Ok becomes more and more uniform. Under these circumstances deviations from the mean that are negligible in comparison with the overall displacement can produce large kurtosises. Instead, we suggest a modified (and simpler) definition of Nc, which better satisfies our objectives.

We first define the overlap of atomic displacements in a different manner. It can be seen from Equation 3.11 that the simultaneous displacement of l atoms is included in each set of overlaps {Ok} with k l. For example, if three atoms move cooperatively then both the {O2} and {O3} sets will include large elements corresponding to these contributions. Another redundancy is present within {Ok}, since values in this set are calculated for all possible subsets of k atoms and the displacement of each atom is therefore considered more than once. However, we can avoid this redundancy by defining a single koverlap, rather than dealing with CNk different values.

Recall that di(j) is the displacement of atom i between frames j 1 and j. The ordering of the atoms is arbitrary but remains the same for each frame number j. We now define Δi(j) as the displacement of atom i in frame j, where index i numbers the atoms in frame j in descending order, according to the magnitude of di(j), e.g. atom 1 in frame 2 is now the atom with the maximum displacement between frames 1 and 2, atom 2 has the second largest displacement etc. As the ordering may vary from frame to frame, the atoms labelled i in different frames can now be different. This relabelling greatly simplifies the notation we are about to introduce. Consider the k-overlap defined as 

Θk = 1 Δtot j=2Nf [Δk(j) Δk+1(j)], (3.12)

where k ranges from 1 to N, Δtot = j=2NfΔ1(j) and ΔN+1(j) is defined to be zero for all j. A schematic illustration of this construct is presented in Figure 3.4. For example, if only two atoms move in the course of the rearrangement, and both are displaced by the same amount (which may vary from frame to frame), the only non-zero overlap will be Θ2.

Figure 3.4:The Θ indices for a hypothetical rearrangement localised on three atoms. For each of these atoms the displacement d as a function of frame number j is shown. The di values in successive frames are connected with dotted (d1), dashed (d2) and solid (d3) lines. The corresponding contributions to Θ1, Θ2, and Θ3 are shown as shaded squares and are labelled accordingly. If the remaining N 3 atoms do not participate, and the area of one square is S, the only non-zero overlaps will be Θ1, Θ2, and Θ3 with values 59, 39 and 19, respectively.

We can now define an index to quantify the number of atoms that move cooperatively as

Nc = k=1NkΘ k. (3.13)

If only one atom moves during the rearrangement then Nc = 1, while if K atoms displace cooperatively during the rearrangement then Nc = K. This definition is independent of the total displacement, the integrated path length, and the number of atoms, which makes it possible to compare Nc indices calculated for different systems.

3.4 Applications to LJ13 and LJ75 Clusters and BLJ256 Liquid

Figure 3.2 shows results for the most cooperative and uncooperative processes we have found for the LJ75 cluster that are localised mainly on two atoms. In these calculations we have used the database of transition states that was found previously as the result of a discrete path sampling calculation conducted for this system [8, 10]. The cooperative rearrangement [Figure 3.2(a,c)] is the one with the maximum two-overlap Θ2. For this pathway Θ2 = 0.7, Np = 3.4, and Nc = 7.7. The values of Np and Nc both reflect the fact that the motion of the two atoms is accompanied by a slight distortion of the cluster core. This example shows that while Np and Nc allow us to quantify localisation and cooperativity, and correctly reflect the number of atoms that participate and move cooperatively in ideal cases, there will not generally be a simple correspondence between their values and the number of atoms that move. This complication is due to the fact that small displacements of atoms in the core will generally occur, no matter how localised the rearrangement is. In addition, the data reduction performed in Equation 3.7 and Equation 3.13 means that a range of pathways can give the same value for Np or Nc. Since the size of the contribution from a large number of small displacements depends on the shape of the displacement distribution function the number of possibilities grows with the size of the index.

The uncooperative rearrangement depicted in the Figure 3.2(b,d) was harder to identify. In principle we could have selected the pathway that maximises Θ1 from all the rearrangements localised on two atoms. However, this approach picks out rearrangements localised on one atom, where distortion of the core accounts for the value of Np > 1. Instead, we first selected from all the rearrangements with Np < 4 those where two atoms move by approximately the same amount, while the displacement of any other atom is significantly smaller. These are the rearrangements that maximise 1γ({d1,d2}), where d1 and d2 are the total displacements of the two atoms that move the most. After this procedure we selected the rearrangement with the maximum value of Θ1. Figure 3.2 (b) shows that this rearrangement features the displacement of one atom at a time, and the atom that moves first also moves last. For this pathway the values of Θ1, Np and Nc are 0.7, 3.8 and 5.3, respectively. Further illustrations and movies of the corresponding rearrangements are available online [181].

Figure 3.2 illustrates several general trends that we have observed for cluster rearrangements. Firstly, we have found that the barrier height is smaller for the cooperative rearrangements [Figure 3.2(c,d)]. Usually atoms that move cooperatively are neighbours. Rearrangements generally involve a change of the environment for the atoms that move. Cooperative motion can reduce this perturbation since for any of the participants the local environment is partly preserved because it moves with the atom in question. Flatter points on the energy profile [circled in Figure 3.2(d)] usually signify a change in the mechanism, i.e. one group of atoms stops moving and another group starts. By comparing (b) and (d) in Figure 3.2 we conclude that flatter points on the energy profile correlate with the most cooperative parts of this rearrangement.

A simple correlation between barrier heights and Np and Nc does not seem to exist. The barrier height is not a function of cooperativity alone, but also of the energetics of the participating atoms. The way the Np and Nc indices have been defined can make them insensitive to details of the rearrangements that will affect the energetics. For instance, neither index depends on the location of the participating atoms or the directionality of their motion. In most cases cooperatively moving particles are adjacent, i.e. localised in space; however, long-distance correlations of atomic displacements also occur. One such case is depicted in Figure 3.5 (a). This path is nearly symmetric with respect to the integrated path length (π = 0.01), but is very asymmetric with respect to the uphill and downhill barrier heights (β = 0.91). This rearrangement has Np = Nc = 10. Interestingly, Np and Nc calculated separately for both sides of the pathway are very similar, i.e. the two steepest-descent pathways cannot be distinguished using these indices. Close inspection of this rearrangement reveals that one side of the pathway involves the rearrangement of two atomic triplets that share a vertex, while the other side involves the drift of all five atoms on the surface of the cluster [see the insets in Figure 3.5 (a)]. Although Nc does not distinguish these cases, the motion in the second side of the path is more cooperative. The participating atoms move together, which results in a significantly lower downhill barrier.

Figure 3.5:Two limitations of the cooperativity index Nc. (a) Nc is not sensitive to the spatial positions of the cooperatively moving atoms, nor to the directionality of their motion. The energy profile is depicted for a rearrangement of the LJ75 cluster, which is very asymmetric with respect to barrier height but has similar integrated path lengths on either side of the transition state. The cooperativity index Nc evaluated separately for the two sides is about 10 in both cases. The motion of the five atoms that displace the most is shown schematically relative to a reference atom (black). (b) The displacement of two (left) and three (right) atoms d as a function of the frame number j is shown schematically for a hypothetical pathway. The rearrangement on the left is more cooperative because two atoms move together over a longer region of the path. However, the current definition of Nc does not distinguish between these two cases.

Np and Nc also describe properties of the whole pathway. A significant number of pathways that we observed were rather non-uniform, i.e. very cooperative phases alternated with uncooperative ones. To distinguish such pathways in the LJ75 database we calculated a set {Np} containing Nps evaluated for each pair of adjacent frames. Then a selection of pathways was made with m2(m1)2 < 0.01, where m2(m1)2 is the moment ratio evaluated for the set {Np}. While this procedure ensured that Np corresponds closely to the number of atoms that moves between any two snapshots of the rearrangement, it did not distinguish cases where different atoms contribute to the value of Np in different frames [see Figure 3.5 (b)]. The average uphill and downhill barriers for this subset of rearrangements are 100 times smaller than the average barriers for the complete LJ75 database (Table 3.1). Figure 3.6 shows that Np and Nc calculated for these rearrangements are highly correlated and span a range of values, implying that widely different pathways are represented. Finally, all the selected pathways are an order of magnitude shorter than the average path length for the whole database, even though this database contains many short rearrangements localised on one atom.

Table 3.1:Average uphill and downhill barriers, average integrated path length for the LJ75 rearrangement pathway database. Values are given for the whole database containing 31,342 paths, and for a subset containing the 57 most cooperative paths. The units of energy and distance are 𝜖 and σ, respectively.

All Cooperative

Uphill barrier 3.03 0.06
Downhill barrier0.97 0.03
Path length 3.08 0.58

Figure 3.6:Nc as a function of Np calculated for the 57 most cooperative rearrangements from the LJ75 pathway database.

Figure 3.7 shows the values of the Np and Nc indices plotted against each other for pathway databases calculated for LJ13, LJ75 and BLJ256. The two databases for BLJ256 labelled as 1 and 12 are taken from Reference [171] and correspond to databases BLJ1 and BLJ12 in that paper. BLJ1 and BLJ12 were obtained using two different sampling schemes intended to provide an overview of a wide range of configuration space and a thorough probe of a smaller region, respectively. These databases were constructed by systematic exploration of the PES, and we refer the reader to the original work for further details [171]. Each BLJ256 database contains 10,000 transition states. The LJ13 and LJ75 databases were obtained in discrete path sampling (DPS) studies [8, 10]. and contain 28,756 and 31,342 transition states, respectively. Figure 3.7 is a density plot where darker shading signifies a higher concentration of data points. The outlying points are connected by a solid line to define the area in which all the points lie. Figure 3.7 shows that as Np grows the allowed range of Nc increases, especially for LJ13. For the LJ75 database rearrangements with Nc > N2 appear to be very rare or poorly sampled. Figure 3.7 also shows that for all these systems rearrangements localised on two or three atoms dominate. This result may be an intrinsic property. However, it may also be due to the geometric perturbation scheme used in producing the starting points for the transition state searches employed in generating these databases. For databases BLJ1 and BLJ12 there are significantly more rearrangements with larger values of Np and Nc compared to LJ75, which suggests that the abundance of very localised rearrangements for clusters may be a surface effect. The apparent absence of cooperative rearrangements in LJ75 database for large values of Np may be due to the fact that only the pathways between compact phases of this cluster were sampled thoroughly, i.e. rearrangements involving liquid-like structures that are expected to have larger values of Nc are probably underrepresented.

Figure 3.7:Nc as a function of Np calculated from four pathway databases for LJ13, LJ75 and BLJ256 systems. Due to the large number of data points we employ a density plot, where the darkest shading corresponds to the highest concentration of points. Outlying points are connected to illustrate the boundaries of the data area. The two BLJ databases are taken from Reference [171].

Figure 3.8 depicts the average barrier as a function of the participation and cooperativity indices. Np and Nc were calculated separately for both sides of the pathway and the corresponding barriers (uphill or downhill) were averaged to produce a density plot of barrier height. Data boundaries do not coincide with those shown in Figure 3.7 due to numerical imprecision in preparing Figure 3.7. Figure 3.8 illustrates that for each system cooperative rearrangements have the lowest barriers, irrespective of the value of Np. For clusters, cooperative rearrangements have lower barriers than uncooperative rearrangements with Np as small as 1 3, while for bulk barriers corresponding to rearrangements with low Np become comparable to these for very cooperative rearrangements with high Np.

Figure 3.8:Average barrier height as a function of Np and Nc calculated for the same LJ13, LJ75 and BLJ256 pathway databases used in Figure 3.7. The indices were calculated separately for the two sides of each path. In this case the darkest shading corresponds to the highest barriers. Outlying points are connected to illustrate the boundaries of the data area. The two BLJ databases are taken from Reference [171].

In further computational experiments we found that attempts to connect the endpoints of uncooperative pathways using the algorithm described in Chapter 2 either required more images and iterations or converged to an alternative pathway. In some cases additional difficulties arose, such as convergence to a higher index saddle instead of a transition state, which can happen if the linear interpolation guess conserves a symmetry plane. Figure 3.9 shows Nc calculated from Equation 3.13 plotted against Np for the LJ75 pathway database. Knowing the integrated path length, s, for each pathway we started doubly nudged elastic band calculations with three images per unit of distance and 30 iterations per image. Most of the points that correspond to runs that failed or converged to an alternative pathways are concentrated in the region of small values for NcNp.

Figure 3.9:Nc as a function of Np calculated for the LJ75 pathway database. For each pathway we conducted DNEB calculations [7] assuming prior knowledge of the path. Every DNEB calculation employed 3s images and 90s iterations in each case, where s is the integrated path length. Out of 31,342 DNEB runs 25,158 yielded a connected pathway while the rest did not (FAILED). Connected pathways are classified further as one-step pathways involving the correct transition state (OK) or an alternative transition state (ALT), or as multi-step pathways (MULTI), which involve more than one transition state. For each set of data points best fit straight lines obtained from linear regression are shown and labelled appropriately.

As can be seen from Figure 3.7, the LJ13 database contains significantly more pathways with large values of Np and Nc compared to LJ75, where most of the rearrangements are localised on two atoms. The fact that the LJ13 database is almost exhaustive then suggests that localised rearrangements either start to dominate as the system size increases or that the sampling scheme used for LJ75 was biased towards such mechanisms. Systematic sampling of the configuration space for stationary points often employs perturbations of every degree of freedom followed by minimisation [162]. The LJ13 database was obtained in this fashion, while the LJ75 database was generated during the DPS approach [8]. In this procedure discrete paths are perturbed by replacing local minima with structures obtained after perturbing all the coordinates and minimising. To investigate whether the perturbation scheme can affect the resulting database of stationary points in more detail we consider the case of LJ13, since nearly all the transition states are known. Figure 3.10 presents the results of two independent runs aimed at locating most of the transition states for this system. Every cycle a perturbation was applied to a randomly selected transition state from the database and the resulting geometry was used as a starting point for a new transition state search using eigenvector-following [19, 2224, 72]. Only distinct permutation-inversion isomers were saved. In the first run (bottom curve) every degree of freedom was perturbed by 0.4x, where x is a random number in the interval [1,1] [162]. For the second run (top curve) we introduced a perturbation scheme including correlation. 2 K N2 atoms out of N were displaced by a vector 0.4(x1,x2,x3), where the components x1, x2 and x3 are again random numbers drawn from [1,1]. The K atoms to be displaced were selected based on their relative positions in the cluster. One atom was first selected at random, while the remaining K 1 were chosen to be its closest neighbours. The top curve was generated from a run with K = 6. Both runs required approximately the same time to produce two nearly identical databases, each containing about 29,000 pathways. However, as can be seen from Figure 3.7, random perturbation of all the degrees of freedom results in uncooperative rearrangements being found first, while employing correlated perturbations has the opposite effect.

Figure 3.10:The average value of Nc for LJ13 pathway databases as new paths are added. 6-atom correlated perturbations (top curve) and random perturbations of every degree of freedom (bottom curve) were used to produce starting points for refinement by eigenvector-following [19222472]. Average values were calculated every time 100 new pathways were added.

3.5 Summary

The most important result of this chapter is probably the introduction of an index to quantify the cooperativity of atomic rearrangements. With this new measure it becomes possible to correlate cooperativity and barrier heights, and to show that cooperative rearrangements generally have lower barriers and shorter path lengths. We hope that these results will shed new light on relaxation mechanisms in complex systems, such as glasses and biomolecules, in future applications. For example, in a peptide or protein a large geometrical change can result from a rearrangement that could be described in terms of a single dihedral angle. In glasses and supercooled liquids an important research goal is to understand how observed dynamical properties, such as atomic diffusion and correlation functions [171, 174, 182, 183], are related to features of the underlying PES. The classification of elementary rearrangements ascage-breakingorcage-preserving [171, 184], and the emergence of structures such asmegabasins [171, 184186] can now be investigated more precisely in terms of localisation and cooperativity.

We have also demonstrated that cooperative rearrangements are relatively easy to characterise using double-ended transition state searching algorithms, since linear interpolation produces an effective initial guess. Uncooperative rearrangements are usually harder to find using such methods, and alternative initial guesses may be helpful in these cases.

Single-ended transition state searching has been used both in conjunction with double-ended methods, and as a way to sample potential energy surfaces for stationary points. Stationary point databases constructed using random perturbations followed by quenching are likely to be biased towards uncooperative rearrangements. We have therefore outlined a strategy for generating initial guesses appropriate to single-ended transition state searching algorithms, which instead favours cooperative rearrangements. This approach also includes a parameter that is likely to influence the degree of localisation.