Atomistic configurational forces in crystalline fracture

Configurational atomistic forces contribute to the configurational mechanics (i.e. non-equilibrium) problem that determines the release of total potential energy of an atomistic system upon variation of the atomistic positions relative to the initial atomic configuration. These forces drive energetically favorable irreversible reorganizations of the material body, and thus characterize the tendency of crystalline defects to propagate. In this work, we provide new expressions for the atomistic configurational forces for two realistic interatomic potentials, i.e. the embedded atom potential (EAM) for metals, and second generation reactive bond order (REBO-II) potential for hydrocarbons. We present a range of numerical examples involving quasistatic fracture for both FCC metals and mono and bi-layer graphene at zero Kelvin that demonstrate the ability to predict defect nucleation and evolution using the proposed atomistic configurational mechanics approach. Furthermore, we provide the contributions for each potential including two-body stretching, three-body mixed-mode stretchingbending, and four-body mixed-mode stretching-bending-twisting terms that make towards defect nucleation and propagation.


Introduction
We use the concept of spatial and material atomistic forces derived for various many-body potentials. In a quasi-static setting considered here, spatial atomistic forces contribute to the classical deformational (equilibrium) problem seeking to minimise the total potential energy of an atomistic system when varying the spatial atomistic positions. By contrast, material atomistic forces contribute to the configurational (nonequilibrium) problem that determines the release of total potential energy of an atomistic system when varying the material (initial) atomistic positions. Thereby, material atomistic forces characterise the tendency of generic atomistic defects to propagate, i.e., they drive energetically favourable re-organisations of the material atomistic configuration. Atomistic configurational forces provide a tool to study configurational changes which inherently occur during a fracture and failure process. Configurational forces are associated with any lattice defect irrespective of its type, i.e., point, line and planar defects which can be for example, vacancies, interstitials, substitutionals, Stone-Wales, 5-8-5, cracks, dislocations, twinning, stacking faults, surfaces, grain or twining boundaries. Taken together, every lattice irregularity results in configurational force driving its evolution relative to the material configuration.
In our prior works [1,2], we established the first principles of atomistic configurational mechanics. Therein, the concept of deformational and configurational atomistic pair, triplet and tetrad forces are introduced. Configurational pair force accounts for the changes in total energy due to the changes in the pair length. Variations in the area formed by a triplet of atoms is captured using the notion of configurational triplet force. In a similar fashion, the energetic variations associated with the variations in the tetrad volume defines configurational tetrad force. Notably, the proposed theoretical framework in [2] compares well to the energy-momentum structure of the Eshelby stress in the continuum configurational mechanics [3]. The continuum setting of configurational mechanics are presented for example in [4][5][6][7][8][9][10] as well as our own contributions to configurational mechanics in [11][12][13][14][15][16]. However, to a large extent we have not yet studied computational examples from the realm of atomistic fracture-and configurational-mechanics.
In this work, we initially elaborate on deformational and configurational formulations applicable on substantive materials. Next, we present a number of diverse numerical examples to demonstrate the applicability of the approach to crystalline fracture. We investigate in the derivation of energy release during the fracture of crystalline structures. We distinguish between spatial and material settings for two major categories of chemical bonds, i.e., metallic, and covalent bonding. We provide a novel framework to study the versatility and applicability of atomistic configurational forces in pure metals representing metallic bonding, and covalently bonded semiconductors. Atomistic configurational forces are practically used in [17] as a crack propagation criterion using interatomic bond deletions.
The above-mentioned cases allow us to assess the physical meaning of the configurational forces during the rearrangement of interatomic bonds i.e., bond dissociation and rotation in the presence of in-/out-ofplane externally imposed stretch, shear, bending and torsion. We use configurational forces to measure energy release and predict crack paths in a mono-/bi-layer graphene as a zero-gap semiconductor as well as in FCC structure of pure Copper. Graphene's exceptional mechanical and physical properties have been investigated in a number of studies over the decades, as well as authors contributions in [18][19][20][21][22]. By modeling graphene using a commonly-used reactive bond order potential [23] enables us to examine the contributions of two-(stretching), three-(bending), and four-(torsion) body interactions to failure and fracture. Similarly, by considering the two-body embedded atom potential for metals [24] provides a comparative configurational study between metallic and covalent interatomic interactions. In doing so, we show in this work that not only bond stretching, but also triplet bending and tetrad torsion play a significant role in brittle crack growth and propagation in crystalline lattices.
We organize this contribution as follows: after the introduction, Section 2 introduces the notation and elaborates on the kinematics of the problem. Thereafter in Sections 3, 4 and 5, we derive the governing equations from a deformational and configurational viewpoint. Next in Section 6, we derive particularly the deformational and configurational forces for an atomistic system modelled via reactive empirical bondorder potential and embedded-atom method. For the sake of demonstration, Section 7 focuses on numerical examples providing insights into the detailed analysis of configurational forces whereby the effectiveness of the methodology in the fracture of a single and double layers graphene as well as FCC lattice of copper are presented. Finally, Section 8 concludes the paper.

Kinematics
We assume a system of atoms described by their position X α occupying material (undeformed) configuration and x α denoting spatial (deformed) configuration of atomistic systems in ambient space E 3 . They are linked by a discrete, i.e., atom-wise motion with a time-like variable S ∈ R + , The vector-valued spatial and material pair direction follow as where spatial and material pair length vectors pointing from atom α to atoms β and the corresponding scalar-valued spatial and material pair lengths as x αβ := |x αβ | and X αβ := The introduction of spatial and material triplet angle is motivated by the classical concept of rotational spring in terms of angles between adjacent atomistic pairs. Here, we reassign triplet angles from the interval [0, π] to the interval [+1, − 1] by the cosine map, thus, ϕ αβγ := n αβ ⋅n αγ ≡ n αγ ⋅n αβ =: ϕ αγβ , and, In a similar fashion, spatial and material tetrad dihedral angle are introduced based on the concept of rotational springs in terms of dihedral angles between planes spanned by adjacent atomistic triplets. The tetrad angles from the interval [0, π] to [+1, − 1] by the cosine map, thus, ϕ αβγδ := n αβγ ⋅n αβδ ≡ n αβδ ⋅n αβγ =: ϕ αβδγ (6) and, We introduce spatial and material pair stretches as λ αβ := x αβ X αβ and Λ αβ := X αβ x αβ (7) Note that the objective and parity symmetric spatial and material pair stretches associated with each atom belonging to a specific pair are the same.
We also introduce a novel concept of spatial and material triplet twist as importantly, the objective and parity symmetric spatial and material triplet twists associated with each atom belonging to a specific atomistic triplets are in general distinct. We also introduce spatial and material tetrad dihedral twists as the objective and parity symmetric tetrad dihedral twits associated with each atom belonging to a specific tetrad are in general distinct. Finally, fully symmetric spatial and material pair projection tensors are defined as, x αβ and P ⊥ αβ := similarly, spatial and material projection tensors take the form of x βαγ and P ⊥ βαγ :=

Deformational and configurational mechanics
In this section, we reiterate briefly the main idea and the compact expressions of spatial and material settings proposed by Steinmann et al. [2] for two-, three-, and four-body atomistic systems.

Total potential energy
The total potential energy E of this atomistic system depends on the sets of spatial and material positions x ϵ and X ϵ , respectively, of all atoms. For the deformational problem, x ϵ denote the unknown variables, whereas X ϵ serve as a given parameterisation. For the configurational problem these roles are reversed, i.e., X ϵ denote the variable, whereas X ϵ serve as parameterisation. The total potential energy E consists of an internal and external contribution: The external potential energy, capturing interaction of finite atomistic system with external world, depends merely on the spatial atomistic position x ϵ of all atoms through the atom-wise external potential energy V α . That is, The internal potential energy, capturing internal interatomic interactions, depends in general on the sets of spatial and material position x ϵ and X ϵ of all atoms. We shall here consider the assumption that an additive expansion into pair, triplet and tetrad potentials is sufficiently accurate to approximate the complex energetic landscape dictated by underlaying quantum mechanics, thus

Equilibrium of spatial atomistic forces
For a quasi-static process, the equilibrium condition for spatial atomistic forces follow from minimising the total potential energy E of an atomistic system considered as a function of the spatial atomistic positions in x ϵ and parameterised in X ϵ . We denote spatial variations of the spatial atomistic positions x ϵ at fixed material atomistic positions X ϵ as D δ x ϵ . Thus, the minimum condition for the total potential energy E under spatial variations reads Consequently, the minimum condition for the total potential energy E of an atomistic system is given by The internal spatial force acting on atom α follow straightforwardly from the pair, triplet and tetrad potentials as Likewise, the external spatial force acting on atom α is given by explicitly by From the above atomistic equilibrium condition, we realize that the sum of all internal spatial forces k int α := k pair α + k triplet α + k tetrad α and external spatial forces k ext α acting on atoms α is zero.

Non-equilibrium of material atomistic forces
Here, a re-parameterisation of the total potential energy E of an atomistic system as a function of the set of material atomistic positions X ϵ , i.e., with X ϵ as the variables and x ϵ as the parameterisation. We denote material variations of the material atomistic positions X ϵ at fixed spatial atomistic positions x ϵ as d δ X ϵ . The material variation of the total potential energy E equates with the virtual energy release R δ that satisfies an inequality constraint in accordance with the second law of Thermodynamics, that is, Thereby, the material variation of the total potential energy E (note that d δ E ext ≡ 0) is given by The internal material force acting on atom α that follow from the pair, triplet, and tetrad potentials, respectively, as The external potential energy is due to the interaction with the world outside the atomistic system, e.g. due to gravity or contact forces at the boundary or electro-magneto fields. Here, it is assumed that the external potential energy depends solely on the spatial atomistic positions, i.e., E ext (x α ), thus there is no configurational force stemming from the variation of the external energy with respect to material changes.

Spatial energy density
The atomistic pair, triplet and tetrad contributions -E pair , E triplet and E tetrad -to the internal potential energy E int is expressed in terms of their densities -W αβ 0 , ω αβγ 0 and ω αβγδ 0 -per unit material pair length X αβ , angle Φ αβγ and dihedral angle Φ αβγδ , respectively. These energy densities are parameterised in terms of spatial pair stretch, triplet and tetrad twists,

Spatial interaction force
The stretch-based parameterisation of pair potential and twist-based parameterisation of the triplet and tetrad potentials lead to the corresponding spatial atomistic force expand as sums of corresponding spatial interaction forces The stretch-based spatial pair interaction force and twist-based spatial triplet and tetrad interaction forces are given by These interaction forces can be stated in terms of the stretch-based spatial pair, twist-based triplet and tetrad interaction force generator g pair αβ , g triplet αβ and g tetrad αβγ in the compact format as We define stretch-based spatial pair interaction force generator as the pair-wise sensitivities of the stretch-based pair potential. Thus, Similarly, we define twist-based spatial triplet and tetrad interaction force generators as pair-and triplet-wise sensitivities of twist-based triplet and tetrad potentials. In the following, we present two examples of such generators (for complete list see Appendix B.2 of [2]). Thus, and, We introduce the signed magnitude of the stretch-based spatial pair and twist-based spatial triplet and tetrad interaction force generators as Note that, the stretch-based spatial pair interaction force k pair αβ is the force exerted on atom α due to its interaction with atom β. It is oriented along the line connecting the pair α, β. The twist-based spatial triplet interaction force k triplet αβγ is the force on atom α due to its interaction with atom β and γ. It lies in the plane spanned by the triplet α, β and γ with g triplet and g triplet [αγ] oriented perpendicular to the lines connecting the pairs α, β and α, γ, respectively. The twist-based spatial tetrad interaction force k tetrad αβγδ is the force exerted on atom α due to its interaction with atoms β, γ and δ. It is assembled from contributions oriented perpendicular to the three planes spanned by the triplets α, β, γ and α, δ, γ as well as α, δ, β.

Material energy density
For the material setting, we consider the atomistic pair, triplet and tetrad contributions -E pair , E triplet and E tetrad -to the internal potential energy E int expanded in terms of their densities -W αβ t , ω αβγ t and ω αβγδ t -per unit spatial length x αβ , angle ϕ αβγ and dihedral angle ϕ αβγδ , respectively.
Observe that the spatial and material potential energy densities relate, for example, as whereby the material pair-stretch, triplet-and tetrad-twists Λ αβ , Ω αβγ and Ω αβγδ , respectively, serve to transform spatial into material densities. The internal energy described in this work includes the interatomic potentials suitable for studying zero Kelvin lattices in Molecular Statics (MS). However, an extension to finite temperature lattices in the context of MS can be provided by incorporating the atomistic frequency capturing the fluctuations occurring at non-zero temperatures. Through changing the variables in the definition of an empirical potential, we can study the effect of the atomistic frequency. A detailed explanation regarding the introduction of effective empirical potentials can be found in [25].

Material interaction force
The material atomistic forces corresponding to the stretch-based pair, twist-based triplet and tetrad interaction forces given by We also introduce the signend magnitudes of the stretch-based material pair and the twist-based material triplet and tetrad interaction force generators, It should also be noted that, for a finite temperature lattice, a further term denoted as thermal configurational force captures the influence of atomistic fluctuations in the variation of the total energy, however it is outside the scope of this study.

Many-body potential energy
In the following, we use the proposed theoretical framework to derive deformational and configurational settings of empirical manybody potentials. Embedded-atom method and reactive bond order model enable to assess configurational forces in the fracture of crystalline lattices e.g., copper and graphene, composed of the two-to fourbody atomistic systems.

Embedded-atom method
Functional forms of many-body potentials such as pair functional, e. g., embedded-atom method, have been developed to study materials e. g., metals, with better computational performance and higher accuracy in comparison with cluster potentials e.g., classical Morse potential. Generally, they take the form We consider embedded-atom method [24], thus ϕ αβ is a pair potential described by a Morse potential function [26] ϕ αβ ( where is a Morse function and H(x αβ ) is a unit step function. Equation (45) includes a cutoff function ψ(x αβ ) defined as The last term in Eq. (45) controls the strength of pairwise repulsion between atoms at short distances. The second term in Eq. (44), F (ρ α ), denotes the electronic energy, that is an embedding energy associated with placing an atom α in its electronic cloud, parameterized by an intermediary function ρ α . That is, the homogenous electron gas density, ρ α , is determined using the contributions from the neighbouring atoms as a linear superposition of electron density from every neighbour through the function ρ β , The electron density ρ β is a function of distance and represents a spherically-averaged electron density field around an isolated atom β, Embedded-atom method in total requires 20 material parameters and coefficients to fit a potential function using experimental data.

Spatial setting of embedded-atom method
The total potential energy E of an atomistic system is determined as a function of the spatial atomistic positions. Here W 0 αβ denotes the internal bond potential density per unit length of the material bond vector. Considering the same unit for energy density, deformational embedding energy density D 0 is defined. It is parameterized in terms of deformational electron density L α 0 . The external potential energies within the domain and at the boundary are denoted as V C α and V ∂C α , respectively, Pair potential density is parameterized using spatial stretch λ αβ , that is, whereby deformational Morse potential density is described as Deformational electron density, L α 0 , per unit material pair distance is defined as the linear sum of deformational electron densities in interaction with the neighboring atoms, Finally the deformational embedding energy density is described by, To derive the unit of energy per material pair distance (eV/Å), we redefine the fitting coefficients as presented (Table 1), Next, we derive pair deformational force k αβ oriented via the spatial pair direction resulting from the variations of pair and embedding energies due to the changes in spatial pair stretch, thus, subsequent functions derivatives to derive pair deformational force are Table 1 Function fitting parameters for deformational embedded-atom potential density of copper. Table 2 Fitting coefficient required to derive configurational embedded-atom potential density of copper. Table 3 Coefficients for REBO-II potential for solid-state carbon-based structures given in [23], which also include the parameters for the fifth-order polynomial spline function g βαγ .  presented in Appendix A.1. Finally, net deformational forces acting on an individual atom is derived as the cluster sum of pair deformational forces. That is,

Material setting of embedded-atom method
The total potential energy in the material setting is defined as the combination of internal and external energy densities, Configurational pair potential density is parameterized using material stretch Λ αβ , that is, where configurational Morse potential density is redefined based on material unit of density and corresponding pair stretch, that is, The configurational embedding energy takes the form of, where local configurational electron density is defined based on its linear summation of pair configurational electron density, such that, ultimately, we define pair configurational electron density, Some coefficients and variables are redefined here to maintain the unit of energy density based on spatial pair distance, such that we need (Table 2), Using material pair direction, configurational pair force acting on an atom α is defined as, subsequent functions derivatives to calculate configurational pair force are presented in Appendix A.2. Finally, net configurational forces acting on an individual atom is derived as the cluster sum of pair configurational forces. That is,

Second-generation reactive empirical bond order potential
The analytical form for the pair potential, V αβ , of the Brenner et al. [23] denoted as second-generation reactive empirical bond order (REBO-II) potential energy for solid-state carbon-based structures is given by the following functions with corresponding parameters listed in Table 3.  4. Illustration of copper FCC bulk deformation under in-plane increasing shear. In a.1a.3, the spatial configurations are presented. The bulk lattice is sliced in y direction to twodimensional layers. Configurational forces appear in 2D outer layers (b.1) takes a different form than an inner layer (c.1). Consecutive increase of shear (b.2) leads to the propagation of crack. Based on analogous configurational forces in c.2, the crack growth can be understood as the stretch of bonds along crack path. b.3 implies a dislocation glide occurring during deformation. Corresponding configurational forces in c.3 indicates rigorously the contributing atoms in crack propagation and dislocation glide.
The repulsive and attractive pairwise terms are defined by f R αβ and f A αβ , respectively, the cutoff function f c i j limits the contributions only to nearest-neighbors Fifth-order polynomial spline function g βαγ describes the energy of a triplet as a function of angle θ βαγ . The coefficients β n are derived via fitting the function to the experimental data.
The energy associated with the quadruplet dihedral angle is described as The ternary and quaternary interactions are involved in b αβ depending on the local coordination of atoms around atom α and the triplet and quadruplet angles between the atoms (α, β, γ, δ), The pi-bond term b σ− π αβ depends on the local coordination of atoms around atom α and the angle between atoms (β, α, γ), The dihedral function b DH αβ includes the quadruplet energy g γαβδ . The cutoff functions ensure that atoms α and β interact only with their nearest-neighbor atoms which in this case includes atoms γ and δ, respectively. This function ultimately describes the torsion between

Spatial setting of REBO-II potential
The term capturing the repulsion between carbon atoms is defined in spatial setting as deformational repulsive energy density, such that, Similar to the case of EAM potential we require to redefine various parameters and coefficients in REBO-II potential to provide an energy in terms of its density, in this way, for Eq. (76) we need, We also introduce the signed magnitude of stretch based interaction force associated with repulsive energy density, The spatial form of attractive force between the atoms depends on spatial pair stretch. Thus, and it is defined as energy density per unit of spatial pair length by modification of the coefficients The signed magnitude of attractive interaction force is derived as, As classical REBO-II potential requires to have a cut-off function limiting the interaction between carbon atoms to first nearest neighbour, we define it as a deformational energy density in the form of  The atoms in the adjacency of displacement-carrying zone experience the highest bending, thus reflected in the magnitude of configurational forces. b) is the corresponding spatial configuration illustrating the maximum curvature tolerated by SLGS before detachment.
The energy between a triplet of atoms depends on their spatial triplet twist ω βαγ , the energy is defined as its density per material triplet angle, Φ βαγ , such that, therefore we need, The signed magnitude associated with the variation of triplet energy due to spatial triplet twists, i.e., bending, is defined as Spatial tetrad twist ω γαβδ is used to define tetrad potential density, that is, The corresponding deformational signed magnitude of four-body interaction force is defined using infinitesimal changes in tetrad energy Deformational σ − π bond energy density identifies a double covalent bond between two carbon atoms. It depends on the triplet twist accounting for the spatial changes in the cosine of the angle bonds between a triplet of atoms. Similar to the classical case, using deformational cut-off function, the energy of triplet is limited to the contributions from nearest neighbour. Thus, the dependency of deformational σ − π energy density on spatial pair stretch and triplet twist requires the derivation of pair and triplet signed magnitudes, respectively, and, Deformational dihedral energy density describes the energy for a spatial rotation about dihedral angles for carbon double bonds per material tetrad twist and pair stretches, assuming only nearest neighbour interactions, we thus derive, (92) Fig. 15. The magnitude of configurational forces at the crack tip predicts the direction of crack growth. As it is illustrated in a, b and c using a black vector, the crack propagates in the direction that the configurational forces are larger. Fig. 16. Illustration of bi-crystalline mono-layer graphene. The second grain is oriented by θ = 9 ∘ from x axis. Fig. 17. Illustration of crack propagation in a bi-crystalline mono-layer graphene. The crack initiates and propagates in a zigzag fashion into the crystal. a depicts the contour plot for magnitude of pair configurational forces, whereby the configurational forces appear at the free surfaces, crack region and grain boundary. b presents that configurational forces not only trace the path of crack growth but also they predict the path of crack growth at the crack tip (see black vector). c presents that crack has chosen the path suggested in b to grow.
Next, we derive signed magnitude of interaction forces as a response to the changes in deformational dihedral energy due to the variations in spatial pair stretches and tetrad twist, namely, Net pair deformational force acting on atom α is defined using the cluster sum of pair deformational forces, where, Next, triplet deformational force takes the form of triplet force of REBO-II potential is defined as Finally, we derive atom-wise tetrad deformational force as, and tetrad deformational force is defined as,

Material setting of REBO-II potential
In the following, we present material atomistic energy and force formulations contributing to configurational problem to determine the release of energy in carbon-based atomistic structures modelled with REBO-II potential experiencing any variations in their material atomistic setting. Configurational repulsive energy density per unit of spatial pair distance is defined based on material pair stretch Next, we define configurational attractive energy density per unit of spatial pair distance, what follows is the derivation of the magnitude of configurational attractive force, presenting the variations in the attractive energy density with respect to the variations in material pair stretch, Similar to spatial setting we require to define configurational cut-off potential density per unit of spatial pair distance corresponding the signed magnitude of configurational cut-off force takes the form of, Configurational triplet energy density as a function of material triplet twist Ω βαγ measures the energy between a triplet of atoms per unit of spatial triplet angle ϕ βαγ , the signed magnitude of configurational triplet force presents the required force to bend the triplet at material configuration where we require to redefine a coefficient in the form of, Configurational tetrad energy density based on material tetrad twist per unit of spatial tetrad angle ϕ γαβδ is defined as the signed magnitude of configurational tetrad force determines the energy release due to the rotations of the atomistic tetrad at the material configuration, that is, Configurational σ − π potential density captures the double covalent indicates the material configuration based on the color mapping of magnitude of configurational forces. We observe the concentration of the tetrad forces in the vicinity of edges at step 1250, indicating significant bendings that are necessary to fold the tube. Moreover, in step 1500 configurational pair forces are focused in the vicinity of bukling zone, thus we can measure the associated energetical variations based on pair stretches.
interactions between carbon atoms whereby the spatial position of triplet is fixed to identify the energy of system at the material configuration, configurational σ − π forces are presented as their magnitude. We derive, and, Configurational dihedral energy density describes the energy for a material rotation about a dihedral angle per spatial tetrad twist and pair stretches, assuming only nearest neighbour interactions we have, Net pair configurational force acting on atom α is defined using the cluster sum of pair configurational forces, where, Next, triplet configurational force takes the form of triplet force of REBO-II potential is defined as Finally, we derive atom-wise tetrad configurational force as, tetrad configurational force is defined as,

Numerical examples
In this section, we study two-and three-dimensional atomistic crystalline models as a way of verifying configurational formulation of Section 6 for a wide range of interatomic potentials representing different types of atomic bonding, i.e., covalent and metallic. Considering the suitability of REBO-II potential [23] for studying the fracture and failure of graphene, we provide various models composed of carbon under quasi-static process. Moreover, we consider fracture of copper as it is originally investigated in the introduction of EAM potential [24]. In Fig. 1 a flowchart for solving following examples is presented. The material (initial) positions of atoms are initially recorded. External energy is then applied to the boundary atoms in the consecutive steps, requiring minimization of total potential energy. Thereby, we obtain spatial atomistic positions to compute configurational forces. Next, the magnitude of pair, triplet and tetrad configurational forces are compared with a threshold allowing us to delete or preserve an interatomic bond using Configurational-Force-Criterion (CFC) proposed by Birang O. and Steinmann [17].

Fracture of copper under tensile loading
In the following example, we study fracture of FCC copper to assess our proposed configurational formulation for embedded-atom method (see Section 6.1). An edge crack of length 27.014 Å along [1 0 0] axis is introduced on (0 0 1) plane of a single crystal FCC copper of size 68.685 × 14.46 × 68.685 Å 3 . The pre-crack is formed by turning off the interaction between yellow and blue atoms. Mode-I loading is applied to the pink atoms by uniform stretch in z direction, whereby other boundaries are traction-free. The lattice constant of FCC copper is 3.615 Å and it is arranged in cubic directions, i.e., x[1 0 0], y[0 1 0] and z[0 0 1] (see Fig. 2). All boundaries are constrained by non-periodic boundary condition. After crack insertion, the lattice is minimised using L-BFGS optimization algorithm. An increasing displacement by ΔZ = 0.005 Å in 1600 steps is then applied. Upon imposing the displacement the pink atoms are frozen to find the atomistic positions at equilibrium. Zero Kelvin temperature of lattice suppresses the thermally activated dislocations.
We evaluate pair configurational forces at every equilibrium configuration as pair forces are applicable to functional potentials such as EAM. We observe that in this model the crack propagates in a brittle manner (Fig. 3 (b)), a similar behavior of pre-defined crack in single crystalline copper under tensile is reported in [27]. The brittle propagation of crack is associated with symmetric geometry of assumed lattice as well as symmetry of imposed boundary conditions. Moreover, the reversibility of crack propagation in such a model due to the preservation of all bonds prevents the observation of ductile crack. The visualization of configurational forces (Fig. 3 (c)) demonstrates that these forces appear at free surfaces, as well as crack surfaces. In Fig. 3 (a), we present the contour plot of z component of these forces, thus they trace crack growth precisely.

Fracture of copper under shear loading
In the previous example, we studied configurational changes for a perfect brittle crack growth in bulk copper. In line with it, in the following example we aim to study the capability of configurational forces to study the mechanism of dislocations under shear loading. A uniform displacement of 0.005 Å in 1600 steps in x direction is applied to the lattice. Initially, the deformation of the lattice leads to the propagation of the pre-defined crack, next we observe an emission of a dislocation from the crack tip as illustrated in Fig. 4 (a.1-a.3). The evaluation of net configurational forces at various 2D slices of 3D cube illustrate that during crack growth configurational fores of each layer are distinctive. In Fig. 4 (b.1-c.1) at step 0, we have chosen an inner and outer layer to compute configurational forces, as it is observed, the outer layer entails configurational forces on all atoms including crack zone as it represents a free boundary of lattice, however in an inner layer configurational forces appear at the pre-defined crack zone and edge boundary atoms. The increment of shear leads to further reorganizations at material configuration, thus new configurational forces emerge in the crack path (see Fig. 4 (c.2)). Ultimately the angle and path of dislocation can be precisely identified using configurational forces, See Fig. 4 (c.3) depicting plainly the slip plane (1,0,1)[0,1,1]. These forces determine accurately the direction of energy release and the contributing atoms in the complex mechanism of dislocation nucleation/glide. This example classifies the capability of atomistic configurational forces to study plasticity occurring at the crack tip to determine subjectively and objectively the sources of local energy release under configurational changes.

Fracture of mono-layer graphene under mixed-mode bendingdisplacement
Graphene monolayers exhibit excellent load-carrying capability, even under out-of-plane displacement, due to its strength under large bending and torsion. Duan and Wang [28] carried out study to assess bending behavior of a circular graphene under point-loading using Molecular Dynamics (MD) and Von-Karman plate theory. Their results from MD and Von-Karman theory show that there is stress concentration near loading region and free boundaries.
We consider a single-layer graphene sheet (SLGS) under out-of-plane displacement to analyze pair, triplet and tetrad configurational forces at a material configuration that corresponds to a specific spatial configuration of SLGS.
A hexagonal sheet including 14,200 atoms is assumed in Fig. 5. The atoms at the boundary depicted as pink circles are constrained in all directions. We study a quasi-statics model thus the atoms do not possess any momenta. Initially at unstretched state, the total energy of the plane is minimized to obtain initial relaxed lattice. Next, step-wise displacement by Δz = 0.25 (Å) in +z direction is applied to the central yellow atoms followed by energy minimization. We use quick-min optimization algorithm [29] to find the minimum energy path. The position of atoms at every displacement step is recorded as the convergence criterion which is the norm of total forces acting on plane is satisfied, i.e., | ∑ α∈I C k α | = 1e − 9. Based on the recorded spatial position of atoms, obtained from the solution of standard MS problem, and using initially specified material position of atoms, we calculate the material pair stretch, triplet and tetrad twists, Λ αβ , Ω βαγ , Ω γαβδ , to evaluate local pair, triplet and tetrad configurational forces, respectively.
In [17], it is presented that configurational forces can identify free boundaries and the location of pre-defined defects in a crystalline lattice involving short-range interactions. Here, we measure configurational forces resulting from pair, triplet and tetrad configurational terms as illustrated in Fig. 6. At initial undeformed configuration (Fig. 6 (A.sp)), configurational forces are observed on the atoms defining fixed boundaries. It can be explained by considering that the edge of the plane in two-dimensional models or surface of a bulk in three-dimensional models behave as a defect. Moreover, during the deformation of SLGS, configurational forces at the boundary experience variations in their energy.
Next, we observe the emergence of new configurational forces on the atoms in the vicinity of load-carrying zone during the interatomic debondings as illustrated in Fig. 6 (B.ma)-(D.ma). It can be inferred from the pictures, the configurational forces are compatible with the spatial configurations illustrated in Fig. 6 (B.sp)-(D.sp). The configurational forces result from stretches and bendings on interface atoms. In this example, we show that the energy release associated with the fracture of SLGS under a mixed mode stretching-bending can be determined individually using pair and triplet configurational forces.
According to Fig. 7 (b), the angle between a triplet of atoms at the interface changes from undeformed configuration (n = 0) to the deformed one (n = 106). Therefore, we observe the emergence of triplet configurational forces before the detachment of load-carrying zone at step 106 in Fig. 8 (a), compatible with displacement illustrated in Fig. 8 (b).

Fracture of mono-and bi-layer graphene under tensile loading
Zhang et al. [30] illustrated an armchair crack morphology in a bilayer graphene that is similar to a mono-layer graphene, however they ignored the influence of interlayer interaction that might influence crack growth path. On the other hand, Jang et al. [31] experimentally presented the asynchronous crack with dissimilar path occurs in fracture of bilayer graphene. In the following example, using configurational approach we aim to examine the fracture of a mono-layer zigzag crack in comparison to a bi-layer zigzag crack entailing interlayer Van der Waals forces. We compare the influence of pair, triplet and tetrad energy release and emergence of subsequent configurational forces in fracture of both layers, particularly their differences at the crack tip. The bi-layer graphene with size by 170.53 × 59.08 Å 2 has a distance by 3.34 Å between two layers. The pink atoms are stretched by 0.02 Å in 200 steps. A pre-defined crack is introduced in both layers by deleting three layers of atoms in zig-zag direction (see Fig. 9). The equilibrium configuration is obtained using L-BFGS optimization algorithm considering energy derivative by 1e − 10 as convergence criterion. For stable crack, we observe wrinkles occurring due to the consideration of van der Waals interlayer interactions (see Fig. 10 (a.ma)). However, in corresponding material ( Fig. 10 (a.sp)) we observe only configurational forces at the free boundaries and crack surfaces, thus the wrinkles do not lead to the emergence of primary pair configurational forces.
The propagation of crack into the graphene layers are presented in Fig. 11 (a.sp) and (b.sp) for various snapshots 251 and 252. Moreover, at material configuration we observe that the magnitude of pair configurational forces are consistent with crack growth at the spatial configuration. Next, we compare the magnitude of pair, triplet and tetrad configurational forces, namely, Fig. 12 a, b and c at various snapshots. At it is comprehensible from the magnitude of contour bar, the largest energy release during crack growth occurs due to bond stretches, however the contribution of bending and torsion are remarkable.
Furthermore, the crack morphology at both layers are compared. As illustrated in Fig. 13, using the magnitude of pair configurational forces we observe that crack length at the top layer is longer than the bottom layer. The dissimilarity of crack paths might stem from highly non-linear deformation of system due to interlayer interactions.
In the following, the fracture of a mono-layer graphene is compared with the bilayer. In this model, we use the geometry of bilayer, while only one layer of it is considered. A stretch by 0.1 Å, in 500 steps is applied to the boundary atoms. From Fig. 14 we observe that for a similar geometry, the crack morphology in a mono-layer graphene is different than a bilayer graphene. Moreover, mono-layer graphene shows a higher fracture toughness in comparison to bilayer graphene, thus larger external energy is required to fracture the mono-layer. This behavior can be associated with the intensive wrinkling of bilayer, such that it makes the system highly brittle. Moreover, as it is illustrated in Fig. 12, bilayer fracture is accompanied with large triplet and tetrad energy release, while these terms do not influence fracture of monolayer.
Finally, we assess the magnitude of pair configurational forces in the prediction of crack path. As it is shown in Fig. 15 at snapshot n = 250, we can see that due to the artificial blunting at the tip, crack choose the path suggested by configurational force (depicted by black vector) as it leads to crack growth in armchair direction; that is crack grows toward top-right direction rather than bottom-right or zigzag direction. Similarly at the crack tip in snapshot 251, we observe that crack propagates in a direction that the magnitudes of configurational forces of an atomistic pair are largest at the tip, i.e., bottom-right direction. Magnified Fig. 15 at n = 252, shows that the prediction suggested at previous snapshot is satisfied.

Fracture of bicrystalline mono-layer graphene under tensile loading
In the following example we aim to examine the ability of configurational approach to predict fracture propagation through a grain boundary where the bonding environment is much more complex than for a single crystalline structure. Figure 16 illustrates a bi-crystalline mono-layer graphene under increasing displacement at the top and bottom boundaries (pink atoms) by 0.046 Å at 1000 steps. The energy of system is minimized using L-BFGS algorithm to derive the configuration at equilibrium. What follows is the evaluation of configurational forces. At the initial material configuration we observe that configurational forces emerge at the edge and the grain boundaries, see the contour plot of pair configurational forces in Fig. 17 (a). Moreover, in Fig. 17 (b) we observe that configurational forces appear by the initiation of crack, as well as predicting the crack propagation path at the grain boundary; see magnified image at the boundary where atoms with largest configurational forces, that is red atoms illustrate the prospective candidates for debonding. Finally, in Fig. 17 (c) we observe that the prediction at step (b) is accurate and crack has chosen the suggested path to propagate, however it should be noted that the prediction is only possible for one lattice unit.
An extension for studying fracture occurring at grain boundaries (GBs) using the configurational mechanics approach might require the consideration of an analytical energy describing in particular GB, for example the Bulatov-Reed-Kumar five-parameter grain boundary energy function for face-centered cubic metals [32]. In the continuum setting, configurational forces indicating the energetical variations occurring at surfaces and interfaces is developed by Steinmann [33] which might help the development of a framework to study configurational mechanics of GBs at atomistic scale.

Buckling of carbon nanotube under torsion
In the following, we study the fracture of a carbon nanotube (CNT) at zero Kelvin using configurational forces. The goal is the determination of the contribution of two-body stretch, three-body bending and four-body torsion under quasi-static torsion to an armchair single-walled CNT as illustrated in Fig. 18.
In recent years, torsional behavior of CNTs is studied in a number of studies at molecular scale (see for example [34,35]). Yu et al. [36] investigated in the mechanical behavior of CNTs under torsion, therein, the torsion of carbon nanotube causes the bukling with dumbbell-shape cross section. Khoei et al. [35] reported the reversible nature of CNTs bukling under torsion. Energy variation of CNTs accompanied with bukling under twist of CNTs are also investigated in other works [37][38][39][40]. In a separate study [41], the influence of temperature and various types of defects on the torsional capacity of CNTs are investigated. The previous studies of CNTs provide a detailed analysis on the structural evolution of CNTs during torsional deformation, however a study on the configurational changes during geometrically nonlinear deformation of CNTs is not yet reported. This motivates the following example characterizing in detail the structural evolution of twisted single-walled CNTs modelled by REBO-II potential and analyses of various emerging configurational forces which correlate well with the structural evolution.
We model a perfect CNT that is twisted by Δθ = 0.9 ∘ in 400 steps, whereby the rotation axis is along z direction that goes through origin point (0,0,0). Upon applying the torque, the pink atoms are constrained in all degrees of freedom to minimize the potential energy with respect to the free atoms using L-BFGS algorithm using the second norm of forces as the convergence criterion by 1e − 10. Figure 19 (a) presents a series of snapshots from the spatial configuration of the armchair CNT. The pictures are plotted in a parallel view to tube axis. The single-walled carbon nanotube during this study does not fail, indicating its extraordinary flexibility under rotational loading. After a rotational angle of 225 ∘ , initial folding occurs all along the tube. As the rotation increases by 262 ∘ , the tube wall folds. Finally, at the rotation of 270 ∘ , we observe the collapse thus buckling occurs at the center of tube, introducing concentrated plasticity. For the sake of clarification, we select a ring of the armchair CNT and plot its evolution above the associated tube, as shown in Fig. 19 (a-c). The ring located in the height by 50 Å from the bottom indicates that the imposed shear to the tube during folding changes the circular cross section to elliptical, revealing large variations in the angle.
For such a highly nonlinear deformation, we present local configurational forces during the progress of twist in Fig. 19 (b). At snapshot 1250 for the rotation angle by 225 ∘ , triplet configurational forces illustrate locally the atoms contributing in the process of twisting and folding. Next, pair configurational forces at the angle 270 ∘ of snapshot 1500 illustrate the exact location of buckling, thus simplifying our analysis on the changes of total potential energy. In Table 4, we present the contribution of pair, triplet and tetrad configurational forces during torsion of CNT. The most significant source of energy release is the pair stretches during folding, however the triplet bending and tetrad torsion are also sufficiently large to be considered in the failure of CNTs.

Conclusion
In our prior works, we established the theoretical foundation to describe the discrete configurational mechanics of crystalline lattices, which interact via many-body potentials. In this work, we briefly reiterated the theory, but focused on examining the usefulness of discrete configurational forces to analyze quasi-static fracture of covalently bonded graphene and metallically bonded FCC copper, which are modelled by REBO-II and EAM potentials, respectively. The chosen empirical potentials capture interatomic stretches, bending and torsion. We utilized our kinematical parameters to describe configurational expansions, whereby the notions of pair, triplet and tetrad configurational forces capture the release of energy due to the changes from the material (initial) positions. These findings confirm that discrete configurational forces offer a means to predict both the mechanics driving the crack nucleation and propagation, as well as the fracture paths, including unpredictable, zigzag crack paths. Furthermore, we demonstrated that not only bond stretching, but also triplet bending and tetrad torsion play a significant role in brittle crack growth and propagation in crystalline lattices. We envision that this approach can be used to study fracture in other crystalline nanostructures, such as nanowires, while also incorporating thermal effects on the fracture mechanics.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
where the derivative of Morse function density with respect to the spatial stretch is derived as, similarly the derivative of cutoff function with respect to the spatial stretch is derived, that is, The infinitesimal changes in the deformational embedding energy by the variation of deformational electron density L α Finally, we evaluate the variation of electron density at the position of an isolated atom with respect to the changes occurring in the spatial stretch

A2.
In the following we explain the required mathematical steps to derive configurational forces in embedded-atom method. The derivative of pair potential density with respect to material bond stretch takes the form of, where the derivative of Morse function density with respect to the material stretch is derived as, similarly the derivative of cutoff function with respect to the material stretch is derived, that is, The infinitesimal changes in the configurational embedding energy by the variation of deformational electron density Finally, we evaluate the variation of electron density at the position of an isolated atom with respect to the changes occurring in the material stretch Zhang et al. [30] measured fracture toughness of centrally pre-cracked mono-and bi-layer graphene using nanomechanical tensile testing. They validated the results via performing MD study at a constant temperature to focus on crack growth only due to mechanical loading. The experimental and numerical results in mono-and bi-layer graphene show that fast brittle fracture occurs which manifests the formation of long, atomically sharp crack. Budarapu et al. [42] carried out MD study of mono-layer graphene to study the variations in interatomic pair length and triplet angle for a pre-crack in zigzag and armchair directions. Their study shows that for zigzag pre-crack, critical pair stretch determines the crack propagation, however in armchair fracture pair stretches and triplet angles are both determining factors in crack growth. Le and Batra [43] worked on the influence of pre-crack length and distribution in crack tip speed for a MD study, their results show that crack tip speed rapidly reaches steady state, thus it can be interpreted as its independency from initial crack length and distribution. However, their efforts for measuring fracture toughness was not fruitful as  Illustration of crack propagation for pre-defined screening edge crack (case 3). The crack morphology is presented at different displacement steps. Associated material configuration for every selected displacement is also shown to present pair configurational forces emerging on the atoms located in the crack path. For the sake of comprehensible illustration, the material configuration is sliced from crack line. Fig. A.4. Spatial configuration of SLGS containing screening central crack is presented in (a). The magnitude of pair configurational forces is presented in (b) consistent with spatial configuration. Largest configurational forces emerge at the free boundaries and crack path. the discrepancy between measured toughness for various far-field strains was large. Here, we present fracture of a graphene strip as illustrated in Fig. A.1. We assume that the strip is subjected to four types of pre-defined cracks individually. We analyze the toughness of SLGS and assess the utility of configurational forces in the determination of crack morphology. To this end, in every case of study, we apply a total displacement of Δx = 0.12 (Å) in 300 steps to the boundary atoms at top and bottom layers in armchair direction. Initially, at undeformed state, the strip is relaxed by minimization of energy using quick-min approach with energy gradient as the convergence criterion by 1e-9. Next, at every step, the energy of the system is minimized upon applying the displacement and keeping the boundary atoms fixed in all degrees of freedom during minimization. The positions of atoms are also recorded step-wise to calculate configurational parameters for the evaluation of configurational forces. We model the system using REBO-II potential, however we assume that minimum and maximum cutoff radi, d and D, respectively, for nearest-neighbor interactions is set to 1.92 (Å) as it provides the brittle propagation of crack in a mono-layer graphene (see Methods in [30]). The geometry of pre-defined cracks, inserted into the strip, are presented in Fig. A.2.
We evaluate pair configurational forces for the case of screening edge crack. The location of crack at different displacement steps are presented in Fig. A.3, whereby at material configuration we observe configurational forces appear at crack edges which are identical due to the model's symmetry. Fig. A.4 presents the contour plot of magnitude of pair configurational forces for the case of central crack. Next, we observe the case of central and edge crack (Fig. A.5), thus the crack is introduced by removing some atomistic layers. Finally in Fig. A.6, the contour plot of pair and triplet configurational forces for edge blunting crack are illustrated.