Elastic bag model for molecular dynamics simulations of solvated systems: application to liquid argon

A new approach is developed to study the dynamics of the localized process in solutions and other condensed phase systems. The approach employs a fluctuating elastic boundary (FEB) model which encloses the simulated system in an elastic bag that mimics the effects of the bulk solvent. This alleviates the need for periodic boundary conditions and allows for a reduction in the number of solvent molecules that need to be included in the simulation. The boundary bag is modeled as a mesh of quasi-particles connected by elastic bonds. The FEB model allows for volume and density fluctuations characteristic of the bulk system, and the shape of the boundary fluctuates during the course of the simulation to adapt to the configuration fluctuations of the explicit solute-solvent system inside. The method is applied to the simulation of a Lennard-Jones model of liquid argon. Various structural and dynamical quantities are computed and compared with those obtained from conventional periodic boundary simulations. The agreement between the two is excellent in most cases, thus validating the viability of the FEB method.     

Nonperiodic boundary conditions for solvated systems

The simulation of charged and/or strongly polar solutes represents a challenge for standard molecular-dynamics techniques. The use of periodic boundary conditions (PBCs) leads to artifacts due to the interaction between two replicas in the presence of the long-range Coulomb forces. A way to avoid these problems is the use of nonperiodic boundary conditions. A possible realization is to consider a finite system, a sphere, embedded in a reaction field described by the method of the images. In the present work the modified image approximation has been implemented in a molecular-dynamics code and optimized for the use of two standard solvents, water and acetonitrile. The methodology has then been applied to investigate the conformational changes in water-solvated alanine dipeptide. The free-energy surface calculated with this method is comparable to that obtained with PBC.     

Classical and quantum mechanical/molecular mechanical molecular dynamics simulations of alanine dipeptide in water: comparisons with IR and vibrational circular dichroism spectra

We have implemented the combined quantum mechanical (QM)/molecular mechanical (MM) molecular dynamics (MD) simulations of alanine dipeptide in water along with the polarizable and nonpolarizable classical MD simulations with different models of water. For the QM/MM MD simulation, the alanine dipeptide is treated with the AM1 or PM3 approximations and the fluctuating solute dipole moment is calculated by the Mulliken population analysis. For the classical MD simulations, the solute is treated with the polarizable or nonpolarizable AMBER and polarizable CHARMM force fields and water is treated with the TIP3P, TIP4P, or TIP5P model. It is found that the relative populations of right-handed alpha-helix and extended beta and P(II) conformations in the simulation trajectory strongly depend on the simulation method. For the QM/MM MD simulations, the PM3/MM shows that the P(II) conformation is dominant, whereas the AM1/MM predicts that the dominant conformation is alpha(R). Polarizable CHARMM force field gives almost exclusively P(II) conformation and other force fields predict that both alpha-helical and extended (beta and P(II)) conformations are populated with varying extents. Solvation environment around the dipeptide is investigated by examining the radial distribution functions and numbers and lifetimes of hydrogen bonds. Comparing the simulated IR and vibrational circular dichroism spectra with experimental results, we concluded that the dipeptide adopts the P(II) conformation and PM3/MM, AMBER03 with TIP4P water, and AMBER polarizable force fields are acceptable for structure determination of the dipeptide considered in this paper.     

Replica exchange with nonequilibrium switches: enhancing equilibrium sampling by increasing replica overlap

We describe a replica exchange strategy where trial swap configurations are generated by nonequilibrium switching simulations. By devoting simulation time to the switching simulations, one can systematically increase an effective overlap between replicas, which leads to an increased exchange acceptance rate and less correlated equilibrium samples. In this paper, we derive our method for a general class of stochastic dynamics, and discuss various strategies for enhancing replica overlap through novel dynamical schemes and prudent choices of switching protocols. We then demonstrate our method on a model system of alanine dipeptide in implicit solvent, characterizing decreases in data correlations and gains in sampling efficiency.     

Free energy simulations: Applications to the study of liquid water,hydrophobic interactions and solvent effects on conformational stability

Free energy simulations using the Metropolis Monte Carlo method and the coupling parameter approach with umbrella sampling are described for several problems of interest in structural biochemistry: the liquid water, the hydrophobic interaction of alkyl and phenyl groups in water, and solvent effects on the conformational stability of the alanine dipeptide and the dimethyl phosphate anion in water. Proximity analysis of the results is employed to identify stabilizing factors. Implications of the result with respect to the structural chemistry of proteins and nucleic acids is considered.     

Infrared and vibrational CD spectra of partially solvated alpha-helices: DFT-based simulations with explicit solvent

Theoretical simulations are used to investigate the effects of aqueous solvent on the vibrational spectra of model alpha-helices, which are only partly exposed to solvent to mimic alpha-helices in proteins. Infrared absorption (IR) and vibrational circular dichroism (VCD) amide I' spectra for 15-amide alanine alpha-helices are simulated using density functional theory (DFT) calculations combined with the property transfer method. The solvent is modeled by explicit water molecules hydrogen bonded to the solvated amide groups. Simulated spectra for two partially solvated model alpha-helices, one corresponding to a more exposed and the other to a more buried structure, are compared to the fully solvated and unsolvated (gas phase) simulations. The dependence of the amide I spectra on the orientation of the partially solvated helix with respect to the solvent and effects of solvation on the amide I' of 13C isotopically substituted alpha-helices are also investigated. The partial exposure to solvent causes significant broadening of the amide I' bands due to differences in the vibrational frequencies of the explicitly solvated and unsolvated amide groups. The different degree of partial solvation is reflected primarily in the frequency shifts of the unsolvated (buried) amide group vibrations. Depending on which side of the alpha-helix is exposed to solvent, the simulated IR band-shapes exhibit significant changes, from broad and relatively featureless to distinctly split into two maxima. The simulated amide I' VCD band-shapes for the partially solvated alpha-helices parallel the broadening of the IR and exhibit more sign variation, but generally preserve the sign pattern characteristic of the alpha-helical structures and are much less dependent on the alpha-helix orientation with respect to the solvent. The simulated amide I' IR spectra for the model peptides with explicitly hydrogen-bonded water are consistent with the experimental data for small alpha-helical proteins at very low temperatures, but overestimate the effects of solvent on the protein spectra at ambient temperatures, where the peptide-water hydrogen bonds are weakened by thermal motion.     

A microscopic view of substitution reactions solvated by ionic liquids

The solvation effect of the ionic liquid 1-N-butyl-3-methylimidazolium hexafluorophosphate on nucleophilic substitution reactions of halides toward the aliphatic carbon of methyl p-nitrobenzenesulfonate (pNBS) was investigated by computer simulations. The calculations were performed by using a hybrid quantum-mechanical/molecular-mechanical (QM/MM) methodology. A semiempirical Hamiltonian was first parametrized on the basis of comparison with ab initio calculations for Cl(-) and Br(-) reaction with pNBS at gas phase. In condensed phase, free energy profiles were obtained for both reactions. The calculated reaction barriers are in agreement with experiment. The structure of species solvated by the ionic liquid was followed along the reaction progress from the reagents, through the transition state, to the final products. The simulations indicate that this substitution reaction in the ionic liquid is slower than in nonpolar molecular solvents proper to significant stabilization of the halide anion by the ionic liquid in comparison with the transition state with delocalized charge. Solute-solvent interactions in the first solvation shell contain several hydrogen bonds that are formed or broken in response to charge density variation along the reaction coordinate. The detailed structural analysis can be used to rationalize the design of new ionic liquids with tailored solvation properties.     

Moving solvated electrons with light: nonadiabatic mixed quantum/classical molecular dynamics simulations of the relocalization of photoexcited solvated electrons in tetrahydrofuran (THF)

Motivated by recent ultrafast spectroscopic experiments [Martini et al., Science 293, 462 (2001)], which suggest that photoexcited solvated electrons in tetrahydrofuran (THF) can relocalize (that is, return to equilibrium in solvent cavities far from where they started), we performed a series of nonequilibrium, nonadiabatic, mixed quantum/classical molecular dynamics simulations that mimic one-photon excitation of the THF-solvated electron. We find that as photoexcited THF-solvated electrons relax to their ground states either by continuous mixing from the excited state or via nonadiabatic transitions, approximately 30% of them relocalize into cavities that can be over 1 nm away from where they originated, in close agreement with the experiments. A detailed investigation shows that the ability of excited THF-solvated electrons to undergo photoinduced relocalization stems from the existence of preexisting cavity traps that are an intrinsic part of the structure of liquid THF. This explains why solvated electrons can undergo photoinduced relocalization in solvents like THF but not in solvents like water, which lack the preexisting traps necessary to stabilize the excited electron in other places in the fluid. We also find that even when they do not ultimately relocalize, photoexcited solvated electrons in THF temporarily visit other sites in the fluid, explaining why the photoexcitation of THF-solvated electrons is so efficient at promoting recombination with nearby scavengers. Overall, our study shows that the defining characteristic of a liquid that permits the photoassisted relocalization of solvated electrons is the existence of nascent cavities that are attractive to an excess electron; we propose that other such liquids can be found from classical computer simulations or neutron diffraction experiments.     

On the proper calculation of electrostatic interactions in solid-supported bilayer systems

Modeling systems that are not inherently isotropic, e.g., extended bilayers, using molecular simulation techniques poses a potential problem. Since these methods rely on a finite number of atoms and molecules to describe the system, periodic boundary conditions are implemented to avoid edge effects and capture long-range electrostatic interactions. Systems consisting of a solvated bilayer adsorbed on a solid surface and exposed to an air/vacuum interface occur in many experimental settings and present some unique challenges in this respect. Here, we investigated the effects of implementing different electrostatic boundary conditions on the structural and electrostatic properties of a quartz/water/vacuum interface and a similar quartz-supported hydrated lipid bilayer exposed to vacuum. Since these interfacial systems have a net polarization, implementing the standard Ewald summation with the conducting boundary condition for the electrostatic long-range interactions introduced an artificial periodicity in the out-of-plane dimension. In particular, abnormal orientational polarizations of water were observed with the conducting boundary condition. Implementing the Ewald summation technique with the planar vacuum boundary condition and calculating electrostatic properties compatible with the implemented electrostatic boundary condition removed these inconsistencies. This formulation is generally applicable to similar interfacial systems in bulk solution.     

Role of solvent in determining conformational preferences of alanine dipeptide in water

Evidence from a variety of spectroscopic probes indicates that (phi, psi) values corresponding to the left-handed polyproline II helix (P(II)) are preferred for short alanine-based peptides in water. On the basis of results from theoretical studies, it is believed that the observed preference is dictated by favorable peptide-solvent interactions, which are realized through formation of optimal hydrogen-bonding water bridges between peptide donor and acceptor groups. In the present study, we address this issue explicitly by analyzing the hydration structure and thermodynamics of 16 low-energy conformers of the alanine dipeptide (N-acetylalanine-N'-methylamide) in liquid water. Monte Carlo simulations in the canonical ensemble were performed under ambient conditions with all-atom OPLS parameters for the alanine dipeptide and the TIP5P model for water. We find that the number of hydrogen-bonded water molecules connecting the peptide group donor and acceptor atoms has no effect on the solvation thermodynamics. Instead, the latter are determined by the work done to fully hydrate the peptide. This work is minimal for conformations that are characterized by a minimal overlap of the primary hydration shells around the peptide donor and acceptor atoms. As a result, peptide-solvent interactions favor "compact" conformations that do not include P(II)-like geometries. Our main conclusion is that the experimentally observed preference for P(II) does not arise due to favorable direct interactions between the peptide and water molecules. Instead, the latter act to unmask underlying conformational preferences that are a consequence of minimizing intrapeptide steric conflicts.     

First principles and classical molecular dynamics simulations of solvated benzene

We have performed extensive ab initio and classical molecular dynamics (MD) simulations of benzene in water in order to examine the unique solvation structures that are formed. Qualitative differences between classical and ab initio MD simulations are found and the importance of various technical simulation parameters is examined. Our comparison indicates that nonpolarizable classical models are not capable of describing the solute-water interface correctly if local interactions become energetically comparable to water hydrogen bonds. In addition, a comparison is made between a rigid water model and fully flexible water within ab initio MD simulations which shows that both models agree qualitatively for this challenging system.     

An efficient parallelization scheme for molecular dynamics simulations with many-body, flexible, polarizable empirical potentials: application to water

An efficient parallelization scheme for classical molecular dynamics simulations with flexible, polarizable empirical potentials is presented. It is based on the standard Ewald summation technique to handle the long-range electrostatic and induction interactions. The algorithm for this parallelization scheme is designed for systems containing several thousands of polarizable sites in the simulation box. Its performance is evaluated during molecular dynamics simulations under periodic boundary conditions with unit cell sizes ranging from 128 to 512 molecules employing two flexible polarizable water models [DC(F) and TTM2.1-F] containing 1 and 3 polarizable sites, respectively. The time-to-solution for these two polarizable models is compared with the one for a flexible, pairwise-additive water model (TIP4F). The benchmarks were performed on both shared and distributed memory platforms. As a result of the efficient calculation of the induced dipole moments, a superlinear scaling as a function of the number of the processors is observed. To the best of our knowledge, this is the first reported results of parallel scaling and performance for simulations of liquid water with a polarizable potential under periodic boundary conditions.     

Efficient multiple time scale molecular dynamics: Using colored noise thermostats to stabilize resonances

Multiple time scale molecular dynamics enhances computational efficiency by updating slow motions less frequently than fast motions. However, in practice, the largest outer time step possible is limited not by the physical forces but by resonances between the fast and slow modes. In this paper we show that this problem can be alleviated by using a simple colored noise thermostatting scheme which selectively targets the high frequency modes in the system. For two sample problems, flexible water and solvated alanine dipeptide, we demonstrate that this allows the use of large outer time steps while still obtaining accurate sampling and minimizing the perturbation of the dynamics. Furthermore, this approach is shown to be comparable to constraining fast motions, thus providing an alternative to molecular dynamics with constraints.     

A systematic study of minima in alanine dipeptide

The alanine dipeptide is a standard system to model dihedral angles in proteins. It is shown that obtaining the Ramachandran plot accurately is a hard problem because of many local minima; depending on the details of geometry optimizations, different Ramachandran plots can be obtained. To locate all energy minima, starting from geometries from MD simulations, 250,000 geometry optimizations were performed at the level of RHF/6-31G*, followed by re-optimizations of the located 827 minima at the level of MP2/6–311++G**, yielding 30 unique minima, most of which were not previously reported in literature. Both in vacuo and solvated structures are discussed. The minima are systematically categorized based on four backbone dihedral angles. The Gibbs energies are evaluated and the structural factors determining the relative stabilities of conformers are discussed. © 2018 Wiley Periodicals, Inc.     

Efficient implementation of periodic boundary conditions in Monte Carlo simulation

Determination of the shortest distances between particles is one of the most time-consuming parts of molecular simulation by the Monte Carlo method. In this work, we demonstrate that the use of signed-integer storage of coordinates in a scaled box allows one to skip multiple conditional statements in realization of periodic boundary conditions in cubic and rectangular boxes, which, in turn, increases the performance. Performance of the improved procedure was tested in NVT Monte Carlo simulations for liquid krypton and water. © 2018 Wiley Periodicals, Inc.     

Independent pairs and Monte-Carlo simulations of the geminate recombination of solvated electrons in liquid-to-supercritical water

Independent pairs (IP) and Monte Carlo (MC) simulations are employed to model experimental femtosecond time-resolved pump-probe spectroscopic data on the geminate recombination dynamics of solvated electrons in liquid-to-supercritical water. The hydrated electron was created by two-photon ionization of the neat fluid with a total ionization energy of 9.3 eV. In both numerical approaches, the ejection length, , (i.e. the distance from the ionization core, at which the electron is thermally and spatially localized) is used as the primary adjustable fitting parameter that can bring both model simulations into quantitative agreement with the ultrafast kinetic experiment. The influence of the thermodynamic conditions on the ejection length and on the recombination mechanism is discussed. Whereas in the compressed liquid associated with a high dielectric constant (ε ≥ 20), the electron recombines predominantly with the OH radical, the dissociative recombination via charge neutralization with the hydronium cation takes over at small dielectric constants (ε < 20). The importance of charge-dipole interactions for Monte-Carlo simulations of the recombination reactions of the hydrated electrons in the low-permittivity region is stressed.     

Stochastic boundary conditions for molecular dynamics simulations of ST2 water

The deformable stochastic boundary method developed previously for treating simple liquids without periodic boundary conditions, is extended to the ST2 model of water. The method is illustrated by a molecular dynamics simulation of a sphere containing 98 water molecules. Comparison with the results of the periodic boundary simulation by Stillinger and Rahman shows very good agreement for structural and dynamic properties.