Search

Enhanced sampling techniques that target a set of collective variables and that use molecular dynamics as the driving engine have seen widespread application in the computational molecular sciences as a means to explore the free-energy landscapes of complex systems. The use of molecular dynamics as the fundamental driver of the sampling requires the introduction of a time step whose magnitude is limited by the fastest motions in a system. While standard multiple time-stepping methods allow larger time steps to be employed for the slower and computationally more expensive forces, the maximum achievable increase in time step is limited by resonance phenomena, which inextricably couple fast and slow motions. Recently, we introduced deterministic and stochastic resonance-free multiple time step algorithms for molecular dynamics that solve this resonance problem and allow ten- to twenty-fold gains in the large time step compared to standard multiple time step algorithms [P. Minary et al., Phys. Rev. Lett. 93, 150201 (2004); B. Leimkuhler et al., Mol. Phys. 111, 3579-3594 (2013)]. These methods are based on the imposition of isokinetic constraints that couple the physical system to Nosé-Hoover chains or Nos'e-Hoover Langevin schemes. In this paper, we show how to adapt these methods for collective variable-based enhanced sampling techniques, specifically adiabatic free-energy dynamics/temperature-accelerated molecular dynamics, unified free-energy dynamics, and by extension, metadynamics, thus allowing simulations employing these methods to employ similarly very large time steps. The combination of resonance-free multiple time step integrators with free-energy-based enhanced sampling significantly improves the efficiency of conformational exploration. [ABSTRACT FROM AUTHOR]

Extended phase-space isokinetic methods in their deterministic [Minary et al., Phys. Rev. Lett. 93, 150201 (2004)] and stochastic forms [Leimkuhler et al., Mol. Phys. 111, 3579 (2013)] have proved tremendously successful in allowing multiple time-scale molecular dynamics simulations to be performed with very large time steps. These methods work by coupling the physical degrees of freedom to a set of Nosé-Hoover chain or Nosé-Hoover Langevin thermostats via an isokinetic constraint, which has the effect of avoiding resonance artifacts that plague multiple time-step algorithms. In this paper, we introduce a new resonance-free approach that achieves the same gains in time step but without the imposition of isokinetic constraints or the introduction of extended phase-space variables. Rather, we modify the physical Hamiltonian that effects the same regulation of resonances achieved by the isokinetic constraints. In so doing, we show that sampling errors can be controlled and performance improvements are possible within a simpler Hamiltonian framework. The method is demonstrated in simulations of the structure of liquid water and, in conjunction with enhanced sampling, in generation of the Ramachandran free-energy surface of the solvated alanine dipeptide. [ABSTRACT FROM AUTHOR]

Higher order factorization schemes are developed for path integral molecular dynamics in order to improve the convergence of estimators for physical observables as a function of the Trotter number. The methods are based on the Takahashi-Imada and Susuki decompositions of the Boltzmann operator. The methods introduced improve the averages of the estimators by using the classical forces needed to carry out the dynamics to construct a posteriori weighting factors for standard path integral molecular dynamics. The new approaches are straightforward to implement in existing path integral codes and carry no significant overhead. The Suzuki higher order factorization was also used to improve the end-to-end distance estimator in open path integral molecular dynamics. The new schemes are tested in various model systems, including an ab initio path integral molecular dynamics calculation on the hydrogen molecule and a quantum water model. The proposed algorithms have potential utility for reducing the cost of path integral molecular dynamics calculations of bulk systems. [ABSTRACT FROM AUTHOR]

Extraction of relaxation times, lifetimes, and rates associated with the transport of topological charge defects in hydrogen-bonded networks from molecular dynamics simulations is a challenge because proton transfer reactions continually change the identity of the defect core. In this paper, we present a statistical mechanical theory that allows these quantities to be computed in an unbiased manner. The theory employs a set of suitably defined indicator or population functions for locating a defect structure and their associated correlation functions. These functions are then used to develop a chemical master equation framework from which the rates and lifetimes can be determined. Furthermore, we develop an integral equation formalism for connecting various types of population correlation functions and derive an iterative solution to the equation, which is given a graphical interpretation. The chemical master equation framework is applied to the problems of both hydronium and hydroxide transport in bulk water. For each case it is shown that the theory establishes direct links between the defect's dominant solvation structures, the kinetics of charge transfer, and the mechanism of structural diffusion. A detailed analysis is presented for aqueous hydroxide, examining both reorientational time scales and relaxation of the rotational anisotropy, which is correlated with recent experimental results for these quantities. Finally, for OH-(aq) it is demonstrated that the 'dynamical hypercoordination mechanism' is consistent with available experimental data while other mechanistic proposals are shown to fail. As a means of going beyond the linear rate theory valid from short up to intermediate time scales, a fractional kinetic model is introduced in the Appendix in order to describe the nonexponential long-time behavior of time-correlation functions. Within the mathematical framework of fractional calculus the power law decay ∼t-σ, where σ is a parameter of the model and depends on the dimensionality of the system, is obtained from Mittag-Leffler functions due to their long-time asymptotics, whereas (stretched) exponential behavior is found for short times. [ABSTRACT FROM AUTHOR]

The problems of ergodicity and internal consistency in the centroid and ring-polymer molecular dynamics methods are addressed in the context of a comparative study of the two methods. Enhanced sampling in ring-polymer molecular dynamics (RPMD) is achieved by first performing an equilibrium path integral calculation and then launching RPMD trajectories from selected, stochastically independent equilibrium configurations. It is shown that this approach converges more rapidly than periodic resampling of velocities from a single long RPMD run. Dynamical quantities obtained from RPMD and centroid molecular dynamics (CMD) are compared to exact results for a variety of model systems. Fully converged results for correlations functions are presented for several one dimensional systems and para-hydrogen near its triple point using an improved sampling technique. Our results indicate that CMD shows very similar performance to RPMD. The quality of each method is further assessed via a new χ2 descriptor constructed by transforming approximate real-time correlation functions from CMD and RPMD trajectories to imaginary time and comparing these to numerically exact imaginary time correlation functions. For para-hydrogen near its triple point, it is found that adiabatic CMD and RPMD both have similar χ2 error. [ABSTRACT FROM AUTHOR]

An efficient computational approach to perform Car–Parrinello ab initio molecular dynamics (CPAIMD) simulations under cluster (free) boundary conditions is presented. The general approach builds upon a recent real-space CPAIMD formalism using discrete variable representation (DVR) basis sets [Y. Liu et al., Phys. Rev. B 12, 125110 (2003); H.-S. Lee and M. E. Tuckerman, J. Phys. Chem. A 110, 5549 (2006)]. In order to satisfy cluster boundary conditions, a DVR based on sinc functions is utilized to expand the Kohn–Sham orbitals and electron density. Poisson’s equation is solved in order to calculate the Hartree potential via an integral representation of the 1/r singularity. Excellent convergence properties are achieved with respect to the number of grid points (or DVR functions) and the size of the simulation cell. A straightforward implementation of the present approach leads to near linear scaling [O(N4/3)] of the computational cost with respect to the system size (N) for the solution of Poisson’s equation. The accuracy and stability of CPAIMD simulations based on sinc DVR are tested for a model problem as well as for N2 and a water dimer. [ABSTRACT FROM AUTHOR]

Dynamical properties of liquid water were studied using Car-Parrinello [Phys. Rev. Lett. 55, 2471 (1985)] ab initio molecular dynamics (AIMD) simulations within the Kohn-Sham (KS) density functional theory employing the Becke-Lee-Yang-Parr exchange-correlation functional for the electronic structure. The KS orbitals were expanded in a discrete variable representation basis set, wherein the complete basis set limit can be easily reached and which, therefore, provides complete convergence of ionic forces. In order to minimize possible nonergodic behavior of the simulated water system in a constant energy (NVE) ensemble, a long equilibration run (30 ps) preceded a 60 ps long production run. The temperature drift during the entire 60 ps trajectory was found to be minimal. The diffusion coefficient [0.055 Å2/ps] obtained from the present work for 32 D2O molecules is a factor of 4 smaller than the most up to date experimental value, but significantly larger than those of other recent AIMD studies. Adjusting the experimental result so as to match the finite-sized system used in the present study brings the comparison between theory and experiment to within a factor of 3. More importantly, the system is not observed to become “glassy” as has been reported in previous AIMD studies. The computed infrared spectrum is in good agreement with experimental data, especially in the low frequency regime where the translational and librational motions of water are manifested. The long simulation length also made it possible to perform detailed studies of hydrogen bond dynamics. The relaxation dynamics of hydrogen bonds observed in the present AIMD simulation is slower than those of popular force fields, such as the TIP4P potential, but comparable to that of the TIP5P potential. [ABSTRACT FROM AUTHOR]

We present a rigorous description of chemical space within a molecular grand-canonical ensemble multi-component density functional theory framework. A total energy density functional for chemical compounds in contact with an electron and a proton bath is introduced using Lagrange multipliers which correspond to the energetic response to changes of the elementary particle densities. From a generalized Gibbs-Duhem equation analog, reactivity indices such as the nuclear hardness and a molecular Fukui function, which couples the grand-canonical electronic and nuclear degrees of freedom, are obtained. Maxwell relations between composition particles, ionic displacements, and the external potential are discussed. Numerical results for the molecular Fukui function are presented as well as finite temperature estimates for the oxidation of ammonia. [ABSTRACT FROM AUTHOR]

Structural properties of liquid water at ambient temperature were studied using Car-Parrinello [Phys. Rev. Lett. 55, 2471 (1985)] ab initio molecular dynamics (CPAIMD) simulations combined with the Kohn-Sham (KS) density functional theory and the BLYP exchange-correlation functional for the electronic structure. Unlike other recent work on the same subject, where plane-wave (PW) or hybrid Gaussian/plane-wave basis sets were employed, in the present paper, a discrete variable representation (DVR) basis set is used to expand the KS orbitals, so that with the real-space grid adapted in the present work, the properties of liquid water could be obtained very near the complete basis set limit. Structural properties of liquid water were extracted from a 30 ps CPAIMD-BLYP/DVR trajectory at 300 K. The radial distribution functions (RDFs), spatial distribution functions, and hydrogen bond geometry obtained from the CPAIMD-BLYP/DVR simulation are generally in good agreement with the most up to date experimental measurements. Compared to recent ab initio MD simulations based on PW basis sets, less significant overstructuring was found in the RDFs and the distributions of hydrogen bond angles, suggesting that previous plane-wave and Gaussian basis set calculations have exaggerated the tendency toward overstructuring. [ABSTRACT FROM AUTHOR]

A multistate empirical valence bond model for proton transport in water, which explicitly includes solvent polarization, is presented. Polarization is included for each valence-bond state via induced point dipoles, and the model is parametrized to be used with an effective path integral derived potential surface, so as to include quantum effects of the transferring proton. The new model is shown to reproduce ab initio geometries and energetics for small protonated clusters. It is also shown that the new model gives a diffusion constant for the excess proton in water, which is in good agreement with experiment, and that the qualitative features of ab initio path integral simulations [D. Marx, M. E. Tuckerman, J. Hutter, and M. Parrinello, Nature (London) 397, 601 (1999)] are well reproduced. [ABSTRACT FROM AUTHOR]

The theoretical principles underpinning the calculation of infrared spectra for condensed-phase systems in the context of ab initio molecular dynamics have been recently developed in literature. At present, most ab initio molecular dynamics calculations are restricted to relatively small systems and short simulation times. In this paper we devise a method that allows well-converged results for infrared spectra from ab initio molecular dynamics simulations using small systems and short trajectories characteristic of simulations typically performed in practice. We demonstrate the utility of our approach by computing the imaginary part of the dielectric constant [variant_greek_epsilon]″(ω) for H2O and D2O in solid and liquid phases and show that it compares well with experimental data. We further demonstrate that maximally localized Wannier orbitals can be used to separate the individual contributions of different molecular species to the linear spectrum of complex systems. The new spectral decomposition method is shown to be useful in present-day ab initio molecular dynamics calculations to compute the magnitude of the “continuous absorption” generated by excess protons in aqueous solutions with good accuracy even when other species present in the solutions absorb strongly in the same frequency window. [ABSTRACT FROM AUTHOR]

The transport of protons through aqueous, partially aqueous, or nonaqueous hydrogen-bonded media is a fundamental process in many biologically and technologically important systems. Liquid methanol is an example of a hydrogen-bonded system that, like water, supports anomalously fast proton transport. Using the methodology of ab initio molecular dynamics, in which internuclear forces are computed directly from electronic structure calculations as the simulation proceeds, we have investigated the microscopic mechanism of the proton transport process in liquid methanol at 300 K. It is found that the defect structure associated with an excess proton in liquid methanol is a hydrogen-bonded cationic chain whose length generally exceeds the average chain length in pure liquid methanol. Hydrogen bonds in the first and second solvation shells of the excess proton are considerably shorter and stronger than ordinary methanol-methanol hydrogen bonds. Along this chain, proton transfer reactions occur in an essentially random manner described by Poisson statistics. Structural diffusion of the defect structure is possible if the proton migrates toward an end of the defect chain, which causes a weakening of the hydrogen bonds at the opposite end. The latter can, therefore, be easily ruptured by ordinary thermal fluctuations. At the end of the chain where the proton resides, new hydrogen bonds are likely to form due to the strong associative nature of the excess proton. It is through this "snake-like" mechanism that the defect structure is able to diffuse through the hydrogen-bond network of the liquid. The estimated activation enthalpy of this proposed mechanism is found to be in reasonable agreement with the experimentally determined activation enthalpy. [ABSTRACT FROM AUTHOR]

A new formalism designed to treat long range interactions on surfaces, systems which are infinitely replicated in two spatial directions but have finite extent in the third, is developed. The new formalism is based in reciprocal space and, thus, permits the facile extension of standard plane-wave based density functional theory, Ewald summation, and smooth particle-mesh Ewald methods to handle surfaces efficiently. The method is tested on both model (body centered cubic lattices) and realistic problems (an ice surface with a defect and the 2 x 1 surface reconstruction of silicon) and found to be accurate, efficient, and a marked improvement on existing formulations in speed, accuracy, and utility. [ABSTRACT FROM AUTHOR]

Quantum time correlation functions (TCFs) involving two states are important for describing nonadiabatic dynamical processes such as charge transfer (CT). Based on a previous single-state method, we propose an imaginary-time open-chain path-integral (OCPI) approach for evaluating the two-state symmetrized TCFs. Expressing the forward and backward propagation on different electronic potential energy surfaces as a complex-time path integral, we then transform the path variables to average and difference variables such that the integration over the difference variables up to the second order can be performed analytically. The resulting expression for the symmetrized TCF is equivalent to sampling the open-chain configurations in an effective potential that corresponds to the average surface. Using importance sampling over the extended OCPI space via open path-integral molecular dynamics, we tested the resulting path-integral approximation by calculating the Fermi's golden rule CT rate constant within a widely used spin-boson model. Comparing with the real-time linearized semiclassical method and analytical result, we show that the imaginary-time OCPI provides an accurate two-state symmetrized TCF and rate constant in the typical turnover region. It is shown that the first bead of the open chain corresponds to physical zero-time and that the endpoint bead corresponds to final time t; oscillations of the end-to-end distance perfectly match the nuclear mode frequency. The two-state OCPI scheme is seen to capture the tested model's electronic quantum coherence and nuclear quantum effects accurately. [ABSTRACT FROM AUTHOR]

Presents a reciprocal space based formalism for treating long range forces in clusters. Incorporation of the formalism into plane-wave based density function theory calculations, Ewald summation calculations and smooth particle-mesh Ewald calculations to yield accurate and numerically efficient descriptions of long range interactions in cluster systems.

Investigates the double proton transfer in the formic acid dimer using ab initio path integral molecular dynamics (PIMD). Description of the electronic structure by Kohn-Sham formulation of the density functional theory; Details on the PIMD approach of Tuckerman; Results of the thermal and quantum fluctuations around the equilibrium structure.

Examines the empirical valence model for proton transfer in aqueous solutions. Properties of the gas-phase H5O2+ complex; Use of gradient-corrected density functional theory in the study; Simulations of the solvated hydronium.

Focuses on the nonequilibrium molecular dynamics generated by the SLLOD algorithm for fluid flow. Existence of a conserved energy for the equations of motion; View that coupling each degree of freedom in the system to a Nose-Hoover chain thermostat as a means of efficiently generating the phase space equilibrium.

In path integral molecular dynamics, efficient sampling of the phase space is not guaranteed due to the stiff harmonic part of the action arising from the quantum kinetic energy. This problem has been eliminated by incorporating a sufficient number of thermostats into the dynamical scheme and by introducing a transformation of the path ‘‘bead’’ variables. In this paper, an efficient Car–Parrinello path integral molecular dynamics algorithm, sufficiently general to include the use of ultrasoft pseudopotentials is introduced. Difficulties encountered when combining thermostats and transformations of the Cartesian ‘‘bead’’ coordinates with the generalized orthonormality condition are circumvented by employing a constrained nonorthogonal orbital method. © 1996 American Institute of Physics. [ABSTRACT FROM AUTHOR]

A new fourth order predictor–corrector integration scheme is presented. The unique feature of the new algorithm and what distinguishes it from a Gear predictor–corrector is that the method is derived from the Trotter decomposition of a specially formulated evolution operator and as such, is both symplectic and reversible. In addition, the method retains the useful property of Gear methods that only one force evaluation per time step is required. The new integrator is tested on a harmonic plus quartic oscillator and the Henon–Heiles system. Comparisons are made to the second order velocity Verlet integrator, the true fourth order Yoshida/Suzuki schemes and fourth order Gear. In all cases, the new method works well, giving energy conservation and trajectories of much better quality than velocity Verlet and of comparable quality to the results of the true fourth order schemes for the same computational cost as velocity Verlet. © 1995 American Institute of Physics. [ABSTRACT FROM AUTHOR]

In this paper, new techniques for integrating the Car–Parrinello equations based on multiple time scale methodology are presented. The formalism of multiple time scale methodology based on operator factorizations of the classical Liouville propagator is reviewed. It is shown how the techniques are applied to Car–Parrinello for use with the velocity Verlet and Gaussian dynamics schemes presented in the preceding paper [M. E. Tuckerman and M. Parrinello, J. Chem. Phys. 101, 1302 (1994)], and a detailed discussion is presented of how a reference system for Car–Parrinello simulations may be chosen. It is shown that the use of such techniques can save up to a factor of 5–10 in cpu time over the standard methods. [ABSTRACT FROM AUTHOR]

In this paper and in a companion paper [M. E. Tuckerman and M. Parrinello, J. Chem. Phys. 101, 1316 (1994)] the problem of integrating the equations of motion in Car–Parrinello simulations is addressed. In this paper, new techniques for treating the constraint problem based on the velocity Verlet integrator and the Gaussian dynamics are presented. Questions of adiabaticity and temperature control are discussed, and it is shown how to combine the new techniques with the recently developed Nosé–Hoover chain thermostat method. All new techniques are described using the formalism of operator factorizations applied to the classical Liouville propagator. In the companion paper, the formalism and application of multiple time scale methodology in Car–Parrinello simulations are discussed. [ABSTRACT FROM AUTHOR]

The reference system propagator algorithm of Tuckerman et al. is applied to a simple model of flexible CO2. The forces are broken up according to three different time scales corresponding to the internal motion, short-range impulsive collisions, and long-range interactions, and three time steps are introduced to integrate the equations of motion. The use of this scheme allows a factor of 8 to 12 to be saved in CPU time over standard integration schemes. [ABSTRACT FROM AUTHOR]

New path integral molecular dynamics (PIMD) and path integral hybrid Monte Carlo (PIHMC) algorithms are developed. It is shown that the use of a simple noncanonical change of variables that naturally divides the quadratic part of the action into long and short wavelength modes and multiple time scale integration techniques results in very efficient algorithms. The PIMD method also employs a constant temperature MD technique that has been shown to give canonical averages even for stiff systems. The new methods are applied to the simple quantum mechanical harmonic oscillator and to electron solvation in fluid helium and xenon. Comparisons are made with PIMC and the more basic PIMD and PIHMC methods. [ABSTRACT FROM AUTHOR]

Rate constants evaluated from (1) the energy-loss turnover theory of Pollak, Grabert, and Hänggi (PGH), (2) the Grote–Hynes extension of Kramers theory (GH), and (3) the microcanonical variational transition state theory for dissipative systems of Tucker and Pollak (μVTST) are compared with rate constants determined from direct computer simulations of generalized Langevin dynamics. The comparisons are made for a cubic oscillator under the influence of a slow bath characterized by a Gaussian friction kernel. In the μVTST calculations, which are based on an effective two degree of freedom Hamiltonian, barrier crossing due to energy transfer from the bath to the effective Hamiltonian is neglected. This neglect is significant only at very strong coupling, where it causes the μVTST results to drop below the simulation results. Both GH and μVTST theories fail (as expected) in the energy diffusion regime, while PGH theory is only moderately successful. The μVTST results agree extremely well with the simulation results in the spatial diffusion regime, providing a significant improvement over the GH results at intermediate coupling strengths and over the PGH results at strong coupling strengths. This improvement is a result of nonlinear effects which are included in the μVTST approach but neglected in the PGH and GH theories. [ABSTRACT FROM AUTHOR]

The generalized Langevin equation (GLE) has been used to model a wide variety of systems in which a subset of the degrees of freedom move on a potential of mean force surface subject to fluctuating forces and dynamic friction. When there is a wide separation in the time scales for motion on the potential surface and for relaxation of the friction kernel, direct integration of the GLE is very costly in CPU time. In this paper we introduce an integrator based on our previous work using numerical analytical propagator algorithm (NAPA) and reference system propagator algorithm (RESPA) that greatly accelerates such simulations. We also discuss sampling methods for the random force. Accuracy of this algorithm is assessed by comparisons with an analytically solvable example. Introducing dynamic friction kernels determined from full molecular dynamics (MD) simulations allows us to compare the accuracy of the GLE simulations with full scale molecular dynamics simulations. [ABSTRACT FROM AUTHOR]

A frequently encountered problem in molecular dynamics is how to treat the long times that are required to simulate condensed systems consisting of particles interacting through long range forces. Standard methods require the calculation of the forces at every time step. Because each particle interacts with all particles within the interaction range of the potential the longer the range of the potential the larger the number forces that must be calculated at each time step. In this note we present a variant of the RESPA (reference system propagator algorithm), which we developed for handling systems with multiple time scales like disparate mass mixtures. This version of RESPA greatly reduces the number of forces that must be computed at each time step and thereby leads to a dramatic acceleration of such simulations. The RESPA method uses ideas similar to NAPA, an algorithm we invented to treat high frequency oscillators interacting with low frequency bath. The method is based on a choice of a reference system in which the particles interact through short range forces. The reference system is numerically integrated for n time steps δt and the error incurred by using short range forces is corrected by solving a rigorous set of equations once every Δt=nδt. This method reduces the cpu time dramatically. It is shown that this approach and suitable generalizations should be very useful for future simulations of quantum and classical condensed matter systems. [ABSTRACT FROM AUTHOR]

A frequently encountered problem in molecular dynamics is how to treat the long times that are required to simulate condensed systems consisting of mixtures of light and heavy particles. Standard methods require the choice of time step sufficiently small to guarantee stable solution for the low mass component with the consequence that these simulations require a very large number of central processing unit cycles to treat the relaxation of the heavier component. In this note, we present a new method that allows one to use a time step appropriate for the heavy particles. This method uses a similar idea to numerical analytical propogator algorithm, an algorithm we invented to treat high frequency oscillators interacting with low frequency baths and is based on a choice of a reference system for the light particle motions. The method is applied to the case of a liquid containing 864 Lennard-Jones spheres, 824 of these particles having a mass, M=100 and 40 spheres picked at random have a mass m=1. It is shown that molecular dynamics using the new algorithm runs seven to ten times faster than standard methods and this approach as well as suitable generalizations should be very useful for future simulations of quantum and classical condensed matter systems. [ABSTRACT FROM AUTHOR]

A frequently encountered problem in molecular dynamics simulations is the long runs required to study condensed systems consisting of both high frequency and low frequency degrees of freedom. Standard integrators require the choice of time step sufficiently small to guarantee stable solution of the highest frequency motion with the consequence that simulations require a very large number of central processing unit (CPU) cycles. In this note we present a new integrator that allows one to use a time step appropriate for the low frequency degrees of freedom without making any approximations related to the separation of time scales. This method is based on a choice of an analytically solvable reference system for the high frequency motion. We show how the analytical solution can be incorporated into a numerical integrator. The method is applied to two cases which are paradigms for this problem and it is shown that this approach and suitable generalizations should be very useful for future simulations of quantum and classical condensed matter systems. [ABSTRACT FROM AUTHOR]

Comments on the scientific paper 'Simple Reversible Molecular Dynamics Algorithms for Nose-Hoover Chain Dynamics,' by S. Jang and J.A. Voth. Criticism on the proposed leap-frog algorithms of Jang and Voth; Failure to consider energy conservation in the Nose-Hoover chain dynamics.

The hydrated morphology of either proton exchange membranes (PEMs) or anion exchange membranes (AEMs) determines many aspects of species transport. The present work seeks to understand the morphology and microstructure of a triblock copolymer, polystyrene-b-poly(ethylene-co-butylene)-b-polystyrene (SEBS), functionalized with alkyl-substituted quaternary ammonium groups. Mesoscale dissipative particle dynamics (DPD) simulations were utilized and parametrized by reproducing the experimental morphology of the SEBS copolymer. It was found that the AEM (i.e., quaternary ammonium-functionalized SEBS) phase separates into a functionalized polystyrene-rich phase that is hydrophilic and a hydrophobic phase consisting of the SEBS mid-blocks. The morphology was controlled by the water content and was transformed from perforated and interconnected lamellae to perfect lamellae and then to disordered bicontinuous domains by increasing the hydration level (λ = H2O/functional head group) from 4 to 20. The hydrophilic phase swelled upon the hydration of the membrane consistent with AFM phase imaging of a similar SEBS-based ionomer. Domains exclusively consisting of water were formed at high levels of hydration (λ = 16 and 20) within the hydrophilic phase. Changing the anion from OH– to Cl– resulted in larger water domains at the highest hydration levels. [ABSTRACT FROM AUTHOR]

When there are high and low frequency motions in systems with long and short range forces a judicious choice of reference system leads to very large accelerations in molecular dynamic (MD) simulations. Building on our previous work where we have developed reference system methods for systems with high frequency oscillators, disparate masses, or long range forces, we present a double reference system method which leads to acceleration of as much as 20 in systems consisting of 864 molecules with Lennard-Jones (12-6) forces. Much larger savings should be achieved when this method is applied to longer range forces and larger systems. [ABSTRACT FROM AUTHOR]

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. [ABSTRACT FROM AUTHOR]

Abstract: Phosphoric acid based polymer electrolyte membranes (e.g. poly(benzimidazole)/phosphoric acid) are the only quasi-anhydrous separator electrolytes for medium temperature fuel cells currently in use. Pure phosphoric acid has the highest known intrinsic proton conductivity of any substance, which is very sensitive towards doping with Brønsted–Lowry bases or even acids. Proton conductivity shows the most dramatic decrease upon addition of basic species, e.g. polymers containing N-heterocyclic groups, and high phosphoric acid contents, are required to ensure high conductivities in these materials. We perform first principles molecular dynamics simulations in order to investigate these processes at the fundamental molecular level. Proton dynamics and transport are studied in a system of phosphoric acid and imidazole having 2:1 molar ratio. This model system represents a first step towards the more detailed understanding of the proton conduction mechanisms in realistic phosphoric acid based polymer electrolyte materials. Our results indicate that imidazole molecules have a much stronger proton affinity (basicity) than phosphoric acid, forming very stable imidazolium–dihydrogenphosphate complexes. We evaluate a number of properties which are important for the proton conduction and compare them to the previously obtained data on pure phosphoric acid. The reduced proton density in the proton conducting H3PO4 phase is shown to have strong effects on the energetics of proton transfer, local proton dynamics, electrostatic proton–proton interactions as well as on the structure of the hydrogen bond network. These effects are suggested to be the reason for the impaired proton conducting properties of this, and similar, phosphoric acid based systems containing N-heterocycles. [Copyright &y& Elsevier]

Neat liquid phosphoric acid (H3PO4) has the highest intrinsic proton conductivity of any known substance and is a useful model for understanding proton transport in other phosphate-based systems in biology and clean energy technologies. Here, we present an ab initio molecular dynamics study that reveals, for the first time, the microscopic mechanism of this high proton conductivity. Anomalously fast proton transport in hydrogen-bonded systems involves a structural diffusion mechanism in which intramolecular proton transfer is driven by specific hydrogen bond rearrangements in the surrounding environment. Aqueous media transport excess charge defects through local hydrogen bond rearrangements that drive individual proton transfer reactions. In contrast, strong, polarizable hydrogen bonds in phosphoric acid produce coupled proton motion and a pronounced protic dielectric response of the medium, leading to the formation of extended, polarized hydrogen-bonded chains. The interplay between these chains and a frustrated hydrogen-bond network gives rise to the high proton conductivity. [ABSTRACT FROM AUTHOR]

Abstract: The properties of the hydrated hydroxide ion (aq) stand in sharp contrast to those of most other aqueous ions. Chief among these is its anomalously high mobility, which is shared only by the aqueous hydronium ion (aq). However, while the transport mechanism of (aq) is now well understood, the details of (aq) diffusion at ambient conditions are just beginning to emerge, and the effects of temperature on the transport mechanism remain largely unelucidated. Here, we undertake an ab initio molecular dynamics study of the effect of temperature on the solvation and transport of (aq). In addition to revealing new details of the transport process, our analysis provides an explanation for the experimentally observed temperature dependence of the (aq) reorientation time. The calculations reveal a suppression of proton transfer events that underly the structural diffusion process caused by a pronounced change in the population of dominant (aq) solvation complexes. [Copyright &y& Elsevier]

When only equilibrium thermodynamics or averages of position-dependent functions in classical molecular dynamics calculations are sought, the choice of the atomic masses in the Hamiltonian becomes irrelevant. Consequently, the masses can be used as free parameters to help improve conformational sampling efficiency in the spirit of Bennett's mass-tensor dynamics (Bennett C. H. J. Comput. Phys. 1975 , 19 , 267.). Here we report on the effect of reducing side-chain and solvent masses on the folding behavior of a 9-residue hairpin-forming peptide in solution. A physical motivation for reducing the solvent mass is an effective reduction in solvent viscosity. Side-chain mass scaling is motivated by the idea of solvent potentials of mean force which are employed in select coarse-grained protein models. In the limit of very large mass differences, the mass reduction effectively creates an adiabatic decoupling between solvent, side-chain, and backbone motions, so that both the backbone and side chains move on the instantaneous solvent potential of mean force (PMF) surface, and the backbone additionally moves on the side-chain PMF. Because of the arbitrariness in the choice of masses, this limit only needs to be reached approximately in practice. In particular, we show that a 10-fold reduction in solvent masses and a side-chain mass scale that is intermediate between the scaled solvent and the backbone lead to a quantitative enhancement in conformational sampling. [ABSTRACT FROM AUTHOR]

Intermolecular enol tautomers of Watson-Crick base pairs could emerge spontaneously via interbase double proton transfer. It has been hypothesized that their formation could be facilitated by thermal fluctuations and proton tunneling, and possibly be relevant to DNA damage. Theoretical and computational studies, assuming classical nuclei, have confirmed the dynamic stability of these rare tautomers. However, by accounting for nuclear quantum effects explicitly through Car-Parrinello path integral molecular dynamics calculations, we find the tautomeric enol form to be dynamically metastable, with lifetimes too insignificant to be implicated in DNA damage. [ABSTRACT FROM AUTHOR]

Adiabatic free energy dynamics (AFED) was introduced by Rosso et al. [J. Chem. Phys.2002, 116, 4389] for computing free energy profiles quickly and accurately using a dynamical adiabatic separation between a set of collective variables or reaction coordinates and the remaining degrees of freedom of a system. This approach has been shown to lead to a significant gain in efficiency versus traditional methods such as umbrella sampling, thermodynamic integration, and free energy perturbation for generating one-dimensional free energy profiles. More importantly, AFED is able to generate multidimensional free energy surfaces efficiently via full sweeps of the surface that rapidly map out the locations of the free energy minima. The most significant drawback to the AFED approach is the need to transform the coordinates into a generalized coordinate system that explicitly contains the collective variables of interest. Recently, Maragliano and Vanden-Eijnden built upon the AFED approach by introducing a set of extended phase-space variables, to which the adiabatic decoupling and high temperature are applied [Chem. Phys. Lett.2006, 426, 168]. In this scheme, which the authors termed "temperature accelerated molecular dynamics" or TAMD, the need for explicit coordinate transformations is circumvented. The ability of AFED and TAMD to generate free energy surfaces efficiently depends on the thermostatting mechanism employed, since both approaches are inherently nonequilibrium due to the adiabatic decoupling. Indeed, Maragliano and Vanden-Eijnden did not report any direct generation of free energy surfaces within the overdamped Langevin dynamics employed by these authors. Here, we show that by formulating TAMD in a manner that is closer to the original AFED approach, including the generalized Gaussian moment thermostat (GGMT) and multiple time-scale integration, multidimensional free energy surfaces for complex systems can be generated directly from the probability distribution function of the extended phase-space variables. The new TAMD formulation, which we term driven AFED or d-AFED, is applied to compare the conformational preferences of small peptides both in gas phase and in solution for three force fields. The results show that d-AFED/TAMD accurately and efficiently generates free energy surfaces in two collective variables useful for characterizing the conformations, namely, the radius of gyration, RG, and number of hydrogen bonds, NH. [ABSTRACT FROM AUTHOR]

The difficulties encountered in sampling of systems with rough energy landscapes using present methodology significantly limit the impact of simulation on molecular biology, in particular protein folding and design. Here, we present a major methodological development based on a promising new technique, the reference potential spatial warping algorithm (REPSWA) [Z. Zhu et al., Phys. Rev. Lett., 88 (2002), pp. 100201-100204], and present applications to several realistic systems. REPSWA works by introducing a variable transformation in the classical partition function that reduces the volume of phase space associated with a priori known barrier regions while increasing that associated with attractive basins. In this way, the partition function is preserved so that enhanced sampling is achieved without the need for reweighting phase-space averages. Here, a new class of transformations, designed to overcome the barriers induced by intermolecular/nonbonded interactions, whose locations are not known a priori, is introduced. The new transformations are designed to work in synergy with transformations originally introduced for overcoming intramolecular barriers. The new transformation adapts to the fluctuating local environment and is able to handle barriers that arise "on the fly." Thus, the new method is referred to as dynamic contact REPSWA (DC-REPSWA). In addition, combining hybrid Monte Carlo (HMC) with DC-REPSWA allows more aggressive sampling to take place. The combined DC-REPSWA-HMC method and its variants are shown to substantially enhance conformational sampling in long molecular chains composed of interacting single beads and beads with branches. The latter topologies characterize the united residue and united side chain representation of protein structures. [ABSTRACT FROM AUTHOR]

Understanding the chemical reaction mechanisms governing how small organic molecules attach to semiconductor surfaces can lead to new strategies for creating specific surface patterns such as single adduct monolayers. In this study, room-temperature ab initio molecular dynamics simulations of one and two 1 ,3-cyclohexadiene (CHD) molecule(s) reacting with the Si(1 OO)-2 xl surface reveal that adducts form via a carbocation-mediated two-step mechanism. Dimer flipping can either promote or prevent bond formation depending on how the CHD approaches. CHDs often travel past several Si dimers before finding the proper local environment. The resulting intermediate can persist for more than 4 ps, allowing the second bond to form with any adjacent Si dimer. The additional reactive site accounts for a large portion of the discrepancy between the predicted thermodynamic and observed experimental product distribution. Surface adducts protect a 5.6 Å region, direct unbound CHD exploration, and can cause adjacent dimers to flip. [ABSTRACT FROM AUTHOR]

A review is presented of the article "Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning" by Frank Noé, Jonas Köhler and others on the use of a machine learning methodology to address a problem with atomistic simulations of complex molecular systems using Boltzmann probability distributions to predict equilibrium properties.

The charge transport mechanism in solid ammonium perchlorate (AP) crystal exposed to an ammonia-rich environment is studied using ab initio molecular dynamics (AIMD). Ammonium perchlorate is an ionic crystal composed of NH4+ and ClO4− units that possesses an orthorhombic phase at T<513 K and a cubic phase at T>513 K. Exposure to an ammonia-rich atmosphere allows ammonia molecules to be absorbed into the crystal at interstitial sites. It has been proposed that these neutral ammonias can form short-lived N2H7+ complexes with the NH4+ ions allowing proton transfer between them, thereby enhancing the conductivity considerably. To date, however, there has been no direct evidence of this proposed mechanism. In this paper, ab initio molecular dynamics techniques are employed to explore this mechanism. By comparing computed infrared spectra of the pure and ammonia-doped crystals, we observe a significant broadening of the NH stretch peak into a lower frequency region, indicating through an experimentally verifiable observation, the formation of hydrogen bonds between NH3 and NH4+ units. This suggestion is confirmed by direction observation of N2H7+ complexes from the trajectory. Comparison of the diffusion constants of NH4+ in the pure and doped crystals yields a ratio that is comparable to the experimentally measured conductivity ratio and clearly shows an enhanced positive charge mobility. Finally, compelling evidence suggesting the possibility of an ammonia umbrella inversion following proton transfer from NH4+ and NH3 is obtained. [Copyright &y& Elsevier]

Mixed quantum mechanical/molecular mechanical (QM/MM) models, in which one part of a system is treated at an ab initio level while the rest is described by an empirical force field, are becoming increasingly popular. In this Letter, we describe a simple QM/MM model for methanol (CH3OH) in which the OH constitutes the ab initio region and is treated within a density functional representation, while the methyl group, the MM region, is treated using the AMBER force field. The model is used to study the structure of the bulk liquid and is shown to compare favorably with the experimental structure. [Copyright &y& Elsevier]

Examines the solution structures and transport mechanisms of hydrated hydroxide. Variation in the transport mechanism for protons and hydroxide ions; Mechanism of transport of hydrated proton in water; Assessment of the role of solvation complexes of hydroxide ions.

Determining collective variables (CVs) for conformational transitions is crucial to understanding their dynamics and targeting them in enhanced sampling simulations. Often, CVs are proposed based on intuition or prior knowledge of a system. However, the problem of systematically determining a proper reaction coordinate (RC) for a specific process in terms of a set of putative CVs can be achieved using committor analysis (CA). Identifying essential degrees of freedom that govern such transitions using CA remains elusive because of the high dimensionality of the conformational space. Various schemes exist to leverage the power of machine learning (ML) to extract an RC from CA. Here, we extend these studies and compare the ability of 17 different ML schemes to identify accurate RCs associated with conformational transitions. We tested these methods on an alanine dipeptide in vacuum and on a sarcosine dipeptoid in an implicit solvent. Our comparison revealed that the light gradient boosting machine method outperforms other methods. In order to extract key features from the models, we employed Shapley Additive exPlanations analysis and compared its interpretation with the "feature importance" approach. For the alanine dipeptide, our methodology identifies ϕ and θ dihedrals as essential degrees of freedom in the C7ax to C7eq transition. For the sarcosine dipeptoid system, the dihedrals ψ and ω are the most important for the cisαD to transαD transition. We further argue that analysis of the full dynamical pathway, and not just endpoint states, is essential for identifying key degrees of freedom governing transitions. [ABSTRACT FROM AUTHOR]

Fuel cell-based anion-exchange membranes (AEMs) and proton exchange membranes (PEMs) are considered to have great potential as cost-effective, clean energy conversion devices. However, a fundamental atomistic understanding of the hydroxide and hydronium diffusion mechanisms in the AEM and PEM environment is an ongoing challenge. In this work, we aim to identify the fundamental atomistic steps governing hydroxide and hydronium transport phenomena. The motivation of this work lies in the fact that elucidating the key design differences between the hydroxide and hydronium diffusion mechanisms will play an important role in the discovery and determination of key design principles for the synthesis of new membrane materials with high ion conductivity for use in emerging fuel cell technologies. To this end, ab initio molecular dynamics simulations are presented to explore hydroxide and hydronium ion solvation complexes and diffusion mechanisms in the model AEM and PEM systems at low hydration in confined environments. We find that hydroxide diffusion in AEMs is mostly vehicular, while hydronium diffusion in model PEMs is structural. Furthermore, we find that the region between each pair of cations in AEMs creates a bottleneck for hydroxide diffusion, leading to a suppression of diffusivity, while the anions in PEMs become active participants in the hydronium diffusion, suggesting that the presence of the anions in model PEMs could potentially promote hydronium diffusion. [ABSTRACT FROM AUTHOR]

Addresses questions concerning the character and anomalously high mobility of the hydrated excess proton in liquid water using computer simulation techniques. Widely held views of the proton's structure and transport; Conclusion that the proton's quantum-mechanical character plays a more important role than previously thought; Implications of results in areas related to aqueous acid-base chemistry, enzymatic proton transfers, proton transport in biological channels and atmospheric proton transfers on ice surfaces.

Explicit reversible integrators, suitable for use in large-scale computer simulations, are derived for extended systems generating the canonical and isothermal-isobaric ensembles. The new methods are compared with the standard implicit (iterative) integrators on some illustrative example problems. In addition, modification of the proposed algorithms for multiple time step integration is outlined. [ABSTRACT FROM AUTHOR]

Looks at the relative influence of thermal and quantum fluctuations on the proton transfer properties of the charged water complexes H5O2+ and H3O2-. Behavior of the shared protons in each complex; Cause of the differing behavior of the proton in each instance.