Tang Qing Yu

Microscopic mechanisms of equilibrium melting of a solid

Melting can follow many pathways Melting involves the loss of order as additional kinetic energy is added to a system. Although simple models of this sort of phase transition exist, it can be very difficult to observe the initial stages either experimentally or using simulations. Samanta et al. developed a robust rareevent sampling technique that makes it possible to examine melting events without needing excessive computing time (see the Perspective by van de Walle). For both copper and aluminum, they observed the formation of defects that act as starting points for the melting process rather than the homogeneous loss of order assumed in classic nucleation theory. Science, this issue p. 729 Multiple competing pathways direct a metallic solid to its molten state. [Also see Perspective by van de Walle] The melting of a solid, like other first-order phase transitions, exhibits an intrinsic time-scale disparity: The time spent by the system in metastable states is orders of magnitude longer than the transition times between the states. Using rare-event sampling techniques, we find that melting of representative solids—here, copper and aluminum—occurs via multiple, competing pathways involving the formation and migration of point defects or dislocations. Each path is characterized by multiple barrier-crossing events arising from multiple metastable states within the solid basin. At temperatures approaching superheating, melting becomes a single barrier-crossing process, and at the limit of superheating, the melting mechanism is driven by a vibrational instability. Our findings reveal the importance of nonlocal behavior, suggesting a revision of the perspective of classical nucleation theory.

Order-parameter-aided temperature-accelerated sampling for the exploration of crystal polymorphism and solid-liquid phase transitions

The problem of predicting polymorphism in atomic and molecular crystals constitutes a significant challenge both experimentally and theoretically. From the theoretical viewpoint, polymorphism prediction falls into the general class of problems characterized by an underlying rough energy landscape, and consequently, free energy based enhanced sampling approaches can be brought to bear on the problem. In this paper, we build on a scheme previously introduced by two of the authors in which the lengths and angles of the supercell are targeted for enhanced sampling via temperature accelerated adiabatic free energy dynamics [T. Q. Yu and M. E. Tuckerman, Phys. Rev. Lett. 107, 015701 (2011)]. Here, that framework is expanded to include general order parameters that distinguish different crystalline arrangements as target collective variables for enhanced sampling. The resulting free energy surface, being of quite high dimension, is nontrivial to reconstruct, and we discuss one particular strategy for performing the free energy analysis. The method is applied to the study of polymorphism in xenon crystals at high pressure and temperature using the Steinhardt order parameters without and with the supercell included in the set of collective variables. The expected fcc and bcc structures are obtained, and when the supercell parameters are included as collective variables, we also find several new structures, including fcc states with hcp stacking faults. We also apply the new method to the solid-liquid phase transition in copper at 1300 K using the same Steinhardt order parameters. Our method is able to melt and refreeze the system repeatedly, and the free energy profile can be obtained with high efficiency.

Locating landmarks on high-dimensional free energy surfaces.

Significance The problem of generating and navigating high-dimensional free energy surfaces is a significant challenge in the study of complex systems. The approach introduced represents an advance in this area, and its ability to generate and organize the key features of a high-dimensional free energy surface, i.e., its landmarks, with high efficiency impacts numerous problems in the materials and biomolecular sciences for which prediction of optimal structures is key. These include polypeptide and nucleic acid structure and crystal design and structure prediction. Moreover, as the algorithm targets the free energy surface, candidate structures can be ranked based on their relative free energies, which is not possible with algorithms that target only the bare potential energy surface. Coarse graining of complex systems possessing many degrees of freedom can often be a useful approach for analyzing and understanding key features of these systems in terms of just a few variables. The relevant energy landscape in a coarse-grained description is the free energy surface as a function of the coarse-grained variables, which, despite the dimensional reduction, can still be an object of high dimension. Consequently, navigating and exploring this high-dimensional free energy surface is a nontrivial task. In this paper, we use techniques from multiscale modeling, stochastic optimization, and machine learning to devise a strategy for locating minima and saddle points (termed “landmarks”) on a high-dimensional free energy surface “on the fly” and without requiring prior knowledge of or an explicit form for the surface. In addition, we propose a compact graph representation of the landmarks and connections between them, and we show that the graph nodes can be subsequently analyzed and clustered based on key attributes that elucidate important properties of the system. Finally, we show that knowledge of landmark locations allows for the efficient determination of their relative free energies via enhanced sampling techniques.

Multiscale implementation of infinite-swap replica exchange molecular dynamics

Significance Efficiently sampling the conformational space of complex molecular systems is a difficult problem that frequently arises in the context of molecular dynamics (MD) simulations. Here we present a modification of replica exchange MD (REMD) that is as simple to implement and parallelize as standard REMD but is shown to perform better in terms of sampling efficiency. The method is used here to investigate the multifunnel structure of the C-terminal β-hairpin of protein G in explicit solvent, which to date has required much longer REMD simulations to produce converged predictions of conformational statistics. In particular, we identify a potentially important folded β structure previously undetected by other importance sampling methods. Replica exchange molecular dynamics (REMD) is a popular method to accelerate conformational sampling of complex molecular systems. The idea is to run several replicas of the system in parallel at different temperatures that are swapped periodically. These swaps are typically attempted every few MD steps and accepted or rejected according to a Metropolis–Hastings criterion. This guarantees that the joint distribution of the composite system of replicas is the normalized sum of the symmetrized product of the canonical distributions of these replicas at the different temperatures. Here we propose a different implementation of REMD in which (i) the swaps obey a continuous-time Markov jump process implemented via Gillespie’s stochastic simulation algorithm (SSA), which also samples exactly the aforementioned joint distribution and has the advantage of being rejection free, and (ii) this REMD-SSA is combined with the heterogeneous multiscale method to accelerate the rate of the swaps and reach the so-called infinite-swap limit that is known to optimize sampling efficiency. The method is easy to implement and can be trivially parallelized. Here we illustrate its accuracy and efficiency on the examples of alanine dipeptide in vacuum and C-terminal β-hairpin of protein G in explicit solvent. In this latter example, our results indicate that the landscape of the protein is a triple funnel with two folded structures and one misfolded structure that are stabilized by H-bonds.

Kinetics of O2 Entry and Exit in Monomeric Sarcosine Oxidase via Markovian Milestoning Molecular Dynamics

The flavoenzyme monomeric sarcosine oxidase (MSOX) catalyzes a complex set of reactions currently lacking a consensus mechanism. A key question that arises in weighing competing mechanistic models of MSOX function is to what extent ingress of O2 from the solvent (and its egress after an unsuccessful oxidation attempt) limits the overall catalytic rate. To address this question, we have applied to the MSOX/O2 system the relatively new simulation method of Markovian milestoning molecular dynamics simulations, which, as we recently showed [ Yu et al. J. Am. Chem. Soc. 2015 , 137 , 3041 ], accurately predicted the entry and exit kinetics of CO in myoglobin. We show that the mechanism of O2 entry and exit, in terms of which possible solvent-to-active-site channels contribute to the flow of O2, is sensitive to the presence of the substrate-mimicking competitive inhibitor 2-furoate in the substrate site. The second-order O2 entry rate constants were computed to be 8.1 × 10(6) and 3.1 × 10(6) M(-1) s(-1) for bound and apo MSOX, respectively, both of which moderately exceed the experimentally determined second-order rate constant of (2.83 ± 0.07) × 10(5) M(-1) s(-1) for flavin oxidation by O2 in MSOX. This suggests that the rate of flavin oxidation by O2 is likely not strongly limited by diffusion from the solvent to the active site. The first-order exit rate constants were computed to be 10(7) s(-1) and 7.2 × 10(6) s(-1) for the apo and bound states, respectively. The predicted faster entry and slower exit of O2 for the bound state indicate a longer residence time within MSOX, increasing the likelihood of collisions with the flavin isoalloxazine ring, a step required for reduction of molecular O2 and subsequent reoxidation of the flavin. This is also indirectly supported by previous experimental evidence favoring the so-called modified ping-pong mechanism, the distinguishing feature of which is an intermediate complex involving O2, the flavin, and the oxidized substrate simultaneously in the cavity. These findings demonstrate the utility of the Markovian milestoning approach in contributing new understanding of complicated enyzmatic function.