Concurrent atomistic-continuum modeling of crystalline materials

Chen, Youping; Shabanov, Sergei; McDowell, David L.

doi:10.1063/1.5099653

In this work, we present a concurrent atomistic-continuum (CAC) method for modeling and simulation of crystalline materials. The CAC formulation extends the Irving-Kirkwood procedure for deriving transport equations and fluxes for homogenized molecular systems to that for polyatomic crystalline materials by employing a concurrent two-level description of the structure and dynamics of crystals. A multiscale representation of conservation laws is formulated, as a direct consequence of Newton's second law, in terms of instantaneous expressions of unit cell-averaged quantities using the mathematical theory of distributions. Finite element (FE) solutions to the conservation equations, as well as fluxes and temperature in the FE representation, are introduced, followed by numerical examples of the atomic-scale structure of interfaces, dynamics of fracture and dislocations, and phonon thermal transport across grain boundaries. In addition to providing a methodology for concurrent multiscale simulation of transport processes under a single theoretical framework, the CAC formulation can also be used to compute fluxes (stress and heat flux) in atomistic and coarse-grained atomistic simulations.

I. INTRODUCTION

Solid state physics describes the structure of all crystals in terms of primitive unit cells. These primitive unit cells, each containing an identical group of atoms, are the building blocks of crystals. When added continuously in space, they form the structure of a crystal and completely fill the space the crystal occupies,¹ cf. Fig. 1.

FIG. 1.

View large Download slide

A two-level structural description of crystals: a continuously distributed lattice + a discrete basis.

Although the motion of an atom in a crystal is an irregularly fluctuating function of time, it turns out that the dynamics of a crystal can be most readily described, not in terms of individual atoms, but in terms of traveling waves, named lattice vibrations by Max Born. These vibrational modes are called phonons and are of two types: acoustic and optical, cf. Fig. 2. With the former, atoms in a unit cell move in the same phase, resulting in the deformation of the lattice; with the latter, atoms undergo relative motion within the lattice cells, leaving the lattice unchanged. Thus, the displacements of atoms in a crystal can be described as the sum of continuous lattice displacements and the internal displacements of the atoms relative to the lattice.²

FIG. 2.

View large Download slide

Transverse atomic motion in (a) an acoustic mode and (b) an optical mode.

It is seen that the solid state physics description of the structure and dynamics of crystals is based on a concurrent atomistic-continuum (CAC) approach: the structure is continuous at the lattice level but discrete at the atomic level; the lattice deformation is continuous until a defect is generated, cf. Fig. 3, while the internal motion is discrete.

FIG. 3.

View large Download slide

Schematics illustrating that (a) the lattice deformation is continuous until (b) structural discontinuities generated by slip and/or twinning.

Historically, atomistic and continuum descriptions offer two fundamentally different approaches to our understanding of crystalline matter. From the atomistic viewpoint, matter consists of discrete particles; while from the continuum viewpoint, matter is infinitely divisible. The continuum view leads to formulations of field equations of conservation laws. Supplemented by constitutive relations, such as Hooke's law and Fourier's law, these conservation equations serve as the governing equations in continuum modeling of materials.

This work introduces the concurrent atomistic-continuum (CAC) method that links and unifies the atomistic and continuum descriptions of crystalline materials. The Irving-Kirkwood procedure for linking the molecular description to hydrodynamical equations is reviewed in Sec. II; the CAC formalism is introduced in Sec. III; the finite element (FE) formulation of CAC as well as the fluxes and temperature in the FE representation is outlined in Sec. IV; numerical examples are presented in Sec. V, followed by a summary in Sec. VI. Detailed derivations of the conservation equations using the mathematical theory of distributions for monatomic crystals are included in Appendix A and that for polyatomic crystals in Appendix B.

II. THE IRVING-KIRKWOOD PROCEDURE

The statistical mechanics formulation to link discrete and continuum descriptions of conservation laws was pioneered by Irving and Kirkwood (IK) in 1950.³ In their landmark paper, IK defined local densities in terms of molecular variables using the Dirac delta for “localized density,”⁴ “since mass or momentum of any molecule may be considered as localized at that molecule.”³ A local density of a dynamic phase function

A (r_{k}, v_{k})

⁠, where

r_{k}

and

v_{k}

are the position and velocity of the kth particle, is defined as an ensemble-averaged point function

\bar{a} (x, t)

as

\bar{a} (x, t) \equiv ⟨ \sum_{k = 1}^{N} A (r_{k}, v_{k}) δ (r_{k} - x) ⟩,

(1)

where the symbol “≡” means “equal by definition.”

Equation is usually interpreted as mapping the phase functions in the 6N-dimensional phase space to the local densities in the 4-dimensions of physical space and time. Based on the local density definition in Eq. , IK derived the rate of change for the densities of mass, momentum, and energy in the form of partial differential equations that relate the density of a conserved quantity

\bar{a} (x, t)

and its flux

\bar{J} (x, t)

as

\frac{\partial}{\partial t} \bar{a} (x, t) = \nabla_{x} \cdot \bar{J} (x, t) .

(2)

Formulas for fluxes (stress and heat flux) are then obtained through comparing the rate equations with the differential equations of hydrodynamics.

The IK formalism has inspired numerous efforts to link molecular variables to field quantities in classical continuum mechanics.^5–8 To understand the motivation of a new formalism, we recall the following noted by Irving and Kirkwood³ and by Kirkwood.⁹

Local densities and fluxes in the IK formulation are point functions. “The point functions, although averaged neither over space nor time, satisfy equations that are identical in form to the equations of hydrodynamics for a single component, single phase system.”³ These equations were called “the hydrodynamical-like equations” in the IK paper.³ It was noted that “to obtain the hydrodynamical equations themselves it is merely necessary to perform the appropriate space and time averages.”³
The densities for mass, momentum, and energy were defined without restriction to a single component single phase system. For fluxes (stress and heat flux), IK stated that “we must impose this restriction” and the flux formulas obtained with their formalism are “only valid for a single component, single phase system.”³
The statistical mechanical theory of transport processes published in a series of 14 papers by Kirkwood and co-workers only considers rigid molecules, ignoring the internal degrees of freedom (DOFs) of the molecules. In the first one of the series, Kirkwood envisioned extension of the formulation to molecules possessing internal degrees of freedom (DOFs).⁹

III. THE CAC FORMALISM

The CAC formalism is an extension of the IK procedure for homogeneous hydrodynamic systems to a two-level structural description of crystals that includes the atomic DOFs within each lattice cell.^10,11 With such an extension, the CAC formulation is applicable to polyatomic crystalline materials and can be used to solve atomic trajectories or quantify instantaneous fluxes, conferring full characteristics of a concurrent atomistic-continuum method.

A. Ensemble-averaged local densities

For a polyatomic crystal of

N_{l}

unit cells, with each unit cell containing

N_{α}

atoms, the density of a phase function

A (r_{k ξ}, v_{k ξ})

at point r in the physical space can be expressed as

\bar{a} (r, t) = ⟨ \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} A (r_{k ξ}, v_{k ξ}) δ (r_{k ξ} - r) ⟩,

(3)

where

r_{k ξ}

and

v_{k ξ}

are the position vector and velocity vector of the ξth atom in the kth unit cell, and the symbol ⟨ ⟩ denotes ensemble averaging that involves

N_{l} N_{a}

fold integrals of the phase function A and a probability distribution function in the phase space.

Recall that the structure and deformation of a crystal are continuous in space at the lattice level but are discrete at the atomic level. It is, therefore, consistent with solid state physics to describe the crystal structure in terms of the positions of lattice cells and the relative positions of atoms inside the lattice cells. For this purpose, we resolve the sum in the definition of the local density in Eq. into contributions from groups of atoms within the same lattice cell and express the locations relative to the unit cell for those atoms, i.e.,

\bar{a} (x, y, t) \equiv ⟨ \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a} - 1} A (r_{k ξ}, v_{k ξ}) δ (x - r_{k}) δ (y - Δ r_{k ξ}) ⟩,

(4)

where

r_{k}

is the position of kth unit cell and

Δ r_{k ξ}

is the relative position of ξth atom in the unit cell k.

The density in Eq. is expressed in terms of the direct product of two deltas, with the two field variables x and y representing a distinction of large- and small-scale variations of the variable r in Eq. . Note that the ensemble average in Eq. also involves

N_{l} N_{a}

fold integrals and, although expressed in terms of different phase variables, the values of a local density defined in Eqs. and can be mathematically equivalent. However, different from that in Eq. , a conservation equation for the density defined in Eq. relates a conserved quantity to two fluxes,

{\bar{J}}_{1}

and

{\bar{J}}_{2}

⁠, owing to the use of large- and small-scale variations in the latter, i.e.,

\frac{\partial}{\partial t} \bar{a} (x, y, t) = \nabla_{x} \cdot {\bar{J}}_{1} (x, t) + \nabla_{y} \cdot {\bar{J}}_{2} (x, y, t) .

(5)

B. Space and time-averaged local densities

Evans and Morriss argued that, since conservation laws hold instantaneously, conservation equations and fluxes should be definable without ensemble averaging.¹² This has been demonstrated recently for monatomic crystals.^13–16 To derive the conservation equations and flux expressions that hold instantaneously based on the two-level structural description of crystals, we consider a unit cell of volume V located at point x in the physical space; within the unit cell, there are

N_{α}

volume elements, each containing one atom and having volume V_α satisfying

\sum_{α} V_{α} = V

⁠. Since the distribution of the unit cells is continuous and homogeneous in space, we may define a local density per unit-cell volume as

a (x, t) \equiv \sum_{k = 1}^{N_{l}} (\sum_{ξ = 1}^{N_{a}} A (r_{k ξ}, v_{k ξ})) δ_{V} (x - r_{k}) = \sum_{α = 1}^{N_{a}} a_{α} (x, t),

(6)

where

a_{α} (x, t) = \sum_{k} A (r_{k α}, v_{k α}) δ_{V} (x - r_{k})

is the contribution of the αth atom to

a (x, t)

and can be expressed as

a_{α} (x, t) = \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} A (r_{k ξ}, v_{k ξ}) δ_{V} (x - r_{k}) δ_{V α} (y - Δ r_{k ξ}) = a (x, y, t)

(7)

and

δ_{V} (x - r_{k}) \equiv \frac{1}{V} {\begin{matrix} 1 & if x - r_{k} \in V, \\ 0 & if x - r_{k} \notin V, \end{matrix}

(8)

δ_{V α} (y - Δ r_{k ξ}) \equiv {\begin{matrix} 1 & if y - Δ r_{k ξ} \in V_{α} or ξ = α, \\ 1 & if y - Δ r_{k ξ} \notin V_{α} or ξ \neq α . \end{matrix}

(9)

Since $\sum_{α = 1}^{N_{α}} δ_{V α} (y - Δ r_{k ξ}) = 1$ ⁠, one can see that Eq. (7) satisfies Eq. (6), and $a (x, y, t) = a_{α} (x, t)$ holds for all $y \in V_{α}$ ⁠.

Note that the box functions defined in Eqs. (8) and (9) have a jump discontinuity at the enclosing surfaces V and V_α, respectively. These functions do not have derivatives in the classical sense. It thus requires a generalization of the concept of a function. The corresponding mathematical theory is known as the theory of distributions.¹⁷ This theory will be used throughout this work as the mathematical tool to derive the conservation laws for physical quantities described by nondifferentiable functions.

In atomistic simulations, the equation of motion is solved in discrete time steps. Accordingly, local densities should be further averaged over a time-interval. The contribution of the αth atom to the mass density ρ_α, momentum density

p_{α} (p_{α} = ρ_{α} v_{α})

⁠, and energy density

E_{α} (E_{α} = ρ_{α} e_{α})

per unit volume can then be defined as

ρ_{α} (x, t) \equiv \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} m_{ξ} {\bar{δ}}_{V} (x - r_{k}) {\bar{δ}}_{V α} (y - Δ r_{k ξ}),

(10)

p_{α} (x, t) \equiv \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} m_{ξ} v_{k ξ} {\bar{δ}}_{V} (x - r_{k}) {\bar{δ}}_{V α} (y - Δ r_{k ξ}),

(11)

E_{α} (x, t) \equiv \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a} - 1} (\frac{1}{2} m_{ξ} v_{k ξ}^{2} + Φ_{k ξ}) {\bar{δ}}_{V} (x - r_{k}) {\bar{δ}}_{V α} (y - Δ r_{k ξ}),

(12)

where

Φ_{k ξ}

is the potential energy of atom kξ;

{\bar{δ}}_{V}

denotes the average of δ_V over a time-interval Δt, i.e.,

{\bar{δ}}_{V} (r_{k} - x) = \frac{1}{Δ t} \int_{0}^{Δ t} δ_{V} [r_{k} (t + τ) - x] d τ

⁠. The time evolution of a conserved quantity, i.e., the conservation equation, can then be expressed in terms of a lattice-level flux J as

\frac{\partial}{\partial t} a (x, t) = \frac{1}{V} \underset{\partial V}{\int \int ◯} J (x + x^{'}, t) \cdot n d S

(13)

or in terms of atomic-level fluxes J₁ and J₂ as

\begin{aligned} \frac{\partial}{\partial t} a_{α} (x, t) & = - \frac{1}{V} \underset{\partial V}{\int \int ◯} J_{1} (x + x^{'}, t) \cdot n d S \\ - \frac{1}{V} \underset{\partial V_{α}}{\int \int ◯} J_{2} (x, y + y^{'}, t) \cdot n d S_{α}, \end{aligned}

(14)

where n is the outward unit normal vector to the surface element, J₁ is the part of flux that flows across the enclosing surface ∂V, and J₂ is the part of flux that flows back and forth inside V and hence only cross the bounding surface ∂V_α, as shown in Fig. 4.

FIG. 4.

2D schematic for J1, the flux across the bounding surface of the unit cell (the rectangle) ∂V, and J2, the flux that does not cross ∂V but rather ∂Vα (the triangle) within the unit cell.

View large Download slide

2D schematic for J₁, the flux across the bounding surface of the unit cell (the rectangle) ∂V, and J₂, the flux that does not cross ∂V but rather ∂V_α (the triangle) within the unit cell.

It is noted that by integrating Eq. (5) over V and invoking the divergence theorem, one obtains an integral equation in the form of Eq. (14). However, Eq. (5) is not valid instantaneously and it requires the fluxes to be differentiable. By contrast, the conservation equation in Eq. (14) holds instantaneously in the distributional sense and the fluxes are not required to be continuous functions.

C. Conservation equations and fluxes

Averaging a local density over a unit cell of volume V leads to the integral form of conservation laws that equate the rate of the change of conserved quantities in the unit cell to their fluxes across the enclosing surface of the unit cell. Following from Eq. , the equation of conservation of mass can be derived in the form of Eq. as

\begin{aligned} \frac{\partial ρ_{α}}{\partial t} & = \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} m_{k ξ} \frac{\partial}{\partial t} {\bar{δ}}_{V} (x - r_{k}) {\bar{δ}}_{V α} (y - Δ r_{k ξ}) \\ + \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} m_{k ξ} {\bar{δ}}_{V} (x - r_{k}) \frac{\partial}{\partial t} {\bar{δ}}_{V α} (y - Δ r_{k ξ}) \\ = - \frac{1}{V} \underset{\partial V}{\int \int ◯} ρ_{α} v \cdot n d S - \frac{1}{V} \underset{\partial V_{α}}{\int \int ◯} ρ_{α} Δ v_{α} \cdot n d S_{α}, \end{aligned}

(15)

where

v = \dot{x}

⁠,

v_{α} = p_{α} / ρ_{α}

⁠, and

Δ v_{α} = v_{α} - v

⁠.

The momentum conservation law can also be derived, using the momentum density defined in Eq. , from Newton's second law, as

\begin{aligned} \frac{\partial (ρ_{α} v_{α})}{\partial t} = & \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} m_{ξ} {\dot{v}}_{k ξ} {\bar{δ}}_{V} (x - r_{k}) {\bar{δ}}_{V α} (y - Δ r_{k ξ}) \\ + \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} m_{ξ} v_{k ξ} \frac{\partial}{\partial t} {\bar{δ}}_{V} (x - r_{k}) {\bar{δ}}_{V α} (y - Δ r_{k ξ}) \\ + \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} m_{ξ} v_{k ξ} \frac{\partial}{\partial t} {\bar{δ}}_{V} (x - r_{k}) \frac{\partial}{\partial t} {\bar{δ}}_{V α} (y - Δ r_{k ξ}) \\ = & \frac{1}{V} \underset{\partial V}{\int \int ◯} (t_{α} - ρ_{α} v_{α} v \cdot n) d S + \frac{1}{V} \underset{\partial V_{α}}{\int \int ◯} (τ_{α} - ρ_{α} v_{α} Δ v_{α} \cdot n) d S_{α}, \end{aligned}

(16)

where t_α and

τ_{α}

are the momentum flux (stress) that cross ∂V and ∂V_α, respectively, as a result of atomic interaction and motion involving atom α. The total stress vector t on a surface element Aⁿ of area A and normal n, located at position y in the unit cell at x, from all the atoms can then be expressed as

t (x, y, n) = \sum_{α = 1}^{N_{a}} (t_{α} + τ_{α}) = t^{pot} (x, y, n) + t^{kin} (x, y, n),

(17)

where

t^{pot}

is the potential part of t that measures the interaction force per unit area transmitted across Aⁿ, i.e.,

\begin{aligned} t^{pot} (x, y, n) & = \sum_{k, l (\neq k)}^{N_{l}} \sum_{ξ, η}^{N_{a}} \frac{\partial Φ_{k ξ}}{\partial r_{l η}} \int_{r_{l η}}^{r_{k ξ}} {\bar{δ}}_{A}^{n} (φ - x) d φ \\ + \sum_{k}^{N_{l}} \sum_{ξ, η \neq ξ}^{N_{a}} V {\bar{δ}}_{V} (x - r_{k}) \frac{\partial Φ_{k ξ}}{\partial r_{k η}} \int_{r_{k η}}^{r_{k ξ}} {\bar{δ}}_{A}^{n} (φ - y) d φ, \end{aligned}

(18)

and

t^{kin}

is the kinetic part of t that measures the flow of momentum per unit area and time across Aⁿ, i.e.,

\begin{aligned} t^{kin} (x, y, n) = & - \sum_{k = 1}^{N_{l}} m {\tilde{v}}_{k} {\tilde{v}}_{k} {\bar{δ}}_{A}^{n} (x - r_{k}) \\ - \sum_{k = 1}^{N_{l}} V {\bar{δ}}_{V} (x - r_{k}) \sum_{α = 1}^{N_{a}} m_{α} Δ {\tilde{v}}_{k α} Δ {\tilde{v}}_{k α} {\bar{δ}}_{A}^{n} (y - Δ r_{k α}) . \end{aligned}

(19)

Here, ${\tilde{v}}_{k} = v_{k} - v$ ⁠, $Δ {\tilde{v}}_{k α} = Δ v_{k α} - Δ v_{α}$ ⁠, $\int_{r_{l η}}^{r_{k ξ}} {\bar{δ}}_{A}^{n} (φ - x) d φ$ ⁠, defined in Eq. (A3), and ${\bar{δ}}_{A}^{n} (y - Δ r_{k α})$ ⁠, defined in Eq. (A4), represent the potential and kinetic flux as a line-plane intersection problem in space and time, respectively.

Similarly, denoting q_α and j_α as the parts of heat flux that cross ∂V and ∂V_α, respectively, the time rate of change of the energy density E_α, can be derived as

\frac{\partial (ρ_{α} e_{α})}{\partial t} = \frac{\partial}{\partial t} \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} E_{k ξ} {\bar{δ}}_{V} (x - r_{k}) {\bar{δ}}_{V α} (y - Δ r_{k ξ}) = \frac{1}{V} \underset{\partial V}{\int \int ◯} (q_{α} + t_{α} \cdot v_{α} - v ρ_{α} e_{α}) \cdot n d S + \frac{1}{V} \underset{\partial V_{α}}{\int \int ◯} (j_{α} + τ_{α} \cdot v_{α} - Δ v_{α} ρ_{α} e_{α}) \cdot n d S_{α} .

(20)

The total heat flux is then given by

q (x, y) = \sum_{α = 1}^{N_{a}} (q_{α} + j_{α}) = q^{p o t} (x, y) + q^{k i n} (x, y) = q (x, y, e^{i}) e^{i},

(21)

where eⁱ (i = 1, 2, 3) are the orthonormal basis at x,

q_{}^{pot} (x, y, e^{i}) = \sum_{k, l \neq k}^{N_{l}} \sum_{α, η}^{N_{a}} \frac{\partial Φ_{k α}}{\partial r_{l η}} \cdot {\tilde{v}}_{k α} \int_{r_{l η}}^{r_{k α}} {\bar{δ}}_{A}^{i} (φ - x) d φ + \sum_{k}^{N_{l}} V {\bar{δ}}_{V} (r_{k} - x) \sum_{α, η \neq α}^{N_{a}} \frac{\partial Φ_{k α}}{\partial r_{k η}} \cdot {\tilde{v}}_{k α} \int_{r_{l η}}^{r_{k α}} {\bar{δ}}_{A}^{i} (φ - y) d φ

(22)

and

q_{}^{kin} (x, y, e^{i}) = - \sum_{k = 1}^{N_{l}} \sum_{α = 1}^{N_{α}} (\frac{1}{2} m_{α} {\tilde{v}}_{k α}^{2} + Φ_{k α}) {\tilde{v}}_{k} {\bar{δ}}_{A}^{i} (x - r_{k}) - \sum_{k = 1}^{N_{l}} V {\bar{δ}}_{V} (r_{k} - x) \sum_{α}^{N_{a}} (\frac{1}{2} m_{α} {\tilde{v}}_{k α}^{2} + Φ_{k α}) Δ {\tilde{v}}_{k α} {\bar{δ}}_{A}^{i} (y - Δ r_{k α}) .

(23)

The link between atomistic and continuum descriptions at this point has fundamentally departed from the IK formalism. As distinct from classical continuum mechanics, this new formulation describes the structure and properties of a crystalline system as a continuous function in x at the lattice level, cf. Figs. 3 and 5, while discrete in α or y at the atomic level; fluxes (stress and heat flux) are composed of components resulting from the motion and deformation of lattice cells (denoted as t and q, respectively) and the rearrangement of atoms within the cells (denoted as τ and j, respectively). Most importantly, instead of linking molecular variables to classical continuum mechanics equations, the CAC formalism leads to a new representation of conservation laws and fluxes.

FIG. 5.

View large Download slide

Schematic of the mass density distribution in space for a crystal with two different atoms in the primitive unit cell: (a) using Eq. (3), [(b) and (c)] using Eq. (7), showing that the mass density is a continuous and homogeneous function at the lattice level, while it is discrete and inhomogeneous at the atomic scale.

D. Temperature

An important difference between atomistic and classical continuum theories is the definition of temperature; in the former, temperature is a derived quantity, while in the latter, temperature is a fundamental state variable. Specifically, the abiding definition of temperature in classical molecular simulations is the kinetic temperature defined by the kinetic theory in terms of the mean kinetic energy of the random motion of atoms, while that in quantum calculations is based on the mean energy of phonons.

The kinetic theory provides a microscopic description of temperature through the classical mechanics of particles. The equipartition theorem of the kinetic theory asserts that the average kinetic energy of a particle in an equilibrium gaseous system with temperature T is

\frac{1}{2} m v_{k α}^{2} = \frac{3}{2} k_{B} T

⁠. Away from thermal equilibrium, the particle velocity is usually replaced by

{\tilde{v}}_{k α}

⁠, the difference between particle velocity and the velocity field. Such a temperature definition is consistent with the physical picture of an ideal-gas thermometer in which the velocity is measured relative to the co-moving frame of the kinetic thermometer.¹⁸ The nonequilibrium temperature is expressed¹⁹ as

\frac{3}{2} k_{B} T^{N} (x) = \sum_{k = 1}^{N_{l}} \frac{1}{2} m_{α} {\tilde{v}}_{k α}^{2} V {\bar{δ}}_{V} (x - r_{k}) .

(24)

This temperature definition, however, has unambiguous meaning only in thermal equilibrium or an approximate local thermal equilibrium condition such as that in a steady state. To elucidate this point, we decompose the kinetic energy density into two parts,

\begin{aligned} E_{α}^{kin} (x, t) & = \frac{1}{2} \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} m_{ξ} {(v_{k ξ})}^{2} {\bar{δ}}_{V} (x - r_{k}) {\bar{δ}}_{V α} (y - Δ r_{k ξ}) \\ = \frac{1}{2} ρ_{α} {(v_{α})}^{2} + \frac{1}{2} \sum_{k = 1}^{N_{l}} m_{α} {({\tilde{v}}_{k α})}^{2} {\bar{δ}}_{V} (x - r_{k}) = k_{α 1} + k_{α 2}, \end{aligned}

(25)

where

k_{α 1}

is the kinetic energy density represented by the velocity field and

k_{α 2}

is the kinetic energy due to the velocity difference between particle velocity and velocity field, usually interpreted as the random motion of particles. Although the definition of total kinetic energy is unambiguous, the velocity field can only be uniquely defined when it is zero (for thermal equilibrium) or a constant (for steady state). This is because the velocity field depends on length and time scales, i.e., the averaging, of the linear momentum. Thus, regardless of whether it is in an atomistic or a multiscale simulation, temperature can only be well defined at thermal equilibrium or steady state. It is a concept describing the state of a system rather than the dynamics of an atom.

Differing from the classical concept, the quantum description of temperature is expressed in terms of phonons. There are more phonons at high temperature and fewer phonons at low temperature. The average phonon number is determined by the equilibrium temperature T through the Planck constant ħ and the Boltzmann constant

k_{B}

⁠. This relationship is called the Bose-Einstein distribution, i.e.,

n (κ, ν) = {[\exp (\frac{ℏ ω (κ, ν)}{k_{B} T}) - 1]}^{- 1},

(26)

where

n (κ, ν)

is the phonon number for branch ν and wave vector κ in thermal equilibrium. The total kinetic energy of an equilibrium system can then be expressed in terms of all the available phonon modes in the systems, i.e.,

K = \frac{1}{2} \sum_{κ}^{N_{l}} \sum_{ν}^{3 N_{a}} ℏ ω (k, ν) [\frac{1}{2} + n (k, ν)] .

(27)

Recall that the number of phonon modes is equal to the number of the particle degrees of freedom of the system and that the total kinetic energy per phonon mode in the classical limit is equal to k_BT/2. The difference between the quantum and classical temperatures can thus be quantified through equating the phonon and classical descriptions of the total kinetic energies,²⁰ i.e.,

\begin{aligned} K = & \frac{1}{2} \sum_{κ}^{N_{l}} \sum_{ν}^{3 N_{a}} ℏ ω (k, ν) [\frac{1}{2} + n (k, ν)] \\ = & \frac{1}{2} \sum_{κ}^{N_{l}} \sum_{ν}^{3 N_{a}} ℏ ω (k, ν) {\frac{1}{2} + {[\exp (\frac{ℏ ω (k, ν)}{k_{B} T_{q u a n t u m}}) - 1]}^{- 1}} \\ = & \frac{3}{2} N_{l} N_{a} k_{B} T_{c l a s s i c a l} . \end{aligned}

(28)

For a given quantum temperature, the corresponding classical temperature can be calculated using Eq. (28). Figure 6 presents such calculations for LiF and SrTiO₃. As can be seen from Fig. 6, the two temperatures are different at low temperature, but the difference becomes negligible as the temperature increases. At 0 K quantum temperature, LiF and SrTiO₃ have a zero-point energy equivalent to 237.1 K and 279.7 K classical temperature, respectively, and hence there are corresponding motions at 0 K, i.e., the zero-point vibrations.

FIG. 6.

Quantum temperature vs classical temperature (blue solid curves) for (a) LiF and (b) SrTiO3, showing that the two temperatures converge at high temperature, indicated as the solid blue curves approaching the dotted green lines.

View large Download slide

Quantum temperature vs classical temperature (blue solid curves) for (a) LiF and (b) SrTiO₃, showing that the two temperatures converge at high temperature, indicated as the solid blue curves approaching the dotted green lines.

IV. FINITE ELEMENT FORMULATION

As shown in Fig. 5, a local density such as the mass density in a polyatomic system is not a continuous function at the atomic scale. By decomposing the atomic position into the position of the lattice cell, x, and its internal position within the unit cell at x, we are able to express a density a_α(x,t) or a(x, y, t) at an equilibrium or steady state system as a homogeneous and continuous function in x. Consequently, the field equations of the conservation laws can be discretized in x and solved using continuum-based numerical methods such as the finite element (FE) method.

A. Temperature in FE formulation

As shown in Eq. , the decomposition of the total kinetic energy into

k_{α 1}

and

k_{α 2}

depends on the length scale of the velocity field. In particular, in a FE model

k_{α 1}

is represented by FE nodal velocities, in which the displacement field in an element is interpolated with FE shape functions

S_{ξ} (x)

⁠, i.e.,

{\hat{u}}_{α} (x, t) = S_{ξ} (x) U_{ξ α} (t),

(29)

where

U_{ξ α} (t)

is the displacement vector of αth atom embedded within the ξth node of the element; the contribution of the αth atom to the kinetic energy in an element of volume V_e is thus given by

\underset{V_{e}}{\iint \int} k_{α 1} d V = \underset{V_{e}}{\iint \int} \frac{1}{2} ρ_{α} {({\dot{u}}_{α})}^{2} d V = \underset{V_{e}}{\iint \int} \frac{1}{2} ρ_{α} {(S_{ξ} (x) {\dot{U}}_{ξ α})}^{2} d V .

(30)

The kinetic energy of an element can also be expressed in terms of phonons using the phonon dispersion relations of the FE model.²¹ However, a FE model employing the usual trilinear shape functions cuts off phonons whose wavelengths are smaller than the element size. The motion of the FEs thus only represents a subset of the phonons of the system. This means that the FE mesh, i.e., the FE displacement approximation, uniquely determines the decomposition of total kinetic energy into $k_{α 1}$ and $k_{α 2}$ ⁠, where $k_{α 1}$ is represented by FE nodal velocities, while $k_{α 2}$ is related to $k_{α 1}$ through the temperature-dependent phonon number $n (κ, ν)$ ⁠.

A special case is for systems at high-temperature thermal equilibrium, in which every phonon mode or every particle DOF has the same energy

k_{B} T / 2

⁠. The total kinetic energy of an 8-node element containing

n_{l} = V_{e} / V

unit cells is thus given by

\underset{V_{e}}{\iint \int} k_{α 1} d V + \underset{V_{e}}{\iint \int} k_{α 2} d V = n_{l} \times \frac{3}{2} k_{B} T,

(31)

with

\underset{V_{e}}{\iint \int} k_{α 1} (x, t) d V = \underset{V_{e}}{\iint \int} \frac{1}{2} ρ_{α} {[S_{ξ} (x) {\dot{U}}_{ξ α}]}^{2} d V = 8 (\frac{3}{2} k_{B} T)

(32)

and

\underset{V_{e}}{\iint \int} k_{α 2} (x, t) d V = \underset{V_{e}}{\iint \int} \frac{1}{2} \sum_{k = 1}^{N_{l}} m_{α} {({\tilde{v}}_{k α})}^{2} {\bar{δ}}_{V} (x - r_{k}) d V = (n_{l} - 8) (\frac{3}{2} k_{B} T) .

(33)

Equations (32) and (33) can be used to prescribe FE nodal velocities for a given initial equilibrium temperature T or to find $k_{α 2}$ from $k_{α 1}$ ⁠.

B. Finite element equations

Substituting the conservation equations of mass into the linear momentum equation, we obtain

ρ_{α} {\ddot{u}}_{α} = f_{α}^{int} (x) + \frac{1}{V} \underset{\partial V}{\int \int ◯} t_{α}^{k i n} d S + \frac{1}{V} \underset{\partial V_{α}}{\int \int ◯} τ_{α}^{k i n} d S_{α} .

(34)

Denoting the sum of the two surface integrals of kinetic stresses as

f_{α}^{T} (x, t)

⁠, noting that the averaged

f_{α}^{T} (x, t)

over a long-time duration can be linked to

k_{α 2}

and hence the local temperature, we have

\frac{1}{t} \int_{0}^{t} f_{α}^{T} (x, τ) d τ \equiv \frac{1}{t} \int_{0}^{t} \frac{1}{V} \underset{\partial V}{\int \int ◯} t_{α}^{k i n} d S + \frac{1}{t} \int_{0}^{t} \frac{1}{V} \underset{\partial V_{α}}{\int \int ◯} τ_{α}^{k i n} d S_{α} = λ \nabla_{x} T,

(35)

where

\nabla_{x} T

and λ can be determined from

\nabla_{x} (k_{α 1} + k_{α 2})

⁠, in which

\nabla_{x} k_{α 1}

can be calculated based on FE nodal velocities, and

\nabla_{x} k_{α 2}

can be determined from

\nabla_{x} k_{α 1}

if we assume local thermal equilibrium and employ the assumption used in the linearized phonon Boltzmann equation, i.e.,

\nabla_{x} n_{k} (x, t) = \nabla_{x} n_{k}^{0} (x, t)

⁠,²² where

n_{k}^{0}

and

n_{k}

are the equilibrium and nonequilibrium phonon mode distribution for wave vector k, respectively.

Equation can be solved using the finite element method, with the weighted residual of Eq. over an element being given by

\underset{V_{e}}{\iint \int} S_{η} (x) (ρ_{α} {\ddot{u}}_{α} - f_{α}^{int} - f_{α}^{e x t} - f_{α}^{T}) d V = 0.

(36)

This is the weak (Galerkin) form of Eq. (34). The integral can be numerically evaluated using Gaussian quadrature, and the equation can be solved for equilibrium or nonequilibrium processes.

For systems at thermal equilibrium, $\nabla_{x} T = 0$ ⁠, $f_{α}^{T}$ is a fluctuating function of time with a zero mean everywhere except at the system's boundary if thermal expansion is constrained. It can thus be modeled as a constant force applied at the boundary. When atomic-scale fluctuations are important to the simulated phenomena such as thermal activation of dislocation motion or phase transition, $f_{α}^{T}$ can be modeled as a periodic or randomly fluctuating force in time with the frequencies of the phonons that are cut off by the FE shape functions.

For nonequilibrium processes, $\nabla_{x} T$ is a constant within an element if the usual linear interpolation is used. $f_{α}^{T}$ can then be modeled as a body force with a constant mean in each element but with periodic or random fluctuation in t. It is noted that, to model it as a random force, $f_{α}^{T}$ needs to be partitioned into a fluctuating part and a damping part, with the latter to be determined by and to balance the fluctuation force, according to the fluctuation-dissipation theorem.²³

C. Fluxes in FEA solutions

Similar to that in MD,^13–15,24 fluxes across a surface element in a CAC model can be computed through calculations of the interaction forces and atomic motions that cross the surface element. Not only can the stress and heat flux at any given surface element be calculated, their spatial distributions can also be plotted as contours for the entire model.

With the usual trilinear shape functions to approximate the displacement field in FE, maximum fluxes occur at the nodes of the elements. Thus, only fluxes at the nodes and a few other locations on the boundaries of the elements are needed for the flux contours, such as the unit cells shown in Fig. 7. In regions of structural disorders or discontinuities, such as those near a dislocation core or crack tip, atomic-scale fluxes can be calculated using Eqs. (18) and (19) and plotted as a contour at atomic resolution.

FIG. 7.

View large Download slide

(a) 2D schematic of an element containing many unit cells, and each unit cell contains two atoms. The shaded unit cells are those that can be used to calculate fluxes for interpolation of the flux fields in the element so as to generate contours of fluxes within the element; (b) shear stress contour near a dislocation core.

V. NUMERICAL EXAMPLES

Currently, there are two public versions of CAC codes, the CAC package in LAMMPS and PyCAC. The general procedure for a CAC simulation can be summarized as follows:

Select an interatomic potential(s) for the system;
Discretize the system based on its morphological length scales, using atomic resolution for atomic-scale structural disorders such as grain boundaries (GBs), while using FEs for single crystal regions with the mesh size determined by, e.g., the smallest possible dislocation spacing;
Relax the system with initial temperature. One may need to systematically change the initial configurations to obtain the dislocation or GB structures that have the minimal energy;
Apply loading and boundary conditions;
Run simulation to solve for and to visualize atomic trajectories, which are either atomically resolved or interpolated through FE nodal displacements, for observations of dynamic processes or emerging phenomena, using software such as Ovito for atomic visualization or Tecplot for FE visualization; and
Calculate fluxes (stress or heat flux) or other transport properties.

The accuracy and efficiency of the CAC method have been tested through one-to-one comparisons with MD in space- and time-resolved simulations of phenomena such as crack initiation and branching,^26,27 phase transitions,²⁸ dislocation nucleation,^29–32 loop formations^33–36 interactions with interfaces and other defects,^37–43 as well as phonon-dislocation interaction^44,45 and phonon-GB interaction.⁴⁶ In general, for an element containing N by N by N unit cells, the ratio of the number of mathematical operations in CAC to that in MD scales as O(N⁻³).⁴² This means that CAC can be very efficient by using large elements, e.g., for polycrystalline materials with micrometer-sized grains or for crystals with low dislocation density, i.e., large dislocation spacing.

To illustrate the applicability of CAC, in this section, we present three sets of prior simulation results. These simulations were performed at low temperature without considering the fluctuations from phonons whose wavelengths are smaller than the element size.

A. Atomic-scale structures of interfaces

CAC has been employed to study the structure and energy of grain boundaries (GBs) and their dynamic interaction with cracks and dislocations in SrTiO₃ polycrystals.^42,43 In these studies, GBs are atomically resolved, while regions away from GBs are coarse-grained. Simulation results show that both the energy and atomic displacements of SrTiO₃ GBs from CAC agree well with those that from MD if the width of the atomically resolved region for the GBs is larger than 12 Å.⁴² The ghost forces reported at the atomic-continuum interface in many multiscale models are not present in the CAC models.⁴² The atomic-scale GB structures from CAC⁴² are shown in Fig. 8 to compare well with those from experimental results.⁴⁷

FIG. 8.

CAC results of (a) ∑5(310), (b) ∑13(510), and (c) ∑5(210) GBs in SrTiO3, and (d) STEM images of ∑5(210) GB;47 the unit structures around each GB plane are outlined with solid lines.42 Reproduced with permission from Yang and Chen, Proc. R. Soc. A Math. Phys. Eng. Sci. 471, 20140758 (2015). Copyright 2015 The Royal Society.

View large Download slide

CAC results of (a) ∑5(310), (b) ∑13(510), and (c) ∑5(210) GBs in SrTiO₃, and (d) STEM images of ∑5(210) GB;⁴⁷ the unit structures around each GB plane are outlined with solid lines.⁴² Reproduced with permission from Yang and Chen, Proc. R. Soc. A Math. Phys. Eng. Sci. 471, 20140758 (2015). Copyright 2015 The Royal Society.

Figure 9 presents CAC simulation results of the misfit dislocation networks in the phase interfaces in heteroepitaxial thin films. CAC is shown to have reproduced the experimentally observed dislocation networks in PbTe/PbSe (001) interfaces.⁴⁸ The dislocation spacing in the system as a function of the epilayer thickness obtained by CAC is also in good agreement with the experimental measurements. In addition, the critical thickness above which misfit dislocations are observed to form in the heteroepitaxial thin film compares very well with experimental observations.⁴⁹

FIG. 9.

View large Download slide

CAC simulation results of (a) the dislocation network in a PbTe/PbSe (001) interface and (b) a comparison between CAC simulation results and experimental measurements of the dislocation spacing as a function of PbTe epilayer thickness; the inset in (a) is a STM image of the dislocation network in PbTe/PbSe (001).⁴⁸ Reproduced with permission from Springholz and Wiesauer, Phys. Rev. Lett. 88, 015507 (2001). Copyright 2001 American Physical Society.

B. Fracture and dislocations

Since the only constitutive relation in CAC is the nonlocal force field, and CAC is of an integral form, continuity between FEs is not required. Consequently, nucleation and propagation of dislocations and cracks can be simulated in the FE region via sliding and separation between elements. Simulation results can be output in terms of elements or in terms of atomic trajectories. This makes CAC a convenient tool for the simulation and visualization of the dynamic processes of defect nucleation, interaction, and propagation.

Figure 10 presents CAC simulation results for the loading of a notched solid normal to the notch depth. The simulations have captured the dynamics processes of crack initiation, propagation, and branching. They also reveal that the interactions of the stress waves that emitted from the propagating crack with those that are reflected back by the boundaries of the specimen lead to a sharp increase of the crack speed, which in turn triggers the dynamic instability and hence the branching of cracks in micrometer-scale specimens.²⁶ In addition, we have measured the critical stress intensity factor and observed a crack speed near the speed of Rayleigh wave before crack branching.

FIG. 10.

View large Download slide

Time sequence of CAC simulations of a brittle crystal showing (a) stress waves emitting from a propagating crack and [(b)–(d)] crack branching as a result of the interaction between waves propagating from the crack tips and those reflected from the specimen boundaries.⁵⁰ Reproduced with permission from Deng et al., Int. J. Plast. 26, 1402–1414 (2010). Copyright 2010 Elsevier Ltd.

CAC has also been used to characterize the complex dynamics of fast moving dislocations, including the energy intensities and the wavelengths of phonons emitted from sonic dislocations, as well as the velocity-dependent stress fluctuations around the core of moving dislocations. Figure 11 presents snapshots comparing dislocation motion simulated using CAC with a uniform mesh and that using MD, with the number of DOFs of the CAC model being 1% of that of MD. The dislocation velocity is measured to be ∼2900 m/s in both CAC and MD simulations. Mach cones in a V-shaped pattern of the phonon wave-fronts are observed in the wake of the sonic dislocations. Analysis of simulation results based on a wavelet transform shows that the faster a dislocation is moving, the longer the emitted phonon wavelength.^29–32

FIG. 11.

View large Download slide

Snapshots of dislocation motion in a V-notched specimen; only the left part of the specimen in MD and the right part of the specimen in CAC are displayed. The atoms in MD and the elements in CAC are color coded in displaying stress component $σ_{y y}$ ⁠. The color, red or blue, respectively, indicates tensile or compressive stress. Reproduced with permission from Xiong et al., Acta Mater. 104, 143–155 (2016). Copyright 2015 Acta Materialia Inc.

Figure 12 presents CAC simulation results of the sequential slip transfer of a₀ ̸2 [110](111) dislocations across a Σ3 (111) coherent twin boundary (CTB) in Ni⁵⁶ under ostensibly quasistatic strain rate conditions. Simulation results shows that CAC can serve as an effective tool for the study of slip transfer reactions of sequences of incoming and outgoing long-range dislocation pileups of mixed dislocation characters. Five embedded atom method (EAM) potentials^51–55 were compared; results demonstrate the uncertainty in computed dislocation-interface reactions associated with the use of a variety of interatomic potentials and suggest that the applicability of dislocation/GB interaction criteria in the literature for slip transfer may have limitations.

FIG. 12.

View large Download slide

Snapshots of dislocation pileup with dominant leading screw character impinging against a CTB. Atoms are colored by adaptive common neighbor analysis: red are of hexagonal-close packed local structures, blue are BCC atoms, and all FCC atoms are deleted. In (a), incoming dislocations approach the CTB. In (b), the leading dislocation is constricted at the CTB, where two Shockley partial dislocations are recombined into a full dislocation. All five potentials produced the same results in (a) and (b). In (c), with Mishin-EAM⁵¹ and Voter-EAM potentials,⁵² the dislocation cross-slips into the outgoing twinned grain via redissociation into two partials. In (d) with Angelo-EAM,⁵³ Foiles-EAM,⁵⁴ and Zhou-EAM,⁵⁵ the redissociated dislocation is absorbed by the CTB, with two partials gliding on the twin plane in opposite directions.⁵⁶ Reproduced with permission from Xu et al., JOM 69, 814–821 (2017). Copyright 2017 The Minerals, Metals & Materials Society.

C. Phonon thermal transport

CAC solves for approximate atomic trajectories. Collective motion of atoms, i.e., phonons, and the nucleation and propagation of defects, as well as phonon scattering by defects and other phonons, naturally emerges in the simulation. It thus provides a unified treatment for phonon transport and defect dynamics, as both are represented by the same atomic motion, without the need to assume the nature of phonon transport or scattering, or the mechanism for defect nucleation, propagation, and interactions.

Figure 13 presents CAC simulation results of a heat pulse propagating across GBs.⁴⁶ A phonon representation of heat pulses composed of spatio-temporal Gaussian wave packets obeying the Bose-Einstein distribution⁵⁷ is used to mimic the coherent lattice excitation by ultrashort lasers. As can be seen from Fig. 13, the simulation has captured the following:

The phenomenon of phonon focusing “caustics,”⁵⁸ manifested as the concentration of energy flux in the sixfold symmetric focusing pattern shown in Fig. 13(a);
The Kapitza resistance, shown in Fig. 13(d) as the discontinuities in the kinetic energy at the GBs;
Simultaneous ballistic and diffusive phonon transport across GBs, and a co-existence of coherent-incoherent phonon scattering by the GBs. These are clearly shown in Figs. 13(b)–13(d) as both wavelike transport of the propagating heat pulse and the breaking of the 6-fold symmetry of the phonon focusing pattern due to diffusive transport and incoherent scattering by the GBs; and
Local reconstruction of the GB structures, as shown in Fig. 13(e), as a result of the phonon-GB interaction during the heat pulse propagation.

FIG. 13.

View large Download slide

Snapshots of kinetic energy distribution in space during the propagation a heat pulse across GBs. GB locations are indicated by white arrows. Reproduced with permission from Chen et al., Acta Mater. 136, 355–365 (2017). Copyright 2017 Acta Materialia Inc.

VI. SUMMARY AND DISCUSSIONS

We have presented the CAC method for simulation of transport processes in crystalline materials. The formalism extends the Irving-Kirkwood statistical mechanical theory of transport processes for single-phase, single-component molecular systems to a concurrently coupled atomic-continuum description of polyatomic crystals. Conservation equations are derived in terms of instantaneous expressions of cell-averaged quantities using the mathematical theory of distributions.⁵⁹ Fluxes are then obtained that quantify the flow of conserved quantities across the lattice cells as well as those that flow back and forth within the cells, as direct consequences of the local density definition and Newton's second law.

Basic features of the CAC formulation, in comparison with classical continuum mechanics and molecular dynamics, may be summarized as follows:

As in classical continuum mechanics, the CAC governing equations are field equations, and hence can be discretized and solved using the finite element method. However, unlike the situation with continuum mechanics, each element node in CAC represents a unit cell that further contains a group of atoms; the atomic interaction and atomic-level crystal structure are thus fully built into the formulation, with no additional constitutive equations other than an interatomic potential being needed.
By including the internal DOFs of each lattice cell associated with atomic movement in the formulation, CAC can resolve the atomic structure of crystalline materials and reproduce both acoustic and optical phonon modes. It also can measure fluxes at the atomic scale, such as stress near a dislocation core or heat flux across an interface or a defect.
Since CAC is derived in a bottom-up manner from the atomistic, it becomes identical to the underlying atomistic model when discretized at the lattice level. CAC thus has the full applicability of MD in predicting materials behavior without the need for a priori assumption for mechanisms. It can also go beyond MD to scale up for mesoscale simulations by utilizing the continuous structure and deformation in regions where structure disorders or other discontinuities are absent or where atomic resolution is not needed.
CAC is thus naturally a concurrent multiscale tool that resolves details to full atomic resolution at regions of interest, while coarse-graining to thousands of atoms per finite element elsewhere; it admits the propagation of dislocations, addresses stacking fault structures in the FE regions, and provides a unified treatment for material microstructures, defects, and phonon transport.

A comparison of representative coarse grained (CG) methods can be found in Ref. 60, and a detailed review of dynamic multiscale methods can be found in Ref. 61. CAC is currently the only concurrently coupled atomistic and continuum method whose continuum description is not based on the classical continuum mechanics equations but a reformulated field equation of conservation laws. As distinct from other CG and multiscale methods, the CAC formulation links and unifies atomistic and continuum description of transport equations. There are two unique capabilities of the CAC formulation:

It can serve as a multiscale theory for concurrent atomistic-continuum simulation of polyatomic crystalline materials under a single set of governing equations, thereby expanding the current atomic-interaction-based simulation capabilities from the nanoscale to the mesoscale crystalline materials.
It provides consistent formulas for calculations of continuum quantities such as stress and heat flux in atomistic or CG atomistic simulations of heterogeneous materials, noting that the majority of flux formulas in the literature or currently implemented in MD software are only applicable to homogenized, single-phase, single-component materials.

Irving and Kirkwood recognized that the classical hydrodynamics is a single-scale theory as a result of a single-scale materials description.³ Beyond the single-scale continuum mechanics or hydrodynamics, we have the generalized continuum mechanics (GCM) that advocate and employ a concurrent two-level structural description of materials.^62–64 Well-established GCM theories include the Micromorphic theory,⁶⁵ the Micropolar theory,⁶⁶ and the Microstructural theory.⁶⁷ Unlike CAC, these GCMs were formulated via a top down approach.⁶⁸ Although the link between molecular dynamics and GCM has been attempted,^69–72 the popular GCMs are two-level continuum-continuum formulations,⁶⁰ lacking a description for the discrete structure or motion at the atomic scale.

We have also presented three sets of simulation results in this tutorial to illustrate the applicability of CAC in reproducing atomic-scale structures of interfaces, dynamic processes of dislocations and fracture, and phonon transport across GBs. These are low-temperature simulations, in which thermal fluctuations with wavelength smaller than the size of the elements in the FE regions are ignored.

For simulation of defect dynamics at finite temperature, the FE mesh size in CAC is determined by the density of defects as well as by the dominant wavelengths of phonons that interact with the defects. In regions where atomic-scale accuracy is necessary to capture interface and defect reconstruction, full atomic resolution is needed.

For low-temperature transport processes in which phonon-phonon scattering has a negligible effect, the usual FE trilinear shape functions may be replaced by those that combine the trilinear shape functions and Bloch wave functions consisting of phonons whose wavelengths are cut off by the shape function. This allows short-wavelength phonons from an atomically resolved region to traverse a coarse-grained FE region and enter another atomically resolved region⁷³ so as to model phonon scattering by dislocations, impurities, GB or phase interfaces with all possible phonons.

For finite-temperature phonon transport when high-frequency short-wavelength phonons are dominated by thermally resistive⁷⁴ or diffusive processes, the surface integrals of the kinetic momentum fluxes, i.e., $f_{α}^{T}$ ⁠, in the momentum equation can be modeled as a body force with a constant mean but fluctuating in time in terms of the frequencies of the missing phonons. This is intended to reproduce the phonon density of state of the underlying atomistic model, thereby approximating the effect of the thermal fluctuations by the missing phonons. Phonon transport represented by the motion and deformation of the elements, whether it is diffusive, ballistic, or simultaneously ballistic and diffusive, will naturally emerge in the simulation as the consequence of the conservation laws, the interatomic potential, and the initial temperature.

Finally, we would like to note that, while the conservation equations and expressions of fluxes and temperature in the CAC formulation are identical to those in an atomistic model of particles, the FE solutions to the conservation equations are approximate numerical solutions. Consequently, simulated phenomena and calculated fluxes using the FE approximations feature a trade-off of efficiency and accuracy. There are thus coarse-graining errors that may increase with the FE mesh size, depending on the nature of the phenomena and the morphological length scales of the materials. Advanced algorithms that can improve the accuracy and efficiency of the simulation, such as adaptive mesh refinement and advanced multiple time step algorithm, are needed. Also, the current CAC formulation is based on the traditional interatomic potentials with well-defined site energies. Future work will be needed to extend the formulation to quantum-accurate, machine-learned potentials.

ACKNOWLEDGMENTS

This material is based upon research supported by the U.S. Department of Energy (DOE), Basic Energy Sciences under Award No. DE-SC0006539. The contributions of D.L.M. and Y.C. are partially supported by the National Science Foundation (NSF) under Grant Nos. CMMI 1761553 and 1761512, respectively. We thank Zexi Zheng for generating the data plotted in Fig. 6.

APPENDIX A: CONSERVATION EQUATIONS FOR MONATOMIC CRYSTALS

The classical formulations of boundary value problems are based on the assumption that their solutions are smooth (i.e., with continuous derivatives) and that the equation is satisfied at every point in the domain of the problems. For problems that involve point mass, point charge, or other discreteness or discontinuities, the smoothness assumption breaks down. The mathematical tool for such problems is the theory of distributions, also known as the theory of generalized functions, developed by Laurent Schwartz in the 1950s.¹⁷ One generalized function frequently used in science and engineering is the Dirac delta-function introduced by Paul Dirac in his research in quantum mechanics and is employed in most statistical mechanics formulations to link a particle description to a field representation. The theory of distributions enables a rigorous solution to the mathematical difficulties in linking discrete and continuum descriptions that have been bypassed by ad hoc constructions or by heuristic arguments.

Appendixes A and B present formal derivations of conservation laws based on the mathematical theory of distributions. A more detailed discussion of the integral form of conservation laws and atomistic formulas for fluxes for monatomic crystals can be found in Refs. 13–16. Due to additional quantities such as temperature discussed in this work, there are notational changes in this paper.

1. Local density definition and the balance equation of linear momentum for monatomic crystals

We define the volumetric density of mass, linear momentum, and energy over a unit cell and a time-interval as

\begin{matrix} ρ (x, t) \equiv \sum_{k} m_{k} {\bar{δ}}_{V} (x - r_{k}), \\ ρ (x, t) v (x, t) \equiv \sum_{k} m_{k} v_{k} {\bar{δ}}_{V} (x - r_{k}), \\ ρ (x, t) e (x, t) \equiv \sum_{k} (\frac{1}{2} m_{k} v_{k}^{2} + Φ_{k}) {\bar{δ}}_{V} (x - r_{k}), \end{matrix}

(A1)

where

{\bar{δ}}_{V}

is the time step-averaged box function

δ_{V}

defined in and is used to serve the role of averaging the IK point-function local density over a unit cell V and a time-interval Δt.

The derivative of

δ_{V}

does not exist everywhere as a classical function, but it exists as a distribution. A distribution f(x) is not required to have a value at any point x, it is defined by its action on a test function φ(x) that infinitely many times differentiable and has a bounded support. For example, the Dirac delta-function δ(x) is defined as (δ, φ) = φ(0). In this sense, the derivative of δ_V is the surface distribution that acts on a test function φ as

(\frac{\partial}{\partial t} δ_{V} (x - r_{k}, φ (x, t)) = {\int \int ◯}_{\partial V} φ v_{k} \cdot n d S

⁠. Any smooth function g(x) can be viewed as a distribution whose action on a test function φ(x) is defined as the integral of the product g(x) φ(x) over the whole range of the variable x. Furthermore, any distribution f(x) is proved to be the limit of a sequence of smooth functions f_ɛ in the sense that the numerical sequence (f_ɛ, φ) converges to the number (f, φ) for any test function as the regularization parameter ɛ approaches 0. In this sense, the derivative of δ_V can formally be written as a definite integral of a δ-function, that is, as the limit of the corresponding integral in which the δ-function is replaced by a regularized distribution δ_ɛ that is a smooth function and converges to δ as

ε \to 0

in the distributional sense,⁵⁹ i.e.,

\frac{\partial}{\partial t} δ_{V} (x - r_{k}) = - lim_{ε \to 0} \frac{1}{V} \iint_{\partial V} v_{k} \cdot n δ_{ε} (x + x^{'} - r_{k}) d S ≜ - \frac{1}{V} \iint_{\partial V} v_{k} \cdot n δ (x + x^{'} - r_{k}) d S .

(A2)

For the rest of the appendixes, any formal integral of a δ-function, like the one in the right-hand side of [Eq. (A2)], is understood as a distribution defined by such a limit. Any relation between distributions obtained in this way is proved to be independent of the choice of regularization.⁵⁹

The internal force density is also a volume density. An important fact is that the total internal force in a volume element V is equal to the total interaction force transmitted across the enclosing surface ∂V, which can be described as a line-plane intersection problem. We denote an oriented surface element that is centered at x with normal n and area A as Aⁿ; the intersection of the line segment L_kl and surface element Aⁿ can be expressed as

\int_{L_{k l}} δ_{A}^{n} (x - φ) d φ \equiv \frac{1}{A} \iint_{A^{n}} \int_{L_{k l}} δ (x + x^{'} - φ) d φ n \cdot d S = \frac{1}{A} {\begin{matrix} 1 & if L_{k l} intersects A^{n} and n \cdot L_{k l} \geq 0, \\ - 1 & if L_{k l} intersects A^{n} and n \cdot L_{k l} < 0, \\ 0 & if L_{k l} does not intersect A^{n} . \end{matrix}

(A3)

It is noted that the time-interval average of $\int_{L_{k l}} δ_{A}^{n} (x - φ) d φ$ over a time step is itself.

The crossing of an atom through a surface element during a time step Δt can also be expressed as the intersection of the path of the atom

r_{k}

and a surface element Aⁿ with area A and normal n as

{\bar{δ}}_{A}^{n} (x - r_{k}) {\tilde{v}}_{k} \equiv \frac{1}{A Δ t} \iint_{A^{n}} \int_{0}^{Δ t} δ (x - r_{k}) {\tilde{v}}_{k} n \cdot d τ d S = \frac{1}{A Δ t} {\begin{matrix} 1 & if r_{k} intersects A^{n} in Δ t and n \cdot {\tilde{v}}_{k} \geq 0, \\ - 1 & if r_{k} intersects A^{n} in Δ t and n \cdot {\tilde{v}}_{k} < 0, \\ 0 & if r_{k l} does not intersect A^{n} . \end{matrix}

(A4)

Expressing the total internal force in the volume element V as the total interaction forces crossing the bounding surface ∂V, we obtain the internal force density in terms of site energy

Φ_{k}

⁠, according to Eq. , as

\begin{aligned} \sum_{k} F_{k} δ_{V} (x - r_{k}) = & \sum_{k, l} \frac{\partial Φ_{k}}{\partial r_{l}} δ_{V} (x - r_{k}) = \sum_{k, l \neq k} (\frac{\partial Φ_{k}}{\partial r_{l}} - \frac{\partial Φ_{l}}{\partial r_{k}}) δ_{V} (x - r_{k}) = \frac{1}{V} \underset{\partial V}{\int \int ◯} \sum_{k \in V, l \notin V} (\frac{\partial Φ_{k}}{\partial r_{l}} - \frac{\partial Φ_{l}}{\partial r_{k}}) \int_{L_{k l}} d φ δ (x + x^{'} - φ) d S \\ = & \frac{1}{V} \underset{\partial V}{\int \int ◯} \sum_{k, l (\neq k)} \frac{\partial Φ_{k}}{\partial r_{l}} \int_{L_{k l}} d φ δ (x + x^{'} - φ) d S . \end{aligned}

(A5)

The time rate of change of the linear momentum within V can be derived, as a result of Newton's second law applying to the particles in V, using Eqs. and , as

\begin{aligned} \frac{\partial}{\partial t} \sum_{k} m_{k} v_{k} {\bar{δ}}_{V} (x - r_{k}) = & \sum_{k} m_{k} {\dot{v}}_{k} {\bar{δ}}_{V} (x - r_{k}) + \sum_{k} m_{k} v_{k} \frac{\partial}{\partial t} {\bar{δ}}_{V} (x - r_{k}) = \sum_{k} F_{k} {\bar{δ}}_{V} (x - r_{k}) - \frac{1}{Δ t} \int_{0}^{Δ t} d τ \frac{1}{V} \underset{\partial V}{\int \int ◯} \sum_{k} m_{k} v_{k} v_{k} δ (x + x^{'} - r_{k}) \cdot n d S \\ = & \frac{1}{Δ t} \int_{0}^{Δ t} d τ \frac{1}{V} \underset{\partial V}{\int \int ◯} \sum_{k, l (\neq k)} \frac{\partial Φ_{k}}{\partial r_{l}} \int_{L_{k l}} d φ δ (x + x^{'} - φ) d S - \frac{1}{Δ t} \int_{0}^{Δ t} d τ \frac{1}{V} \underset{\partial V}{\int \int ◯} {\sum_{k} m_{k} {\tilde{v}}_{k} {\tilde{v}}_{k} δ (x + x^{'} - r_{k}) + ρ v v} \cdot n d S \\ = & \underset{\partial V}{\int \int ◯} t d^{2} x^{'} - \underset{\partial V}{\int \int ◯} ρ v v \cdot n d S, \end{aligned}

(A6)

where the stress vector t on the surface element of normal n is expressed as a line-plane intersection problem, cf. Eqs. and , as

\begin{aligned} t (x, n) & = \sum_{k, l (\neq k)} \frac{\partial Φ_{k}}{\partial r_{l}} \int_{L_{k l}} {\bar{δ}}_{A}^{n} (φ - x) d φ - \sum_{k} m_{k} {\tilde{v}}_{k} {\tilde{v}}_{k} {\bar{δ}}_{A}^{n} (x - r_{k}) \\ = \sum_{k < l} F_{k l} \int_{L_{k l}} {\bar{δ}}_{A}^{n} (φ - x) d φ - \sum_{k} m_{k} {\tilde{v}}_{k} {\tilde{v}}_{k} {\bar{δ}}_{A}^{n} (x - r_{k}), \end{aligned}

(A7)

where

F_{k l} ≜ - (\partial Φ_{l} / \partial r_{k} - \partial Φ_{k} / \partial r_{l})

is a measure of the strength of the interaction between atom k and atom l.

2. Conservation equation of energy for monatomic crystals

Similarly, the time rate of change of the energy density can also be expressed in terms of atomic variables as

\frac{\partial (ρ e)}{\partial t} = \frac{\partial}{\partial t} \sum_{k} E_{k} {\bar{δ}}_{V} (x - r_{k}) = \sum_{k} {\dot{E}}_{k} {\bar{δ}}_{V} (x - r_{k}) + \sum_{k} E_{k} \frac{\partial}{\partial t} {\bar{δ}}_{V} (x - r_{k}) .

(A8)

The first term in the right-hand side of Eq. can be expressed as

\begin{aligned} \sum_{k} {\dot{E}}_{k} {\bar{δ}}_{V} (x - r_{k}) = & \sum_{k} \frac{d}{d t} (\frac{1}{2} m_{k} v {_{k}^{2}}_{} + Φ_{k}) {\bar{δ}}_{V} (x - r_{k}) = \sum_{k} (m_{k} v_{k} \cdot {\dot{v}}_{k} + {\dot{Φ}}_{k}) {\bar{δ}}_{V} (x - r_{k}) = \sum_{k} (F_{k} \cdot v_{k} + \sum_{l} \frac{\partial Φ_{k}}{\partial r_{l}} \cdot {\dot{r}}_{l}) {\bar{δ}}_{V} (x - r_{k}) \\ = & \sum_{k, l} (- \frac{\partial Φ_{l}}{\partial r_{k}} \cdot v_{k} + \frac{\partial Φ_{k}}{\partial r_{l}} \cdot v_{l}) {\bar{δ}}_{V} (x - r_{k}) = \frac{1}{Δ t} \int_{0}^{Δ t} d τ \frac{1}{V} \underset{\partial V}{\int \int ◯} \sum_{k, l (\neq k)} \frac{\partial Φ_{k}}{\partial r_{l}} \cdot v_{l} \int_{L_{k l}} d φ δ (x + x^{'} - φ) d S, \end{aligned}

(A9)

and the second term, according to Eq. , as

\begin{aligned} \sum_{k} E_{k} \frac{\partial}{\partial t} {\bar{δ}}_{V} (x - r_{k}) & = - \frac{1}{Δ t} \int_{0}^{Δ t} d τ \underset{\partial V}{\int \int ◯} \sum_{k} \frac{1}{V} E_{k} v_{k} δ (x - x^{'} - r_{k}) \cdot n d S \\ = - \frac{1}{Δ t} \int_{0}^{Δ t} d τ \underset{\partial V}{\int \int ◯} \sum_{k} \frac{1}{V} (\frac{1}{2} m_{k} {\tilde{v}}_{k}^{2} + m_{k} {\tilde{v}}_{k} \cdot v + \frac{1}{2} m_{k} v^{2} + Φ_{k}) {\tilde{v}}_{k} δ (x + x^{'} - r_{k}) \cdot n d S - \frac{1}{V} \underset{\partial V}{\int \int ◯} ρ e v \cdot n d S . \end{aligned}

(A10)

Substituting Eqs. and into Eq. , we obtain the integral form of the energy conservation law as

\frac{\partial (ρ e)}{\partial t} = \underset{\partial V}{\int \int ◯} (q_{pot} + σ_{pot} \cdot v) \cdot n d S + \underset{\partial V}{\int \int ◯} (q_{kin} + σ_{kin} \cdot v - ρ e v) \cdot n d S

(A11)

and the heat flux vector in the indicial notation as

q^{i} (x, t) = - \sum_{k} (\frac{1}{2} m_{k} {\tilde{v}}_{k}^{2} + Φ_{k}) {\tilde{v}}_{k} {\bar{δ}}_{A}^{i} (x - r_{k}) + \sum_{k, l (\neq k)} \frac{\partial Φ_{k}}{\partial r_{l}} \cdot {\tilde{v}}_{l} \int_{L_{k l}} {\bar{δ}}_{A}^{i} (φ - x) d φ .

(A12)

APPENDIX B: CONSERVATION EQUATIONS FOR POLYATOMIC CRYSTALS

1. Balance equation of linear momentum for polyatomic crystals

The rate of change of the density of linear momentum defined in Eq. can be expressed as

\begin{aligned} \frac{\partial (ρ_{α} v_{α})}{\partial t} & = \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} m_{ξ} {\dot{v}}_{k ξ} {\bar{δ}}_{V} (r_{k} - x) {\bar{δ}}_{V α} (Δ r_{k ξ} - y) \\ + \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} m_{ξ} v_{k ξ} \frac{\partial}{\partial t} {\bar{δ}}_{V} (r_{k} - x) {\bar{δ}}_{V α} (Δ r_{k ξ} - y) \\ + \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} m_{ξ} v_{k ξ} {\bar{δ}}_{V} (r_{k} - x) \frac{\partial}{\partial t} {\bar{δ}}_{V α} (Δ r_{k ξ} - y) . \end{aligned}

(B1)

The first term in the right-hand side of Eq. can be derived according to Newton's second law and, in the absence of a body force, it is equal to the internal force density

f_{α}^{int} (x)

⁠. Note that the translational symmetry of potential energies requires that

\frac{\partial Φ_{k ξ}}{\partial r_{k ξ}} = - \sum_{l (\neq k)}^{N_{l}} \sum_{η}^{N_{a}} \frac{\partial Φ_{k ξ}}{\partial r_{l η}} - \sum_{η (\neq ξ)}^{N_{a}} \frac{\partial Φ_{k ξ}}{\partial r_{k η}} .

(B2)

This leads to the expression of the total force acting on atom kξ due to all the other atoms as

\begin{aligned} \frac{\partial Φ}{\partial r_{k ξ}} = & \sum_{l}^{N_{l}} \sum_{η}^{N_{a}} \frac{\partial Φ_{l η}}{\partial r_{k ξ}} = \sum_{l (\neq k)}^{N_{l}} \sum_{η}^{N_{a}} \frac{\partial Φ_{l η}}{\partial r_{k ξ}} + \sum_{η \neq ξ}^{N_{a}} \frac{\partial Φ_{k η}}{\partial r_{k ξ}} + \frac{\partial Φ_{k ξ}}{\partial r_{k ξ}} \\ = & \sum_{l (\neq k)}^{N_{l}} \sum_{η}^{N_{a}} \frac{\partial Φ_{l η}}{\partial r_{k ξ}} + \sum_{η \neq ξ}^{N_{a}} \frac{\partial Φ_{k η}}{\partial r_{k ξ}} - \sum_{l (\neq k)}^{N_{l}} \sum_{η}^{N_{a}} \frac{\partial Φ_{k ξ}}{\partial r_{l η}} - \sum_{η (\neq ξ)}^{N_{a}} \frac{\partial Φ_{k ξ}}{\partial r_{k η}} \\ = & \sum_{l (\neq k)}^{N_{l}} \sum_{η}^{N_{a}} (\frac{\partial Φ_{l η}}{\partial r_{k ξ}} - \frac{\partial Φ_{k ξ}}{\partial r_{l η}}) + \sum_{η \neq ξ}^{N_{a}} (\frac{\partial Φ_{k η}}{\partial r_{k ξ}} - \frac{\partial Φ_{k ξ}}{\partial r_{k η}}) . \end{aligned}

(B3)

Figure 14 shows that the total interaction force on atom α in the unit cell at point x can be decomposed into the interaction forces on atom α due to the interactions between atoms in different unit cells and that between atoms in the same unit cell. This then leads to the definitions of surface tractions,

t_{α}^{p o t}

and

τ_{α}^{p o t}

⁠. The internal force density

f_{α}^{int} (x)

can then be written, according to Eq. , using the line-plane intersection defined in Eq. , as

\begin{aligned} f_{α}^{int} (x) & = \frac{1}{V} \underset{\partial V}{\int \int ◯} \sum_{k \in V, l \notin V}^{N_{l}} \sum_{η}^{N_{a}} (\frac{\partial Φ_{l η}}{\partial r_{k α}} - \frac{\partial Φ_{k α}}{\partial r_{l η}}) \int_{r_{l η}}^{r_{k α}} δ (x + x^{'} - φ) d φ d S \\ + \underset{\partial V_{α}}{\int \int ◯} \sum_{k}^{N_{l}} {\bar{δ}}_{V} (x - r_{k}) \sum_{η \neq α}^{N_{a}} (\frac{\partial Φ_{k η}}{\partial r_{k α}} - \frac{\partial Φ_{k α}}{\partial r_{k η}}) \int_{r_{k η}}^{r_{k α}} δ (y + y^{'} - φ) d φ d S_{α} ≜ \frac{1}{V} \underset{\partial V}{\int \int ◯} t_{α}^{p o t} d S + \frac{1}{V} \underset{\partial V_{α}}{\int \int ◯} τ_{α}^{p o t} d S_{α}, \end{aligned}

(B4)

where

t_{α}^{p o t} (x, n) = \sum_{k \in V, l \notin V}^{N_{l}} \sum_{η}^{N_{a}} (\frac{\partial Φ_{l η}}{\partial r_{k α}} - \frac{\partial Φ_{k α}}{\partial r_{l η}}) \int_{r_{l η}}^{r_{k α}} {\bar{δ}}_{A}^{n} (x - φ) d φ,

(B5)

τ_{α}^{p o t} (x, y, n) = \sum_{k}^{N_{l}} V {\bar{δ}}_{V} (x - r_{k}) \sum_{η \neq α}^{N_{a}} (\frac{\partial Φ_{k η}}{\partial r_{k α}} - \frac{\partial Φ_{k α}}{\partial r_{k η}}) \int_{r_{k η}}^{r_{k α}} {\bar{δ}}_{A}^{n} (y - φ) d φ = \sum_{k \in V}^{N_{l}} \sum_{η \neq α}^{N_{a}} (\frac{\partial Φ_{k η}}{\partial r_{k α}} - \frac{\partial Φ_{k α}}{\partial r_{k η}}) \int_{r_{k η}}^{r_{k α}} {\bar{δ}}_{A}^{n} (y - φ) d φ .

(B6)

FIG. 14.

View large Download slide

(a) The volume element V of the unit cell at point x; (b) two contributions to $f_{α}^{int} (x)$ ⁠: (1) the interaction forces that cross ∂V and (2) the interaction forces that only cross $\partial V_{α}$ ⁠; (c) the volume element $V_{α}$ ⁠.

The second term in Eq. can be expressed according to Eq. , noting that all the cross terms involving the velocity differences vanish over a time-interval average and

\sum_{α = 1}^{N_{α}} δ_{V α} (y - Δ r_{k ξ}) = 1

⁠, as

\begin{aligned} \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} m_{ξ} v_{k ξ} \frac{\partial}{\partial t} {\bar{δ}}_{V} (x - r_{k}) {\bar{δ}}_{V α} (y - Δ r_{k ξ}) = & - \frac{1}{V} \underset{\partial V}{\int \int ◯} \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} m_{ξ} v_{k ξ} v_{k} \cdot n \bar{δ} (x + x^{'} - r_{k}) {\bar{δ}}_{V α} (y - Δ r_{k ξ}) d S \\ = & - \frac{1}{V} \underset{\partial V}{\int \int ◯} \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} m_{ξ} ({\tilde{v}}_{k} + Δ {\tilde{v}}_{k ξ} + v_{α}) ({\tilde{v}}_{k} + v) \cdot n \bar{δ} (x + x^{'} - r_{k}) {\bar{δ}}_{V α} (y - Δ r_{k ξ}) d S \\ = & - \frac{1}{V} \underset{\partial V}{\int \int ◯} \sum_{k = 1}^{N_{l}} m_{α} {\tilde{v}}_{k} {\tilde{v}}_{k} \cdot n \bar{δ} (x + x^{'} - r_{k}) d S - \frac{1}{V} \underset{\partial V}{\int \int ◯} ρ_{α} v_{α} v \cdot n d S ≜ \frac{1}{V} \underset{\partial V}{\int \int ◯} [t_{α}^{k i n} - ρ_{α} v_{α} v \cdot n] d S \end{aligned}

(B7)

and the third term as

\begin{aligned} \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} m_{ξ} v_{k ξ} {\bar{δ}}_{V} (x - r_{k}) \frac{\partial}{\partial t} {\bar{δ}}_{V α} (y - Δ r_{k ξ}) = & - \underset{\partial V_{α}}{\int \int ◯} \sum_{k = 1}^{N_{l}} δ_{V} (x - r_{k}) \sum_{ξ = 1}^{N_{a}} m_{ξ} v_{k ξ} Δ v_{k ξ} δ (y + y^{'} - Δ r_{k ξ}) \cdot n d S_{α} \\ = & - \sum_{k = 1}^{N_{l}} δ_{V} (x - r_{k}) \underset{\partial V_{α}}{\int \int ◯} \sum_{ξ = 1}^{N_{a}} m_{ξ} Δ {\tilde{v}}_{k ξ} Δ {\tilde{v}}_{k ξ} δ (y + y^{'} - Δ r_{k ξ}) \cdot n d S_{α} \\ + \frac{1}{V} \underset{\partial V_{α}}{\int \int ◯} ρ v_{α} Δ v_{α} \cdot n d S_{α} ≜ \frac{1}{V} \underset{\partial V_{α}}{\int \int ◯} [τ_{α}^{k i n} - ρ v_{α} Δ v_{α} \cdot n] d S_{α}, \end{aligned}

(B8)

with

t_{α}^{k i n} (x, n) = - \sum_{k = 1}^{N_{l}} m_{α} {\tilde{v}}_{k} {\tilde{v}}_{k} {\bar{δ}}_{A}^{n} (x - r_{k}),

(B9)

τ_{α}^{k i n} (x, y, n) = - \sum_{k}^{N_{l}} V {\bar{δ}}_{V} (x - r_{k}) m_{α} Δ {\tilde{v}}_{k α} Δ {\tilde{v}}_{k α} {\bar{δ}}_{A}^{n} (y - Δ r_{k α}) = - \sum_{k \in V}^{N_{l}} m_{α} Δ {\tilde{v}}_{k α} Δ {\tilde{v}}_{k α} {\bar{δ}}_{A}^{n} (y - Δ r_{k α}) .

(B10)

Substituting Eqs. , and into Eq. , we obtain the conservation equation of linear momentum as

\frac{\partial (ρ_{α} v_{α})}{\partial t} = \frac{1}{V} \underset{\partial V}{\int \int ◯} (t_{α} - ρ_{α} v_{α} v \cdot n) d S + \frac{1}{V} \underset{\partial V_{α}}{\int \int ◯} (τ_{α} - ρ_{α} v_{α} Δ v_{α} \cdot n) d S_{α},

(B11)

where

t_{α} = t_{α}^{p o t} (x, n) + t_{α}^{k i n} (x, n), and τ_{α} = τ_{α}^{p o t} (x, y, n) + τ_{α}^{k i n} (x, y, n) .

(B12)

The total stress vector on a surface element at point y in the unit cell located at x due to the interaction and motion of all the atoms is then expressed as a sum of potential and kinetic parts as

t (x, y, n) = \sum_{α = 1}^{N_{a}} (t_{α} + τ_{α}) = \sum_{α = 1}^{N_{a}} (t_{α}^{p o t} + t_{α}^{k i n} + τ_{α}^{p o t} + τ_{α}^{k i n}) = t^{p o t} (x, y, n) + t^{k i n} (x, y, n),

(B13)

where

\begin{aligned} t^{p o t} (x, y, n) = & \sum_{α = 1}^{N_{a}} (t_{α}^{p o t} + τ_{α}^{p o t}) = \sum_{k \in V, l \notin V}^{N_{l}} \sum_{α, η}^{N_{a}} (\frac{\partial Φ_{l η}}{\partial r_{k α}} - \frac{\partial Φ_{k α}}{\partial r_{l η}}) \int_{r_{l η}}^{r_{k α}} {\bar{δ}}_{A}^{n} (x - φ) d φ + \sum_{k \in V}^{N_{l}} \sum_{α, η \neq α}^{N_{a}} (\frac{\partial Φ_{k η}}{\partial r_{k α}} - \frac{\partial Φ_{k α}}{\partial r_{k η}}) \int_{r_{k η}}^{r_{k α}} {\bar{δ}}_{A}^{n} (y - φ) d φ \\ = & \sum_{k, l (\neq k)}^{N_{l}} \sum_{α, η}^{N_{a}} \int_{r_{l η}}^{r_{k α}} \frac{\partial Φ_{k α}}{\partial r_{l η}} {\bar{δ}}_{A}^{n} (x - φ) d φ + \sum_{k}^{N_{l}} \sum_{α, η \neq α}^{N_{a}} V {\bar{δ}}_{V} (x - r_{k}) \int_{r_{k η}}^{r_{k α}} \frac{\partial Φ_{k α}}{\partial r_{k η}} {\bar{δ}}_{A}^{n} (y - φ) d φ \end{aligned}

(B14)

and

t^{k i n} (x, y, n) = \sum_{α = 1}^{N_{a}} (t_{α}^{k i n} + τ_{α}^{k i n}) = - \sum_{α = 1}^{N_{a}} \sum_{k = 1}^{N_{l}} m_{α} {\tilde{v}}_{k} {\tilde{v}}_{k} {\bar{δ}}_{A}^{n} (x - r_{k}) - \sum_{α = 1}^{N_{a}} \sum_{k = 1}^{N_{l}} m_{α} Δ {\tilde{v}}_{k α} Δ {\tilde{v}}_{k α} V {\bar{δ}}_{V} (x - r_{k}) {\bar{δ}}_{A}^{n} (y - Δ r_{k α}) .

(B15)

The stress tensor can then be expressed in terms of the stress vector as

\begin{aligned} σ^{i j} (x, y) = t^{i} (x, y, e^{j}) or σ_{α}^{i j} (x, y) = t_{α}^{i} (x, y, e^{j}), i, j = 1, 2, 3. \end{aligned}

(B16)

2. Conservation equation of energy for polyatomic crystals

The rate of change of the energy density can be written as

\begin{aligned} \frac{\partial (ρ_{α} e_{α})}{\partial t} = & \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} {\dot{E}}_{k ξ} {\bar{δ}}_{V} (x - r_{k}) {\bar{δ}}_{V α} (y - Δ r_{k ξ}) + \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} E_{k ξ} \frac{\partial}{\partial t} {\bar{δ}}_{V} (x - r_{k}) {\bar{δ}}_{V α} (y - Δ r_{k ξ}) + \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} E_{k ξ} {\bar{δ}}_{V} (x - r_{k}) \frac{\partial}{\partial t} {\bar{δ}}_{V α} (y - Δ r_{k ξ}) . \end{aligned}

(B17)

With

\begin{aligned} {\dot{E}}_{k ξ} & = m_{k ξ} v_{k ξ} \cdot {\dot{v}}_{k ξ} + {\dot{Φ}}_{k ξ} = - \sum_{l}^{N_{l}} \sum_{η}^{N_{a}} (\frac{\partial Φ_{l η}}{\partial r_{k ξ}} \cdot v_{k ξ} - \frac{\partial Φ_{k ξ}}{\partial r_{l η}} \cdot v_{l η}) \\ = - \sum_{l (\neq k)}^{N_{l}} \sum_{η}^{N_{a}} \frac{\partial Φ_{l η}}{\partial r_{k ξ}} v_{k ξ} - \sum_{η \neq ξ}^{N_{a}} \frac{\partial Φ_{k η}}{\partial r_{k ξ}} v_{k ξ} - \frac{\partial Φ_{k ξ}}{\partial r_{k ξ}} v_{k ξ} + \sum_{l (\neq k)}^{N_{l}} \sum_{η}^{N_{a}} \frac{\partial Φ_{k ξ}}{\partial r_{l η}} \cdot v_{l η} + \sum_{η \neq ξ}^{N_{a}} \frac{\partial Φ_{k ξ}}{\partial r_{k η}} v_{k η} + \frac{\partial Φ_{k ξ}}{\partial r_{k ξ}} v_{k ξ} \\ = - \sum_{l (\neq k)}^{N_{l}} \sum_{η}^{N_{a}} (\frac{\partial Φ_{l η}}{\partial r_{k ξ}} v_{k ξ} - \frac{\partial Φ_{k ξ}}{\partial r_{l η}} \cdot v_{l η}) - \sum_{η \neq ξ}^{N_{a}} (\frac{\partial Φ_{k η}}{\partial r_{k ξ}} v_{k ξ} - \frac{\partial Φ_{k ξ}}{\partial r_{k η}} v_{k η}), \end{aligned}

(B18)

the first term in the right-hand side of Eq. can be expressed, using Eq. and the line-plane intersection theorem, as

\begin{aligned} \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} {\dot{E}}_{k ξ} {\bar{δ}}_{V} (x - r_{k}) {\bar{δ}}_{V α} (y - Δ r_{k ξ}) = & \frac{1}{V} \underset{\partial V}{\int \int ◯} \sum_{k \in V, l \notin V}^{N_{l}} \sum_{η}^{N_{a}} (\frac{\partial Φ_{l η}}{\partial r_{k α}} v_{k α} - \frac{\partial Φ_{k α}}{\partial r_{l η}} \cdot v_{l η}) \int_{r_{l η}}^{r_{k α}} δ (x + x^{'} - φ) d φ d S \\ + \sum_{k}^{N_{l}} {\bar{δ}}_{V} (x - r_{k}) \underset{\partial V_{α}}{\int \int ◯} \sum_{η \neq α}^{N_{a}} (\frac{\partial Φ_{k η}}{\partial r_{k α}} v_{k α} - \frac{\partial Φ_{k α}}{\partial r_{k η}} v_{k η}) \int_{r_{k η}}^{r_{k α}} δ (y + y^{'} - φ) d φ d S_{α} . \end{aligned}

(B19)

The second term can be derived according to the distributional derivative of the box function expressed in Eq. as

\begin{aligned} \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} E_{k ξ} \frac{\partial}{\partial t} {\bar{δ}}_{V} (x - r_{k}) {\bar{δ}}_{V α} (y - Δ r_{k ξ}) = & \frac{1}{V} \underset{\partial V}{\int \int ◯} \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} E_{k ξ} v_{k} \cdot n δ (x + x^{'} - r_{k}) {\bar{δ}}_{V α} (y - Δ r_{k ξ}) d S \\ = & - \frac{1}{V} \underset{\partial V}{\int \int ◯} \sum_{k = 1}^{N_{l}} (\frac{1}{2} m_{α} {\tilde{v}}_{k α}^{2} + {\tilde{v}}_{k α} v_{α} + {(v_{α})}^{2} + Φ_{k α}) {\tilde{v}}_{k} \cdot n δ (x + x^{'} - r_{k}) d S - \underset{\partial V}{\int \int ◯} ρ_{α} e_{α} v \cdot n d S \\ = & - \frac{1}{V} \underset{\partial V}{\int \int ◯} \sum_{k = 1}^{N_{l}} (\frac{1}{2} m_{α} {\tilde{v}}_{k α}^{2} + Φ_{k α}) {\tilde{v}}_{k} δ (x + x^{'} - r_{k}) \cdot n d S \\ - \frac{1}{V} \frac{1}{Δ t} \int_{0}^{Δ t} d t \underset{\partial V}{\int \int ◯} \sum_{k = 1}^{N_{l}} {\tilde{v}}_{k α} v_{α} {\tilde{v}}_{k} \bar{δ} (x + x^{'} - r_{k}) \cdot n d S - \underset{\partial V}{\int \int ◯} ρ_{α} e_{α} v \cdot n d S \\ ≜ & \frac{1}{V} \underset{\partial V}{\int \int ◯} (q_{α}^{k i n} + t_{α}^{k i n} \cdot v_{α} - ρ_{α} e_{α} v) \cdot n d S \end{aligned}

(B20)

and the third term as

\begin{aligned} \sum_{k = 1}^{N_{l}} \sum_{ξ = 1}^{N_{a}} E_{k ξ} {\bar{δ}}_{V} (x - r_{k}) \frac{\partial}{\partial t} {\bar{δ}}_{V α} (y - Δ r_{k ξ}) = & \sum_{k = 1}^{N_{l}} {\bar{δ}}_{V} (x - r_{k}) \underset{\partial V_{α}}{\int \int ◯} \sum_{ξ = 1}^{N_{a}} E_{k ξ} Δ v_{k ξ} \cdot n δ (y + y^{'} - Δ r_{k ξ}) d S_{α} \\ = & - \sum_{k = 1}^{N_{l}} {\bar{δ}}_{V} (x - r_{k}) \underset{\partial V_{α}}{\int \int ◯} \sum_{ξ = 1}^{N_{a}} (\frac{1}{2} m_{ξ} {\tilde{v}}_{k ξ}^{2} + Φ_{k ξ} + {\tilde{v}}_{k ξ} v_{α} + {(v_{α})}^{2}) Δ {\tilde{v}}_{k ξ} \cdot n δ (y - Δ r_{k ξ}) d S_{α} \\ - \underset{\partial V V_{α}}{\int \int ◯} ρ_{α} e_{α} Δ v_{α} \cdot n d S ≜ \frac{1}{V} \underset{\partial V_{α}}{\int \int ◯} (j_{α}^{k i n} + τ_{α}^{k i n} \cdot v_{α} - ρ_{α} e_{α} Δ v_{α}) \cdot n d S_{α} . \end{aligned}

(B21)

This leads to

\frac{\partial (ρ_{α} e_{α})}{\partial t} = \frac{1}{V} \underset{\partial V}{\int \int ◯} (q_{α} + t_{α} \cdot v_{α} - ρ_{α} e_{α} v) \cdot n d S + \frac{1}{V} \underset{\partial V_{α}}{\int \int ◯} (j_{α} + τ_{α} \cdot v_{α} - ρ_{α} e_{α} Δ v_{α}) \cdot n d S_{α} .

(B22)

The total heat flux is thus expressed as

q (x, y) = \sum_{α = 1}^{N_{a}} (q_{α} + j_{α}) = q^{i} e^{i},

(B23)

with

q^{i} (x, y) = \sum_{k, l \neq k}^{N_{l}} \sum_{α, η}^{N_{a}} \frac{\partial Φ_{k α}}{\partial r_{l η}} \cdot {\tilde{v}}_{k α} \int_{r_{l η}}^{r_{k α}} {\bar{δ}}_{A}^{i} (x - φ) d φ - \sum_{k = 1}^{N_{l}} \sum_{α = 1}^{N_{α}} (\frac{1}{2} m_{α} {\tilde{v}}_{k α}^{2} + Φ_{k α}) {\tilde{v}}_{k} {\bar{δ}}_{A}^{i} (x - r_{k})

(B24)

and

\begin{aligned} j^{i} (x, y) = & \sum_{k}^{N_{l}} V {\bar{δ}}_{V} (x - r_{k}) \sum_{α, η \neq α}^{N_{a}} \frac{\partial Φ_{k α}}{\partial r_{k η}} \cdot {\tilde{v}}_{k α} \int_{r_{k η}}^{r_{k α}} {\bar{δ}}_{A}^{i} (y - φ) d φ \\ - \sum_{k = 1}^{N_{l}} V {\bar{δ}}_{V} (r_{k} - x) \sum_{α}^{N_{a}} (\frac{1}{2} m_{α} {\tilde{v}}_{k α}^{2} + Φ_{k α}) \\ \times Δ {\tilde{v}}_{k α} {\bar{δ}}_{A}^{i} (y - Δ r_{k α}) . \end{aligned}

(B25)

Irving and Kirkwood noted that a definition for the site energies of particles is required for definition of energy density as well as for the derivation of a field representation of energy conservation law in terms of particle variables.³ The conservation equations and the expressions for fluxes derived in the appendixes are valid for any analytical many-body potentials with well-defined site energies. A detailed comparison and discussion of the fluxes defined as a surface density presented in Appendix A and the fluxes formulas expressed as a volume density can be found in Ref. 13. Specifically, analytical and MD simulation results presented in Ref. 13 show that replacing the Dirac δ with a finite-sized volume weighting function changes the fundamental nature of fluxes as a surface density. Consequently, the dynamic balance between the rate of change of the total momentum or energy in a volume element and the fluxes across the surface boundary cannot be established, leading to the failure of the volume averaged flux formulas in satisfying the momentum and energy conservation laws as well as typical transport boundary conditions. The equivalence between flux formulas expressed as a pair density derived using the kinetic theory^75,76 and that expressed as a single density derived using statistical mechanics,³ as ensemble averaged point functions or as surface-averages, can be found in Ref. 16.

REFERENCES

1.

C.

Kittel

,

Introduction to Solid State Physics

(

John Wiley & Sons, Inc.

,

New York

,

1956

).

Google Scholar

2.

M.

Born

and

K.

Huang

,

Dynamical Theory of Crystal Lattices

(

Clarendon Press

,

Oxford

,

1954

).

Google Scholar

3.

J.

Irving

and

J. G.

Kirkwood

, “

The statistical mechanical theory of transport processes. IV. The equations of hydrodynamics

,”

J. Chem. Phys.

18

,

817

–

829

(

1950

).

https://doi.org/10.1063/1.1747782

Google Scholar

Crossref

4.

P.

Dirac

,

Principles of Quantum Mechanics

(Clarendon Press, Oxford,

1930

).

Google Scholar

5.

B.

Todd

,

P. J.

Daivis

, and

D. J.

Evans

, “

Heat flux vector in highly inhomogeneous nonequilibrium fluids

,”

Phys. Rev. E

51

,

4362

(

1995

).

https://doi.org/10.1103/PhysRevE.51.4362

Google Scholar

Crossref

6.

B. D.

Todd

,

D. J.

Evans

, and

P. J.

Daivis

, “

Pressure tensor for inhomogeneous fluids

,”

Phys. Rev. E

52

,

1627

–

1638

(

1995

).

https://doi.org/10.1103/PhysRevE.52.1627

Google Scholar

Crossref

7.

Y.

Chen

, “

Local stress and heat flux in atomistic systems involving three-body forces

,”

J. Chem. Phys.

124

,

054113

(

2006

).

https://doi.org/10.1063/1.2166387

Google Scholar

Crossref

PubMed

8.

E.

Smith

,

D.

Heyes

,

D.

Dini

, and

T.

Zaki

, “

Control-volume representation of molecular dynamics

,”

Phys. Rev. E

85

,

056705

(

2012

).

https://doi.org/10.1103/PhysRevE.85.056705

Google Scholar

Crossref

9.

J. G.

Kirkwood

, “

The statistical mechanical theory of transport processes. I. General theory

,”

J. Chem. Phys.

14

,

180

–

201

(

1946

).

https://doi.org/10.1063/1.1724117

Google Scholar

Crossref

10.

Y.

Chen

, “

Reformulation of microscopic balance equations for multiscale materials modeling

,”

J. Chem. Phys.

130

,

134706

(

2009

).

https://doi.org/10.1063/1.3103887

Google Scholar

Crossref

PubMed

11.

Y.

Chen

and

J.

Lee

, “

Atomistic formulation of a multiscale field theory for nano/micro solids

,”

Philos. Mag.

85

,

4095

–

4126

(

2005

).

https://doi.org/10.1080/14786430500362595

Google Scholar

Crossref

12.

D. J.

Evans

and

G.

Morriss

,

Statistical Mechanics of Nonequilibrium Liquids

(

Cambridge University Press

,

Cambridge

,

2008

).

Google Scholar

Crossref

13.

Y.

Chen

and

A.

Diaz

, “

Physical foundation and consistent formulation of atomic-level fluxes in transport processes

,”

Phys. Rev. E

98

,

052113

(

2018

).

https://doi.org/10.1103/PhysRevE.98.052113

Google Scholar

Crossref

14.

Y.

Chen

and

A.

Diaz

, “

Local momentum and heat fluxes in transient transport processes and inhomogeneous systems

,”

Phys. Rev. E

94

,

053309

(

2016

).

https://doi.org/10.1103/PhysRevE.94.053309

Google Scholar

Crossref

PubMed

15.

Y.

Chen

, “

The origin of the distinction between microscopic formulas for stress and Cauchy stress

,”

Europhys. Lett.

116

,

34003

(

2016

).

https://doi.org/10.1209/0295-5075/116/34003

Google Scholar

Crossref

16.

A.

Diaz

,

D.

Davydov

, and

Y.

Chen

, “

On the equivalence of the two foundational formulations for atomistic flux in inhomogeneous transport processes

,”

Proc. R. Soc. A

475

,

20180688

(

2019

).

https://doi.org/10.1098/rspa.2018.0688

Google Scholar

Crossref

17.

L.

Schwartz

,

Théorie des Distributions

(

Hermann

,

Paris

,

1950

).

Google Scholar

18.

W. G.

Hoover

,

Computational Statistical Mechanics

(

Elsevier

,

New York

,

1991

).

Google Scholar

19.

S.

Root

,

R. J.

Hardy

, and

D. R.

Swanson

, “

Continuum predictions from molecular dynamics simulations: Shock waves

,”

J. Chem. Phys.

118

,

3161

–

3165

(

2003

).

https://doi.org/10.1063/1.1537243

Google Scholar

Crossref

20.

J.

Turney

,

A.

McGaughey

, and

C.

Amon

, “

Assessing the applicability of quantum corrections to classical thermal conductivity predictions

,”

Phys. Rev. B

79

,

224305

(

2009

).

https://doi.org/10.1103/PhysRevB.79.224305

Google Scholar

Crossref

21.

Y.

Li

et al, “

Phonon spectrum and phonon focusing in coarse-grained atomistic simulations

,”

Comput. Mater. Sci.

162

,

21

–

32

(

2019

).

https://doi.org/10.1016/j.commatsci.2019.02.020

Google Scholar

Crossref

22.

A.

Cantarero

and

F. X.

Àlvarez

,

Nanoscale Thermoelectrics

(

Springer

,

New York

,

2014

), pp.

1

–

39

.

Google Scholar

23.

R.

Kubo

, “

The fluctuation-dissipation theorem

,”

Rep. Prog. Phys.

29

,

255

(

1966

).

https://doi.org/10.1088/0034-4885/29/1/306

Google Scholar

Crossref

24.

J.

Rigelesaiyin

,

A.

Diaz

,

W.

Li

,

L.

Xiong

, and

Y.

Chen

, “

Asymmetry of the atomic-level stress tensor in homogeneous and inhomogeneous materials

,”

Proc. R. Soc. A Math. Phys. Eng. Sci.

474

,

20180155

(

2018

).

https://doi.org/10.1098/rspa.2018.0155

Google Scholar

Crossref

25.

S.

Xu

et al, “

PyCAC: The concurrent atomistic-continuum simulation environment

,”

J. Mater. Res.

33

,

857

–

871

(

2018

).

https://doi.org/10.1557/jmr.2018.8

Google Scholar

Crossref

26.

Q.

Deng

,

L.

Xiong

, and

Y.

Chen

, “

Coarse-graining atomistic dynamics of brittle fracture by finite element method

,”

Int. J. Plast.

26

,

1402

–

1414

(

2010

).

https://doi.org/10.1016/j.ijplas.2010.04.007

Google Scholar

Crossref

27.

Q.

Deng

and

Y.

Chen

, “

A coarse-grained atomistic method for 3D dynamic fracture simulation

,”

J. Multiscale Comput. Eng.

11

,

227

–

237

(

2013

).

https://doi.org/10.1615/IntJMultCompEng.2013005442

Google Scholar

Crossref

28.

L.

Xiong

and

Y.

Chen

, “

Coarse-grained simulations of single-crystal silicon

,”

Model. Simul. Mater. Sci. Eng.

17

,

035002

(

2009

).

https://doi.org/10.1088/0965-0393/17/3/035002

Google Scholar

Crossref

29.

L.

Xiong

,

Q.

Deng

,

G.

Tucker

,

D. L.

McDowell

, and

Y.

Chen

, “

A concurrent scheme for passing dislocations from atomistic to continuum domains

,”

Acta Mater.

60

,

899

–

913

(

2012

).

https://doi.org/10.1016/j.actamat.2011.11.002

Google Scholar

Crossref

30.

L.

Xiong

,

G.

Tucker

,

D. L.

McDowell

, and

Y.

Chen

, “

Coarse-grained atomistic simulation of dislocations

,”

J. Mech. Phys. Solids

59

,

160

–

177

(

2011

).

https://doi.org/10.1016/j.jmps.2010.11.005

Google Scholar

Crossref

31.

S.

Yang

,

L.

Xiong

,

Q.

Deng

, and

Y.

Chen

, “

Concurrent atomistic and continuum simulation of strontium titanate

,”

Acta Mater.

61

,

89

–

102

(

2013

).

https://doi.org/10.1016/j.actamat.2012.09.032

Google Scholar

Crossref

32.

L.

Xiong

et al, “

Coarse-grained elastodynamics of fast moving dislocations

,”

Acta Mater.

104

,

143

–

155

(

2016

).

https://doi.org/10.1016/j.actamat.2015.11.037

Google Scholar

Crossref

33.

L.

Xiong

,

Q.

Deng

,

G. J.

Tucker

,

D. L.

McDowell

, and

Y.

Chen

, “

Coarse-grained atomistic simulations of dislocations in Al, Ni and Cu crystals

,”

Int. J. Plast.

38

,

86

–

101

(

2012

).

https://doi.org/10.1016/j.ijplas.2012.05.002

Google Scholar

Crossref

34.

L.

Xiong

,

D. L.

McDowell

, and

Y.

Chen

, “

Nucleation and growth of dislocation loops in Cu, Al and Si by a concurrent atomistic-continuum method

,”

Scr. Mater.

67

,

633

–

636

(

2012

).

https://doi.org/10.1016/j.scriptamat.2012.07.026

Google Scholar

Crossref

35.

S.

Xu

,

L.

Xiong

,

Y.

Chen

, and

D. L.

McDowell

, “

Edge dislocations bowing out from a row of collinear obstacles in Al

,”

Scr. Mater.

123

,

135

–

139

(

2016

).

https://doi.org/10.1016/j.scriptamat.2016.06.018

Google Scholar

Crossref

36.

S.

Xu

,

L.

Xiong

,

Y.

Chen

, and

D. L.

McDowell

, “

An analysis of key characteristics of the Frank-Read source process in FCC metals

,”

J. Mech. Phys. Solids

96

,

460

–

476

(

2016

).

https://doi.org/10.1016/j.jmps.2016.08.002

Google Scholar

Crossref

37.

S.

Xu

,

L.

Xiong

,

Y.

Chen

, and

D. L.

McDowell

, “

Sequential slip transfer of mixed-character dislocations across Σ3 coherent twin boundary in FCC metals: A concurrent atomistic-continuum study

,”

npj Comput. Mater.

2

,

15016

(

2016

).

https://doi.org/10.1038/npjcompumats.2015.16

Google Scholar

Crossref

38.

L.

Xiong

,

S.

Xu

,

D. L.

McDowell

, and

Y.

Chen

, “

Concurrent atomistic–continuum simulations of dislocation–void interactions in fcc crystals

,”

Int. J. Plast.

65

,

33

–

42

(

2015

).

https://doi.org/10.1016/j.ijplas.2014.08.002

Google Scholar

Crossref

39.

S.

Xu

,

L.

Xiong

,

Y.

Chen

, and

D. L.

McDowell

, “

Shear stress- and line length-dependent screw dislocation cross-slip in FCC Ni

,”

Acta Mater.

122

,

412

–

419

(

2017

).

https://doi.org/10.1016/j.actamat.2016.10.005

Google Scholar

Crossref

40.

S.

Xu

,

L.

Xiong

,

Y.

Chen

, and

D.

McDowell

, “

Validation of the concurrent atomistic-continuum method on screw dislocation/stacking fault interactions

,”

Crystals

7

,

120

(

2017

).

https://doi.org/10.3390/cryst7050120

Google Scholar

Crossref

41.

S.

Yang

and

Y.

Chen

, “Concurrent atomistic-continuum simulation of defects in polyatomic ionic materials,” in Multiscale Materials Modeling for Nanomechanics, edited by C. Weinberger and G. Tucker (Springer, Cham, Switzerland) Vol. 245, pp. 261–296.

42.

S.

Yang

and

Y.

Chen

, “

Concurrent atomistic and continuum simulation of bi-crystal strontium titanate with tilt grain boundary

,”

Proc. R. Soc. A Math. Phys. Eng. Sci.

471

,

20140758

(

2015

).

https://doi.org/10.1098/rspa.2014.0758

Google Scholar

Crossref

43.

S.

Yang

,

N.

Zhang

, and

Y.

Chen

, “

Concurrent atomistic–continuum simulation of polycrystalline strontium titanate

,”

Philos. Mag.

95

,

2697

–

2716

(

2015

).

https://doi.org/10.1080/14786435.2015.1076178

Google Scholar

Crossref

44.

X.

Chen

,

L.

Xiong

,

D. L.

McDowell

, and

Y.

Chen

, “

Effects of phonons on mobility of dislocations and dislocation arrays

,”

Scr. Mater.

137

,

22

–

26

(

2017

).

https://doi.org/10.1016/j.scriptamat.2017.04.033

Google Scholar

Crossref

45.

L.

Xiong

,

D. L.

McDowell

, and

Y.

Chen

, “

Sub-THz phonon drag on dislocations by coarse-grained atomistic simulations

,”

Int. J. Plast.

55

,

268

–

278

(

2014

).

https://doi.org/10.1016/j.ijplas.2013.11.004

Google Scholar

Crossref

46.

X.

Chen

et al, “

Ballistic-diffusive phonon heat transport across grain boundaries

,”

Acta Mater.

136

,

355

–

365

(

2017

).

https://doi.org/10.1016/j.actamat.2017.06.054

Google Scholar

Crossref

47.

H.-S.

Lee

et al, “

Defect energetics in SrTiO₃ symmetric tilt grain boundaries

,”

Phys. Rev. B

83

,

104110

(

2011

).

https://doi.org/10.1103/PhysRevB.83.104110

Google Scholar

Crossref

48.

G.

Springholz

and

K.

Wiesauer

, “

Nanoscale dislocation patterning in PbTe/PbSe (001) lattice-mismatched heteroepitaxy

,”

Phys. Rev. Lett.

88

,

015507

(

2001

).

https://doi.org/10.1103/PhysRevLett.88.015507

Google Scholar

Crossref

49.

Y.

Li

,

Z.

Fan

,

W.

Li

,

D. L.

McDowell

, and

Y.

Chen

, “

A multiscale study of misfit dislocations in PbTe/PbSe(001) heteroepitaxy

,”

J. Mater. Res.

,

34

,

2306

–

2314

(

2019

).

https://doi.org/10.1557/jmr.2019.69

Google Scholar

Crossref

50.

Q.

Deng

, “Coarse-graining atomistic dynamics of fracture by finite element method formulation, parallelization and applications,” Ph.D. dissertation (

University of Florida

,

2011

).

51.

Y.

Mishin

,

D.

Farkas

,

M.

Mehl

, and

D.

Papaconstantopoulos

, “

Interatomic potentials for monoatomic metals from experimental data and ab initio calculations

,”

Phys. Rev. B

59

,

3393

(

1999

).

https://doi.org/10.1103/PhysRevB.59.3393

Google Scholar

Crossref

52.

A. F.

Voter

and

S. P.

Chen

, “

Accurate interatomic potentials for Ni, Al and Ni₃Al

.”

MRS Online Proc. Libr. Arch.

82,

175 (

1986

).

https://doi.org/10.1557/PROC-82-175

Google Scholar

53.

J. E.

Angelo

,

N. R.

Moody

, and

M. I.

Baskes

, “

Trapping of hydrogen to lattice defects in nickel

,”

Model. Simul. Mater. Sci. Eng.

3

,

289

(

1995

).

https://doi.org/10.1088/0965-0393/3/3/001

Google Scholar

Crossref

54.

S. M.

Foiles

and

J.

Hoyt

, “

Computation of grain boundary stiffness and mobility from boundary fluctuations

,”

Acta Mater.

54

,

3351

–

3357

(

2006

).

https://doi.org/10.1016/j.actamat.2006.03.037

Google Scholar

Crossref

55.

X.

Zhou

,

R.

Johnson

, and

H.

Wadley

, “

Misfit-energy-increasing dislocations in vapor-deposited CoFe/NiFe multilayers

,”

Phys. Rev. B

69

,

144113

(

2004

).

https://doi.org/10.1103/PhysRevB.69.144113

Google Scholar

Crossref

56.

S.

Xu

,

L.

Xiong

,

Y.

Chen

, and

D. L.

McDowell

, “

Comparing EAM potentials to model slip transfer of sequential mixed character dislocations across Two symmetric tilt grain boundaries in Ni

,”

JOM

69

,

814

–

821

(

2017

).

https://doi.org/10.1007/s11837-017-2302-1

Google Scholar

Crossref

57.

X.

Chen

,

A.

Chernatynskiy

,

L.

Xiong

, and

Y.

Chen

, “

A coherent phonon pulse model for transient phonon thermal transport

,”

Comput. Phys. Commun.

195

,

112

–

116

(

2015

).

https://doi.org/10.1016/j.cpc.2015.05.008

Google Scholar

Crossref

58.

J. P.

Wolfe

,

Imaging Phonons: Acoustic Wave Propagation in Solids

(

Cambridge University Press

,

Cambridge

,

1998

).

Google Scholar

Crossref

59.

V. S.

Vladimirov

,

A.

Jeffrey

, and

F. E.

Schroeck

, “

Equations of mathematical physics

,”

Am. J. Phys.

39

,

1548

–

1548

(

1971

).

https://doi.org/10.1119/1.1976723

Google Scholar

Crossref

60.

Y.

Chen

,

J.

Zimmerman

,

A.

Krivtsov

, and

D.

McDowell

, “

Assessment of atomistic coarse-graining methods

,”

Int. J. Eng. Sci.

49

,

1337

–

1349

(

2011

).

https://doi.org/10.1016/j.ijengsci.2011.03.018

Google Scholar

Crossref

61.

A.

Diaz

,

D.

McDowell

, and

Y.

Chen

,

Generalized Models and Non-Classical Approaches in Complex Materials 2

(

Springer

,

New York

,

2018

), pp.

55

–

77

.

Google Scholar

Crossref

62.

A. C.

Eringen

,

Microcontinuum Field Theories: Foundations and Solids

(

Springer

,

New York

,

1999

), Vol. 487.

Google Scholar

Crossref

63.

Y.

Chen

,

J. D.

Lee

, and

A.

Eskandarian

, “

Atomistic viewpoint of the applicability of microcontinuum theories

,”

Int. J. Solids Struct.

41

,

2085

–

2097

(

2004

).

https://doi.org/10.1016/j.ijsolstr.2003.11.030

Google Scholar

Crossref

64.

Y.

Chen

,

J. D.

Lee

, and

A.

Eskandarian

, “

Examining the physical foundation of continuum theories from the viewpoint of phonon dispersion relation

,”

Int. J. Eng. Sci.

41

,

61

–

83

(

2003

).

https://doi.org/10.1016/S0020-7225(02)00141-6

Google Scholar

Crossref

65.

A. C.

Eringen

,

Mechanics of Micromorphic Continua

(

Springer

,

New York

,

1968

).

Google Scholar

Crossref

66.

A. C.

Eringen

,

Microcontinuum Field Theories

(

Springer

,

New York

,

1999

), pp.

101

–

248

.

Google Scholar

Crossref

67.

R. D.

Mindlin

, “

Micro-structure in linear elasticity

,”

Arch. Rational Mech. Anal.

16

,

51

–

78

(

1964

).

https://doi.org/10.1007/BF00248490

Google Scholar

Crossref

68.

J. D.

Lee

and

Y.

Chen

, “

Constitutive relations of micromorphic thermoplasticity

,”

Int. J. Eng. Sci.

41

,

387

–

399

(

2003

).

https://doi.org/10.1016/S0020-7225(02)00240-9

Google Scholar

Crossref

69.

Y.

Chen

,

J.

Lee

, and

A.

Eskandarian

, “

Atomistic counterpart of micromorphic theory

,”

Acta Mech.

161

,

81

–

102

(

2003

).

https://doi.org/10.1007/s00707-002-0981-2

Google Scholar

Crossref

70.

Y.

Chen

and

J. D.

Lee

, “

Connecting molecular dynamics to micromorphic theory. (I). Instantaneous and averaged mechanical variables

,”

Phys. A Stat. Mech. Appl.

322

,

359

–

376

(

2003

).

https://doi.org/10.1016/S0378-4371(02)01921-0

Google Scholar

Crossref

71.

Y.

Chen

and

J. D.

Lee

, “

Connecting molecular dynamics to micromorphic theory. (II). Balance laws

,”

Phys. A Stat. Mech. Appl.

322

,

377

–

392

(

2003

).

https://doi.org/10.1016/S0378-4371(02)01922-2

Google Scholar

Crossref

72.

Y.

Chen

and

J. D.

Lee

, “

Determining material constants in micromorphic theory through phonon dispersion relations

,”

Int. J. Eng. Sci.

41

,

871

–

886

(

2003

).

https://doi.org/10.1016/S0020-7225(02)00321-X

Google Scholar

Crossref

73.

X.

Chen

,

A.

Diaz

,

L.

Xiong

,

D. L.

McDowell

, and

Y.

Chen

, “

Passing waves from atomistic to continuum

,”

J. Comput. Phys.

354

,

393

–

402

(

2018

).

https://doi.org/10.1016/j.jcp.2017.10.038

Google Scholar

Crossref

74.

B. H.

Armstrong

, “

N processes, the relaxation-time approximation, and lattice thermal conductivity

,”

Phys. Rev. B

32

,

3381

(

1985

).

https://doi.org/10.1103/PhysRevB.32.3381

Google Scholar

Crossref

75.

J. G.

Kirkwood

,

F. P.

Buff

, and

M. S.

Green

, “

The statistical mechanical theory of transport processes. III. The coefficients of shear and bulk viscosity of liquids

,”

J. Chem. Phys.

17

,

988

–

994

(

1949

).

https://doi.org/10.1063/1.1747099

Google Scholar

Crossref

76.

M.

Born

and

H.

Green

, “

A general kinetic theory of liquids. III. Dynamical properties

,”

Proc. R. Soc. Lond. A

190

,

455

–

474

(

1947

).

https://doi.org/10.1098/rspa.1947.0088

Google Scholar

Crossref

2019

Author(s)

Concurrent atomistic-continuum modeling of crystalline materials

I. INTRODUCTION

II. THE IRVING-KIRKWOOD PROCEDURE

III. THE CAC FORMALISM

A. Ensemble-averaged local densities

B. Space and time-averaged local densities

C. Conservation equations and fluxes

D. Temperature

IV. FINITE ELEMENT FORMULATION

A. Temperature in FE formulation

B. Finite element equations

C. Fluxes in FEA solutions

V. NUMERICAL EXAMPLES

A. Atomic-scale structures of interfaces

B. Fracture and dislocations

C. Phonon thermal transport

VI. SUMMARY AND DISCUSSIONS

ACKNOWLEDGMENTS

APPENDIX A: CONSERVATION EQUATIONS FOR MONATOMIC CRYSTALS

1. Local density definition and the balance equation of linear momentum for monatomic crystals

2. Conservation equation of energy for monatomic crystals

APPENDIX B: CONSERVATION EQUATIONS FOR POLYATOMIC CRYSTALS

1. Balance equation of linear momentum for polyatomic crystals

2. Conservation equation of energy for polyatomic crystals

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Concurrent atomistic-continuum modeling of crystalline materials

I. INTRODUCTION

II. THE IRVING-KIRKWOOD PROCEDURE

III. THE CAC FORMALISM

A. Ensemble-averaged local densities

B. Space and time-averaged local densities

C. Conservation equations and fluxes

D. Temperature

IV. FINITE ELEMENT FORMULATION

A. Temperature in FE formulation

B. Finite element equations

C. Fluxes in FEA solutions

V. NUMERICAL EXAMPLES

A. Atomic-scale structures of interfaces

B. Fracture and dislocations

C. Phonon thermal transport

VI. SUMMARY AND DISCUSSIONS

ACKNOWLEDGMENTS

APPENDIX A: CONSERVATION EQUATIONS FOR MONATOMIC CRYSTALS

1. Local density definition and the balance equation of linear momentum for monatomic crystals

2. Conservation equation of energy for monatomic crystals

APPENDIX B: CONSERVATION EQUATIONS FOR POLYATOMIC CRYSTALS

1. Balance equation of linear momentum for polyatomic crystals

2. Conservation equation of energy for polyatomic crystals

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only