The flow of a compressible fluid with slip through a cylinder with an asymmetric local constriction has been considered both numerically, as well as analytically. For the numerical work, a particle-based method whose dynamics is governed by the multiparticle collision (MPC) rule has been used together with a generalized boundary condition that allows for slip at the wall. Since it is well known that an MPC system corresponds to an ideal gas and behaves like a compressible, viscous flow on average, an approximate analytical solution has been derived from the compressible Navier–Stokes equations of motion coupled to an ideal gas equation of state using the Karman–Pohlhausen method. The constriction is assumed to have a polynomial form, and the location of maximum constriction is varied throughout the constricted portion of the cylinder. Results for centerline densities and centerline velocities have been compared for various Reynolds numbers, Mach numbers, wall slip values and flow geometries. Full article

A good deal of molecular dynamics simulations aims at predicting and quantifying rare events, such as the folding of a protein or a phase transition. Simulating rare events is often prohibitive, especially if the equations of motion are high-dimensional, as is the case in molecular dynamics. Various algorithms have been proposed for efficiently computing mean first passage times, transition rates or reaction pathways. This article surveys and discusses recent developments in the field of rare event simulation and outlines a new approach that combines ideas from optimal control and statistical mechanics. The optimal control approach described in detail resembles the use of Jarzynski’s equality for free energy calculations, but with an optimized protocol that speeds up the sampling, while (theoretically) giving variance-free estimators of the rare events statistics. We illustrate the new approach with two numerical examples and discuss its relation to existing methods. Full article

Explicit or implicit expressions of potential energy surfaces (PES) represent the basis of our ability to simulate condensed matter systems, possibly understanding and sometimes predicting their properties by purely computational methods. The paper provides an outline of the major approaches currently used to approximate and represent PESs and contains a brief discussion of what still needs to be achieved. The paper also analyses the relative role of empirical and ab initio methods, which represents a crucial issue affecting the future of modeling in chemical physics and materials science. Full article

Quantum Monte Carlo methods are among the most accurate algorithms for predicting properties of general quantum systems. We briefly introduce ground state, path integral at finite temperature and coupled electron-ion Monte Carlo methods, their merits and limitations. We then discuss recent calculations using these methods for dense liquid hydrogen as it undergoes a molecular/atomic (metal/insulator) transition. We then discuss a procedure that can be used to assess electronic density functionals, which in turn can be used on a larger scale for first principles calculations and apply this technique to dense hydrogen and liquid water. Full article

Rare, but important, transition events between long-lived states are a key feature of many molecular systems. In many cases, the computation of rare event statistics by direct molecular dynamics (MD) simulations is infeasible, even on the most powerful computers, because of the immensely long simulation timescales needed. Recently, a technique for spatial discretization of the molecular state space designed to help overcome such problems, so-called Markov State Models (MSMs), has attracted a lot of attention. We review the theoretical background and algorithmic realization of MSMs and illustrate their use by some numerical examples. Furthermore, we introduce a novel approach to using MSMs for the efficient solution of optimal control problems that appear in applications where one desires to optimize molecular properties by means of external controls. Full article

Malliavin weight sampling (MWS) is a stochastic calculus technique for computing the derivatives of averaged system properties with respect to parameters in stochastic simulations, without perturbing the system’s dynamics. It applies to systems in or out of equilibrium, in steady state or time-dependent situations, and has applications in the calculation of response coefficients, parameter sensitivities and Jacobian matrices for gradient-based parameter optimisation algorithms. The implementation of MWS has been described in the specific contexts of kinetic Monte Carlo and Brownian dynamics simulation algorithms. Here, we present a general theoretical framework for deriving the appropriate MWS update rule for any stochastic simulation algorithm. We also provide pedagogical information on its practical implementation. Full article

Quantum time correlation functions are often the principal objects of interest in experimental investigations of the dynamics of quantum systems. For instance, transport properties, such as diffusion and reaction rate coefficients, can be obtained by integrating these functions. The evaluation of such correlation functions entails sampling from quantum equilibrium density operators and quantum time evolution of operators. For condensed phase and complex systems, where quantum dynamics is difficult to carry out, approximations must often be made to compute these functions. We present a general scheme for the computation of correlation functions, which preserves the full quantum equilibrium structure of the system and approximates the time evolution with quantum-classical Liouville dynamics. Several aspects of the scheme are discussed, including a practical and general approach to sample the quantum equilibrium density, the properties of the quantum-classical Liouville equation in the context of correlation function computations, simulation schemes for the approximate dynamics and their interpretation and connections to other approximate quantum dynamical methods. Full article

We review a selection of methods for performing enhanced sampling in molecular dynamics simulations. We consider methods based on collective variable biasing and on tempering, and offer both historical and contemporary perspectives. In collective-variable biasing, we first discuss methods stemming from thermodynamic integration that use mean force biasing, including the adaptive biasing force algorithm and temperature acceleration. We then turn to methods that use bias potentials, including umbrella sampling and metadynamics. We next consider parallel tempering and replica-exchange methods. We conclude with a brief presentation of some combination methods. Full article

This paper invites the reader to learn more about time integrators for Molecular Dynamics simulation through a simple MATLAB implementation. An overview of methods is provided from an algorithmic viewpoint that emphasizes long-time stability and finite-time dynamic accuracy. The given software simulates Langevin dynamics using an explicit, second-order (weakly) accurate integrator that exactly reproduces the Boltzmann-Gibbs density. This latter feature comes from adding a Metropolis acceptance-rejection step to the integrator. The paper discusses in detail the properties of the integrator. Since these properties do not rely on a specific form of a heat or pressure bath model, the given algorithm can be used to simulate other bath models including, e.g., the widely used v-rescale thermostat. Full article

We analyze the time reversible Born-Oppenheimer molecular dynamics (TRBOMD) scheme, which preserves the time reversibility of the Born-Oppenheimer molecular dynamics even with non-convergent self-consistent field iteration. In the linear response regime, we derive the stability condition, as well as the accuracy of TRBOMD for computing physical properties, such as the phonon frequency obtained from the molecular dynamics simulation. We connect and compare TRBOMD with Car-Parrinello molecular dynamics in terms of accuracy and stability. We further discuss the accuracy of TRBOMD beyond the linear response regime for non-equilibrium dynamics of nuclei. Our results are demonstrated through numerical experiments using a simplified one-dimensional model for Kohn-Sham density functional theory. Full article

Computing quantum dynamics in condensed matter systems is an open challenge due to the exponential scaling of exact algorithms with the number of degrees of freedom. Current methods try to reduce the cost of the calculation using classical dynamics as the key ingredient of approximations of the quantum time evolution. Two main approaches exist, quantum classical and semi-classical, but they suffer from various difficulties, in particular when trying to go beyond the classical approximation. It may then be useful to reconsider the problem focusing on statistical time-dependent averages rather than directly on the dynamics. In this paper, we discuss a recently developed scheme for calculating symmetrized correlation functions. In this scheme, the full (complex time) evolution is broken into segments alternating thermal and real-time propagation, and the latter is reduced to classical dynamics via a linearization approximation. Increasing the number of segments systematically improves the result with respect to full classical dynamics, but at a cost which is still prohibitive. If only one segment is considered, a cumulant expansion can be used to obtain a computationally efficient algorithm, which has proven accurate for condensed phase systems in moderately quantum regimes. This scheme is summarized in the second part of the paper. We conclude by outlining how the cumulant expansion formally provides a way to improve convergence also for more than one segment. Future work will focus on testing the numerical performance of this extension and, more importantly, on investigating the limit for the number of segments that goes to infinity of the approximate expression for the symmetrized correlation function to assess formally its convergence to the exact result. Full article

As shown by Jarzynski, free energy differences between equilibrium states can be expressed in terms of the statistics of work carried out on a system during non-equilibrium transformations. This exact result, as well as the related Crooks fluctuation theorem, provide the basis for the computation of free energy differences from fast switching molecular dynamics simulations, in which an external parameter is changed at a finite rate, driving the system away from equilibrium. In this article, we first briefly review the Jarzynski identity and the Crooks fluctuation theorem and then survey various algorithms building on these relations. We pay particular attention to the statistical efficiency of these methods and discuss practical issues arising in their implementation and the analysis of the results. Full article

Although substantial progress has been made in recent years in research onsheared granular matter, relatively few studies concentrate on the behavior of materials withvery strong polydispersity. In this paper, shear deformation of a two-dimensional granularmaterial composed of frictional disk-shaped grains with power-law size distribution isanalyzed numerically with a finite-difference model. The analysis of the results concentrateson those aspects of the behavior of the modeled system that are related to its polydispersity. Itis demonstrated that many important global material properties are dependent on the behaviorof the largest grains from the tail of the size distribution. In particular, they are responsiblefor global correlation of velocity anomalies emerging at the jamming transition. They alsobuild a skeleton of the global contact and force networks in shear-jammed systems, leadingto the very open, “sparse” structure of those networks, consisting of only ~ 35% of all grains.The details of the model are formulated so that it represents fragmented sea ice moving ona two-dimensional sea surface; however, the results are relevant for other types of stronglypolydisperse granular materials, as well. Full article

Coarse-grained models of soft matter are usually combined with implicit solvent models that take the electrostatic polarizability into account via a dielectric background. In biophysical or nanoscale simulations that include water, this constant can vary greatly within the system. Performing molecular dynamics or other simulations that need to compute exact electrostatic interactions between charges in those systems is computationally demanding. We review here several algorithms developed by us that perform exactly this task. For planar dielectric surfaces in partial periodic boundary conditions, the arising image charges can be either treated with the MMM2D algorithm in a very efficient and accurate way or with the electrostatic layer correction term, which enables the user to use his favorite 3D periodic Coulomb solver. Arbitrarily-shaped interfaces can be dealt with using induced surface charges with the induced charge calculation (ICC*) algorithm. Finally, the local electrostatics algorithm, MEMD(Maxwell Equations Molecular Dynamics), even allows one to employ a smoothly varying dielectric constant in the systems. We introduce the concepts of these three algorithms and an extension for the inclusion of boundaries that are to be held fixed at a constant potential (metal conditions). For each method, we present a showcase application to highlight the importance of dielectric interfaces. Full article

Ab initio molecular dynamics simulations were performed to investigate the elasticity of cubic CaSiO₃ perovskite at high pressure and temperature. All three independent elastic constants for cubic CaSiO₃ perovskite, C₁₁, C₁₂, and C₄₄, were calculated from the computation of stress generated by small strains. The elastic constants were used to estimate the moduli and seismic wave velocities at the high pressure and high temperature characteristic of the Earth’s interior. The dependence of temperature for sound wave velocities decreased as the pressure increased. There was little difference between the estimated compressional sound wave velocity (V_P) in cubic CaSiO₃ perovskite and that in the Earth’s mantle, determined by seismological data. By contrast, a significant difference between the estimated shear sound wave velocity (V_S) and that in the Earth’s mantle was confirmed. The elastic properties of cubic CaSiO₃ perovskite cannot explain the properties of the Earth’s lower mantle, indicating that the cubic CaSiO₃ perovskite phase is a minor mineral in the Earth’s lower mantle. Full article

A methodology for the determination of the solid-fluid contact angle, to be employed within molecular dynamics (MD) simulations, is developed and systematically applied. The calculation of the contact angle of a fluid drop on a given surface, averaged over an equilibrated MD trajectory, is divided in three main steps: (i) the determination of the fluid molecules that constitute the interface, (ii) the treatment of the interfacial molecules as a point cloud data set to define a geometric surface, using surface meshing techniques to compute the surface normals from the mesh, (iii) the collection and averaging of the interface normals collected from the post-processing of the MD trajectory. The average vector thus found is used to calculate the Cassie contact angle (i.e., the arccosine of the averaged normal z-component). As an example we explore the effect of the size of a drop of water on the observed solid-fluid contact angle. A single coarse-grained bead representing two water molecules and parameterized using the SAFT-γ Mie equation of state (EoS) is employed, meanwhile the solid surfaces are mimicked using integrated potentials. The contact angle is seen to be a strong function of the system size for small nano-droplets. The thermodynamic limit, corresponding to the infinite size (macroscopic) drop is only truly recovered when using an excess of half a million water coarse-grained beads and/or a drop radius of over 26 nm. Full article

Molecular modeling is an important subdomain in the field of computational modeling, regarding both scientific and industrial applications. This is because computer simulations on a molecular level are a virtuous instrument to study the impact of microscopic on macroscopic phenomena. Accurate molecular models are indispensable for such simulations in order to predict physical target observables, like density, pressure, diffusion coefficients or energetic properties, quantitatively over a wide range of temperatures. Thereby, molecular interactions are described mathematically by force fields. The mathematical description includes parameters for both intramolecular and intermolecular interactions. While intramolecular force field parameters can be determined by quantum mechanics, the parameterization of the intermolecular part is often tedious. Recently, an empirical procedure, based on the minimization of a loss function between simulated and experimental physical properties, was published by the authors. Thereby, efficient gradient-based numerical optimization algorithms were used. However, empirical force field optimization is inhibited by the two following central issues appearing in molecular simulations: firstly, they are extremely time-consuming, even on modern and high-performance computer clusters, and secondly, simulation data is affected by statistical noise. The latter provokes the fact that an accurate computation of gradients or Hessians is nearly impossible close to a local or global minimum, mainly because the loss function is flat. Therefore, the question arises of whether to apply a derivative-free method approximating the loss function by an appropriate model function. In this paper, a new Sparse Grid-based Optimization Workflow (SpaGrOW) is presented, which accomplishes this task robustly and, at the same time, keeps the number of time-consuming simulations relatively small. This is achieved by an efficient sampling procedure for the approximation based on sparse grids, which is described in full detail: in order to counteract the fact that sparse grids are fully occupied on their boundaries, a mathematical transformation is applied to generate homogeneous Dirichlet boundary conditions. As the main drawback of sparse grids methods is the assumption that the function to be modeled exhibits certain smoothness properties, it has to be approximated by smooth functions first. Radial basis functions turned out to be very suitable to solve this task. The smoothing procedure and the subsequent interpolation on sparse grids are performed within sufficiently large compact trust regions of the parameter space. It is shown and explained how the combination of the three ingredients leads to a new efficient derivative-free algorithm, which has the additional advantage that it is capable of reducing the overall number of simulations by a factor of about two in comparison to gradient-based optimization methods. At the same time, the robustness with respect to statistical noise is maintained. This assertion is proven by both theoretical considerations and practical evaluations for molecular simulations on chemical example substances. Full article

A reduction of the cost for long-range interaction calculation is essential for large-scale molecular systems that contain a lot of point charges. Cutoff methods are often used to reduce the cost of long-range interaction calculations. Molecular dynamics (MD) simulations can be accelerated by using cutoff methods; however, simple truncation or approximation of long-range interactions often offers serious defects for various systems. For example, thermodynamical properties of polar molecular systems are strongly affected by the treatment of the Coulombic interactions and may lead to unphysical results. To assess the truncation effect of some cutoff methods that are categorized as the shift function method, MD simulations for bulk water systems were performed. The results reflect two main factors, i.e., the treatment of cutoff boundary conditions and the presence/absence of the theoretical background for the long-range approximation. Full article

Molecular dynamics simulations are used to study the evaporation of water droplets containing either dissolved LiCl, NaCl or KCl salt in a gaseous surrounding (nitrogen) with a constant high temperature of 600 K. The initial droplet has 298 K temperature and contains 1,120 water molecules, 0, 40, 80 or 120 salt molecules. The effects of the salt type and concentration on the evaporation rate are examined. Three stages with different evaporation rates are observed for all cases. In the initial stage of evaporation, the droplet evaporates slowly due to low droplet temperature and high evaporation latent heat for water, and pure water and aqueous solution have almost the same evaporation rates. In the second stage, evaporation rate is increased significantly, and evaporation is somewhat slower for the aqueous salt-containing droplet than the pure water droplet due to the attracted ion-water interaction and hydration effect. The Li⁺-water has the strongest interaction and hydration effect, so LiCl aqueous droplets evaporate the slowest, then NaCl and KCl. Higher salt concentration also enhances the ion-water interaction and hydration effect, and hence corresponds to a slower evaporation. In the last stage of evaporation, only a small amount of water molecules are left in the droplet, leading to a significant increase in ion-water interactions, so that the evaporation becomes slower compared to that in the second stage. Full article

In the last few decades, computer simulations have become a fundamental tool in the field of soft matter science, allowing researchers to investigate the properties of a large variety of systems. Nonetheless, even the most powerful computational resources presently available are, in general, sufficient to simulate complex biomolecules only for a few nanoseconds. This limitation is often circumvented by using coarse-grained models, in which only a subset of the system’s degrees of freedom is retained; for an effective and insightful use of these simplified models; however, an appropriate parametrization of the interactions is of fundamental importance. Additionally, in many cases the removal of fine-grained details in a specific, small region of the system would destroy relevant features; such cases can be treated using dual-resolution simulation methods, where a subregion of the system is described with high resolution, and a coarse-grained representation is employed in the rest of the simulation domain. In this review we discuss the basic notions of coarse-graining theory, presenting the most common methodologies employed to build low-resolution descriptions of a system and putting particular emphasis on their similarities and differences. The AdResS and H-AdResS adaptive resolution simulation schemes are reported as examples of dual-resolution approaches, especially focusing in particular on their theoretical background. Full article

In this review, we discuss the Dynamical approach to Non-Equilibrium Molecular Dynamics (D-NEMD), which extends stationary NEMD to time-dependent situations, be they responses or relaxations. Based on the original Onsager regression hypothesis, implemented in the nineteen-seventies by Ciccotti, Jacucci and MacDonald, the approach permits one to separate the problem of dynamical evolution from the problem of sampling the initial condition. D-NEMD provides the theoretical framework to compute time-dependent macroscopic dynamical behaviors by averaging on a large sample of non-equilibrium trajectories starting from an ensemble of initial conditions generated from a suitable (equilibrium or non-equilibrium) distribution at time zero. We also discuss how to generate a large class of initial distributions. The same approach applies also to the calculation of the rate constants of activated processes. The range of problems treatable by this method is illustrated by discussing applications to a few key hydrodynamic processes (the “classical” flow under shear, the formation of convective cells and the relaxation of an interface between two immiscible liquids). Full article

Performing molecular dynamics in electronically excited states requires the inclusion of nonadiabatic effects to properly describe phenomena beyond the Born-Oppenheimer approximation. This article provides a survey of selected nonadiabatic methods based on quantum or classical trajectories. Among these techniques, trajectory surface hopping constitutes an interesting compromise between accuracy and efficiency for the simulation of medium- to large-scale molecular systems. This approach is, however, based on non-rigorous approximations that could compromise, in some cases, the correct description of the nonadiabatic effects under consideration and hamper a systematic improvement of the theory. With the help of an in principle exact description of nonadiabatic dynamics based on Bohmian quantum trajectories, we will investigate the origin of the main approximations in trajectory surface hopping and illustrate some of the limits of this approach by means of a few simple examples. Full article

Since most experimental observations are performed at constant temperature and pressure, the isothermal-isobaric (NPT) ensemble has been widely used in molecular simulations. Nevertheless, the NPT ensemble has only recently been placed on a rigorous foundation. The proper formulation of the NPT ensemble requires a “shell” particle to uniquely identify the volume of the system, thereby avoiding the redundant counting of configurations. Here, we review our recent work in incorporating a shell particle into molecular dynamics simulation algorithms to generate the correct NPT ensemble averages. Unlike previous methods, a piston of unknown mass is no longer needed to control the response time of the volume fluctuations. As the volume of the system is attached to the shell particle, the system itself now sets the time scales for volume and pressure fluctuations. Finally, we discuss a number of tests that ensure the equations of motion sample phase space correctly and consider the response time of the system to pressure changes with and without the shell particle. Overall, the shell particle algorithm is an effective simulation method for studying systems exposed to a constant external pressure and may provide an advantage over other existing constant pressure approaches when developing nonequilibrium molecular dynamics methods. Full article

In this work, we make an attempt to answer the question of what a multiscale problem is in Molecular Dynamics (MD), or, more in general, in Molecular Simulation (MS). By introducing the criterion of separability of scales, we identify three major (reference) categories of multiscale problems and discuss their corresponding computational strategies by making explicit examples of applications. Full article