Baron Peters

A growing string method for determining transition states: Comparison to the nudged elastic band and string methods

Interpolation methods such as the nudged elastic band and string methods are widely used for calculating minimum energy pathways and transition states for chemical reactions. Both methods require an initial guess for the reaction pathway. A poorly chosen initial guess can cause slow convergence, convergence to an incorrect pathway, or even failed electronic structure force calculations along the guessed pathway. This paper presents a growing string method that can find minimum energy pathways and transition states without the requirement of an initial guess for the pathway. The growing string begins as two string fragments, one associated with the reactants and the other with the products. Each string fragment is grown separately until the fragments converge. Once the two fragments join, the full string moves toward the minimum energy pathway according to the algorithm for the string method. This paper compares the growing string method to the string method and to the nudged elastic band method using the alanine dipeptide rearrangement as an example. In this example, for which the linearly interpolated guess is far from the minimum energy pathway, the growing string method finds the saddle point with significantly fewer electronic structure force calculations than the string method or the nudged elastic band method.

Obtaining reaction coordinates by likelihood maximization

We present a new approach for calculating reaction coordinates in complex systems. The new method is based on transition path sampling and likelihood maximization. It requires fewer trajectories than a single iteration of existing procedures, and it applies to both low and high friction dynamics. The new method screens a set of candidate collective variables for a good reaction coordinate that depends on a few relevant variables. The Bayesian information criterion determines whether additional variables significantly improve the reaction coordinate. Additionally, we present an advantageous transition path sampling algorithm and an algorithm to generate the most likely transition path in the space of collective variables. The method is demonstrated on two systems: a bistable model potential energy surface and nucleation in the Ising model. For the Ising model of nucleation, we quantify for the first time the role of nuclei surface area in the nucleation reaction coordinate. Surprisingly, increased surface area increases the stability of nuclei in two dimensions but decreases nuclei stability in three dimensions.

Molecular-level origins of biomass recalcitrance: decrystallization free energies for four common cellulose polymorphs.

Cellulose is a crystalline polymer of β1,4-D-glucose that is difficult to deconstruct to sugars by enzymes. The recalcitrance of cellulose microfibrils is a function of both the shape of cellulose microfibrils and the intrinsic work required to decrystallize individual chains, the latter of which is calculated here from the surfaces of four crystalline cellulose polymorphs: cellulose Iβ, cellulose Iα, cellulose II, and cellulose III(I). For edge chains, the order of decrystallization work is as follows (from highest to lowest): Iβ, Iα, ΙΙΙ(Ι), and II. For cellulose Iβ, we compare chains from three different locations on the surface and find that an increasing number of intralayer hydrogen bonds (from 0 to 2) increases the intrinsic decrystallization work. From these results, we propose a microkinetic model for the deconstruction of cellulose (and chitin) by processive enzymes, which when taken with a previous study [Horn et al. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 18089] identifies the thermodynamic and kinetic attributes of enzyme and substrate engineering for enhanced cellulose (or chitin) conversion. Overall, this study provides new insights into the molecular interactions that form the structural basis of cellulose, which is the primary building block of plant cell walls, and highlights the need for experimentally determining microfibril shape at the nanometer length scale when comparing conversion rates of cellulose polymorphs by enzymes.

Homogeneous nucleation of methane hydrates: unrealistic under realistic conditions.

Methane hydrates are ice-like inclusion compounds with importance to the oil and natural gas industry, global climate change, and gas transportation and storage. The molecular mechanism by which these compounds form under conditions relevant to industry and nature remains mysterious. To understand the mechanism of methane hydrate nucleation from supersaturated aqueous solutions, we performed simulations at controlled and realistic supersaturation. We found that critical nuclei are extremely large and that homogeneous nucleation rates are extremely low. Our findings suggest that nucleation of methane hydrates under these realistic conditions cannot occur by a homogeneous mechanism.

Extensions to the likelihood maximization approach for finding reaction coordinates.

This paper extends our previous work on obtaining reaction coordinates from aimless shooting and likelihood maximization. We introduce a simplified version of aimless shooting and a half-trajectory likelihood score based on the committor probability. Additionally, we analyze and compare the absolute log-likelihood score for perfect and approximate reaction coordinates. We also compare the aimless shooting and likelihood maximization approach to the earlier genetic neural network (GNN) approach of Ma and Dinner [J. Phys. Chem. B 109, 6769 (2005)]. For a fixed number of total trajectories in the GNN approach, the accuracy of the transition state ensemble decreases as the number of trajectories per committor probability estimate increases. This quantitatively demonstrates the benefit of individual committor probability realizations over committor probability estimates. Furthermore, when the least squares score of the GNN approach is applied to individual committor probability realizations, the likelihood score still provides a better approximation to the true transition state surface. Finally, the polymorph transition in terephthalic acid demonstrates that the new half-trajectory likelihood scheme estimates the transition state location more accurately than likelihood schemes based on the probability of being on a transition path.

Transmission Coefficients, Committors, and Solvent Coordinates in Ion-Pair Dissociation

From a hypothetical perfect dividing surface, all trajectories commit to opposite basins in forward and backward time without recrossing, transition state theory is exact, the transmission coefficient is one, and the committor distribution is perfectly focused at 1/2. However, chemical reactions in solution and other real systems often have dynamical trajectories that recross the dividing surface. To separate true dynamical effects from effects of a nonoptimal dividing surface, the dividing surface and/or reaction coordinate should be optimized before computing transmission coefficients. For NaCl dissociation in TIP3P water, we show that recrossing persists even when the 1/2-committor surface itself is used as the dividing surface, providing evidence that recrossing cannot be fully eliminated from the dynamics for any configurational coordinate. Consistent with this finding, inertial likelihood maximization finds a combination of ion-pair distance and two solvent coordinates that improves the committor distribution and increases the transmission coefficient relative to those for ion-pair distance alone, but recrossing is not entirely eliminated. Free energy surfaces for the coordinates identified by inertial likelihood maximization show that the intrinsic recrossing stems from anharmonicity and shallow intermediates that remain after dimensionality reduction to the dynamically important variables.

Competing nucleation pathways in a mixture of oppositely charged colloids: Out-of-equilibrium nucleation revisited

Recent simulations of crystal nucleation from a compressed liquid of oppositely charged colloids show that the natural Brownian dynamics results in nuclei of a charge-disordered FCC (DFCC) solid whereas artificially accelerated dynamics with charge swap moves result in charge-ordered nuclei of a CsCl phase. These results were interpreted as a breakdown of the quasiequilibrium assumption for precritical nuclei. We use structure-specific nucleus size coordinates for the CsCl and DFCC structures and equilibrium based sampling methods to understand the dynamical effects on structure selectivity in this system. Nonequilibrium effects observed in previous simulations emerge from a diffusion tensor that dramatically changes when charge swap moves are used. Without the charge swap moves diffusion is strongly anisotropic with very slow motion along the charge-ordered CsCl axis and faster motion along the DFCC axis. Kramers-Langer-Berezhkovskii-Szabo theory predicts that under the realistic dynamics, the diffusion anisotropy shifts the current toward the DFCC axis. The diffusion tensor also varies with location on the free energy landscape. A numerical calculation of the current field with a diffusion tensor that depends on the location in the free energy landscape exacerbates the extent to which the current is skewed toward DFCC structures. Our analysis confirms that quasiequilibrium theories based on equilibrium properties can explain the nonequilibrium behavior of this system. Our analysis also shows that using a structure-specific nucleus size coordinate for each possible nucleation product can provide mechanistic insight on selectivity and competition between nucleation pathways.

Reaction coordinates, one-dimensional Smoluchowski equations, and a test for dynamical self-consistency

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.

Nucleation in a Potts lattice gas model of crystallization from solution

Nucleation from solution is important in many pharmaceutical crystallization, biomineralization, material synthesis, and self-assembly processes. Simulation methodology has progressed rapidly for studies of nucleation in pure component and implicit solvent systems; however little progress has been made in the simulation of explicit solvent systems. The impasse stems from the inability of rare events simulation methodology to be combined with simulation techniques which maintain a constant chemical potential driving force (supersaturation) for nucleation. We present a Potts lattice gas (PLG) to aid in the development of new simulation strategies for nucleation from solution. The PLG captures common crystallization phase diagram features such as a eutectic point and solute/solvent melting points. Simulations of the PLG below the bulk solute melting temperature reveal a competition between amorphous and crystalline nuclei. As the temperature is increased toward the bulk melting temperature, the nucleation pathway changes from a one step crystalline nucleation pathway to a two step pathway, where an amorphous nucleus forms and then crystallizes. We explain these results in terms of classical nucleation theory with different size-dependant chemical potentials for the amorphous and crystalline nucleation pathways. The two step pathway may be particularly important when crystallization is favored only at postcritical sizes.

Optimizing Nucleus Size Metrics for Liquid–Solid Nucleation from Transition Paths of Near-Nanosecond Duration

We determine the mechanism for the liquid-solid phase transition in the Lennard-Jones fluid close to coexistence with aimless shooting and likelihood maximization. The reaction coordinate for this process is a product of a structural descriptor and the size of the nascent solid nucleus and is quantitatively verified with the committor probability histogram test. This study identifies the first accurate scalar reaction coordinate for the liquid-solid nucleation process in Lennard-Jonesium, which will likely extend to nucleation processes in other spherically symmetric fluids. On the basis of our results, we propose a structural correction factor for the commonly cited nucleus size reaction coordinate from classical nucleation theory that enables connection of simulation data to stochastic models of nucleation kinetics. In addition, we show that aimless shooting is able to obtain reasonable acceptance rates for transitions with highly diffusive characteristics, which has been problematic for transition path sampling methods for diffusive processes such as nucleation and macromolecular transitions.