Notes for Pair-Potential Code Argon ver 2.1

Ju Li

April 3, 1999

1 Basic Setup

We will be concerned with particle systems that interact via short-ranged pair potential, with cutoff distance r_{cut}. Given any function W (r), such a potential can be constructed by

V (r) ≡ \{ \begin{matrix}W (r) − W (r_{cut}) − W′(r_{cut})(r − r_{cut}), & 0 < r < r_{cut} \\ 0, & r ≥ r_{cut} \end{matrix}  (1)

W (r) is then smooth in both value and first derivative at r = r_{cut}. Binary components are supported in the code, but they should only interact via scaled forms of a single V (r),

V_{ij}(r) = E_{ij}V (\frac{r}{L_{ij}} ) , 0 < r < ∞,  (2)

where the ij-pair is one of the three possibilities: Ar-Ar(0), Ar-NonAr(1), NonAr-NonAr(2), and E_{ij}, L_{ij} are dimensionless constants describing the differences with pure Ar-Ar interaction potential.

The supercell is designed to be a parallelepiped, which if desired can then be stacked infinitely in space. The three edges in the Cartesian frame (row vectors) are^{1}

a_{1} = (H_{11}, H_{12}, H_{13}), a_{2} = (H_{21}, H_{22}, H_{23}), a_{3} = (H_{31}, H_{32}, H_{33}),  (3)

with H_{ij} forming a matrix H. The position of a particle is specified by a row vector

s = (s_{1}, s_{2}, s_{3})  (4)

with s_{i}’s usually satisfying

0 ≤ s_{i} < 1, i = 1, 2, 3,  (5)

and the Cartesian coordinate of this particle x, also a row vector, is

x = s_{1}a_{1} + s_{2}a_{2} + s_{3}a_{3} = sH,  (6)

or in subscript notation,

x_{j} = ∑_{i=1}^{3} s_{i}H_{ij}, j = 1, 2, 3.  (7)

s_{i} has the geometrical meaning of the fraction of the ith edge to make up x.

The volume of the supercell is

Ω = | det H| = |a_{1} · (a_{2} × a_{3})| > 0.  (8)

The inverse of the H matrix B ≡ H^{−1}, such that

I = HB = BH.  (9)

If we define row vectors b_{i} = H_{ji},

b_{1} = (B_{11}, B_{21}, B_{31}), b_{2} = (B_{12}, B_{22}, B_{32}), b_{3} = (B_{13}, B_{23}, B_{33}),  (10)

then (9) is equivalent to

a_{i} · b_{j} ≡ a_{i}b_{j}^{T} = δ_{ij}.  (11)

1This labelling is literally followed in all my C (1−3 becomes 0−2) and Fortran source codes, irrespective of internal storage arrangements. That is, I may have sacrificed some efficiency for clarity in Fortran, where columns are stored contiguously) .

Since b_{1} is perpendicular to both a_{2} and a_{3}, it must be collinear with the normal vector n of the a_{2}/a_{3} plane,

b_{1} = |b_{1}|n,  (12)

and so by (11),

1 = a_{1} · b_{1} = a_{1} · (|b_{1}|n) = |b_{1}|(a_{1} · n),  (13)

but |a_{1} · n| is nothing other than the thickness of the entire supercell along the a_{1} edge. Thus, the thicknesses (distances between two parallel surfaces) of the supercell are

d_{1} = \frac{1}{|b_{1}|}, d_{2} = \frac{1}{|b_{2}|}, d_{3} = \frac{1}{|b_{3}|}.  (14)

This result will be useful to determine how large the bins should be in order for a particle to interact with particles in just 27 adjacent bins.

We shall use the Lagrangian strain \underline{\underline{η}} as a measure of deformation of the supercell, where a reference state H_{0} is chosen at a prior moment,

dx = (ds)H_{0},  (15)

and with ds kept constant, the material point is carried to dx′ by changing H_{0} → H,

dx′ = (ds)H.  (16)

The Lagrangian strain is defined by the change of differential line length,

dl^{′2} = dx′dx^{′T} ≡ dx(I + 2\underline{\underline{η}})dx^{T} ,  (17)

which, by plugging in

dx′ = (ds)H = (dx)H_{0}^{−1}H,  (18)

can be easily shown to be

\underline{\underline{η}} = \frac{1}{2} (H_{0}^{−1}HH^{T} H_{0}^{−T} − I) .  (19)

Because \underline{\underline{η}} is a symmetric 3 × 3 matrix, it always has three mutually orthogonal eigendirections, x_{1}, x_{2}, x_{3},

x_{1}\underline{\underline{η}} = x_{1}η_{1}, x_{2}\underline{\underline{η}} = x_{2}η_{2}, x_{3}\underline{\underline{η}} = x_{3}η_{3},  (20)

and along those directions, the line segment lengths change by

dl′ = dl\sqrt{1 + 2η_{i}},  (21)

which can be shown to achieve extrema among all directions as three poles of an ellipsoid. Thus, as long as η_{1}, η_{2} and η_{3} oscillate between (−η_{max}, η_{max}) for some chosen constant η_{max}, any line segment at H_{0} lengthens by no more than \sqrt{1 + 2η_{max}} and shortens by no less than \sqrt{1 − 2η_{max}}. This condition defines a strain session. At the moment when the criterion

|η_{1}| < η_{max}, |η_{2}| < η_{max}, |η_{3}| < η_{max}  (22)

is violated, the strain session terminates. And we usually pick the current H as the new H_{0} to start off a new session, which is called a re-mesh incident. During any one session, distance between any two points with fixed Δs can change by no more than factors of

(f^{−1}, f),  where f ≡ \sqrt{\frac{1 + 2η_{max}}{1 − 2η_{max}}} .  (23)

And the supercell thicknesses shorten by no less than \sqrt{1 − 2η_{max}} from the values at H_{0}.

2 Maintaining a Neighbor List with O(N) Computa-

tion

Since interaction between atom i and j is truncated at L_{ij}r_{cut}, ideally we would have to do a calculation only when

r_{ij} ≡ |x_{ij}| < L_{ij}r_{cut}.  (24)

However, without compartmentalizing the supercell, checking out the above condition is itself an O(N^{2}) calculation. The idea is that we should subdivide the supercell into m_{1} × m_{2} × m_{3} smaller bins: each bin maintains a list of the particles inside, and the bin size is chosen to approximate L_{max}r_{cut}, such that a particle only needs to be concerned with other particles in nearby 27 bins, as possible candidates for interaction.

Upon it, we further build and maintain a smaller, “personal” list for each particle. That is, a particle does not need to look over 27 bins at each step, which is O(N) effort but nonetheless burdensome, but only needs to do that quite infrequently at so-called re-anchor events; in all other steps it only needs to consult its own personal list. It is exactly at those re-anchor events that a personal list is modified. Re-anchor events are usually asynchronous, meaning different particles do it at different times; except at the re-mesh incident when all particles gets re-anchored, essentially resetting the system. Depending on whether the supercell undergoes significant deformation (η_{max} is usually chosen to be ∼ 0.02), re-mesh incidents may not happen at all during a simulation.

To make things easier, we define a unique “caretaker” for the ij-interaction between atoms i or j. For load-balancing purposes, it is best if the responsibilities are evenly distributed. We use the pseudo-random choice of “if i + j is odd and i > j, or if i + j is even and i < j, then i is the caretaker; otherwise it is j.” Because i > j is “uncorrelated” with whether i + j is even, there is no significant load concentration for most geometrical configurations. A particle’s personal list only keeps possible interactions that are under its own care.

A re-anchor event works like this: particle i refreshes its own anchor position s_{i}^{a} – a reduced coordinate (instead of Cartesian one) which stays constant at other times,

s_{i}^{a} = s_{i}, at reanchor,  (25)

and finds out (in nearby 27 bins) all particles with anchor position s_{j}^{a} satisfying

|R_{ij}^{∗}| = |Δs_{ij}^{a}H^{∗}| ≤ L_{ij}r_{list},  (26)

where the superscript ∗ is to denote this particular re-anchor moment, and r_{list} is a constant to be determined. If this ij-interaction belongs to i, it is put in i’s list; otherwise it is added to j’s list. As for “obsolete” interactions, those used to be taken care by i are liquidated at this particular re-anchor; but those taken care by old j’s become zombies, which in this version of the code is not cleared at i’s re-anchor for time-efficiency concerns, since it will be cleared at those j’s re-anchors.

The role of particle i’s anchor, s_{i}^{a}, is to localize that particle,

|d_{i}| = |(s_{i} − s_{i}^{a})H| ≤ L_{i}^{min}βr_{cut},  (27)

within distance L_{i}^{min}βr_{cut}, where

L_{i}^{min} ≡ min (L_{ij}) = min (L_{i,Ar}, L_{i,NonAr}) = min (L_{ii}, L_{Ar,NonAr}) .  (28)

And β is a constant usually chosen to be ∼ 0.08. Whenever condition (27) is violated, it triggers a re-anchor event for particle i.

Now we will prove the central result of this section, that is, as long as the twin conditions (26) and (27) are satisfied during one strain session, there exist a choice of r_{list} that guarantees all pair interactions with

r_{ij} ≡ |Δs_{ij}H| ≤ L_{ij}r_{cut}  (29)

are recorded in its caretaker’s list. For there is,

Δs_{ij}H = Δs_{ij}^{a}H + (s_{i} − s_{i}^{a})H − (s_{j} − s_{j}^{a})H ≡ R_{ij} + d_{i} − d_{j},  (30)

i.e., the current displacement between particles i and j are the sum of their current anchor displacement R_{ij}, and current drifts d_{i}, d_{j} from individual anchors. Because (26) is satisfied at the moment of anchoring (∗), those not in particle i’s list must have

|R_{ij}^{∗}| > L_{ij}r_{list},  (31)

and by virtue of (23), during the entire strain session must satisfy

|R_{ij}| ≥ f^{−1}|R_{ij}^{∗}| > f^{−1}L_{ij}r_{list}.  (32)

Also, since L_{i}^{min}, L_{j}^{min} ≤ L_{ij} for a particular ij, by (27) there are,

|d_{i}| ≤ L_{i}^{min}βr_{cut} ≤ L_{ij}βr_{cut},  (33)
|d_{j}| ≤ L_{j}^{min}βr_{cut} ≤ L_{ij}βr_{cut}.  (34)

Thus, applying the inequality

r_{ij} = |R_{ij} + d_{i} − d_{j}| ≥ |R_{ij}| − |d_{i}| − |d_{j}|,  (35)

and canceling common factor L_{ij} from both sides in (29), one obtains

r_{ij} > L_{ij}r_{cut}  (36)

during the entire strain session for those unrecorded particles at ∗ if only

f^{−1}r_{list} − 2βr_{cut} ≥ r_{cut},  (37)

or, finally,

r_{list} ≥ f(1 + 2β)r_{cut} ≡ \sqrt{\frac{1 + 2η_{max}}{1 − 2η_{max}}} (1 + 2β)r_{cut}  (38)

for any sensible choice of η_{max} and β. The conditions (26) and (27) are initialized at the start of each strain session by re-anchoring all particles.

Let us next look at how to choose the bin partition m_{1} × m_{2} × m_{3} at every re-mesh incident, such that conditions (26) can be implemented efficiently for each particle during the entire strain session, by just looking into nearby 27 bins; since a particular bin partition cannot change during a strain session. It would work if the bin thicknesses t_{i} in all three directions are always ≥ L_{ij}r_{list} during the session, which is equivalent to demanding that

t_{1}, t_{2}, t_{3} ≥ L_{max}r_{list}  (39)

where

L_{max} = max (L_{Ar−Ar}, L_{Ar−NonAr}, L_{NonAr−NonAr}) ,  (40)

if there are non-Ar atoms in the system; and t_{i}’s are the current thicknesses of bins. So, the bin thicknesses should be still larger at the re-mesh incident as the supercell might shrink during the strain session. There should be

t_{1}^{0}, t_{2}^{0}, t_{3}^{0} ≥ \frac{L_{max}r_{list}}{\sqrt{1 − 2η_{max}}}  (41)

at the re-mesh incident to guarantee that (39) always hold. And so, one should partition the supercell using

m_{i} ≤ \frac{\sqrt{1 − 2η_{max}}}{|b_{i}^{0}|L_{max}r_{list}} , i = 1, 2, 3,  (42)

based on the supercell thickness formula (14).

The above derivations are fully rigorous in continuous time, but a tiny error may be introduced when time-stepping, where interacting particles may re-anchor at the same step. There are sufficient leeway built into the derivation that makes such a problem highly improbable and insignificant. On the other hand, if we did make an error in either the derivation or the coding, the error will manifest itself since the many particles are always testing these inequalities (possible equals) inexorably during the simulation.

For calculating the Radial Distribution Function g(r), we often demand greater r ranges than merely r_{cut} or r_{list}. There, we can utilize the property (39) of bins, which always holds during the entire strain session, and defensively choose the number of relevant bins in order to accomplish r_{gcut}. However, one should avoid the pitfall where the same bin are selected twice if the system turns out to be small.

3 Longwave Instability in RVE: an Upper Bound to

Material Strength

Imagine we gradually increase the load on a homogeneous material without dislocation density, either under constant displacement, constant traction or constant stress (whatever it means) – the effect is the buildup of elastic energy in the system, which it dearly wants to get rid of, either by^{2} nucleation and growth of void(s), extension of existent crack(s), or, if all fails, the initiation of long-wave instability (LWI). Since LWI does not demand any prerequisites on the material, such as the presence of a crack or a void (the nucleation of which requires a critical r_{c} and thermally activated transition of a barrier), it can and will happen in a homogeneous material, thus putting an upper bound on the RVE’s strength.

2There are other modes whereby the elastic energy can be released, for instance emission of dislocation loops by Frank-Read source, GB sliding etc.; for clarity of this discussion we shall omit them.

That is, even if the good material is lucky enough not to have a crack, and not to nucleate a void during the loading, it cannot escape being ruined by LWI at high enough stress, usually on the order of 0.025-0.1E, which is the same order as σ_{IC}, the interfacial decohesion stress from UBER.

We shall demonstrate that LWI is intrinsically related to σ_{IC}, and also the apex of so-called τ-γ curve of “shear transformation element”, important for inelastic deformation. They, viewed from the perspective of LWI, are the two faces of the same mirror: one comes from unstable longitudinal waves that do not preserve volume, the other from unstable transverse waves that produces shear. What we strive to achieve in this writeup are not new insights into the grand theory of materials deformation – features of UBER or τ-γ curves are wellknown and well-theorized – but emphasizing on practical issues like how we can calculate such saddle-points/apexes from first principles in an atomistic simulation. We are going to do a rigorous analysis at finite temperature (unlike UBER) using Lagrangian strain, as the loading is no longer small (indeed it is huge); and obtain exact particle formulas (usually called fluctuation formulas) for τ_{ij} and elastic constants C_{ijkl}, defined based on such, which can be readily evaluated in a microcanonical-ensemble MD simulation. That is, we are able to not only observe LWI directly in a simulation, but are also able to fully characterize it, and verify our little theory.

The “long-wave” is somewhat a misnomer – the instability has been observed in MD simulations in 15Ȧ cells, consisting of 200 atoms, and the critical stress one obtains does not differ significantly with a 100Ȧ simulation. Like the Navier-Stokes equation or elasticity theory, it might be pushed down to amazingly small scales. So, LWI is a valid theory for “activation volume elements” – a tiny microscopic region that is under very high stress. It is quite important: I think it is one of the two ultimate, smallest-scale, progenitors of material deformation (the other one is Frank-Read source), since the end product of LWI, besides carrying strain, are also observed to result in voids and shear faults, which then might be turned into micro-cracks and dislocations. Also, this method may be generalizable to inhomogeneous systems such as grain boundaries, or hard particle interfaces, although the analysis would be much harder.

The story is a simple one. Consider a homogeneous material at configuration X with constant temperature T: if one slightly distort the configuration by affine transformation J,

X → Y = JX, J_{ij} ≡ δ_{ij} + u_{i,j}, u_{i,j} ≡ \frac{∂u_{i}(x)}{∂x_{j}},  (43)

the Lagrangian strain of state Y with respect to X are defined as

η ≡ \frac{1}{2}(J^{T} J − 1) = \frac{1}{2}(u_{i,j} + u_{j,i} + u_{k,i}u_{k,j}).  (44)

The reasons for using η instead of the more commonly used ϵ

ϵ_{ij} ≡ \frac{1}{2}(u_{i,j} + u_{j,i})  (45)

are,

• η_{ij}, unlike ϵ_{ij}, is a genuine second-order tensor satisfying tensor transformation law. That remains true even when u_{i,j} is no longer small.

• The volume element’s free energy depends only on η, which has six degrees of freedom, instead of J, which has nine.

• We could obtain explicit finite temperature particle formulas for τ_{ij} and C_{ijkl}, defined based on η.

Since the Helmholtz free energy of state Y, F (Y), is rotationally invariant,

F (Y) = F (RY),  (46)

where R is an arbitrary rotation, F (Y) is a function of just η,

F (Y) = F (η, X),  (47)

and one can always do a Taylor expansion on F (η, X) with respect to η when η is small:

F (Y) = F (X) + V (X) (τ_{ij}(X)η_{ij} + \frac{1}{2}C_{ijkl}(X)η_{ij}η_{kl} + ...) .  (48)

Now let us consider a volume element of originally homogeneous material with fixed outer boundary X. The materials inside remain homogeneous at low stress levels, not because “it choose to”, but because either it costs more energy to generate certain inhomogeneities inside, such as a wave, or because the transition barrier to creating a singularity inside such as a void (if it does reach there, F will be lowered) is too high at the current T for it to happen. Suppose the second option is suppressed, the volume element would remain homogeneous up to a point where generating a wave costs no energy. At first glance of (48), one may suggest that it depends on whether the matrix C_{ijkl} is positive definite or not. However, because of the compatibility relation, one cannot freely assign any strain wave he desires – it must come from a genuine displacement wave, u(x). And since η in (44) contains second-order terms of u, we see that we has slightly more work to do. The final positive definiteness relation will involve both τ and C.

By long-wave, we assume that the system free energy is additive, i.e., it is the algebraic sum of free energies of constituent parts, each undergoing “homogeneous deformation” in its own differential volume. In effect it is to replace V (X) in (48) by an integral,

V (X) ⟶ ∫_{V (X)} d^{3}x,  (49)

with spatially-varying η(x) generated by u(x). If we take u(x) to be a plane wave,

u(x) = we^{ik·r}  (50)

with wave-vector k satisfying the Born-von Karman boundary condition

k_{i} = \frac{2πn_{i}}{L_{i}} , i = 1, 2, 3,  (51)

where L_{i} is the dimension of the system in the ith direction, u(x) will vanish on the boundary of the sample, thus still satisfying the constraint X. However the material is no longer homogeneous. If the energy cost of creating such a wave is equal to or less than zero, then the material is unstable at X, and the waves will develop into a catastrophic event.

For the initiation of such an instability, |u| can be taken to be very small, and it suffices to do expansion to the lowest effective order,

ΔF_{wave} = ∫_{V (X)} d^{3}x (τ_{jl}(X)η_{jl}(x) + \frac{1}{2}C_{ijkl}(X)η_{ij}(x)η_{kl}(x)) + O(u^{3})  (52)

When we plug in

η_{jl}(x) = \frac{1}{2} (u_{j,l}(x) + u_{l,j}(x) + u_{k,j}(x)u_{k,l}(x))  (53)

and do the integration, first-order contributions of u vanish because the integral of any sinusoidal function over its period must be zero. We are thus left with u^{2} contributions, and keeping to that order of accuracy, we obtain

ΔF_{wave} = \frac{1}{2} ∫_{V (X)} d^{3}x (τ_{jl}(X)u_{k,j}(x)u_{k,l}(x) + C_{ijkl}(X)u_{i,j}(x)u_{k,l}(x)) + O(u^{3}),  (54)

which can be re-written as

ΔF_{wave} = \frac{1}{2} ∫_{V (X)} d^{3}xD_{ijkl}(X)u_{i,j}(x)u_{k,l}(x),  (55)

where

D_{ijkl}(X) ≡ C_{ijkl}(X) + τ_{jl}(X)δ_{ik}.  (56)

D is a somewhat strange-looking rank-4 tensor: it has ij ↔ kl symmetry, but not i ↔ j or k ↔ l symmetry we are used to in linear elasticity. A close examination suggests that it is caused by the asymmetric positions of displacements (k) and derivatives (i, j) in the second-order u_{k,i}u_{k,j} term in (44). This means that at finite stress, a wave of w = a and k = b does not in general have the same energy as a wave of w = b and k = a, and this asymmetry is the more apparent when we approach the instability, since in all other loading conditions, τ ≪ C.

Inserting (50) into (54), we see that the stability criterion is

Γ(w, k) = (C_{ijkl}w_{i}w_{k} + τ_{jl})k_{j}k_{l} > 0, |w| ≡ 1, |k| ≡ 1,  (57)

for all possible orientations of w and k. One method of finding out whether (57) holds is to numerically tabulate the solid angles of w in θ = (0, π) and ϕ = (0, π) (since w ↔ −w in (57)), evaluate the following 3 × 3 symmetric matrix,

γ_{jl}(w) ≡ C_{ijkl}w_{i}w_{k} + τ_{jl},  (58)

and calculate the determinant or the smallest eigenvalue of γ_{jl}(w). When that surface touches 0, the homogeneous volume element is no longer stable.

We can look into some special cases analytically. For isotropic (such as amorphous) materials, the elastic constant keeps the form

C_{ijkl} = λδ_{ij}δ_{kl} + μδ_{ik}δ_{jl} + μδ_{il}δ_{jk},  (59)

even at high stress, and there is, by (58),

γ_{jl}(w) = (λ + μ)w_{j}w_{l} + μδ_{jl} + τ_{jl}.  (60)

In general, for that material to keep being isotropic, the applied stress must be hydrostatic

τ_{jl} = tδ_{jl},  (61)

and so the stability criterion is easily seen to be

λ + 2μ + t > 0 (longitudinal), μ + t > 0 (transverse).  (62)

One see that for failure under tension, the elastic constant have to go negative, whereas linear elastic theories predicts that λ = 0 or μ = 0 would be sufficient.