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é, 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 the 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.

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 and Kirkwood1 in a paper that appeared in 1950. Their path was completed a few years later by Noll2 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 is conserved.

An innovative and effective technique for accessing different statistical ensembles was introduced by Andersen3 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 and Rahman4 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é,5 their papers revolutionized the approach to investigating structural phase transitions. A few years later, Nosé,6 inspired once again by Andersen’s paper, conceived the idea of an extended Hamiltonian befitting the MD simulations consistent with a constant-temperature canonical ensemble (NVT), in the sense that “the canonical distribution is realized in a physical system.”5 Simulations of extended Hamiltonian dynamics were also considered within the more general class of non-Hamiltonian MD in important articles on their statistical properties.7 A development that, to our knowledge, went unnoticed in the literature is that Noll,2 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.5,6,8 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 and Rahman,4 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,3 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 Murdoch9 (who proposes an identification procedure alternative to that used in the IKN procedure along lines anticipated by Hardy10), Admal and Tadmor,11 Seguin and Fried,12 and others.13–20 

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 and Fried12 and Davydov and Steinmann.20 Nevertheless, in the present treatment, we restrict attention to the 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,21,22 which makes it directly relevant to certain recently proposed multiscale computational schemes.23–25 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 and Morriss.26 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 Sec. 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 Sec. 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 Sec. IV.

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 Z. The evolution of the system is given by Hamilton’s equations

ζ̇(t)={z,H}|z=ζ(t),
(1)

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

{A,B}:=AqBpApBq.
(2)

Given a Hamiltonian system, it is readily checked that

divz{z,H}=0,
(3)

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 Z is locally volume-preserving.

The incompressibility condition expressed by (3) has an important consequence: any probability density function ρ̃ defined on Z at the initial time t0 is convected along the Hamiltonian trajectories in the sense that

ρ̃̇(ζ(t),t)=0,
(4)

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

tρ̃(z,t){H(z,t),ρ̃(z,t)}=0,
(5)

where

[]:={H(z,t),}
(6)

is the Liouville operator for a system of Hamiltonian H. Liouville’s equation evolves a given initial value ρ̃(z,t0)=ρ̃0(z) of the probability measure at any given point of Z into its current value ρ̃(z,t) at that same point. The kernel of the Liouville operator 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.

If we consider a decomposition Z=Zs×Ze of the phase space Z, where Zs is associated with a subset of Hamiltonian coordinates, denoted by (qs, ps), and Ze 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

{A,B}Zs:=AqsBpsApsBqs;
(7)

the complementary bracket {A,B}Ze has an analogous expression. With such definitions, we obtain an additive decomposition of the Poisson bracket:

{A,B}={A,B}Zs+{A,B}Ze.
(8)

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

ρ̃t={H,ρ̃}Zs+{H,ρ̃}Ze.
(9)

A microscopic observable is a mapping zô(z) defined over the phase space Z. At this stage, the tensorial order of ô can be left arbitrary. The corresponding macroscopic observable is the expected value of ô at time t with respect to the chosen probability density ρ̃(z,t), namely

o(t):=Zρ̃(z,t)ô(z)dz.
(10)

An important application of Liouville’s equation is to compute the time derivative of the expected value o of a microscopic observable ô, which is given by

ot=Z{H,ρ̃}ôdz.
(11)

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

ot=Z{H,ρ̃}Zsôdz+Z{H,ρ̃}Zeôdz.
(12)

Irving and Kirkwood1 (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).

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

q=(r1,,rN)andp=(m1ṙ1,,mNṙN),
(13)

where rk is the current position vector, with respect to a fixed origin, of the k-th material point, the mass of which is mk. 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

ô(z,r)=ko^k(z)δ(rkr),
(14)

where 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

o(r,t)=ko^k(z)δ(rkr)=Zô(z,r)ρ̃(z,t)dz
(15)

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.

The first system to which we apply the IKN procedure is the extended Hamiltonian system introduced by Nosé,5,6 which forms the basis for the thermostatting strategy proposed by Hoover.8 We shall see that 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.

1. The extended Hamiltonian

Nosé’s extended Hamiltonian can be written as

HNH(q,p;Q,T):=k|pk|22mks2+ps22Q+V(q)+(3N+1)kBT(lns1),
(16)

the underlying coordinates being

q=(q,s)andp=(p1,,pN,ps).
(17)

The string q is the standard string of Lagrangian coordinates in (13)1, augmented by a dimensionless scalar coordinate s; ps, the momentum conjugate to s, augments a string of virtual momenta pk,k{1,,N}. For the physical velocity vk of the k-th material point, the relationship between the virtual momentum pk and the physical momentum mkvk is mediated by s as follows:

pks=mkvk.
(18)

The role assigned by Nosé to s is to scale the physical time t in terms of the virtual time τ involved in the Hamiltonian evolution associated with HNH according to

t=0τdαs(α).
(19)

This interpretation of s, which depends on the virtual time, is consistent with Nosé’s intention to achieve temperature control by controlling particle velocities and granted the understanding that the temperature of a system is related to the average of its kinetic energy.5 For this formulation to make sense, s must be positive. Consistent with this observation, we stipulate that the probability density function ρ̃(q,p,τ) must vanish identically for s0.

The expression of HNH 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 kB is Boltzmann’s constant. Specifically, Q modulates the new kinetic term proportional to ps2 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)kBT(lns1) 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 = Vi + Ve, where the internal contribution

Vi(q):=12jkVjk(|rjrk|),
(20)

accounts for pairwise interactions between material points and

Ve(q):=kVke(rk)
(21)

is the unary potential of external forces.

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 τ; in what follows, the inversion of relation (19) is implied. Specifically, we define the mass density

ρ(r,t)kmkδ(rkr),
(22)

the linear momentum density

ρv(r,t)kpksδ(rkr)=kmkvkδ(rkr),
(23)

the kinetic energy density

𝜖K(r,t)k|pk|22mks2δ(rkr)=kmk2|vk|2δ(rkr),
(24)

and the potential energy density

𝜖V(r,t)jk12Vjk(|rjrk)|δ(rkr)+kVke(rk)δ(rkr).
(25)

The presence of additional degrees of freedom and energy terms in HNH 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 r, and corresponding distributed observables (indicated with superposed bars) of the IKN form.

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

𝜖ps(t)1ωrefps22Q,
(26)
𝜖s(t)Aωref(lns1),
(27)

where we have introduced A:= (3N + 1)kBT and ω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

k1Nδ(rkr)
(28)

of the number density

n(r,t)k1Nδ(rkr).
(29)

Following that prescription, we define the distributed observables

𝜖¯ps(r,t)ps22Qk1Nδ(rkr),
(30)
𝜖¯s(r,t)A(lns1)k1Nδ(rkr).
(31)

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.

3. The mass and momentum balances

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

ρ̃τ=kHNHrkρ̃pkρ̃rkHNHpk+HNHsρ̃psρ̃sHNHps,
(32)

we can deduce the macroscopic balances of the defined observables.

The time derivative of the mass density is

ρt=Zsρ̃τkmkδ(rkr)dz=divrZρ̃kpksδ(rkr)dz+Zρ̃psQkmkδ(rkr)dz=divr(ρv)+σρ,
(33)

where the source term is given by

σρ(r,t)psQkmkδ(rkr).
(34)

As for the linear momentum density, we find that

(ρv)t=Zρ̃kVrkδ(rkr)dzdivrZρ̃kmkvkvkδ(rkr)dz=divr(ρvvT)+fe,
(35)

where the Cauchy stress tensor T is the sum of the kinetic stress tensor

TK(r,t)kmk(vkv)(vkv)δ(rkr)
(36)

and the stress tensor associated with internal interactions, namely,

TV(r,t)12jkR3xx|x|Vjk(|x|)J(x,r,t)dx,
(37)

with

J(x,r,t):=01δ(rjrαx)δ(rkr+(1α)x)dα,
(38)

and where the external force field is

f e(r,t)kVkerkδ(rkr).
(39)

The detailed derivation of formulae (37) and (38) was provided originally by Noll,2 who more recently generalized and refined his calculations.27 Note that the very presence of the source term σρ in the mass balance, a term that has no classical counterpart, is dictated by the form (16) of the extended Hamiltonian HNH, which however does not give rise to specific new terms in the balance of linear momentum.

4. Energy balances

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

𝜖=𝜖K+𝜖V+𝜖ps+𝜖s
(40)

and

𝜖¯=𝜖K+𝜖V+𝜖¯ps+𝜖¯s.
(41)

In both instances, the time derivatives of 𝜖K and 𝜖V produce the counterpart of the continuum mechanical heat flux q, which involves the sum of a kinetic contribution

qK(r,t)kmk2|vkv|2(vkv)δ(rkr),
(42)

an interaction contribution

qV(r,t)12jkR3x|x|Vjk(|x|)x0112(vj+vk)×δ(rjrαx)δ(rkr+(1α)x)dαdx,
(43)

and a transport contribution

qT(r,t)jk12Vjk(|rjrk|)(vkv)δ(rkr)+kVke(rk)(vkv)δ(rkr).
(44)
a. Energy balance with collective observables.

The time derivatives of the collective energy densities 𝜖s and 𝜖ps defined in (26) and (27) are

d𝜖psdt=ps32Q2ωref+psQωrefkmk|vk|2Aσ𝜖ps
(45)

and

d𝜖sdt=psQωrefAlnsσ𝜖s;
(46)

these derivatives induce two uniform source terms, σ𝜖ps and σ𝜖s.

All in all, the balance of the energy density 𝜖 defined in (40) takes the form

𝜖t=divr(q+𝜖Kv+𝜖VvTTv)+σ𝜖,
(47)

where q and σ𝜖 are given by

q=qK+qV+qT,
(48)
σ𝜖=σ𝜖0+σ𝜖K+σ𝜖V+σ𝜖ps+σ𝜖s.
(49)

In addition to the time-dependent but spatially uniform sources σ𝜖ps and σ𝜖s, σ𝜖 features the classical term

σ𝜖0(r,t)kVkerkvkδ(rkr)
(50)

and two other new terms, namely

σ𝜖K(r,t)psQkmk2|vk|2δ(rkr)
(51)

and

σ𝜖V(r,t)psQjk12Vjk(|rjrk)|δ(rkr)+psQkVke(rk)δ(rkr).
(52)

While the time evolution of 𝜖K and 𝜖V produces spatially dependent diffusive fluxes, convective terms, and sources, the time evolution of 𝜖ps and 𝜖s generates only uniform sources.

b. Energy balance with distributed observables.

To deduce the balance of the energy density 𝜖¯ defined in (41), we compute the time derivatives of the distributed observables 𝜖¯ps and 𝜖¯s defined in (30) and (31). Specifically, we find that

𝜖¯pst=divr(qps+𝜖¯psv)+σ𝜖¯ps,
(53)

with

qps(r,t)ps22Qk1N(vkv)δ(rkr)
(54)

and

σ𝜖¯ps(r,t)ps32Q2k1Nδ(rkr)+psQjmj|vj|2Ak1Nδ(rkr).
(55)

Additionally, we find that

𝜖¯st=divr(qs+𝜖¯sv)+σ𝜖¯s,
(56)

with

qs(r,t)A(lns1)k1N(vkv)δ(rkr)
(57)

and

σ𝜖¯s(r,t)psQAlnsk1Nδ(rkr).
(58)

All in all, the balance of the energy density 𝜖¯ reads

𝜖¯t=divr(q+𝜖¯vTTv)+σ𝜖¯,
(59)

where q and σ𝜖¯ are given by

q=qK+qV+qT+qps+qs,
(60)
σ𝜖¯=σ𝜖0+σ𝜖K+σ𝜖V+σ𝜖¯ps+σ𝜖¯s.
(61)

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

In this section, we apply our generalized IKN procedure to an extended Hamiltonian system akin to but more general than that considered by Parrinello and Rahman4 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 schemes24,25 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 and Rahman under specific conditions analyzed by Podio-Guidugli.21 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 AB, with a single centered dot, the scalar product defined through the trace operator by tr(ATB).

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 F, with positive determinant, which maps the referential position vector sk of the k-th material point into its current position vector rk:

rk=Fsk.
(62)

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

ṙk=Fṡk+Ḟsk=vk.
(63)

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:

KL:=12kmk|Fṡk+Ḟsk|2=12FTFkmkṡkṡk+12ḞTḞkmksksk+14(ḞTF+FTḞ)kmk(ṡksk+skṡk).
(64)

The Lagrangian kinetic energy KL is a symmetric and positive semi-definite quadratic form in the variables ((ṡk)k=1N,Ḟ), which depends parametrically on q. The corresponding Hamiltonian kinetic energy K is the Legendre–Fenchel transform of KL. Although KL depends on sk and ṡk, it is independent of rk and, therefore, so is K. We thus have

pk=KLṡk,G=KLḞ,ṡk=Kpk,and Ḟ=KG.
(65)

The conjugate momentum pk and the physical momentum mkvk of the k-th material point are linked by the transformation

pk=mkFTṙk=mkFT(Fṡk+Ḟsk),
(66)

from which we deduce the alternative expression FTpk for the physical momentum of the k-th material point.

The APR extended Hamiltonian that we consider is

HAPR(q,p)=K(q,p)+V(q)ωrefPF,
(67)

with p=((pk)k=1N,G). Here, as usual, the potential energy depends on the pairwise interactions of material points in the positions they currently occupy:

V(q):=jk12Vjk(|F(sjsk)|).
(68)

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), ωrefPF, which is of an enthalpic character. In that term, ωref>0 denotes the constant reference volume of the computational cell, while the time-independent control parameter 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 ρ̃(q,p,t) is determined by Liouville’s equation (9), which takes the following explicit form:

ρ̃t=kHAPRskρ̃pkρ̃skHAPRpk+HAPRFρ̃Gρ̃FHAPRG.
(69)

2. Definition of the observables

The discussion in Sec. III B 1 provides us with precise indications about the appropriate form for the various observables. We define the mass density

ρ(r,t)kmkδ(Fskr),
(70)

the linear momentum density

ρv(r,t)kFTpkδ(Fskr)=kmkvkδ(Fskr),
(71)

the kinetic energy density

𝜖K(r,t)k12mk|FTpk|2δ(Fskr)=kmk2|vk|2δ(Fskr),
(72)

and the potential energy density

𝜖V(r,t)jk12Vjk(|F(sjsk)|)δ(Fskr).
(73)

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

𝜖P(r,t)PF=PF,
(74)

and its distributed counterpart

𝜖¯P(r,t)ωref(PF)k1Nδ(Fskr).
(75)

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 T = TK + TV is obtained from (36)–(38) modulo the substitutions

rkFsk,
(76)
pksmk(Fṡk+Ḟsk)=FTpk=mkvk.
(77)

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

𝜖t=divr(q+𝜖Kv+𝜖VvTTv)+σ𝜖,
(78)

with

𝜖=𝜖K+𝜖V+𝜖P,
(79)
q=qK+qV+qT,
(80)

and

σ𝜖𝜖̇P.
(81)

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 ωrefF constitutes an extensive variable and σ𝜖 provides a reckoning of the power expenditure associated with changing the shape of the computational cell under the influence of the Piola stress P.28,29

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

𝜖¯t=divr(q+𝜖¯vTTv),
(82)

where

𝜖¯=𝜖K+𝜖V+𝜖¯p
(83)

and

q=qK+qV+qT+qP.
(84)

The stress-related terms 𝜖¯Pv and qP are both due to the enthalpic term in the extended Hamiltonian. The former has a convective nature; the latter, whose form is

qPωref(PF)kvkvNδ(Fskr),
(85)

accounts for the heat flux generated by the action of the prescribed macroscopic Piola stress P.

4. About the use of the Parrinello–Rahman kinetic energy

The Lagrangian kinetic energy introduced by Parrinello and Rahman4 is

KPR:=12FTFkmkṡkṡk+12W|Ḟ|2,W>0.
(86)

Since KPR is a symmetric and positive definite quadratic form in the variables ((ṡk)k=1N,Ḟ), the above analysis involving KL is applicable also in this case, with similar results. What changes is not the form of the resulting continuum balances but rather is the physical significance of the terms in those balances that depend directly on the microscopic evolution equations derived from the choice of KPR 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 KL and there is no need for inertial parameters, whose value is regarded as adjustable in more than one way.3,4 In this connection, we also observe that the difference between KL and KPR vanishes if two constraints are imposed on the motion of the computational cell and the particles therein,21 one being holonomic and of the form

kmksksk=WI,
(87)

where I is the identity tensor, and the other being non-holonomic and of the form

ḞF1=(ḞF1)T.
(88)

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

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

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.

1.
J. H.
Irving
and
J. G.
Kirkwood
, “
The statistical mechanical theory of transport processes. IV. The equations of hydrodynamics
,”
J. Chem. Phys.
18
,
817
829
(
1950
).
2.
W.
Noll
, “
Die herleitung der grundgleichungen der thermomechanik der kontinua aus der statistischen mechanik
,”
Indiana. Univ. Math. J.
4
,
627
646
(
1955
);
R. B.
Lehoucq
and
A.
Von Lilienfeld-Toal
, “
Translation of Walter Noll’s derivation of the fundamental equations of continuum thermodynamics from statistical mechanics
,”
J. Elasticity
100
,
5
24
(
2010
).
3.
H. C.
Andersen
, “
Molecular dynamics simulations at constant pressure and/or temperature
,”
J. Chem. Phys.
72
,
2384
2393
(
1980
).
4.
M.
Parrinello
and
A.
Rahman
, “
Polymorphic transitions in single crystals: A new molecular dynamics method
,”
J. Appl. Phys.
52
,
7182
7190
(
1981
).
5.
S.
Nosé
, “
Constant temperature molecular dynamics methods
,”
Prog. Theor. Phys. Suppl.
103
,
1
46
(
1991
).
6.
S.
Nosé
, “
A molecular dynamics method for simulations in the canonical ensemble
,”
Mol. Phys.
52
,
255
262
(
1984
);
S.
Nosé
, “
A unified formulation of the constant temperature molecular dynamics methods
,”
J. Chem. Phys.
81
,
511
519
(
1984
).
7.
M. E.
Tuckerman
,
C. J.
Mundy
, and
G. J.
Martyna
, “
On the classical statistical mechanics of non-Hamiltonian systems
,”
Europhys. Lett.
45
,
149
(
1999
);
M. E.
Tuckerman
,
Y.
Liu
,
G.
Ciccotti
, and
G. J.
Martyna
, “
Non-Hamiltonian molecular dynamics: Generalizing Hamiltonian phase space principles to non-Hamiltonian systems
,”
J. Chem. Phys.
115
,
1678
1702
(
2001
).
8.
W. G.
Hoover
, “
Canonical dynamics: Equilibrium phase-space distributions
,”
Phys. Rev. A
31
,
1695
1697
(
1985
);
W. G.
Hoover
, “
Nosé–Hoover nonequilibrium dynamics and statistical mechanics
,”
Mol. Simul.
33
,
13
19
(
2007
).
9.
A. I.
Murdoch
, “
On molecular modelling and continuum concepts
,”
J. Elasticity
100
,
33
61
(
2010
);
A. I.
Murdoch
,
Physical Foundations of Continuum Mechanics
(
Cambridge University Press
,
Cambridge
,
2012
), p.
xiv+423
.
10.
R. J.
Hardy
, “
Formulas for determining local properties in molecular dynamics simulations: Shock waves
,”
J. Chem. Phys.
76
,
622
628
(
1982
).
11.
N. C.
Admal
and
E. B.
Tadmor
, “
A unified interpretation of stress in molecular systems
,”
J. Elasticity
100
,
63
143
(
2010
);
N. C.
Admal
and
E. B.
Tadmor
, “
Stress and heat flux for arbitrary multibody potentials: A unified framework
,”
J. Chem. Phys.
134
,
184106
(
2011
);
[PubMed]
N. C.
Admal
and
E. B.
Tadmor
, “
Material fields in atomistics as pull-backs of spatial distributions
,”
J. Mech. Phys. Solids
89
,
59
76
(
2016
).
12.
B.
Seguin
and
E.
Fried
, “
Statistical foundations of liquid-crystal theory. I: Discrete systems of rod-like molecules
,”
Arch. Ration. Mech. Anal.
206
,
1039
1072
(
2012
);
[PubMed]
B.
Seguin
and
E.
Fried
, “
Statistical foundations of liquid-crystal theory II: Macroscopic balance laws
,”
Arch. Ration. Mech. Anal.
207
,
1
37
(
2013
).
[PubMed]
13.
E.
Wajnryb
,
A. R.
Altenberger
, and
J. S.
Dahler
, “
Uniqueness of the microscopic stress tensor
,”
J. Chem. Phys.
103
,
9782
9787
(
1995
).
14.
C.
Denniston
and
M. O.
Robbins
, “
Mapping molecular models to continuum theories for partially miscible fluids
,”
Phys. Rev. E
69
,
021505
(
2004
);
C.
Denniston
and
M. O.
Robbins
, “
General continuum boundary conditions for miscible binary fluids from molecular dynamics simulations
,”
J. Chem. Phys.
125
,
214102
(
2006
).
[PubMed]
15.
I.
Goldhirsch
, “
Stress, stress asymmetry and couple stress: From discrete particles to continuous fields
,”
Granular Matter
12
,
239
252
(
2010
).
16.
Y. A.
Humenyuk
and
M. V.
Tokarchuk
, “
Extension of hydrodynamic balance equations for simple fluids
,”
J. Stat. Phys.
142
,
1052
1084
(
2011
).
17.
R. B.
Lehoucq
and
M. P.
Sears
, “
Statistical mechanical foundation of the peridynamic nonlocal continuum theory: Energy and momentum conservation laws
,”
Phys. Rev. E
84
,
031112
(
2011
).
18.
J. Z.
Yang
,
X.
Wu
, and
X.
Li
, “
A generalized Irving–Kirkwood formula for the calculation of stress in molecular dynamics models
,”
J. Chem. Phys.
137
,
134104
(
2012
).
19.
G.
Capriz
and
P. M.
Mariano
, “
Objective fluxes in a multi-scale continuum description of sparse medium dynamics
,”
Phys. A
415
,
354
365
(
2014
).
20.
D.
Davydov
and
P.
Steinmann
, “
Reviewing the roots of continuum formulations in molecular systems. Part I: Particle dynamics, statistical physics, mass and linear momentum balance equations
,”
Math. Mech. Solids
19
,
411
433
(
2014
);
D.
Davydov
and
P.
Steinmann
, “
Reviewing the roots of continuum formulations in molecular systems. Part II: Energy and angular momentum balance equations
,”
Mathematics Mech. Solids
19
,
852
867
(
2014
);
D.
Davydov
and
P.
Steinmann
, “
Reviewing the roots of continuum formulations in molecular systems. Part III: Stresses, couple stresses, heat fluxes
,”
Mathematics Mech. Solids
20
,
1153
1170
(
2015
).
21.
P.
Podio-Guidugli
, “
On (Andersen–)Parrinello–Rahman molecular dynamics, the related metadynamics, and the use of the Cauchy–Born rule
,”
J. Elasticity
100
,
145
153
(
2010
).
22.
P.
Podio-Guidugli
, “
On microscopic and macroscopic notions of stress
,” in
Elasticity and Inelasticity
(
Moscow University Supercomputing Center
,
2011
), pp.
279
285
.
23.
D.
Davydov
,
J.-P.
Pelteret
, and
P.
Steinmann
, “
Comparison of several staggered atomistic-to-continuum concurrent coupling strategies
,”
Comput. Methods Appl. Mech. Eng.
277
,
260
280
(
2014
).
24.
M. H.
Ulz
, “
A multiscale molecular dynamics method for isothermal dynamic problems using the seamless heterogeneous multiscale method
,”
Comput. Methods Appl. Mech. Eng.
295
,
510
524
(
2015
).
25.
S.
Li
and
Q.
Tong
, “
A concurrent multiscale micromorphic molecular dynamics
,”
J. Appl. Phys.
117
,
154303
(
2015
);
Q.
Tong
and
S.
Li
, “
From molecular systems to continuum solids: A multiscale structure and dynamics
,”
J. Chem. Phys.
143
,
064101
(
2015
).
[PubMed]
26.
D. J.
Evans
and
G.
Morriss
,
Statistical Mechanics of Nonequilibrium Liquids
, 2nd ed. (
Cambridge University Press
,
2008
).
27.
W.
Noll
, “
Thoughts on the concept of stress
,”
J. Elasticity
100
,
25
32
(
2010
).
28.
A.
McLellan
,
The Classical Thermodynamics of Deformable Materials
, Cambridge Monographs on Physics (
Cambridge University Press
,
1980
).
29.
C.-S.
Man
, “
Classical thermodynamics of elastic solids as open systems
,”
J. Elasticity
126
,
271
280
(
2017
).