^{1,2}, Peter G. Bolhuis

^{3}, Ryan G. Mullen

^{1}and Joan-Emma Shea

^{2}

### Abstract

We propose a method for identifying accurate reaction coordinates among a set of trial coordinates. The method applies to special cases where motion along the reaction coordinate follows a one-dimensional Smoluchowski equation. In these cases the reaction coordinate can predict its own short-time dynamical evolution, i.e., the dynamics projected from multiple dimensions onto the reaction coordinate depend only on the reaction coordinate itself. To test whether this property holds, we project an ensemble of short trajectory swarms onto trial coordinates and compare projections of individual swarms to projections of the ensemble of swarms. The comparison, quantified by the Kullback-Leibler divergence, is numerically performed for each isosurface of each trial coordinate. The ensemble of short dynamical trajectories is generated only once by sampling along an initial order parameter. The initial order parameter should separate the reactants and products with a free energy barrier, and distributions on isosurfaces of the initial parameter should be unimodal. The method is illustrated for three model free energy landscapes with anisotropic diffusion. Where exact coordinates can be obtained from Kramers-Langer-Berezhkovskii-Szabo theory, results from the new method agree with the exact results. We also examine characteristics of systems where the proposed method fails. We show how dynamical self-consistency is related (through the Chapman-Kolmogorov equation) to the earlier isocommittor criterion, which is based on longer paths.

This work was supported by NSF CAREER (Grant No. 0955502) and by a NSF Computational Discovery and Innovation (CDI) (Grant No. 1125235). We thank Professor Scott Shell whose work stimulated our thinking about the Kullback-Leibler divergence. We thank the reviewers for helpful suggestions. Finally, we thank David Wu, Amadeu Sum, Gregg Beckham, Valeria Molinero, and Frank Noe for stimulating discussions.

INTRODUCTION

THEORY

METHOD

ALGORITHM

RESULTS

Example I: Nucleation with anisotropic structure formation kinetics

The model

Simulations

Optimization of the reaction coordinate

Example II: A narrow tube landscape

The model

The effect of the choice of λ

Example III: When isosurfaces of λ give a bimodal distribution

CONCLUSIONS

### Key Topics

- Free energy
- 54.0
- Diffusion
- 35.0
- Anisotropy
- 19.0
- Chemical reaction theory
- 17.0
- Nucleation
- 14.0

## Figures

Schematic depicting a projection of short trajectory swarm data (gray) from an initial point **x** _{0} (black) to specific coordinates **q** _{1} and **q** _{2} at later times.

Schematic depicting a projection of short trajectory swarm data (gray) from an initial point **x** _{0} (black) to specific coordinates **q** _{1} and **q** _{2} at later times.

For a reaction coordinate whose dynamics follow a one-dimensional Smoluchowski equation, swarms of trajectories from different individual configurations on each isosurface will drift and diffuse similarly. Therefore, for each isosurface, the projection of the individual swarms should each resemble the combined projection of all swarms (depicted to the right of the free energy surface). In the case depicted above, the individual swarms behave differently from each other, and therefore differently from their combined projection.

For a reaction coordinate whose dynamics follow a one-dimensional Smoluchowski equation, swarms of trajectories from different individual configurations on each isosurface will drift and diffuse similarly. Therefore, for each isosurface, the projection of the individual swarms should each resemble the combined projection of all swarms (depicted to the right of the free energy surface). In the case depicted above, the individual swarms behave differently from each other, and therefore differently from their combined projection.

If the dynamics of *q* follow a one-dimensional Smoluchowski equation, dynamical self-consistency should apply at all isosurfaces of *q*. (Left) For a poor reaction coordinate *q*(**x**) the value of *q* alone is not a good predictor of drift or diffusion from the configuration **x**. (Right) If *q* shows dynamical self-consistency for all isosurfaces of *q* from reactant A to product B, then *q* is an accurate reaction coordinate.

If the dynamics of *q* follow a one-dimensional Smoluchowski equation, dynamical self-consistency should apply at all isosurfaces of *q*. (Left) For a poor reaction coordinate *q*(**x**) the value of *q* alone is not a good predictor of drift or diffusion from the configuration **x**. (Right) If *q* shows dynamical self-consistency for all isosurfaces of *q* from reactant A to product B, then *q* is an accurate reaction coordinate.

The solid curves are contours of the actual free energy landscape *βF*(*q* _{1}, *q* _{2}) with a saddle point at the round dot. The coordinate *q* _{1} has been used as the initial coordinate, i.e., λ(**x**) = *q* _{1}(**x**). The dotted contours are curves of constant Λ. Note that Λ is peaked where the original landscape had a saddle point.

The solid curves are contours of the actual free energy landscape *βF*(*q* _{1}, *q* _{2}) with a saddle point at the round dot. The coordinate *q* _{1} has been used as the initial coordinate, i.e., λ(**x**) = *q* _{1}(**x**). The dotted contours are curves of constant Λ. Note that Λ is peaked where the original landscape had a saddle point.

The free energy landscape for a model of nucleation where the nucleus size can change along either fast (*n* _{ F }) and slow (*n* _{ S }) mobility directions.

The free energy landscape for a model of nucleation where the nucleus size can change along either fast (*n* _{ F }) and slow (*n* _{ S }) mobility directions.

The free energy as a function of the initial coordinate λ = *n* _{ F } + *n* _{ S }.

The free energy as a function of the initial coordinate λ = *n* _{ F } + *n* _{ S }.

Histograms of endpoints of 1000 trajectories 2.0τ after initiation. Each figure shows nine swarms initiated at nine different points on the free energy landscape. In (a) the anisotropy is *s* = 1.0, and in (b) the anisotropy is *s* = 0.1. Differences in the way the swarms drift are difficult to visually discern, but the dynamical self-consistency test can detect differences.

Histograms of endpoints of 1000 trajectories 2.0τ after initiation. Each figure shows nine swarms initiated at nine different points on the free energy landscape. In (a) the anisotropy is *s* = 1.0, and in (b) the anisotropy is *s* = 0.1. Differences in the way the swarms drift are difficult to visually discern, but the dynamical self-consistency test can detect differences.

The trial reaction coordinate *q* can be represented by the angle θ between the direction of progress along *q* and the *n* _{ F }-axis.

The trial reaction coordinate *q* can be represented by the angle θ between the direction of progress along *q* and the *n* _{ F }-axis.

exp[−Λ] weighted Kullback-Leibler divergences for different trial reaction coordinates (represented by θ) and for different values of the diffusion anisotropy, *s* = 1.0, 0.3, 0.1, 0.03.

exp[−Λ] weighted Kullback-Leibler divergences for different trial reaction coordinates (represented by θ) and for different values of the diffusion anisotropy, *s* = 1.0, 0.3, 0.1, 0.03.

Illustrating three mechanistic regimes that prevail for different degrees of diffusion anisotropy. When *s* ≈ 1, the diffusion tensor is isotropic and most pathways follow the minimum free energy path. When *s* < 1, but not extremely small, a two step nucleation mechanism prevails with initial motion along the slow coordinate before escape in the *n* _{ F }-direction. Finally, when *s* is extremely small the Berezhkovskii-Zitserman (BZ) regime prevails. In the BZ-regime, trajectories can escape in the *n* _{ F }-direction with the *n* _{ S } degree of freedom frozen.

Illustrating three mechanistic regimes that prevail for different degrees of diffusion anisotropy. When *s* ≈ 1, the diffusion tensor is isotropic and most pathways follow the minimum free energy path. When *s* < 1, but not extremely small, a two step nucleation mechanism prevails with initial motion along the slow coordinate before escape in the *n* _{ F }-direction. Finally, when *s* is extremely small the Berezhkovskii-Zitserman (BZ) regime prevails. In the BZ-regime, trajectories can escape in the *n* _{ F }-direction with the *n* _{ S } degree of freedom frozen.

Narrow tube type free energy landscape as a function of fast (*n* _{ F }) and slow (*n* _{ S }) coordinates. The saddle point is at (*n* _{ F }, *n* _{ S }) = (20, 20).

Narrow tube type free energy landscape as a function of fast (*n* _{ F }) and slow (*n* _{ S }) coordinates. The saddle point is at (*n* _{ F }, *n* _{ S }) = (20, 20).

Projections of the free energy onto different initial coordinates. (a) The initial coordinate was λ = *n* _{ S }. (b) The initial coordinate was λ = *n* _{ F }.

Projections of the free energy onto different initial coordinates. (a) The initial coordinate was λ = *n* _{ S }. (b) The initial coordinate was λ = *n* _{ F }.

⟨Δ_{KL}[q]⟩_{Λ} for different coordinates (represented by the angle θ as shown in Figure 8 ). (a) The initial coordinate was λ = *n* _{ S }. (b) The initial coordinate was λ = *n* _{ F }. In both cases the optimal coordinate (the minimum) rotates toward *n* _{ S } as the diffusion tensor becomes more anisotropic.

⟨Δ_{KL}[q]⟩_{Λ} for different coordinates (represented by the angle θ as shown in Figure 8 ). (a) The initial coordinate was λ = *n* _{ S }. (b) The initial coordinate was λ = *n* _{ F }. In both cases the optimal coordinate (the minimum) rotates toward *n* _{ S } as the diffusion tensor becomes more anisotropic.

Contour plots of (a) the EVB potential V(x, y) with the minimum energy path between the reactant and product minima, (b) the potential Λ(x, y) constructed with λ chosen as the end-to-end direction along the minimum energy path, and (c) the potential Λ(x, y) constructed with the ideal energy gap coordinate for λ. Contour spacings are 5k_{B}T in all plots.

Contour plots of (a) the EVB potential V(x, y) with the minimum energy path between the reactant and product minima, (b) the potential Λ(x, y) constructed with λ chosen as the end-to-end direction along the minimum energy path, and (c) the potential Λ(x, y) constructed with the ideal energy gap coordinate for λ. Contour spacings are 5k_{B}T in all plots.

⟨Δ_{KL}[q]⟩_{Λ} for different linear trial coordinates q(x, y) represented by the angle θ. The initial coordinate λ(x, y) gives a hysteretic free energy F_{λ}(λ) if sampled imperfectly. Coordinates similar to λ clearly do not have dynamical self-consistency, but other coordinates are not clearly distinguished in ⟨Δ_{KL}[q]⟩_{Λ.}

⟨Δ_{KL}[q]⟩_{Λ} for different linear trial coordinates q(x, y) represented by the angle θ. The initial coordinate λ(x, y) gives a hysteretic free energy F_{λ}(λ) if sampled imperfectly. Coordinates similar to λ clearly do not have dynamical self-consistency, but other coordinates are not clearly distinguished in ⟨Δ_{KL}[q]⟩_{Λ.}

⟨Δ_{KL}[q]⟩_{Λ} for different linear trial coordinates q(x, y) represented by the angle θ. The initial coordinate λ(x, y) is the vertical energy gap between diabatic states of the EVB model.

⟨Δ_{KL}[q]⟩_{Λ} for different linear trial coordinates q(x, y) represented by the angle θ. The initial coordinate λ(x, y) is the vertical energy gap between diabatic states of the EVB model.

## Tables

Comparison between reaction coordinates from the dynamical self-consistency test and from KLBS theory. Coordinates identified by dynamical self-consistency are summarized for two different initial coordinates. For narrow tube potential energy landscapes the dynamical self-consistency test can correctly identify accurate reaction coordinates even for an inaccurate initial coordinate λ.

Comparison between reaction coordinates from the dynamical self-consistency test and from KLBS theory. Coordinates identified by dynamical self-consistency are summarized for two different initial coordinates. For narrow tube potential energy landscapes the dynamical self-consistency test can correctly identify accurate reaction coordinates even for an inaccurate initial coordinate λ.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content