Adaptive free energy sampling in multidimensional collective variable space using boxed molecular dynamics. Faraday Discussions, 195,

The past decade has seen the development of a new class of rare event methods in which molecular con ﬁ guration space is divided into a set of boundaries/interfaces, and then short trajectories are run between boundaries. For all these methods, an important concern is how to generate boundaries. In this paper, we outline an algorithm for adaptively generating boundaries along a free energy surface in multi-dimensional collective variable (CV) space, building on the boxed molecular dynamics (BXD) rare event algorithm. BXD is a simple technique for accelerating the simulation of rare events and free energy sampling which has proven useful for calculating kinetics and free energy pro ﬁ les in reactive and non-reactive molecular dynamics (MD) simulations across a range of systems, in both NVT and NVE ensembles. Two key developments outlined in this paper make it possible to automate BXD, and to adaptively map free energy and kinetics in complex systems. First, we have generalized BXD to multidimensional CV space. Using strategies from rigid-body dynamics, we have derived a simple and general velocity-re ﬂ ection procedure that conserves energy for arbitrary collective variable de ﬁ nitions in multiple dimensions, and show that it is straightforward to apply BXD to sampling in multidimensional CV space so long as the Cartesian gradients V CV are available. Second, we have modi ﬁ ed BXD to undertake on-the-ﬂ y statistical analysis during a trajectory, harnessing the information content latent in the dynamics to automatically determine boundary locations. Such automation not only makes BXD considerably easier to use; it also guarantees optimal boundaries, speeding up convergence. We have tested the multidimensional adaptive BXD procedure by calculating the potential of mean force for a chemical reaction recently investigated using both experimental and computational approaches – i.e. , F + CD 3 CN / DF + D 2 CN in both the gas phase and a strongly coupled explicit CD 3 CN solvent. The results obtained using multidimensional adaptive BXD agree well with previously published experimental and computational results, providing good evidence for its reliability


Introduction
The solution to a wide range of problems that can be addressed with molecular simulation consists fundamentally of determining rate coefficients.For example, biochemical systems rely on a delicate balance of rate coefficients within larger coupled kinetic networks. 1Similarly, bulk oxidation timescales in atmospheric chemistry 2 and combustion 3 (required to predict pollutant lifetimes, or to optimize an engine) are linked to detailed kinetic networks comprised of a wide range of elementary kinetic steps. 4With developments in both statistical mechanics and electronic structure theory, it is now possible to identify the important stationary points on a molecular potential energy surface (PES), 5 and carry out accurate calculations of the energy and partition function at each point.This enables extremely accurate calculations of the rates at which small molecules undergo structural changes, in either canonical or microcanonical ensembles. 6owever, calculating accurate rate coefficients for larger molecules (e.g., enzymes, long-chain hydrocarbon fuels, unsaturated volatile organic pollutants, etc.) remains an outstanding challenge for a number of reasons: (1) it remains difficult to calculate an accurate PES along a given path, (2) there is a combinatorial explosion in the number of paths with increasing system dimensionality, and (3) the conformational exibility inherent in larger molecular systems makes it very difficult to calculate accurate partition functions.Particularly as a result of the latter two challenges, the calculation of rate coefficients in complex systems tends to not to focus on stationary points, but rather on free energy surfaces along a particular path between states, typically dened in terms of a small set of collective variables (CVs).In cases where it is a good assumption that the full system dynamics along a particular path is mostly associated with changes in a small set of CVs, then the maximum on the free energy surface may be utilized to calculate rate coefficients in the Eyring equation. 7In cases where this is not a good assumption, an additional correction in the form of the so-called 'recrossing coefficient' is typically applied. 8,9n this paper, we present a relatively simple adaptive algorithm for discovering minimum free energy pathways between states in a multidimensional space of CVs, which can then be used to calculate rate coefficients in complex systems.There is strong evidence within computational complexity theory that problems of this sort are NP-complete [10][11][12] i.e., it is possible to verify (within polynomial time) whether any proposed solution is indeed a solution, but there is no known polynomial time algorithm to nd a solution in the rst place.This has rather profound consequences for how we think about free energy path sampling in complex molecular systems: the emphasis is less on nding an algorithm which is well-suited to every type of rare event problem, but rather on having access to a exible range of methods which can be practically used to tackle different conformational search problems.
'Boxed Molecular Dynamics' (BXD), [13][14][15][16] a method we have been actively involved in developing over the last few years, allows one to obtain both thermodynamic and kinetic information from the same run, producing data that produces a Markov master equation. 1,4,14,17BXD can be formulated so as to conserve energy, accelerating NVE simulations as well as NVT simulations.22]24 The fundamental idea in BXD is to accelerate dynamics simulations by introducing a set of hard boundaries within the hyperdimensional conguration space of the system being simulated.When a trajectory passes a boundary, those components of the velocity vector that take the trajectory across the boundary are reected.The statistics of reections at the boundary of the box are subsequently used to renormalize the results.Within BXD, 'boxes' refer to the conguration space domain between a particular set of boundaries.In principle, it is possible to implement boundaries which depend on the 6n dimensional phase space of Cartesian coordinates and momenta (where n is the number of atoms); however, in practice the original implementations of BXD utilized one-dimensional CVs in conguration space.
BXD falls within a class of sampling methods in which molecular conguration space is divided into a set of boundaries (also called interfaces or hypersurfaces), and short trajectories are run between boundaries.These methods (e.g., milestoning, 31,32 forward ux sampling, 33,34 transition interface sampling, 35 nonequilibrium umbrella sampling, 36 and others [37][38][39] ) have yet to displace umbrella sampling 40 as the most widely used method to determine free energies (or potentials of mean force), but in fact they have a number of features which we believe make them more attractive than umbrella sampling: (1) because they do not require modication of the potential energy function, they perturb the dynamics far less than umbrella sampling; (2) they allow for exact renormalization of the results in each box (unlike the iterative numerical WHAM scheme typically utilized to renormalize umbrella sampling results); (3) they require specication of fewer parameters than umbrella sampling (i.e., BXD only requires specifying a boundary location; umbrella sampling requires specifying the umbrella position and force constant); (4) they can provide both thermodynamic (free energy) and kinetic (rate) data simultaneously; (5) unlike umbrella sampling, they provide results which are in fact dynamically meaningful; and (6) it is possible to rigorously dene the regimes in which the accelerated dynamics they provide map onto the results that would have been obtained using standard unbiased simulations with standard initial conditions sampling strategies.
An important concern with these methods is how to generate boundaries (analogous to the umbrella sampling issue of how best to choose 'umbrella' potentials).In a broad range of molecular simulation studies, boundaries (or umbrella potentials) are located along a particular set of CVs which align with the intuition of the investigator (i.e., "user").For example, in enzyme catalysis, it is usually possible to highlight a few key bonds as being particularly important; similarly in a drug binding study, it is oen possible to identify a few key motions as particularly important to binding.Such user intuition is not a panacea: it may in fact fail to identify important CVs, and there are potential pitfalls 41 owing to the fact that it is oen extremely difficult to nd good CVs. 42Nevertheless, for understanding dynamics in hyperdimensional systems, user 'intuition' as to the important CVs usually constitutes an important guess as to where to initiate sampling and make practical progress in a simulation study.
With BXD's implementation in the CHARMM molecular simulation package, 43 it has found application to a range of chemical systems.These applications have highlighted two important issues: (1) the BXD velocity reection procedure must be generalized to deal with a wider range of CVs than the relatively small subset with which it is currently compatible; (2) with the implementation of a wider range of CVs, BXD must be formulated in a way that can automatically identify optimal boundaries in multi-dimensional collective variable space.The reason for the latter point is that BXD in multi-dimensional CV space requires specifying a large number of parameters.The number of parameters required scales as N CV Â N B , where N CV is the number of collective variables, and N B is the number of boundaries (N B is typically between 10 and 100 in systems studied so far).For relatively small systems where N CV ¼ 1 and which require no more than $10 boundaries, a user can typically keep track of the number of parameters requiring specication; however, for larger systems where N CV > 1, the number of parameters which requires specication rapidly expands beyond what even an expert can keep track of, becoming extremely tedious (if not altogether impossible).By automating the boundary selection scheme outlined in this paper using an adaptive algorithm, we avoid these problems entirely, and we also guarantee the specication of optimal boundaries.Adaptive sampling strategies have been previously explored in the context of umbrella sampling, 44,45 force biasing, 46 weighted ensemble sampling, 39 transition interface sampling, 47 accelerated molecular dynamics, 48 metadynamics 49,50 and steered MD. 51 BXD's robustness arises in part from the fact that it generates free energy proles which are largely insensitive to the location of boundaries, so long as the typical transit time from one boundary of the box to the other is larger than the system's characteristic decorrelation timescale. 13,14This is in fact the only 'hardand-fast' rule which must be satised in order for BXD to yield physically meaningful results: the average time between boundary reections in any given box must be larger than the system's characteristic dynamical decorrelation timescale in that region of the free energy surface. 14This rule places a lower limit on the allowed distance between any box's boundaries; otherwise, ballistic reection between box boundaries will occur, and the results are meaningless.So long as the boundaries are far enough apart to avoid problems related to dynamical decorrelation, then the choice of box boundaries is exible, and the BXD results do not depend on boundary location.
However, the computational efficiency of BXD (i.e., the speed at which it converges a free energy or a rate calculation) does depend on the boundary placement.For maximum efficiency, the boundaries should be placed close enough together so that a typical trajectory will visit the boundaries of any given box in a reasonable amount of time.Optimally placed boundaries will result in faster convergence.This is an issue that has become particularly apparent as we have attempted to use BXD to accelerate dynamics obtained from on-the-y electronic structure theory, and also in condensed phase reactions, where force evaluations are very expensive.Our experience to date has shown that 'userselected' box boundaries are oen far from ideal, and can result in wasted clock cycles, a point which is easily understood from Scheme 1.In regions with a large gradient, boxes should be smaller, given that an unbiased trajectory free to sample the box is more likely to get trapped downhill rather than travel uphill, while in atter regions that have a small gradient, the boxes can be larger, given that an unbiased trajectory will more readily sample wide regions of the conguration space.Scheme 1 therefore allows us to understand how clock cycles are wasted as a result of two common boundary-selection pitfalls: (1) large boxes in a region of the free energy surface with steep gradients, or (2) small boxes in a relatively at region of the free energy surface.In the former case, the trajectory will rarely visit high free energy congurations within the box, and convergence will be slow.In the latter case, clock cycles are wasted on constraining sampling in at regions of the free energy surface that the trajectory would have naturally visited anywayi.e., the boundaries actually slow down an intrinsic sampling rate which was already satisfactory.Scheme 1 highlights a nal important pointi.e., sampling on any given free energy path oen requires boxes of varying sizes, with the size of the box inversely related to the gradient of the free energy surface along a particular coordinate, which is generally unknown in advance.Box boundary placement is also sensitive to the local friction regime in which the dynamical process of interest takes place (a point discussed in further detail below).In general, efficient sampling in high friction environments (e.g., a chemical reaction occurring in a solvent) requires closely spaced boundaries, while boundaries in low friction environments (e.g., a chemical reaction in the gas phase) are farther apart.
In this paper, we outline an extension of BXD to multidimensional CV space, and an automated procedure that adaptively generates optimal BXD hypersurfaces to sample dynamical pathways within a user-specied multidimensional CV space.The underlying idea guiding this approach is simple, and exploits one of the key advantages of BXD compared to a method like umbrella sampling: because the underlying dynamics are in fact meaningful, 'on-the-y' analysis of their information content is in fact the most reliable guide to boundary placement.This philosophy allows us to use BXD for generating optimal boundaries in multidimensional applications, which may be subsequently used to accelerate rare events or carry out free energy sampling.We also report on results using this multi-dimensional adaptive BXD scheme to accelerate free energy sampling along the F + CD 3 CN / DF + CD 2 CN reactive pathway, in both CD 3 CN solvent and in the gas phase.This system constitutes a stringent test of the methodology, owing to the extreme asymmetry of the PES either side of the transition state (TS)e.g., similar to that shown in Scheme 1.The results are in good agreement with previous experimental and modelling studies, providing good evidence for the reliability of our extended BXD algorithm.We believe that the adaptive scheme described in this article may be useful to other methods that rely on sampling between conguration space interfaces.
2 Theoretical framework

BXD along a single collective variable
BXD is an exact extension of transition state theory, 13,14 with origins in Intramolecular Dynamics Diffusion Theory (IDDT), [52][53][54][55][56] which describes the motion of a trajectory along a reaction coordinate in terms of a diffusional equation or equivalent Langevin equation.BXD was initially formulated in order to accelerate dynamics by introducing a series of constraints along a one-dimensional collective variable, which provide a series of 'boxes' within which to lock the trajectory, as illustrated in Fig. 1.The region dened by the collective variable r is split into m boxes by the introduction of m + 1 user dened constraints.The trajectory is constrained within each box, which allows one to sample regions that would otherwise be visited only rarely.The trajectory constraint procedure involves an elastic collision procedure applied at the boundaries, which works as follows: whenever the next time step in the dynamics would result in the trajectory crossing the boundary, the trajectory is reset to the previous step, and a velocity inversion (i.e., reection) procedure is applied to those atoms that contribute to the denition of the collective variable.For a given box i bounded by r i and r iÀ1 , the rate coefficient for transfer from box i to i À 1 is determined by the inverse of the mean rst passage time (MFPT) hsi.The simplest way to compute this is to keep track of the number of times the trajectory has undergone velocity transformation at each boundary, h i,iÀ1 , along with the total amount of time, t i , that the trajectory spends within box i.This gives the rate coefficient for transfer from box i to box i À 1 as follows: Fig. 1 Illustration of the original one-dimensional BXD scheme along some collective variable r.The equilibrium constant between box i and box i À 1 may then be obtained from equilibrium statistical mechanics as where DG iÀ1,i is the free energy difference between box i and box i À 1. Eqn (2) allows us to obtain a full set of box-to-box free energy differences.Dening some arbitrary zero G 0 , the set of box-to-box free energy differences may then be summed appropriately to obtain DG i , the free energy of any given box relative to G 0 .This allows calculation of p i (the probability of the residing in box i) as follows: Having determined the probability of residing in any specic box according to (3), it is then possible to determine p(r) to arbitrary resolution by renormalizing the statistics within each box using histogram binning.Letting p i (r) be the probability of a particular value of r observed in box i, estimated by histogram binning from a sample within the box, then the probability of residing anywhere along the reaction coordinate dened by the boxes is given by Since only the box-to-box rate coefficients need to be computed, the length of time the trajectory needs to spend in each box is only determined by how long it takes for these rate coefficients to converge.The BXD method of partitioning the conguration space allows regions that are poorly sampled in standard MD trajectories to be isolated within a box and sampled independently, which lends itself well to parallelisation on modern cluster architectures.Alternatively, it is easy to formulate the BXD algorithm so that a given trajectory -aer a specied number of reection events at a particular boundaryis allowed to proceed to the next box, as illustrated in Fig. 1.Such a 'box-to-box' strategy allows trajectories to scan over adjacent boxes until convergence is achieved.

Extending BXD to multidimensional collective variable space
In this section, we present a generalisation of BXD to multidimensional collective variables.For a system of N atoms, we dene r(t) ˛R3N to be the vector of Cartesian coordinates of atoms in the system, and ṽ(t) ˛R3N to be the vector of corresponding velocities.A collective variable at some time t is a function s(t) of r(t) andṽ(t).In cases where one wants to characterize the dynamics of a molecular system at some time t using M collective variables, then the CV space may be represented as an M-dimensional vector s(t) ¼ [s 1 (t), s 2 (t), ., s M (t)], where M is generally much less than N.In the simplest case, where M ¼ 1, s(t) is oen referred to as a reaction coordinate.In its original implementation, BXD partitioned a one-dimensional collective variable space into an ordered set of zerodimensional points along the reaction coordinate.An intuitive route to generalising BXD is thus to partition the M-dimensional CV space into a series of (M À 1) dimensional boundaries, which is a strategy that follows naturally from BXD's origins in transition state theory. 54,57For example, a two-dimensional CV space may be partitioned by an ordered set of lines, a three-dimensional CV space by an ordered set of planes, and so onto the general case of hyperplanes.Within an Mdimensional collective variable space s(t) ¼ [s 1 (t), s 2 (t), ., s M (t)], any given BXD boundary B j may be dened as a plane in Hessian normal formi.e., in terms of a unit norm ñ ¼ [n 1 , n 2 , ., n M ] and a constant D j : Using the notation outlined above, Fig. 2 schematically illustrates a set of BXD boundaries that one might choose in order to partition a system dened in terms of two collective variables.

General velocity reection procedure in multidimensional collective variable space
Having specied a set of boundaries which partition the space of collective variables into smaller regions, a standard MD trajectory is performed within boundaries B j and B jÀ1 at every step, the collective variable vectors(t) is computed, and the velocities and positions of the previous time step are stored.For times t where the trajectory crosses either boundary B j or B jÀ1 , a velocity reection procedure is applied to constrain the trajectory so that it does not cross the boundary.In what follows, we focus on the velocity reection procedure to be used for reecting off multi-dimensional boundaries of the sort dened in eqn (5), generalizing the one-dimensional velocity reection procedure outlined in our previous BXD papers to multidimensional collective variable space.
Within the space of collective variables, eqn (5) species that a BXD boundary B j is dened in terms of a unit norm ñj ˛RM , which lies a distance D j from the origin.The function f(r(t)) ¼s(t)$ñj + D j provides a measure of how far the system Fig. 2 Schematic illustration of BXD boundaries that one could choose to partition a multi-dimensional system with two potential energy surface (PES) wells.The potential energy isosurface in the figure is projected into the collective variable s is from a particular boundary at time t, with changes in the sign of f(r(t)) indicating that the system has crossed B j .In order to constrain dynamics so that they lie to a specic side of a particular boundary B j , we wish to satisfy the following inequality: The inequality in this equation gives it the form of a so-called "unilateral constraint" 58 i.e., a constraint which is enforced only at times when the inequality is unsatised.For example, consider a case where f(r(t)) $ 0 at time t, and f(r(t + Dt)) < 0 at the next timestep t + Dt.In this case, the BXD procedure species that we revert back to r(t), and invert the velocities to give new velocities ṽ 0 (t), propagation according to which ensures that the constraints are satised at timestep t + Dt.By the chain rule, the time derivative of the constraint may be written as the projection of the atomic velocities onto the gradient of f(r(t)): To ensure that the constraint will be satised at time t + Dt, the inverted velocities must satisfy the following: In the general case of a system of K constraints, Vf is a matrix of K rows by 3N columns, but here we are restricting ourselves to the case of a single constraint, and therefore Vf in eqn (8) represents a row vector.The inverted velocities are related to the unbiased ones through eqn (8) in order to ensure a fully elastic reection of the velocities normal to B j .This procedure is in contrast to the sorts of holonomic constraints typically employed in molecular dynamics (e.g., SHAKE 59 and RATTLE 60 ), in which velocities normal to the constraint are set to zero in order to constrain the dynamics.The equation of motion for dynamics 58,[60][61][62] under a single constraint may be written as: where M ˛R3NÂ3N is a diagonal matrix of atomic masses, ã ˛R3N is the vector of accelerations, F is the force vector from the MD energy function, and G are the forces due to the constraint, given by where l is a time-dependent Lagrangian multiplier, and f T represents the transpose of f.Rather than applying the constraint directly as an acceleration, the constraint is enforced upon the inverted velocities as follows: By substituting eqn (11) into eqn (8) and rearranging for l we have The Lagrangian multiplier and subsequent impulse is only computed and applied for time steps in which an unaltered velocity would result in the constraint being unsatised, similar to the strategy used in the original BXD velocity reection algorithm.Dening BXD boundaries as hyperplanes ensures that the derivatives of f in eqn ( 12) may be computed by combining the derivatives of the components of s as follows: Eqn (13) means that the reection procedure can easily be constructed from a linear combination of derivatives of collective variables.This allows for straightforward combination of arbitrary reaction coordinates for which gradients are dened.The appendix to this paper includes an illustrative example of how to implement a velocity inversion procedure in the space of two CVs.

Adaptively generated boundaries in multidimensional CV space
The extension of BXD to multidimensional collective-variable space raises interesting questions as to where initial boundaries should be placed.For studies involving only a single collective variable (i.e., a 1d case), determining those boundary placements which most efficiently partition a reaction coordinate (either for rare event acceleration or free energy sampling) has generally been undertaken through some combination of 'user intuition' and trial and error.In multidimensional collective variable space, such a strategy quickly becomes unfeasible owing to the fact that the number of variables required to specify a boundary increases with the dimensionality of the CV space.In this section we present an automated adaptive path sampling procedure, in which optimal boundaries are generated through on-the-y statistical analysis carried out during a trajectory.
2.4.1 Overall adaptive scheme.Whereas previous implementations of BXD required a list of box boundaries, the adaptive implementation of BXD requires the user to provide the following input data, all of which are schematized in Fig. 2: (1) Specication of the CVs which the user wishes to adaptively sample along with a pair of limits that bound the sampling within a particular CV.In many cases, one of the CV limits (e.g., B 0 in Fig. 2) helps dene the extremum for what can be considered a reactant state, and the other CV limit (e.g., B n in Fig. 2) helps dene the extremum for what can be considered a product state.
(3) A 'target' or 'product' geometry (characterized by a set of 'target' CVs).We dene G 3 R M to be the region of CV space dened by two boundaries B R and B P (in Fig. 2, B R h B 0 and B P h B n ), and B i to be some arbitrary boundary that lies within G.The region G 1 3 R M lies between B R and B i , while the region G 2 3 R M lies between B i and B P , with G 1 + G 2 ¼ G.The approach of adaptive BXD is to carry out a single sampling run that makes two passes over the CV spacei.e., from B R to B P , and then to reverse direction and go from B P to B R .Along the way, statistical analysis determines the most efficient location at which to place the next bound.Aer the placement of a bound, the BXD velocity reection procedure is used to enhance the sampling of the next region.Passes in both directions are generally required so that barriers on the energy landscape are sampled in both directions (i.e., what may not require any acceleration going downhill will in fact require acceleration when going uphill).
The overall adaptive BXD procedure is illustrated by the owchart in Fig. 3, which assumes that the adaptive procedure has been initialised near B R , so that the rst pass involves generating boundaries en route to B P .At the start of the trajectory, B i ) B R , and B End ) B P (i.e., B i and B End initially enclose the region G h G 2 , with G 1 ¼ 0).Aer n steps of dynamics, samplings within G 2 (constrained through application of the BXD velocity reection procedure), there are two possible outcomes: (1) velocity reections against B End were observed, implying that the path from B i to B End requires no additional acceleration, or (2) velocity reections against B End were not observed, which means that an additional bound B new is required according to the procedure outlined below in section 2.4.2.Dynamics are then run until the system crosses B new , at which point B i ) B new .The dynamics in this hitherto unexplored space are then restricted in the region of G 2 through the application of the BXD velocity reection procedure.The sampling procedure is repeated until the dynamics reach B End , at which point B i Fig. 3 Flowchart illustrating the adaptive BXD boundary generation procedure.
) B P and B End ) B R , and the dynamics sweep back for a second pass in the opposite direction.Upon completion of the reverse pass, the Fig. 3 schematic arrives at the "Stop" point, and G will have been partitioned into a set of boxes with bounds B R , B 1 , B 2 , ., B n , B P for subsequent use in BXD runs to generate free energy surfaces to a specied degree of convergence.
2.4.2Procedure for adaptively generating a new boundary.An important aspect of the adaptive BXD scheme is 'on-the-y' analysis of the statistics collected during the sampling procedure for generation of a new boundary, B new .Let S RnÂM be the set of the sampled values of s, illustrated as blue circles in Fig. 4A, and let R ˛Rn be the vector of distances r from B i to each sampled values in S. The vector R provides information on how far from B i the next boundary should be placed, the location of which is determined as follows: (1) A normalized histogram of R is computed to give p(r), a probability density function representing the distances from B i that a trajectory samples between reections, as shown in Fig. 4B.
(2) From p(r), we calculate the cumulative distribution function pðrÞ.We then identify a histogram bin b max in p(r) with a bin centre r max chosen so that P(r max ) $ (1 À 3). 3 ˛(0, 1) is a parameter which species the "probability threshold" at which to place a new boundary (the value of 3 is specied by the user, and typically ranges from 0.01-0.1).We then identify smax (the mean value of the sampled values in S that fall within bin b max ), illustrated in Fig. 4B, as the point at which to place a new boundary B new .
(3) To determine the orientation of B new as an (M À 1)-dimensional plane, we use a simple strategy consistent with BXD's origins in transition state theory (TST) i.e., B new should be more or less orthogonal to the path of the observed dynamics. 13,14With b min dened as the rst bin in the histogram of R (see Fig. 4B), we calculates min (the mean value of the sampled values in S that fall within b min ).With this denition, smin represents the average value of s immediately prior to and aer reection against B i , i.e. the mean crossing point through B i .Similarly, smax represents the average crossing point through B new .The vector from smin to smax thus serves as an approximate dynamical pathway through the box, and we dene the unit norm for B new as The unit norm in eqn ( 14), combined withs max , allow us to fully dene the new boundary B new , as illustrated in Fig. 4C.The next time the trajectory crosses B new , it becomes enforced as a constraint (i.e., B i ) B new ), and an identical analysis will be carried out to determine the next B new .
As discussed above, adaptive boundary generation in this fashion will eventually lead to reection against B P .Once a trajectory reaches the barrier via adaptive boundary generation on the reactant side of the barrier, reection against B P generally follows rapidly without any need for boundaries on the product side of the barrier.To generate boundaries on the product side, a second adaptive sweep from B P to B R is required.To do this, the direction of sampling is reversed, and the adaptive boundary generation process is repeated going the opposite way.The only difference is thatbecause adaptive boundaries are already in place on the reactant side of the barrierthe reactant region is unlikely to require any more boundaries on the second sweep.For example, consider a BXD trajectory on its second sweep which is passing through the reactant region enclosed by boundaries B i and B iÀ1 (both of which were adaptively generated in the rst pass).It is likely for reection events against B iÀ1 to be observedi.e., sampling within this region is already suitably accelerated by BXD, and the trajectory can move on to the region dened by boundaries B iÀ1 and B iÀ2 .Should we observe that the trajectory has not inverted against B iÀ1 aer n steps, then an additional boundary is adaptively generated as described above.

Multidimensional adaptive sampling of chemical reactions in liquids
As an initial test of the multidimensional adaptive BXD scheme outlined above, we investigated F + CD 3 CN / DF + CD 2 CN in CD 3 CN solvent.This system has recently been the subject of both ultrafast transient IR spectroscopy experiments and corresponding non-equilibrium MD simulations. 63,64As such, it provides an excellent test case for investigating the algorithms described above, and also for evaluating their performance and accuracy.The reaction, which takes place in deuterated acetonitrile solvent (CD 3 CN), consists of deuterium abstraction from acetonitrile by the uorine atom, snapshots of which are illustrated in Fig. 5. Reactive molecular dynamics are possible using a customized version of the CHARMM molecular dynamics soware suite, using a parallel implementation of the multi-state empirical valence bond (MS-EVB) method.The simulation includes a single F radical embedded in a periodic box of 62 CD 3 CN solvent molecules.With a total of 64 MS-EVB states parallelized across 64 CPU cores, our simulations are able to treat the reactive process leading to DF as well as transient deuterium transfer from the nascent DF to the nitrile group on the other solvent molecules. 64The MS-EVB coupling elements were t to explicitly correlated CCSD(T)-F12 electronic structure theory extrapolated to the innite basis set limit (the contours of this PES are shown in Fig. 6).This procedure yields an accurate reactive PES, which is critical to understanding non-equilibrium energy deposition for reactions of this sort.
A one-dimensional implementation of BXD was previously used in these simulations to restrict the distance between the F radical and the reactive deuterium between 1.5 Å and 1.8 Å.This prevented the F radical diffusing away from its co-reactant during pre-production equilibration sampling runs.In the production NVE runs, the lower bound was removed.This allowed the reaction to occur, so that we could obtain an accurate measurement of energy deposition and relaxation in the nascent reaction products.At the time these studies were published, it was not to use BXD to generate a free energy surface for this reaction given that reversible reactive sampling requires the use of at least two CVs: the distance between the F radical and deuterium (F-D distance), and the distance between transferring deuterium and the carbon atom to which it is bonded (C-D distance).Constraint of the F-D distance accelerates abstraction over a relatively early barrier, and constraint of the C-D distance prevents the product DF from immediately diffusing away from its co-product.In addition to BXD sampling of the condensed phase reaction, we also carried out BXD sampling of the gas phase reaction, which included only three EVB states: the reactant F + CD 3 CN state, the co-product DF + CD 2 CN state, and the [CD 3 CND] + /[F] À state.Unless stated otherwise, all the results presented herein were run with a time step of 0.1 fs, using a Langevin thermostat at 300 K with a friction coef-cient of 20 ps À1 .

Adaptively generated BXD boundaries along the F + CD 3 CN reaction path
We applied the adaptive boundary generation procedure described in section 2.4 to sample this reaction and create BXD boundaries that could be used to accelerate the calculation of a free energy surface.This constitutes an interesting and particularly stringent test of our adaptive BXD procedure because of the large change in gradients along the reaction pathway: the gradients on the reactant side of the TS are very at, while those on the product side are very steep.Application of adaptive BXD to this system also enables us to comment on an outstanding experimental questionnamely, to what extent does the free energy surface of the gas phase chemical reaction resemble the free energy surface of the reaction in a strongly coupled solvent like CD 3 CN?In the gas phase, the 0 K reaction enthalpy is À37 kcal mol À1 , most of which is potentially available for deposition into the nascent DF product.Measurements carried out using ultrafast transient IR spectroscopy in solution showed deposition of substantial vibrational energy (i.e., at least v ¼ 2) in the stretching motion of the nascent DF product for the reaction taking place in solution.This value places a rm lower limit on exothermicity of the reaction free energy; however, a detailed analysis of the free energy proles in both the gas phase and in solvent is beyond experimental reach.As outlined above, adaptive BXD free energy sampling was undertaken in a CV space comprised of the F-D and C-D distances: an F-D distance of 2.7 Å was used to dene B R , and a C-D distance of 3.5 Å was used to dene B P .Adaptive sampling times of 100 ps (in the gas phase) and 30 ps (in solution phase) per box were used to determine the placement of new boundaries, with 3 ¼ 0.01 (guaranteeing that new boundaries are placed at a location visited no more than 1% of the time).Fig. 6 shows a series of snapshots taken during the automated boundary location procedure, illustrating how the adaptive algorithm works.Beginning from an initial point sampled near B R , BXD adaptively generates a boundary, which allows it to sample regions near the TS.Once the dynamics arrive at the TS, the system rapidly descends toward products, and quickly arrives at B P .At this point, BXD begins sweeping back in the opposite direction, adaptively generating boundaries which eventually return it back to the rst box bounded by B R .Fig. 7 shows the nal set of adaptively generated boundaries used to sample the free energy along the reaction pathway in both solution phase and in the gas phase.The plots also show the dynamical traces in CV space used to construct the BXD boundaries.
There are some important points to note with respect to Fig. 6 and 7: (1) the adaptively generated boundaries generally follow the route taken by the dynamics along the reaction pathway, with orientations that are roughly orthogonal to the dynamical pathway through CV space; (2) the spacing between boundaries varies as a function of the steepness of the free energy surface (the gradient) of the underlying PES and the corresponding free energy prolei.e., steep regions with large gradients require several boundaries, whereas less steep regions with smaller gradients require fewer boundaries; and (3) the reaction pathway in solvent has more adaptively generated BXD boundaries than the corresponding gas phase pathway, as a result of solvent friction effects that do not occur in the gas phase.Placing such a large number of BXD boundaries by user trial and error would be an extremely labour intensive process.

Free energy sampling within the adaptive boundaries
Having adaptively generated boundaries for both the solution and gas phase reactions, the standard BXD sampling procedure could then be applied.For the gas phase, the system is small enough that it was possible to gather all the required statistics with a single 100 ns trajectory, where the trajectory was sequentially restrained within each box until 100 reection events had occurred on either side of the boundary before being allowed through to the next box (another way of deciding how long to remain in the box is to monitor the point at which the MFPT reaches a user-specied convergence criterion).Given the larger size of the solution phase system along with the increased computational requirements that result from the 64 EVB states, we exploited the trivial parallelism of BXD to run trajectories in each box until meeting a user-specied convergence criterion (i.e., that the box-to-box MFPTs did not change by more than 0.1% with increased sampling), giving a total of 12 ns of dynamics across all boxes.Fig. 8 shows examples of the sampled values of the CVs obtained in the solution phase simulations, and demonstrates the sort of statistics obtained in two different regions along the free energy prole: (A) in the vicinity of the transition state, and (B) along a steep 'post-transition' state region aer DF has formed.
Once sampling was completed within each box (generating statistics similar to those shown in Fig. 7), MFPTs were calculated as described in section 2.1, and the results used to generate a 'box-averaged' free energy prole and corresponding 'box-averaged' probability spanning B R to B P .A higher-resolution free energy prole was obtained placing the statistics for a particular box into histogram bins and then using eqn (4) to renormalize by the box-averaged probabilities.Fig. 9A shows the smaller histogram bins into which we partitioned the statistics in each box to accomplish this.In the 1d case, high-resolution partitioning along the dynamical pathway is straightforward; in this case (and more generally for higherdimensional cases), our strategy is as follows: (1) Dene a path r which passes through the average dynamical crossing points through each boundary (i.e., smax ), and spans B R to B P .
(2) Each region between a set of boundaries is then partitioned into a series of bisecting hyperplanes, to an arbitrary user-specied resolution.The regions between these bisecting hyperplanes constitute the high-resolution histogram bins.The centre of each bin is chosen to be the point along r which is equidistant from the hyperplanes that bound the bin.
The blue line in Fig. 9A shows the path r which spans B R to B P , and which was used to generate ner histogram bins for plotting the high-resolution free energy prole.Fig. 9B shows the corresponding high-resolution BXD free energy proles for the gas phase reaction.The overlapping curves in this plot show how the free energy proles change with increasing sampling time, giving some indication of how quickly the BXD free energy prole converges in this particular system.Fig. 10 shows a comparison of the reactive free energy surfaces obtained in both the gas phase and in solvent.In the vicinity of the reactants and transition   states, the proles are very similar; however, they show considerable differences in the post-reaction region.The reason for this difference arises from post-reaction hydrogen bonding complexes formed by the nascent DF.In the gas phase BXD free energy sampling, the DF rotates around the backside of its CD 2 CN coproduct and nds a stable hydrogen-bonding complex with the nitrile moiety.In solvent, such interactions are possible with any of a wide range of nearby solvent molecules, and therefore no distinct minima can be observed along the free energy prole.In terms of understanding the DF energy deposition observed in the previously published experimental and MD results, the key quantity in Fig. 10 is the free energy difference between: (1) the maximum value observed near the transition state region, and (2) the minimum observed near the product state region.This quantity places an upper bound on the amount of energy which may be deposited into the nascent DF: for the reaction taking place in solvent the value is 27.6 kcal mol À1 , and for the reaction taking place in the gas phase the value is 31.1 kcal mol À1 .Both of these values are in good agreement with previous experimental and modelling studies.Our previously published experimental and MD studies indicated the prompt deposition of $23 kcal mol À1 into the stretching motion of the nascent DF prior to relaxation. 63,64While gas phase experiments of F + CD 3 CN are not available for direct comparison to our free energy results, experiments examining gas-phase energy deposition into HF in the F + CH 3 CN reaction have been performed, 65 and suggest that the nascent diatomic product in solvent contains slightly less excitation than in the gas phase. 63This is consistent with the results in Fig. 10, which indicate that more energy is available to the products in the gas phase reaction than in the solvent reaction.

Conclusions
In this paper, we have outlined an adaptive and automated procedure for generating boundaries in a multi-dimensional space of CVs.Our automated algorithm reduces the user effort required to carry out both rare event and free energy sampling in both one-dimensional and multi-dimensional cases; it generates box boundaries which are far enough apart to avoid any problems related to dynamical decorrelation, but which afford optimal acceleration.The extension of BXD to multidimensional collective variables provides an effective way to sample increasingly complex systems, but retains much of the simplicity and original properties of the 1-dimensional BXD implementation.The adaptive BXD scheme tested in this paper has been implemented in CHARMM, and will soon be available in the release version (we have also made initial efforts toward a BXD implementation in the TeraChem 66 ab initio dynamics package).The tests reported in this paper were carried out using the CHARMM implementation, in conjunction with parallelizable MS-EVB machinery also available in CHARMM. 64This framework allowed us to map free energy along a deuterium abstraction reaction pathway in both gas and solution phases.The results we obtained are in agreement with previously published experimental and modelling studies, providing good evidence for the reliability of our adaptive multidimensional BXD implementation.We believe that the adaptive scheme outlined in this paper, which allows us to generate hyperplanes in multi-dimensional collective variable space, may be more broadly useful to a wide range of techniques which rely on splicing up conguration space into a set of interfaces or boundaries.In the future, we will explore rigorous methods for estimating the error bars of free energy surfaces generated using BXD. 67We also plan to explore extensions of the adaptive BXD scheme in systems with CV spaces that have dimensionalities of three and highere.g., enzyme reactions and conformational dynamics, 68 drug binding, 69 and chemical reactions at surfaces and in liquids. 70As shown in eqn (7), implementation of BXD in multi-dimensional CV space requires denitions of the gradient in CV space, Vf, a wide library of which are available in the PLUMED 71 package.We are presently working on writing the BXD algorithm as a portable, and mostly 'standalone' plugin that may be easily interfaced with a wide range of molecular dynamics packages, a similar philosophy to that which has been adopted by PLUMED. 71Implementation of adaptive BXD in a package of this sort should allow it to be used in a wide range of contexts.
5 Appendix: velocity reflection in twodimensional CV space In this section we give details on the calculations required to perform velocity reection for a simple but illustrative case.Consider a system of atoms A, B and C where our collective variables are the distances AB and BC.This style of collective variable is useful in many situations, including the acceleration of abstraction reactions as discussed in the main document.
For the sake of brevity we restrict ourselves to 2 spatial coordinates.Letr ¼ [a x , a y , b x , b y , c x , c y ] be the coordinates and ṽ ¼ [V x a , V y a , V x b , V y b , V x c , V y a boundary being crossed, and we require a velocity reection.In order to perform the velocity reection using a Lagrangian multiplier, we need to compute Vf, the transpose of which is given by The expression above demonstrates how it is simple to construct the reection procedure from the gradients of the individual collective variables.With Vf in hand, it is a simple matter to determine the Lagrangian multiplier l with which the velocities may be inverted.From eqn (12) we may compute l via l ¼ À2Vf$ṽ VfM À1 Vf T ; (A.5) and then subsequently use it to compute inverted velocities from eqn (11) as ṽ0 (t) ¼ ṽ(t) + lM À1 Vf T .In the resulting velocities, the components normal to the boundary are inverted, and thus the constraint is satised.

Scheme 1
Scheme 1 Illustration of the relationship between a system's characteristic dynamics [red lines] in a given region of the free energy surface G(r) [black line] sampled along some CV r.Optimal boundaries are shown by grey lines.In steep regions of G(r), optimal boundaries are closely spaced; in flatter regions of G(r), optimal boundaries are farther apart.

Fig. 4
Fig. 4 Illustration of the adaptive boundary generation procedure.Panel A shows sampling of values within a 2d collective variable space, with an existing boundary B i ; panel B shows histogram binning of the distances with respect to the existing boundary B i .smin ands max , located within histogram bins b min and b max , are both shaded in red.Panel C shows generation of a new bound B new , where the norm defined in eqn (14) is illustrated by the purple arrow.

Fig. 5
Fig. 5 Snapshots from a molecular dynamics simulation of F + CD 3 CN in an explicit solvent of 62 CD 3 CN molecules.The images show: (1) approach of F to a CD 3 CN coreactant, (2) passage over the abstraction TS, (3) the nascent DF and its CD 2 CN coproduct, and (4) formation of a hydrogen-bonded complex between DF and another solvent molecule.

Fig. 6
Fig.6Time series illustrating the dynamical sequence that generates adaptive boundaries along the F + CD 3 CN reaction path in the gas phase.The grey dots indicate points in CV space that have already been sampled, and the black x indicates the position of the system at the time when the snapshot for each respective panel was taken.Snapshot 1 shows initial sampling near B R and snapshot 2 shows generation of the first boundary.Snapshot 3 shows the state of the system immediately following transition state passage and rapid downhill transit toward B P .Snapshots 4-6 show adaptive boundary placement as the system attempts to find its way back to the first box (i.e., that which is bounded by B R ).

Fig. 7
Fig. 7 Grey lines show the final set of adaptively generated BXD boundaries along the F + CD 3 CN reaction path in the gas phase, and in solution.The black traces show those values of the CVs which were sampled during the dynamics used to construct the boundaries.The contours indicate the underlying 0 K MS-EVB potential energy profile, and are provided for reference.The 0 K reaction enthalpy is À37 kcal mol À1 .

Fig. 8
Fig. 8 2D histograms of observed values of the collective variables from BXD sampling in solution.Panel A shows statistics sampled in the vicinity of the transition state, while panel B shows observed values on a steep 'post-transition state' region of the PES after DF has formed.The colors indicate the CV sampling frequency: dark red indicates a very high frequency, deep blue indicates a lower frequency, and white indicates zero frequency.

Fig. 9
Fig.9Panel A shows the BXD boundaries (black lines), and the high-resolution histogram bins (light grey lines) generated using the procedure outlined in the text for the gas phase reaction path.The blue line shows the path through the average dynamical boundary crossing points.Panel B shows the corresponding high-resolution BXD free energy profile for gas phase CD 3 CN.The overlapping curves in this plot show how the free energy profiles change with increasing sampling time, giving some indication of the rate of convergence for this particular system.

Fig. 10
Fig. 10 Reaction free energy profile in both gas and solution phases.
n 1 ða x À b x Þ=r AB n 1 À a y À b y Á r AB n 1 ða x À b x Þ=r AB þ n 2 ðb x À c x Þ=r BC n 1 À a y À b y Á r AB þ n 2 À b y À c y Á r BC n 2 ðc x À b x Þ=r BC