Continuum balances from extended Hamiltonian dynamics

TL;DR

This paper extends the Irving-Kirkwood-Noll continuum balance laws to systems governed by extended Hamiltonian dynamics, accounting for auxiliary degrees of freedom that influence continuum properties and are crucial for multiscale modeling.

Contribution

It provides a physically consistent method to derive continuum balance equations from extended Hamiltonian systems, incorporating additional fluxes and sources.

Findings

Derivation of modified continuum balance laws with extra terms

Identification of auxiliary degrees of freedom effects on fluxes

Implications for multiscale computational schemes

Abstract

The classical procedure devised by Irving and Kirkwood in 1950 and completed slightly later by Noll produces counterparts of the basic balance laws of standard continuum mechanics starting from an ordinary Hamiltonian description of the dynamics of a system of material points. Post-1980 molecular dynamics simulations of the time evolution of such systems use extended Hamiltonians such as those introduced by Andersen, Nos\'e, and Parrinello and Rahman. The additional terms present in these extensions affect the statistical properties of the system so as to capture certain target phenomenologies that would otherwise be beyond reach. We here propose a physically consistent application of the Irving-Kirkwood-Noll procedure to extended Hamiltonian systems of material points. Our procedure produces balance equations at the continuum level featuring non-standard terms, because the presence of auxiliary degrees of freedom gives rise to additional fluxes and sources that influence the thermodynamic and transport properties of the continuum model. Being aware of the additional contributions may prove crucial when designing multiscale computational schemes in which information is exchanged between the atomistic and continuum levels.

Equations

$$\dot{\zeta }(t)=\{z,H\}{|}_{z=\zeta (t)},$$

$$\{A,B\}:=\frac{\partial A}{\partial q}\cdot \frac{\partial B}{\partial p}-\frac{\partial A}{\partial p}\cdot \frac{\partial B}{\partial q}.$$

$${\mathrm{div}}_z\{z,H\}=0,$$

$$\dot{\tilde{\rho }}(\zeta (t),t)=0$$

$$\frac{\partial }{\partial t}\tilde{\rho }(z,t)-\{H(z,t),\tilde{\rho }(z,t)\}=0,$$

$$\ell [\cdot ]:=-\{H(z,t),\cdot \}$$

$$\{A,B{\}}_{{\mathcal{Z}}_{\mathrm{s}}}:=\frac{\partial A}{\partial q_{\mathrm{s}}}\cdot \frac{\partial B}{\partial p_{\mathrm{s}}}-\frac{\partial A}{\partial p_{\mathrm{s}}}\cdot \frac{\partial B}{\partial q_{\mathrm{s}}};$$

$$\{A,B\}=\{A,B{\}}_{{\mathcal{Z}}_{\mathrm{s}}}+\{A,B{\}}_{{\mathcal{Z}}_{\mathrm{e}}}.$$

$$\frac{\partial \tilde{\rho }}{\partial t}=\{H,\tilde{\rho }{\}}_{{\mathcal{Z}}_{\mathrm{s}}}+\{H,\tilde{\rho }{\}}_{{\mathcal{Z}}_{\mathrm{e}}}.$$

$$⟨o⟩(t):={\int }_{\mathcal{Z}}\tilde{\rho }(z,t)\hat{o}(z)dz.$$

$$\frac{\partial ⟨o⟩}{\partial t}={\int }_{\mathcal{Z}}\{H,\tilde{\rho }\}\hat{o}dz.$$

$$\frac{\partial ⟨o⟩}{\partial t}={\int }_{\mathcal{Z}}\{H,\tilde{\rho }{\}}_{{\mathcal{Z}}_{\mathrm{s}}}\hat{o}dz+{\int }_{\mathcal{Z}}\{H,\tilde{\rho }{\}}_{{\mathcal{Z}}_{\mathrm{e}}}\hat{o}dz.$$

$$q=({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N)\text{and}p=(m_1{\dot{\mathbf{r}}}_1,\dots ,m_N{\dot{\mathbf{r}}}_N),$$

$$\hat{o}(z,\mathbf{r})=\sum _k{\hat{o}}_k(z)\delta ({\mathbf{r}}_k-\mathbf{r}),$$

$$⟨o⟩(\mathbf{r},t)$$

$$={\int }_{\mathcal{Z}}\hat{o}(z,\mathbf{r})\tilde{\rho }(z,t)dz$$

$$\begin{gathered} H_{\mathrm{NH}}(q^{\prime },p^{\prime };Q,T):=\sum _k\frac{|{\mathbf{p}}_k{|}^2}{2m_k{s}^2}+\frac{p_s^2}{2Q}+V(q) \\ +(3N+1)k_{\mathrm{B}}T(\ln s-1), \end{gathered}$$

$$q^{\prime }=(q,s)\text{and}{p}^{\prime }=({\mathbf{p}}_1,\dots ,{\mathbf{p}}_N,p_s).$$

$$\frac{{\mathbf{p}}_k}{s}=m_k{\mathbf{v}}_k.$$

$$t={\int }_0^{\tau }\frac{d\alpha }{s(\alpha )}.$$

$$V^{\mathrm{i}}(q):=\frac{1}{2}\sum _{j\ne k}{V}_{jk}(|{\mathbf{r}}_j-{\mathbf{r}}_k|)$$

$$V^{\mathrm{e}}(q):=\sum _k{V}_k^{\mathrm{e}}({\mathbf{r}}_k),$$

$$\rho(\mathbf{r},t)\equiv\Big{\langle}\sum_{k}m_{k}\delta(\mathbf{r}_{k}-\mathbf{r})\Big{\rangle},$$

$$\rho \mathbf{v}(\mathbf{r},t)$$

$$\displaystyle\mbox{}=\Big{\langle}\sum_{k}m_{k}\mathbf{v}_{k}\delta(\mathbf{r}_{k}-\mathbf{r})\Big{\rangle},$$

$${ϵ}_K(\mathbf{r},t)$$

$$\displaystyle\mbox{}=\Big{\langle}\sum_{k}\frac{m_{k}}{2}|\mathbf{v}_{k}|^{2}\delta(\mathbf{r}_{k}-\mathbf{r})\Big{\rangle},$$

$$\epsilon_{V}(\mathbf{r},t)\equiv\mbox{}\Big{\langle}\sum_{j\neq k}\frac{1}{2}V_{jk}(|\mathbf{r}_{j}-\mathbf{r}_{k})|)\delta(\mathbf{r}_{k}-\mathbf{r})\Big{\rangle}\\ \mbox{}+\Big{\langle}\sum_{k}V^{\mathrm{e}}_{k}(\mathbf{r}_{k})\delta(\mathbf{r}_{k}-\mathbf{r})\Big{\rangle}.$$

$${ϵ}_{p_s}(t)$$

$${ϵ}_s(t)$$

$$\sum _k\frac{1}{N}\delta ({\mathbf{r}}_k-\mathbf{r})$$

$$n(\mathbf{r},t)\equiv\Big{\langle}\sum_{k}\frac{1}{N}\delta(\mathbf{r}_{k}-\mathbf{r})\Big{\rangle}.$$

Full text

I Introduction

Whereas classical continuum mechanics rests on the description of the deterministic evolution of macroscopic quantities, classical statistical mechanics studies probabilistic features of microscopic observables. Continuum mechanics provides phenomenological descriptions typically aimed at determining over a time interval and a matter-comprised space region the values taken by a list of fields, such as mass density, motion, stress, and temperature. This is done by solving initial/boundary-value problems governed by sets of partial-differential evolution equations embodying such fundamental laws as the balances of mass, linear momentum, and total energy, specialized for a constitutive class of choice. Statistical mechanics provides expressions for equations of state and forces based on a chosen potential of internal interactions in combination with the ensemble properties. Their different natures notwithstanding, these edifices can supply consistent and complementary information when applied to one and the same physical system with the objective of building a multiscale method. When the statistical information retrieved at the atomistic level is used to obtain deterministic predictions at the continuum level, or vice versa, there are two main consistency issues to settle: one concerning macroscopic constitutive responses to deformations and microscopic interatomic potentials; the other concerning the evolution equations in the two frameworks. In this paper, we deal with the latter issue.

To close the atomistic-to-continuum scale gap to the extent of associating a set of deterministic balances with a probabilistic account of the Newtonian evolution of a system of material points, we take the bottom-up path pioneered by Irving & Kirkwood Irving and Kirkwood (1950) in a paper that appeared in 1950. Their path was completed a few years later by Noll,Noll (1955); *NT10 who derived closed-form integral expressions relating microscopic observables to the continuum mechanical notions of stress and heat flux.

The Irving–Kirkwood–Noll (IKN) procedure was derived for the classical Hamiltonian dynamics of a collection of material points. This represents a significant limitation, because the statistical information that can be retrieved via Molecular Dynamics (MD) from a classical Hamiltonian system refers only to the microcanonical ensemble (NVE), in which the number of particles, the volume, and the energy are conserved.

An innovative and effective technique for accessing different statistical ensembles was introduced by AndersenAndersen (1980) in 1980. Andersen’s idea was to view the MD computational cell as a carrier of mesoscopic information, conveyed by an extended Hamiltonian, in which the introduction of an additional dynamical variable—the volume of the cell in his case—allows the system to explore the phase space along trajectories for which the value of the extended Hamiltonian, but not the physical energy, is conserved. As a consequence, while the extended system retains statistics pertinent to the extended microcanonical ensemble, the statistical ensemble induced on the physical degrees of freedom is modified to obtain the conservation of two observables—the pressure applied to the cell and the related enthalpy. Prompted by Andersen’s paper, Parrinello & RahmanParrinello and Rahman (1980); *PR2 let not only the volume but also the shape of the computational cells change. They showed that the availability of shape-related degrees of freedom allows for simulating stress-induced phase transitions in solids; to paraphrase the words of Nosé,Nosé (1991) their papers revolutionized the approach to investigating structural phase transitions. A few years later, Nosé,Nosé (1984a); *Nos84 inspired once again by Andersen’s paper, conceived the idea of an extended Hamiltonian befitting MD simulations consistent with a constant-temperature canonical ensemble (NVT), in the sense that “the canonical distribution is realized in a physical system”.Nosé (1991) Simulations of extended Hamiltonian dynamics were also considered within the more general class of non-Hamiltonian MD in important articles on their statistical properties.Tuckerman, Mundy, and Martyna (1999); *TucLiu01 A development that, to our knowledge, went unnoticed in the literature is that Noll,Noll (1955) whose treatment allows for the presence of non-conservative external forces, actually considered the extension to a special class of non-Hamiltonian evolutions in as early as 1955.

In the present paper, we couple the IKN procedure—in itself a microscopic-to-macroscopic scale-bridging procedure—with two extended Hamiltonians. The first models the Nosé–Hoover (NH) thermostat.Nosé (1984a); *Nos84; Nosé (1991); Hoover (1985); *Hoo07 The second is a Hamiltonian of Andersen–Parrinello–Rahman (APR) type, which includes a complete accounting of the kinetic energy. The extended Lagrangian used by Parrinello & Rahman,Parrinello and Rahman (1980); *PR2 which featured a simplified accounting of the kinetic energy, was a generalization of that previously devised by Andersen with the objective of achieving an isobaric ensemble,Andersen (1980) whence our use of the acronym APR.

As is well-known, the IKN procedure tells us what statistical observables to associate with the continuum mechanical notions of mass, linear momentum, and energy density. On evolving these observables à la Liouville with respect to a chosen Hamiltonian, three consequences of the Newtonian motion of a system of material points are derived. These can be associated with three balance laws basic to continuum mechanics, namely mass conservation, balance of linear momentum, and balance of energy. The association involves convenient identifications of other statistical observables with the continuum mechanical notions of stress and heat flux. In recent years, there has been a renewed interest in these matters; with no pretense to exhaustivity, we mention the works of MurdochMurdoch (2010); *Mur12 (who proposes an identification procedure alternative to that used in the IKN procedure along lines anticipated by HardyHardy (1982)), Admal & Tadmor,Admal and Tadmor (2010); *AdmTad11; *AdmTad16 Seguin & Fried,Seguin and Fried (2012); *SegFri13 and others.Wajnryb, Altenberger, and Dahler (1995); Denniston and Robbins (2004); *DenRob06; Goldhirsch (2010); Humenyuk and Tokarchuk (2011); Lehoucq and Sears (2011); Yang, Wu, and Li (2012); Capriz and Mariano (2014); Davydov and Steinmann (2014a); *DavSte14b; *DavSte15

We remark that in the case of microstructured continua a microscopic definition of angular momentum is often necessary and would significantly affect the continuum balance of the density of angular momentum, as recently discussed by Seguin & FriedSeguin and Fried (2012); *SegFri13 and Davydov & Steinmann.Davydov and Steinmann (2014a); *DavSte14b; *DavSte15 Nevertheless, in the present treatment, we restrict attention to continua associated with systems of material points that do not possess rotational degrees of freedom. In such a context, the continuum balance of angular momentum is equivalent to the symmetry of the Cauchy stress tensor, a property guaranteed under fairly general assumptions concerning the nature of the internal interactions.

Our goal is to implement the IKN procedure for both the NH and the APR Hamiltonians. Our main expectation is to provide an indication as to which continuum equations (and of what kind) should be coupled to MD simulations when building consistent multiscale schemes. Our interest in the NH thermostat is motivated by the assumption of isothermal evolution, posited for many continuum mechanical models. On the other hand, the APR approach has the peculiarity of introducing a dynamic generalization of the Cauchy–Born rule,Podio-Guidugli (2010, 2011) which makes it directly relevant to certain recently proposed multiscale computational schemes.Davydov, Pelteret, and Steinmann (2014); Ulz (2015); Li and Tong (2015); *TongLi15 A discussion of how thermostatting techniques can affect the evaluation of transport coefficients in linear response theory (which are related to the balance of linear momentum) is presented by Evans & Morriss.Evans and Morriss (2008) Our goal, however, is to obtain a full set of continuum balances featuring a more accurate account of the macroscopic effects induced by the extended Hamiltonian dynamics.

When setting up the IKN procedure for an extended Hamiltonian system, that is to say, a particle system whose Hamiltonian $H$ is nonstandard, two key observations are that (i) certain microscopic observables must be defined by combining ordinary physical variables with other possibly “unphysical” variables and that (ii) all macroscopic observables, whatever their microscopic antecedents, consist of ensemble averages weighted with respect to probability densities defined over the extended phase space. This in itself guarantees a basic statistical-continuum consistency in the derived balances of macroscopic observables. However, we find that an extended Hamiltonian evolution always entails the presence of new terms in the balance of energy. This signals that the choice of a peculiar probabilistic dynamics at the microscopic scale modifies the thermodynamic and transport properties of the system at the macroscopic scale. Awareness of this denouement is crucial when devising computational schemes in which information is exchanged between the macroscopic and microscopic levels.

Our paper is organized as follows. In Section II, we recall a few well-known notions regarding Liouville’s evolution of probability densities and observables for Hamiltonian particle systems; we also specify the changes in a format necessary to deal with extended Hamiltonian systems. In Section III, we discuss the application to extended Hamiltonian systems of the IKN procedure for obtaining macroscopic balance equations; in particular, we analyze in detail the cases of the NH thermostat and of APR dynamics. The main results of our study, as well as their relevance to multiscale computational schemes, are finally summarized in Section IV.

II Liouville’s evolution of observables

Consider a system defined by a Hamiltonian $H$ depending on Lagrangian coordinates $q$ and conjugate momenta $p$, with $z=(q,p)$ denoting a generic point in the phase space $\mathcal{Z}$. The evolution of the system is given by Hamilton’s equations

[TABLE]

where $\{\cdot,\cdot\}$ denotes the Poisson bracket defined, for any pair $A$ and $B$ of functions of $(z,t)=(q,p,t)$, as

[TABLE]

II.1 Liouville’s theorem and Liouville’s equation

Given a Hamiltonian system, it is readily checked that

[TABLE]

which is a statement of Liouville’s theorem. A geometrical statement of Liouville’s theorem based on (1) and (3) is that any Hamiltonian evolution of a region in the phase space $\mathcal{Z}$ is locally volume-preserving.

The incompressibility condition expressed by (3) has another important consequence: any probability density function $\tilde{\rho}$ defined on $\mathcal{Z}$ at the initial time $t_{0}$ is convected along Hamiltonian trajectories in the sense that

[TABLE]

whenever $\zeta(t)$ is a solution of (1). This can be equivalently expressed by Liouville’s equation

[TABLE]

where

[TABLE]

is the Liouville operator for a system of Hamiltonian $H$. Liouville’s equation evolves a given initial value $\tilde{\rho}(z,t_{0})=\tilde{\rho}_{0}(z)$ of the probability measure at any given point of $\mathcal{Z}$ into its current value $\tilde{\rho}(z,t)$ at that same point. The kernel of the Liouville operator $\ell$ consists of all the stationary probability measures, that is, those which are independent of time $t$. Notice that the discussion of the present section applies also to time-dependent Hamiltonians, even though, in the applications discussed subsequently, we only encounter time-independent Hamiltonians.

II.2 Extended form of Liouville’s equation

If we consider a decomposition $\mathcal{Z}=\mathcal{Z}_{\mathrm{s}}\times\mathcal{Z}_{\mathrm{e}}$ of the phase space $\mathcal{Z}$, where $\mathcal{Z}_{\mathrm{s}}$ is associated with a subset of Hamiltonian coordinates, denoted by $(q_{\mathrm{s}},p_{\mathrm{s}})$, and $\mathcal{Z}_{\mathrm{e}}$ is associated with additional degrees of freedom that extend the set of coordinates, we can then define the restricted Poisson bracket for any pair $A$ and $B$ of scalar functions of $(z,t)$ by

[TABLE]

the complementary bracket $\{A,B\}_{\mathcal{Z}_{\mathrm{e}}}$ has an analogous expression. With such definitions, we obtain an additive decomposition of the Poisson bracket:

[TABLE]

Moreover, we can exploit (8) to express Liouville’s evolution of probability densities for an extended Hamiltonian system as

[TABLE]

II.3 Extended evolution of observables

A microscopic observable is a mapping $\{z\mapsto\hat{o}(z)\}$ defined over the phase space $\mathcal{Z}$. At this stage, the tensorial order of $\hat{o}$ can be left arbitrary. The corresponding macroscopic observable is the expected value of $\hat{o}$ at time $t$ with respect to the chosen probability density $\tilde{\rho}(z,t)$, namely

[TABLE]

An important application of Liouville’s equation is to compute the time derivative of the expected value $\langle o\rangle$ of a microscopic observable $\hat{o}$, which is given by

[TABLE]

For extended Hamiltonian systems, substituting the additive decomposition (8) in (11) gives

[TABLE]

Irving & KirkwoodIrving and Kirkwood (1950) (in a footnote on their p. 822) commented on the possible presence of additional degrees of freedom, leading to (12). However, they considered only observables that are independent of the additional degrees of freedom. At variance with them, we will consider observables that depend also on such coordinates, entailing non-vanishing contributions from both terms in the right-hand side of (12).

III The IKN procedure for extended Hamiltonian systems

We begin by briefly recapping the classical IKN procedure, which is applied to a collection of $N$ material points moving in three-dimensional space. The Hamiltonian coordinates are

[TABLE]

where $\mathbf{r}_{k}$ is the current position vector, with respect to a fixed origin, of the $k\textrm{-th}$ material point, the mass of which is $m_{k}$. The fundamental presumption is that all macroscopic properties of matter are deducible from a microscopic discrete picture; the idea is to exploit the similarity in format of continuum balance equations and statistical Liouville flows.

The procedure consists of three steps: (i) The densities of mass, linear momentum, kinetic energy, and internal energy are associated with the ensemble averages of certain microscopic observables, the common form of which is

[TABLE]

where $\mathbf{r}$ is the position vector of the typical space point. (ii) Each of the above ensemble averages is evolved, in the order given, à la Liouville. (iii) The so-obtained Liouville flows are identified term-by-term with their continuum mechanical counterparts. An all important feature of the procedure is that the representation (14) of each microscopic observable—a finite sum of Dirac distributions, each supported at one of the material points constituting the system—implies that the associated macroscopic observable

[TABLE]

is a time-dependent spatial field.

When applying the IKN procedure to an extended Hamiltonian system, it is important that the basic macroscopic observables maintain their physical meaning. This poses restrictions on the definitions of the relative microscopic observables, where possibly ‘unphysical’ additional Hamiltonian coordinates inevitably enter. Moreover, the additional Hamiltonian coordinates affect the structure of both the phase space and the probability measure entering ensemble averages. Altogether, it can happen that certain terms of the continuum balances produced by the IKN procedure turn out to have forms different from the usual ones. In the following, we exemplify these developments for two paradigmatic extended Hamiltonian systems.

III.1 Application to the NH thermostat

The first system to which we apply the IKN procedure is the extended Hamiltonian system introduced by Nosé,Nosé (1984a); *Nos84; Nosé (1991) which forms the basis for the thermostatting strategy proposed by Hoover.Hoover (1985); *Hoo07 We shall see, a physically sound definition of observables produces macroscopic balances that are consistent with classical results. Nevertheless, additional source terms that stem from the extended Hamiltonian are present both in the mass balance and in the energy balance. At the continuum level, those additional terms can alter the spatial and temporal distributions of the fields. At the microscopic level, it is, however, the ensemble that changes.

III.1.1 The extended Hamiltonian

Nosé’s extended Hamiltonian can be written as

[TABLE]

the underlying coordinates being

[TABLE]

The string $q^{\prime}$ is the standard string of Lagrangian coordinates in $\eqref{IKobs1}_{1}$, augmented by a dimensionless scalar coordinate $s$; $p_{s}$, the momentum conjugate to $s$, augments a string of virtual momenta $\mathbf{p}_{k},\;k\in\{1,\ldots,N\}$. For the physical velocity $\boldsymbol{v}_{k}$ of the $k$-th material point, the relationship between the virtual momentum $\mathbf{p}_{k}$ and the physical momentum $m_{k}\mathbf{v}_{k}$ is mediated by $s$ as follows:

[TABLE]

The role assigned by Nosé to $s$ is to scale the physical time $t$ in terms of the virtual time $\tau$ involved in the Hamiltonian evolution associated with $H_{\mathrm{NH}}$ according to:

[TABLE]

This interpretation of $s$, which depends on the virtual time, is consistent with Nosé’s intention to achieve temperature control by controlling particle velocities, granted the understanding that the temperature of a system is related to the average of its kinetic energy.Nosé (1991) For this formulation to make sense, $s$ must be positive. Consistent with this observation, we stipulate that the probability density function $\tilde{\rho}(q^{\prime},p^{\prime},\tau)$ must vanish identically for $s\leq 0$.

The expression of $H_{\mathrm{NH}}$ in (16) features two parameters: $Q$, which enters the new kinetic term; and $T$, which enters the new potential energy, where $N$ is the number of material points and $k_{\mathrm{B}}$ is Boltzmann’s constant. Specifically, $Q$ modulates the new kinetic term proportional to $p_{s}^{2}$ and plays the role of an effective thermal inertia, associated with the characteristic frequency of thermal fluctuations, and $T$ is the target temperature. The role of the ‘entropic’ potential $(3N+1)k_{\mathrm{B}}T(\ln s-1)$ is to induce the constant-temperature ensemble characterizing the equilibrium properties of the system.

Finally, in (16) the potential energy $V$ has the standard form $V=V^{\mathrm{i}}+V^{\mathrm{e}}$, where the internal contribution

[TABLE]

accounts for pairwise interactions between material points and

[TABLE]

is the unary potential of external forces.

III.1.2 Definition of the observables

We now identify the statistical counterparts of the fields entering the basic continuum balances as macroscopic observables under the form of ensemble averages of suitable physically relevant microscopic observables. We do so by applying the general IKN prescription recalled in the introductory paragraph of this section. By construction, such ensemble averages depend on the virtual time $\tau$; in what follows, the inversion of relation (19) is implied. Specifically, we define the mass density

[TABLE]

the linear momentum density

[TABLE]

the kinetic energy density

[TABLE]

and the potential energy density

[TABLE]

The presence of additional degrees of freedom and energy terms in $H_{\mathrm{NH}}$ calls for the definition of further observables. Since the additional degrees of freedom are not explicitly ascribable to individual material points, it is not immediately obvious whether the corresponding observables should be given the form (14), entailing generally nontrivial variations with spatial position. These observables could also be viewed as collective properties of the system and be described by uniform spatial fields. To investigate the implications of those two perspectives, we define both collective observables, which are independent of the space point $\mathbf{r}$, and corresponding distributed observables (indicated with superposed bars) of IKN form.

The collective observables associated with the extra kinetic energy density and the entropic energy density are

[TABLE]

where we have introduced $A:=(3N+1)k_{\mathrm{B}}T$ and $\omega_{\mathrm{ref}}$ is used to denote the volume of the computational cell.

The construction of an IKN observable for a collective quantity can always be achieved by multiplying that quantity by the microscopic precursor

[TABLE]

of the number density

[TABLE]

Following that prescription, we define the distributed observables

[TABLE]

The foregoing definitions and the corresponding evolution terms represent an important novelty of the continuum theory associated with an extended Hamiltonian system via the IKN procedure. We anticipate that the main difference in choosing collective or distributed observables is that the collective ones, being spatially uniform, are not convected with the motion of material points. For this reason, they produce only uniform source terms in the energy balance, whereas distributed observables produce both convective and diffusive terms in that balance.

III.1.3 The mass and momentum balances

By applying Liouville’s equation (9), which in this case takes the form

[TABLE]

we can deduce the macroscopic balances of the defined observables.

The time derivative of the mass density is

[TABLE]

where the source term is given by

[TABLE]

As for the linear momentum density, we find that

[TABLE]

where the Cauchy stress tensor ${\mathbf{T}}$ is the sum of the kinetic stress tensor

[TABLE]

and the stress tensor associated with internal interactions, namely

[TABLE]

with

[TABLE]

and where the external force field is

[TABLE]

The detailed derivation of formulae (37) and (38) was provided originally by Noll,Noll (1955); *NT10 who more recently generalized and refined his calculations.Noll (2010b) Note that the very presence of the source term $\sigma_{\rho}$ in the mass balance, a term that has no classical counterpart, is dictated by the form (16) of the extended Hamiltonian $H_{\mathrm{NH}}$, which however does not give rise to specific new terms in the balance of linear momentum.

III.1.4 Energy balances

In view of our previous definitions of collective and distributed energy observables, we must balance two energy densities, namely

[TABLE]

and

[TABLE]

In both instances, the time derivatives of $\epsilon_{K}$ and $\epsilon_{V}$ produce the counterpart of the continuum mechanical heat flux $\mathbf{q}$, which involves the sum of a kinetic contribution

[TABLE]

an interaction contribution

[TABLE]

and a transport contribution

[TABLE]

Energy balance with collective observables. The time derivatives of the collective energy densities $\epsilon_{s}$ and $\epsilon_{p_{s}}$ defined in (26) and (27) are

[TABLE]

and

[TABLE]

these derivatives induce two uniform source terms, $\sigma_{\epsilon}^{p_{s}}$ and $\sigma_{\epsilon}^{s}$.

All in all, the balance of the energy density $\epsilon$ defined in (40) takes the form

[TABLE]

where $\mathbf{q}$ and $\sigma_{\epsilon}$ are given by

[TABLE]

In addition to the time-dependent but spatially uniform sources $\sigma_{\epsilon}^{p_{s}}$ and $\sigma_{\epsilon}^{s}$, $\sigma_{\epsilon}$ features the classical term

[TABLE]

and two other new terms, namely

[TABLE]

and

[TABLE]

While the time evolution of $\epsilon_{K}$ and $\epsilon_{V}$ produces spatially dependent diffusive fluxes, convective terms, and sources, the time evolution of $\epsilon_{p_{s}}$ and $\epsilon_{s}$ generates only uniform sources. Energy balance with distributed observables. To deduce the balance of the energy density $\bar{\epsilon}$ defined in (41), we compute the time derivatives of the distributed observables $\bar{\epsilon}_{p_{s}}$ and $\bar{\epsilon}_{s}$ defined in (30) and (31). Specifically, we find that

[TABLE]

with

[TABLE]

and

[TABLE]

Additionally, we find that

[TABLE]

with

[TABLE]

and

[TABLE]

All in all, the balance of the energy density $\bar{\epsilon}$ reads

[TABLE]

where $\mathbf{q}$, and $\sigma_{\bar{\epsilon}}$ are given by

[TABLE]

In this case, each of the energy densities produce flux and convective terms influencing the energy balance (59).

III.2 Application to the APR extended Hamiltonian

In this section, we apply our generalized IKN procedure to an extended Hamiltonian system akin to but more general than that considered by Parrinello & Rahman,Parrinello and Rahman (1980); *PR2 whose successful predictions of stress-induced displacive phase transitions in certain crystalline materials initiated an irreversible change in format in the MD simulation of such phenomena. We reiterate that we chose this application not only for its intrinsic importance, but also stimulated by recently proposed multiscale numerical schemesUlz (2015); Li and Tong (2015); *TongLi15 for coupling atomistic and continuum models. Importantly, those numerical schemes are based on an APR approach, and we view determining the form of the implied continuum balances to be an endeavor of particular importance and interest. In a departure from previous treatments, we do so by employing the general and exact form of the kinetic energy, a form which reduces to that postulated by Parrinello & Rahman under specific conditions analyzed by Podio-Guidugli.Podio-Guidugli (2010) We regard maintaining this level of generality as important for obtaining a consistent and unprejudiced physical description at different scales.

In what follows, we denote by a simple juxtaposition the contraction of a single tensorial index between a tensor and a vector or between two tensors. For both tensors and vectors, we denote by $\boldsymbol{A}\cdot\boldsymbol{B}$, with a single centered dot, the scalar product defined through the trace operator by $\mathrm{tr}(\boldsymbol{A}^{\mathsf{T}}\mskip-3.0mu\boldsymbol{B})$.

III.2.1 The extended Hamiltonian

The additional degrees of freedom that extend the ordinary Lagrangian within an APR framework are those used to describe the deformations of the computational cell; they comprise an invertible second-order tensor $\mathbf{F}$, with positive determinant, which maps the referential position vector $\mathbf{s}_{k}$ of the $k$-th material point into its current position vector $\mathbf{r}_{k}$:

[TABLE]

Accordingly, the Lagrangian coordinates of this system are given by $q=((\mathbf{s}_{k})_{k=1}^{N},\mathbf{F})$. The time-dependence of $\mathbf{F}$ renders (62) a generalization of the Cauchy–Born rule, from which it follows that

[TABLE]

As we shall see, (63) is crucial for a correct definition of physical observables. Importantly, it must be used to compute the Lagrangian form of the kinetic energy:

[TABLE]

The Lagrangian kinetic energy $K_{L}$ is a symmetric and positive semi-definite quadratic form in the variables $((\dot{\mathbf{s}}_{k})_{k=1}^{N},\dot{\mathbf{F}})$, which depends parametrically on $q$. The corresponding Hamiltonian kinetic energy $K$ is the Legendre–Fenchel transform of $K_{\mathrm{L}}$. Although $K_{\mathrm{L}}$ depends on $\mathbf{s}_{k}$ and $\dot{\mathbf{s}}_{k}$, it is independent of $\mathbf{r}_{k}$ and, therefore, so is $K$. We thus have

[TABLE]

The conjugate momentum $\mathbf{p}_{k}$ and the physical momentum $m_{k}\mathbf{v}_{k}$ of the $k$-th material point are linked by the transformation

[TABLE]

from which we deduce the alternative expression $\mathbf{F}^{-\mathsf{T}\mskip-3.0mu}\mathbf{p}_{k}$ for the physical momentum of the $k$-th material point.

The APR extended Hamiltonian that we consider is

[TABLE]

with $p=((\mathbf{p}_{k})_{k=1}^{N},\mathbf{G})$. Here, as usual, the potential energy depends on the pairwise interactions of material points in the positions they currently occupy:

[TABLE]

Given our present purposes, the addition of an external potential would not add anything of importance. What makes an approach of APR type interesting is the last term of (67), $-\omega_{\mathrm{ref}}{\mathbf{P}}\cdot\mathbf{F}$, which is of an enthalpic character. In that term, $\omega_{\mathrm{ref}}>0$ denotes the constant reference volume of the computational cell, while the time-independent control parameter ${\mathbf{P}}$ specifies the prescribed macroscopic Piola stress in the referential configuration of the computational domain throughout the MD simulation.

In this extended Hamiltonian system, the evolution of any probability density $\tilde{\rho}(q,p,t)$ is determined by Liouville’s equation (9), which takes the following explicit form:

[TABLE]

III.2.2 Definition of the observables

The discussion in the previous sections provides us with precise indications about the appropriate form for the various observables. We define the mass density

[TABLE]

the linear momentum density

[TABLE]

the kinetic energy density

[TABLE]

and the potential energy density

[TABLE]

In keeping with the treatment of Nosé’s extended Hamiltonian (16), we can further define the collective observable

[TABLE]

and its distributed counterpart

[TABLE]

III.2.3 The continuum balances

In the case of the APR extended Hamiltonian, the classical form of the mass balance is preserved, as is the form of the linear momentum balance, in which the stress tensor ${\mathbf{T}}={\mathbf{T}}_{K}+{\mathbf{T}}_{V}$ is obtained from (36)–(38) modulo the substitutions

[TABLE]

Choosing the collective form (74) for the enthalpic energy density generates the energy balance

[TABLE]

with

[TABLE]

and

[TABLE]

In these relations, the terms with suffixes $K$, $V$, and $T$ have the form of the corresponding terms in the application of the IKN procedure to the NH extended Hamiltonian. If the computational cell is viewed as a homogeneously deforming elastic body, then $\omega_{\mathrm{ref}}\mathbf{F}$ constitutes an extensive variable and $\sigma_{\epsilon}$ provides a reckoning of the power expenditure associated with changing the shape of the computational cell under the influence of the Piola stress ${\mathbf{P}}$.McLellan (1980); Man (2017)

If we choose instead the distributed form (75) for the enthalpic energy density, the energy balance reads

[TABLE]

where

[TABLE]

and

[TABLE]

The stress-related terms $\epsilon_{{\mathbf{P}}}\mathbf{v}$ and $\mathbf{q}_{{\mathbf{P}}}$ are both due to the enthalpic term in the extended Hamiltonian. The former has a convective nature; the latter, whose form is

[TABLE]

accounts for the heat flux generated by the action of the prescribed macroscopic Piola stress ${\mathbf{P}}$.

III.2.4 About the use of the Parrinello–Rahman kinetic energy

The Lagrangian kinetic energy introduced by Parrinello & RahmanParrinello and Rahman (1980); *PR2 is

[TABLE]

Since $K_{\mathrm{PR}}$ is a symmetric and positive definite quadratic form in the variables $((\dot{\mathbf{s}}_{k})_{k=1}^{N},\dot{\mathbf{F}})$, the above analysis invoving $K_{\mathrm{L}}$ is applicable also in this case, with similar results. What changes is not the form of the resulting continuum balances but the physical significance of the terms in those balances that depend directly on the microscopic evolution equations derived from the choice of $K_{\mathrm{PR}}$ for the kinetic energy.

It is worth observing that the inertial parameter $W$ plays a role akin to that of $Q$ in the Nosé–Hoover Hamiltonian. This marks an important distinction between the extension pertinent to the Nosé–Hoover Hamiltonian, which would be paralleled by one based on the Parrinello–Rahman Lagrangian, and the extension realized in the APR Hamiltonian discussed in this paper. When the physical kinetic energy is translated in terms of the extended variables, the outcome is $K_{L}$ and there is no need for inertial parameters, whose value is regarded as adjustable in more than one way.Andersen (1980); Parrinello and Rahman (1981) In this connection, we also observe that the difference between $K_{\mathrm{L}}$ and $K_{\mathrm{PR}}$ vanishes if two constraints are imposed on the motion of the computational cell and the particles therein,Podio-Guidugli (2010) one being holonomic and of the form

[TABLE]

where $\mathbf{I}$ is the identity tensor, and the other being non-holonomic and of the form

[TABLE]

It would be of some interest to introduce an extended Hamiltonian in which ad hoc Lagrangian multipliers would mediate the presence of these constraints.

IV Conclusions

We discussed the application of the Irving–Kirkwood–Noll procedure to extended Hamiltonian systems employed in certain molecular dynamics simulations. This procedure makes it possible to obtain macroscopic balance laws of continuum mechanics that can be regarded as specific collective counterparts of the extended Hamiltonian dynamics of a system of material points.

Our analysis demonstrates that the presence of auxiliary degrees of freedom, needed to control the statistical properties at the atomistic level, affects the structure of the balance equations at the continuum level. Precisely, the modifications induced at the microscopic scale influence the thermodynamic and transport properties at the macroscopic scale by way of novel source and flux terms in the continuum balances for the mass and energy densities.

We considered in detail two important examples of extended Hamiltonian systems. First, we discussed the application of the Irving–Kirkwood–Noll procedure to the extended Hamiltonian underlying the Nosé–Hoover thermostat. Notably, the effective thermal inertia present in that model produces source terms in the balances of both mass and energy, source terms which do not feature in standard continuum models for isothermal processes. Second, we discussed the application of the Irving–Kirkwood–Noll procedure to an extended Hamiltonian that we designated after Andersen, Parrinello, and Rahman because they were first to introduce extended Lagrangians to control the pressure and stress imposed on a collection of material points during a molecular dynamics simulation. In this context, it is possible to identify a contribution to the heat flux that is absent in the classical context. The term in question encodes the effects on the macroscopic balance of the energy needed to keep a constant the stress applied to the system.

Our findings are particularly relevant when considering computational schemes aiming at connecting the atomistic and continuum scales. Indeed, we view it as essential to correctly subsume the statistical properties imposed on the microscopic model on the thermodynamic and transport properties of the coupled continuum. When developing multiscale numerical schemes that aim at coupling, in a physically consistent manner, statistical information retrieved from MD with continuum modeling, this goal can be achieved by adopting continuum mechanical balances of the type we derived. We surmise that a failure to do so may lead to undesirable inconsistencies and computational artifacts.

Acknowledgments

The authors wish to acknowledge the support of CECAM – Centre Européen de Calcul Atomique et Moléculaire through the organization of a workshop that fostered the development of the present study. E. F. and G. G. G. gratefully acknowledge support from the Okinawa Institute of Science and Technology Graduate University with subsidy funding from the Cabinet Office, Government of Japan.