Search

This letter illustrates the opinion of the molecular dynamics (MD) community on the need to adopt a new FAIR paradigm for the use of molecular simulations. It highlights the necessity of a collaborative effort to create, establish, and sustain a database that allows findability, accessibility, interoperability, and reusability of molecular dynamics simulation data. Such a development would democratize the field and significantly improve the impact of MD simulations on life science research. This will transform our working paradigm, pushing the field to a new frontier. We invite you to support our initiative at the MDDB community (https://mddbr.eu/community/)

The conformational dynamics of single-stranded nucleic acids are fundamental for nucleic acid folding and function. However, their elementary chain dynamics have been difficult to resolve experimentally. Here we employ a combination of single-molecule F\"orster resonance energy transfer, nanosecond fluorescence correlation spectroscopy, fluorescence lifetime analysis, and nanophotonic enhancement to determine the conformational ensembles and rapid chain dynamics of short single-stranded nucleic acids in solution. To interpret the experimental results in terms of end-to-end distance dynamics, we utilize the hierarchical chain growth approach, simple polymer models, and refinement with Bayesian inference of ensembles to generate structural ensembles that closely align with the experimental data. The resulting chain reconfiguration times are exceedingly rapid, in the 10-ns range. Solvent viscosity-dependent measurements indicate that these dynamics of single-stranded nucleic acids exhibit negligible internal friction and are thus dominated by solvent friction. Our results provide a detailed view of the conformational distributions and rapid dynamics of single-stranded nucleic acids.

In molecular dynamics simulations in the NPT ensemble at constant pressure, the size and shape of the periodic simulation box fluctuate with time. For particle images far from the origin, the rescaling of the box by the barostat results in unbounded position displacements. Special care is thus required when a particle trajectory is unwrapped from a projection into the central box under periodic boundary conditions to a trajectory in full three-dimensional space, e.g., for the calculation of diffusion coefficients. Here, we review and compare different schemes in use for trajectory unwrapping. We also specify the corresponding rewrapping schemes to put an unwrapped trajectory back into the central box. On this basis, we then identify a scheme for the calculation of meaningful diffusion coefficients, which is a primary application of trajectory unwrapping. In this scheme, the wrapped and unwrapped trajectory are mutually consistent and their statistical properties are preserved. We conclude with advice on best practice for the consistent unwrapping of constant-pressure simulation trajectories and the calculation of accurate translational diffusion coefficients., Comment: 30 pages, 8 figures

Markov-chain Monte Carlo algorithms rely on trial moves that are either rejected or accepted based on certain criteria. Here, we provide an efficient algorithm to generate random rotation matrices in four dimensions (4D) covering an arbitrary pre-defined range of rotation angles. The matrices can be combined with Monte Carlo methods for the efficient sampling of the SO(4) group of 4D rotations. The matrices are unbiased and constructed such that repeated rotations result in uniform sampling over SO(4). 4D rotations can be used to optimize the mass partitioning for stable time integration in coarse-grained molecular dynamics simulations and should find further applications in the fields of robotics and computer vision., Comment: 9 pages, 2 figures, 1 table

Disordered proteins and nucleic acids play key roles in cellular function and disease. Here we review recent advances in the computational exploration of the conformational dynamics of flexible biomolecules. We focus on hierarchical chain growth (HCG) from fragment libraries built with atomistic molecular dynamics simulations. HCG combines chain fragments in a statistically reproducible manner into ensembles of full-length atomically detailed biomolecular structures. The input fragment structures are typically collected from molecular dynamics simulations, but could also come from structural databases. Experimental data can be integrated during and after chain assembly. Applications to the neurodegeneration-linked proteins $\alpha$-synuclein, tau, and TDP-43, including as condensate, illustrate the use of HCG. We conclude by highlighting the emerging connections to AI-based structural modeling., Comment: 21 pages, 4 figures

Analysis of bond rupture data from single-molecule force spectroscopy experiments commonly relies on the strong assumption that the bond dissociation process is irreversible. However, with increased spatiotemporal resolution of instruments it is now possible to observe multiple unbinding-rebinding events in a single pulling experiment. Here, we augment the theory of force-induced unbinding by explicitly taking into account rebinding kinetics, and provide approximate analytic solutions of the resulting rate equations. Furthermore, we use a short-time expansion of the exact kinetics to construct numerically efficient maximum likelihood estimators for the parameters of the force-dependent unbinding and rebinding rates, which pair well with and complement established methods, such as the analysis of rate maps. We provide an open-source implementation of the theory, evaluated for Bell-like rates, which we apply to synthetic data generated by a Gillespie stochastic simulation algorithm for time-dependent rates., Comment: 14 pages, 5 figures, 1 table. The article has been submitted to Physical Review Research

Among other improvements, the Martini 3 coarse-grained force field provides a more accurate description of the solvation of protein pockets and channels through the consistent use of various bead types and sizes. Here, we show that the representation of Na$^+$ and Cl$^-$ ions as "tiny" (TQ5) beads limits the accessible time step to 25 fs. By contrast, with Martini 2, time steps of 30-40 fs were possible for lipid bilayer systems without proteins. This limitation is relevant for, e.g., phase separating lipid mixtures that require long equilibration times. We derive a quantitative kinetic model of time-integration instabilities in molecular dynamics (MD) as a function of time step, ion concentration and mass, system size, and simulation time. With this model, we demonstrate that ion-water interactions are the main source of instability at physiological conditions, followed closely by ion-ion interactions. We show that increasing the ionic masses makes it possible to use time steps up to 40 fs with minimal impact on static equilibrium properties and on dynamical quantities such as lipid and ion diffusion coefficients. Increasing the size of the bead representing the ions (and thus changing their hydration) also permits longer time steps. The use of larger time steps in Martini 3 simulations results in a more efficient exploration of configuration space. The kinetic model of MD simulation crashes can be used to determine the maximum allowed time step whenever sampling efficiency is critical., Comment: Manuscript: 28 pages, 6 figures, 6 tables. Supporting Information: 16 pages, 10 figures, 7 tables. The article has been submitted for publication in The Journal of Chemical Physics

Simulations with an adaptive time-dependent bias, such as metadynamics, enable an efficient exploration of the conformational space of a system. However, the dynamic information of the system is altered by the bias. With infrequent metadynamics it is possible to recover the transition rate of crossing a barrier, if the collective variables are ideal and there is no bias deposition near the transition state. Unfortunately, for simulations of complex molecules, these conditions are not always fulfilled. To overcome these limitations, and inspired by single-molecule force spectroscopy, we developed a method based on Kramers' theory for calculating the barrier-crossing rate when a time-dependent bias is added to the system. We assess the quality of the bias parameter by measuring how efficiently the bias accelerates the transitions compared to ideal behavior. We present approximate analytical expressions of the survival probability that accurately reproduce the barrier-crossing time statistics, and enable the extraction of the unbiased transition rate even for challenging cases, where previous methods fail.

Molecular self-organization driven by concerted many-body interactions produces the ordered structures that define both inanimate and living matter. Understanding the physical mechanisms that govern the formation of molecular complexes and crystals is key to controlling the assembly of nanomachines and new materials. We present an artificial intelligence (AI) agent that uses deep reinforcement learning and transition path theory to discover the mechanism of molecular self-organization phenomena from computer simulations. The agent adaptively learns how to sample complex molecular events and, on the fly, constructs quantitative mechanistic models. By using the mechanistic understanding for AI-driven sampling, the agent closes the learning cycle and overcomes time-scale gaps of many orders of magnitude. Symbolic regression condenses the mechanism into a human-interpretable form. Applied to ion association in solution, gas-hydrate crystal formation, and membrane-protein assembly, the AI agent identifies the many-body solvent motions governing the assembly process, discovers the variables of classical nucleation theory, and reveals competing assembly pathways. The mechanistic descriptions produced by the agent are predictive and transferable to close thermodynamic states and similar systems. Autonomous AI sampling has the power to discover assembly and reaction mechanisms from materials science to biology.

Single-molecule localization microscopy allows practitioners to locate and track labeled molecules in biological systems. When extracting diffusion coefficients from the resulting trajectories, it is common practice to perform a linear fit on mean-square-displacement curves. However, this strategy is suboptimal and prone to errors. Recently, it was shown that the increments between observed positions provide a good estimate for the diffusion coefficient, and their statistics are well-suited for likelihood-based analysis methods. Here, we revisit the problem of extracting diffusion coefficients from single-particle tracking experiments subject to static and dynamic noise using the principle of maximum likelihood. Taking advantage of an efficient real-space formulation, we extend the model to mixtures of subpopulations differing in their diffusion coefficients, which we estimate with the help of the expectation-maximization algorithm. This formulation naturally leads to a probabilistic assignment of trajectories to subpopulations. We employ the theory to analyze experimental tracking data that cannot be explained with a single diffusion coefficient. We test how well a dataset conforms to the assumptions of a diffusion model and determine the optimal number of subpopulations with the help of a quality factor of known analytical distribution. To facilitate use by practitioners, we provide a fast open-source implementation of the theory for the efficient analysis of multiple trajectories in arbitrary dimensions simultaneously., Comment: 19 pages, 8 figures, 1 table. The article has been accepted for publication in The Journal of Chemical Physics

Translational diffusion coefficients are routinely estimated from molecular dynamics simulations. Linear fits to mean squared displacement (MSD) curves have become the de facto standard, from simple liquids to complex biomacromolecules. Nonlinearities in MSD curves at short times are handled with a wide variety of ad hoc practices, such as partial and piece-wise fitting of the data. Here, we present a rigorous framework to obtain reliable estimates of the diffusion coefficient and its statistical uncertainty. We also assess in a quantitative manner if the observed dynamics is indeed diffusive. By accounting for correlations between MSD values at different times, we reduce the statistical uncertainty of the estimator and thereby increase its efficiency. With a Kolmogorov-Smirnov test, we check for possible anomalous diffusion. We provide an easy-to-use Python data analysis script for the estimation of diffusion coefficients. As an illustration, we apply the formalism to molecular dynamics simulation data of pure TIP4P-D water and a single ubiquitin protein. In a companion paper [J. Chem. Phys. XXX, YYYYY (2020)], we demonstrate its ability to recognize deviations from regular diffusion caused by systematic errors in a common trajectory "unwrapping" scheme that is implemented in popular simulation and visualization software., Comment: 15 pages, 6 figures. The following article has been accepted for publication at The Journal of Chemical Physics

In molecular dynamics simulations under periodic boundary conditions, particle positions are typically wrapped into a reference box. For diffusion coefficient calculations using the Einstein relation, the particle positions need to be unwrapped. Here, we show that a widely used heuristic unwrapping scheme is not suitable for long simulations at constant pressure. Improper accounting for box-volume fluctuations creates, at long times, unphysical trajectories and, in turn, grossly exaggerated diffusion coefficients. We propose an alternative unwrapping scheme that resolves this issue. At each time step, we add the minimal displacement vector according to periodic boundary conditions for the instantaneous box geometry. Here and in a companion paper [J. Chem. Phys. XXX, YYYYY (2020)], we apply the new unwrapping scheme to extensive molecular dynamics and Brownian dynamics simulation data. We provide practitioners with a formula to assess if and by how much earlier results might have been affected by the widely used heuristic unwrapping scheme., Comment: 6 pages, 5 figures. The following article has been accepted for publication at The Journal of Chemical Physics

Cryo-electron microscopy is a revolutionary technique that can provide 3D density maps at near-atomic resolution. However, map validation is still an open issue in the field. Despite several efforts from the community, it is possible to overfit the reconstructions to noisy data. Here, inspired by modern statistics, we develop a novel methodology that uses a small independent particle set to validate the 3D maps. The main idea is to monitor how the map probability evolves over the control set during the refinement. The method is complementary to the gold-standard procedure, which generates two reconstructions at each iteration. We low-pass filter the two reconstructions for different frequency cutoffs, and we calculate the probability of each filtered map given the control set. For high-quality maps, the probability should increase as a function of the frequency cutoff and of the refinement iteration. We also compute the similarity between the probability distributions of the two reconstructions. As higher frequencies are added to the maps, more dissimilar are the distributions. We optimized the BioEM software package to perform these calculations, and tested the method on several systems, some which were overfitted. Our results show that our method is able to discriminate the overfitted sets from the non-overfitted ones. We conclude that having a control particle set, not used for the refinement, is essential for cross-validating cryo-EM maps.

In single-molecule force spectroscopy experiments, a biomolecule is attached to a force probe via polymer linkers, and the total extension -- of molecule plus apparatus -- is monitored as a function of time. In a typical unfolding experiment at constant force, the total extension jumps between two values that correspond to the folded and unfolded states of the molecule. For several biomolecular systems the committor, which is the probability to fold starting from a given extension, has been used to extract the molecular activation barrier (a technique known as "committor inversion"). In this work, we study the influence of the force probe, which is much larger than the molecule being measured, on the activation barrier obtained by committor inversion. We use a two-dimensional framework in which the diffusion coefficient of the molecule and of the pulling device can differ. We systematically study the free energy profile along the total extension obtained from the committor, by numerically solving the Onsager equation and using Brownian dynamics simulations. We analyze the dependence of the extracted barrier on the linker stiffness, molecular barrier height, and diffusion anisotropy, and thus, establish the range of validity of committor inversion. Along the way, we showcase the committor of 2-dimensional diffusive models and illustrate how it is affected by barrier asymmetry and diffusion anisotropy., Comment: 22 pages, 7 figures

Exascale computing holds great opportunities for molecular dynamics (MD) simulations. However, to take full advantage of the new possibilities, we must learn how to focus computational power on the discovery of complex molecular mechanisms, and how to extract them from enormous amounts of data. Both aspects still rely heavily on human experts, which becomes a serious bottleneck when a large number of parallel simulations have to be orchestrated to take full advantage of the available computing power. Here, we use artificial intelligence (AI) both to guide the sampling and to extract the relevant mechanistic information. We combine advanced sampling schemes with statistical inference, artificial neural networks, and deep learning to discover molecular mechanisms from MD simulations. Our framework adaptively and autonomously initializes simulations and learns the sampled mechanism, and is thus suitable for massively parallel computing architectures. We propose practical solutions to make the neural networks interpretable, as illustrated in applications to molecular systems., Comment: 11 pages, 5 figures, supporting information

By performing molecular dynamics simulations with up to 132 million coarse-grained particles in half-micron sized boxes, we show that hydrodynamics quantitatively explains the finite-size effects on diffusion of lipids, proteins, and carbon nanotubes in membranes. The resulting Oseen correction allows us to extract infinite-system diffusion coefficients and membrane surface viscosities from membrane simulations despite the logarithmic divergence of apparent diffusivities with increasing box width. The hydrodynamic theory of diffusion applies also to membranes with asymmetric leaflets and embedded proteins, and to a complex plasma-membrane mimetic.

We describe and implement iMapD, a computer-assisted approach for accelerating the exploration of uncharted effective Free Energy Surfaces (FES), and more generally for the extraction of coarse-grained, macroscopic information from atomistic or stochastic (here Molecular Dynamics, MD) simulations. The approach functionally links the MD simulator with nonlinear manifold learning techniques. The added value comes from biasing the simulator towards new, unexplored phase space regions by exploiting the smoothness of the (gradually, as the exploration progresses) revealed intrinsic low-dimensional geometry of the FES.

In cryo-electron microscopy (EM), molecular structures are determined from large numbers of projection images of individual particles. To harness the full power of this single-molecule information, we use the Bayesian inference of EM (BioEM) formalism. By ranking structural models using posterior probabilities calculated for individual images, BioEM in principle addresses the challenge of working with highly dynamic or heterogeneous systems not easily handled in traditional EM reconstruction. However, the calculation of these posteriors for large numbers of particles and models is computationally demanding. Here we present highly parallelized, GPU-accelerated computer software that performs this task efficiently. Our flexible formulation employs CUDA, OpenMP, and MPI parallelization combined with both CPU and GPU computing. The resulting BioEM software scales nearly ideally both on pure CPU and on CPU+GPU architectures, thus enabling Bayesian analysis of tens of thousands of images in a reasonable time. The general mathematical framework and robust algorithms are not limited to cryo-electron microscopy but can be generalized for electron tomography and other imaging experiments.

We present a new method that enables the identification and analysis of both transition and metastable conformational states from atomistic or coarse-grained molecular dynamics (MD) trajectories. Our algorithm is presented and studied by using both analytical and actual examples from MD simulations of the helix-forming peptide Ala5, and of a larger system, the epidermal growth factor receptor (EGFR) protein. In all cases, our method identifies automatically the corresponding transition states and metastable conformations in an optimal way, with the input of a set of relevant coordinates, by capturing accurately the intrinsic slowest relaxation rate. Our approach provides a general and easy to implement analysis method that provides unique insight into the molecular mechanism and the rare but crucial rate limiting conformational pathways occurring in complex dynamical systems such as molecular trajectories.

We describe different Bayesian ensemble refinement methods, examine their interrelation, and discuss their practical application. With ensemble refinement, the properties of dynamic and partially disordered (bio)molecular structures can be characterized by integrating a wide range of experimental data, including measurements of ensemble-averaged observables. We start from a Bayesian formulation in which the posterior is a functional that ranks different configuration space distributions. By maximizing this posterior, we derive an optimal Bayesian ensemble distribution. For discrete configurations, this optimal distribution is identical to that obtained by the maximum entropy "ensemble refinement of SAXS" (EROS) formulation. Bayesian replica ensemble refinement enhances the sampling of relevant configurations by imposing restraints on averages of observables in coupled replica molecular dynamics simulations. We show that the strength of the restraint should scale linearly with the number of replicas to ensure convergence to the optimal Bayesian result in the limit of infinitely many replicas. In the "Bayesian inference of ensembles" (BioEn) method, we combine the replica and EROS approaches to accelerate the convergence. An adaptive algorithm can be used to sample directly from the optimal ensemble, without replicas. We discuss the incorporation of single-molecule measurements and dynamic observables such as relaxation parameters. The theoretical analysis of different Bayesian ensemble refinement approaches provides a basis for practical applications and a starting point for further investigations., Comment: Paper submitted to The Journal of Chemical Physics (15 pages, 4 figures); updated references; expanded discussions of related formalisms, error treatment, and ensemble refinement with and without replicas; appendix

We calculate the pair diffusion coefficient D(r) as a function of the distance r between two hard-sphere particles in a dense monodisperse suspension. The distance-dependent pair diffusion coefficient describes the hydrodynamic interactions between particles in a fluid that are central to theories of polymer and colloid dynamics. We determine D(r) from the propagators (Green's functions) of particle pairs obtained from discontinuous molecular dynamics simulations. At distances exceeding 3 molecular diameters, the calculated pair diffusion coefficients are in excellent agreement with predictions from exact macroscopic hydrodynamic theory for large Brownian particles suspended in a solvent bath, as well as the Oseen approximation. However, the asymptotic 1/r distance dependence of D(r) associated with hydrodynamic effects emerges only after the pair distance dynamics has been followed for relatively long times, indicating non-negligible memory effects in the pair diffusion at short times. Deviations of the calculated D(r) from the hydrodynamic models at short distances r reflect the underlying many-body fluid structure, and are found to be correlated to differences in the local available volume. The procedure used here to determine the pair diffusion coefficients can also be used for single-particle diffusion in confinement with spherical symmetry., Comment: 6 pages, 5 figures

We present a recently developed one-dimensional dipole lattice model that accurately captures the key properties of water in narrow nanopores. For this model, we derive three equivalent representations of the Hamiltonian that together yield a transparent physical picture of the energetics of the water chain and permit efficient computer simulations. In the charge representation, the Hamiltonian consists of nearest-neighbor interactions and Coulomb-like interactions of effective charges at the ends of dipole ordered segments. Approximations based on the charge picture shed light on the influence of the Coulomb-like interactions on the structure of nanopore water. We use Monte Carlo simulations to study the system behavior of the full Hamiltonian and its approximations as a function of chemical potential and system size and investigate the bimodal character of the density distribution occurring at small system sizes.

We describe a reverse integration approach for the exploration of low-dimensional effective potential landscapes. Coarse reverse integration initialized on a ring of coarse states enables efficient "navigation" on the landscape terrain: escape from local effective potential wells, detection of saddle points, and identification of significant transition paths between wells. We consider several distinct ring evolution modes: backward stepping in time, solution arc--length, and effective potential. The performance of these approaches is illustrated for a deterministic problem where the energy landscape is known explicitly. Reverse ring integration is then applied to "noisy" problems where the ring integration routine serves as an outer "wrapper" around a forward-in-time inner simulator. Three versions of such inner simulators are considered: a system of stochastic differential equations, a Gillespie--type stochastic simulator, and a molecular dynamics simulator. In these "equation-free" computational illustrations, estimation techniques are applied to the results of short bursts of "inner" simulation to obtain the unavailable (in closed form) quantities (local drift and diffusion coefficient estimates) required for reverse ring integration; this naturally leads to approximations of the effective landscape., Comment: article submitted to J Chem Phys

We study the diffusive dynamics of a hard-sphere fluid confined between parallel smooth hard walls. The position-dependent diffusion coefficient normal to the walls is larger in regions of high local packing density. High density regions also have the largest available volume, consistent with the fast local diffusivity. Indeed, local and global diffusivities as a function of the Widom insertion probability approximately collapse onto a master curve. Parallel and average normal diffusivities are strongly coupled at high densities and deviate from bulk fluid behavior.

We show how accurate kinetic information, such as the rates of protein folding and unfolding, can be extracted from replica-exchange molecular dynamics (REMD) simulations. From the brief and continuous trajectory segments between replica exchanges, we estimate short-time propagators in conformation space and use them to construct a master equation. For a helical peptide in explicit water, we determine the rates of transitions both locally between microscopic conformational states and globally for folding and unfolding. We show that accurate rates in the ~1/(100 ns) to ~1/(1 ns) range can be obtained from REMD with exchange times of 5 ps, in excellent agreement with results from long equilibrium molecular dynamics., Comment: 4 pages, 4 figures; v2: updated abstract and Fig. 3

We use computer simulations to study the kinetics and mechanism of proton passage through a narrow-pore carbon-nanotube membrane separating reservoirs of liquid water. Free energy and rate constant calculations show that protons move across the membrane diffusively in single-file chains of hydrogen-bonded water molecules. Proton passage through the membrane is opposed by a high barrier along the effective potential, reflecting the large electrostatic penalty for desolvation and reminiscent of charge exclusion in biological water channels. At neutral pH, we estimate a translocation rate of about 1 proton per hour and tube., Comment: 4 pages, 4 figures

Using a Coarse-grained Molecular Dynamics (CMD) approach we study the apparent nonlinear dynamics of water molecules filling/emptying carbon nanotubes as a function of system parameters. Different levels of the pore hydrophobicity give rise to tubes that are empty, water-filled, or fluctuate between these two long-lived metastable states. The corresponding coarse-grained free energy surfaces and their hysteretic parameter dependence are explored by linking MD to continuum fixed point and bifurcation algorithms. The results are validated through equilibrium MD simulations., Comment: 4 pages, 3 figures; accepted version

We present a ``coarse molecular dynamics'' approach and apply it to studying the kinetics and thermodynamics of a peptide fragment dissolved in water. Short bursts of appropriately initialized simulations are used to infer the deterministic and stochastic components of the peptide motion parametrized by an appropriate set of coarse variables. Techniques from traditional numerical analysis (Newton-Raphson, coarse projective integration) are thus enabled; these techniques help analyze important features of the free-energy landscape (coarse transition states, eigenvalues and eigenvectors, transition rates, etc.). Reverse integration of (irreversible) expected coarse variables backward in time can assist escape from free energy minima and trace low-dimensional free energy surfaces. To illustrate the ``coarse molecular dynamics'' approach, we combine multiple short (0.5-ps) replica simulations to map the free energy surface of the ``alanine dipeptide'' in water, and to determine the ~ 1/(1000 ps) rate of interconversion between the two stable configurational basins corresponding to the alpha-helical and extended minima., Comment: The article has been submitted to "The Journal of Chemical Physics."

A theoretical treatment based upon the quasi-chemical theory of solutions predicts the most probable number of water neighbors in the inner shell of a Li+ ion in liquid water to be four. The instability of a six water molecule inner sphere complex relative to four-coordinated structures is confirmed by an `ab initio' molecular dynamics calculation. A classical Monte Carlo simulation equilibrated 26 water molecules with a rigid six-coordinated Li(H2O)6+ complex with periodic boundary conditions in aqueous solution. With that initial configuration for the molecular dynamics, the six-coordinated structure relaxed into four-coordinated arrangements within 112 fs and stabilized. This conclusion differs from prior interpretations of neutron and X-ray scattering results on aqueous solutions., Comment: 9 pages, 2 figures

Ewald summation and physically equivalent methods such as particle-mesh Ewald, kubic-harmonic expansions, or Lekner sums are commonly used to calculate long-range electrostatic interactions in computer simulations of polar and charged substances. The calculation of pressures in such systems is investigated. We find that the virial and thermodynamic pressures differ because of the explicit volume dependence of the effective, resummed Ewald potential. The thermodynamic pressure, obtained from the volume derivative of the Helmholtz free energy, can be expressed easily for both ionic and rigid molecular systems. For a system of rigid molecules, the electrostatic energy and the forces at the atom positions are required, both of which are readily available in molecular dynamics codes. We then calculate the virial and thermodynamic pressures for the extended simple point charge (SPC/E) water model at standard conditions. We find that the thermodynamic pressure exhibits considerably less system size dependence than the virial pressure. From an analysis of the cross correlation between the virial and thermodynamic pressure, we conclude that the thermodynamic pressure should be used to drive volume fluctuations in constant-pressure simulations., Comment: RevTeX, 19 pages, 2 EPS figures; in press: Journal of Chemical Physics, 15-August-1998

Recent progress on molecular theories of hydration of nonpolar solutes in aqueous solution has led to new ways of thinking about the old issue of free, or "available," volume in liquids. This article surveys the principal new results with particular attention to general issues of packing in liquids., Comment: 14 pages included 3 figures; submitted to the proceedings of the NATO Advanced Study Institute: New Approaches to Old and New Problems in Liquid State Theory, Patti Marina (Messina) Italy, July 1998 replacement: minor typographical corrections

Recent developments in molecular theories and simulation of ions and polar molecules in water are reviewed. The hydration of imidazole and imidazolium solutes is used to exemplify the theoretical issues. The treatment of long-ranged electrostatic interactions in simulations is discussed extensively. It is argued that the Ewald approach is an easy way to get correct hydration free energies in the thermodynamic limit from molecular calculations; and that molecular simulations with Ewald interactions and periodic boundary conditions can also be more efficient than many common alternatives. The Ewald treatment permits a conclusive extrapolation to infinite system size. The picture that emerges from testing of simple models is that the most prominent failings of the simplest theories are associated with solvent proton conformations that lead to non-gaussian fluctuations of electrostatic potentials. Thus, the most favorable cases for the second-order perturbation theories are monoatomic positive ions. For polar and anionic solutes, continuum or gaussian theories are less accurate. The appreciation of the specific deficiencies of those simple models have led to new concepts, multistate gaussian and quasi-chemical theories, that address the cases for which the simpler theories fail. It is argued that, relative to direct dielectric continuum treatments, the quasi-chemical theories provide a better theoretical organization for the computational study of the electronic structure of solution species., Comment: invited J. Phys. Chem. feature article, 30 pages including 10 figures, minor stylistic revisions

Free energies of ionic solvation calculated from computer simulations exhibit a strong system size dependence. We perform a finite-size analysis based on a dielectric-continuum model with periodic boundary conditions. That analysis results in an estimate of the Born ion size. Remarkably, the finite-size correction applies to systems with only eight water molecules hydrating a sodium ion and results in an estimate of the Born radius of sodium that agrees with the experimental value., Comment: 2 EPS figures

Theories of solvation free energies often involve electrostatic potentials at the position of a solute charge. Simulation calculations that apply cutoffs and periodic boundary conditions based on molecular centers result in center-dependent contributions to electrostatic energies due to a systematic sorting of charges in radial shells. This sorting of charges induces a surface-charge density at the cutoff sphere or simulation-box boundary that depends on the choice of molecular centers. We identify a simple solution that gives correct, center-independent results, namely the radial integration of charge densities. Our conclusions are illustrated for a Lennard-Jones solute in water. The present results can affect the parameterization of force fields., Comment: 10 pages, 3 Encapsulated Postscript figures, LaTeX, uses RevTeX macro

An information theory model is used to construct a molecular explanation why hydrophobic solvation entropies measured in calorimetry of protein unfolding converge at a common temperature. The entropy convergence follows from the weak temperature dependence of occupancy fluctuations for molecular-scale volumes in water. The macroscopic expression of the contrasting entropic behavior between water and common organic solvents is the relative temperature insensitivity of the water isothermal compressibility. The information theory model provides a quantitative description of small molecule hydration and predicts a negative entropy at convergence. Interpretations of entropic contributions to protein folding should account for this result., Comment: Phys. Rev. Letts. (in press 1996), 3 pages, 3 figures

A closed form of the electrostatic potential of a homogeneously charged cube is derived by integration. The exact result is compared with multipole expansions for the exterior and interior of the cube. The electrostatic potential of a homogeneously charged square in two-dimensional electrostatics is also determined., Comment: 5 pages, 1 Postscript figure, uses revtex.sty

We study the chemical potential of water as a function of charge based on perturbation theory. By calculating the electrostatic-energy fluctuations of two states (fully charged and uncharged) we are able to determine accurate values for the dependence of the chemical potential on charge. We find identical results for the chemical-potential difference of fully charged and uncharged water from overlapping-histogram and acceptance-ratio methods and by smoothly connecting the curves of direct exponential averages. Our results agree with those of Rick and Berne (J. Am. Chem. Soc., 1994, 116, 3949) with respect to both the chemical-potential difference and its dependence on the charge coupling parameter. We observe significant deviations from simple Gaussian-fluctuation statistics. The dependence on the coupling parameter is not quadratic, as would be inferred from linear continuum models of electrostatics., Comment: 10 pages, 4 Postscript figures, uses revtex.sty

The hydration free energies of ions exhibit an approximately quadratic dependence on the ionic charge, as predicted by the Born model. We analyze this behavior using second-order perturbation theory. This provides effective methods to calculating free energies from equilibrium computer simulations. The average and the fluctuation of the electrostatic potential at charge sites appear as the first coefficients in a Taylor expansion of the free energy of charging. Combining the data from different charge states allows calculation of free-energy profiles as a function of the ionic charge. The first two Taylor coefficients of the free-energy profiles can be computed accurately from equi- librium simulations; but they are affected by a strong system-size dependence. We apply corrections for these finite-size effects by using Ewald lattice sum- mation and adding the self-interactions consistently. Results are presented for a model ion with methane-like Lennard-Jones parameters in SPC water. We find two very closely quadratic regimes with different parameters for positive and negative ions. We also studied the hydration free energy of potassium, calcium, fluoride, chloride, and bromide ions. We find negative ions to be solvated more strongly compared to positive ions of equal size, in agreement with experimen- tal data. We ascribe this preference of negative ions to their strong interac- tions with water hydrogens, which can penetrate the ionic van der Waals shell without direct energetic penalty in the models used. We also find a positive electrostatic potential at the center of uncharged Lennard-Jones particles in water, which favors negative ions. Regarding finite-system-size effects, we show that even using only 16 water molecules it is possible to calculate accu- rately the hydration free energy of sodium if self-interactions are considered., Comment: 17 pages, 6 Postscript figures, uses revtex.sty

Ewald summation is widely used to calculate electrostatic interactions in computer simulations of condensed-matter systems. We present an analysis of the errors arising from truncating the infinite real- and Fourier-space lattice sums in the Ewald formulation. We derive an optimal choice for the Fourier-space cutoff given a screening parameter $\eta$. We find that the number of vectors in Fourier space required to achieve a given accuracy scales with $\eta^3$. The proposed method can be used to determine computationally efficient parameters for Ewald sums, to assess the quality of Ewald-sum implementations, and to compare different implementations., Comment: 6 pages, 3 figures (Encapsulated PostScript), LaTex

Recent molecular simulation and integral equation results alkali-halide ion pair potentials-of-mean-force in water are discussed. Dielectric model calculations are implemented to check that these models produce that characteristic structure of contact and solvent-separated minima for oppositely charged ions in water under physiological thermodynamic conditions. Comparison of the dielectric model results with the most current molecular level information indicates that the dielectric model does not, however, provide an accurate description of these potentials-of-mean-force. We note that linear dielectric models correspond to modelistic implementations of second-order thermodynamic perturbation theory for the excess chemical potential of a distinguished solute molecule. Therefore, the molecular theory corresponding to the dielectric models is second-order thermodynamic perturbation theory for that excess chemical potential. The second-order, or fluctuation, term raises a technical computational issue of treatment of long-ranged interactions similar to the one which arises in calculation of the dielectric constant of the solvent. It is contended that the most important step for further development of dielectric models would be a separate assessment of the first-order perturbative term (equivalently the {\it potential at zero charge} ) which vanishes in the dielectric models but is generally nonzero. Parameterization of radii and molecular volumes should then be based of the second-order perturbative term alone. Illustrative initial calculations are presented and discussed., Comment: 37 pages and 8 figures. LA-UR-93-4205