A numerical-continuation-enhanced flexible boundary condition scheme applied to Mode I and Mode III fracture
AA numerical-continuation-enhanced flexible boundary condition schemeapplied to Mode I and Mode III fracture
Maciej Buze ∗ School of Mathematics, Cardiff University, Senghennydd Road, Cardiff, CF24 4AG, United Kingdom
James R. Kermode † Warwick Centre for Predictive Modelling, School of Engineering,University of Warwick, Coventry CV4 7AL, United Kingdom (Dated: September 1, 2020)Motivated by the inadequacy of conducting atomistic simulations of crack propagation using static boundaryconditions that do not reflect the movement of the crack tip, we extend Sinclair’s flexible boundary conditionalgorithm [
Philos. Mag. , 647–671 (1975)], enabling full solution paths for cracks to be computed withpseudo-arclength continuation, and present a method for incorporating more detailed far-field information intothe model for next to no additional computational cost. The new algorithms are ideally suited to study details oflattice trapping barriers to brittle fracture and can be easily incorporated into density functional theory and mul-tiscale quantum/classical QM/MM calculations. We demonstrate our approach for Mode III fracture with a 2Dtoy model, and for Mode I fracture of silicon using realistic interatomic potentials, highlighting the superiorityof the new approach over employing a corresponding static boundary condition. In particular, the inclusion ofnumerical continuation enables converged results to be obtained with realistic model systems containing a fewthousand atoms, with very few iterations required to compute each new solution. We also introduce a methodto estimate the lattice trapping range of admissible stress intensity factors K − < K < K + very cheaply anddemonstrate its utility on both the toy and realistic model systems. I. INTRODUCTION
Atomistic modelling of brittle fracture in crystals goes backto the pioneering work carried out by Sinclair and coworkersin the 1970s [1–3]; for a recent review of contributions madeto our understanding of fracture from atomistic simulationssee Refs. 4 and 5. The principal distinction from continuummodels is the discreteness of the atomic lattice, which leadsto the concept of lattice trapping, first identified by Thom-son in 1971 [6]. A consequence of lattice trapping is thatcracks remain stable over a range of stress intensity factors K − < K < K + . Lattice trapping can lead to anisotropy inpropagation directions [7], and the associated energy barri-ers imply that cleavage does not necessary produce smoothfracture surfaces at low energies [8]. The phenomenon hasa dynamical analogue, the velocity gap, which is a forbiddenband of crack velocities at low temperatures [9]. The velocitygap vanishes at larger temperatures; thermal activation overlattice trapping barriers has been proposed as an explanationfor observations of low speed crack propagation on the ( ) cleavage plane in silicon [10].Detailed investigation of these phenomena are currently ex-tremely challenging for two interconnected reasons. Firstly,realistic interatomic potentials capable of describing the veryhigh strains near crack tips are very hard to construct [11]. Atthe same time, the requirement for large model systems andthe strong coupling between lengthscales associated with frac-ture make the applicaton of quantum mechanical techniquessuch as density functional theory (DFT) extremely challeng- ∗ [email protected] † [email protected] ing, despite the considerable success such techniques have en-joyed elsewhere in materials science [12]. Even when accu-rate atomistic models are available, determining the relevantstable crack tip configurations and the energy pathways thatlink them is extremely challenging because of the high dimen-sionality of the atomistic configuration space [10]. The pictureis further complicated if the modelled crack propagates, caus-ing an effective shift of origin of the entire strain field, whichis often not reflected in the supplied boundary condition.In the mathematical sciences, numerical continuation tech-niques have been previously applied to study crack propaga-tion in Refs. 13 and 14. Relatedly, a mathematically rigorousnumerical analysis of domain size effects for a static boundarycondition scheme coupled with numerical continuation tech-niques has been conducted in Ref. 15. It provides a basicframework in which proving convergence rates to the infinitelimit is possible and is the principal motivation for the currentwork.Here, we return to the flexible boundary condition (FBC)approach introduced by Sinclair [3]. In this class of ap-proaches, a localised atomic core region is coupled to a linearelastic far-field. The method has been developed and appliedextensively to model dislocations [16–18], but applications tofracture have received comparatively little attention. We aug-ment the FBC method, and address the challenge of identify-ing and analysing stable and unstable crack tip configurationsby combining it with numerical continuation techniques. Wedemonstrate our ideas firstly for Mode III fracture in a toymodel of a 2D crystal, considered a useful stepping stone forour theory, as it has readily calculable exact Hessians and per-mits a mathematically rigorous analysis. This is then followedby the more realistic and much-studied example of Mode Ifracture of silicon on the ( ) cleavage plane, using bondorder potentials that have been modified to extend the interac- a r X i v : . [ phy s i c s . c o m p - ph ] A ug tion range and introduce screening to provide a qualitativelycorrect description of bond-breaking processes [19]. II. METHODOLOGYA. Discrete Kinematics
For the purpose of describing the method, we first considera simplified system consisting of a two-dimensional infinitecrystal of atoms forming a triangular lattice and interacting viaa known interatomic potential with a finite interaction radius,with a crack forming along the horizontal axis and a cracktip located at ( α , ) ∈ R . The position of the i th atom isdenoted by x ( i ) = ( x ( i ) , x ( i ) , x ( i )) ∈ R , which is alwaysof the form x ( i ) = ˆ x ( i ) + Y ( i ) , (1)where ˆ x ( i ) is the crystalline lattice position and Y ( i ) = ( Y ( i ) , Y ( i ) , Y ( i )) is the displacement from thecrystalline lattice.The theory will be presented for two crack modes:pure Mode III, described by out-of-plane displacements( Y = Y = Y = φ , which, to avoid unnecessarytechnicalities is initially taken to be a pair potential with totalenergy of the form E = ∑ i (cid:54) = j φ ( r i j ) , where r i j = | x ( i ) − x ( j ) | . (2)The restriction to pair potentials is not needed for the analy-sis, and will be lifted in the numerical examples considered inSection III B.Following the ideas of Sinclair [3], the system is dividedinto three regions. Region I, also known as the defect core isa finite collection of atoms, labelled from i = i = N , inthe vicinity of the crack tip. Each atom in Region I is freeto move, thus, in Mode I, there are 2 N degrees of freedomassociated with Region I and, in Mode III, there are N degreesof freedom.Atoms in Region II, known as the interface , all contain atleast one atom from Region I in its interaction range and arelabelled from i = N + i = N . Region III is the far field .The total energy is divided into two parts E = E ( ) ( { x ( i ) } N i = ) + E ( ) ( { x ( i ) } i > N ) . (3)The far-field part represented by E ( ) is known to be un-bounded for a body containing a crack opening. In practicalapplications the quantity of interest is thus the energy differ-ence between a suitably chosen initial configuration { x ( i ) } and a relaxed configuration { x ( i ) } , which is denoted by E − E = E ( { x ( i ) } ) − E ( { x ( i ) } ) . There are different ways of specifying the behaviour ofatoms in Region II and III. In what follows we first review two approaches, namely a simple static boundary conditionspecified by continuum linearised elasticity and a simplifiedversion of the flexible boundary scheme due to Sinclair [3].Subsequently, we show how the flexible boundary schemeleads to a simple equation to check for admissible values ofthe stress intensity factor for which equilibria exist, which mo-tivates defining an alternative version of the flexible boundaryscheme with improved accuracy.This is then followed by a discussion about applying nu-merical continuation techniques to both formulations and theresulting bifurcation diagrams capturing crack propagation,energy barriers and the phenomenon of lattice trapping [6].
B. Static boundary scheme
Prescribing a simple static far-field boundary conditionconsists of constraining atoms in both Region II and RegionIII to be displaced according to continuum linear elasticity(CLE) equations arising from the mode of crack considered,the Cauchy-Born relation and the interatomic potential em-ployed. Crucially, these equations are derived for the cracktip fixed at the origin, i.e. with α = ˆ x ( i ) = r i ( cos θ i , θ i ) employed, the anti-plane CLE displacement for Mode IIIcrack is given by U III
CLE ( i ) = √ r i ( , , sin ( θ i / ) , (4)whereas the isotropic in-plane CLE displacement for Mode Icrack is given by U I CLE ( i ) = √ r i (cid:16) ( θ i / ) − cos ( θ i / ) , (5)5 sin ( θ i / ) − sin ( θ i / ) , (cid:17) . The anisotropic Mode I case will be considered in Sec-tion III B, where we also allow out-of-plane relaxation as wellas in-plane. The superscript on U CLE distinguishing differentcrack modes is dropped whenever a distinction is not needed.The displacement fields we consider are in fact of the form { K U
CLE ( i ) } , where K ∈ R is the stress intensity factor.A suitable way of encoding this far-field behaviour is toconsider configurations { x ( i ) } in the form x ( i ) = ˆ x ( i ) + K U
CLE ( i ) + U ( i ) , (6)where, as in (1), ˆ x ( i ) is the crystalline lattice position and U ( i ) the atomistic correction of the i th atom, accounting forthe fact that atoms within Region I are free to relax under theinteratomic potential. This correction is constrained to satisfy U ( i ) = for i > N , which ensures that atoms outside the coreremain fixed at the CLE displacement field.In this framework the initial configuration { x ( i ) } againstwhich the energy difference is computed corresponds to set-ting U ( i ) = i . As a result, it trivially holds that E ( { x ( i ) } i > N ) = E ( { x ( i ) } i > N ) , and hence E − E = E ( { x ( i ) } N i = ) − E ( { x ( i ) } N i = ) , which is a finite quantity.If one defines a function of ( { U ( i ) } , K ) given by F = ( { f ( i ) } N i = ) , (7)where f ( i ) = − ∂ E ∂ x ( i ) is the force acting on the i th atom, thenan equilibrium configuration can be found by solving F = .For a fixed K , in Mode III, this is a system of N equationsfor N variables, and, in Mode I there are 2 N equations for2 N variables. In a fully 3D case, to be considered in numer-ical tests in Section III B, the system considered consists of3 N equations for 3 N variables.The remaining difficulty is the interplay between the choiceof K and the crack tip position. This will be addressed in Sec-tion II D with the help of numerical continuation, in particu-lar highlighting fundamental limitations of the static boundaryscheme. C. Flexible boundary scheme
1. Standard formulation
The central idea of the flexible boundary scheme describedby Sinclair [3] is to allow the crack tip position ( α , ) to vary.The displacements can be shifted to account for the currentcrack tip position by redefining the polar coordinates used in(4) and (5) so that r i ( cos θ i , sin θ i ) = ˆ x ( i ) − ( α , ) . The configurations considered are, similarly to (6), of the form x ( i ) = ˆ x ( i ) + K U α CLE ( i ) + U ( i ) , where the CLE displacement is now written as U α CLE to em-phasise the dependence on α throught the shift of the polarcoordinate system. The chosen initial unrelaxed configuration { x ( i ) } is obtained by setting α = U ( i ) = i .The effect that varying α has on the system can be capturedby considering the notion of a generalised forcef ∞ α = − ∂ E ∂ α = ∑ i − ∂ E ∂ x ( i ) · ∂ x ( i ) ∂ α = ∑ i f ( i ) · K V α ( i ) , (8)where V α = − ∂ U α CLE . As stated, this is an infinite sum, whichcan be shown to be convergent since U CLE solves the CLEequation [21].Somewhat arbitrarily, Sinclair assumes that in Region IIIthe crystal is fully ’linearly elastic’ [3], in the sense that the continuum CLE displacement is an equilibrium by itself,meaning that f ( i ) = , for i > N , (9)for any choice of K and α , effectively truncating the infinitesum in (8).The Sinclair scheme can be formalised by defining a func-tion of ( { U ( i ) } , K , α } given by F = ( { f ( i ) } N i = , f α ) , (10)where f α = N ∑ i = f ( i ) · K V α ( i ) . (11)An equilibrium configuration in this scheme is then obtainedby solving F = .Notably, the summation in (11) is effectively over { i } N i = N + , since at an equilibrium f ( i ) = i ≤ N .Obviously, (9) is only true in an approximate sense, whichleaves open to interpretation whether the truncation enforcedthrough (9) is the optimal choice.It is further worth noting that in the limit when N → ∞ , thegeneralised force f α in (11) is null at any equilibrium, hencethe extra equation f α = E − E in the newscheme, we follow the procedure described in Ref. 3,Appendix 1, with the far-field contribution to the energy E ( ) − E ( ) from (3) approximated as E ( ) − E ( ) = − N ∑ i = N + (cid:16) f ( ) ( i ) + f ( ) ( i ) (cid:17) · ( x ( i ) − x ( i )) , where f ( ) ( i ) = − ∂ E ( ) ∂ x ( i ) .Making sense of the arguably unsubstantiated far-field ap-proximation in (9) as well as addressing the question of con-vergence leads to several interesting realisations that will beaddressed in the next section.
2. Predicting admissible stress intensity factors
The strain fields associated with atomistic corrections { U ( i ) } are known to decay more quickly away from the corethan the strain fields associated to { U CLE ( i ) } (proven rigor-ously in a simplified setup in Ref. 22), meaning that theircontributions are effectively negligible beyond a small regionaround the crack tip. It can thus be conjectured that a reason-able approximation to the flexible boundary scheme condition f α =
0, defined in (11), is to allow the unrelaxed configuration { x ( i ) } to depend on α and look at the generalised force atthe unrelaxed configuration, namely − ∂ E ∂ α = ∑ i − ∂ E ∂ x ( i ) · ∂ x ( i ) ∂ α = ∑ i f ( i ) · K V α ( i ) . Employing the same truncation as in (9), we can postulate acondition f α =
0, where f α = N ∑ i = f ( i ) · K V α ( i ) . (12)With the unrelaxed configuration { x ( i ) } , now determinedsolely by K and α , verifying whether f α = f α = K values nearly perfectlyoscillating around a fixed interval of admissible values K − < K < K + .With the numerical tests indicating that the predicted in-terval is strongly dependent on the size of the computationaldomain, it seems plausible that changing the far-field trunca-tion rule from (9) can have a drastic effect on the computedsolution path. This will be investigated in the next section.
3. Effect of changing far-field truncation rule
The truncation in (9) is equivalent to stating that the atom-istic information associated with atoms in Region III, whichconceptually is an infinite far field, is completely disregarded.One can provide the flexible boundary scheme with moreatomistic input from Region III by changing the truncationin (9) to f ( i ) = , for i > N , (13)where N > N is much larger, e.g. N = N . As a result, anew condition is ˜ f α =
0, where˜ f α = N ∑ i = f ( i ) · K V α ( i ) . Thus, the new approach testing the effect of different trunca-tion can be formalised by defining a function of ( { U ( i ) } , K , α ) given by F = ( { f ( i ) } N i = , ˜ f α ) , (14)with an equilibrium configuration obtained by solving F = F = F = f α from (11) satisfies f α = N ∑ i = N + f ( i ) · K V α ( i ) . (15) The right-hand side admits input only from atoms in Re-gion III, whose displacements are determined solely by α and K , since by design, U ( i ) = i ≥ N , which highlights thegeneral rationale behind this formulation: the far-field regionwithin the computational domain is vastly enlarged, but onlytwo degrees of freedom remain attached to it, meaning thatin practice there is virtually no additional computational cost,apart from the ability to compute the right-hand side of (15).It will be shown through numerical tests presented in Sec-tion III A 4 that the new scheme results in much improvedaccuracy for small sizes of the core region, implying that inpractice the new scheme is numerically preferable, enablingincreased accuracy to be achieved with decreased numericalcost. D. Pseudo-arclength numerical continuation
The basic premise of numerical continuation applied to theproblem at hand is as follows. Suppose we have identifiedsome K for which some equilibrium configuration x exists.Can we use this knowledge to quickly find another equilib-rium for K + δ K , for some small δ K ? Such an approachwill work well if there exists a continuous path of solutions K (cid:55)→ { U ( i ) K } N i = (and K (cid:55)→ α K in the case of the flexibleboundary scheme). Such a path can be shown to exist, cour-tesy of Implicit Function Theorem [23], in the neighbourhoodof K if the associated Hessian operator is invertible at K .A more sophisticated version, which is particularly usefulfor the problem at hand, is known as the pseudo-arclengthcontinuation. It postulates that the quantities involved all aresmooth functions of an arclength parameter s . The questionthus changes to: given some triplet ( K s , { U s ( i ) } , α s ) (in thecase of static boundary α ≡ { x s ( i ) } , can we find a new tripletfor s : = s + δ s , for some small δ s , which gives us a newequilibrium { x s ( i ) } ? The key advantage of this approach isthat it can handle index-1 saddle points, which makes it a use-ful tool for studying energy barriers and the phenomenon oflattice trapping.Numerical continuation can be incorporated into the frame-work by including K as a variable in the systems of equations F j = j = , , K as a variable renders each system ofequations F j = f K = f K = N ∑ i = ( U s ( i ) − U s ( i )) · ˙ U s ( i ) + ( α s − α s ) ˙ α s (16) + ( K s − K s ) ˙ K s − δ s . Here (cid:0) ˙ K s , { ˙ U s ( i ) } , ˙ α s (cid:1) refer to derivatives with respect to s evaluated at s . Note that in the static boundary scheme wesimply have α s i =
0, and hence ˙ α s i =
0, for all i .If there indeed exist a continuous path of solutions, it can be X s n X s n +1 X s n s FIG. 1. Suppose X s n = ( K s n , { U s n ( i ) } , α s n ) belongs to the solu-tion path, which is smooth (solid green line). The schematic two-dimensional plot indicates that imposing an extra constraint in theform of f K = X s n ,represented by a purple arrow, with δ s determining how far along ˙ X s n we choose to travel. Choosing δ s small enough ensures that the al-gorithm can safely traverse the solution path even along folds. It isfurther clear from the plot that, for δ s sufficiently small, X s n + δ s ˙ X s n (orange dot) is a very good initial guess for subsequent steps of theNewton iteration. shown [24] that for a step size δ s small enough, a simple New-ton iteration will converge to a new solution of ( F j , f K ) = .The remaining difficulty is to compute (cid:0) ˙ K s , { ˙ U s ( i ) } , ˙ α s (cid:1) .This can be achieved by first noting that with s being an ar-clength parameter, by definition it has to hold that N ∑ i = ˙ U s ( i ) · ˙ U s ( i ) + ( ˙ K s ) + ( ˙ α s ) = . (17)This eliminates one degree of freedom. The remaining de-grees of freedom can be eliminated by differentiating bothsides of F j = s , which ispossible under the assumption of there existing a smooth pathof solutions. Details are presented in the Appendix.The resulting pseudo-arclength continuation algorithms as-sociated with both schemes are presented as Algorithm 1 andAlgorithm 2. Algorithm 1
Static boundary condition pseudo-arclength con-tinuation Given δ s ; given a stable equilibrium configuration solving ( F , f K ) = determined by ( K s , { U s ( i ) } ) ; compute ( ˙ K s , { ˙ U s ( i ) } ) using (A.3); compute a new stable equilibrium configuration ( K s , { U s ( i ) } ) by solving ( F , f K ) = using Newton iteration with initialguess ( K s + δ s ˙ K s , { U s ( i ) + δ s ˙ U s ( i ) } ) ; for n > do given ( K s n , { U s n ( i ) } ) and ( ˙ K s n − , { ˙ U s n − ( i ) } ) ; compute ( ˙ K s n , { ˙ U s n ( i ) } ) by solving linear system (A.6); compute a new equilibrium configuration ( K s n + , { U s n + ( i ) } ) by solving ( F , f K ) = using Newton iteration with initialguess ( K s n + δ s ˙ K s n , { U s n ( i ) + δ s ˙ U s n ( i ) } ) . end for Algorithm 2
Flexible boundary condition pseudo-arclengthcontinuation Given δ s ; given a stable equilibrium configuration solving ( F , f K ) = determined by ( K s , { U s ( i ) } , α s ) ; compute ( ˙ K s , { ˙ U s ( i ) } , ˙ α s ) using (A.5); compute a new stable equilibrium configuration ( K s , { U s ( i ) } , α s ) by solving ( F , f K ) = using Newton iter-ation with initial guess ( K s + δ s ˙ K s , { U s ( i ) + δ s ˙ U s ( i ) } , α s + δ s ˙ α s ) ; for n > do given ( K s n , { U s n ( i ) } , α s n ) and ( ˙ K s n − , { ˙ U s n − ( i ) } , ˙ α s n − ) ; compute ( ˙ K s n , { ˙ U s n ( i ) } , ˙ α s n ) by solving linear system (A.7); compute a new equilibrium configuration ( K s n + , { U s n + ( i ) } , α s n + ) by solving ( F , f K ) = using Newton iteration with initial guess ( K s n + δ s ˙ K s n , { U s n ( i ) + δ s ˙ U s n ( i ) } , α s n + δ s ˙ α s n ) . end for Bearing in mind that most realistic interatomic potentialsonly provide analytic forces but not Hessians, meaning thatAlgorithm 2 cannot be readily used, as it requires a computa-tion of the Hessian while differentiating F j = ( ˙ K s n , ˙ U s n , ˙ α s n ) , we also propose a simple finite-differencebased approximate scheme as a Hessian-free alternative.The method consists of first computing two stable equi-librium configurations determined by ( K s , { U s ( i ) } , α s ) and ( K s , { U s ( i ) } , α s ) , which crucially satisfy K s ≈ K s (andalso α s ≈ α s ). In the first step the tangent ( ˙ K s , { ˙ U s ( i ) } , ˙ α s ) can be approximated as˙ U s ( i ) = K s − K s (cid:0) U s ( i ) − U s ( i ) (cid:1) , (18a)˙ α s = K s − K s (cid:0) α s − α s (cid:1) , (18b)˙ K s = , (18c)with the last line a direct consequence of K being the effectivecontinuation parameter in the first step, since it is K that isvaried to obtain two stable equilibrium configurations.With the tangents computed, one can now assemble the ex-tended system and solve ( F , f K ) = to obtain an equilibriumdetermined by ( K s , { U s ( i ) } , α s ) . The switch to the extendedsystem entails that now the arclength s is the continuation pa-rameter and in particular s − s = δ s , with δ s fixed through-out. As a result, subsequent tangent approximations are com-puted, for n = , . . . , as˙ U s n ( i ) = δ s (cid:0) U s n ( i ) − U s n − ( i ) (cid:1) , (19a)˙ α s n = δ s (cid:0) α s n − α s n − (cid:1) , (19b)˙ K s n = δ s (cid:0) K s n − K s n − (cid:1) . (19c)The details of this approximate scheme are summarised asAlgorithm 3 below. In practical applications, to avoid pos-sible numerical artefacts, the finite-difference approach couldbe substituted by the automatic differentiation approach [25]. Algorithm 3
Hessian-free approximate flexible boundarycondition pseudo-arclength continuation Given δ s ; given two stable equilibrium configurations solving ( F , f K ) = , determined by ( K s , { U s ( i ) } , α s ) and ( K s , { U s ( i ) } , α s ) ,and satisfying K s ≈ K s and α s ≈ α s ; compute an approximate ( ˙ K s , { ˙ U s ( i ) } , ˙ α s ) using (18); compute a new stable equilibrium configuration ( K s , { U s ( i ) } , α s ) by solving ( F , f K ) = using Newton iter-ation with initial guess ( K s + δ s ˙ K s , { U s ( i ) + δ s ˙ U s ( i ) } , α s + δ s ˙ α s ) ; for n > do given ( K s n , { U s n ( i ) } , α s n ) and ( K s n − , { U s n − ( i ) } , α s n − ) ; compute an approximate ( ˙ K s n , { ˙ U s n ( i ) } , ˙ α s n ) using (19); compute a new equilibrium configuration ( K s n + , { U s n + ( i ) } , α s n + ) by solving ( F , f K ) = using Newton iteration with initial guess ( K s n + δ s ˙ K s n , { U s n ( i ) + δ s ˙ U s n ( i ) } , α s n + δ s ˙ α s n ) . end for III. RESULTS In this section we discuss numerical tests based around ap-plying the pseudo-arclength continuation to both the static andflexible boundary schemes.We begin by directly comparing the static boundary schemeand the flexible boundary scheme when applied to a simple toymodel, highlighting the superiority of the latter. This is thenfollowed by a study of fracture on the ( ) cleavage plane insilicon with two interatomic potentials. A. Mode III toy model
We first consider a toy model of anti-plane Mode III frac-ture posed on a triangular lattice with lattice constant equal tounity and atoms interacting according to a nearest neighbourpair potential. The total energy is thus of the form E = ∑ i (cid:54) = j | ˆ x ( i ) − ˆ x ( j ) | = φ ( r i j ) , where r i j = | x ( i ) − x ( j ) | , Quantity Value a µ γ K G a , shear modulus µ , sur-face energy γ and Griffith stress intensity factor K G computed for theMode III toy model. where φ ( r ) = (cid:0) − exp ( − r ) (cid:1) . The resulting material properties are reported in Table I, in-cluding the shear modulus, the surface energy and the Griffithprediction for the critical stress intensity factor K G .To investigate domain size effects, we consider computa-tional domains of different sizes, each geometrically repre-sented by ball of radius R around the origin. The three choiceof radii are (1) R =
32, (2) R =
64 and (3) R = | ˆ x ( i ) | < R − R out − R φ , (20)where R out = . R φ = . N = N = N =
1. Pseudo-arclength continuation with static boundary scheme
Algorithm 1 is first employed to compute solution pathspresented in Figure 2. With no knowledge of the actual cracktip position, the y -axis was chosen to represent the Euclideannorm of { U s ( i ) } .The plot confirms the intuitively clear notion that | U s | willbe smallest when there is no mismatch between the predictedcrack tip position (in the static boundary scheme fixed at α =
0) and the actual crack tip position. Periodic wiggles furtherindicate a repeating bond-breaking behaviour.The solution paths are heavily tilted, with no clear range ofstress intensity factors for which equilibria exist, as K growsto effectively compensate for α being fixed. In particular, nounstable equilibria are found and the energy is monotonicallyincreasing in K , implying that no study of energy barriers ispossible.An ad-hoc post-processing way of estimating actual valuesof α and K is to findmin α , K (cid:32) N ∗ ∑ i = | K s U CLE ( i ) + U s ( i ) − KU α CLE ( i ) | (cid:33) , (21)where, in order to avoid boundary effects, { i } N ∗ i = correspondsto all atoms such that | ˆ x ( i ) | < R . The resulting plots of α against K are shown with dashed lines in Figure 3. Stress intensity factor K / K G E u c li d e a n n o r m o f t h e c o rr e c t o r U R = 32 R = 64 R = 128 FIG. 2. Solution paths obtained for Mode III toy model using Algo-rithm 1 for three choices of domain size.
2. Pseudo-arclength continuation with flexible boundary scheme
Algorithm 2 is now employed to compute solution paths ofthe toy model for three different domain sizes, as described inSection III A.A direct comparison of both scheme is shown in Figure3, revealing that the flexible scheme is superior to the post-processing of the static scheme in terms of predicting therange of the stress intensity factors for which equilibria exist.In particular, the flexible scheme employed on a core regionwith radius R =
32 is as accurate as the post-processed staticscheme employed on a core region with radius R = Stress intensity factor K / K G H o r i z o n t a l c r a c k t i p p o s i t i o n R = 32 flexible R = 64 flexible R = 64 static R = 128 flexible R = 128 static FIG. 3. Comparison of solutions paths obtained for Mode III toymodel using Algorithm 2 (solid lines) and Algorithm 1 (dotted lines,post-processed via (21)) for three domain sizes (static R =
32 casetoo far away to the right to include).
With unstable equilibrium configurations correspondingto index-1 saddle points captured in the flexible boundaryscheme, a study of energy barriers is now feasible, as show-cased in Figure 4 and later in Figure 6. H o r i z o n t a l c r a c k t i p p o s i t i o n Stress intensity factor K / K G R e s c a l e d e n e r g y g a i n FIG. 4. A study of energy barriers in the Mode III toy model forcomputational domain of radius R = ( E − E ∗ ) / E ∗ , where E ∗ is the energy of the bottom left configu-ration. The dashed parts of the solution path denote index-1 saddlepoints, which correspond to energetic cost of crack propagation at agiven value of K , which can be seen by observing in the lower plotthat dashed lines lie above their neighbouring solid lines. A nudgedelastic band calculation further confirming this being the case is pre-sented in Figure 6. The point where stable parts of the solution pathscross corresponds to the critical stress intensity factor K c , notablynot quite matching the Griffith stress intensity factor K G . This phe-nomenon is elaborated upon in Section III A 5.
3. Predicting the admissible range for K
The ideas developed in Section II C 2 are now checked nu-merically for the toy model presented in Section III A, againemploying three domain sizes. The results are presented inFigure 5.The prediction of the range of admissible values of thestress intensity factor based on the CLE displacements only isshown to be fairly accurate, with the magnitude for K match-ing, while the predicted length of the interval considerablylarger than in reality. Importantly, the prediction correctlyshifts with the changing domain size, indicating that the rangeof admissible values for K is to a considerable extent deter-mined by the far-field behaviour only, thus strongly motivat-ing the new formulation of the flexible scheme presented inSection II C 3, which will be tested numerically in the nextsection. Stress intensity factor K / K G H o r i z o n t a l c r a c k t i p p o s i t i o n R = 32 full R = 32 approx R = 64 full R = 64 approx R = 128 full R = 128 approx FIG. 5. Solutions paths computed for the Mode III toy model withAlgorithm 2 (solid lines) for the three domain sizes, plotted againsta corresponding path of approximate solutions obtained by solving f α =
4. Pseudo-arclength continuation with flexible boundary schemewith extended far-field region
To test the effect of extending the far-field region discussedin Section II C 3, we consider a computational domain in theform of a ball of radius R =
128 with varying sizes of Region I.As before, the core region is chosen to consists of all atomssatisfying (20), this time with (1) R =
8, (2) R =
16, (3) R = R =
64. Region II is again an annulus of width R φ = . R − R + R out , as opposed to just R out = . N = N =
292 (2) N = N = N = K against α arepresented in the middle panel of Figure 7, which also includethe solution path computed with the standard flexible schemewith R =
256 for comparison.The extension of the far-field region drastically increasesthe accuracy of the flexible boundary scheme, with a tiny fullyatomistic region required to have a very accurate predictionfor the admissible range of values for the stress intensity fac-tor. This is demonstrated quantitatively in the error analysis inSection III A 5.Despite the large far-field region, the system of nonlinearequations associated with the new scheme when R = R = R =
32, with details presented in Figure 6.
Stress intensity factor K / K G H o r i z o n t a l c r a c k t i p p o s i t i o n Crack advance R e s c a l e d e n e r g y g a i n FIG. 6. A study of energy barriers in the flexible boundary schemewith extended far-field with R =
32 and R = K are chosen and in each casean initial minimum energy path (MEP) is formed by linear interpo-lation between the first stable equilibrium, the saddle in between andthe second stable equilibrium. The path is then optimised using thenudged elastic band method [26]. The resulting MEPs are shown inthe bottom figure, together with larger dots representing the equilib-ria computed with pseudo-arclength continuation, thus confirmingthat the middle equilibrium is indeed a saddle and also confirminglack of other critical points along a given path.
5. Error analysis
To conclude the numerical investigation of the toy model,a brief error analysis is presented in Figure 7. The referencesolution path, imitating the infinite limit N → ∞ is obtainedwith the standard flexible boundary scheme, as described inSection III A, with R = R = , , , R =
128 and R = , , , O ( R − ) , whichimproves upon a known rate of convergence O ( R − / ) of thestatic scheme proven in Ref. 15. A mathematically rigorousproof of the improved rate of convergence will be a subjectof further study. Notably, the error analysis together with thestudy of energy barriers presented in Figure 4 and 6 clearlyshow that the Griffith prediction for the critical stress inten-sity factor K G is only valid in the limit N → ∞ .Secondly, the extended far field flexible boundary schemeremains as accurate as the outer radius, which in the currentstudy is fixed at R = R , confirming the intuitionbehind this reformulation of the flexible boundary scheme.The underlying reasons for this are also to be explored in afuture work. B. Mode I fracture of silicon on the ( ) cleavage plane We next test our new algorithms on a more complex prob-lem: fracture of silicon on the ( ) cleavage plane in the [ ] propagation direction. This is known to be the preferredlow-energy cleavage orientation, but the precise details of thelattice trapping barriers to brittle fracture remain elusive forthe reasons outlined in the introduction, making this a prob-lem of scientific interest as well as an interesting test case.We consider two interatomic potentials known to give aqualitatively correct description of brittle fracture for this sys-tem: modifications of the Tersoff [28] and Kumagai [29] po-tentials, with the interaction length increased and additionalscreening functions introduced to improve the description ofbond-breaking processes [19]. Without these modifications,neither potential predicts brittle behaviour. The modified po-tentials have been shown to predict lattice trapping ranges K − < K < K + for the ( ) cleavage plane in reasonableagreement with DFT, albeit restricted to a small model systemwith static boundary conditions [19], and we thus use themhere as a proxy for a fully description of interatomic bondingin silicon. The potentials have not previously been applied tostudy fracture on the ( ) plane, in part because of the com-plexities introduced by the played by surface reconstructionssuch as the Pandey 2 × π -bonded chain [12, 30], which wedo not study here.A number of small modifications to the FBC method de-scribed above are needed. Since analytical Hessians arenot readily available for these potentials, we use the finite-difference reformulation of the scheme outlined in Algorithm3. For comparison with the static case, we also consider aHessian-free version of Algorithm 1, which can be obtainedfrom Algorithm 3 by fixing α = α = U α CLE and its derivative V α are alsoredefined to account for the anisotropy of the silicon crys- Quantity Tersoff+S Kumagai+S r c / Å 6.0 6.0 a / Å 5.432 5.429 C / GPa 143 165 C / GPa 75 65 C / GPa 69 77 γ ( ) / Jm − K G / MPa √ m 1.07 0.97TABLE II. Values of the cutoff radius r c , lattice constant a , cubicelastic constants C , C , C , surface energy γ ( ) and Griffithstress intensity factor K G computed with the screened versions ofthe Tersoff and Kumagai interatomic potentials. tal using the near field solution for a crack in a rectilinearanisotropic elastic medium (noting that V α can convenientlybe obtained from the xx and xy elements of the deformationtensor) [31]. The CLE solutions are expanded from two tothree dimensions using plane strain loading conditions appro-priate for a simulation cell periodic along the crack front line,i.e. Y =
0, with the atomistic corrector U ( i ) for each atomalso becoming three dimensional. In place of the Newton iter-ation, we solve ( F , f K ) = LGMRES package [32].For large systems, it is necessary to precondition the solver.We used a general purpose preconditioner for materials sys-tems [33], augmented by a diagonal rescaling of the f α and f K components of the preconditioner to balance their mag-nitudes with that of the atomic forces f ( i ) (as suggested bySinclair [3]). Finally, the crack tip force f α is now computedby summing only over atoms in region II (or regions II andIII for the extended far-field variant); as discussed after (11)this does not affect the equilibria obtained. A software imple-mentation of the algorithm is available with the framework ofthe Atomic Simulation Environment (ASE) [34] as part of theopen source matscipy package [35].To setup the simulations, the lattice and elastic constantsand the surface energy of the ( ) plane are computed foreach potential and are reported in Table II, along with theGriffith prediction for the critical stress intensity factor K G ,obtained using the relaxed surface energy γ ( ) .Similar to the toy model above, we consider three domainradii (1) R =
32 Å, (2) R =
64 Å and (3) R =
128 Å, with theradius of the fully atomistic region I chosen to consider atomswith crystal positions | ˆ x ( i ) | < R − R out − R φ where now we take R out = r c =
12 Å as the width of theannulus of atoms defining Region III and R φ = r c = . r c is added to ensure the forces on atoms inRegion III are unaffected by the presence of the outer surface.The corresponding numbers of atoms in region I are (1) N = N = N = N + N + Stress intensity factor K / K G H o r i z o n t a l c r a c k t i p p o s i t i o n Flexible scheme R = 8 R = 16 R = 32 R = 64 R = 256
Stress intensity factor K / K G Extended far field flexible scheme R = 8 ( R = 128) R = 16 ( R = 128) R = 32 ( R = 128) R = 64 ( R = 128) R = 256
10 20 30 50
Region I radius R D i s t a n c e t o t h e r e f e r e n c e s o l u t i o n p a t h R Convergence of solutions paths flexflex ext
FIG. 7. Solution paths and error analysis for flexible and extended far field flexible boundary schemes. The solution path computed with theflexible boundary scheme with R =
256 serving as a reference in the error analysis. C r a c k a d v a n c e / Å R = 32 Å static R = 64 Å static R = 128 Å static0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 Stress intensity factor K / K G C r a c k a d v a n c e / Å R = 32 Å flexible R = 64 Å flexible R = 128 Å flexible FIG. 8. Comparison of solution paths obtained with Algorithm 3 fora Si ( )[ ] crack modelled with the screened Kumagai potentialusing flexible (solid lines) and static (dashed lines) boundary condi-tions, for three choices of domain size. For the static cases α and K are obtained by a post-processing fit to the CLE solution.
1. Pseudo-arclength continuation with the Static and FlexibleBoundary Conditions
We first perform arc-length continuation calculations withthe Kumagai potential for three choices of domain radii, usingboth static and flexible boundary conditions. The results areshown in Figure 8. For the static cases, we employ the post-processing fit for K and α given in (21), leading to the resultsshown with dashed lines in the figure. For small domain sizes, the static solutions are highly tilted, while the flexible solu-tions show the correct periodic behaviour even at the smallestdomain size. We note that the absolute value of α is some-what arbitrary, since shifts of the solution path along the crackpropagation direction can be incorporated in different choicesof the corrector field U .The high accuracy of the flexible scheme allows a care-ful comparison of the lattice trapping predicted by differentchoices of interatomic potential to be made, as shown in Fig-ure 9. At a domain size of 128 Å the solution paths are al-ready very close to periodic in the crack propagation direc-tion. The energy differences computed with (3) illustrated inthe lower panels confirm that there is a critical stress intensityfactor K c for which the total energy of the atomistic plus con-tinuum system is equal at all stable energy minima, i.e. beforeand after crack advance. While the range of lattice trapping K − < K < K + predicted by the two potentials is similar, forboth potentials K c is less than the Griffith equilibrium value K G . The values of K G used here were computed from the re-laxed ( ) surface energy, indicating, as well as remainingfinite size effects, some of the difference could be attributedto local modifications of the surface energy close the cracktip — a discrepancy that could be further exacerbated by thepresence of more complex surface features such as the Pandey2 × K = . K G .The inset schematics illustrate how this feature arises: movingalong the stable path from (A) to (B), the bond at the crack tipremains intact as the centre of the continuum field α advances.The bond gradually opens as we move towards point (C) in theunstable part of the solution path, while between (C) and (D)it opens more rapidly as the atoms ‘snap’ apart. We postu-late that this sharp feature is associated with the finite cutoffof the potential, which, despite the screening terms that makefracture simulations feasible, is still a modelling assumption.In future work we aim to compute solution paths with DFTto remove the uncertainty associated with the use of simpli-1 C r a c k a d v a n c e / Å Kumagai+S, R = 128 Å A BCD
Tersoff+S, R = 128 Å Stress intensity factor K / K G E n e r g y d i ff e r e n c e EE / e V Stress intensity factor K / K G A BCD (A) = 5.7 K / K G = 0.85(B) = 7.1 K / K G = 1.17(C) = 7.7 K / K G = 0.98(D) = 8.3 K / K G = 0.95 FIG. 9. Comparison of lattice trapping of Si ( )[ ] cracks predicted by the screened Kumagai (left panels) and Tersoff (middle andright panels) potentials. Upper panels: solution paths obtained with Algorithm 3, including stable parts (solid lines) corresponding to energyminima and unstable parts (dashed lines) corresponding to saddle points. Lower panels: energy difference with respect to the CLE solutionwith α = , K = K G . Right insets: near-tip atomic positions corresponding to marked points A, B, C, D on the Tersoff solution and energypaths, with the opening bond highlighted in red. fied potentials: this remains out of reach for the present since,despite the considerable improvements in accuracy affordedby the flexible scheme, converged solution paths still requirea large number of force evaluations on systems comprisingseveral thousand atoms.
2. Predicting the admissible range of K
The admissible range of K is now predicted by finding rootsof equation (12), i.e. where f α ( K , α ) =
0, leading to the pre-dictions shown with the dashed lines in Figure 10. Here, K isfound numerically for each value of α in a 200-element grid.For both potentials, the admissible range of K is in reason-able agreement with that computed in the full solution paths,suggesting that our approach provides a useful way to esti-mate the stable range of K for the cost of a fixed number of force evaluations on the full domain.
3. Pseudo-arclength continuation with an extended far-field region
To conclude the numerical tests, we also apply the extendedfar-field scheme of (14) to the Si ( )[ ] crack system,modelled using the screened Tersoff potential. The overalldomain size is fixed at ¯ R =
128 Å, and two choices of radiifor Region I are considered: R I =
14 Å and R I =
46 Å, cho-sen since these lead to problems with the same numbers ofdegrees of freedom as the R =
32 Å and R =
64 Å flexiblemodels considered earlier.The results are illustrated in Figure 11. Although there isan improvement over the standard flexible scheme in conver-gence towards the reference R =
128 Å solution path, par-ticularly for the smallest Region I size, these results do not2
Stress intensity factor K / K G C r a c k a d v a n c e / Å Kumagai+S, R = 128 Å Full solutionApprox solution 0.9 1.0 1.1 1.2
Stress intensity factor K / K G Tersoff+S, R = 128 Å Full solutionApprox solution
FIG. 10. Comparison of full solution paths obtained by arc-lengthcontinuation with Algorithm 3 (solid lines) and corresponding ap-proximate solution paths for f α = ( )[ ] crack modelled with the screened Kumagai and Tersoff potentials,using a domain radius of 128 Å. C r a c k a d v a n c e / Å R I = 14 Å ( R = 32 Å) flex R I = 46 Å ( R = 64 Å) flex R I = 110 Å ( R = 128 Å) flex0.8 0.9 1.0 1.1 1.2 1.3 1.4 Stress intensity factor K / K G C r a c k a d v a n c e / Å R I = 110 Å ( R = 128 Å) flex R I = 14 Å ( R = 128 Å) flex ext R I = 46 Å ( R = 128 Å) flex ext FIG. 11. Comparison of original (denoted ‘flex’, upper panel, solidlines) and extended far field (denoted ‘flex ext’, lower panel, dashedlines) variants of the flexible boundary condition approaches topseudo-arclength continuation in an Si ( )[ ] crack system mod-elled with the screened Tersoff potential. The R =
128 Å ‘flex’ resultis shown in both panels to allow comparison. Calculations with theextended scheme show improved accuracy at smaller domain radii. provide convincing evidence that a larger far-field region sig-nificantly enhances the accuracy of the scheme. This in incontrast to the results obtained with the toy model, suggestingthat an enhanced CLE predictor that improves the match withthe atomistic model is needed to further increase accuracy.
IV. CONCLUSIONS
In this work we have reported an extension of Sinclair’sflexible boundary condition algorithm to allow full solutionpaths for cracks to be computed using pseudo-arclength con-tinuation. We have also introduced an extension of the FBCalgorithm which allows information to be incorporated froma larger far-field region, and which also provides a steppingstone towards putting the method on a more rigorous math-ematical footing. We demonstrated the approach for ModeIII fracture with a 2D toy model, and for Mode I fracture ofsilicon using realistic interatomic potentials that give a quali-tatively correct description of fracture.In future, our approach will enable a detailed study of lat-tice trapping barriers to brittle fracture to be carried out us-ing increasingly realistic models of interatomic bonding, go-ing beyond the screened bond-order potentials demonstratedhere, for example by using machine-learning interatomic po-tentials [36] or DFT directly. This could help to resolve ques-tions such as the role of blunt-sharp-blunt crack tip recon-struction observed during fracture in the Si ( )[ ] cracksystem [37], where NEB calculations demonstrated the crackis bluntened at stable minima and sharp at the unstable tran-sition states. Moreover, the new approach could be ex-panded to study crack path selection, known to exhibit com-plex phenomenon in anisotropic materials [38], or the dy-namics of three dimensional crack fronts, going beyond pre-vious work that was limited to simple interatomic potentialsand small model systems [10]. Truly accurate predictions ofcritical stress intensity factors and lattice trapping ranges re-quire a quantum mechanical approach, at least near the cracktip. Hybrid schemes such as QM/MM (quantum mechan-ics/molecular mechanics), previously applied to dynamic frac-ture [12], could be combined with the algorithms introducedhere to make quantitative fracture toughness calculations ac-curate and affordable. Established routes could then be fol-lowed to produce atomistically informed continuum mod-els [39–41].Before this can be done, however, further work is needed toassess finite-size effects. For the silicon fracture application,we have demonstrated that the flexible scheme is superior tostatic boundaries, but not yet quantified the convergence rate,meaning that the new algorithms cannot yet be used for pre-dictive materials science. Ultimately, it is hoped that the flex-ible boundary scheme and numerical continuation techniquescan be combined with higher-order far-field predictions to in-crease accuracy in a quantifiable manner.Finally, we note that the pseudo-arclength continuationused here would also be applicable to other defects such asdislocations by replacing the stress intensity factor K as a bi-furcation parameter with the applied shear stress, which alsoenters as a prefactor in front of the CLE solution. ACKNOWLEDGMENTS
We thank Christoph Ortner and Lars Pastewka for use-ful discussions. We acknowledge funding from the EP-3SRC under grant numbers EP/R012474/1, EP/R043612/1 andEP/S028870/1. Additional support was provided by the Lev-erhulme Trust under grant RPG-2017-191 and the Royal Soci-ety under grant number RG160691. The authors would like toacknowledge the University of Warwick Scientific ComputingResearch Technology Platform for assistance in the researchdescribed in this paper.
Appendix: Computation of tangents in the pseudo-arclengthcontinuation scheme
In the static boundary scheme given by F =
0, differenti-ating both sides with respect s yields = H s ˙ U s + ˙ K s b Ks , (A.1)where (cid:0) H s ˙ U s (cid:1) ( i ) = N ∑ j = H s ( i , j ) · ˙ U s ( j ) . (A.2)Here H s ( i , j ) is ( i , j ) -th entry of the Hessian operator evalu-ated at s . In an infinite crystal, the Hessian evaluated at s n isan infinite block matrix with H s n ( i , j ) = ∂ E ∂ x s n ( i ) ∂ x s n ( j ) , with a short-hand notation H s n used to denote its part relatedto atoms in Region I, which is thus a N × N block matrix.The other term on the right-hand side of (A.1) is given by b Ks ( i ) = N ∑ j = H s ( i , j ) · U α s CLE ( j ) , where crucially the summation here is over both the core andthe interface regions (thus the Hessian operator here is effec-tively a rectangular block matrix of size N × N ), whereas in(A.2) the summation is only over the core region.It follows from (17) (with ˙ α ≡ U s = − ˙ K s (cid:0) H − s b Ks (cid:1) , ˙ K s = ± (cid:0) | H − s b Ks | + (cid:1) − / , (A.3)provided the square block matrix H s is invertible. The casewhen it is not invertible is known as a bifurcation point and itwill be discussed below.In the flexible boundary scheme given by F = s implies (cid:40) = H s ˙ U s + ˙ K s b Ks + ˙ α s b α s = b α s · ˙ U s + ˙ K s C α , Ks + ˙ α s C α , α s , (A.4) where b α s ( i ) = N ∑ j = H s ( i , j ) · K V α s ( j ) , C α , Ks = N ∑ i = (cid:16) b Ks ( i ) · K V α s ( i ) + f ( i ) · V α s ( i ) (cid:17) , C α , α s = N ∑ i = (cid:16) b α s ( i ) · K V α s ( i ) + f ( i ) · K V ( ) α s ( i ) (cid:17) , with V ( ) α = − ∂ V α .Note that (A.4) applies to the newly formulated scheme F = as well, except that the sums defining b Ks ( i ) , b α s ( i ) , C α , Ks and C α , α s should be over i = , . . . , N .If H s is invertible, then equations (17) and (A.4) togetherimply that˙ K s = ± (cid:18) | A | + A A + (cid:19) − / , ˙ U s = ˙ K s A , ˙ α s = ˙ K s A A , (A.5)where A = N ∑ i = − (cid:0) H − s b α s ( i ) (cid:1) · b α s ( i ) + C α , α s , A = N ∑ i = (cid:0) H − s b Ks ( i ) (cid:1) · b α s ( i ) − C α , Ks , A = − A A H − s b α s − H − s b Ks . With ( K s , { U s ( i ) } , α s ) and ( ˙ K s , { ˙ U s ( i ) } , ˙ α s ) known, astandard Newton iteration with initial guess ( K s , { U s ( i ) } , α s ) + δ s ( ˙ K s , { ˙ U s ( i ) } , ˙ α s ) , is guaranteed to converge to a new solution ( K s , { U s ( i ) } , α s ) satisfying ( F i , f K ) = provided δ s is small enough (see Figure 1 for visual insight behind this).Furthermore, the s derivative ( ˙ K s , { ˙ U s ( i ) } , ˙ α s ) can nowbe handily computed with an approximate finite-difference-like scheme, given, for F = , by = H s ˙ U s + ˙ K s b Ks , (A.6a)1 = N ∑ i = ˙ U s ( i ) · ˙ U s ( i ) + ˙ K s ˙ K s (A.6b)and, for F = (and also for F = after adjusting limits ofsummation), by = H s ˙ U s + ˙ K s b Ks + ˙ α s b α s , (A.7a)0 = b α s · ˙ U s + ˙ K s C α , Ks + ˙ α s C α , α s , (A.7b)1 = N ∑ i = ˙ U s ( i ) · ˙ U s ( i ) + ˙ K s ˙ K s , (A.7c)It is a standard assertion of bifurcation theory [42] that the lin-ear systems of equations given by (A.6) or (A.7) remain solv-able even at the points where stability change, correspondingto cases when H s is not invertible, thus allowing us to traversefull the full bifurcation diagram.4 REFERENCES [1] J. E. Sinclair and B. R. Lawn, Int. J. Fract. Mech. , 125 (1972).[2] J. E. Sinclair and B. R. Lawn, “An atomistic study of cracks inDiamond-Structure crystals,” (1972).[3] J. E. Sinclair, Philos. Mag. , 647 (1975).[4] E. Bitzek, J. R. Kermode, and P. Gumbsch, Int. J. Fract. ,13 (2015).[5] M. Marder, Int. J. Fract. , 169 (2016).[6] R. Thomson, C. Hsieh, and V. Rana, J. Appl. Phys. , 3154(1971).[7] R. Perez and P. Gumbsch, Phys. Rev. Lett. , 5347 (2000).[8] P. Gumbsch and R. M. Cannon, MRS Bull. , 15 (2000).[9] L. I. Slepyan, Soviet Physics - Doklady , 538 (1981).[10] J. R. Kermode, A. Gleizer, G. Kovel, L. Pastewka, G. Csányi,D. Sherman, and A. De Vita, Phys. Rev. Lett. , 135501(2015).[11] D. Holland and M. P. Marder, “Erratum: Ideal brittle fractureof silicon studied with molecular dynamics [phys. rev. lett. 80,746 (1998)],” (1998).[12] J. R. Kermode, T. Albaret, D. Sherman, N. Bernstein, P. Gumb-sch, M. C. Payne, G. Csányi, and A. De Vita, Nature , 1224(2008).[13] X. Li, The European Physical Journal B , 258 (2013).[14] X. Li, Journal of Applied Physics , 164314 (2014).[15] M. Buze, T. Hudson, and C. Ortner, ESAIM: MathematicalModelling and Numerical Analysis , 1821 (2020).[16] J. E. Sinclair, P. C. Gehlen, R. G. Hoagland, and J. P. Hirth, J.Appl. Phys. , 3890 (1978).[17] J. A. Yasi and D. R. Trinkle, Phys. Rev. E , 066706 (2012).[18] A. M. Z. Tan and D. R. Trinkle, Phys Rev E , 023308 (2016).[19] L. Pastewka, A. Klemenz, P. Gumbsch, and M. Moseler, Phys.Rev. B: Condens. Matter Mater. Phys. , 205410 (2013).[20] M. Ostoja-Starzewski, Probabilistic engineering mechanics ,112 (2006).[21] V. Ehrlacher, C. Ortner, and A. V. Shapeev, Arch. Ration.Mech. Anal. , 1217 (2016).[22] M. Buze, T. Hudson, and C. Ortner, Mathematical Models andMethods in Applied Sciences , 2469 (2019).[23] S. Lang, Fundamentals of Differential Geometry, GraduateTexts in Mathematics (New York Springer, 1999).[24] W.-J. Beyn, A. Champneys, E. Doedel, W. Govaerts, Y. A.Kuznetsov, and B. Sandstede, in
Handbook of Dynamical Sys-tems III: Towards Applications (2002). [25] R. D. Neidinger, SIAM review , 545 (2010).[26] S. Makri, C. Ortner, and J. R. Kermode, J. Chem. Phys. ,094109 (2019).[27] R. T. Rockafellar and R. J.-B. Wets, Variational analysis (Springer, 1998).[28] J. Tersoff, “New empirical approach for the structure and en-ergy of covalent systems,” (1988).[29] T. Kumagai, S. Izumi, S. Hara, and S. Sakai, Comput. Mater.Sci. , 457 (2007).[30] D. Fernandez-Torre, T. Albaret, and A. De Vita, Phys. Rev.Lett. , 185502 (2010).[31] G. C. Sih, P. C. Paris, and G. R. Irwin, Int. J. Fract. Mech. ,189 (1965).[32] A. H. Baker, E. R. Jessup, and T. Manteuffel, SIAM J. MatrixAnal. Appl. , 962 (2005).[33] D. Packwood, J. Kermode, L. Mones, N. Bernstein, J. Wool-ley, N. Gould, C. Ortner, and G. Csányi, J. Chem. Phys. ,164109 (2016).[34] A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli,R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Ham-mer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen,J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaas-bjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen,L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt,M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Wal-ter, Z. Zeng, and K. W. Jacobsen, J. Phys. Condens. Matter ,273002 (2017).[35] J. Kermode and L. Pastewka, “matscipy package,” (2020), https://github.com/libAtoms/matscipy .[36] A. P. Bartók, J. Kermode, N. Bernstein, and G. Csányi, Phys.Rev. X , 041048 (2018).[37] T. D. Swinburne and J. R. Kermode, Phys. Rev. B , 144102(2017).[38] A. Mesgarnejad, C. Pan, R. M. Erb, S. J. Shefelbine, andA. Karma, Phys Rev E , 013004 (2020).[39] J. J. Möller, A. Prakash, and E. Bitzek, Modell. Simul. Mater.Sci. Eng. , 055011 (2013).[40] A. M. Tahir, R. Janisch, and A. Hartmaier, Modell. Simul.Mater. Sci. Eng. , 075005 (2013).[41] J. J. Möller, E. Bitzek, R. Janisch, H. u. Hassan, and A. Hart-maier, J. Mater. Res. , 3750 (2018).[42] K. A. Cliffe, A. Spence, and S. J. Tavener, Acta Numerica9