Lianqing Zheng

Random walk in orthogonal space to achieve efficient free-energy simulation of complex systems.

In the past few decades, many ingenious efforts have been made in the development of free-energy simulation methods. Because complex systems often undergo nontrivial structural transition during state switching, achieving efficient free-energy calculation can be challenging. As identified earlier, the “Hamiltonian” lagging, which reveals the fact that necessary structural relaxation falls behind the order parameter move, has been a primary problem for generally low free-energy simulation efficiency. Here, we propose an algorithm by achieving a random walk in both the order parameter space and its generalized force space; thereby, the order parameter move and the required conformational relaxation can be efficiently synchronized. As demonstrated in both the alchemical transition and the conformational transition, a leapfrog improvement in free-energy simulation efficiency can be obtained; for instance, (i) it allows us to solve a notoriously challenging problem: accurately predicting the pKa value of a buried titratable residue, Asp-66, in the interior of the V66E staphylococcal nuclease mutant, and (ii) it allows us to gain superior efficiency over the metadynamics algorithm.

Melting of Cu under hydrostatic and shock wave loading to high pressures

Molecular dynamics simulations are performed to investigate hydrostatic melting and shock-induced melting of single crystal Cu described by an embedded-atom method potential. The thermodynamic (equilibrium) melting curve obtained from our simulations agrees with static experiments and independent simulations. The planar solid–liquid interfacial energy is found to increase with pressure. The amount of maximum superheating or supercooling is independent of pressure, and is 1.24 ± 0.01 and 0.68 ± 0.01 at a heating or cooling rate of 1 K ps−1, respectively. We explore shock loading along three main crystallographic directions: , and . Melting along the principal Hugoniot differs considerably from and , possibly due to different extents of solid state disordering. Along , the solid is superheated by about 20%, before it melts with a pronounced temperature drop. In contrast, melting along and is quasi-continuous, and premelting (~7%) is observed.

Practically Efficient and Robust Free Energy Calculations: Double-Integration Orthogonal Space Tempering.

The orthogonal space random walk (OSRW) method, which enables synchronous acceleration of the motions of a focused region and its coupled environment, was recently introduced to enhance sampling for free energy simulations. In the present work, the OSRW algorithm is generalized to be the orthogonal space tempering (OST) method via the introduction of the orthogonal space sampling temperature. Moreover, a double-integration recursion method is developed to enable practically efficient and robust OST free energy calculations, and the algorithm is augmented by a novel θ-dynamics approach to realize both the uniform sampling of order parameter spaces and rigorous end point constraints. In the present work, the double-integration OST method is employed to perform alchemical free energy simulations, specifically to calculate the free energy difference between benzyl phosphonate and difluorobenzyl phosphonate in aqueous solution, to estimate the solvation free energy of the octanol molecule, and to predict the nontrivial Barnase-Barstar binding affinity change induced by the Barnase N58A mutation. As demonstrated in these model studies, the DI-OST method can robustly enable practically efficient free energy predictions, particularly when strongly coupled slow environmental transitions are involved.

Homogeneous nucleation and growth of melt in copper

Molecular dynamics simulations are conducted to investigate homogeneous nucleation and growth of melt in copper described by an embedded-atom method (EAM) potential. The accuracy of this EAM potential for melting is validated by the equilibrium melting point obtained with the solid-liquid coexistence method and the superheating-supercooling hysteresis method. We characterize the atomistic melting process by following the temperature and time evolution of liquid atoms. The nucleation behavior at the extreme superheating is analyzed with the mean-first-passage-time (MFPT) method, which yields the critical size, steady-state nucleation rate, and the Zeldovich factor. The value of the steady-state nucleation rate obtained from the MFPT method is consistent with the result from direct simulations. The size distribution of subcritical nuclei appears to follow a power law similar to three-dimensional percolation. The diffuse solid-liquid interface has a sigmoidal profile with a 10%-90% width of about 12 A near the critical nucleation. The critical size obtained from our simulations is in reasonable agreement with the prediction of classical nucleation theory if the finite interface width is considered. The growth of melt is coupled with nucleation and can be described qualitatively with the Johnson-Meh-Avrami law. System sizes of 10(3)-10(6) atoms are explored, and negligible size dependence is found for bulk properties and for the critical nucleation.

Molecular dynamics simulations of melting and the glass transition of nitromethane

Molecular dynamics simulations have been used to investigate the thermodynamic melting point of the crystalline nitromethane, the melting mechanism of superheated crystalline nitromethane, and the physical properties of crystalline and glassy nitromethane. The maximum superheating and glass transition temperatures of nitromethane are calculated to be 316 and 160 K, respectively, for heating and cooling rates of 8.9 x 10(9) Ks. Using the hysteresis method [Luo et al., J. Chem. Phys. 120, 11640 (2004)] and by taking the glass transition temperature as the supercooling temperature, we calculate a value of 251.1 K for the thermodynamic melting point, which is in excellent agreement with the two-phase result [Agrawal et al., J. Chem. Phys. 119, 9617 (2003)] of 255.5 K and measured value of 244.73 K. In the melting process, the nitromethane molecules begin to rotate about their lattice positions in the crystal, followed by translational freedom of the molecules. A nucleation mechanism for the melting is illustrated by the distribution of the local translational order parameter. The critical values of the Lindemann index for the C and N atoms immediately prior to melting (the Lindemann criterion) are found to be around 0.155 at 1 atm. The intramolecular motions and molecular structure of nitromethane undergo no abrupt changes upon melting, indicating that the intramolecular degrees of freedom have little effect on the melting. The thermal expansion coefficient and bulk modulus are predicted to be about two or three times larger in crystalline nitromethane than in glassy nitromethane. The vibrational density of states is almost identical in both phases.

Practically Efficient QM/MM Alchemical Free Energy Simulations: The Orthogonal Space Random Walk Strategy

The difference between free energy changes occurring at two chemical states can be rigorously estimated via alchemical free energy (AFE) simulations. Traditionally, most AFE simulations are carried out under the classical energy potential treatment; then, accuracy and applicability of AFE simulations are limited. In the present work, we integrate a recent second-order generalized ensemble strategy, the orthogonal space random walk (OSRW) method, into the combined quantum mechanical/molecular mechanical (QM/MM) potential based AFE simulation scheme. Thereby, within a commonly affordable simulation length, accurate QM/MM alchemical free energy simulations can be achieved. As revealed by the model study on the equilibrium of a tautomerization process of hydrated 3-hydroxypyrazole and by the model calculations of the redox potentials of two flavin derivatives, lumichrome (LC) and riboflavin (RF) in aqueous solution, the present OSRW-based scheme could be a viable path toward the realization of practically efficient QM/MM AFE simulations.

The relation between shock-state particle velocity and free surface velocity : A molecular dynamics study on single crystal Cu and silica glass

We investigate the ratio Rrp of the free surface velocity to the shock-state particle velocity during shock wave loading with molecular dynamics simulations on two representative solids, single crystal Cu, and silica glass. The free surface velocity is obtained as a function of the particle velocity behind the shock front (or shock stress) for loading on Cu along ⟨100⟩, ⟨110⟩, and ⟨111⟩, and on the isotropic glass. Rrp≥1 for Cu and Rrp<1 for silica glass, and it increases with shock strength; the simulations agree well with the experimental results. For supported shock loading of silica glass at 30–90 GPa, the SiIV–SiVI transition occurs upon shock, inducing substantial densification and thus small Rrp (0.65–0.78). For single crystal Cu, Rrp deviates from 1 near the Hugoniot elastic limit and reaches ∼1.2 at 355 GPa for ⟨100⟩ shock. Rrp is anisotropic, e.g., it is about 1.02, 1.08, and 1.06 for shock loading to about 80 GPa along ⟨100⟩, ⟨110⟩, and ⟨111⟩, respectively. Such an anisotropy is mostly due to t...

Molecular dynamics simulations of melting of perfect crystalline hexahydro-1,3,5-trinitro-1,3,5-s-triazine

The melting mechanism of superheated perfect crystalline hexahydro-1,3,5-trinitro-1,3,5-s-triazine (alpha-RDX) has been investigated using molecular dynamics simulations with the fully flexible force field developed by Smith and Bharadwaj [J. Phys. Chem. B 103, 3570 (1999)]. Sequential 50 ps equilibration simulations of the constant stress-constant temperature ensemble were performed at 10 K intervals over the range of 300-650 K, corresponding to a heating rate of 2.0 x 10(11) Ks. A solid-solid phase transition is observed between 480 and 490 K, followed by melting, which occurs between 500 and 510 K. The solid-solid phase transition, both displacive and rotational, is characterized by an abrupt decrease in the lengths of the unit cell edges a and b and an increase of the length of edge c. The molecular conformation in the new phase is AAE, although the axial nitro groups have different changes: one shift is more axial and the other is more equatorial. Phases other than alpha-RDX have been observed experimentally, however, there are insufficient data for comparisons to ascertain that the new phase observed here corresponds to a real phase. At the high heating rate (2.0 x 10(11) Ks) used in the simulations, the melted RDX reaches full orientational disorder at about 540 K and translational freedom at around 580 K. If the simulation at the melting temperature (510 K) is run sufficiently long complete rotational freedom is achieved in a few hundreds of picoseconds, while complete translational freedom requires much longer. These results show that given a sufficiently high heating rate, the system can exist for significant periods of time in a near-liquid state in which the molecules are not as free to rotate and diffuse as in the true liquid state. The bond lengths and bond angles undergo little change upon melting, while there are significant changes in the dihedral angles. The molecular conformation of RDX changes from AAE to EEE upon melting. The ramification of this for formulating force fields that accurately describe melting is that it is important that the torsional motions are accurately described.

Enhancing QM/MM molecular dynamics sampling in explicit environments via an orthogonal-space-random-walk-based strategy.

Accurate prediction of molecular conformations in explicit environments, such as aqueous solution and protein interiors, can facilitate our understanding of various molecular recognition processes. Most computational approaches are limited as a result of their compromised choices between the underlying energy model and the sampling length. Taking advantage of a recent second-order generalized ensemble scheme [e.g., the orthogonal space random walk (OSRW) strategy], which can synergistically accelerate the motion of a focused region and its coupled environmental response, we are presenting a QM/MM (combined quantum mechanical/molecular mechanical)-based molecular dynamics sampling technique to explore molecular conformational landscapes in explicit environments. The present QM/MM potential scaling-based OSRW sampling scheme is employed to study the binding of DMSO to the FKBP12 protein, the conformation distribution of a novel mercaptosulfonamide inhibitor in aqueous solution, and its binding poses in zinc-containing matrix metalloproteinase-9 (MMP-9). As demonstrated, the present QM/MM second-order generalized ensemble sampling technique enables feasible usage of the QM/MM model to sample molecular conformations in condensed environments.

Essential energy space random walks to accelerate molecular dynamics simulations: Convergence improvements via an adaptive-length self-healing strategy

Recently, accelerated molecular dynamics (AMD) technique was generalized to realize essential energy space random walks so that further sampling enhancement and effective localized enhanced sampling could be achieved. This method is especially meaningful when essential coordinates of the target events are not priori known; moreover, the energy space metadynamics method was also introduced so that biasing free energy functions can be robustly generated. Despite the promising features of this method, due to the nonequilibrium nature of the metadynamics recursion, it is challenging to rigorously use the data obtained at the recursion stage to perform equilibrium analysis, such as free energy surface mapping; therefore, a large amount of data ought to be wasted. To resolve such problem so as to further improve simulation convergence, as promised in our original paper, we are reporting an alternate approach: the adaptive-length self-healing (ALSH) strategy for AMD simulations; this development is based on a recent self-healing umbrella sampling method. Here, the unit simulation length for each self-healing recursion is increasingly updated based on the Wang-Landau flattening judgment. When the unit simulation length for each update is long enough, all the following unit simulations naturally run into the equilibrium regime. Thereafter, these unit simulations can serve for the dual purposes of recursion and equilibrium analysis. As demonstrated in our model studies, by applying ALSH, both fast recursion and short nonequilibrium data waste can be compromised. As a result, combining all the data obtained from all the unit simulations that are in the equilibrium regime via the weighted histogram analysis method, efficient convergence can be robustly ensured, especially for the purpose of free energy surface mapping.