## Abstract

In this paper, we discuss the formulation, recent developments and findings obtained from a mesoscale mechanics technique called phase field dislocation dynamics (PFDD). We begin by presenting recent advancements made in modelling face-centred cubic materials, such as integration with atomic-scale simulations to account for partial dislocations. We discuss calculations that help in understanding grain size effects on transitions from full to partial dislocation-mediated slip behaviour and deformation twinning. Finally, we present recent extensions of the PFDD framework to alternative crystal structures, such as body-centred cubic metals, and two-phase materials, including free surfaces, voids and bi-metallic crystals. With several examples we demonstrate that the PFDD model is a powerful and versatile method that can bridge the length and time scales between atomistic and continuum-scale methods, providing a much needed understanding of deformation mechanisms in the mesoscale regime.

## 1. Introduction

Over several decades, the dynamics of dislocations in metallic crystals or polycrystals have been examined using a variety of mechanics-based simulation tools that differ in length and time scales. At the lowest of these scales, atomic-scale simulation methods have allowed for close study of one to a few dislocations. However, computational limitations prevent such methods from achieving length and time scales associated with basic experimental mechanical tests. At the other end of the spectrum are crystal plasticity simulations, where the plastic deformation of single crystals and polycrystalline materials is realized by the activation of a number of crystallographic slip systems. Dislocations are appreciated by the crystallographic shear they enable. Length scales are assumed sufficiently large such that it is not necessary to model dislocations discretely. Consequently, much of the dislocation processes occurring within strained crystals is homogenized, and the outcomes of dislocation motion are reflected in rate laws for the evolution of a homogeneous density of dislocations.

Thus, it has long been realized that dislocation-based phenomena at an intermediate time and length scale regime, called the mesoscale, are critical to understanding and predicting the response of metals to mechanical loads. At this scale, plasticity is realized as a collection of individual dislocations, closely interacting with each other and internal material boundaries, such as free surfaces, grain boundaries (GBs) and heterophase interfaces. The few simulation methods developed to operate at this scale have the potential to take essential atomic-driven behaviour seen at lower length scales and link it to the crystallographic plastic slip taken into account at the macroscale.

In this article, we present the development and findings of a mesoscale mechanics technique called phase field dislocation dynamics (PFDD). PFDD determines the motion and configuration of individual dislocations by minimizing the total system energy. Interactions between dislocations are based on continuum linear elastic dislocation theory. Hence, atomic motions and interactions are not explicitly calculated and much larger crystal sizes and longer time scales (of the order of seconds) can be assessed. In efforts to recover some atomic-scale physics, PFDD has been extended to permit incorporation of the energetics of atomic-scale dislocation cores from density functional theory (DFT) calculations [1,2], which is particularly important for partial dislocation motion and interactions. Most recent model advancements include additional energetics accounting for differences in moduli in two-phase materials and dislocation interactions with bi-phase interfaces [3].

We begin by reviewing the basics of the PFDD formulation and integration of results from atomic-scale calculations. To demonstrate its capability, we will first focus on modelling extended dislocation configurations in face-centred cubic (fcc) metals. We then discuss a number of problems and key findings related to understanding the variables that govern transitions from full dislocation slip to partial dislocation slip and deformation twinning in nano-sized fcc crystals. Such studies help define the dominant deformation mechanisms active for a given material and loading condition and determine the impact of such transitions on the macroscopic stress–strain response. We continue discussion with problems addressing alternative crystal structures (body-centred cubic (bcc) metals) and two-phase materials, which includes bi-metallic crystals as well as free surfaces and voids. When possible, we compare the results with experimental data reported in the literature. The PFDD model has great potential to be a powerful method to provide an understanding of deformation mechanisms that are needed to inform physically based continuum models.

## 2. Phase field dislocation dynamics

In general, phase field frameworks describe physical behaviour by tracking and evolving one or more scalar order parameters. These order parameters are related to some quantity of interest, which, in the case of PFDD, are line defects in the crystal lattice. In the phase field framework, the system is then evolved through the minimization of the total energy. All physics are, therefore, described using an energy functional that is dependent on the order parameters.

The PFDD model assumes that plasticity is mediated by the motion and interaction of dislocations. In this section, we limit our discussion of the PFDD framework to perfect dislocations in order to introduce the overall framework. The model extensions required to account for partial dislocation motion and interactions are discussed later in §3.

Considering only perfect dislocations, the phase field variables (or order parameters), *ζ*_{α}(**x**,*t*), in the PFDD model track the number and sign of the dislocations that glide on each active slip system *α*. For instance, for most fcc metals, there are 12 independently oriented slip systems that belong to the family {111}〈110〉. Hence, there can be up to 12 different phase field variables active in an fcc unit cell. Extensions for other crystal structures, such as bcc, are addressed later in §6.

Because we assume that plasticity is mediated by the motion and interaction of dislocations, the plastic strain, *ϵ*^{p}, will be directly proportional to the number of gliding dislocations. Hence, the plastic strain can be expressed in terms of the order parameters and can be written as [1,3,4]
2.1
where *N* is the number of slip systems, *b* is the magnitude of the Burgers vector, *s*^{α} is the slip direction, *m*^{α} is the slip plane normal associated with slip system *α*, and *δ*_{n} is a Dirac distribution supported on the slip plane, *n*.

As mentioned above, the order parameters (and also equation (2.1)) are used in the formulation for the energy functional that describes the needed physics for system evolution. In the PFDD model, we have two key energy terms [1,4,5]:
2.2
The first term, *E*^{strain}, is the strain energy associated with the system. This term accounts for both short- and long-range interactions between dislocations, in addition to interactions between an applied external stress and the dislocations in the system. The second term is the core energy, *E*^{core}, which describes the energy barrier encountered as the dislocation moves through the crystal lattice (i.e., the Peierls barrier).

The strain energy stored within the system can be expressed as an additive decomposition of internal and external interaction terms:
2.3
where *C*_{ijkl} is the elastic moduli tensor, *ϵ*^{e} is the elastic strain and *σ*^{appl}_{ij} is the externally applied stress. To evaluate the elastic strain in the first term in equation (2.3), it is assumed that the total strain can be additively decomposed into elastic and plastic components . Using stress equilibrium (*σ*_{ij,j}=0), the total strain can be solved in terms of . This translation is greatly simplified by transforming into Fourier space, where the internal strain energy can be expressed as [4–6]
2.4
where , a superposed (^{∧}) denotes the Fourier transform, the superscript symbol (*) indicates complex conjugation, *k*_{i} is the wavenumber vector, is the Fourier transform of the Green's tensor of linear elasticity [7], and the denotes the principal value of the integral. With this, the strain energy can now be written in terms of the phase field variables as
2.5

Thus far, the energy terms presented describe how dislocations will be impacted by their environment, whether internally (via dislocation–dislocation interactions) or externally (through the applied load). What has not yet been taken into account is the energy expended through the motion of a single dislocation through the crystal lattice by breaking and reforming atomic bonds. Because the phase field formulation does not explicitly account for individual atoms, the energetic contribution due to the motion of the dislocation core is described using a material-dependent energy landscape, called the *γ*-surface, which can be calculated with atomistic methods. For perfect dislocations, the core energy *E*^{core} can be modelled by the Peierls potential, which exhibits a periodicity that follows that of the crystal lattice. It can be expressed appropriately as the integral of a Fourier sine series, because this series mimics the sinusoidal nature of the Peierls potential in a regular cubic lattice [5,8]:
2.6
where *B* defines the magnitude of the energy barrier that must be overcome to activate slip. Accordingly, integer values of the phase field parameter represent perfect Burgers vector translations and correspond to zero core energy. These regions represent places that have been slipped by one or more dislocations (or not slipped at all in the case that the order parameter equals zero), hence the atomic bonds have already been broken and reformed to propagate any dislocations. Conversely, at the dislocation line where the dislocation core lies, the phase field parameter transitions from zero to one. This region represent areas where atomic bonds are being broken and reformed. It is also where the core energy is non-zero and directly related to *B*.

In the PFDD formulation, the dislocation dynamics are governed by the minimization of the total energy, which is dictated by a time-dependent Ginzburg–Landau (TDGL) kinetic equation. This equation relates the change in the phase field variables in time to the variation in the total system energy with respect to the phase field variables [1,4,5] as follows:
2.7
where *L* is a kinetic coefficient that defines the time scale of the simulation. Once the total energy from equation (2.2) is substituted into equation (2.7), the variational derivative with respect to the phase field, *ζ*_{α}(*x*,*t*), can be taken. This action produces a set (one for each active slip system) of coupled integro-differential equations that describe the dislocation dynamics for the entire system. The solution to these equations provides the location of each dislocation in the system.

For the simulations discussed in the sections below, we assume quasi-static loading. Hence, we are mainly concerned with the equilibrium state at the end of each time step, even though a time-dependent equation is used for system evolution. Equilibrium is reached when the left-hand side of equation (2.7) becomes zero. Thus, the equilibrium state is independent of the parameter *L*.

## 3. Extended dislocations

In the foregoing §2, the PFDD formulation for perfect dislocations belonging to the {111}〈110〉 systems was presented. In this section, extensions required to account for the motion and interactions of Shockley partial dislocations belonging to the {111}〈112〉 systems are discussed. Because the phase field variables are defined in terms of perfect dislocation slip directions, linear combinations of these order parameters must be used in order to account for partials [1,6]. Hence, at a minimum, in order to model partial dislocations three phase field variables must be active on the same glide plane (for fcc metals). These partials have a smaller Burgers vectors than full dislocations. Thus, as another consequence, the phase field variables associated with the passage of partial dislocations are non-integers. When leading and trailing partials cross the same glide plane, the combined displacement corresponds to that of a full dislocation and the appropriate phase field parameter will achieve an integer value, whereas the other two active order parameters return to zero. Equation (2.5) will still calculate the appropriate interactions between partial dislocation lines and between the partial dislocations and the externally applied stress.

Conversely, for partial dislocations, equation (2.6) no longer applies. In fact, to accurately model the stacking fault region bounded by leading and trailing partial dislocations, a better approximation of the material *γ*-surface must be included. To this end, we replace *E*^{core} above, with the generalized stacking fault energy (GSFE), *E*^{gsfe}. Similar to *E*^{core}, *E*^{gsfe} also accounts for the energy required to move the dislocation core through the crystal lattice [9–12]. However, to enable the formation, motion, and interaction of partial dislocations in an fcc metal, this term is formulated to account for the entire three-dimensional fcc material *γ*-surface.

For simplicity and without loss of generality, we reduce the problem to treat dislocations gliding on a single glide plane so that *E*^{gsfe} becomes an explicit function of all three phase variables *ζ*_{1},*ζ*_{2}, and *ζ*_{3} belonging to this plane. *E*^{gsfe} is expressed with the following complex Fourier series [1,9,13,14]:
3.1
The coefficients *c*_{0} to *c*_{4} and *a*_{1},*a*_{3} can be determined from either DFT or molecular dynamics (MD) simulations and are discussed in more detail in [2]. They are defined as [13]
3.2
where *T*,*T*_{1},*G*,*G*_{1},*G*_{2},*G*_{3} are certain energetic maxima and minima taken directly from the *γ*-surface. For example, points *G* and *G*_{2} represent the intrinsic and unstable SFEs, respectively. In the present calculations, these points will be taken directly from the *γ*-surface or GSFE curves generated via atomistic methods [1,2,9,13,14,]. Figure 1*a* shows the GSFE curves for a wide range of fcc metals as calculated with DFT [2]. Figure 1*b* presents a parametrized *γ*-surface for Ni as calculated with equation (3.1). The points marked on this *γ*-surface indicate the energy values taken from the GSFE curves and used to calculated the parameters in equations (3.2).

The total energy shown as equation (2.2) can be rewritten as 3.3 This new equation is then substituted into equation (2.7) to evolve the system, which now includes partial dislocations and extended core regions.

In the following sections, examples are presented that include calculations on the evolution of partial dislocations in fcc metals. For these, we employ *ab initio* DFT to calculate the parameters needed to generate the full material *γ*-surface. For consistency, these calculations for all metals examined in this work were carried out with the Vienna *ab initio* simulation package [15,16] using the projector augmented wave method within the Perdew–Burke–Ernzerhof (PBE) approximation [17,2]. For the elastic moduli of all metals examined here, we also use the results from DFT calculations with the PBE approximation, which are reported in Hunter *et al.* [2].

### (a) Case study: stacking fault widths in face-centred cubic metals

In the PFDD simulation, we model a periodic cell, which is a block of material containing an infinitely long perfect dislocation with either edge or screw orientation. Without any externally applied stress, the perfect dislocation will dissociate into leading and trailing Shockley partial dislocations bounding an intrinsic stacking fault region. The partials will spread apart until the repulsion between two like-signed dislocations is balanced by the energy penalty required to expand the stacking fault. Once an equilibrium configuration is reached, the distance between the leading and trailing partials is the equilibrium stacking fault width, *w*_{0}. Because the *γ*-surface varies widely from material to material (i.e. figure 1*a*), the energy penalty to expand the stacking fault is large in some materials and small in others, which correlates to small and large equilibrium core widths, respectively.

Equilibrium stacking fault widths, *w*_{0}, have been calculated for nine fcc metals, namely Ag, Al, Au, Cu, Ni, Ir, Pd, Pt and Rh using the PFDD model with the extension to partial dislocations [2] and are presented in table 1. The objective of these calculations is to provide insights into the material parameters, particularly those from the material *γ*-surface, that control dislocation core properties. Over the years, defining the relationships between core properties and material parameters that apply to a wide range of material properties has been difficult. Experimental equilibrium stacking fault widths exhibit large variations owing to dislocation line curvature, orientation and nearby unknown internal stress fields [18–20]. The classical model from dislocation theory for stacking fault widths does not consistently fit the experimental data for all metals [9,21]. MD is the current numerical choice for core calculations; however, results are dependent on the interatomic potential used in the calculations [9,22–25].

Figure 2*a* shows the relationship between *w*_{0}, as calculated with PFDD, in relation to the normalized intrinsic SFE *γ*_{I}/*μb*. Clearly, *w*_{0} exhibits a strong inverse relationship with *γ*_{I}/*μb*. There is, however, one notable outlier, Pt. To obtain a deeper understanding of the core structure of Pt, we use a metric Δ to better visualize the partial dislocations. Δ is the vector projection of crystallographic slip displacements on the glide plane in the direction of the initial edge dislocation. By this metric, Δ=0 corresponds to no displacement, Δ=1 represents the passage of a perfect dislocation, and Δ=0.5 represents a compact Shockley partial. Other values of Δ are associated with fractional dislocations.

Figure 2*b* shows the dΔ/d*x* profile for Pt and Al. These metals were selected for comparison, because both metals have essentially the same normalized intrinsic SFE, yet their intrinsic SFEs correspond to very different depths in the *γ*-surface. In other words, the difference Δ*γ*=*γ*_{U}−*γ*_{I} is not the same, being notably very small in Pt relative to that in Al. The profile in figure 2*b* for Al shows two well-defined peaks in dΔ/d*x*, each representing a partial dislocation. The distance between these two peaks defines the width of the stacking fault region. Conversely, the profile for Pt shows that the core region is ‘smeared’, or comprised of a distribution of several small partial dislocations and no clear intrinsic stacking fault region. The PFDD simulation suggests that unlike in Al, in Pt, a diffuse distribution of defects spread throughout the core region is the energetically favoured equilibrium structure. Therefore, the core of the Pt dislocation does not exhibit the classical picture of an extended dislocation core comprised of two distinct Shockley partials bounding a well-defined stacking fault region. This may explain why Pt exhibits a wider equilibrium stacking fault width than expected, as shown in figure 2*a*. Most interestingly, it implies a relationship between the depth of the energy well in the *γ*-surface and the defect configuration of the extended dislocation core.

## 4. Intragranular dislocation mechanics

In coarse-grained metals, dislocations can nucleate from both the grain interior and the grain boundaries (GBs). As grain sizes in polycrystals decrease to nanoscale dimensions they nucleate primarily from the GBs. First, the leading partial emits and extends into the grain forming a stacking fault. If the trailing partial subsequently emits after the leading partial, then the stacking fault is erased and the trailing and leading partials combine into a full dislocation. The crystal is considered to deform by full-dislocation-mediated slip. If the trailing partial does not emit and predominantly leading partials accommodate slip in the grain, then the crystal is said to deform alternatively by partial-mediated slip. Whether partial-mediated slip or conventional slip operates greatly affects the propensity for basic deformation mechanisms, such as twinning, cross-slip, dislocation storage (strain hardening) and GB mobility (grain growth) [18,20,26–28].

To date, the key variables that determine the transition from partial- to full-dislocation mediated slip have yet to be clarified. Experimental and atomistic studies have observed that grain size and SFE can influence the crossover from partial-mediated to full-dislocation-mediated slip [24,26,29–33]. As the PFDD model can treat both partial and full dislocation slip in the same simulation, it is well suited to uncover the influential variables governing this transition. Below, we apply the PFDD formulation with DFT-calculated *γ*-surfaces to study single-loop nucleation and expansion within a nano-sized fcc grain. Later in §5, we discuss an extension of this study to a polycrystal.

### (a) Case study: nucleation from grain boundaries

We model a grain as a finite-sized, cubic grain as shown in figure 3. The side-length of the cube corresponds to the grain size *D*. The lattice orientation of the cube is , , and [111] in the *x*,*y* and *z* directions, respectively. We employ periodic boundary conditions in all three directions. Interactions between the central dislocation and its image are expected to affect the dislocation speed as it nears the opposing boundary and for this reason, dislocation speed is not discussed. The GBs are modelled as an array of pinned obstacles of width *w*∼1.25 nm. In one GB, a ledge with a given length *L* and step height *w*_{L} is placed. The effective Burgers vector of this GB defect is equal to that of a full edge dislocation. Points within the boundary and ledge are not evolved using the TDGL equation and hence not allowed to change value and have zero mobility. The crystal is then deformed under a shear stress *σ*_{xz}. It is increased at a quasi-static rate until a critical value is reached, sufficient to just produce a Shockley partial from the ledge. Once this threshold value *σ*_{th} is defined, it is kept constant throughout the entire simulation.

Figure 4 shows the evolution of a loop under application of *σ*_{th}. The material is Ni and the grain size is *D* = 31.87 nm. In the first snapshot, a leading Shockley partial dislocation loop has emitted from the GB ledge into the crystal. It has expanded some distance into the grain creating a stacking fault, marked by the green area. In the next frame, the stacking fault has become large enough that it is no longer energetically favourable to continue to expand the fault region. Rather the trailing Shockley partial emits from the grain boundary defect. Just prior to trailing partial emission, the stacking fault width is at its maximum. We will refer to this critical stacking fault width as . Once the trailing partial has been emitted, it glides at a significantly faster speed than the leading partial to remove the faulted region. The two Shockley partials maintain a separation nearly equal to *w*_{0} as they glide together across the grain. The two partials bow out in different directions owing to their differing Burgers vector and assume an asymmetric shape in order to maximize the length of the dislocation with screw character, because it has the lower self-energy.

Next, the effects of *L* and *D* on *σ*_{th} are systematically studied. For this purpose, the range of *D* examined is made sufficiently large, so that local fluctuations in stress at the GBs do not influence the results. Figure 5*a* shows *σ*_{th} as a function of *D* for Ni, Pd, Al and Au, but with varying ledge sizes. The results show that the variation with *D* is negligible. However, for a given *L*, *σ*_{th} can be related to the nucleation stress. Prior works have suggested that the nucleation of a Shockley partial is related to *γ*_{U}, whereas others to the difference Δ*γ*=*γ*_{U}−*γ*_{I} [34–36]. With this broad range of metals, we also examine the effects of the *γ*-surface on *σ*_{th}. Figure 5*b* plots *σ*_{th} versus *γ*_{U}. We observe a strong linear scaling. The same is not found for Δ*γ*, which is seen if one compares Al and Au. Both *σ*_{th} and *γ*_{U} are larger for Al than Au, but they have nearly the same Δ*γ*. These results indicate that *γ*_{U}, and not Δ*γ*, is the relevant energetic barrier for nucleating a partial dislocation.

### (b) Case study: grain size effects on single loop behaviour

The extent to which the leading partial can glide from its source in the GB and into the grain before the trailing partial is emitted is the maximum width of the stacking fault (see the second frame in figure 4). The value of relative to *D* has been used as a measure of the transition from full to partial dislocation-mediated plasticity [26,32].

To explore the effect of grain size *D* on , we repeated the Ni simulation in figure 4 for larger grain sizes using the same *σ*_{th}. As shown in figure 6, we see that increases as *D* increases. Eventually, it saturates to a value of ∼20 nm at larger *D*. For much smaller grain sizes, such as *D* = 7.97 nm, the partial emerges and is able to glide across the grain without emission of the trailing partial (hence this is not shown in figure 6), forming a monolayer stacking fault that spans the entire grain cross section. The same occurs for *D*<7.97 nm but not larger; therefore, the maximum grain size for full traversal of the leading partial, *D*_{SF}, is 7.97 nm, for Ni.

Taken together, we find that there are two apparent bounds on the variation of the stacking fault width with grain size, depicted in figure 7. As grain size *D* increases, eventually saturates, corresponding to an upper bound on . As *D* decreases, eventually coincides with *D*. The maximum grain size *D*_{SF}, below which the stacking fault covers the grain cross section is an apparent lower bound.

Using the same set-up, we vary the fcc material studied in order to examine the effects of SFEs, key variables expected to influence . In these simulations, the ledge size is fixed at 5*b* and therefore, for each material, the *σ*_{th} is fixed, whereas *D* is varied. The results suggest a strong dependency of on the intrinsic SFE *γ*_{I}. While this is not surprising, we do find it intriguing that the scaling is nonlinear, such that small differences in *γ*_{I} lead to large differences in . The normalized *γ*_{I}/*μb* are very similar for Ni and Au, yet there is at least an order of magnitude difference in their . Likewise, *γ*_{I}/*μb* for Cu is slightly smaller than that for Ni, yet is more than 10 times larger than that for Ni.

To summarize, the transition from partial to full dislocation-mediated slip can be associated with grain sizes at or below *D*_{SF}. This critical grain size decreases as the intrinsic SFE increases and notably is not influenced by the unstable SFE. The unstable SFE, however, does affect the stress to nucleate the partial dislocations from the boundary more than the extent to which they glide thereafter.

### (c) Case study: deformation twinning

In fcc metals, the minimum number of layers for a stable nascent twin has been calculated to vary from two to three [38,39]. However, the kinetic pathway via which two or three adjacent faults can construct the twin is unknown. Proposals made over the years can be classified as either monotonic activation of partials (MAP) or random activation of partials (RAPs). In MAP, a twin is built by sequential glide of twinning dislocations with like Burgers vectors on adjacent planes [33,40]. Consequently, this pathway leaves ever increasing steps in the GB, and unless the boundary can reconstruct itself, backstresses will develop as the twin expands. In RAP, which can also be called a zero-strain twinning mechanism, the twinning dislocations also glide on adjacent planes, but have different Burgers vectors that sum to zero [41]. This mechanism leaves no possibility for backstresses to develop, but the twin does not accommodate strain.

In both MAP and RAP, twins form from partial dislocations emitted from the same GB. In contrast, a third mechanism has been proposed, called alternate emission (AE) in which twins are built from partials emitting from GBs located on opposing sides of the grain. This mechanism has been seen in atomic-scale simulation and indirectly in *in situ* transmission electron microscopy [42]. PFDD reveals that this AE mechanism can be just as, if not more, energetically favourable than MAP and RAP [6,43]. Upon applying the threshold stress *σ*_{th}, each ledge emits a Shockley partial. They are of opposite sign and they glide towards one another. They eventually traverse the grain and form a stable fault. The final structure is a two-layer twin.

Observations from experiment suggest that the propensity for twinning in fcc nanocrystalline materials depends, at least, on grain size and fault energies [31,44–47]. For twinning, the grain size must be sufficiently small for partial-mediated slip. In this model, it would appear that the collaborative effect of these two criteria correspond to *D*<*D*_{SF}, the maximum grain size for which a leading partial nucleated from the GB can traverse the grain without nucleating the trailing partial. As seen earlier, *D*_{SF} depends on fault energies, generally being larger for smaller intrinsic SFEs. Therefore, we can expect that smaller grain sizes and lower *γ*_{I} promote twinning.

To test this idea, the initial configuration shown in figure 3 is extended to include another GB ledge on the opposing GB and an adjacent plane. Figure 8 shows the formation of a two-layer twin via the AE mechanism in a 37.6 nm grain in Ag. This two-ledge configuration was deformed for many other fcc metals (Al, Au, Cu, Ir, Ni, Pd, Rh) for *D*<*D*_{SF}. Interestingly, we find that not all *D*<*D*_{SF} cases are able to twin via the AE mechanism. Only Ag, Cu and Rh developed twins via the AE mechanism, whereas the rest of them did not. Therefore, another criterion, apart from a sufficiently small grain size, *D*<*D*_{SF}, and suitably low *γ*_{I}, must affect the propensity for twinning.

To expose it, we examine closely Au, a non-twinning case (also shown in figure 8), because it has nearly the same low *γ*_{I} as Ag, which does twin. These two materials also have very similar shear moduli and lattice parameters, but notably different *γ*_{U}. We begin with the same crystal with two GB ledges on opposing sides. Under the *σ*_{th} for Au, we see that as before in Ag, Shockley partials of opposite sign are emitted from each ledge and driven towards one another. Once the partial loops get close enough to interact, they expand so as to maximize their screw components while continuing to accommodate strain. When this action eventually becomes energetically unfavourable, a trailing partial is emitted from either one or both ledges. Perfect dislocations form as a result. They glide across the grain, accommodating strain. Thus, deformation is accommodated by full-dislocation slip and not twinning. This is the case in Au in grain sizes as small at 10 nm.

Evidently, for AE twinning, the activation barrier for trail nucleation must be larger than the attractive interaction energy between two approaching, oppositely signed partials. This criterion could be translated to mean that the energy difference for nucleation of a trail over creating the stacking fault, i.e. Δ*γ*=(*γ*_{U}−*γ*_{I})/*μb* must be sufficiently high. Considering all criteria discussed thus far, twinning would depend on both the intrinsic and unstable SFEs, becoming more likely as *γ*_{I}/*μb* decreases and Δ*γ* increases. The map in figure 9 plots Δ*γ* versus *γ*_{I}/*μb* and locates all the fcc metals studied here with respect to their twinning likelihood, provided that their grain size is sufficiently small for partial-mediated slip. Twinning is more probable for metals lying in the top left hand corner, and as we see, Ag, Cu and Rh lie in this region and are separated from the others.

## 5. Polycrystalline dislocation mechanics

Thus far, we have applied PFDD to study dislocation motion within a strained crystal. PFDD can also be used to study similar phenomenon in nanocrystalline materials. In this section, we discuss recent results for polycrystals found through application of the PFDD model.

### (a) Case study: transition from full dislocation to partial dislocation slip

In recent work, PFDD with DFT-calculated *γ*-surfaces was employed to study the nucleation of partial and full dislocations from GBs in nanocrystalline structures [48]. Two materials were compared, Ni and Al, with the aim of elucidating the effects of SFEs (figure 1). The material properties used in simulation are summarized in [48].

The basic nanocrystalline structure used in the PFDD simulations is shown in figure 10. This grain structure was constructed from atomic coordinates using a Voronoi algorithm [49]. More specifically, the coordination number of each atom is calculated and non-coordinated atoms are regarded as GBs. The GBs in the atomic-scale model are then mapped onto grids that are then used in the PFDD simulations. These GBs are coloured in black in figure 10. As shown, they have a finite volume. Each GB is assigned an initial dislocation structure represented by a set of GB dislocations. In this way, under stress, dislocations will nucleate from these GB defects.

As in all simulations presented thus far in this paper, the grid size is the Burgers vector. Simulations were repeated for nanocrystalline structures ranging from 10 to 50 nm. An initial dislocation structure is assigned to each GB. When mobile dislocations in the grains interact with these GB dislocations, the GBs can emit new dislocations, annihilate absorbed dislocations, transmit the dislocation across the grain boundary or stop its motion.

In deformation, in both materials and all grain sizes, the GBs nucleate leading partial dislocations, which are then followed by trailing partials. To explore transitions to and from partial- and full-dislocation-mediated slip, we study the evolution of areas slipped by partial dislocations, *A*_{p}, and by extended full dislocations, *A*_{f}. The areas created by the glide of leading partial dislocations without the trailing partial count towards *A*_{p}. After the trailing partials nucleate, the stacking fault is removed as *A*_{p} and counts towards *A*_{f}. If the trailing partial nucleates immediately after the leading partial, then an extended full dislocation with a stacking fault ribbon smaller than the grain size is formed. The area slipped by this full dislocation is also added to *A*_{f}.

Figure 11*a*,*b* shows the strain evolution of *A*_{p} and *A*_{f} for a 46 nm Al grain and a 40 nm grain size in Ni during quasi-static deformation. Figure 11*c*,*d* shows an example of the dislocation evolution in the Al and Ni grains, respectively. In Ni, *A*_{p} stabilizes, because a balance between the addition and subtraction of faults is reached. In contrast, in Al, the amount of *A*_{f} greatly exceeds *A*_{p} even at small strains. This is not surprising as Al has the higher intrinsic SFE and a relative lower unstable SFE (figure 1). Thus, compared with Ni, stacking faults in Al have a higher energetic penalty and trailing partials in Al have a lower nucleation barrier.

Two more interesting outcomes emerged from these nanocrystalline simulations. First, we see that the amount of *A*_{p} increased as the grain size decreased. However, the crossover from full to partial dislocation-mediated slip with a reduction in grain size is not abrupt. Second, the contribution of partial and full dislocation-mediated slip is not constant with strain. As strain increases, the involvement of full dislocation-mediated slip tends to increase.

## 6. Alternative crystal structures

The PFDD model was originally developed to simulate the evolution of dislocations in fcc materials. Fundamentally, fcc metals deform by one slip mode, unlike bcc metals and many low-symmetry metals that deform by multiple slip modes. Moreover, the motion of dislocations in other metals can be more complex than those in fcc metals [50–53]. In this section, we present an extension of the PFDD framework to bcc crystal structures.

### (a) Case study: loop expansion in body-centred cubic metals

Dislocations in bcc metals can glide using three slip modes: {110}〈111〉, {112}〈111〉 and {123}〈111〉. Prior atomic scale calculations find that these three modes have different core structures and Peierls barriers [54–59], meaning that the relative values of the characteristic stress required to activate dislocation motion can be different for each of these three slip modes. The relative differences in these critical stresses can greatly influence macroscopic flow stress and texture evolution [60–62].

Accounting for multiple slip modes within the phase field framework is relatively straightforward. The phase field variables themselves are scalar and carry no direction or orientation. The orientations of the different planes are accounted for through the Burgers vectors and slip plane normals included in the description of the plastic strain given earlier in equation (2.1). Once the appropriate slip modes are assigned, this modification naturally carries through in the strain energy (equation (2.3)). The slip modes most often studied and used in larger-scale crystal plasticity codes for bcc materials are the {110} and {112} slip modes [63,64]. In a similar fashion, the results presented here will focus on the {110}〈111〉 slip mode.

Accounting for the dependence of the Peierls barrier on dislocation character in the PFDD model is a more complex problem. As mentioned above, for dislocations in a bcc crystal, it is important to identify the screw/edge character of the dislocation at all points along the loop. It has been found via measurements that the screw dislocations move much slower than the edge, with all else being the same [50–53]. Because this is an effect that is not dependent on the internal or external stress states, it must be accounted for in the core energy term (equation (2.6)). The equation for the core energy itself remains the same sinusoidal function, because perfect dislocations rather than partial dislocations are modelled. The key difference is that for bcc metals the amplitude *B* is no longer a material-dependent constant as for fcc metals. Rather *B* varies with the character of the dislocation line such that it reaches a minimum where the dislocation line is pure edge, representing a low energy barrier for slip. Conversely, *B* achieves a maximum where the dislocation has a pure screw orientation. Maximum and minimum values for these energetic barriers are material-dependent. To distinguish between the edge and screw parts, a geometrical criterion based on the angle between the line normal to the dislocation line and Burgers vector was applied. When this angle is 0°, then the dislocation is of pure edge character. When equal to 90°, then it is pure screw character. When lying in between, it is mixed. The amplitude, *B*, of the core energy varies accordingly.

To illustrate the PFDD model extensions for bcc crystal structures, we compare the expansion of a single dislocation loop in an fcc metal (Cu) and a bcc metal (Nb). Figure 12*a* shows the initial set-up for a loop in Nb. The set-up for Cu is the same; however, the orientation of the cell is as shown in figure 3. The grain size in both the Cu and Nb simulation is ∼20 nm. Figure 12*b,c* shows how the dislocation loop expands in Cu and Nb, respectively. In the case of Cu, the loop expands evenly, because the barrier for slip is isotropic. Conversely, an oval shape quickly emerges in the Nb loop. The edge-oriented dislocation segments glide through the crystal lattice much faster than the screw-oriented segments. This distinction can have a non-negligible effect on the deformation response of the crystal [50–53].

## 7. Bi-phase materials

We now present model extensions to the PFDD framework for heterogeneous or bi-phase materials. Heterogeneous materials encompass a wide range of materials, such as polycrystalline materials, materials with cracks or voids, and composite materials with bi-material interfaces. Whether the heterogeneity is a misoriented grain, a dissimilar material, a void or free surface, the differences in elastic constants create an internal stress field that causes image forces. These forces can impact the motion of dislocations.

In the inhomogeneous PFDD model, recently presented by Lei *et al.* [65], the image forces caused by the moduli mismatch are calculated by making use of Eshelby's equivalent inclusion method [7]. In this method, an inhomogeneous system consisting of a matrix containing inhomogeneities is replaced with a homogeneous system comprised of the matrix and an eigen- or virtual strain, *ϵ*^{v} that is non-zero in the regions where the inhomogeneity had existed. Notationally, we will denote quantities in the matrix material with a superscript (1) and those in the heterogeneous region with a superscript (2).

For a system that has inhomogeneities, *E*^{strain} first presented in equation (2.3) can be rewritten in terms of the energy of the equivalent, homogeneous system, *E*^{eq}, plus the difference in energy, Δ*E*, between the inhomogeneous and homogeneous systems. The complete *E*^{strain} can be expressed as [65]
7.1
where *ϵ*_{ij}(**x**) is the total strain. The first integral for *E*^{eq} is taken over the whole domain, and the second integral for Δ*E* is taken only over the region that contains the inhomogeneity. In equation (7.1), the strain is defined as
7.2
In addition, we define and as the stiffness tensors in the matrix material and inhomogeneous region, respectively. The quantity is only defined in the region with the inhomogeneous region. Using stress equilibrium and the principle of superposition, we have the following expression for the total strain *ϵ*_{ij}(**x**) [65]:
7.3
where is the average stress-free strain, *V* is the volume of the computational domain, and is the compliance tensor in the matrix material. Substituting equations (7.2) and (7.3) into equation (7.1) allows us to express *E*^{strain} in terms of the plastic and virtual strains [65]:
7.4
where . Equation (7.4) is similar to equation (2.3); however, the strain energy now accounts for differences in elastic moduli through the virtual strain. The first three terms constitute the internal strain energy. Similar to equation (2.3), the first integral describes elastic interactions including dislocation–dislocation interactions, and the second describes the interaction with the externally applied stress *σ*^{appl}_{ij}. The key difference is that equation (7.4) not only describes these interactions with the plastic component of the system, but also with the virtual strain. The final integral in equation (7.4) accounts for the elastic strain energy due to the presence of the inhomogeneity, which generates an internal stress state as a consequence of to the difference in elastic moduli. This term was not included in the earlier formulation.

At every time increment, both the phase field variables and virtual strain are determined by the minimization of the total energy in the system, *E*, which is now composed of equation (7.4) and equation (2.6) for perfect dislocations. As before, minimization is carried out with the TDGL kinetic equations; however, now there is an additional evolution equation for the virtual strain [4,5,65]. The two equations are given as
7.5
where *K* is a material constant related to the inhomogeneous region. The phase field variables are associated with plasticity and these evolve in the entire domain in order to account for plasticity throughout the system. In contrast, the virtual strain evolves only within the region of the inhomogeneity. Outside of this region, the virtual strain is zero. As before we assume quasi-static loading for the results discussed in the following sections.

### (a) Case study: voids and free surfaces

This extension of the PFDD model to account for inhomogeneities has been used by Lei *et al.* [65] to study dislocation evolution close to a void surface. Simulations consisted of a periodic array of through-thickness cylindrical holes in nickel under an increasing externally applied stress. Under this applied stress, the dislocations glide following the first of equations (7.5). The second part of equation (7.5) accounts for eigenstrains associated with the presence of the void. Voids with radii of 16*b* and 32*b* were considered in addition to a homogeneous material without voids. For the sake of simplicity, dislocation glide was restricted to the (111) slip system.

Simulation results showed that the steady state dislocation structure is not evenly distributed throughout the domain. Rather, the dislocation density distribution varies owing to the presence of an inhomogeneous distribution of shear stress caused by the void. Lei *et al.* also examined the stress–strain response of the configurations with voids and the homogeneous medium. The calculated stress–strain curves show that the domain with the largest void radius (32*b*) has the smallest elastic slope and the lowest yield stress. We note that Lei *et al.* [65] have also applied this model to investigate crazing in glassy polymers.

### (b) Bi-material interfaces

In this section, we extend the inhomogeneous PFDD model to account for a bi-material interface. In application of the PFDD inhomogeneous model, this bi-material system consists of material one, the matrix material, and material two, the inhomogeneous region. Properties of materials one and two will be indicated by superscripts (1) and (2). Here, we discuss the addition of two local interface effects to the formulation: (i) the presence of misfit strains due to differences in the lattice parameter and (ii) the energy required to form a residual dislocation in the interface, should a dislocation transmit across it.

For demonstration, we apply the model to a bicrystal comprised of two fcc crystals that differ in their elastic moduli and lattice parameters and are joined at a single common interface. Figure 13 shows the basic model geometry of the two-phase bicrystal. The interface in this model is assigned a crystallographic character, described by the crystallographic plane of each crystal corresponding to the interface plane and the orientation relationship of the two crystals. In the examples to follow, the interface will either be coherent or semi-coherent with a cube-on-cube orientation relationship and a (001) interface plane. With a cube-on-cube orientation relationship, the same axes in crystals (1) and (2) are parallel and the slip planes and slip directions are continuous across the interface. As a consequence, the same 12 phase field variables corresponding to the 12 slip systems in fcc metals are used for the entire bicrystal system. For instance, the phase field variable associated with the slip system is the same variable whether in material (1) or material (2). It is the same order parameter but evaluated at two different positions *x*. With the phase field variables known, equation (2.1) can be used to calculate *ϵ*^{p} in each crystal.

#### (i) Misfit strains

In the bicrystal model, the lattice parameters *a*^{(1)} and *a*^{(2)} of material (1) and material (2) will generally not be equal. A lattice parameter mismatch at the interface generates a localized misfit strain as the atoms must distort in the region near the interface in order to attain a coherent interface. The misfit stress corresponding to this strain is equal in magnitude but opposite in sign in material (1) and material (2) in order to maintain equilibrium. The material with the larger lattice parameter is in compression, whereas the material with the smaller lattice parameter is in tension.

In the PFDD model individual atomic positions are not modelled, and therefore, to account for the lattice mismatch we add the misfit strain to the formulation of the total system energy. It should be noted that the formulation we select to use applies to a bi-material system with a low lattice parameter mismatch (approx. <3.5%) such that the interfaces are either coherent or semi-coherent with widely spaced misfit dislocations and coherent regions in between. For larger mismatches between lattice parameter, a different formulation would be required.

First, consider a local coordinate system (*x*_{1},*x*_{2},*x*_{3}) lying in the interface plane. For an interface with a cube-on-cube orientation relationship, the assumption of plane stress in the plane of the interface can be appropriately adopted. This leads to the following basic formulation for the misfit strain [66]:
7.6

In the local coordinate system, the misfit strain components in equation (7.6) can be written as
7.7
where *C*^{(c)}_{eff}=*E*^{(c)}/(1−*ν*^{(c)}), *E*^{(c)} is the Young's modulus, and *ν*^{(c)} is the Poisson's ratio for *c*=1 and 2. The and are the in-plane normal strain components in the *x*_{1} and *x*_{2} directions, and is the out-of-plane normal strain component, normal to the interface plane in the *x*_{3} direction. The elastic moduli near the interface are set equal to the bulk elastic moduli. We note that the PFDD model does allow for variation in the elastic moduli. Hence, *C*^{(c)}_{eff} could be replaced by a better estimate of the moduli within the interface region should such a measurement or calculation become available.

Once the local misfit strain has been transformed into the global coordinate system, it can be incorporated into the total strain under an assumption of additive decomposition of strain: 7.8 We note that equation (7.8) is a modification of equation (7.3), which shows the total strain for the inhomogeneous biphase system. The strain energy is revised by inserting the modified total strain in equation (7.8) into equation (7.1), which gives 7.9 The misfit strain results in additional contributions to the energy, which are found in the second and third lines above. The second and third integral terms in the second line describe the strain energy associated with the misfit strain. The additional integral term in the third line accounts for the interaction between the external applied stress and misfit strain. We note that the misfit strain is non-zero only within a localized volume close to the interface.

#### (ii) Energy associated with forming an interfacial residual dislocation

In this section, we further extend the PFDD formulation to include the formation energy of the residual dislocation that is left behind after the dislocation crosses the bi-material interface. Formation of a residual dislocation will add an additional energy term, *E*^{res}, to the total system energy giving
7.10
where *E*^{res} is non-zero in the event that a dislocation transmits from material (1) to material (2) across the interface.

When a dislocation from the donor material (material (1)) transmits to the recipient material (material (2)) across an interface, the Burgers vectors must be conserved. Hence, a residual dislocation is left at the interface with a Burgers vector **b**^{r}=**b**^{(2)}−**b**^{(1)} after the dislocation has transmitted into material (2). In the bicrystal model, the interface region is a volume comprised of two planes, the interface plane contributed by material (1) and the interface plane contributed by material (2), and has surface area *S* and width *w*. Let **x**^{(I1)} be the subset of points **x** located along the interface plane contributed by material (1) and **x**^{(I2)} be the subset of points along the interface plane contributed by material (2). The phase field variables in the interface region are denoted as *ζ*(**x**^{(I1)}) and *ζ*(**x**^{(I2)}) (or *ζ*^{(1)} and *ζ*^{(2)} in short). For clarity, we highlight that these are not separate or different phase field variables from *ζ* introduced earlier.

When the residual dislocation has formed, the net displacement caused by the residual dislocation **u**_{r} is given by
7.11

The displacement **u**_{r} will give rise to tractions *τ*^{(1)} and *τ*^{(2)} on both the material (1) and material (2) side of the interface, respectively. The energy *E*^{res} associated with formation of the residual is given by
7.12

The distortions in each material caused by the residual dislocation are given by
7.13
and the stress generated by the plastic distortion is approximated by
7.14
where is the symmetric part of . We have made a very simple approximation that is directly related to through a nonlinear modulus for the interface . Values for are not generally known for the coherent and semi-coherent interfaces considered here. Consistent with an assumption we made earlier with the misfit strains in equation (7.7), we equate with the elastic moduli of the donor and recipient materials. Finally, the tractions are given by
7.15
where **n** is the interface normal.

### (c) Case study: slip transmission across bi-metal interfaces

In this section, we apply the model to calculate the critical stress *τ*_{crit} required to transmit a perfect dislocation across the bi-material interface in figure 13. A series of composite material systems are studied in an attempt to understand the role of lattice and moduli mismatch on *τ*_{crit}. Of these, we first consider a bicrystal made of Cu and Ni. These materials differ in lattice parameter by 2.5%, and in effective Young's moduli and shear moduli by approximately 50%. Other material parameters used in the Cu–Ni PFDD simulations are given in table 1 presented previously in §2a.

We apply a uniform shear stress *σ*_{xz} (for edge orientation), which is sufficient to drive the edge dislocation to the Cu–Ni interface. Once it reaches the interface, *σ*_{xz} is incrementally increased until it is just large enough to cause the dislocation to transmit through it. This gives the critical resolved shear stress *τ*_{crit} required for transmission.

Figure 14*a*–*c* shows snapshots of the PFDD simulation for an edge dislocation transmitting through a Cu–Ni interface. The dislocation line direction ( for an edge dislocation) points out of the plane of the page. The colour represents the value of the phase field variable, where green corresponds to a value of 1 and are areas slipped by the perfect edge dislocation. The critical stress for transmission is calculated to be 2.58 GPa. In figure 14*c*, we see that the dislocation has fully traversed the recipient material leaving a residual dislocation in the interface (red point), which is an edge dislocation and is expressed with a non-integer value of the phase field variable. We note that because the dislocation is moving from the material with the larger Burgers vector (Cu) to a material with a smaller Burgers vector (Ni), the fractional component of the phase field variable that represents the residual dislocation is greater than one. Conservation of the Burgers vector requires that a residual dislocation **b**^{r} forms with a Burgers vector magnitude given by *b*^{Cu}=*b*^{r}+*b*^{Ni}, where *b*^{r} should equal 0.0064 nm when using the values presented in table 1. The PFDD value of the residual Burgers vector *b*^{r} in figure 14*c* is 0.0064 nm, as expected.

We repeated the simulations for the reverse path, with the dislocation being driven to transmit from Ni to Cu. In this case, the recipient material has the smaller shear modulus and larger Burgers vector (or lattice parameter). Figure 14*d*–*f* shows the evolution of an edge dislocation crossing from Ni to Cu. We find that *τ*_{crit}=5.36 GPa is significantly higher than in the previous simulation of an edge dislocation transmitting from Cu to Ni. As expected from the conservation of the Burgers vector, the magnitude of the residual Burgers vector left in the interface is the same as that left by an edge dislocation: 0.0064 nm.

For this case, we see an asymmetry in *τ*_{crit} with respect to the direction of slip transmission. Transmission from Cu to Ni is easier than from Ni to Cu. Transmission from the softer to the stiffer material and from the larger to smaller Burgers vector is favoured. Repeating the calculations for transmission of a screw dislocation shows the same trends [3].

The asymmetry is not related to the value of the residual dislocation because it is the same regardless of the direction taken. It also does not arise from changes in self-energy of the dislocation. With the strain energy per unit length of a dislocation proportional to *Gb*^{2}, dislocations in Cu have the lower self-energy. However, in simulation, transmission from Ni to Cu is not favoured compared with the reverse. The asymmetry apparently does not arise solely from changes in moduli. From our calculations, the Koehler stress (or the virtual stress calculated here due to the moduli mismatch at the interface) has a weak influence on *τ*_{crit}. The Koehler stress impacts the stress required for the dislocation to approach the interface and less so the stress required for the dislocation to cross it. To explore the effects of lattice and moduli mismatch on this apparent asymmetry in *τ*_{crit}, we carried out similar simulations to those shown in figure 14 but for different fcc/fcc composites. The material parameters for all systems studied here are summarized in table 1. These two-phase systems studied have a low lattice parameter mismatch (approx. ≤3.5%), such that the bi-material interface is either coherent or semi-coherent. Among them, Al–Pt has the greatest and Ag–Au the least lattice and moduli mismatch.

As seen in the Cu–Ni simulations, the PFDD slip transmission simulations for these composite systems show an asymmetry in *τ*_{crit}. The results are tabulated in table 2. In all cases, *τ*_{crit} is higher when moving from a material with a smaller *a* to one with a larger *a*. The *τ*_{crit} is higher when transmission takes place from the stiffer to softer material, with the exception of Ag–Au. Also, the greater the lattice mismatch, the larger *τ*_{crit} is, with the exception of Cu–Ni and Pt–Al. In the light of these exceptions, it is clear that both the lattice and moduli mismatch collaborate to contribute to the observed asymmetry in *τ*_{crit}.

Based on the PFDD results and foregoing analysis, it is natural to postulate that *τ*_{crit} is governed by the formation energy of the residual dislocation. Like *τ*_{crit}, this energy barrier would be path-dependent, sensitive to whether the dislocation transmits from the softer to stiffer material and/or from the material with the smaller to larger lattice parameter.

To isolate the effects of the residual dislocation formation energy, we use the PFDD formulation to develop an analytical expression for the energetic barrier required to form the residual dislocation [3]. Starting with equation (7.12) for *E*^{res}, the peak value of *E*^{res} with respect to *ζ*^{(2)} is given by
7.16
and the critical resolved shear stress on the transmitting dislocation *τ*_{loc} associated with *E*^{res}_{barrier} is given by
7.17

The final expression for *τ*_{loc} suggests that *τ*_{crit} ought to scale with a bi-material factor, . Figure 15 compares *τ*_{crit} as calculated with the PFDD model for all systems studied here against this factor. As shown, all calculated *τ*_{crit} values follow well this scaling.

This scaling suggests a path dependence in the formation energy of the residual dislocation. This energy penalty results in an asymmetry in *τ*_{crit}, in which dislocation transmission is harder when going from a stiffer material with a smaller lattice parameter to a softer material with a larger lattice parameter.

## 8. Near future extensions

In this article, we have discussed a mesoscale technique, called PFDD, for modelling mechanically deformed crystals. The fundamental framework of the PFDD model is energy-based; hence it determines the motion and configuration of individual dislocations by minimizing the total system energy. As we demonstrated through several examples, PFDD permits modelling of discrete dislocations in finite-sized single crystals and polycrystals. To recover some physics at the atomic scale, the PFDD formulation has been extended to incorporate the energetics of dislocation cores from DFT calculations. This extension has enabled studies of partial dislocation nucleation and deformation twinning [1,2,6,43]. We have also presented modifications of the PFDD to treat two-phase metals and dislocation–interface interactions. Future extensions of the PFDD technique include treating polycrystals with texture or fcc/bcc and fcc/fcc bi-materials with partial dislocations.

## Authors' contributions

I.J.B. and A.H. equally contributed to the preparation of the manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

Los Alamos National Laboratory, an affirmative action equal opportunity employer, is operated by Los Alamos National Security, LLC, for the National Nuclear Security Administration of the US Department of Energy under contract no. DE-AC52-06NA25396.

## Acknowledgements

I.J.B. acknowledges support from the Los Alamos National Laboratory Directed Research and Development (LDRD) Programme through project no. 20140348ER. A.H. acknowledges support from the Los Alamos National LDRD Programme through project 20160156ER.

## Footnotes

One contribution of 11 to a theme issue ‘Trends and challenges in the mechanics of complex materials’.

- Accepted November 26, 2015.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.