Effect of coarse graining in water models for the study of kinetics and mechanisms of clathrate hydrates nucleation and growth.

Clathrate hydrates are crystalline inclusion compounds wherein a water framework encages small guest atoms/molecules within its cavities. Among the others, methane clathrates are the largest fossil fuel resource still available. They can also be used to safely transport gases and can also form spontaneously under suitable conditions plugging pipelines. Understanding the crystallization mechanism is very important, and given the impossibility of experimentally identifying the atomistic path, simulations played an important role in this field. Given the large computational cost of these simulations, in addition to all-atom force fields, scientists considered coarse-grained water models. Here, we have investigated the effect of coarse-graining, as implemented in the water model mW, on the crystallization characteristics of methane clathrate in comparison with the all-atom TIP4P force field. Our analyses revealed that although the characteristics directly depending on the energetics of the water models are well reproduced, dynamical properties are off by the orders of magnitude. Being crystallization a non-equilibrium process, the altered kinetics of the process results in different characteristics of crystalline nuclei. Both TIP4P and mW water models produce methane clathrate nuclei with some amount of the less stable (in the given thermodynamic conditions) structure II phase and an excess of pentagonal dodecahedral cages over the tetrakaidecahedral ones regarding the ideal ratio in structure I. However, the dependence of this excess on the methane concentration in solution is higher with the former water model, whereas with the latter, the methane concentration in solution dependence is reduced and within the statistical error.

Jarzynski equality on work and free energy: Crystal indentation as a case study.

Mathematical relations concerning particle systems require knowledge of the applicability conditions to become physically relevant and not merely formal. We illustrate this fact through the analysis of the Jarzynski equality (JE), whose derivation for Hamiltonian systems suggests that the equilibrium free-energy variations can be computational or experimentally determined in almost any kind of non-equilibrium processes. This apparent generality is surprising in a mechanical theory. Analytically, we show that the quantity called "work" in the Hamiltonian derivation of the JE is neither a thermodynamic quantity nor mechanical work, except in special circumstances to be singularly assessed. Through molecular dynamics simulations of elastic and plastic deformations induced via nano-indentation of crystalline surfaces that fall within the formal framework of the JE, we illustrate that the JE cannot be verified and that the results of this verification are process dependent.

Dynamical nonequilibrium molecular dynamics reveals the structural basis for allostery and signal propagation in biomolecular systems.

ABSTRACT: A dynamical approach to nonequilibrium molecular dynamics (D-NEMD), proposed in the 1970s by Ciccotti et al., is undergoing a renaissance and is having increasing impact in the study of biological macromolecules. This D-NEMD approach, combining MD simulations in stationary (in particular, equilibrium) and nonequilibrium conditions, allows for the determination of the time-dependent structural response of a system using the Kubo-Onsager relation. Besides providing a detailed picture of the system's dynamic structural response to an external perturbation, this approach also has the advantage that the statistical significance of the response can be assessed. The D-NEMD approach has been used recently to identify a general mechanism of inter-domain signal propagation in nicotinic acetylcholine receptors, and allosteric effects in ß -lactamase enzymes, for example. It complements equilibrium MD and is a very promising approach to identifying and analysing allosteric effects. Here, we review the D-NEMD approach and its application to biomolecular systems, including transporters, receptors, and enzymes.

Quantum Trajectories for the Dynamics in the Exact Factorization Framework: A Proof-of-Principle Test.

In the framework of the exact factorization of the time-dependent electron-nuclear wave function, we investigate the possibility of solving the nuclear time-dependent Schrödinger equation based on trajectories. The nuclear equation is separated in a Hamilton-Jacobi equation for the phase of the wave function, and a continuity equation for its (squared) modulus. For illustrative adiabatic and nonadiabatic one-dimensional models, we implement a procedure to follow the evolution of the nuclear density along the characteristics of the Hamilton-Jacobi equation. Those characteristics are referred to as quantum trajectories, since they are generated via ordinary differential equations similar to Hamilton's equations, but including the so-called quantum potential, and they can be used to reconstruct exactly the quantum-mechanical nuclear wave function, provided infinite initial conditions are propagated in time.

Adiabatic motion and statistical mechanics via mass-zero constrained dynamics.

In recent work [Coretti et al., J. Chem. Phys., 2018, 149, 191102], a new algorithm to solve numerically the dynamics of the shell model for polarization was presented. The approach, broadly applicable to systems involving adiabatically separated dynamical variables, employs constrained molecular dynamics to strictly enforce the condition that the force on the fast degrees of freedom, modeled as having zero mass, is null at each time step. The algorithm is symplectic and fully time reversible, and results in stable and efficient propagation. In this paper we complete the discussion of the mechanics of mass-zero constrained dynamics by showing how to adapt it to problems where the fast degrees of freedom must satisfy additional conditions. This extension includes, in particular, the important case of first principles molecular dynamics. We then consider the statistical mechanics of the mass-zero constrained dynamical system demonstrating that the marginal probability sampled by the dynamics in the physical phase space recovers the form of the Born-Oppenheimer probability density. The effectiveness of the approach and the favorable scaling of the algorithm with system size are illustrated in test calculations of solid Na via orbital-free density functional dynamics.

A General Mechanism for Signal Propagation in the Nicotinic Acetylcholine Receptor Family.

Nicotinic acetylcholine receptors (nAChRs) modulate synaptic activity in the central nervous system. The α7 subtype, in particular, has attracted considerable interest in drug discovery as a target for several conditions, including Alzheimer's disease and schizophrenia. Identifying agonist-induced structural changes underlying nAChR activation is fundamentally important for understanding biological function and rational drug design. Here, extensive equilibrium and nonequilibrium molecular dynamics simulations, enabled by cloud-based high-performance computing, reveal the molecular mechanism by which structural changes induced by agonist unbinding are transmitted within the human α7 nAChR. The simulations reveal the sequence of coupled structural changes involved in driving conformational change responsible for biological function. Comparison with simulations of the α4ß2 nAChR subtype identifies features of the dynamical architecture common to both receptors, suggesting a general structural mechanism for signal propagation in this important family of receptors.

Ab initio accelerated molecular dynamics study of the hydride ligands in the ruthenium complex: Ru(H2)2H2(P(C5H9)3)2.

The dihydrogen complex Ru(H2)2H2(P(C5H9)3)2 has been investigated, via ab initio accelerated molecular dynamics, to elucidate the H ligands dynamics and possible reaction paths for H2/H exchange. We have characterized the free energy landscape associated with the H atoms positional exchange around the Ru centre. From the free energy landscape, we have been able to estimate a barrier of 6 kcal mol-1 for the H2/H exchange process. We have also observed a trihydrogen intermediate as a passing state along some of the possible reaction pathways.

On computing the solubility of molecular systems subject to constraints using the extended Einstein crystal method.

A method to compute solubilities for molecular systems using atomistic simulations, based on an extension of the Einstein crystal method, has recently been presented [Li et al., J. Chem. Phys. 146, 214110 (2017)]. This methodology is particularly appealing to compute solubilities in cases of practical importance including, but not limited to, solutions where the solute is sparingly soluble and molecules of importance for the pharmaceutical industry, which are often characterized by strong polar interactions and slow relaxation time scales. The mathematical derivation of this methodology hinges on a factorization of the partition function which is not necessarily applicable in the case of a system subject to holonomic molecular constraints. We show here that, although the mathematical procedure to derive it is slightly different, essentially the same mathematical relation for calculating the solubility can be safely applied for computing the solubility of systems subject to constraints, which are the majority of the systems used for practical molecular simulations.

Solubility of paracetamol in ethanol by molecular dynamics using the extended Einstein crystal method and experiments.

Li and co-workers [Li et al., J. Chem. Phys. 146, 214110 (2017)] have recently proposed a methodology to compute the solubility of molecular compounds from first principles, using molecular dynamics simulations. We revise and further explore their methodology that was originally applied to naphthalene in water at low concentration. In particular, we compute the solubility of paracetamol in an ethanol solution at ambient conditions. For the simulations, we used a force field that we previously reparameterized to reproduce certain thermodynamic properties of paracetamol but not explicitly its solubility in ethanol. In addition, we have determined the experimental solubility by performing turbidity measurements using a Crystal16 over a range of temperatures. Our work serves a dual purpose: (i) methodologically, we clarify how to compute, with a relatively straightforward procedure, the solubility of molecular compounds and (ii) applying this procedure, we show that the solubility predicted by our force field (0.085 ± 0.014 in mole ratio) is in good agreement with the experimental value obtained from our experiments and those reported in the literature (average 0.0585 ± 0.004), considering typical deviations for predictions from first principle methods. The good agreement between the experimental and the calculated solubility also suggests that the method used to reparameterize the force field can be used as a general strategy to optimize force fields for simulations in solution.

The physics of open systems for the simulation of complex molecular environments in soft matter.

Molecular dynamics (MD) has become one of the most powerful tools of investigation in soft matter. Despite such success, simulations of large molecular environments are mostly run using the approximation of closed systems without the possibility of exchange of matter. Due to the molecular complexity of soft matter systems, an optimal simulation strategy would require the application of concurrent multiscale resolution approaches such that each part of a large system can be considered as an open subsystem at a high resolution embedded in a large coarser reservoir of energy and particles. This paper discusses the current capability and the future perspectives of multiscale adaptive resolution MD methods to satisfy the conceptual principles of open systems and to perform simulations of complex molecular environments in soft matter.

Communication: Constrained molecular dynamics for polarizable models.

A new algorithm to solve numerically the evolution of empirical shell models of polarizable systems is presented. It employs constrained molecular dynamics to satisfy exactly, at each time step, the crucial condition that the gradient of the potential with respect to the shell degrees of freedom is null. The algorithm is efficient, stable, and, contrary to the available alternatives, it is symplectic and time reversible. A proof-of-principle calculation on a polarizable model for NaCl is presented to illustrate its properties in comparison with the current method, which employs a conjugate-gradient procedure to enforce the null gradient condition. The proposed algorithm is applicable to other cases where a minimum condition on a function of an auxiliary set of driven dynamical variables must be satisfied.

Nucleation of Molecular Crystals Driven by Relative Information Entropy.

Simulating nucleation of molecular crystals is extremely challenging for all but the simplest cases. The challenge lies in formulating effective order parameters that are capable of driving the transition process. In recent years, order parameters based on molecular pair-functions have been successfully used in combination with enhanced sampling techniques to simulate nucleation of simple molecular crystals. However, despite the success of these approaches, we demonstrate that they can fail when applied to more complex cases. In fact, we show that order parameters based on molecular pair-functions, while successful at nucleating benzene, fail for paracetamol. Hence, we introduce a novel approach to formulate order parameters. In our approach, we construct reduced dimensional distributions of relevant quantities on the fly and then quantify the difference between these distributions and selected reference distributions. By computing the distribution of different quantities and by choosing different reference distributions, it is possible to systematically construct an effective set of order parameters. We then show that our new order parameters are capable of driving the nucleation of ordered states and, in particular, the form I crystal of paracetamol.

Thermal Diffusion in Binary Mixtures: Transient Behavior and Transport Coefficients from Equilibrium and Nonequilibrium Molecular Dynamics.

Equilibrium and nonequilibrium molecular dynamics simulations are combined to compute the full set of coefficients that appear in the phenomenological equations describing thermal transport in a binary mixture subject to a constant thermal gradient. The Dynamical Non-Equilibrium Molecular Dynamics approach (D-NEMD) is employed to obtain the microscopic time evolution of the density and temperature fields, together with that of the mass and energy fluxes. D-NEMD enables one to study not only the steady state, but also the evolution of the fields during the transient that follows the onset of the thermal gradient, up to the establishment of the steady state. This makes it possible to ensure that the system has indeed reached a stationary condition, and to analyze the transient mechanisms and time scales of the mass and energy transport. A local time averaging procedure is applied to each trajectory contributing to the calculation to improve the signal-to-noise ratio in the estimation of the fluxes and to obtain a clear signal with the, relatively limited, statistics available.

Time-reversal symmetry for systems in a constant external magnetic field.

The time-reversal properties of charged systems in a constant external magnetic field are reconsidered in this paper. We show that the evolution equations of the system are invariant under a new symmetry operation that implies a new signature property for time-correlation functions under time reversal. We then show how these findings can be combined with a previously identified symmetry to determine, for example, null components of the correlation functions of velocities and currents and of the associated transport coefficients. These theoretical predictions are illustrated by molecular dynamics simulations of superionic AgI.

On metastability and Markov state models for non-stationary molecular dynamics.

Unlike for systems in equilibrium, a straightforward definition of a metastable set in the non-stationary, non-equilibrium case may only be given case-by-case-and therefore it is not directly useful any more, in particular in cases where the slowest relaxation time scales are comparable to the time scales at which the external field driving the system varies. We generalize the concept of metastability by relying on the theory of coherent sets. A pair of sets A and B is called coherent with respect to the time interval [t1, t2] if (a) most of the trajectories starting in A at t1 end up in B at t2 and (b) most of the trajectories arriving in B at t2 actually started from A at t1. Based on this definition, we can show how to compute coherent sets and then derive finite-time non-stationary Markov state models. We illustrate this concept and its main differences to equilibrium Markov state modeling on simple, one-dimensional examples.

On the properties of a bundle of flexible actin filaments in an optical trap.

We establish the statistical mechanics framework for a bundle of Nf living and uncrosslinked actin filaments in a supercritical solution of free monomers pressing against a mobile wall. The filaments are anchored normally to a fixed planar surface at one of their ends and, because of their limited flexibility, they grow almost parallel to each other. Their growing ends hit a moving obstacle, depicted as a second planar wall, parallel to the previous one and subjected to a harmonic compressive force. The force constant is denoted as the trap strength while the distance between the two walls as the trap length to make contact with the experimental optical trap apparatus. For an ideal solution of reactive filaments and free monomers at fixed free monomer chemical potential µ1, we obtain the general expression for the grand potential from which we derive averages and distributions of relevant physical quantities, namely, the obstacle position, the bundle polymerization force, and the number of filaments in direct contact with the wall. The grafted living filaments are modeled as discrete Wormlike chains, with F-actin persistence length ℓp, subject to discrete contour length variations ±d (the monomer size) to model single monomer (de)polymerization steps. Rigid filaments (ℓp = ∞), either isolated or in bundles, all provide average values of the stalling force in agreement with Hill's predictions Fs (H)=NfkBTln(ρ1/ρ1c)/d, independent of the average trap length. Here ρ1 is the density of free monomers in the solution and ρ1c its critical value at which the filament does not grow nor shrink in the absence of external forces. Flexible filaments (ℓp < ∞) instead, for values of the trap strength suitable to prevent their lateral escape, provide an average bundle force and an average trap length slightly larger than the corresponding rigid cases (few percents). Still the stalling force remains nearly independent on the average trap length, but results from the product of two strongly L-dependent contributions: the fraction of touching filaments ∝〈L〉(O.T.) (2) and the single filament buckling force ∝〈L〉(O.T.) (-2).

Spontaneous chiral symmetry breaking in collective active motion.

Chiral symmetry breaking is ubiquitous in biological systems, from DNA to bacterial suspensions. A key unresolved problem is how chiral structures may spontaneously emerge from achiral interactions. We study a simple model of active swimmers in three dimensions that effectively incorporates hydrodynamic interactions. We perform large-scale molecular dynamics simulations (up to 10(6) particles) and find long-lived metastable collective states that exhibit chiral organization although the interactions are achiral. We elucidate under which conditions these chiral states will emerge and grow to large scales. To explore the complex phase space available to the system, we perform nonequilibrium quenches on a one-dimensional Lebwohl-Lasher model with periodic boundary conditions to study the likelihood of formation of chiral structures.

Probabilistic Derivation of Spatiotemporal Correlation Functions in the Hydrodynamic Limit.

In this paper, we use probability theory to prove in suitable conditions the equivalence of equilibrium time correlation functions of microscopic density fields with the time correlation functions of local macroscopic density fields evolved by hydrodynamics in (approximate) phenomenological continuum theories of matter. We further discuss a useful and rigorous numerical algorithm, derived from this framework, to compute macroscopic space- and time-dependent behaviors (such as the hydrodynamical one) via molecular dynamics simulations.

Particle-Based Modeling of Living Actin Filaments in an Optical Trap.

We report a coarse-grained molecular dynamics simulation study of a bundle of parallel actin filaments under supercritical conditions pressing against a loaded mobile wall using a particle-based approach where each particle represents an actin unit. The filaments are grafted to a fixed wall at one end and are reactive at the other end, where they can perform single monomer (de)polymerization steps and push on a mobile obstacle. We simulate a reactive grand canonical ensemble in a box of fixed transverse area A, with a fixed number of grafted filaments N f , at temperature T and monomer chemical potential µ 1 . For a single filament case ( N f = 1 ) and for a bundle of N f = 8 filaments, we analyze the structural and dynamical properties at equilibrium where the external load compensates the average force exerted by the bundle. The dynamics of the bundle-moving-wall unit are characteristic of an over-damped Brownian oscillator in agreement with recent in vitro experiments by an optical trap setup. We analyze the influence of the pressing wall on the kinetic rates of (de)polymerization events for the filaments. Both static and dynamic results compare reasonably well with recent theoretical treatments of the same system. Thus, we consider the proposed model as a good tool to investigate the properties of a bundle of living filaments.