Impact of twin boundaries on bulk elastic constants: Density-functional theory data for Young׳s modulus of Ag

Experimental and theoretical studies on nanowires have reported a size-dependence of the Young׳s modulus in the axial direction, which has been attributed to the increasing influence of surface stresses with decreasing wire diameter. Internal interfaces and their associated interface stresses could lead to similar changes in the elastic properties. In Kobler et al. [1], however, we reported results from atomistic calculations which showed for Ag that twin boundaries have a negligible effect on the Young׳s modulus. Here, we present data of density-functional theory calculations of elastic constants and Young׳s modulus for defect-free bulk Ag as well as for bulk Ag containing dense arrays of twin boundaries. It is shown that rigorous convergence tests are required in order to be able to deduce changes in the elastic properties due to bulk defects in a reliable way.


a b s t r a c t
Experimental and theoretical studies on nanowires have reported a size-dependence of the Young's modulus in the axial direction, which has been attributed to the increasing influence of surface stresses with decreasing wire diameter. Internal interfaces and their associated interface stresses could lead to similar changes in the elastic properties. In Kobler et al. [1], however, we reported results from atomistic calculations which showed for Ag that twin boundaries have a negligible effect on the Young's modulus.
Here, we present data of density-functional theory calculations of elastic constants and Young's modulus for defect-free bulk Ag as well as for bulk Ag containing dense arrays of twin boundaries. It is shown that rigorous convergence tests are required in order to be able to deduce changes in the elastic properties due to bulk defects in a reliable way.  Subject area  Materials Science  More specific  Elastic properties of defect-free  subject area  crystals and with twin boundaries  Type of data  Tables and graphs  How data  The density of the k-point mesh and the width σ of the broadening scheme, which is used to determine the occupation numbers, are the most crucial parameters.

Computational methods
Density-functional theory (DFT) calculations were carried out with the plane-wave code PWscf of the Quantum Espresso software package [2], using the Perdew-Burke-Ernzerhof PBE exchangecorrelation functional [3], Vanderbilt ultrasoft pseudopotentials [4] and a plane wave kinetic energy cutoff of 30 Ry. k-point meshes for Brillouin zone integrations were generated by the Monkhorst-Pack scheme [5], and the fractional occupation numbers of the electronic states were determined by a Gaussian broadening [6]. Atomic positions and lattice parameters were relaxed by minimizing the atomic forces and the stress tensor, with a convergence threshold for the largest residual force and stress component of 5 meV/Å and 0.1 kbar, respectively.
Bulk second order elastic constants (SOEC) were calculated by finite deformations of the conventional face-centered cubic (fcc) unit cell with 4 atoms, Cartesian basis vectors and lattice constant a. The lattice is distorted by applying a small strain ε, which transforms the basis vectors a 1 , a 2 , a 3 to the new vectors e a 1 ; e a 2 ; e a 3 À where I is the 3 Â 3 identity matrix. The three independent elastic constants of the cubic lattice C 11 , C 12 and C 44 were determined by using the deformation strains proposed by Mehl et al. [7]: The first deformation is a homogeneous volume change, which changes the total energy of the unit cell by V 0 is the equilibrium volume of the unit cell and B is the bulk modulus. For a cubic lattice B is related to the elastic constants by The second deformation is a volume conserving orthorhombic strain with energy change and the third deformation is a volume conserving monoclinic shear with energy change From the bulk modulus B and the difference ΔC ¼ C 11 ÀC 12 the two elastic constants C 11 and C 12 are given by For all three deformations a series of total energy calculations was performed for a small set of finite strain values x. The energy values as function of x were fitted to a polynomial and the elastic constants were determined from the second derivative of the polynomial at the minimum. Specifically, the bulk modulus B was calculated by varying the lattice constant a between 7.64 Bohr and 8.04 Bohr in steps of 0.04 Bohr, which corresponds to strain values between À 2.4% and þ2.7%. A polynomial of degree 4 was fitted to the 11 data points. The minimum of the curve gives the equilibrium lattice constant, which was used as starting point for the calculation of C 11 À C 12 and C 44 . The deformations proposed by Mehl have the advantage that the energy is an even function of the strain x. Thus, only positive values for x have to be considered. For the calculation of C 11 À C 12 the strain parameter x was changed between 0 and þ0.084 in steps of þ0.012 (8 values) and for the calculation of C 44 we used 7 values of x between 0 and þ0.12 in steps of þ 0.02. An even polynomial of degree 6 was fitted to the calculated total energy values. The Young's modulus E ½hkl was calculated by a similar quastistatic approach as the SOEC. The unit cell is chosen in such a way that one axis is parallel to the direction ½hkl of the applied strain and the two other axes are perpendicular to the first one. Then a set of finite tensile and compressive strains x is applied and for each strain x the cell vectors perpendicular to the strain direction and the atomic positions are relaxed (Poisson contraction). From the energy change This cell, however, does not allow to introduce twin boundaries. Therefore, the calculation of E ½110 was repeated for a second, 6-atom orthorhombic unit cell with cell vectors along [101], [121] and [111]: The cell contains three atomic layers with ABC stacking along [111], with 2 atoms per plane in the unit cell. This cell can also be used for calculating Young's modulus E ½112 . In a cubic crystal, however, E ½110 and E ½112 are equal due to symmetry. It is important to note that for the orthorhombic unit cell the Poisson contraction gives rise to a tilt of the basis vectors a 2 and a 3 , if a strain is applied in the direction of a 1 . Thus, the angle between a 2 and a 3 has to be relaxed in order to get the correct value for the Young's modulus E ½110 .
Unit cells for the bulk crystal with twin boundaries are derived from the orthorhombic cell, but with a different number of atomic planes in the [111] direction (a 3 axis). Periodic boundary conditions require that always 2 twin boundaries are included in one unit cell. We used cells with 6, 8 and 10 atomic layers (thus containing 12, 16 and 20 atoms) with stacking sequences of ABC BAC, ABCA CBAC and ABCAB ACBAC, respectively, in which the twin boundaries are separated by 3a= ffiffiffi 3 p ¼ 7:18 Å, 9.57 Å, 11.96 Å. In all cases, the Young's modulus E ½110 was calculated by changing the length of a 1 between À0.024 and þ0.030 in steps of 0.06 in units of a= ffiffiffi 2 p . This corresponds to applied strains between À3.4% and þ4.2%. A polynomial of degree 6 was fitted to the 10 total energy values.

Results of the DFT benchmark calculations
Since elastic constants are second derivatives of the total energy, DFT calculations have to be very well converged in order to be able to extract elastic constants in a reliable way (within the limits of the accuracy of the chosen functional). First, we thoroughly tested the influence of the plane wave basis set and density cut-off energies. For our choice of these cut-off energies (30 Ry and 120 Ry, respectively), the elastic constants are well converged within 70.1 GPa. Much more crucial is the convergence with respect to k-point density and Gaussian smearing parameter σ. An additional problem arises from the fact that different unit cells have to be used for the calculation of the bulk value of the Young's modulus and of structures containing twin boundaries, since k-point meshes will not be equivalent. Therefore, to be able to identify changes in elastic properties due to twin boundaries, absolute values of elastic constants have to be converged as good as possible. To have an estimate on how well our calculations are converged and what accuracy for changes in the elastic constants can be expected, we applied the following 3 step strategy: (1) First, the bulk second-order elastic constants (SOEC) are calculated and we examine their convergence with respect to the k-point mesh and the Gaussian smearing parameter σ.
(2) Then we perform direct calculations of the bulk Young's modulus (using both, the tetragonal and orthorhombic unit cell) by quasistatic tensile tests. The results of both calculations are compared to the analytical value of the Young's modulus as given by the SOFCs. In cubic crystals, the Young's moduli E ½110 and E ½112 are identical and are related to the SOEC by The variation between the three approaches indicates how well our calculations are converged. The deviation will be used as an estimate for the error margin that we have to take into account for our results of the Young's modulus. The results for step (1) are shown in Fig. 1 and Table 1. k-point grids were always of the Monkhorst-Pack type with divisions of (n, n, n) of the three basis vectors of the reciprocal lattice. Fig. 1 shows the typical behavior for the convergence of properties with k-point density n: the larger the Gaussian smearing parameter σ, the faster the k-point convergence, but the k-point-converged result depends on the value of σ. This is most obvious for the elastic constant C 12 . However, for our choices of σ of 0.005, 0.010 and 0.015 Ry the deviation from the σ¼0 limit is less than 0.1 GPa for all elastic constants. Table 1 summarizes the results of the elastic constants for the different smearing values σ together with the kpoint density, which is required for obtaining convergence within an accuracy of 70.1 GPa. Obviously, for elastic properties a much denser k-point mesh than for the calculation of energy differences (for example, surface and interface energies) are required.
In Table 1 also the experimental values for the SOEC at 0 K and 300 K are given. The DFT calculations underestimate the elastic constant by up to 20%. This is typical for the PBE functional and other functionals based on the generalized-gradient approximation (GGA). In the local density approximation (LDA), on the other hand, elastic constants are overestimated by up to 20% [11]. However, this uncertainty of DFT calculations (the dependence of the results for SOEC on the choice of functional) is not crucial for our aim: we are not interested in the absolute values of elastic properties, but only in the change of the Young's modulus after the introduction of twin boundaries.
The results of the direct, quasistatic calculations of the bulk value of the Young's modulus E ½110 in step (2) are summarized in Table 2. k-points for the tetragonal and orthorhombic unit cells were also chosen according to the Monkhorst-Pack scheme. In Table 2 the number of divisions n, however, refers to an approximately equivalent (n, n, n) k-point mesh for the conventional 4-atom cubic unit cell with lattice constant a. The lattice vectors of the tetragonal and orthorhombic unit cells have lengths of a= ffiffiffi 2 p , For a Gaussian smearing of σ¼0.015 Ry the SOEC were well converged for a k-point density of 28 and 24 (see Fig. 1). For this σ-value and the k-point mesh with n¼28, however, the Young's modulus E ½110 calculated by quasistatic tensile test still deviates by 1.8 GPa from the analytic value calculated from the SOEC (see Table 2). For a better agreement between direct calculation and analytic result from SOEC for E ½110 , σ has to be lowered to 0.010 Ry and the k-point density has to be increased to n¼32. With the settings of σ¼ 0.015 Ry and the k-point mesh of n¼28 the error margin of the DFT calculations is more than 1 GPa, for σ¼0.010 Ry and the k-point mesh with n¼32 the numerical uncertainty is reduced to about 0.3 GPa.
In Table 3 the calculated Young's moduli for the orthorhombic unit cells with twin boundaries are compared to their corresponding values for the defect-free bulk. For the calculations with σ¼0.015 Ry and the k-point mesh with n ¼28 the results for E ½110 for cells with and without twin boundaries are within our estimated error margin. For these settings the unit cells with twin boundaries were also strained in [112]-direction, which is the actual orientation of the nanowires in the experiments of Ref. [1], but no significant difference between E ½110 and E ½112 is observed.
For the more accurate settings with σ¼0.010 Ry and a k-point mesh with n¼32, the Young's modulus E ½110 is systematically smaller for cells with twin boundaries than the bulk value, and E ½110 decreases systematically with increasing twin boundary density. However, even for the highest density of twin boundaries, in which twin boundaries are separated by only 3 atomic planes, the change in E ½110 is less than 2% compared to the bulk. This can be taken as an upper limit for the modification of Young's modulus due to twin boundaries in Ag. Table 3 Results for the direct quasistatic calculation of the Young's moduli E ½110 and E ½112 for the orthorhombic unit cell with two twin boundaries.