Dynamical Non-Equilibrium Molecular Dynamics

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).


Introduction
The most widespread use of Molecular Dynamics (MD) [1,2], in the same spirit of Monte Carlo (MC) [3,4], is to compute the thermodynamic or statistical behavior of molecular systems at equilibrium.This means that, starting from the assumption of the validity of the ergodic hypothesis, dynamical (MD) or fictitious-time (MC) trajectories are used to sample the equilibrium distribution in phase space (MD) or in configurational space (MC)."Time" averages over the generated trajectories will thereafter provide the statistical properties of the system.
At variance with Monte Carlo, the dynamical approach of Molecular Dynamics can be directly extended to sample distributions corresponding to stationary non-equilibrium conditions, where there exists a stationary distribution but, at variance with equilibrium, its expression is not explicitly known.However, the statistical problem of sampling a time-dependent ensemble cannot be solved by generating states along a single dynamical non-equilibrium trajectory, as long as time cannot be taken as homogeneous and averages over time make no sense.
Generally, to compute macroscopic dynamical behaviors, as, e.g., in hydrodynamics, the assumption of time-scale separation is made and rigorous ensemble averages are substituted with short-time averages equivalent to local smoothing.This may not be the case, sometimes.Moreover, the statistical error implied by this procedure cannot be made as small as desirable and possible.These difficulties can be faced and solved.
In the nineteen-thirties, Lars Onsager [5] observed that an induced (non-equilibrium) relaxation towards equilibrium could be obtained by studying the regression of the corresponding spontaneous fluctuations at equilibrium.Later, in the nineteen-fifties, Kubo [6] provided a mathematical formulation of Onsager's ideas by showing how the (linear) response of a system, initially at equilibrium, to a time-dependent (external) physical perturbation could be obtained by convoluting it with an appropriate equilibrium time-correlation function [7][8][9].Kubo also derived the formal expression for the complete (linear and nonlinear) response.
In the case of Kubo's procedure one does not need to make reference to an initial equilibrium state, but can, rather, refer to an arbitrary initial distribution at time t 0 = 0 of the system.This result has an important consequence for Molecular Dynamics simulations, since it allows one to separate the problem of dynamical evolution from the problem of sampling the initial condition.
Starting from the mid-nineteen-seventies, the direct numerical simulation of the response was used in conjunction with a sample of initial conditions extracted from an equilibrium trajectory [10,11].In this context, the problem of achieving a reasonable signal-to-noise ratio, even for weak perturbations, was solved for short times by introducing the so-called subtraction technique [12], which permitted one to verify, with surprising results [13], the range of the validity of linearity.Some time later on, it was realized that the same approach could be used to calculate dynamical properties for rare events (e.g., transmission coefficients) by averaging the dynamical response over time-dependent trajectories started from initial conditions sampled from a constrained/conditional equilibrium ensemble [14][15][16][17][18].
Quite recently, finally, the idea of creating a large sample of non-equilibrium trajectories starting from a given initial distribution has been extended to cover whatever distribution that can be sampled starting from an equilibrium or a non-equilibrium, but stationary, dynamics.In particular stationary non-equilibrium ensembles can be generated by suitably restraining standard MD simulations.
In particular, we will illustrate the approach by reporting the results of a study of the time evolution of classical fields, including the onset of convective cells and the relaxation of hydrodynamic interfaces in simple liquids.In this context, we will also briefly address a conceptual difficulty of the approach, due to the possible existence of more than one macroscopic state associated with specific perturbations.In particular cases the problem can be circumvented.
The structure of the paper is as follows.In Section 2 we derive the general framework and specify the possible forms for the initial ensemble.In Section 3 we present a few successful applications of the method.Finally, in Section 4 we try to assess the situation and sketch an outlook.
2. Dynamical Approach to Non-Equilibrium: Theoretical Background

General Formulation
We start considering, in a very general way, a (classical) dynamical system with n degrees of freedom, whose time evolution is described by a set of first order differential equations in a phase space of dimension 2n.We will refer to the phase space variables in a collective way with the vector formalism Γ = {q 1 , p 1 , q 2 , p 2 , . . ., q n , p n }, where the q's and the p's reduce to the usual coordinate-momentum pairs for Hamiltonian dynamics.The equations of motion can be written in the compact form The above equations could be the usual Hamiltonian equations of motion for an isolated system of N particles [19], contain a number of holonomic constraints [20] or represent the more general case of an "extended" system, possibly non-Hamiltonian [21,22], including couplings of the system to a thermal and/or pressure bath by means of a few extra degrees of freedom, so that, in general, n > 3N (see, also, [23]).We will only assume that the dynamics described by Equation (1) are ergodic, i.e., if we wait long enough, all regions of the phase space available to the system, in accord with the imposed conditions, will be explored by the dynamic evolution.With this in mind, the statistical mechanics description of the system requires the introduction of the invariant measure dµ( Γ, d 2n Γ) in phase space [23].We start by introducing the generator of time translations in terms of the Liouville operator, so that the equations of motion can be rephrased in the operator form and formally solved.As the Liouville operator depends explicitly on time, integrating Equation (2) from some initial time t 0 to time t, one obtains an implicit integral equation that can be solved by iteration for each j = 1, 2, . . ., 2n The results can be expressed in closed "operatorial" form with the introduction of the evolution operator where T is the time-ordering operator.Time evolution in phase space can be alternatively expressed in term of the Jacobian J Γ(t), Γ(t 0 ) of the time transformation from Γ(t 0 ) to Γ(t).The phase space element d 2n Γ(t 0 ) at time t 0 transforms into the volume element d 2n Γ(t) = J Γ(t), Γ(t 0 ) d 2n Γ(t 0 ) at time t, where J = det J obeys the differential equation [23] dJ and the phase space compressibility κ is defined by For a Hamiltonian system the compressibility κ vanishes, J Γ(t), Γ(t 0 ) = 1 and the dynamics preserves volume in phase space (Liouville Theorem).More generally, when κ does not vanish, d 2n Γ is no longer a dynamical invariant and one needs to introduce a metric factor to define the invariant measure of the phase space under the dynamical evolution.Starting from the general expression for the Jacobian determinant, one gets where w is the indefinite time integral of κ and Z( Γ(t); t) = exp −w( Γ(t); t) .The dynamically invariant volume element in phase space can be defined as Consider, now, an ensemble of systems whose dynamical evolution is defined by Equation (1).The statistical mechanics is described by the time-dependent probability distribution function in phase space f ( Γ; t) which must obey the global conservation law for probabilities The corresponding local, differential, conservation law can be derived by transforming the integral back to the phase space element d 2n Γ, by using and introducing the phase space density ρ( Γ; t) = Z( Γ; t)f ( Γ; t).The continuity equation to be satisfied is which when expressed in terms of the Liouville operator L and the phase space compressibility κ becomes the "generalized" Liouville equation and reduces to the more "familiar" equation for the probability density f ( Γ; t) However, we must point out that this last equation may lead to confusion if one does not keep in mind that, while the Liouville operator L defines the dynamical evolution of the time-dependent probability density in phase space f , the not-vanishing compressibility κ, hidden in the phase space invariant volume, defines the time evolution of the phase space volume d 2n Γ.
The solution of Equation ( 12) can be retrieved along the same lines followed for Equation (2) and the results can be formally written in closed "operatorial" form where we have introduced the adjoint Ŝ † (t, t 0 ) of the previously defined time evolution operator Ŝ(t, t 0 ) acting on the phase space variables Γ and the phase density ρ 0 = ρ( Γ; t 0 ) at the initial time t 0 .The average over the (non-)equilibrium ensemble of a physical observable O(t) = Ô( Γ) t or, more generally, of a macroscopic field O( x, t) = Ô( x, Γ) t = j Ô( Γ)δ( x − R j ) t (the sum is over the particles) can be defined as We can make the time evolution explicit by means of the adjoint time evolution operator ρ( Γ; t) = Ŝ † (t, t 0 )ρ( Γ; t 0 ) and then, by taking advantage of the fact that Ŝ † is the adjoint of the dynamics, we can transfer the effect of time evolution to the physical observables where Ô( Γ; t) = Ŝ(t, t 0 ) Ô( Γ), i.e., the time evolution along the dynamical trajectory of the system starting from the initial condition Γ(t 0 ) at time t 0 .We have introduced the shorthand notation, • • • ρ 0 , for the averages over the ensemble described by the space density ρ 0 at the initial time t 0 .
Despite the apparent complexity of the time evolution operator Ŝ(t, t 0 ) in Equation ( 5), its action is a task that can be simply accomplished by MD, i.e., by the numerical integration of the evolution defined by Equation (1).Note that all this is possible thanks to the fact that the Liouville equation can be integrated by the method of characteristics.
In the following, we will deal with fluid systems where the relevant macroscopic fields are [24] the density field ( x, t), the velocity field v( x, t) and the temperature field T( x, t): where N is the number of particles and the factor f , usually equal to 3N , counts the number of degrees of freedom in the presence of constraints.

Ensembles at t 0
Equations ( 17) and ( 18) express what we like to call the Onsager-Kubo relations and state that we can obtain the time evolution of a macroscopic observable or of a macroscopic field as the average of the time evolved corresponding microscopic expression over the initial-time-ensemble described by the phase space density ρ 0 = ρ( Γ; t 0 ).
If the ensemble at the initial time t 0 can be simulated by a dynamical system in stationary conditions, then such a probability density function can be sampled by MD, generating a set of (possibly independent) phase space points distributed according to ρ 0 .From each of these points, one can then start an independent dynamical trajectory along which the observables Ô( Γ; t) and Ô( x, Γ; t) can be computed.Finally, by averaging over all the trajectories, the values of the involved observables at time t, one can obtain the macroscopic time-dependent behavior of the system as visualized in Figure 1.
In order to use MD to sample the appropriate initial ensemble at time t 0 , one needs to define, for any specific problem, the dynamical evolution, Equation (1), and the auxiliary conditions to which the systems is subjected.Sometimes, but not always, this will be possible within the Hamiltonian formulation of the dynamics.
Figure 1.Phase space representation of the ensemble of dynamical side-trajectories providing the non-equilibrium statistical averages: in blue, the Molecular Dynamics (MD) trajectory sampling the ensemble at time t 0 ; in black, the individual non-equilibrium trajectories sampling the Non-Equilibrium Molecular Dynamics (D-NEMD) ensemble, over which one can average the time behavior of the observable Ô, as a function of the time t.

D-NEMD Selected Applications
We will now list a number of cases, which will later be illustrated with the corresponding application.Transport properties, like viscosity, thermal conductivity, etc., have been computed and their linearity range investigated by non-equilibrium MD since the 1970s [10][11][12][25][26][27][28][29][30][31][32][33][34].These results were obtained by measuring on a computer the mechanical response when switching on the external (at the beginning Hamiltonian and later on, more generally, also non-Hamiltonian) perturbation applied to a model system initially at equilibrium.In other words, we identify in the present case the ensemble at time t 0 with the statistical mechanics equilibrium ensemble, while the dynamical trajectories are carried out under the influence of an external (time-dependent) force field.
More generally, we can generate (and sample) initial ensembles by less trivial procedures, e.g., in the case of the formation of convective cells, gravity is considered as the external perturbation to be applied on a system initially in a steady state under the effect of a thermal gradient.The ensemble at time t 0 no longer corresponds to the equilibrium one, but it is set up by introducing a stationary boundary perturbation which, in the specific case, is just an ad hoc boundary condition, which models a thermal wall stochastically.Moreover, a confining wall, present in the form of an external field acting at the boundary on each particle, confines the system in the simulation box.This boundary condition is perfectly compatible with the presence of a gravity field.
Another possible case we will consider is the relaxation to equilibrium of an interface between two immiscible liquids, starting from an imposed, non-equilibrium, condition in which the curvature of the interface is maintained by a macroscopic restraint fixing the shape of the initial interface.The ensemble at time t 0 is described by a conditional probability density in which an ad hoc restraint is imposed on a field-like observable.The sample is generated by using an advanced MD sampling technique, where the dynamical trajectory evolves under the effect of a suitable restraining potential, from which we can extract an unbiased sample of the conditional probability density function.Time-dependent averages are then taken over dynamical trajectories generated according to the un-restrained dynamics of the systems.The different situations described are summarized in Figure 2.

Transport and Linear Response
Linear Response Theory is a nice result of the nineteen-fifties in the theory of irreversible processes [6], where well-defined microscopic expressions for all transport coefficients have been derived in terms of a properly chosen perturbation [7,8,35,36].In the Dynamical approach to Non-Equilibrium Molecular Dynamics (D-NEMD) framework it has been possible to investigate the linear and, more generally, the non-linear response by making reference to the canonical ensemble for sampling the initial conditions at time t 0 .

Hamiltonian Perturbations
For a system of particles in three dimensions described by the usual set of Cartesian coordinates and momenta, { R j , P j , j = 1, 2, . . .}, the perturbation can be put in Hamiltonian form by choosing a physical property Â( x| Γ) = j A j ( Γ)δ( x − R j ) that describes the coupling of the system to the applied external local field ψ( x, t) = ϕ( x)χ(t), whose time-dependent intensity χ(t) can be constant or periodic or even arbitrary, generating corresponding flux conditions.Especially important are the cases in which the perturbation is either a step function θ(t − t 0 ) (θ(t > t 0 ) = 1 , θ(t < t 0 ) = 0) or a Dirac delta impulse δ(t − t 0 ), at t = t 0 , after which the system is left free to relax.In the linear regime, the general response can be computed as the superposition of these impulsive responses.One then derives the equations of motion using the standard Hamiltonian route, where we start by separating in the Hamiltonian H( Γ, t) = H 0 ( Γ) + H p ( Γ, t) the time-dependent perturbation term where the Hamiltonian H 0 is the equilibrium Hamiltonian to which one can possibly add the coupling to a thermostat or a barostat, something that can be done in a variety of ways that we do not need to specify here.Indicating generically the possible presence of such couplings to different baths with ellipses, the equations of motion for particle j can be written The structure of the equations of motion can be broken into the two terms of the Liouville operator defined in Equation ( 2), ı L( Γ; t) = ı L0 ( Γ) + ı Lp ( Γ; t), with the partial Liouville operator ı L0 defining the dynamical evolution in phase space for the sampling of the ensemble at time t 0 .Accordingly, the corresponding evolution operator for the stationary dynamics will be called Ŝ0 (t).The dynamics of the time-dependent trajectories will be generated by the t 0 -(time dependent) evolution operator Ŝ(t, t 0 ), obeying the (usual) Dyson equation which (if of interest) can be taken as the basis to develop the perturbative approach, whose first term leads to the Linear Response Theory approach.However, in many cases of interest, for example for constrained systems with a Hamiltonian or non-Hamiltonian structure, it becomes very difficult, if not impossible, to carry out the standard manipulations leading to the correlation function expressions for the linear response [17,18].Nevertheless, a linear (or non-linear) response can always be computationally investigated using the procedure defined by Equations ( 17) and ( 18), as outlined in Figure 1.

Non-Hamiltonian Perturbations
A more general scheme has also been used for bulk perturbations, where the new equations of motion, which cannot be derived from a time-dependent Hamiltonian in a way that remains consistent with applied (periodic) boundary conditions, are obtained from Equation ( 23) by substituting the terms derived from the Hamiltonian perturbation h p , with two sets of "ad hoc" phase space functions A specific, notable, example is the one known under the name of "SLLOD tensor" dynamics [37], where C j χ(t) = −( R j • κ) χ(t) and D j χ(t) = ( P j • κ) χ(t) are coupled with specific, synchronized, Lees-Edwards periodic boundary conditions [38] (see Figure 3), which are needed to establish the tensor κ expressing the desired velocity gradient in the non-equilibrium simulation of viscous flows by molecular dynamics [39][40][41][42].

v v
In the typical setup for a planar Couette flow, one establishes a gradient of the x-component of the velocity along the y-axis of the simulation and measures the response using as an observable the xy component of the pressure tensor σ xy , which can be written, for a system where the potential U is given by a sum of pair interactions, as where In the D-NEMD approach, if the external field term is switched on with a step function perturbation in time at t = t 0 = 0, one can measure the viscous time-dependent response η(t) = − σ xy (t) ρ 0 /γ, where γ is the applied shear rate and the asymptotic value η at long times of η(t) gives the viscosity of the fluid.
For the purpose of illustrating the method in the original applications, when the ensemble at the initial time t 0 is an equilibrium ensemble, we will restrict ourselves to the simple case of shear (Couette) flow.We would like to mention, however, that also elongational flows [41,[43][44][45][46] and, later on, mixed shear-elongational flows [47][48][49] have been simulated both in atomic and molecular fluids.In these cases, it becomes technically much more difficult to maintain for an indefinite length of time the periodic boundary conditions and, for that, we refer the interested reader to [50,51].
In Panel (A) of Figure 4, we show the results of a calculation [52] with a step function perturbation on a Lennard-Jones (LJ) fluid at the triple point, = 0.8442 and k B T p = 0.725 in reduced LJ units, i.e., ε for energy, σ for distances and the particle mass m for masses.The temperature of a system of 2,048 particles was controlled using a Nosé-Hoover thermostat [22], both on the long equilibrium trajectory, which samples the (independent) initial conditions from a canonical ensemble at temperature T p and on the non-equilibrium trajectories to handle the heat produced, especially at high shear rates.The behavior of the time-dependent viscous response for the case of a t 0 = 0 impulsive perturbation with a δ(t − t 0 ) term was used to investigate the range of validity of the Linear Response Theory for very small shear rates by comparison with the running time integral of the corresponding stress autocorrelation at equilibrium [53].[54] and Reference [55] respectively.The solid line is the Lorentzian best-fit to the data and the dashed line is the Ree-Eyring-Eu prediction [56]; Panel (B).The running-time integral (solid line) of the D-NEMD viscous dynamical response to a δ(t − t 0 ) perturbation with γ = 10 −4 , averaged over 4000 trajectories versus the running-time integral (dashed line) of the stress autocorrelation function shows the agreement of D-NEMD results with the Green-Kubo linear reponse theory [52].The error bars, extrapolated using the mean square fluctuations over the 4000 trajectories, increase with time restricting the time range over which the response can be computed.(nb: the same kind of time dependent behavior for η(t) is observed directly when using a step function perturbation).

Non-Equilibrium (Steady State) Initial Conditions at Time t 0
The D-NEMD approach can be used also to follow the transient evolution of a system, which, starting from an out-of-equilibrium state under the effect of a stationary thermodynamic field, reaches a final (different) non-equilibrium state in response to an additional external perturbation.Below, we illustrate the approach with a case worked out in [57].This is the case of the build up of a convective roll in a two-dimensional (2D) model fluid kept in an out-of-equilibrium condition by the presence of a thermal gradient when an external gravity field is (instantaneously) switched on.
The 2D system is composed of N = 5, 401 identical particles in a square box of size L in the xz plane with periodic boundary conditions along the x direction and a pair of confining walls along z obtained by means of an external field ψ(z), acting at the top and the bottom of the simulation box to avoid the drifting away of the particles, which interact with each other via a purely repulsive (2D) Weeks-Chandler-Andersen (WCA) [58] pair potential obtained by truncating the Lennard-Jones potential at its minimum r m = 2 1/6 σ and shifting its value by ε in such a way that both the force and potential are continuous and equal to zero for r r m .
The size of the MD box is L = 84.9 in reduced LJ units, which leads to a density f = 0.75, on average.The confining potential V wall is constructed as the result of a (2D) LJ fluid with continuous constant density w filling the two half planes above and below the periodically replicated MD boxes and has a 10-4 power dependence with parameters defined in [57]: where z is the distance from the box edge, z m is the value at which U (z) has its minimum and V n = 2πn!/ (2 n+1 (n/2)!n), obtaining, in analogy with WCA, a purely repulsive wall.The thermal gradient along the x-direction is obtained by means of two stochastic reservoirs, which are implemented in the two stripe regions at the x-extremities of the MD box (see Panel (A) of Figure 5).The velocities of each particle located in these two stripes are sampled from a 2D Maxwellian distribution f ( v) = e −m v 2 /(2k B T i ) /(2πmk B T i ) at the temperature T i of the stripe, with i = 1, 2 labeling the two reservoirs.Periodic boundary conditions along the x-direction mean that a particle can actually travel from the first (cold) to the second (hot) reservoir.To avoid that a non-thermalized particle in the inner region interacts with both reservoirs, a stripe thickness ∆x T = 1.68 > r m was chosen.
For the system in stationary conditions, each reservoir contains an average of around 100 particles.The reservoir temperatures were chosen to be T 1 = 1.5 and T 2 = 9.9, corresponding to a thermal gradient ∇T = 0.1.Using these conditions, in the absence of gravity, a long stationary trajectory was generated and, from it, a set of 1,000 initial conditions at time t 0 was sampled.Time-dependent trajectories have been generated and then suitable properties averaged at times 0 s t, switching on a gravity field with acceleration g = 0.1 in LJ reduced units (while huge compared with Earth gravity, this is a very small value when compared to the accelerations coming from the interatomic interactions).The behavior of the system was analyzed by coarse graining the MD box into a 15 × 15 mesh of square cells of sides = 5.66 that was used to compute local macroscopic fields.Coarse graining is applied by approximating with the value 1/ 2 for particles inside the cell labeled by (j, k) and centered on the mesh point (x j , z k ) and zero otherwise.The velocity field is calculated as an average over the D-NEMD trajectories where m , the mass density, is given in terms of the average of the number of particles n jk (t) inside the cell (j, k) at time t ( A) The collective behavior of the velocity field can be monitored by calculating its circulation C(t) = P v(x(s), z(s); t)ds on a closed path P , as a function of time.Its evolution is shown in Panel (B) of Figure 5 as a function of the time t after the ignition of the gravity field along a path located in the bulk region, but far enough from the center of the box.The circulation starts from zero at t 0 = 0; its value grows with time and, after a small overshooting, reaches its plateau, stationary, value at t ≈ 200.This is the time at which also both temperature and density fields become stationary.The transient is characterized by correlated oscillations of temperature and density with a very similar period τ = 18 in all cells, but with opposite phases for the cells at the bottom of the box with respect to the ones at the top.The velocity field, shown in Figure 6, is initially null, in Panel (A), acquires first at t ≈ 4.5, in Panel (B), an almost uniform downward component in the direction of the force as a consequence of the ignition of gravity, then it is almost null again at t ≈ 9.0; it shows a maximum reaction to compression at t ≈ 13.5, in Panel (C), and almost vanishes again at t ≈ 18.The cycle restarts and at t ≈ 22.5, in Panel (D), the field is again predominantly in the downward direction, although one can start to see the building up of a convective flow, which is shown in its stationary condition at t ≈ 205, in Panel (E).
We have seen how D-NEMD can be used to illustrate the build up of a convective roll when a gravity field is instantaneously switched on in a system where a stationary (non-equilibrium) thermal gradient was already present.This is not the only case in which a convective roll can be observed.Indeed, keeping the same geometry for the system, i.e., with the gravity field orthogonal to the thermal gradient (Panel A of Figure 5), one could alternatively start from initial conditions in which the system is at equilibrium in the presence of the gravity field and, then, follow the dynamics when the thermal gradient is instantaneously switched on or even start from a homogeneous fluid at equilibrium and instantaneously switch on both the gravity field and the thermal gradient [57].In all these cases, although following different paths, the system eventually reaches the same (macroscopic) final steady state with the formation of a clockwise rotating convective roll, centered at the center of the box.Complications arise if the system is setup with a different geometry, e.g., in the case of Rayleigh-Bénard convection when the direction of the thermal gradient is parallel to the direction of the gravity field.The system has a higher symmetry and rolls rotating both in the counterclockwise and clockwise directions are possible.This implies that the D-NEMD averages cannot be carried out directly, as in the case we have described so far.In fact, now, in the ensemble of individual trajectories, one samples with equal probability the initial conditions leading to clockwise or to counterclockwise rolls.Performing ensemble averages without paying attention to the direction of rotation would give a wrong result.One needs either to enforce a mechanism that breaks this symmetry or to weight trajectories differently, according to the direction of rotation of the convective roll.The latter was the choice applied in [57], where it was simply impossible to fix a priori, by tweaking the initial conditions, the final rotation direction of each trajectory.Instead, the direction of the convective roll rotation was analyzed in the steady-state part of each trajectory by time averaging the velocity field over the last 10,000 steps.The ensemble averages were consistently computed, afterwards, by taking the specular image of the fields when the roll rotation was opposite to the one (arbitrarily) chosen as the reference.

Sampling from Conditional Distributions at Time t 0
Probably the most interesting application of the D-NEMD procedure is when the non-equilibrium dynamical trajectories start from states corresponding to a very unlikely fluctuation and we want to follow dynamically the way in which the system relaxes back to equilibrium.An efficient sampling of the points in phase space representing the initial condition cannot be achieved just by waiting long enough for the desired event to occur during a standard MD trajectory, but more advanced methods are required to enhance the sampling.Whenever the conditions can be described using an "order parameter", i.e., an appropriate phase space function or field, one can define a macroscopic constraint that applied to the system will allow one to explore the interesting, but unlikely, region of phase space.A viable method in many cases is the well-known Blue Moon approach [15,59], where the conditional probability density is constructed by augmenting the dynamical system with a set of holonomic constraints that force the system to explore states on the specific hypersurface of interest in phase space.The points sampled on a dynamical trajectory subject to constraints, however, cannot be directly used to start the time-dependent non-equilibrium trajectories, because of the unphysical additional conditions enforced on velocities to keep the dynamics on the constrained hypersurface.In the Blue Moon approach, such caveats are overcome, the correct procedure is outlined in [15], requiring, first, an apt resampling of the velocities and, then, an appropriate reweighting when computing the time-dependent averages.The Blue Moon approach was initially devised and successfully applied to the calculation of rate constants for activated processes.In particular, the rate constants are defined in terms of the product of two terms: a transition state theory term and a "transmission coefficient" (i.e., the plateau value, in the intermediate time scale, reached following the transient behavior of the time-dependent values of the reaction coordinate(s) [14,60]).The calculation of the transmission coefficient is a task that can be accomplished exactly along the same lines of the proposed D-NEMD approach.
As a more advanced illustration of D-NEMD when sampling from a conditional probability density, we describe the case of the hydrodynamic relaxation to equilibrium of the interface between two immiscible liquids [61].The relaxation can be described by following the time evolution of the difference of the density fields of the two species A and B, ∆ ( x; t) = A ( x; t) − B ( x; t) , and the associated velocity field v( x; t).The distribution at time t 0 corresponds to the stationary conditions of the system subject to a macroscopic restraint that forces a non-equilibrium geometry for the interface.This requires the implementation of a method, like the Blue Moon one, that allows one to sample the conditional probability density associated with the constraint.However, using the Blue Moon approach for vector or field-like constraints can become considerably cumbersome and rather inconvenient in practice, especially for molecular systems where constraints are already used in the force field to impose molecular geometries.A much more practical alternative is to use restrained MD, where one substitutes the constraint with an equivalent restraining potential in terms of an additional coupling parameter, asymptotically reproducing unbiased constrained conditions.Let us summarize the restraint MD approach for the case in which the constraint is imposed on a field-like observable, as for the density difference ∆ ( x, t).Constraining the shape of the interface S corresponds to specifying the set of spatial points { x S } where ∆ ( x S , t) = 0.
According to Irving and Kirkwood [24], we associate with this macroscopic field a microscopic observable on which we need to impose the condition ∆ ( x, t 0 ) = ∆˜ ( x) = 0 on all points { x} in the domain corresponding to the desired geometrical surface at the time t 0 .However, in the numerical approach, one cannot deal directly with a continuous (vector) variable x, therefore, the volume available to the system needs to be discretized over a mesh.With the choice of subdividing the volume in elementary cubic cells, one introduces the space discretization { x α , α = 1, 2, . ..},where the reference point x α coincides with the center of the α − th cell and the microscopic observable field at this point x α is defined as the average F of ∆ˆ over the volume Ω α of the α − th cell: on which we now need to impose the condition F ( x α , Γ) = F ( x α ) = 0 at each of the m points { x α , α = 1, 2, . . ., m}, which correspond to the subset of cells that make up the discretized representation of the chosen interface between the two immiscible liquids.
Consider, now, a system described by the Hamiltonian where H( Γ) is the Hamiltonian of the unconstrained system and k is a tunable parameter that defines the strength of the (harmonic) restraining potential, i.e., the last term on the right-hand side of Equation ( 32).This Hamiltonian can be used to drive either an MC simulation or an MD simulation at a fixed temperature T generating trajectories, stochastic or dynamic, which sample the phase space of the system according to the canonical probability density ρ where are the canonical partition functions and We see, then, more explicitly, that, thanks to the restraint potential, at a given value k of the tunable coupling parameter, we are sampling the conditional probability density of Γ given F , whose limit for k → ∞ is just the ensemble associated with that given fluctuation.The idea of using a biasing potential to sample unlikely points in configuration space was pioneered by Torrie and Valleau for MC simulation [62], then presented in this form in [63].However, while in their case of umbrella/window sampling, the bias is tuned in such a way to sample, in a statistically significant manner, a wider portion of the configuration space, in the restraint MD approach, one considers high enough values of the tunable parameter k with the aim of sampling the conditional probability associated with the portion of phase space representing a rare region of interest.Indeed, using that lim , one recovers, in the limit βk → ∞, the joint probability density where, in the normalizing factor, the probability density of the "condition" { F } is given by The choice of a restraining potential, which depends only on the coordinates of the particles, as in this case, does not influence the probability density in the momentum space, which remains the Maxwellian (equilibrium) distribution and, at variance with the Blue Moon approach, independent points along the stationary restrained MD trajectory can be directly taken as initial configurations representative of the probability density at time t 0 .Moreover, if needed, the restrained MD approach can be further generalized to enforce a more general macroscopic constraint affecting also the momenta of the particles, for example, coupling it with the ad hoc boundary conditions and the localized velocity sampling described in Section 3.2 to impose a non-uniform macroscopic velocity field or a temperature gradient in the system.
The definition of the microscopic field in Equation ( 31) has still one important drawback, which can prompt major issues in particular conditions.In fact, because of the presence of δ-functions in the definition, as a particle crosses the border between one cell and the neighboring one, the integrals that define the macroscopic field at the two corresponding points in space change by, plus or minus, respectively, a finite value introducing discontinuities in the restraining potential in Equation (32).This results in the (highly undesirable) appearance of impulsive terms in the forces on the atoms.Consider the δ( x − R j ) function contribution to the integral in Equation (31) for a specific particle j.This is given by the products of three terms, corresponding to the orthogonal directions in space, where each of them is the difference of the values of the cumulative distribution at the edges of the α − th cell.The cumulative distribution for the delta function δ(ξ) is the step function θ(ξ), {θ(ξ < 0) = 0 , θ(ξ > 0) = 1}, so that each of the above three terms can only be either zero or one.In order to smooth the restraining potential, we need to replace the step function with a smoother function, like a sigmoid, resulting in a continuous variation, between zero and one, of each term in the product.This is equivalent to giving a finite extension to the particle size, resulting in the possibility that a particle contributes, fractionally, to the density field of more than one cell at the same time.In this way, the restraining potential changes smoothly with the motion of the particles in time, without discontinuities when particles cross the cell borders.One possible choice for such function is given by the error function, which corresponds to replacing the delta function by the equivalent Gaussian, e − ξ 2 2a / √ 2πa −→ δ(ξ), in the limit a −→ 0, where the parameter a gives the order of magnitude for the (1D) size of the particle.Within such an approximation, Equation (31) becomes where the function Θ(a, x α , R j ) is the product of three terms corresponding to the integrals along the three spatial components (x 1 ≡ x, x 2 ≡ y and x 3 ≡ z) each involving the evaluation of two values of the error function relative to the border of the cubic cell of length : The two fluids, A and B, are modeled using identical Lennard-Jones particles with mass m and (unique) parameters σ = σ AA = σ BB = σ AB and ε = ε AA = ε BB = ε AB for the LJ potential.The immiscibility is obtained by removing the attractive term for the pair interactions between a particle of type A and a particle of type B, keeping only the purely repulsive part, i.e., taking 12 .The simulation was performed at a fixed temperature k B T = 1.5ε on a system totaling 171,500 particles, of which 88,889 for Fluid A and 82,611 for Fluid B, in a rectangular parallelepiped box with the same width and height w = h ≈ 44 and double length d ≈ 88, corresponding to an average particle density n = 1.024,where all figures are in reduced LJ units.The density and temperature are in the fluid region of the phase space of a pure LJ fluid.In order to follow the behavior of the density and velocity fields, the space was discretized using 5,488 cubic cells arranged on a 14 × 14 × 28 grid.In this way, local field values are obtained averaging out on roughly 30 particles in each cell.For the ensemble at time t 0 , the initial configuration for the interface between Fluid A and Fluid B is defined by selecting the m cells (centered on the mesh points, x α , α = 1, 2, . . ., m), which are cut across by the ideal cylindrical surface, S, where A = 50 is the amplitude that determines the curvature of the surface, which is approximatively placed halfway along the z-direction in the simulation box.The restraint potential is completely defined by the choice of the coupling parameter k = 0.004 in LJ units and the imposed values F ( x α ) = 0, α = 1, 2, . . ., m.The initial configuration was first prepared with the equilibration of a pure Type A fluid at the target temperature and density and then identifying, within the simulation box, as particles of Type A all the particles that are on the red side of the surface S, as particles of Type B all the particles that are on the blue side (see the left panel in Figure 7), taking care of having exactly half of the particles of Type A and half of Type B in the m cells that make up the discretized interface at time t 0 .Periodic boundary conditions are applied in all directions, so that a second flat interface (the condition that minimizes the surface tension) is created at the same time at the sides of the box along the z-direction.Then, the system was equilibrated running restrained MD with a time step δt = 4.56 • 10 −4 in LJ units.Such a rather small value for δt was used to ensure a proper numerical integration of the "stiff" restraining forces.A typical snapshot of the isosurface ∆ ( x) = 0 is shown in the left panel of Figure 7.
A long, 10 6 MD restrained trajectory was then carried out, with that same time step, taking out, at regular intervals of 25,000 steps, the configuration of the system in phase space.A set of 40 independent initial conditions was collected in this way, and from each of them was started, now with a regular time step δt = 4.56 • 10 −3 , a 25,000-steps unrestrained MD trajectory at constant energy, i.e., using the equations of motion derived only from the Hamiltonian H( Γ), given in Equation (32).The D-NEMD averaging procedure was used with those 40 trajectories to compute the time-dependent behavior of the macroscopic density and velocity fields, discretized on the previously mentioned cubic mesh.
In the right panel of Figure 7, the dynamical behavior of the interface is shown by plotting the isosurfaces ∆ ( x) = 0 for four successive times.
One can see that the interface curvature diminishes progressively towards the flat, equilibrium condition, while maintaining (approximatively) both the initial uniformity along the direction of the y-axis and the initial mirror symmetry with respect to the middle yz-plane.Small deviations are present, as expected, considering that this macroscopic field results from averaging over a relatively small sample of 40 independent trajectories.They are compatible with the expected amplitudes of the equilibrium fluctuations of the interface.Relaxation of the initially curved surface reaches the flat equilibrium condition fully in a time lapse of approximatively 20,000 steps, i.e., something just short of 100 LJ time units.One can use this information to estimate the order of the relaxation time and, from it, a maximum value v max ≈ 0.5 in LJ units for the average velocity field at the mid-point along x of the interface, i.e., in the region corresponding to the maximum displacement at time t 0 of the interface.If one takes the LJ parameters of argon (for which σ = 3.405 • 10 −10 m and the unit of time corresponds to τ = 2.156 • 10 −12 s), this value translates to an experimentally convincing velocity of ≈ 80 m • s −1 .The second interface on the sides of the MD box remains flat, with even smaller deviations, all along the 25,000 time steps.There seem to be no significant effects on it as a result of the relaxation process, which takes place in the middle of the box.The reason for this will become evident after looking at the time-dependent behavior of the velocity field.We have shown, in fact, how the D-NEMD approach provides very detailed information on the hydrodynamic behavior of the system and unravels the underlying physical mechanisms.
Starting from Equation ( 20), the discretized velocity field is calculated accordingly to: where we made explicit use of the fact that the particles have identical masses m, regardless of whether they are of Type A or B, and the D-NEMD average is taken over the 40 trajectories with initial conditions sampled from the ρ 0 probability density along the restrained MD trajectory.The results of the calculation are shown in the left panels of Figure 8.The direction of the projections of the velocity field in each cell on the xz-plane (after averaging along the translational symmetry axis y) is represented by a small arrow whose length is proportional to its modulus.The results emphasize how the latter approach returns a much flatter picture of the velocity field, exposing features of the relaxation mechanism that are in marked contrast with the underlying symmetries of the process.(reproduced from [61] with permission from the Physical Chemistry Chemical Physics Owner Societies.) Snapshots at three successive times are shown.In the snapshot taken 500 time steps after t 0 , one can notice the build up of some coherence in the velocity field, which becomes structured in a region that is, along the z-direction, about twice the size of the width of the curved interface, the projection of which is represented by the contour line at the value ∆ = 0.For the sake of clarity, the symmetry yz plane is as well highlighted by a dashed straight line.One can distinguish a quasi-symmetric two-tail profile, where the push towards the edges of the interface appears to be more pronounced than the one, in the opposite direction, in the region near the center of the interface (Panel A).This behavior marks the initial build up of a more stable two-roll velocity profile, of roughly the same width across the interface region in the A and B fluids, which becomes very evident after 3,750 steps (Panel B) and is still neatly visible and qualitatively unchanged, even after 7,500 steps from t 0 , when the curvature of the interface is significantly reduced (Panel C).
One can notice that the size and the intensity the field has decreased, but the profile remains highly symmetrical along the mirror symmetry plane.All along, one can notice also that the field in the extreme sides of the simulation box remains essentially unperturbed, which explains why the second flat interface at the boundary is not affected in any significant way and remains stable during the whole relaxation process of the curved interface.It is very interesting to compare these insights on the hydrodynamic processes underlying the interface relaxation, as given by the time-dependent behavior of the dynamical response calculated using the D-NEMD procedure, with the standard approach used on single trajectory simulations of hydrodynamic processes.Starting from the local equilibrium hypothesis of hydrodynamic theory, based on the assumption of a time scale separation between the fast microscopic motion of the particles and the slower hydrodynamic processes, the macroscopic fields are computed as local time averages on the short time scale τ of the atomistic processes [64].
By applying this approach to the time evolution from a single initial condition, one obtains the results shown in the right panels of Figure 8.At first glance, the velocity field appears to be much smoother than the D-NEMD result, which is relatively noisy due to the limited size (40) of the sample of initial conditions used.However, the picture that is returned is quite different, with the physical mechanism exposed by the D-NEMD procedure effectively washed out by the local time averaging, which also presents a velocity field that does seem to violate the mirror plane symmetry initially imposed on the system at time t 0 , contrary to the more convincing evidence, given from the D-NEMD results, of a relaxation mechanism satisfying, on average, such symmetry at all times along the dynamical trajectory.

Conclusions and Perspectives
In this paper, we have presented a dynamical approach to non-equilibrium MD, which makes it possible to compute, numerically, but, otherwise, rigorously, time-dependent non-equilibrium responses, i.e., to observe directly transient responses in non-stationary regimes.We have shown that using a proper simulation setup, it is possible to go beyond the usual situation of initial equilibrium conditions to treat interesting cases in which the initial condition is either a stationary non-equilibrium or a constrained equilibrium condition dictated by means of a macroscopic constraint, which can be expressed, in a general way, either as a scalar or field-like observable, outlining also the connections of the D-NEMD approach to the Blue Moon method [15] to compute the transmission coefficient contribution to the rate constants of activated processes.
We illustrated a few applications of the method starting from the early, historical, approach to the calculation of transport properties in the lines of Linear Response Theory and beyond, to a couple of recent atomistic simulations of hydrodynamic processes: the establishing of a convective cell, when gravity is switched on in the presence of a stationary thermal gradient, and the relaxation of an initially curved interface between two immiscible liquids.We have shown that the method generates rigorous time-dependent non-equilibrium averages, providing valuable insights on the mechanisms of hydrodynamic processes that can be missed using a method like the local time average, which cannot have a rigorous justification, presents a statistical error that cannot be reduced at will and, finally, as we have seen above, can bias the statistical response.A word of caution is needed, though.The time-dependent ensemble averages are meaningful only if the thermodynamical response is unique [65].

Figure 2 .
Figure2.We distinguish three different classes for the sampling of the initial distribution: equilibrium, direct stationary non-equilibrium simulations and advanced conditional sampling.They are shown to be associated with the corresponding sampling techniques and test-case applications.

Figure 3 .
Figure 3.The Lees-Edwards periodic boundary conditions (Panel A) used to establish a stationary Couette flow (Panel B).In the case of a step function perturbation, periodic images above and below the reference MD cell are translated by an amount ±vδt at each time step, starting from time t 0 .Periodic boundary conditions can be effectively imposed using the equivalent non-orthogonal reference cell, highlighted in red (the actual inclination increases uniformly with time).

Figure 4 .
Figure 4. Panel (A).Comparison of shear viscosity values as a function of the shear rate for the planar Couette flow: (a) D-NEMD asymptotic values from Reference [52]; (b) and (c) average values from stationary non-equilibrium calculations from Reference[54] and Reference[55] respectively.The solid line is the Lorentzian best-fit to the data and the dashed line is the Ree-Eyring-Eu prediction[56]; Panel (B).The running-time integral (solid line) of the D-NEMD viscous dynamical response to a δ(t − t 0 ) perturbation with γ = 10 −4 , averaged over 4000 trajectories versus the running-time integral (dashed line) of the stress autocorrelation function shows the agreement of D-NEMD results with the Green-Kubo linear reponse theory[52].The error bars, extrapolated using the mean square fluctuations over the 4000 trajectories, increase with time restricting the time range over which the response can be computed.(nb: the same kind of time dependent behavior for η(t) is observed directly when using a step function perturbation).

Figure 5 .
Figure 5.The simulation setup (Panel A) for sampling the initial distribution showing the regions where the confining field ψ(z) and the two temperature reservoirs act on the particles.In Panel (B) the average evolution of the circulation of the velocity field is shown after averaging over 200 independent initial conditions (in the inset, we show the path along which the circulation was calculated).

Figure 6 .
Figure 6.The build up of the convective flow is shown by visualizing the local velocity field averaged over 1,000 independent initial configurations as a function of time: (A) t = 0; (B) t ≈ 4.5; (C) t ≈ 13.5; (D) t ≈ 22.5; (E) t ≈ 205.

Figure 7 .
Figure 7. (Left panel) A sampled initial condition for the S isosurface ∆ ( x) = 0 separating the two liquids, A and B (the second planar interface at the long edge of the simulation box is not shown).(Right panel) The evolution in time of the initially curved interface (purple) towards the relaxed planar condition (green).The snapshots are the D-NEMD results averaged over a sample of 40 initial conditions.

Figure 8 .
Figure 8.The behavior of the velocity field in the two liquid regions as a function of time: comparison of the results obtained using the D-NEMD approach averaging over 40 initial conditions (Panels A-C) and the local time averaging procedure (Panels D-F).The results emphasize how the latter approach returns a much flatter picture of the velocity field, exposing features of the relaxation mechanism that are in marked contrast with the underlying symmetries of the process.(reproduced from[61] with permission from the Physical Chemistry Chemical Physics Owner Societies.)