' Ω(X){Cijkl ηkl − Ukk τij + Uik τkj + τik Ujk }1' Ω(X){ Cijkl (Ukl + Ulk ) − Ukk τij + Uik τkj + τik Ujk }2= Ω(X){Cijkl Ukl − Ukk τij + Uik τkj + τik Ujk }(2.269)So∗
gij∗ = Ω(X)Bijkl
Ukl + ...

(2.270)∗Bijkl= Cijkl − τij δkl + τjl δik + τil δjk(2.271)

withSimilar to the arguments in symmetric deformation space, stability behavior ingeneral deformation space is governed byA∗ = (B ∗ + (B ∗ )T )/2(2.272)

det |A∗ | = 0(2.273)

Whenfor the first time, the system become unstable. Note that A∗ is a 9 × 9 matrix.
As an observation, when a configuration is stable in general deformation space,i.e., none of the eigenvalues of A∗ are negative, the configuration would be stable inall deformation spaces.

2.4.2

Deformation of a Periodic Simulation Cell

Basic setup
We will simulate particle systems that interact via short-ranged potentials of cutoff
radius rcut . Given any radial function w(r), this can be done by redefining a w̃(r),w̃(r) ≡ w(r) − w(rcut ) − w 0 (rcut )(r − rcut ),0 < r < rcut
 0,r ≥ rcut

87,

(2.274)

which2 ensures continuity in both energy and forces as a particle crosses r = rcut . Inαβ
the case of multiple-component systems, rcut is generalized to a matrix rcut
, whereα ≡ c(i), β ≡ c(j) are the chemical types of atoms i,j. As a convention, i suggests
a “host” atom at the “frame origin”, whereas j suggests a “client” atom. We thendefine
xij ≡ xj − xi ,rij ≡ |xij |,³x̂ij ≡xij.
rij(2.275)´

Quantities such as the pair force fij ≡ −Ṽ 0 (rij ) x̂ij are understood as the force on
j due to i. fij should be parallel to xij when the potential is repulsive, a mnemonicdevice.
The supercell is a parallelepiped, which can be tiled in space indefinitely if desired.
The three edges are row vectors
a1 = (H11 , H12 , H13 ), a2 = (H21 , H22 , H23 ), a3 = (H31 , H32 , H33 ),(2.276)in Cartesian coordinates, with Hµν forming a 3×3 matrix H. 3 The position of particle
i is specified by a row vector, si = (si1 , si2 , si3 ), with siµ ’s usually satisfying
0 ≤ siµ < 1,µ = 1..3,(2.277)and the Cartesian coordinate of this particle, xi , also a row vector, isxi = si1 a1 + si2 a2 + si3 a3 = si H,(2.278)where siµ has the geometrical interpretation of the fraction of the µth edge to build
2An alternative is to define w̃(r) ≡ w(r) exp(rs /(r − rcut )) which has all derivatives continu¡ ¢12 ¡ σ ¢6− r +ous. Anotherefficient scheme for the LJ6-12 potential is w(r) ≡ 4ε[ σrµ ³at r´= rcut¶´12 ¡ ¢³´12´6³³186σσσσr
2 rcut− 3 rcut− rcut+ 2 rcut], which expands in r6 instead of r, and avoids
σusing sqrt or exp.
3This labelling scheme is literally followed in both my C (1 − 3 becomes 0 − 2) and Fortran source
codes, irrespective of internal storage arrangements. That is, I may sacrifice efficiency for clarity
in Fortran (where columns are stored contiguously), to achieve easy-to-read correspondence with C
source codes, where this arrangement is computationally more efficient for edge vector operations.

88

xi .
The volume of the supercell isΩ = | det H| = |a1 · (a2 × a3 )| > 0.

(2.279)The inverse of the H matrix B ≡ H−1 satisfies
I = HB = BH.

(2.280)If we define row vectors
b1 ≡ (B11 , B21 , B31 ), b2 ≡ (B12 , B22 , B32 ), b3 ≡ (B13 , B23 , B33 ),(2.281)then (2.280) is equivalent toai · bj ≡ ai bTj = δij .

(2.282)Since b1 is perpendicular to both a2 and a3 , it must be collinear with the normal
direction n of the a2 /a3 plane: b1 ≡ |b1 |n. And so by (2.282),
1 = a1 · b1 = a1 · (|b1 |n) = |b1 |(a1 · n).

(2.283)

But |a1 · n| is nothing other than the thickness of the supercell along the a1 edge.
Therefore, the thicknesses (distances between two parallel surfaces) of the supercell
are,
d1 =1,|b1 |d2 =1,|b2 |d3 =1.
|b3 |(2.284)

The general design of the simulation should allow for deformation of H that includes rotational components, even though one may choose to impose the constraintof symmetric deformation later, whose dynamics is derived in section 2.4.2. In general one should use the Lagrangian strain η, a true rank-2 tensor under coordinatetransformation, to measure the deformation of a MD supercell as it is unlikely to be89

infinitesimal. To define η, one needs a reference H0 of a previous time, with x0 = sH0and dx0 = (ds)H0 , and then imagine that with s and ds fixed, dx0 is transformed todx = (ds)H, under H0 → H ≡ H0 K.
The Lagrangian strain is defined by the change in the differential line length,dl2 = dxdxT ≡ dx0 (I + 2η)dx0T ,(2.285)

where by plugging in dx = (ds)H = (dx0 )H−1
0 H = (dx0 )K, η is seen to beη=´´
1 ³ −1
1³T
H0 HHT H−TI
KK−=−I.
022(2.286)

Because η is a symmetric matrix, it always has three mutually orthogonal eigendirections x1 η = x1 η1 , x2 η = x2 η2 , x3 η = x3 η3 . Along those directions, the line√√√
lengths are changed by factors 1 + 2η1 , 1 + 2η2 , 1 + 2η3 , which achieve extrema among all line directions. Thus, as long as η1 , η2 and η3 oscillate between[−ηbound , ηbound ] for some chosen ηbound , any line segment at H0 can be lengthened by√√
no more than 1 + 2ηbound and shortened by no less than 1 − 2ηbound . That is, if
we define length measureL(∆s, H) ≡√∆sHHT ∆sT ,(2.287)

then so long as η1 , η2 , η3 oscillate between [ηmin , ηmax ], there isqq

1 + 2ηmin L(∆s, H0 ) ≤ L(∆s, H) ≤1 + 2ηmax L(∆s, H0 ).

(2.288)One can use the above result to define a strain session, which begins with H0 ≡ Hand during which no line segment is allowed to shrink by less than a threshold fc ≤ 1,
compared to its length at H0 . This is equivalent to requiring that,qf≡

1 + 2 (min(η1 , η2 , η3 ))≤fc .

(2.289)Whenever this condition is violated, the session terminates and a new session starts
90with the current H as the new H0 . In my implementation of the O(N ) moleculardynamics program, this is associated with a repartitioning of the supercell into equalsize bins, and is called a strain-induced bin repartitioning.
The purpose of a bin partition and a strain session is the following: it can be avery demanding task to determine if atoms i,j are within rcut or not, for all possible
ij combinations. 4 Formally, this requires checking
rij ≡ L(∆sij , H) ≤ rcut .

(2.290)Because si , sj and H are mobile – they differ from step to step, it appears that wehave to do this at each step. This O(N 2 ) complexity would indeed be the case but
for the observation that, in most MD, MC and static minimization applications, si ’sand H often change only a little from the previous step. Therefore, once we ensured
that (2.290) hold at some previous step, we can devise a sufficient condition to test
whether (2.290) still must hold now or not, at a much smaller cost. Only when thissufficient condition breaks down, which is taken to be less frequent, do we resort to amore complicated search and check in the fashion of (2.290).
My implementation of the above idea is as follows: I associate each si with asemi-mobile reduced coordinate sai called atom i’s anchor. At each step, I check if
L(si − sai , H), that is, the current distance between i and its anchor, is greater than0a certain rdrift ≥ rdrift
or not. If it is not, then sia does not change; if it is, then I

redefine sai ≡ si at this step, which is called atom i’s flash incident. At i’s flash, atom
4
It is often more efficient to count pairs if the potential function allows for easy use of such “halflists”, such as pair- or EAM potentials. In these scenarios we pick a unique “caretaker” among i andj to store the information about the ij-pair, that is, a particle’s “personal list” only keeps possible
pairs that are under its own care. For load-balancing it is best if the responsibilities are distributed
evenly among particles. We use a pseudo-random choice of “if i + j is odd and i > j, or if i + j iseven and i < j, then i is the caretaker; otherwise it is j.” As i > j is “uncorrelated” with whether
i + j is even or odd, significant “load imbalance” is unlikely to occur even if the indices correlatestrongly with the atoms’ positions.

91

i is required to update records 5 of all atoms whose anchors satisfy
L(saj − sai , H0 ) ≤ rlist ≡0rcut + 2rdrift.
fc(2.291)

Note that the distance is between anchors instead of atoms 6 , and the length ismeasured by H0 , not the current H. (2.291) nominally takes O(N ) work per flash,but I accelerate it to O(1) per flash by partitioning the supercell into m1 × m2 × m3
bins at the start of the session, whose thicknesses by H0 (see (2.284)) are required to
be greater than or equal to rlist :d1 (H0 ) d2 (H0 ) d3 (H0 ),,≥ rlist .
m1m2m3(2.292)

The bins deform with H and remains commensurate with it, that is, its s-partition
1/m1 , 1/m2 , 1/m3 remains fixed during a strain session. Each bin keeps an updated
list7 of all anchors inside. Then, if at the time of i’s flash two anchors are separated
by more than one bin, there would beL(saj − sai , H0 ) >d1 (H0 ) d2 (H0 ) d3 (H0 ),,≥ rlist ,m1
m2
m3(2.293)

and they cannot possibly satisfy (2.291). Therefore we only need to test (2.291)
for anchors within adjacent 27 bins. To synchronize, all atoms flash at the start of astrain session. From then on, atoms flash individually whenever L(si −sai , H) > rdrift .8
We see that to maintain anchor lists that captures all solutions to (2.291) among thelatest anchors, it takes only O(N ) work per step, and the pre-factor of which is also0
small because flash events happen quite infrequently for a tolerably large rdrift.

The central claim of the scheme is that if j is not in i’s anchor records (suppose5

Parts of the records may be stored in j’s if pairs are counted and j happens to be the caretakerof the ij pair.
6 asi = si , though.
7When atom i flashes, it also updates the bin-anchor list if necessary.
8If two anchors flash at the same step in a loop, the first flash may get it wrong – that is, missing
the second anchor, but the second flash will correct the mistake. The important thing here is to not
lose an interaction.

92

i’s last flash is more recent than j’s), which was created some time ago in the strain
session, then rij > rcut . The reason is that the current separation between the anchor
0i and anchor j, L(saj − sia , H), is greater than rcut + 2rdrift, since by (2.288), (2.289)

and (2.291),L(sja − sia , H) ≥ f · L(saj − sai , H0 ) > f · rlist ≥ fc · rlist = fc ·0
rcut + 2rdrift. (2.294)
fc

In fact, we see that the rij > rcut conclusion maintains if neither i or j currently
drifts more thanrdrift ≡f · rlist − rcut0≥ rdrift,2(2.295)

from respective anchors. Put it another way, when we design rlist in (2.291), we take
into consideration both atom drifts and H shrinkage which both may bring ij closer
than rcut , but since the current H shrinkage has not yet reached the designed criticalvalue, we can convert it to more leeway for the atom drifts.
For multi-component systems, we define
αβrlistαβ0rcut+ 2rdrift,≡fc(2.296)00
where both fc and rdrift
can be thought
are species-independent constants, and rdriftof as putting a lower bound on rdrift , so flash events cannot occur too frequently – aself-protection mechanism. At each bin repartitioning, we would required1 (H0 ) d2 (H0 ) d3 (H0 )αβ,,≥ max rlist.
α,βm1
m2m3(2.297)

And during the strain session, f ≥ fc , we have
"αrdrift

Ã!

Ãαββα
αββα
f · rlistf · rlist− rcut− rcut≡ min min, minββ
22!#,(2.298)

a time- and species-dependent atom drift bound that controls whether to flash.

93