Skip to main content Accessibility help


  • Access
  • Open access



      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Fracture ab initio: A force-based scaling law for atomistically informed continuum models
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Fracture ab initio: A force-based scaling law for atomistically informed continuum models
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Fracture ab initio: A force-based scaling law for atomistically informed continuum models
        Available formats
Export citation


In fracture mechanics, established methods exist to model the stability of a crack tip or the kinetics of crack growth on both the atomic and the macroscopic scale. However, approaches to bridge the two scales still face the challenge in terms of directly converting the atomic forces at which bonds break into meaningful continuum mechanical failure stresses. Here we use two atomistic methods to investigate cleavage fracture of brittle materials: (i) we analyze the forces in front of a sharp crack and (ii) we study the bond breaking process during rigid body separation of half crystals without elastic relaxation. The comparison demonstrates the ability of the latter scheme, which is often used in ab initio density functional theory calculations, to model the bonding situation at a crack tip. Furthermore, we confirm the applicability of linear elastic fracture mechanics in the nanometer range close to crack tips in brittle materials. Based on these observations, a fracture mechanics model is developed to scale the critical atomic forces for bond breaking into relevant continuum mechanical quantities in the form of an atomistically informed scale-sensitive traction separation law. Such failure criteria can then be applied to describe fracture processes on larger length scales, e.g., in cohesive zone models or extended finite element models.



Current address: Materials Modeling, Fraunhofer Institute for Mechanics of Materials, Wöhlerstr. 11, 79108 Freiburg, Germany.

This paper has been selected as an Invited Feature Paper.


Resistance to crack propagation is undoubtedly one of the most important properties of structural materials. However, our current mechanistic understanding of the energy-dissipation processes involved on the various scales and of their interplay is not sufficiently advanced to predict the fracture behavior of materials ab initio.1,2 As fracture is ultimately caused by the breaking of atomic bonds,3 there is a high demand for atomistically informed fracture criteria. An overview of ongoing attempts to realize them in hierarchical and concurrent schemes can be found, e.g., in Refs. 4 and 5. Here we propose a hierarchical scheme that meets two common challenges: to extract information on the atomic scale in an efficient manner and to coarse-grain this information to obtain practical and meaningful failure criteria for a finite element (FE) implementation on the mesoscale. Our approach is demonstrated by making use of cohesive zones (CZs), but is not limited to this model.

A so-called cohesive zone model is a seemingly natural way to incorporate atomistic information in a continuum model of fracture. In a traditional CZ model, the fracture process zone (i.e., the region in which nonlinear processes occur) ahead of the crack tip is represented by cohesive surfaces (respectively a narrow band of vanishing thickness) which experience cohesive tractions described by a constitutive law, the so-called traction-separation (TS) law. Crack growth occurs when the separation at the tail of the CZ, i.e., at the physical crack tip, reaches a critical value at which the cohesive traction vanishes. The strip yield model of Dugdale,6 which was originally proposed for estimating the size of the plastic zone around a crack tip, is regarded as one of the first CZ-type models. Barenblatt7 also proposed a cohesive fracture concept with the aim to eliminate the crack tip stress singularity in classical linear elastic fracture mechanics (LEFM). For fracture of perfectly brittle materials, the descriptions by LEFM and the CZ model approach are equivalent if the material behavior in the zone around the crack tip is determined by the stress-intensity (K) field and if the cohesive energy density is identical to the critical energy release rate.8 The first condition is met when the fracture process zone is smaller than the K-dominated zone (the region with a significant gradient in the stress), which is the case in brittle materials.9,10

A crucial point to consider in practical applications within a finite element method (FEM) implementation is that CZ-based numerical techniques require cracks to propagate along the element boundaries. Therefore, an inherent nonvanishing amount of mesh dependency arises, which can become dominant for three-dimensional simulations in a homogeneous material.11 For instance, no converged solutions from the simulation of progressive delamination were obtained above a critical mesh size.12 As a general rule, Hillerborg et al.13 suggested a maximum characteristic length ${l_{{\rm{coh}}}} \approx {{{G_{\rm{c}}}} / {\sigma _{\rm{c}}^2}}$ for the CZ under mode I loading, based on the cohesive energy G c and the critical stress to initiate crack advance σc. This characteristic length is always smaller than the radius of the K-dominated zone at critical loading, but also larger than the fracture process zone. In our approach, we couple the mesh size to the critical stress, which is the object of a scaling procedure explained in Sec. III. Furthermore, we do not only treat the fracture process zone as a CZ, but sample the complete K-dominated region on the cleavage plane with cohesive elements. This way, we obtain mesh-independent results for the fracture toughness under mode I loading.

In principle, the TS law can capture the physics of deformation and failure on the microscale, provided its shape is chosen carefully and that it includes important parameters such as maximum strength, work of separation, and critical displacement.14 Within the hierarchical modeling approaches, there are basically two philosophies about how to incorporate this information in TS laws, which differ in the way in which dissipative processes enter the modeling framework.

TS laws derived from atomistic simulations of homogeneous deformation of a single crystal, a bicrystal, or even model sections of polycrystalline microstructures can produce a stress–strain curve which reflects both elastic and dissipative processes during fracture.1521 However, dissipating the energy which is released during crack initiation and growth can be severely constrained in microstructures22 so that the fracture toughness becomes dependent on length scale and geometry. This means that the TS law would have to be adapted for each specific microstructure. Instead, we want to promote an approach in which the process of material separation by breaking of atomic bonds is separated from energy dissipative processes, such as plastic deformation accommodated by generation and motion of dislocations. The description of bond breaking within the CZ can be combined with different levels of description of the plastic deformation in the material adjacent to the CZ. On the smallest scale, plasticity is described based on discrete dislocation dynamics models in two2325 or three dimensions.26 Further developments include the formulation of the virtual internal bond model of cohesive fracture and its application to nanomaterials.27,28

Following the failure criterion by Griffith, the TS law for the bond-breaking process will have to capture the work of separation of the material, but to predict reliably the trends in cohesive behavior, the strength of the bonds is required as well.29 The strength of interatomic bonds is reflected in the theoretical strength of a material, which can be determined via ab initio density functional theory (DFT) calculations for single crystals30,31 as well as interfaces with and without impurities.3234 However, this procedure assumes homogeneous expansion perpendicular to, or cleavage along, an infinite crystallographic plane. During such a process, all bonds across the anticipated cleavage plane are elongated by the same amount and rupture at the same time. By contrast, individual pairs of atoms across the cleavage plane in the vicinity of the crack tip experience different environments, due to the gradient of the elastic displacement ahead of a crack tip, as shown in Fig. 1. Thus, the equivalence of the stresses obtained from such a rigid body separation (RBS)-type calculation with the critical stress needed to advance a crack tip still has to be demonstrated. This is one of the aims of the present paper (see Sec. II).

FIG. 1. Simulation setups for rigid-body separations (a) and crack-tip configurations (b). Note that not all interacting nearest-neighbor (NN) pairs (green; color online) are shown here due to the two-dimensional visualization. In fact the number of interacting NN pairs is twice the number of interacting 2nd NN pairs.

Embedding an atomistics-based TS law into a fracture model via CZs in a FE scheme has also practical issues because the theoretical strength is of the order of GPa and the critical separation is of the order of Å.35 The stresses on the order of the theoretical strength occur very locally in the extreme stress gradient in front of the crack tip, i.e., usually within the first few layers of atoms. In contrast to this, any continuum modeling is based on describing average stresses in some finite volume such that typical stresses occurring in FE models are much smaller. Therefore, Serebrinsky et al. first used a scaling relation for TS laws36 derived from the results of ab initio calculations based on the concepts introduced in Refs. 37 and 38. It assumes that the critical stress needed to separate two half crystals depends on the amount of elastic energy that is stored and released from the material next to the cleavage plane, which in turn scales with the system size. This can easily be observed in atomistic supercell calculations in which two half crystals are shifted and the atomic positions are relaxed afterwards. If the theoretical strength is simply identified with the maximum of the derivative of the energy–displacement curve, it decreases with an increasing number of planes parallel to the cleavage plane.37 However, the amount of energy that localizes and is subsequently released during the fracture process is in fact finite,39 and the seemingly size-dependent fracture stress can be normalized to obtain the unique strength of the material.40 Thus, a scaling approach via the elastic energy is practical, but not well justified, and a physically meaningful coarse-graining scheme providing an unambiguous connection between the macroscopic and microscopic stresses is still needed.

With the work at hand, we want to promote tensile tests carried out by RBSs and translations in combination with a scaling approach based on pairwise forces. RBS schemes can be used to characterize single crystals and grain boundaries on the same footing,41 provide a three-dimensional description of single crystal42 and grain boundary separation,43 and can be extended to include impurity effects.35 Most importantly, the derived critical stresses are size independent without any normalization. Our scaling approach based on interatomic forces will allow us to relate this fundamental atomistic quantity to the required mesh-sensitive model parameters in a rigorous way.

In the following, the results of atomistic simulations of RBS are compared with those of simulations of K-dominated cracks to demonstrate the equivalence of the interatomic forces in both situations; see Sec. II. This is the starting point for our novel approach to scale the atomic-level properties into parameters for continuum models in a physically motivated way, which is presented in Sec. III. This novel approach is tested in the framework of CZ modeling of cleavage fracture in an FE model, which demonstrates its advantages and shortcomings. Our findings are summarized in Sec. IV.


A. Interatomic potentials and simulation setups

We chose tungsten as a brittle, elastically isotropic model material to study cleavage fracture under static loading conditions. The atomic interactions were described either in the framework of ab initio DFT calculations using the local density approximation (LDA) or by using different potentials of the embedded atom method (EAM) type. DFT calculations were carried out using the ABINIT electronic structure code,4446 employing the ultra-soft pseudopotential of Hartwigsen, Goedecker and Hutter.47 For the RBSs described below, a k-point mesh of 8 × 8 × 2 k-points was used, which, after Fourier transformation, corresponds to a maximum k-point spacing of 47.45 a.u. (25 Å). The energy cut-off for the plane-wave basis set was 16 Ha (435 eV).

The EAM potentials used are the classical Ackland–Thetford–Finnis–Sinclair (ATFS) potential from the 1980s48,49 and the more recent potential by Wang.50 Their fracture-relevant properties and cutoff radii r cut are presented in Table I. All atomistic simulations were performed with the molecular dynamics software package IMD51,52 using the MIK53 and FIRE54 relaxators.

TABLE I. Summary of fracture-relevant potential properties of the EAM potentials used for tungsten (ATFS48,49 and Wang50): the lattice parameter a 0 at 0 K, the bulk equilibrium energy E 0, and the elastic constants C ij, the (100) surface energy, and the cut-off radius.

The simulation setups for the corresponding atomistic simulations are shown in Fig. 1. Figure 1(a) shows the case of rigid-body separations, which is often used to determine the theoretical strength via ab initio tensile tests. In this setup, two crystal halves are iteratively separated without allowing the atoms to relax their positions. As mentioned in the introduction, this strategy prevents size effects, which would occur if the crystallographic layers were allowed to relax, and furthermore enables a unified description of single crystals and grain boundaries.41 Note, however, that the strength values determined from methods that include atomic relaxations and lateral contraction of the sample can differ significantly, particularly in the presence of defects, such as grain boundaries.55 The implications will be discussed below. Periodic boundary conditions (PBC) are used in all directions. In the case of a grain boundary calculation, this actually requires the inclusion of two cleavage planes as there would be two grain boundaries in the simulation cell. This is not necessary in the single crystal, but to be able to use the same setup, we did so for consistency.

The size of the simulation box in the present study was L x × L y × L z ≈ 10 × 60 × 10 Å3. The separation distance Δy was varied from −0.5 Å to 2.5 Å in incremental steps of dΔy = 0.025 Å. The normal stress σat acting between the atoms of the two crystal halves at a certain separation distance Δy was then calculated from the average potential energy per atom, E pot, as

(1)$${\sigma _{{\rm{at}}}} = {{{N_{{\rm{atoms}}}}} \over {2{A_{xz}}}}{{{\rm{d}}{E_{{\rm{pot}}}}} \over {{\rm{d}}{\Delta _y}}}\quad ,$$

where N atoms is the number of atoms in the configuration and A xz is the cross-sectional area of the created surface. The factor of 2 results from the PBC.

Figure 1(b) shows the setup used to determine the fracture toughness and structure of an atomistically sharp crack. Although this setup can also be used to study cracks with DFT,56 it is usually used in classical atomistic simulations with interatomic potentials. In this setup, the linear-elastic solution of the crack-tip displacement field under plane-strain conditions, u x,y(r, θ) is directly applied to a cylindrical atomistic configuration with PBC.57,58 The displacements are determined for a given stress intensity factor K I and depend on the distance from the central crack tip, r, and the inclination angle θ of the direction of r and the crack plane. The critical stress intensity factor for the onset of crack propagation, K Ic, was determined by iteratively increasing K I, the corresponding Δu x,y, and the subsequent relaxation of the atoms in the interior of the configuration (with a minimum distance of 2r cut from the outer surface). The value of K I at which the crack propagates by brittle bond breaking was then defined as K Ic (we will denote this quantity as $K_{{\rm{Ic}}}^{{\rm{at}}}$ hereafter to emphasize that it was obtained by atomistic simulations). For more details on the setup and simulation methodology, see Refs. 59 and 60. In the present study, a cylinder of radius 150 Å and length L z ≈ 10 Å was used to study a crack on the (100) crack plane along the crack-front direction [011], which is the natural crack system in tungsten.61

B. Results

1. Simulations of rigid-body separation

The dependence of the binding energy per unit area, i.e., the energy difference of a configuration with respect to the fully separated crystal halves, on the separation distance is shown in Fig. 2 for the two interatomic potentials in comparison to DFT results. From this, the cleavage energy G c is obtained as the difference between the energy at infinite and that at zero separation. The values are summarized in Table II. By comparing with Table I, it can be seen that the cleavage energy approximately equals 2γs(100), which would be the (relaxed) work of separation. The work of separation (WoS) obtained from the DFT calculations (not shown in Table I) equals 8.58 J/m2, which is in good agreement with other DFT values.62,63 The values for the cleavage energies obtained with the EAM potentials are significantly different from the DFT data. Such deviations are not uncommon: see, e.g., Ref. 64 for a review of the applicability of various classes of potentials for bcc transition metals to fracture-like situations. However, for the aim of the present study, it is sufficient to show that the atomistically determined fracture toughness can be reproduced independently of the mesh size of our CZ model. In addition, we will show how the DFT results from the simpler RBS calculations can be used in the same approach, opening up the possibility of a more accurate determination of K Ic without performing the actual DFT crack-tip simulations.

FIG. 2. (a) Binding energy versus separation of crystal halves for two different tungsten potentials and as obtained by DFT. (b) Derivatives of the binding curves, i.e., traction versus separation of crystal halves.

TABLE II. Summary of the characteristic values for rigid-body separations (cleavage energy G c, threshold stress σth, critical separation δcrit, final separation δf) and crack-tip configurations $\left( {K_{{\rm{Ic}}}^{{\rm{at}}}} \right)$ obtained with the ATFS48,49 and Wang50 potentials as well as by DFT calculations.

The TS relationships for the different potentials were readily obtained by applying Eq. (1), i.e., by derivation of the binding-energy versus separation-distance data. The corresponding TS plots are shown in Fig. 2(b). Three characteristic values can be obtained from these plots: the threshold stress, σth, the critical separation distance δcrit (where σth is reached), and the ‘final’ distance δf (where the traction across the surface is zero). All characteristic values are summarized in Table II and compared to values obtained by DFT calculations.

Based on the experience with previous studies on RBS,35,41 the shape of the TS relationships of the Wang potentials appears to be somewhat unexpected as it shows two local maxima and a local minimum at around Δy ≈ 0.65 Å. The origin lies in the not-completely-smooth energy–displacement curves [see Fig. 2(a)]. Both potentials were optimized to reproduce several material properties, but an interplanar separation along the [001] axis was not included. However, as we will see in Sec. III, it is sufficient if the TS law includes the peak stress, and this is reasonably well reproduced by both potentials.

2. Crack-tip simulations

The fracture toughness $\left( {K_{{\rm{Ic}}}^{{\rm{at}}}} \right)$ values determined by atomistic simulations are summarized in Table II. According to the Griffith criterion65 ${K_{\rm{G}}} = \sqrt {2{\gamma _{\rm{s}}}\left( {100} \right) \cdot E}$ (with E being the appropriate orientation-dependent elastic constant), critical values in the range of 1.61–1.64 MPa m1/2 are expected for both potentials. These theoretically expected values using the potential parameters are markedly smaller than the atomistically determined ones; see Sec. II. This well-known behavior66 is generally attributed to the lattice-trapping effect.67 The atomistic crack-tip configurations showing crack propagation by cleavage on the original (100) plane are provided in Fig. 3.

FIG. 3. Atomistic crack-tip configurations in the (100)[011] crack system for the ATFS (a) and Wang (b) potentials at the critical stress intensity factors $K_{{\rm{Ic}}}^{{\rm{at}}}$; see Table II. The original crack-tip atoms are colored black.

3. Comparison of interatomic forces

There are different ways to calculate the force acting across the cleavage plane from the present simulation results: $F_{{\rm{ct}}}^{{\rm{at}}}$ is the y-component of the individual atomic (superscript ‘at’) force vectors calculated during the simulation of crack-containing configurations (subscript ‘ct’). To compare the evolution of the crack-tip forces with the forces during rigid-body separations, the crack-tip distance r ij has to be related to the corresponding separation in the [100] direction, Δy, as follows:

(2)$${\Delta _y} = {{ - {a_0}} / 2} + \sqrt {r_{ij}^2 - {{a_0^2} / 2}} \quad ,$$

where a 0 denotes the bulk lattice parameter. Strictly, this correction is only valid under the assumption that the distance between next-nearest (NN) neighbors (and for cracks on the (100) plane NN neighbor bonds are the critical bonds) changes with Δy, as follows:

(3)$${r_{ij}} = \sqrt {2{{\left( {{{{a_0}} / 2}} \right)}^2} + {{\left( {{{\left( {{{{a_0}} / 2}} \right)}^2} + {\Delta _y}} \right)}^2}} \quad ,$$

i.e., when applying a homogeneous uniaxial strain without allowing Poisson contraction, as is the case in the present study.

$F_{{\rm{rb}}}^{{\rm{at}}}$ is the y-component of the individual atomic force vectors calculated during the simulation of rigid-body separations (subscript ‘rb’). Since there are N surface interacting NN pairs and ${{{N_{{\rm{surface}}}}} \over 2}$ 2nd NN pairs for each (100) surface [Fig. 1(a)], the total force $F_{{\rm{rb}}}^{{\rm{at}}}$ per surface atom is obtained from the individual forces between NN and 2nd NN neighbors, $F_{{\rm{NN}}}^{{\rm{at}}}$ and $F_{2{\rm{NN}}}^{{\rm{at}}}$, by

(4)$$F_{{\rm{rb}}}^{{\rm{at}}} = F_{{\rm{rb,NN}}}^{{\rm{at}}} + {1 / 2}F_{{\rm{rb,2NN}}}^{{\rm{at}}}\quad ,$$

$F_{{\rm{rb}}}^{{\rm{avg}}}$ is calculated from the change of the average potential energy per atom in rigid-body separations:

(5)$$F_{{\rm{rb}}}^{{\rm{avg}}} = {{{N_{{\rm{atoms}}}}} \over {2{N_{{\rm{surface}}}}}}{{{\rm{d}}{E_{{\rm{pot}}}}} \over {{\rm{d}}\Delta }}\quad ,$$

which corresponds to the calculation of the tractions in Eq. (1) but normalized to the number of surface atoms rather than the surface area.

F eff is the derivative of the effective pair potential.68 In this case, the same correction regarding r ij and Δy has to be made as for the crack-tip forces; see Eq. (2).

The evolution of these forces as a function of lateral separation Δy is shown in Figs. 4(a) and 4(b) for the two potentials. The truly remarkable result of this comparison is how well the force laws from the by-far-simpler model systems (rigid-body separations and effective pair potential) agree with the actual forces at the crack tips, considering that the environments in which the forces have been determined (or calculated as in the case of F eff) are completely different.

FIG. 4. Comparison of the relationship between per-atom forces on crack-tip atoms $\left( {F_{{\rm{ct}}}^{{\rm{at}}}} \right)$, the NN and 2NN contributions to the forces on atoms across the cleavage plane during rigid-body separations $\left( {F_{{\rm{rb}}}^{{\rm{at}}}} \right)$, the derivative of the average potential energy during rigid-body separations $\left( {F_{{\rm{rb}}}^{{\rm{avg}}}} \right)$, and the derivative of the effective pair potential (F eff) on the separation distance Δy (a) for the ATFS potential and (b) for the Wang potential.

Second, it should be noted that the unexpected camel-hump shape of the Ty) curve for the Wang potential shown in Fig. 2(b), which is essentially the per-area instead of per-atom versions of $F_{{\rm{rb}}}^{{\rm{avg}}}\left( {{\Delta _y}} \right)$ and thereby $F_{{\rm{rb}}}^{{\rm{at}}}\left( {{\Delta _y}} \right)$, can now be understood by the separation of the total acting force into contributions from NN and 2nd NN neighbors.

Third, the individual force–separation relationships ($F_{{\rm{rb}},{\rm{NN}}}^{{\rm{at}}}$ and $F_{{\rm{rb}},2{\rm{NN}}}^{{\rm{at}}}$) are very similar to each other, and their shape can be qualitatively understood by comparing it to the shape of the derivative of the effective pair potential, F eff. This implicitly shows that the correction made in Eq. (2) is applicable even if it does not account for changes in the electron density due to uniaxial straining, which might, however, be the reason for the minor differences in the magnitudes of $F_{{\rm{rb}},{\rm{NN}}}^{{\rm{at}}}$ and F eff.

An important conclusion which can be drawn from these results is that the forces necessary to break the bond between the crack-tip atoms show only a small—if any at all – dependence on the local atomic surrounding at the crack tip. Even if our observations cannot prove that all EAM potentials will show the same behavior or that crack-tip bonds in metals in general do not show a strong dependence on their local bonding environment, this is to our knowledge the first time that it has been shown that the forces to separate crack-tip atoms can indeed be estimated from simulations of rigid-body displacements. In particular, the comparison between $F_{{\rm{ct}}}^{{\rm{at}}}$ and $F_{{\rm{rb}}}^{{\rm{avg}}}$ in Fig. 4 shows that the characteristics of TS laws derived from simple rigid-body separations, i.e., the maximum force and the critical separation distance δcrit, are in very good agreement with the actual crack-tip forces. This has important consequences for multiscale modeling, as it validates the long-standing use of DFT calculations to determine traction separation laws. Whether the correspondence of forces from rigid-body separation and fracture simulations still holds for systems with complex chemical compositions or grain- and interphase boundaries would still need to be shown. However, an averaging of forces over one structural or chemical unit would be a possible way to extend our approach.


A. Derivation of a scalable TS law

In this section, a method is derived to convert the atomic forces at the crack tip into meaningful continuum quantities, i.e., TS laws describing the cohesive behavior of the material in front of the crack tip. To accomplish this, we consider a chain of N pairs of atoms in front of the crack tip on the crack plane (see Fig. 5). The atom pairs under consideration are those that are affected by the stress concentration of the crack tip and that lie in the crack plane.

FIG. 5. Schematic drawing of the strained atom pairs in the CZ in front of the crack tip.

Formally, the atomic forces F (i), defined in Sec. II.B.2, can be converted into atomic stresses by normalization with the atomic area A at as

(6)$${\sigma _{{\rm{at}}}}\left( {{x^{\left( i \right)}}} \right) = {{{F^{\left( i \right)}}} \over {{A_{{\rm{at}}}}}}\quad ,$$

where x (i) = ia cd is the distance of the atom pair i = 1, 2, …, N to the crack tip and A at = a cfa cd is the projected atomic area on the crack plane, given by the interplanar distance a cf along the crack front and the interplanar distance a cd in crack direction, respectively. The crack itself is defined by the pairs of atoms that are separated by a distance of more than the critical displacement δf, such that the position of the crack tip is located at the last pair of atoms with a larger separation than δf. It is noted here that the elastic distortion of the atomic distances along the crack front and in the crack direction is neglected to obtain a closed-form solution for the scaling relation.

From LEFM, it is known that the stress field in front of a crack tip decays as ${1 / {\sqrt {{x^{\left( i \right)}}} }}$. Figure 6(a) shows that this linear-elastic theory holds excellently at the atomic scale, even very close to the crack tip, in agreement with the findings of Ref. 9. The deviation for the crack-tip atom is partly due to the 20% increased Voronoi volume compared to the atoms in front of the crack tip. In contrast to the atomic stresses, the atomic forces are largest for the crack-tip atom pair. The atomic stresses along the cleavage plane in front of the crack tip can be reasonably expected to follow

(7)$${\sigma _{{\rm{at}}}}\left( {{x^{\left( i \right)}}} \right) = {\sigma _{{\rm{at}}}}\left( {{x^{\left( 1 \right)}}} \right)\sqrt {{{{x^{\left( 1 \right)}}} \over {{x^{\left( i \right)}}}}} = {{{F^{\left( 1 \right)}}} \over {{A_{{\rm{at}}}}}}\sqrt {{{{x^{\left( 1 \right)}}} \over {{x^{\left( i \right)}}}}} = {{{K_{\rm{I}}}} \over {\sqrt {2\pi {x^{\left( i \right)}}} }}\quad ,$$

with the stress intensity factor

(8)$${K_{\rm{I}}} = {{{F^{\left( 1 \right)}}\sqrt {2\pi {a_{{\rm{cd}}}}} } \over {{A_{{\rm{at}}}}}}\quad ,$$

fully defined by atomistic quantities.

FIG. 6. (a) Atomic stresses (black diamonds) in front of the crack tip in the atomistic simulation using the Wang potential at a load of K I = 1.67 MPa m1/2 together with the corresponding analytic LEFM solution, which results in a finite stress at the boundary of our model. The blue points (color online) show the corresponding atom positions (right ordinate axis). (b) Schematic drawing of the stress, with step-wise functions, which indicate the average stress per element in FEM simulations with two different mesh sizes, FE1 and FE2. The region with a significant gradient in the average stress is called the K-dominated region. At the end of the K-dominated zone, the stress gradient decays into a constant far-field stress σ0, which needs to be considered for finite crack sizes; see Eq. (15).

The fundamental relation between the atomistic and the continuum description is demonstrated by calculating the fracture toughness K Ic from Eq. (8), as

(9)$${K_{{\rm{Ic}}}} = {{{F_{\rm{c}}}\sqrt {2\pi {a_{{\rm{cd}}}}} } / {{A_{{\rm{at}}}}}} = \sigma _{\rm{c}}^{{\rm{at}}}\sqrt {2\pi {a_{{\rm{cd}}}}} \quad ,$$

with the critical force F c at which the bond between the first pair of atoms breaks; i.e., this pair of atoms is separated beyond the critical distance δf, and consequently the crack tip advanced by one atomic distance. This force corresponds to the theoretical strength defined in Sec. II.B.1 and can thus be conveniently calculated with atomistic methods. With the values obtained from atomistic calculations with the Wang potential (see Table II), a fracture toughness of K Ic = 1.70 MPa m1/2 is calculated from the critical force, $F_{{\rm{rb}},{\rm{max}}}^{{\rm{avg}}}$, and K Ic = 1.84 MPa m1/2 from the theoretical strength σth with ${a_{{\rm{cd}}}} = {{{a_0}} / {\sqrt 2 }}$. The force-based value corresponds very well to the K Ic values obtained from direct atomistic simulation ($K_{{\rm{Ic}}}^{{\rm{at}}} = 1.68$ MPa m1/2) and from the Griffith criterion ($K_{{\rm{Ic}}}^{\rm{G}} = 1.64$ MPa m1/2). The value from the critical strength is roughly 8% higher as can be expected from the difference in $F_{{\rm{eff}}}^{{\rm{max}}}$ and $F_{{\rm{rb}},{\rm{max}}}^{{\rm{avg}}}$ (see Fig. 4).

However, for most continuum fracture models, a fracture energy or WoS is required in addition to the fracture strength ${\tilde{\sigma }_{\rm{c}}}$ and a critical separation ${\tilde{\delta }_{\rm{c}}}$. These parameters represent the main quantities defining the TS law describing the cohesive behavior of a material in front of a crack tip [see Fig. 6(b)]. While the WoS is scale independent, the fracture strength and the critical separation both depend on the spatial resolution of the applied method. The fundamental material constants are obtained from atomistic methods, but a proper scaling must be developed for the scale-dependent model parameters that takes an appropriate length scale or mesh size into account, as outlined in the introduction.

Since all values used in continuum methods represent average values over a finite volume or a FE, we start by defining the critical fracture strength as average stress over the K I-dominated region

(10)$${\tilde{\sigma }_{\rm{c}}} = {{\sum\limits_{i = 1}^N {{F^{\left( i \right)}}} } \over A}\quad ,$$

as the sum over all atomic forces in this region normalized by the corresponding area. From relation Eq. (7), we also conclude that the atomic forces decay as

(11)$${F^{\left( i \right)}} = {F^{\left( 1 \right)}}\sqrt {{{{x^{\left( 1 \right)}}} \over {{x^{\left( i \right)}}}}} = {{{F^{\left( 1 \right)}}} \over {\sqrt i }}\quad .$$

With Eq. (11), the relation A = NA at and the condition that F (1) = F c, because the onset of crack advance is considered, we obtain

(12)$${\tilde{\sigma }_{\rm{c}}} = {{{F_{\rm{c}}}\sum\limits_{i = 1}^N {{1 \over {\sqrt i }}} } \over {N{A_{{\rm{at}}}}}} = \sigma _{\rm{c}}^{{\rm{at}}}{1 \over N}\sum\limits_{i = 1}^N {{1 \over {\sqrt i }}} \quad .$$

The sum on the right hand side of this equation can be approximated by

(13)$$\tilde{S}\left( N \right) = {1 \over N}\sum\limits_{i = 1}^N {{1 \over {\sqrt i }}} \approx {2 \over {\sqrt N }} - {{1.4} \over N}\quad .$$

Numerical evaluation of the approximation shows that the relative error is less than 10−4 for values of N between 105 and 106; for the latter value, the sum amounts to $\tilde{S}\left( N \right) = 0.002$. The precision of the approximation gets better for larger values of N, which is the usual case for applications. Since the sum is finite, a precise numerical evaluation is always possible, although the simple approximation function is still easier to handle.

To calculate a properly scaled final separation of atoms ${\tilde{\delta }_{\rm{f}}}$, we consider a simple bilinear TS law [see Fig. 7(b)] such that $\Gamma = {1 / {2\sigma _{\rm{c}}^{{\rm{at}}}}}\delta _{\rm{f}}^{{\rm{at}}} = {1 / 2}{\tilde{\sigma }_{\rm{c}}}{\tilde{\delta }_{\rm{f}}}$. It follows

(14)$${\tilde{\delta }_{\rm{f}}} = {{2\Gamma } \over {{{\tilde{\sigma }}_{\rm{c}}}}} = {{\delta _{\rm{f}}^{{\rm{at}}}} \over {\tilde{S}\left( N \right)}}\quad .$$

FIG. 7. (a) Schematic drawing of the region around a crack-opening sampled by elements of different sizes r FE. At each node, the crack opening is characterized by a δ. (b) Bilinear TS law in which critical stress and final displacement depend on the mesh size.

In the final step, a reasonable estimate for the length of the K-dominated region and thus the number N of atom pairs to consider is developed. The criterion for the length of this region is the decay of the stress to a value of σ0, which can be either the yield strength of a plastic material or the far-field stress ${\sigma _0} = {{{K_{\rm{I}}}} / {\sqrt {\pi {a_{{\rm{crack}}}}} }}$ considered in linear-elastic fracture mechanics, where a crack is the crack length. Hence, N is defined by the requirement that

(15)$${{{F^{\left( N \right)}}} \over {{A_{{\rm{at}}}}}} = {{{F^{\left( 1 \right)}}} \over {{A_{{\rm{at}}}}\sqrt N }} = {\sigma _0}\quad .$$

For the critical condition, when the crack starts to propagate (F (1) = F c), we obtain

(16)$${N_{\rm{c}}} = {\left( {{{{F_{\rm{c}}}} \over {{\sigma _0}{A_{{\rm{at}}}}}}} \right)^2} = {\left( {{{\sigma _{\rm{c}}^{{\rm{at}}}} \over {{\sigma _0}}}} \right)^2}\quad .$$

Noting that the critical atomic stresses are typically a 1000 times larger than the yield strength or far-field stress, a reasonable estimate for the number of atom pairs in the K-dominated region is N c ≈ 106. This yields a size of K-dominated zone r coh = N ca cd ≈ 200 μm.

In summary, a scaling relation for the atomistic strengths and separations has been derived, in the form

(17)$${\tilde{\sigma }_{\rm{c}}} = \kappa \sigma _{\rm{c}}^{{\rm{at}}}\quad ,$$
(18)$${\tilde{\delta }_{\rm{f}}} = {{\delta _{\rm{f}}^{{\rm{at}}}} \over \kappa }\quad ,$$

with the scaling factor

(19)$$\kappa = \tilde{S}\left( N \right) = {2 \over {\sqrt N }} - {{1.4} \over N} = 2{{{\sigma _0}} \over {\sigma _{\rm{c}}^{{\rm{at}}}}} - 1.4{\left( {{{{\sigma _0}} \over {\sigma _{\rm{c}}^{{\rm{at}}}}}} \right)^2}\quad .$$

For practical applications, N is replaced by the ratio of the size of each element to the lattice parameter in crack propagation direction, and Eq. (19) becomes

(20)$${\kappa _{{\rm{FE}}}} = 2\sqrt {{{{a_{{\rm{cd}}}}} \over {{r_{{\rm{FE}}}}}}} - 1.4{{{a_{{\rm{cd}}}}} \over {{r_{{\rm{FE}}}}}}\quad .$$

This results in a critical stress and a final displacement, which are both adjusted to the element size, as shown in Fig. 7(b).

It is noted that the directional dependence of the model parameters in Eq. (20) is artificially introduced by the scaling of the FE size with the atomic interplanar distance in the crack propagation direction. A proper directional dependence could be introduced by using a direction-sensitive energy release rate as WoS.

B. Application and validation

As an application of the novel scaling relation, we conducted FE fracture simulations for a linear elastic material representing tungsten to calculate the fracture toughness K Ic. The numerical simulations have been performed using the commercial FE software ABAQUS. All simulations have been performed using an implicit time integration scheme. The surface-based CZ method is used to simulate fracture. The geometry of the model is the same as in the atomistic simulation [see Fig. 1(b)], whereas the ratio between the radius of the sample and the mesh at the crack tip is kept constant at 10,000. This is done to avoid the variation of results due to mesh-size variation in comparison to the sample size. This mesh size is applied to the area around crack tip up to 1% of the total size of the model and then increased gradually toward the outer boundary of the specimen to avoid too much computational cost. The geometry has been meshed with fully integrated 4-node quadratic elements with linear shape functions (CPE4) for the plane strain case. The mechanical boundary conditions have been applied as K I-displacement field in a similar way as in the atomistic simulations. This is achieved by selecting the outermost nodes of the FE model and applying the boundary conditions to those nodes. Isotropic elastic constants, Young’s modulus E and Poisson ratio ν, were determined according to the elastic constants for the Wang potential in Table I, which yields E = 410 GPa and ν = 0.281.

A series of FE simulations has been performed for K I loading using scaled parameters for the TS law according to different mesh sizes, i.e., the scaled critical stress ${\tilde{\sigma }_{\rm{c}}}$, and the displacement at failure ${\tilde{\delta }_{\rm{f}}}$ according to Eq. (20). The value for the initial slope of the TS law has been kept constant at 1.634 × 105 GPa for all cases. The fracture toughness K Ic is defined as the loading at which the first cohesive surface element breaks (δ ≥ δf).

In Fig. 8(a), it can be seen that the values of scaled critical stresses decrease over several orders of magnitude with the increasing mesh size. The values of final displacements (not shown) increases correspondingly to keep the cleavage energy (WoS) constant.

FIG. 8. (a) Scaling of the model parameter of the TS law with the mesh size. (b) Fracture toughness resulting from a FE simulation for a linear elastic material law with the properly scaled TS law (left ordinate axis) and with constant TS law (right ordinate axis) for the cohesive behavior.

Figure 8(b) shows the obtained values of K Ic for different mesh sizes. The value K Ic varies between 1.35 and 1.49 MPa m1/2 over variation of mesh size by 5 orders of magnitude. The average of the K Ic values over the range of mesh sizes is 1.45 MPa m1/2, which is close to the atomistic value $K_{{\rm{Ic}}}^{{\rm{at}}} = 1.68$ MPa m1/2. This is remarkable because there is no fitting of parameters, but rather a consistent scaling of the atomistic quantities according to the mesh size over a range of five orders of magnitude. It can be seen that there is no pronounced scatter around the average value in the sense of a pronounced mesh-size dependence. To emphasize the importance of our scaling procedure, the calculation of K Ic with a constant TS law was also performed. In this case, the scaling is done only for the 10−7 m mesh and the same parameters for the TS law are used for the calculation of K Ic with the 10−6 m mesh (4.59 MPa m1/2) and the 10−4 m mesh (45.82 MPa m1/2). It can be seen in Fig. 8(b) that in this case the critical stress intensity factor scales with the square root of the mesh size (or in other words, the size of the sample), in agreement with the observation of Bazant69 of K Ic variation for self-similar geometries but different sizes.


In the strong stress gradient ahead of the crack tip, the definition of mechanical stress is necessarily scale-dependent. Thus, only fundamental atomistic methods provide a rigorous way to define stress-related material properties, e.g., the fracture strength.

In this work, it is shown that the widely used approach to derive the critical bond separation stress from ab initio calculations of RBSs is well justified. To accomplish this, we have compared forces occurring in atomistic simulations of RBS with those at crack tips under K I loading. The results demonstrate that they are indeed equivalent, which means that they are less environment-dependent than commonly expected. As a consequence, we can derive meaningful TS laws from much simpler model systems than from actual fracture simulations. This enables the use of computationally expensive, but quantitatively reliable methods such as ab initio DFT calculations.

Since the mesh sizes commonly used in continuum methods do not resolve the high atomistic stresses and the strong stress gradient at crack tips, we have introduced a physically meaningful approach to scale atomistic TS laws with the mesh size. It is based on the insight that the atomistically determined stresses ahead of the crack tip agree well with the prediction of LEFM. The approach defined in this work is based on properly averaging atomic forces over cohesive surface elements to convert them in mechanical stresses. This enables the use of fundamental atomistic material constants, such as theoretical strength, theoretical final separation, and work of separation as parameters in continuum fracture models.

It has also been demonstrated that the fracture toughness resulting from the thus-parameterized continuum model is independent of the mesh size, although numerical fluctuations occur. In this way, the atomistic data can be used directly to calculate TS parameters for continuum fracture models, such as CZ models or XFEM.

In our approach, the fracture model describes only the bond-breaking processes during cleavage; any additional (plastic) energy dissipation needs to be taken into account in the constitutive model of the material surrounding the crack. Therefore, the scaling method is ideally used in conjunction with micromechanical plasticity models or discrete dislocation dynamics simulations.

As ab initio tensile tests can be carried out using complex chemical compositions and including grain or interphase boundaries, the scheme enables the material-specific atomistically informed continuum modeling of fracture of engineering materials.


EB gratefully acknowledges the funding from European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme through the project “microKIc—Microscopic Origins of Fracture Toughness” (grant agreement No. 725483). AH and JJM gratefully acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG) through Projects C3 and C4 of the collaborative research center SFB/Transregio 103 superalloy single crystals under grant number INST 213/747-2.


1.Hutchinson, J.W. and Evans, A.G.: Mechanics of materials: Top-down approaches to fracture. Acta Mater. 48, 125 (2000).
2.Needleman, A. and Van der Giessen, E.: Micromechanics of fracture: Connecting physics to engineering. MRS Bull. 26, 211 (2001).
3.Bitzek, E., Kermode, J.R., and Gumbsch, P.: Atomistic aspects of fracture. Int. J. Fract. 191, 13 (2015).
4.Shiari, B. and Miller, R.E.: Multiscale modeling of crack initiation and propagation at the nanoscale. J. Mech. Phys. Solids 88, 35 (2016).
5.Rabczuk, T.: Computational methods for fracture in brittle and quasi-brittle solids: State-of-the-art review and future perspectives. ISRN Appl. Math. 2013, 849231 (2013).
6.Dugdale, D.S.: Yielding of steel sheets containing slits. J. Mech. Phys. Solids 8, 100 (1960).
7.Barenblatt, G.I.: The mathematical theory of equilibrium cracks in brittle fracture. Adv. Appl. Mech. 7, 55 (1962).
8.Rice, J.R.: Mathematical analysis in the mechanics of fracture. In Fracture, Vol. 2, Liebowitz, H., ed. (Academic Press, New York, 1968); chap. 3.
9.Singh, G., Kermode, J.R., De Vita, A., and Zimmerman, R.W.: Validity of linear elasticity in the crack-tip region of ideal brittle solids. Int. J. Fract. 189, 103 (2014).
10.Shimada, T., Ouchi, K., Chihara, Y., and Kitamura, T.: Breakdown of continuum fracture mechanics at the nanoscale. Sci. Rep. 5, 8596 (2015).
11.Zhou, F. and Molinari, J.F.: Dynamic crack propagation with cohesive elements: A methodology to address mesh dependency. Int. J. Numer. Methods Eng. 59, 1 (2004).
12.Turon, A., Dávila, C.G., Camanho, P.P., and Costa, J.: An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models. Eng. Fract. Mech. 74, 1665 (2007).
13.Hillerborg, A., Modéer, M., and Petersson, P-E.: Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem. Concr. Res. 6, 773 (1976).
14.Chandra, N., Li, H., Shet, C., and Ghonem, H.: Some issues in the application of cohesive zone models for metal–ceramic interfaces. Int. J. Solids Struct. 39, 2827 (2002).
15.Spearot, D.E., Jacob, K.I., and McDowell, D.L.: Non-local separation constitutive laws for interfaces and their relation to nanoscale simulations. Mech. Mater. 36, 825 (2004).
16.Yamakov, V., Saether, E., Phillips, D.R., and Glaessgen, E.H.: Molecular-dynamics simulation-based cohesive zone representation of intergranular fracture processes in aluminum. J. Mech. Phys. Solids 54, 1899 (2006).
17.Coffman, V.R., Sethna, J.P., Heber, G., Liu, M., Ingraffea, A., Bailey, N.P., and Barker, E.I.: A comparison of finite element and atomistic modelling of fracture. Modell. Simul. Mater. Sci. Eng. 16, 65008 (2008).
18.Zhou, X.W., Moody, N.R., Jones, R.E., Zimmerman, J.A., and Reedy, E.D.: Molecular-dynamics-based cohesive zone law for brittle interfacial fracture under mixed loading conditions: Effects of elastic constant mismatch. Acta Mater. 57, 4671 (2009).
19.Coffman, V.R., Saetna, J.P., Ingraffea, A.R., Bozek, J.E., Bailey, N.P., Barker, E.I., Sethna, J.P., Ingraffea, A.R., Bozek, J.E., Bailey, N.P., and Barker, E.I.: Challenges in continuum modelling of intergranular fracture. Strain 47, 99 (2011).
20.Lloyd, J.T., Zimmerman, J.A., Jones, R.E., Zhou, X.W., and McDowell, D.L.: Finite element analysis of an atomistically derived cohesive model for brittle fracture. Modell. Simul. Mater. Sci. Eng. 19, 65007 (2011).
21.Barrows, W., Dingreville, R., and Spearot, D.: Traction-separation relationships for hydrogen induced grain boundary embrittlement in nickel via molecular dynamics simulations. Mater. Sci. Eng., A 650, 354 (2016).
22.Kumar, S. and Curtin, W.A.: Crack interaction with microstructure. Mater. Today 10, 34 (2007).
23.Cleveringa, H.H.M., Van der Giessen, E., and Needleman, A.: Discrete dislocation analysis of mode I crack growth. J. Mech. Phys. Solids 48, 1133 (2000).
24.Deshpande, V.S., Needleman, A., and Van der Giessen, E.: Discrete dislocation modeling of fatigue crack propagation. Acta Mater. 50, 831 (2002).
25.Broedling, N.C., Hartmaier, A., and Gao, H.: A combined dislocation—Cohesive zone model for fracture in a confined ductile layer. Int. J. Fract. 140, 169 (2006).
26.Déprés, C., Prasad Reddy, G.V., Robertson, C., and Fivel, M.: An extensive 3D dislocation dynamics investigation of stage-I fatigue crack propagation. Philos. Mag. 94, 4115 (2014).
27.Zhang, P., Klein, P., Huang, Y., Gao, H., and Wu, P.D.: Numerical simulation of cohesive fracture by the virtual-internal-bond model. Comput. Model. Eng. Sci. 3, 263 (2002).
28.Ji, B. and Gao, H.: A study of fracture mechanisms in biological nano-composites via the virtual internal bond model. Mater. Sci. Eng., A 366, 96 (2004).
29.Tahir, A.M., Janisch, R., and Hartmaier, A.: Hydrogen embrittlement of a carbon segregated Σ5 (310) [001] symmetrical tilt grain boundary in α-Fe. Mater. Sci. Eng., A 612, 462 (2014).
30.Ogata, S., Li, J., and Yip, S.: Ideal pure shear strength of aluminum and copper. Science 298, 807 (2002).
31.Černý, M. and Pokluda, J.: The theoretical shear strength of fcc crystals under superimposed triaxial stress. Acta Mater. 58, 3117 (2010).
32.Šob, M., Friák, M., Legut, D., Fiala, J., and Vitek, V.: The role of ab initio electronic structure calculations in studies of the strength of materials. Mater. Sci. Eng., A 387–389, 148 (2004).
33.Ogata, S., Umeno, Y., and Kohyama, M.: First-principles approaches to intrinsic strength and deformation of materials: Perfect crystals, nano-structures, surfaces and interfaces. Modell. Simul. Mater. Sci. Eng. 17, 13001 (2009).
34.Gibson, M.A. and Schuh, C.A.: Segregation-induced changes in grain boundary cohesion and embrittlement in binary alloys. Acta Mater. 95, 145 (2015).
35.Tahir, A.M., Janisch, R., and Hartmaier, A.: Ab initio calculation of traction separation laws for a grain boundary in molybdenum with segregated C impurities. Modell. Simul. Mater. Sci. Eng. 21, 16 (2013).
36.Serebrinsky, S., Carter, E.A., and Ortiz, M.: A quantum-mechanically informed continuum model of hydrogen embrittlement. J. Mech. Phys. Solids 52, 2403 (2004).
37.Nguyen, O. and Ortiz, M.: Coarse-graining and renormalization of atomistic binding relations and universal macroscopic cohesive behavior. J. Mech. Phys. Solids 50, 1727 (2002).
38.Hayes, R.L., Ortiz, M., and Carter, E.A.: Universal binding-energy relation for crystals that accounts for surface relaxation. Phys. Rev. B 69, 1 (2004).
39.Elsner, B.A.M. and Müller, S.: Size effects and strain localization in atomic-scale cleavage modeling. J. Phys.: Condens. Matter 27, 345002 (2015).
40.Enrique, R.A. and Van der Ven, A.: Decohesion models informed by first-principles calculations: The ab initio tensile test. J. Mech. Phys. Solids 107, 494 (2017).
41.Janisch, R., Ahmed, N., and Hartmaier, A.: Ab initio tensile tests of aluminum bulk crystals and grain boundaries: Universality of mechanical behavior. Phys. Rev. B 81, 184108 (2010).
42.Sun, Y.M., Beltz, G.E., and Rice, J.R.: Estimates from atomic models of tension shear coupling in dislocation nucleation from a crack-tip. Mater. Sci. Eng., A 170, 67 (1993).
43.Pang, X.Y., Janisch, R., and Hartmaier, A.: Interplanar potential for tension-shear coupling at grain boundaries derived from ab initio calculations. Modell. Simul. Mater. Sci. Eng. 24, 15007 (2015).
44.Gonze, X., Beuken, J-M., Caracas, R., Detraux, F., Fuchs, M., Rignanese, G-M., Sindic, L., Verstraete, M., Zerah, G., Jollet, F., Torrent, M., Roy, A., Mikami, M., Ghosez, P., Raty, J-Y., and Allan, D.C.: First-principles computation of material properties: The ABINIT software project. Comput. Mater. Sci. 25, 478 (2002).
45.Gonze, X.: A brief introduction to the ABINIT software package. Z. Kristallogr. 220, 558 (2005).
46.Gonze, X., Amadon, B., Anglade, P.M., Beuken, J.M., Bottin, F., Boulanger, P., Bruneval, F., Caliste, D., Caracas, R., Côté, M., Deutsch, T., Genovese, L., Ghosez, P., Giantomassi, M., Goedecker, S., Hamann, D.R., Hermet, P., Jollet, F., Jomard, G., Leroux, S., Mancini, M., Mazevet, S., Oliveira, M.J.T., Onida, G., Pouillon, Y., Rangel, T., Rignanese, G.M., Sangalli, D., Shaltaf, R., Torrent, M., Verstraete, M.J., Zerah, G., and Zwanziger, J.W.: ABINIT: First-principles approach to material and nanosystem properties. Comput. Phys. Commun. 180, 2582 (2009).
47.Hartwigsen, C., Goedecker, S., and Hutter, J.: Relativistic separable dual-space Gaussian pseudopotentials from H to Rn. Phys. Rev. B 58, 3641 (1998).
48.Finnis, M.W. and Sinclair, J.E.: A simple empirical N-body potential for transition metals. Philos. Mag. A 50, 45 (1984).
49.Ackland, G.J. and Thetford, R.: An improved N-body semi-empirical model for body-centred cubic transition metals. Philos. Mag. A 56, 15 (1987).
50.Wang, J., Zhou, Y.L., Li, M., and Hou, Q.: A modified W interatomic potential based on ab initio calculations. Modell. Simul. Mater. Sci. Eng. 22, 15004 (2014).
51.Stadler, J., Mikulla, R., and Trebin, H-R.: IMD: A software package for molecular dynamics studies on parallel computers. Int. J. Mod. Phys. C 8, 1131 (1997).
52.IMD: The ITAP Molecular Dynamics Program. (2003).
53.Beeler, J.R.: Radiation Effects Computer Experiments (North-Holland, Amsterdam, 1983).
54.Bitzek, E., Koskinen, P., Gähler, F., Moseler, M., and Gumbsch, P.: Structural relaxation made simple. Phys. Rev. Lett. 97, 170201 (2006).
55.Černý, M., Šesták, P., Řehák, P., Všianská, M., and Šob, M.: Ab initio tensile tests of grain boundaries in the fcc crystals of Ni and Co with segregated sp-impurities. Mater. Sci. Eng., A 669, 218 (2016).
56.Perez, R. and Gumbsch, P.: An Ab initio study of the cleavage anisotropy in silicon. Acta Mater. 48, 4517 (2000).
57.Sih, G.C., Paris, P.C., and Irwin, G.R.: On cracks in rectilinearly anisotropic bodies. Int. J. Fract. 1, 189 (1965).
58.Sih, G.C. and Liebowitz, H.: Mathematical theories of brittle fracture. In Fracture, Vol. 2, Liebowitz, H., ed. (Academic Press, New York, 1968); chap. 2.
59.Möller, J.J. and Bitzek, E.: Comparative study of embedded atom potentials for atomistic simulations of fracture in α-iron. Modell. Simul. Mater. Sci. Eng. 22, 045002 (2014).
60.Möller, J.J. and Bitzek, E.: Fracture toughness and bond trapping of grain boundary cracks. Acta Mater. 73, 1 (2014).
61.Gumbsch, P., Riedle, J., Hartmaier, A., and Fischmeister, H.F.: Controlling factors for the brittle-to-ductile transition in tungsten single crystals. Science 282, 1293 (1998).
62.Vitos, L., Ruban, A.V., Skriver, H.L., and Kollar, J.: The surface energy of metals. Surf. Sci. 411, 186 (1998).
63.Piazza, Z.A., Ajmalghan, M., Ferro, Y., and Kolasinski, R.D.: Saturation of tungsten surfaces with hydrogen: A density functional theory study complemented by low energy ion scattering and direct recoil spectroscopy data. Acta Mater. 145, 388 (2018).
64.Möller, J.J., Mrovec, M., Bleskov, I., Neugebauer, J., Hammerschmidt, T., Drautz, R., Elsässer, C., Hickel, T., and Bitzek, E.: On {110} planar faults in strained bcc metals: Origins and implications of a commonly observed artifact of classical potentials. Phys. Rev. Mater. 2, 093606 (2018).
65.Griffith, A.A.: The phenomena of rupture and flow in solids. Philos. Trans. R. Soc., A 221, 163 (1921).
66.Gumbsch, P.: Atomistische Modellierung zweidimensionaler Defekte in Metallen: Risse, Phasengrenzflächen. PhD Dissertation, Universität Stuttgart, Stuttgart, 1991.
67.Thomson, R., Hsieh, C., and Rana, V.: Lattice trapping of fracture cracks. J. Appl. Phys. 42, 3154 (1971).
68.Johnson, R.A.: Alloy models with the embedded-atom method. Phys. Rev. B 39, 12554 (1989).
69.Bažant, Z.P., ed.: Fracture mechanics of concrete structures. In Proc. FraMCoS1-lnt. Conf. On Fracture Mechanics of Concrete Structures (Beaver Run Resort, Breckenridge, Colorado, 1992); pp. 623.