Development of a Modified Embedded Atom Force Field for Zirconium

Jul 7, 2016 - The best set of MEAM parameters are determined by employing a ... The Journal of Physical Chemistry Letters 2016 7 (19), 3752-3759...
0 downloads 0 Views 5MB Size
Download PDF
Article pubs.acs.org/JPCC

Development of a Modified Embedded Atom Force Field for Zirconium Nitride Using Multi-Objective Evolutionary Optimization Badri Narayanan,† Kiran Sasikumar,† Zhi-Gang Mei,‡ Alper Kinaci,† Fatih G. Sen,† Michael J. Davis,§ Stephen K. Gray,† Maria K. Y. Chan,† and Subramanian K. R. S. Sankaranarayanan*,† †

Center for Nanoscale Materials, ‡Nuclear Engineering Division, and §Chemical Sciences and Engineering Division, Argonne National Laboratory, Lemont Illinois 60439 United States S Supporting Information *

ABSTRACT: Zirconium nitride (ZrN) exhibits exceptional mechanical, chemical, and electrical properties, which make it attractive for a wide range of technological applications, including wear-resistant coatings, protection from corrosion, cutting/shaping tools, and nuclear breeder reactors. Despite its broad usability, an atomic scale understanding of the superior performance of ZrN, and its response to external stimuli, for example, temperature, applied strain, and so on, is not well understood. This is mainly due to the lack of interatomic potential models that accurately describe the interactions between Zr and N atoms. To address this challenge, we develop a modified embedded atom method (MEAM) interatomic potential for the Zr−N binary system by training against formation enthalpies, lattice parameters, elastic properties, and surface energies of ZrN (and, in some cases, also Zr3N4) obtained from density functional theory (DFT) calculations. The best set of MEAM parameters are determined by employing a multiobjective global optimization scheme driven by genetic algorithms. Our newly developed MEAM potential accurately reproduces structure, thermodynamics, energetic ordering of polymorphs, as well as elastic and surface properties of Zr−N compounds, in excellent agreement with DFT calculations and experiments. As a representative application, we employed molecular dynamics simulations based on this MEAM potential to investigate the atomic scale mechanisms underlying fracture of bulk and nanopillar ZrN under applied uniaxial strains, as well as the impact of strain rate on their mechanical behavior. These simulations indicate that bulk ZrN undergoes brittle fracture irrespective of the strain rate, while ZrN nanopillars show quasi-plasticity owing to amorphization at the crack front. The MEAM potential for Zr−N developed in this work is an invaluable tool to investigate atomic-scale mechanisms underlying the response of ZrN to external stimuli (e.g, temperature, pressure etc.), as well as other interesting phenomena such as precipitation.

■

INTRODUCTION Transition metal nitrides such as zirconium nitride (ZrN) are of broad technological importance owing to their excellent physical and mechanical properties, such as high hardness, high melting point, high electrical conductivity, to name a few. Such properties have enabled their use in a range of different applications from coating materials for surface protection of cutting and surgical tools to automotive and aerospace components and other parts subject to high wear and corrosive environments. In coating applications, interface properties, such as adhesion strength and interfacial diffusivity between the metal substrate and coating, dictate the structural stability.1,2 Thus, it is vital to understand, predict, and control the interface properties of coatings. Furthermore, the structural characteristics of the coatings have a significant bearing on mechanical and physical properties of the tools underneath. For example, with proper grain structure it is possible to prevent erosion and wear in precision medical tools such as scalpels, scissors, elevators, curets, vice grips, and so on that are used for surgical purposes.3−5 While several experimental studies have been performed in order to understand the relationships between the processing parameters, the film structure, and the mechanical © XXXX American Chemical Society

properties, a fundamental understanding of the structure− property relationship at the atomistic level remains scarce. In addition to its applications as coating material, zirconium nitride is also a promising diluent material for actinide nitride targets in accelerator driven systems (ADS), a type of nuclear reactors. Nitride fuels used as advanced fuel materials possess high density, which generates excess neutrons and has a higher potential to transmute the long-lived fission products. Considering the breeding ratio, thermophysical properties (high thermal conductivity, high melting point, high fuel density), chemical compatibility with the Na coolant, and reprocessing feasibility, actinide nitrides appear to be a compromise between oxide and metal fuels. ZrN is added as a diluent, and the nitride fuels typically form single-phase solid solutions under normal irradiation conditions. Understanding the thermophysical properties of ZrN is important for developing the technologies for nuclear fuel cycle based on nitrides. Received: May 25, 2016 Revised: July 6, 2016

A

DOI: 10.1021/acs.jpcc.6b05296 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

GA18 followed by a local optimization method (simplex)19,20 to optimize the values of the independent parameters for the MEAM potential model. In all local and global optimization schemes, the goodness of fit and hence, the predictive capability of the force field, depends on the weighting scheme selected between different observables. The selection of weights is difficult because the errors often have different orders of magnitude and originate from different number of measurements. In order to obtain a more reliable and accurate fit, a more systematic way of weighting is required. One way of removing the effects of weighting is to use a multiobjective GA optimization algorithm, which can give a set of solutions at the Pareto front and enable selection of force field parameters according to desired properties.21 We therefore utilize a multiobjective GA optimization algorithm to parametrize our MEAM potential model for ZrN. Where possible, we also make comparisons with available experimental information or firstprinciples calculations. Note the experimental/theory data employed for comparison was not used in the training set. This paper is organized as follows: Methods describes potential model formalism, that is, the MEAM functional forms as well as the training data set, and the fitting procedure and parametrization strategies we employ for potential model development. Results reports our results on the multiobjective evolutionary optimization, the optimized MEAM parameters obtained in this study, the structural and energetic properties predicted by these parameters, and their success in describing interatomic interactions in ZrN. Performance of Newly Developed MEAM provides representative examples for applying our newly developed MEAM to compare the mechanical deformation of ZrN nanostructures with bulk. Finally, Conclusions summarizes the key findings and provides concluding remarks.

Experimental measurements of thermophysical and interface properties are often difficult. With the exponential growth in computing power, atomistic simulations are now emerging as a viable approach to calculate such properties. While firstprinciples calculations of material properties would be ideal, there are still limitations to the system sizes (∼1000 atoms) that are tractable with quantum calculations. Classical molecular dynamics (MD) simulations present a promising alternative to first-principles methods for the study of thermophysical properties as well as interface dynamics. MD simulations with semiempirical potentials can indeed be used to model systems with sizes up to billions of atoms. The accuracy of MD calculations, however, relies strongly on the quality of the interatomic potentials employed. Typically, these potentials are parameterized to accurately reproduce the values of various measurable material properties, for example, lattice parameters, cohesive energies, elastic constants, thermal conductivity, surface energies, and so on. Despite the clear need for largescale atomistic simulations, the amount of effort devoted to development of interatomic potentials for nitrides, especially ZrN, has been meager. To date, there is only one semiempirical potential function for ZrN that was developed by Adachi and co-workers.6 This interatomic potential has a simplified pairwise form consisting of Morse and Busing-Ida terms.6,7 The parameters for this potential were semiempirically determined by fitting to a limited number of experimental values that included mainly the variation of lattice parameter with temperature and pressure. Binary systems like Zr−N are particularly difficult to describe via interatomic potentials since the constituent elements Zr and N are very different from each other; this makes accounting for Zr−Zr, N−N, and Zr−N interactions within one mathematical formalism (i.e., potential function) extremely challenging. For a given training data set consisting of thermodynamic, structural, and elastic properties, it is important to (a) choose a proper functional form of the potential model to adequately describe the various many body interactions, and (b) utilize a parametrization strategy that can optimize multiple objectives (e.g., elastic constants and surface energies). Here, we develop the first many-body interatomic potential model for Zr−N that successfully captures thermophysical, structural, and surface properties. In terms of the potential formalism, we choose the modified embedded-atom method (MEAM) interatomic potential since it can simultaneously describe a wide range of crystal structures (body-centered cubic (bcc), face-centered cubic (fcc), hexagonal close-packed (hcp), diamond and even gaseous elements). MEAM was originally developed by Baskes and represents a modification of the EAM to include the directionality of bonding. While the original MEAM includes only the first nearest-neighbor interactions,8 a recently developed enhanced version, or the 2NN formalism,9,10 considers second nearest-neighbor interactions; in this work, we adopted the 2NN-MEAM formalism, henceforth referred to as MEAM. This MEAM potential formalism has been previously used for a wide range of systems: pure bcc, fcc, and hcp metals; carbon; gases (H, N, O),11−13 and binary compounds including nitrides and carbides (Fe−C,14 Fe−N,15 Ti−C,12 and Ti−N12). The determination of force-field parameters that have a complex functional form requires a systematic approach, which can enable effective and efficient sampling of the parameter space. Motivated by our recent success in employing global optimization methods based on genetic algorithms (GA) for fitting force fields,16,17 we use

■

METHODS Modified Embedded Atom Method. In the framework of MEAM, the total energy E of the system is expressed as8−10 E=

⎡ 1 ∑ ⎢⎢Fi(ρi ) + 2 i ⎣

⎤ ϕ S ( R ) ∑ ij ij ij ⎥⎥ ⎦ j≠i

(1)

where Fi is the embedding function for atom i, which is associated with a background electron density ρi , while Sij is the screening function and ϕij(Rij) is the pair interaction energy between atoms i and j that are a distance Rij apart. The embedding function F for a given background electron density ρ is defined as ρ ⎛ρ⎞ F(ρ ̅ ) = AEc ̅ ln⎜⎜ ̅ ⎟⎟ ρ0 ⎝ ρ0 ⎠

(2)

where A is an adjustable parameter, while Ec and ρ0 are, respectively, cohesive energy and electron density for a reference structure (i.e., a configuration in which all the atoms are on ideal lattice sites). For single elements, the most stable (equilibrium) crystalline phase under ambient conditions is a typical choice for the reference structure. In this study, the hexagonal-close packed structure is chosen as reference for Zr following ref 22; for nitrogen, however, the N2 molecule is chosen as reference structure similar to ref 15. To account for bond directionality effects, the background electron density ρi at each site i is evaluated by combining contributions from a B

DOI: 10.1021/acs.jpcc.6b05296 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

In addition, we determine all the MEAM parameters (Ec, re, α, ρ0, and screening parameters) for describing interactions between Zr and N atoms. Training Data Set. To develop accurate and transferable interatomic potentials, it is imperative to employ a training data set that amply samples the potential energy landscape. In the MEAM formalism, there are certain parameters for the Zr−N interactions, namely, Ec, re, and α that hold direct relationship with the physical properties of the reference structure, as explained in the previous section. We derived the values of these parameters directly from density functional theory calculations of material properties of rocksalt ZrN, that is, cohesive energy, equilibrium nearest-neighbor distance, and bulk modulus. For the remaining parameters, we employed a DFT-derived training set containing (a) equation of state, (b) elastic constants, and (c) energies of low index surfaces [(001), (110), and (111)] of rocksalt ZrN, as well as heat of formation and lattice parameters of another commonly observed compound, Zr3N4. All the DFT calculations are performed using the projectoraugmented wave (PAW) formalism, as implemented in the Vienna Ab Initio Simulation Package (VASP).24,25 We employ PAW atom potentials for Zr and N, without semicore p-states, supplied with the VASP package.24,25 The exchange correlation is described by the Perdew−Burke−Ernzerhof (PBE) generalized gradient approximation (GGA) functional,26 which is known to describe metal nitrides well.27 The plane wave energy cutoff is set at 500 eV, while Brillouin zone is sampled by a Γcentered Monkhorst−Pack grid. For calculation of heat of formation, equation of state, structural optimization, and elastic constants, we employ a computational supercell containing one conventional unit cell of rocksalt ZrN (containing 8 atoms) or orthorhombic Zr3N4 (containing 28 atoms). A k-point grid of 6 × 6 × 6 is employed. A conjugate gradient method is used to perform atomic relaxations until the force components on any atom are less than 0.01 eV/Å. In terms of low-index surfaces of rocksalt ZrN, we calculated surface energies for three orientations, namely, (001), (110), and (111). The surface energy for a slab S with given surface orientation is obtained by subtracting the bulk energy Eb from the total energy of the slab 1 Es, that is, S = 2 (Es − NE b), where N is the number of atoms present in the slab. All surface calculations are performed on the smallest possible surface unit cell; periodic boundary conditions are imposed in the plane of the slab, while a vacuum of 12 Å is included along the direction normal to the slab. For (001) and (110) surfaces, we used slab configurations consisting of seven layers. Since ZrN has rocksalt structure, each layer of (111) slab contains a single element, that is, it has alternating layers of Zr and N. To maintain the stoichiometry of the (111) slab, and consequently enable the use of bulk energy of rocksalt ZrN as the chemical potential, we employed a (111) slab configuration with eight layers, which results in an overall dipole. Dipole corrections are therefore applied for (111) in order to remove interactions among periodic slabs. The value of ZrN(111) surface energy calculated using such a configuration provides the average value between two possible surface terminations, i.e., Zr and N. Note that our DFT derived lattice parameters, cohesive energies, and elastic constants are within 1, 2, and 10% of experimental values, similar to expected errors in DFT-PBE framework.28 The experimental values of surface energies of ZrN are not available. However, we note that

spherically symmetric term, as well as several angular terms with weight factors t(h) i (h = 1−3). Each of these partial electron density terms are related to the atomic configuration, and atomic electron density ρa(h) from all the atoms j (≠i) that lie j within the second-nearest neighbor distance of atom i. The atomic electron density ρa(h) from an atom j at a distance Rij j from atom i is given by ρja(h) = ρ0 e−β (h)(R ij / re− 1)

(3)

where ρ0 is the scaling factor for atomic electron density, β(h) are decay lengths, and re is the equilibrium distance in the reference structure. As shown in eq 2, the embedding function has a specific functional form. In contrast, the pair interaction ϕij is obtained from the zero-temperature equation of state Eu (Rose et al.; ref 23) for a uniform expansion/contraction in the reference structure. Eu is written as * Eu(R ) = −Ec(1 + a* + da*3 )e−a

(4)

where d is an adjustable parameter, R is the nearest-neighbor distance, and a* is evaluated as

a* = α(R /re − 1)

(5)

⎛ 9BΩ ⎞1/2 α=⎜ ⎟ ⎝ Ec ⎠

(6)

Here B is the bulk modulus and Ω is the atomic volume of the reference structure. Furthermore, in this work, the modified version of MEAM that accounts for second-nearest neighbor interactions (i.e., 2NN-MEAM) is employed; in 2NN-MEAM formalism, the contributions of second-neighbor interactions is screened less severely. The details of these multi-body screening functions Sij as well as 2NN-MEAM formalism are provided elsewhere.9,10,22 To extend the MEAM formalism to binary systems, such as Zr−N, it is necessary to determine the interactions between pairs of dissimilar elements. This can be achieved via a technique similar to that described for single elements explained above, as described in recent work on developing MEAM for Ti−N and Ti−C systems.12 For the Zr−N system, we employ the rocksalt structure as the reference structure, whose total energy per atom can be written as12 u EZrN (R ) =

1⎡ ⎢FZr( ρZr ) + FN( ρN ) + Z1ϕZrN(R ) 2⎣ ⎤ Z + 2 (SZrϕZrZr(aR ) + S NϕNN(aR ))⎥ ⎦ 2

(7)

where Z1 and Z2 are number of first- and second-nearest neighbors, a is the ratio between second- and first-nearest neighbor distances in the reference structure, while SZr and SN are screening factors (whose values depend on Cmin and Cmax) for second neighbor interactions between Zr and N atoms; for rocksalt structure, the values of Z1 and Z2 are 6 and 12, respectively. The universal equation of state EuZrN(R) can also be obtained via eqs 4 and 5 using cohesive energy, atomic volume, and bulk modulus of the reference rocksalt ZrN structure obtained from first-principles calculations or experiments. This provides a way of estimating interaction between Zr and N atoms ϕZrN(R). In this work, the MEAM parameters for Zr is taken from ref 22. For N, we reparameterize A, t(1), t(2), and t(3) while retaining the rest of the parameters as reported in ref 15. C

DOI: 10.1021/acs.jpcc.6b05296 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

parameters of Zr3N4, and (b) Δ(2), which is the sum of squared errors in predictions of surface energies of rocksalt ZrN. For each individual in the initial randomly chosen population, we evaluate all the material properties in the training set using molecular dynamics package LAMMPS30 and determine the objective function values Δ(1) and Δ(2). Thereafter, we perform various genetic operations on these individuals (parents) to use their traits and generate new parameter sets (offsprings). These operations include tournament selection without replacement, crossover using the simulated binary method, or mutation via a polynomial of order 20.18,31,32 The probability of a crossover operation is set to 0.9, while that of mutation is 0.1. The mutations are necessary to introduce sufficient diversity into the population, and avoid premature convergence of the GA run. After the genetic operations, both the old (parents) and the newly created (offsprings) individuals are ranked using a nondominated sorting algorithm,18,31 and the best Np individuals are retained in the subsequent generation. This routine is repeated for sufficient number of generations; at each generation, this technique results in a Pareto front containing Np possible solutions. As the multiobjective GA optimization evolves, that is, as the generation number increases, the Pareto front progressively moves toward lower values of both objective functions (Δ(1) and Δ(2)). Eventually, beyond a certain number of generations, the Pareto front does not advance anymore leading to the convergence of the GA run (See ref 33 for detailed description of Pareto front).

typically DFT-PBE surface energies are within 10% of experimental values.16,28 Evolutionary Strategy To Optimize MEAM Parameters. Using the training set described above, we employ global optimization via genetic algorithms (GA) with multiple objectives to optimize the 15 independent MEAM parameters, namely, (a) A, t(1), t(2), and t(3) for N−N interactions, and (b) scaling factor for atomic electron density ρN0 /ρZr 0 , and screening NZrN ZrZrN ZrNN NNZr ZrN−Zr parameters: CZrNZr , CNZrN min , Cmin , Cmin , Cmin , Cmin , Cmax max , ZrZrN ZrNN NNZr Cmax , Cmax , and Cmax for Zr−N interactions. We employ the single/multiple objective genetic algorithm toolbox developed by Sastry and co-workers to perform all the genetic operations;29 the overall procedure is outlined in Figure 1.

■

RESULTS Determination of MEAM Parameters. Using multiobjective evolutionary optimization technique detailed in Figure 1, we determine the MEAM parameters for the Zr−N binary system by fitting against structural, thermodynamic, elastic, and surface properties of rocksalt ZrN as well as structure/energetics of Zr3N4. Figure 2 shows the evolution of

Figure 1. Schematic representation of multiobjective evolutionary optimization strategy employed in this work to parametrize MEAM interatomic potential for ZrN.

The optimization scheme starts with a population of Np = 100 different MEAM parameter sets, whose values are set randomly within physically meaningful ranges. Each of these parameter sets that constitute the population is called an individual. The fitness (suitability) of a particular individual is assessed by the accuracy of its predicted properties as compared to the target value (training data set). Typically, a weighted sum of errors in predicted properties termed objective function is used to describe fitness of an individualthe higher the objective function value, the lower its fitness/suitability. This technique, however, faces an inherent challenge−the evolution of GA optimization is strongly influenced by the weight factors, whose choice is arbitrary and is largely driven by intuition. For instance, it is often difficult to achieve an optimal parameter set that describes both bulk as well as surface properties accurately using a single objective function. This is primarily related to difficulty in assigning relative weights to the relative contributions of errors in predicted bulk/surface properties. To mitigate this challenge, we employed two objective functions (instead of one), namely, (a) Δ(1), which is the sum of squared errors in predictions of bulk properties, that is, elastic constants of ZrN, heat of formation and lattice

Figure 2. Evolution of the Pareto front with generation number g during a typical multiobjective evolutionary optimization run for two objective functions, namely, the sum of squared errors in the prediction of structural and elastic properties Δ(1), and that for surface energies Δ(2). The values of Δ(2) are scaled by a factor of 100 to enhance clarity. The final optimized parameter set obtained in this study is represented by the orange star.

the Pareto front as a function of generations during our multiobjective GA run; at each generation, the Pareto front provides a collection of non-dominated MEAM parameter sets in the objectives Δ(1) and Δ(2). Initially (g = 0), the population of MEAM parameter sets is randomly chosen within a physically allowable domain in the parameter space. As D

DOI: 10.1021/acs.jpcc.6b05296 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C Table 1. MEAM Parameters for Pure Zr and Na Zr N

Ec

re

α

A

β0

β1

β2

β3

t(1)

t(2)

t(3)

Cmin

Cmax

d

6.36 4.88

3.2 1.1

4.45 5.96

0.68 1.71

2.45 2.75

1.00 4.0

3.00 4.0

2.00 4.0

6.30 0.49

−3.30 1.23

−10.00 −1.97

1.0 2.0

1.44 2.8

0.00 0.00

a The cohesive energy Ec and equilibrium nearest-neighbor distance re are provided in eV/atom and Å, respectively. All other parameters are dimensionless. The parameters for Zr are taken from earlier work.22 For N, the values for A, t(1), t(2), and t(3) are reparameterized in this study, while the remaining parameters are taken from an earlier report.12

expected, this results in parameter sets that exhibit extremely high objective function values (poor suitability) as illustrated in Figure 2. As the optimization run proceeds, fitter MEAM parameter sets (with lower objective values) are sampled, which causes the Pareto front to advance gradually toward lower values of Δ(1) and Δ(2) (Figure 2). Eventually, the GA run converges and the Pareto front stops advancing at ∼400 generations, giving a population of Np best candidate parameter sets. Among these parameter sets in the converged Pareto front (g = 400), we identify one set near the knee33 of the front, which best reproduces the training set, as well as ensures that (a) the ZrN structure is stable at finite temperatures up to its melting point and that (b) no other polymorph of ZrN is more stable than the reference state, that is, rocksalt structure. The optimized MEAM parameter set is provided in Tables 1 and 2.

Table 3. Lattice Parameters and Heat of Formation of Zirconium Nitride Compounds Using 2NN-MEAM Developed in This Worka MEAM

value

Ec re α d CZrNZr min CNZrN min CZrZrN min CZrNN min CNNZr min CZrNZr max CNZrN max CZrZrN max CZrNN max CNNZr max ρN0 /ρZr 0

7.387 2.305 4.475 0.0 0.46 0.71 1.57 1.46 0.24 1.25 2.80 6.25 1.58 1.97 11.565

experiments

4.61 −1.77

4.61,b 4.585c 1.79b

9.91 10.90 3.30 −1.58

9.79d 10.85d 3.30d

ZrN a (Å) ΔH (eV/atom)

4.61 −1.77

a (Å) b (Å) c (Å) ΔH (eV/atom)

10.01 11.10 3.37 −1.51

Zr3N4

a

The corresponding values from our DFT calculations, and previous experiments (when available) are also provided for comparison. bRef 34. cRef 35. dRefs 36 and 37.

Table 2. MEAM Parameters for the Binary Zr−N System Obtained in This Work parametera

DFT

ZrN. Consequently, our MEAM parameters could exactly reproduce the DFT values for heat of formation ΔH and lattice parameter for rocksalt ZrN. The MEAM predicted lattice parameter for ZrN is within ∼0.5% of previous experimental reports; similarly, the value of ΔH for rocksalt ZrN predicted by MEAM is in excellent agreement with experiments (within ∼1%). The newly obtained MEAM parameters provide an accurate description of structure and energetics of Zr3N4 as well. As shown in Table 3, the MEAM predicted lattice parameters for Zr3N4 are within ∼2% of both DFT values as well as experiments. Similarly, the ΔH for Zr3N4 obtained from MEAM is in near-perfect agreement with DFT calculations (Table 3). In addition, the equation of state (i.e., energyvolume dependence) for ZrN and Zr3N4 predicted by MEAM is in good accordance with those obtained from DFT calculations (Figure 3). Table 4 compares the elastic constants Cij and surface energies for low index surfaces of ZrN (Figure 4) predicted by MEAM with those obtained from DFT calculations and experiments (when available). The elastic stiffness values for ZrN predicted by our newly developed MEAM potential are in excellent agreement with DFT calculations and experiments (within ∼10%). These MEAM parameters accurately reproduce the dependence of surface energy on orientation as compared to DFT calculations; in the framework of both MEAM and DFT, (100) surface is the most stable followed by (110) and (111) surfaces. Furthermore, the absolute values of the surface energies obtained from MEAM show an error of ∼15% with respect to their DFT values (Table 4). Such a close match is particularly difficult to achieve with interatomic potentials, including MEAM. This clearly demonstrates the power and applicability of multi objective evolutionary optimization to parametrize accurate interatomic potentials. Fracture Behavior of Bulk ZrN and Nanopillar. To demonstrate the ability of our newly developed MEAM potential to capture atomic scale dynamical phenomena in Zr−N systems, we investigate the response of bulk and

a

The cohesive energy Ec and equilibrium nearest-neighbor distance re are provided in eV/atom and Å, respectively. All other parameters are dimensionless.

Performance of Newly Developed MEAM. We assess the accuracy of our newly developed MEAM parameters (Tables 1 and 2) by calculating the physical properties of Zr−N compounds and comparing them with our DFT calculations as well as previous experiments.34−39 The 2NN-MEAM formalism employed in this work includes interactions up to secondnearest neighbor interactions. Therefore, we employ a radial cutoff distance of 5.26 Å, which lies between the second- and third-nearest neighbor distances of Zr.22 Table 3 compares the MEAM predicted lattice parameters and heat of formation for ZrN and Zr3N4 with the values obtained from our DFT calculations and experiments. As aforementioned, the potential parameters Ec and re for Zr−N interactions are obtained by fitting exactly to the DFT-derived values of cohesive energy and nearest-neighbor Zr−N distance for reference rocksalt phase of E

DOI: 10.1021/acs.jpcc.6b05296 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

nanopillar ZrN to applied strain at different strain rates using molecular dynamics (MD) simulations. Orthorhombic computational supercells are used to simulate both (a) bulk (4.6 nm × 4.6 nm × 23 nm) and (b) nanopillar (6.8 nm × 6.8 nm × 13.7 nm) configurations. For bulk ZrN, periodic boundary conditions are employed in all directions, while for nanopillar, only the axial direction (z) is considered to be periodic. The surfaces normal to the nonperiodic directions (i.e., x and y) in the nanopillars are the crystallographic (100) planes, which is the lowest energy surface for rocksalt ZrN (Table 4). Furthermore, the system sizes employed here (for both bulk and nanopillars) possess a cross-sectional area of at least 10 × 10 unit cells, which is sufficiently large (∼8× the radial cutoff distance of MEAM) to avoid any spurious effects from periodic images. All the MD simulations in this work are performed in LAMMPS using a time step of 0.5 fs. First, both the bulk and nanopillar ZrN are equilibrated under ambient conditions (300 K, 1 atm) using a Nosé−Hoover thermostat and barostat.40 These thermalized configurations are then subjected to uniaxial tensile loading along the axial (z) direction at uniform strain rates. Three different strain rates, 108, 109, and 1010 s−1, are considered for both bulk and nanopillar systems. Note that while these rates are 2−3 orders of magnitude faster than that in typical dynamic/shock loading experiments, computational limitations restrict slowest strain rates accessible via MD simulations to 108 s−1 (e.g., ref 41). Figure 5 shows the stress−strain curves for both the bulk and rectangular nanopillar systems at different strain rates. For all

Figure 3. Equation of state for (a) rocksalt ZrN and (b) orthorhombic Zr3N4 calculated using MEAM parameters developed in this study (blue filled circles) and DFT calculations (red open circles). The heat of formation (ΔH) for each crystal (provided in eV per unit formula, i.e., eV/u.f.) is relative to its enthalpy of formation at equilibrium (ΔH0) as evaluated by the corresponding level of theory. Also, the crystal volumes are normalized by the value at equilibrium in the framework of the corresponding level of theory.

Table 4. Comparison of Elastic Constants and Energies of Low-Index Surfaces of ZrN Predicted by the 2NN-MEAM Potential Developed in This Work with DFT Calculations and Experimentsa MEAM C11 C12 C44 (100) (110) (111)

DFT

Elastic Constants 480 512 83 114 125 113 Surface Energies 1.26 1.48 2.29 2.78 3.39 2.95

experiments 471,b 454c 88,b 121c 138,b 124c

a

The elastic constants are reported in GPa, while the surface energies are provided in J/m2. bRef 38. cRef 39. Figure 5. Stress−strain curves for (a) bulk and (b) rectangular nanopillar ZrN for three different strain rates: 108, 109, and 1010 s−1.

strain rates, both bulk and nanopillar ZrN show elastic behavior until 5% strain and a yield stress of 15 GPa. Bulk ZrN, and to a certain extent the nanopillar ZrN, then undergo quasiplastic stress relaxations. The yield region seen beyond the elastic limit but before brittle fracture has been previously reported for ceramics, particularly at high strain rates.41−44 Bulk ZrN undergoes brittle failure at the fracture strain with a sudden drop in the stress. The fracture strain is also seen to increase with increasing strain rate, similar to previously reported studies on ceramics.41,42 The stress−strain behavior of a nanopillar is fundamentally different from the bulk. The stress relaxation for a strain rate of 108 s−1 is seen to be a continuous decrease from the yield stress, while the two higher strain rates show a sharper

Figure 4. Atomic structures (top view) of ZrN surfaces with different orientations, (a) 001, (b) 110, and (c) 111. The Zr atoms are depicted using green spheres, while N are shown as blue spheres. For clarity, the surface unit cell (whose edges are shown in white) is repeated twice along the principal directions in the plane of the surface.

F

DOI: 10.1021/acs.jpcc.6b05296 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

For all strain rates, as reported previously for brittle ceramics, the crack is initiated by the formation of voids.41 The crack propagates rapidly and within 5−10 ps the sample separates along a near-flat (100) crack plane, as expected from a rocksalt structure brittle ceramic. As shown in Figure 5, the stress−strain behavior of a nanopillar is fundamentally different from the bulk wherein the stress relaxes gradually beyond an intrinsic strength of the material. To understand the atomistic origin of the continuous stress relaxation, we depict atomic snapshots at selected times from a representative MD run on ZrN nanopillars (i.e., at 109 s−1 strain rate) in Figure 7. The gradual stress relaxation is seen

drop, followed by a regime of continuous stress relaxation. Such gradual stress relaxations have been previously reported for metallic nanowires.45 Table 5 summarizes the results of the tensile loading simulations on bulk and nanopillar ZrN for all the strain rates. Table 5. Strain Rate Dependence of Young’s Modulus Y, Fracture Stress σf, Fracture Strain ϵf, and Intrinsic Strength for Bulk and Nanopillar ZrN Obtained from Our MD Simulations bulk strain rate (s−1) 108 109 1010

nanopillar

Y (GPa)

σf (GPa)

ϵf

Y (GPa)

intrinsic strength (GPa)

406 406 407

21.7 22.8 22.4

0.099 0.110 0.124

371 397 400

16.0 18.0 21.5

We see that for the bulk ZrN, the Young’s modulus is independent of the strain rate and is equal to the elastic constant predicted by the interatomic potential. The fracture stress does not significantly depend on strain rate, while the fracture strain increases with increasing strain rate.41−44 For the rectangular nanopillar, both the Young’s modulus and intrinsic strength, the stress before the gradual stress relaxation stage, increase with increasing strain rate. Figure 6 depicts snapshots from the MD simulation of tensile loading of bulk ZrN at a strain rate of 108 s−1. The stress

Figure 7. Atomistic representation of brittle fracture of nanopillar ZrN. The corresponding strains (with reference to the unstrained length of ZrN) are shown for each atomistic snapshot. The molecular dynamics simulation for the above representative case was performed at a strain rate of 109 s−1. The atoms are colored by the per-atom potential energy. Unlike the bulk ZrN, significant amorphization is seen near the crack plane. After brittle fracture (at ∼10%), the crack propagation proceeds gradually with the formation of a narrow connected neck and eventual separation of the cracked ZrN pillars at ∼50% strain.

Figure 6. Atomistic representation of brittle fracture of bulk ZrN. The corresponding strains (with reference to the unstrained length of ZrN) are shown for each atomistic snapshot. The molecular dynamics simulation for the above representative case was performed at a strain rate of 108 s−1. The atoms are colored by the per-atom potential energy. Between 5 and 7% strain, stress relaxation is seen via the formation of small angle tilt boundaries. Subsequent strain hardening occurs before fracture at ∼10% strain.

to be associated with slow crack propagation where the two nanopillars retain contact after fracture via the formation of a gradually thinning neck. Such local continuous amorphization in the nanopillar has similarities with previously reported studies for metallic nanowires.45 They reported the possibility of continuously transforming a homogeneous perfect crystal to an amorphous metal by the application of sufficiently high strain rates. The stress−strain behavior for a ceramic nanowire shows similar tendencies with a continuous stress relaxation. However, in contrast to a metallic nanowire, the amorphization is restricted locally to the fracture plane. By and large, the fracture is still of brittle nature. In Figure 8, we show the spatiotemporal variation of the density of ZrN for both the bulk and nanowire cases strained with a strain rate of 108 s−1 along the z direction. The plots are colored based on local density in 1D bins along the z direction

relaxation beyond the elastic limit (5% strain) is seen to be due to the formation of small angle tilt boundaries. This is accompanied by dislocation-like propagation of local lowcoordination defects in the bulk (see Supporting Information, Figure S1) and a propagation of two stress fronts that meet at the eventual crack plane (see Supporting Information, Figure S2). No necking is observed indicative of the brittle nature of the ceramic. However, local amorphization is possible near the crack plane before fracture, particularly at the high strain rate of 1010 s−1 (see Supporting Information, Figures S3 and S4). This is not surprising since at higher strain rates the sample has less time to relieve its strain by structural rearrangement of atoms. G

DOI: 10.1021/acs.jpcc.6b05296 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

These simulations indicated that bulk ZrN undergoes brittle fracture, regardless of the strain rate; the strain at fracture increases as the strain-rate increases, while the fracture stress, and the Young’s modulus remain largely unchanged. On the other hand, ZrN nanopillars exhibit quasi-plasticity, which is facilitated by localized amorphization similar to previous studies on metallic nanowires. These MD simulations demonstrate the capability of the MEAM parameters developed in this work to capture dynamical atomic-scale events in Zr−N systems that occur in response to external stimuli, for example, temperature and applied pressure. Finally, the Zr−N potential developed here can be easily extended to include other metals (whose MEAM parameters are already available) to investigate interesting atomic-scale phenomena in alloy nitrides, such as precipitation and mechanical behavior.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b05296. Atomic snapshots at selected times during fracture of bulk ZrN and nanopillars at different strain rates (PDF).

Figure 8. Density fluctuation along strained direction (y) vs strain (x) for (top) bulk ZrN and (bottom) nanopillar ZrN at 300 K. In both cases, the strain rate is 108 s−1 along the z-direction.

■

of the simulation cell. The vertical axis of the color maps represents the spatial response along the strained direction, while the horizontal axis is the strain. Since we apply a constant strain rate in our simulations, the horizontal axis is also a proxy for time. For the bulk case, the density of ZrN decreases homogeneously during elastic deformation. A sharp transition can be seen at fracture, indicative of brittle failure. The density in the cracked ZrN sample on either side of the crack plane is inhomogeneous and fluctuates discontinuously. This is indicative of shock propagation in ZrN soon after brittle fracture in the ceramic. In contrast, the nanopillar case shows a smooth and homogeneous transition up to and after fracture. The underlying mechanisms for fracture for the bulk or the nanopillar case are fundamentally similar for the three simulated strain rates (see Supporting Information, Figure S5). The success of our newly developed MEAM in capturing atomistic mechanisms that occur in ZrN in response to external environment (e.g., applied strains) opens up avenues for employing this tool to gain atomistic insights into various interesting material phenomena, such as shock loading, atomic transport near grain boundaries, and so on. Furthermore, the multiobjective evolutionary framework for parametrizing interatomic potentials presented in this work can be extended to any material, and establishes a systematic protocol to develop force fields for multicomponent systems.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Use of the Center for Nanoscale Materials was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-AC0206CH11357. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC0205CH11231. We gratefully acknowledge the computing resources provided on Blues and Fusion, high-performance computing clusters operated by the Laboratory Computing Resource Center at Argonne National Laboratory.

■

REFERENCES

(1) Holleck, H. Material Selection For Hard Coatings. J. Vac. Sci. Technol., A 1986, 4, 2661−2669. (2) Sproul, W. D. New Routes in the Preparation of Mechanically Hard Films. Science 1996, 273, 889−892. (3) Krishnan, V.; Krishnan, A.; Remya, R.; Ravikumar, K. K.; Nair, S. A.; Shibli, S. M. A.; Varma, H. K.; Sukumaran, K.; Kumar, K. J. Development and Evaluation of Two PVD-Coated Beta-Titanium Orthodontic Archwires For Fluoride-Induced Corrosion Protection. Acta Biomater. 2011, 7, 1913−1927. (4) Shieu, F. S.; Cheng, L. H.; Sung, Y. C.; Huang, J. H.; Yu, G. P. Microstructure and Coating Properties of Ion-Plated TiN on Type 304 Stainless Steel. Thin Solid Films 1998, 334, 125−132. (5) Jeongwon, P.; Duk-Jae, K.; Yoon-Kee, K.; Keun-Ho, L.; Ki-Hoon, L.; Saeyoung, A. Improvement of the Biocompatibility and Mechanical Properties of Surgical Tools with TiN Coating by PACVD. Thin Solid Films 2003, 435, 102−7. (6) Adachi, J.; Kurosaki, K.; Uno, M.; Yamanaka, S. A Molecular Dynamics Study of Zirconium Nitride. J. Alloys Compd. 2005, 396, 260−263. (7) Ida, Y. Interionic Repulsive Force and Compressibility of Ions. Phys. Earth Planet. Inter. 1976, 13, 97−104.

■

CONCLUSIONS In summary, we have developed a 2NN-MEAM interatomic potential for Zr−N binary system by training against (a) lattice parameters and cohesive energies of ZrN and Zr 3 N 4 compounds, as well as (b) elastic constants and (c) surface energies of ZrN. To optimize MEAM parameters, we employed a global optimization scheme empowered by multiobjective genetic algorithms. Our newly obtained set of MEAM parameters accurately describe the structure, energetics, dynamics, and elastic properties of bulk ZrN, as well as its surface properties. Using the Zr−N MEAM developed in this work, we investigated the fracture behavior of ZrN in bulk and nanopillar configurations at atomic scale using MD simulations. H

DOI: 10.1021/acs.jpcc.6b05296 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C (8) Baskes, M. I. Modified Embedded Atom Potentials for Cubic Materials and Impurities. Phys. Rev. B: Condens. Matter Mater. Phys. 1992, 46, 2727−2742. (9) Byeong-Joo, L.; Baskes, M. I.; Kim, H.; Yang Koo, C. Second Nearest-Neighbor Modified Embedded Atom Method Potentials for BCC Transition Metals. Phys. Rev. B: Condens. Matter Mater. Phys. 2001, 64, 184102/1−11. (10) Lee, B. J.; Shim, J. H.; Baskes, M. I. Semiempirical Atomic Potentials for the FCC Metals Cu, Ag, Au, Ni, Pd, Pt, Al, and Pb Based on First and Second Nearest-Neighbor Modified Embedded Atom Method. Phys. Rev. B: Condens. Matter Mater. Phys. 2003, 68, 144112. (11) Zhou, X. W.; Wadley, H. N. G.; Filhol, J.-S.; Neurock, M. N. Modified charge transfer embedded atom method potential for metalmetal oxide systems. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 065402. (12) Young-Min, K.; Byeong-Joo, L. Modified embedded-atom method interatomic potentials for the Ti-C and Ti-N binary systems. Acta Mater. 2008, 56, 3481−9. (13) Lee, B.-J.; Jang, J.-W. A Modified Embedded-Atom Method Interatomic Potential for the Fe-H System. Acta Mater. 2007, 55, 6779−6788. (14) Byeong-Joo, L. A. Modified Embedded Atom Method Interatomic Potential for the Fe-C System. Acta Mater. 2006, 54, 701−11. (15) Byeong-Joo, L.; Tae-Ho, L.; Sung-Joon, K. A Modified Embedded Atom Method Interatomic Potential for the Fe-N System: A Comparative Study with the Fe-C System. Acta Mater. 2006, 54, 4597−607. (16) Narayanan, B.; Kinaci, A.; Sen, F. G.; Davis, M. J.; Gray, S. K.; Chan, M. K. Y.; Sankaranarayanan, S. K. Describing the Diverse Geometries of Gold From Nanoclusters to Bulk - a First-Principles Based Hybrid Bond Order Potential. J. Phys. Chem. C 2016, 120, 13787−13800. (17) Sen, F. G.; Kinaci, A.; Narayanan, B.; Gray, S. K.; Davis, M. J.; Sankaranarayanan, S. K. R. S.; Chan, M. K. Y. Towards Accurate Prediction of Catalytic Activity in IrO2 Nanoclusters via First Principles-Based Variable Charge Force Field. J. Mater. Chem. A 2015, 3, 18970−18982. (18) Deb, K. Multi-Objective Optimization Using Evolutionary Algorithms; John Wiley and Sons: Chichester, U.K., 2001. (19) Nelder, J. A.; Mead, R. A Simplex Method for Function Minimization. Computer Journal 1965, 7, 308−313. (20) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes: The Art of Scientific Computing, 3rd ed.; Cambridge University Press: New York, U.S.A., 2007. (21) Jaramillo-Botero, A.; Naserifar, S.; Goddard, W. A. General Multiobjective Force Field Optimization Framework, with Application to Reactive Force Fields for Silicon Carbide. J. Chem. Theory Comput. 2014, 10, 1426−1439. PMID: 26580361. (22) Young-Min, K.; Byeong-Joo, L.; Baskes, M. I. Modified embedded-atom method interatomic potentials for Ti and Zr. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 74, 14101−1−12. (23) Rose, J. H.; Smith, J. R.; Guinea, F.; Ferrante, J. Universal Features of the Equation of State of Metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1984, 29, 2963−2969. (24) Kresse, G.; Furthmuller, J. Efficiency of Ab-initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. Comput. Mater. Sci. 1996, 6, 15−50. (25) Kresse, G.; Furthmuller, J. Efficient Iterative Schemes for Abinitio Total Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169−11186. (26) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (27) Marlo, M.; Milman, V. Density-Functional Study of Bulk and Surface Properties of Titanium Nitride Using Different ExchangeCorrelation Functionals. Phys. Rev. B: Condens. Matter Mater. Phys. 2000, 62, 2899−2907. (28) Lejaeghere, K.; van Speybroek, V.; van Oost, G.; Cottenier, S. Error Estimates for Solid-State Density-Functional Theory Predic-

tions: An Overview by Means of the Ground-State Elemental Crystals. Crit. Rev. Solid State Mater. Sci. 2014, 39, 1−24. (29) Sastry, K.; Goldberg, D. E.; Kendall, G. In Search Methodologies: Introductory Tutorials in Optimization and Decision Support Techniques; Burke, E. K., Kendall, G., Eds.; Springer: Berlin, 2005; Chapter 4, pp 97−125. (30) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (31) Deb, K.; Pratap, S. M. T.; Agarwal, A. A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II. Trans. Evol. Comp 2002, 6, 182−197. (32) Deb, K.; Agrawal, R. B. Simulated Binary Crossover for Continuous Search Space. Complex Systems 1995, 9, 115−148. (33) Bechikh, S.; Ben Said, L.; Ghédira, K. Searching for Knee Regions of the Pareto Front using Mobile Reference Points. Soft Computing 2011, 15, 1807−1823. (34) B, P. W., Ed. In Structure Reports; International Union of Crystallography: Oosthoek, Scheltema, and Holkema, Utrecht, 1913− 1993. (35) Stampfl, C.; Mannstadt, W.; Asahi, R.; Freeman, A. J. Electronic Structure and Physical Properties of Early Transition Metal Mononitrides: Density-Functional Theory LDA, GGA, and Screened-Exchange LDA FLAPW Calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2001, 63, 155106. (36) Lerch, M.; Füglein, E.; Wrba, J. Synthesis, Crystal Structure, and High Temperature Behavior of Zr3N4. Z. Anorg. Allg. Chem. 1996, 622, 367−372. (37) Baur, W. H.; Lerch, M. On Deciding Between Space Groups Pnam and Pna21 for the Crystal Structure of Zr3N4. Z. Anorg. Allg. Chem. 1996, 622, 1729−1730. (38) Chen, X. J.; Struzhkin, V. V.; Wu, Z. G.; Somayazulu, M.; Qian, J.; Kung, S.; Christensen, A. N.; Zhao, Y. S.; Cohen, R. E.; Mao, H. K.; et al. Hard Superconducting Nitrides. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 3198−3201. (39) Christensen, A. N.; Dietrich, O. W.; Kress, W.; Teuchert, W. D. Phonon Anomalies In Transition-Metal Nitrides: ZrN. Phys. Rev. B: Condens. Matter Mater. Phys. 1979, 19, 5699−5703. (40) Evans, D. J.; Holian, B. L. The Nose-Hoover Thermostat. J. Chem. Phys. 1985, 83, 4069−4074. (41) Pedone, A.; Malavasi, G.; Menziani, M. C.; Segre, U.; Cormack, A. N. Molecular Dynamics Studies of Stress-Strain Behavior of Silica Glass under a Tensile Load. Chem. Mater. 2008, 20, 4356−4366. (42) Swiler, T. P.; Simmons, J. H.; Wright, A. C. Molecular Dynamics Study of Brittle Fracture in Silica Glass and Cristobalite. J. Non-Cryst. Solids 1995, 182, 68−77. (43) Fischer-Cripps, A.; Lawn, B. Indentation Stress-Strain Curves for Quasi-Ductile Ceramics. Acta Mater. 1996, 44, 519−527. (44) Kim, B.-N.; Hiraga, K.; Morita, K.; Sakka, Y. A High-Strain-Rate Superplastic Ceramic. Nature 2001, 413, 288−291. (45) Ikeda, H.; Qi, Y.; Ç agin, T.; Samwer, K.; Johnson, W. L.; Goddard, W. A. Strain Rate Induced Amorphization in Metallic Nanowires. Phys. Rev. Lett. 1999, 82, 2900−2903.

I

DOI: 10.1021/acs.jpcc.6b05296 J. Phys. Chem. C XXXX, XXX, XXX−XXX