Ab initio based empirical potential applied to tungsten at high pressure

TL;DR

This paper develops an ab initio derived empirical potential for tungsten that accurately predicts mechanical and structural properties under high pressure, validated against experiments and used to study phase transitions.

Contribution

It introduces a spline-based empirical potential combining Stillinger-Weber and embedded-atom models derived from DFT data for tungsten.

Findings

Accurately predicts elastic constants and phonons up to 100 GPa

Replicates dislocation core structures and deformation twinning

Models high-pressure bcc to fcc phase transition

Abstract

Density-functional theory forces, stresses and energies comprise a database from which the optimal parameters of a spline-based empirical potential combining Stillinger-Weber and modified embedded-atom forms are determined. Accuracy of the potential is demonstrated by predictions of ideal shear, stacking fault, vacancy migration, elastic constants and phonons all between 0 and 100 GPa. Consistency with existing models and experiments is demonstrated by application to screw dislocation core structure and deformation twinning in a tungsten nanorod. Lastly, the potential is used to study the high-pressure bcc to fcc phase transition.

Figures

Tables

Table 1. Table 1: Zero-pressure elastic constants of bcc tungsten in GPa.

	B	C11	C12	C44
MEAMa	319	550	204	147
GGAa	304	513	199	142
F-Sb	309	520	204	161
LDAc	320	552	204	149
F-Sd	310	525	203	159
F-Se	310	522	204	161
EAMf	308	520	202	159
BOPg	310	522	204	161
Expt.h	308–314	501–521	199-207	151–160

Table 2. Table 2: Table of vacancy and self-interstitial formation energies (in eV) for bcc tungsten. Entries with angle-brackets indicate that the defect in question relaxes to the dumbell configuration shown.

Defect	MEAMa	GGAa	F-Sb	GGAc	F-Sd	F-Se	EAMf	GGAg
Vac. Formation	2.99	3.17	3.58	3.11	3.56	3.63	3.57	3.56
Vac. Migration	1.73†	1.70†	1.43	1.66	2.07	1.44	2.98†	1.78
Vac. Activation	4.72†	4.87†	5.01	4.77	5.63	5.07	6.55†	5.34
$⟨$001$⟩$ Dumbell	11.15	$⟨$111$⟩$	11.53	11.74	11.51	9.82	12.20	11.49
$⟨$011$⟩$ Dumbell	9.98	10.64	9.86	10.10	9.84	9.64	9.704	9.84
$⟨$111$⟩$ Dumbell	9.73	10.31	9.58	9.82	9.55	9.82	10.56	9.55
Octahedral	11.76	12.42	11.72	11.99	11.71	10.02	12.03	11.68
Tetrahedral	10.54	$⟨$111$⟩$	10.93	11.64	11.00	10.00	$⟨$011$⟩$	11.05

Table 3. Table 3: Energy and structural relaxation of low-index free surfaces in bcc tungsten. Surface energies E 𝐸 E are given in meV/Å 2 and the changes Δ 12 subscript Δ 12 \Delta_{12} in inter-planar spacing between the first two planes of the surface, where available, are shown in parenthesis in units of percent.

	MEAMa	GGAa	F-Sb	GGAc	F-Se	GGAf	AMEAMf	BOPg
$E_{001}$$\left({\Delta }_{12}^{001}\right)$	233(-5.7)	245(-11.5)	186(-0.9)	289	183(-0.7)	487	373	237(-2.5)
$E_{011}$$\left({\Delta }_{12}^{011}\right)$	198(-3.8)	200(-3.8)	159(-1.1)	249	161(-0.5)	398	353	163(-1.0)
$E_{111}$$\left({\Delta }_{12}^{111}\right)$	204(-18.9)	216(-21.6)	—	278	—	449	314	—

Equations

$$\begin{array}{rl}& V_{tot}=\sum _{i>j}\varphi (r_{ij})+\sum _{i}U(n_i)\\ & +\sum _{\begin{array}{c}i,j>k\\ j\ne i\end{array}}p(r_{ij})p(r_{ik})q(\cos {\theta }_{jik}),\end{array}$$

$$\begin{array}{r}{n}_i=\sum _{j\ne i}\rho (r_{ij})+\sum _{\begin{array}{c}j>k\\ j\ne i\end{array}}f(r_{ij})f(r_{ik})g(\cos {\theta }_{jik}),\end{array}$$

Full text

Interatomic potentials, tungsten, dislocation core, deformation twinning

I Introduction

Tungsten is an exceptional transition metal exhibiting the highest tensile strength, melting point, and elastic modulus of any pure metal and has important applications in aerospace, energy and armament industries. Much interest has been focused on $\alpha$-W (bcc) and $\beta$-W (A15) nanostructures including nanorodsKarabacak et al. (2003, 2005); Wang et al. (2015), nanoparticlesOzawa et al. (2001); Lei et al. (2007); Shibuta and Suzuki (2008); Harilal et al. (2013), and thin filmsChopra et al. (1967); Haghiri-Gosnet et al. (1989); Weerasekera et al. (1994). Due to the technological importance of tungsten, classical interatomic potentials of various forms have been developed to study this metalAckland and Thetford (1987); Baskes (1992); Zhou et al. (2004); Derlet et al. (2007); Mrovec et al. (2007); Wang et al. (2004a); Murphy et al. (2015).

Classical potentials – popular for their favorable scaling compared to first-principles methods ($\mathcal{O}(N)$ vs. $\sim\mathcal{O}(N^{3})$) – were traditionally developed by choosing analytic functional forms with a handful of free parameters determined by fitting directly to experimental bulk data such as cohesive energy, lattice and elastic constants. The force-matching method of Ercolessi and AdamsErcolessi and Adams (1994) has facilitated the development of interatomic potentials based on ab initio calculations of relaxed crystallographic defects, metastable structures and other non-equilibrium configurations. However, many such potentials still possess a small number of free parameters and analytic functional forms which limit their transferability, requiring researchers to take great caution in choosing a potential suitable for their region of interest. It is thus desirable to produce a single potential which can be employed to study the range of materials physics, here in tungsten.

To meet this challenge we develop an unique semi-empirical potential based on a robust database of ab initio calculations that samples much of the potential-energy landscape. Our model combines the Stillinger-Weber (SW)Stillinger and Weber (1985) form with the modified embedded atom methodDaw and Baskes (1984) (MEAM) form with functions parameterized by quintic splines. Section II describes the functional form of the model, the density-functional theory (DFT) calculations comprising the large fitting database, and the genetic algorithm optimization scheme. Accuracy of the fitted potential is demonstrated in Section III by comparing MEAM to DFT for the various structural and elastic properties to which it was fit. Given that the potential is fit directly to important crystallographic defects, structural properties and elastic constants, transferability is demonstrated in Section IV by examining MEAM predictions for \frac{1}{2}$$\langle111$\rangle$ screw dislocation core structure, deformation twinning and detwinning of a nanorod, and dynamics of bcc and fcc tungsten at high pressure. Conclusions are given in Section V.

II Optimization of the empirical potential with first-principles calculations

We present a spline-based empirical potential fit to a large database of highly-converged density functional theory (DFT) calculations using a global optimization scheme based on an evolutionary algorithm.

II.1 Empirical extension of the MEAM potential

The embedded-atom (EAM)Daw and Baskes (1984); Daw et al. (1993) and MEAMBaskes (1987); Baskes et al. (1989); Baskes (1992) methods have been applied to many systems including semiconductorsDaw and Baskes (1984); Baskes (1987); Baskes et al. (1989); Daw et al. (1993); Lenosky et al. (2000); Lee (2007) and transition metalsBaskes (1992); Mishin et al. (1999); Hennig et al. (2008); Fellinger et al. (2010); Park et al. (2012). The original MEAM formalism involves a parametrized analytical functional form which accounts for bond-bending through angular functions with explicit $s$-, $p$-, $d$- and $f$-orbital characteristics. Lenosky *et al.*Lenosky et al. (2000) first parametrized the MEAM formalism through the use of cubic splines for the study of defects in Si. The use of splines for parameterizing empirical potentials increases model flexibility and efficiency, and has been successfully applied to the study of martensitic transformations in pure titaniumHennig et al. (2008), shock-loading in niobiumFellinger et al. (2010); Zhang et al. (2011) and dislocation dynamics in molybdenumPark et al. (2012). SW potentials, initially developed for the modeling of cubic-diamond Si, have been successfully applied to amorphous SiVink et al. (2001) as well as GeZ. Jian et al. (1990) and other systems.

Model flexibility is paramount when constructing empirical potentials. The addition of distinct terms to the model can improve flexibility, as demonstrated by the success of MEAM over EAM Cherne and Baskes (2001); Vella et al. (2015). In the present work, we propose an empirical extension of MEAM which includes a SW-type three-body term in the energy. Our model employs functions parameterized by quintic splines which improve performance for properties requiring continuous third derivatives of the interatomic potential, e.g. $C_{ij}$ vs. $P$ relations, over cubic splines. The total potential energy is that of SW with the addition of an embedding term $U$,

[TABLE]

where the “electronic density” $n_{i}$ at atom $i$ involves two- and three-body contributions

[TABLE]

and $\theta_{jik}$ is the angle of the triplet centered on atom $i$. This is the simplest extension of MEAM which does not include four-body terms and cannot be gauge-transformed back to the original model. It contains as special cases the SW ($U=0$), MEAM ($q=0$), and EAM ($g=q=0$) forms. This flexibility gives us the ability to fit to a large ab-initio database described in the next section.

II.2 Density-functional theory database

DFT calculations performed with vaspKresse and Hafner (1993, 1994); Kresse and Furthmüller (1996a, b) using a projector-augmented planewave basisBlöchl (1994) and Perdew-Burke-Erzenhof (PAW-PBE)Perdew et al. (1996, 1997) generalized-gradient exchange correlation approximation comprise a database of forces, stresses and energies for fitting via the force-matching method. In addition to 6$s$ and 5$d$, the 5$p$ electrons are treated as valence to improve accuracy at high pressures, where the overlap of these semi-core states is not necessarily negligible. A planewave energy cutoff of 600 eV and first-order Methfessel-Paxton smearing width of 0.1 eV are used for all calculations; k-points are sampled on a $\Gamma$-centered 40$\times$40$\times$40 mesh in the bcc brillouin zone. These quantities are chosen to ensure convergence of the total energy to within 0.1 meV/atom. Additional computational details are presented in FellingerFellinger (2013).

The ab-initio fitting database contains 596 configurations with a total of 16,860 atoms and thus 54,752 force components, stress components and energies to be fit. Configurations in the database include volumetric strains of bcc, fcc, hcp, $\beta$-W (A15), $\beta$-Ta ($\beta$-U) and $\omega$-Ti. Tetragonal strains are included for hcp and $\omega$-Ti structures to ensure accurate $c/a$ values. The database also contains elastic constants of the bcc phase at pressures between 0 and 100 GPa, in increments of 25 GPa, using volume-conserving orthorhombic and monoclinic strains of 0.5 % for $C^{\prime}=\frac{1}{2}(C_{11}-C_{12})$ and $C_{44}$, respectively. At zero pressure, configurations with orthorhombic strains up to 10 % and monoclinic strains up to 40 % are added. Unrelaxed symmetry-inequivalent configurations of $\langle$110$\rangle$ and $\langle$112$\rangle$ $\gamma$-surfaces, ideal shear strain and vacancy migration are included at five equally spaced pressures between 0 and 100 GPa. Relaxed zero-pressure structures containing a vacancy at the lattice site and halfway along the $\langle$111$\rangle$ migration path are also added. A 7$\times$7$\times$7 bcc supercell with a single atom displaced by 0.006 $\mathrm{\AA}$ is included to promote accurate force-constants and phonons via the small-displacement methodKresse et al. (1995); Alfé et al. (2001). Using a supercell of this size reduces the interaction of the displaced atom with its images across periodic boundaries and thus improves the accuracy of calculated force-constants and phonon dispersions. Relaxed low-index free surfaces as well as crowdion, octahedral, $\langle$111$\rangle$-split and $\langle$110$\rangle$-split self-interstitial configurations are included. Ab-initio MD snapshots of 125-atom bcc supercells at 1620 K, 2960 K and a liquid tungsten at 6730 K are added to improve performance for simulations at high temperature. A 36-atom hcp supercell at 100 K is also included. A Lastly, a mesh of 36 points on the BainBain and Dunkirk (1924) (bcc-fcc) and BurgersBurgers (1934) (bcc-hcp) energy surfaces at pressures of 0 GPa and 700 GPa in addition to 600 GPa for the Bain path and 800 GPa for the Burgers path are included to ensure that the potential can be used to explore properties of these close-packed phases at high pressure.

II.3 Genetic algorithm optimization

Development of the optimized potential is an iterative process of fitting, testing and database refinement. Ten to twenty fits are performed simultaneously and the resultant potentials are tested for accuracy on a range of properties. The fitting database is refined based on the results of these tests: new structures are added to correct spurious behavior and structures are removed when over-fitting is suspected. Re-tuning of algorithm inputs and/or error weights often accompanies this refinement.

A global optimization scheme combining a genetic algorithm (GA) with a local downhill optimizer provides a method for determining the spline parameters of the potential. At each iteration of the GA all potentials in the population of ten are locally optimized with 60 steps of a PowellPowell (1964) conjugate direction algorithm, then the population is sorted and bred according to total weighted least-squares error. For the presented potential forces, stresses and energies are equally weighted. Breeding is done by a stochastic combination of spline knots from two parent potentials. The following constraints are enforced by introducing a “punishment” error when not satisfied: (i) $\left|\mathrm{max}[f(r)]\right|=\left|\mathrm{max}[p(r)]\right|=1$ and (ii) $n_{i}$ lies within the domain of $U(n)$ for all $i$. If the latter constraint is violated, the embedding function is evaluated at its nearest endpoint. At each GA step, for every potential in the population, there is a 10% chance for the embedding function domain and total density (Eq. 2) to be rescaled by a transformation where $U(n)\rightarrow U(n/\alpha)$ and $n\rightarrow\alpha n$, where $\alpha$ is determined by the minimum and maximum densities at the current step. When this occurs, additional gauge symmetries in the three-body terms of Eqs. 1 and 2 are exploited so that the maximal knot values of $\left|f\right|$ and $\left|p\right|$ are equal to 1.

While forces and energies are invariant with respect to these transformations, the total error is not because constraints (i) and (ii) are always satisfied after rescaling. Because during fitting the embedding function is not extrapolated but rather evaluated at the nearest endpoint when densities lie outside the domain, energies for such configurations are not invariant under the aforementioned rescaling. Furthermore, spline functions are in general not invariant under such a transformation of their argument. Thus, performing this rescaling serves as a genetic mutation of the potential. The algorithm is exited when between successive steps the change in total error for every potential in the population is less than $10^{-3}$. Parameters for the fitted potential and plots of the seven functions are presented in the Supplementary Information, along with more detailed descriptions of the algorithms used. Henceforth all references to MEAM will pertain to the present empirically extended potential.

III Accuracy of the fitted potential

We demonstrate the accuracy of the fitted MEAM potential through the energetics of non-equilibrium structures, crystallographic defects, thermodynamic properties and phonon dispersion. All MEAM calculations in this work (other than those necessary for fitting) are performed in the Large-scale Atomic/Molecular Massively Parallel Simulator (lammps)Plimpton (1995). If at any step during an MD run the density seen by an atom exceeds the embedding function domain, the embedding energy is linearly extrapolated from the nearest endpoint.

III.1 Energetics and elastic properties

Figure 1(a) shows the GGA-DFT and MEAM energy-volume relations for six distinct phases including A15 $\beta$-W and high-energy close-packed structures. MEAM accurately predicts energies of all six phases relative to the ground state. An equilibrium bcc lattice constant of 3.189 $\mathrm{\AA}$ is predicted by both GGA-DFT and MEAM, compared to the published experimental values between 3.15 and 3.165 $\mathrm{\AA}$Davey (1925); Hartmann et al. (1931); Taylor and Doyle (1967). It is well known now that GGA tends to overestimate the lattice constant of metals; the reason for this is discussed in Wang *et al.*Wang et al. (2004b) as well as Favot and Dal CorsoFavot and Corso (1999) and references therein.

Figure 1(b) compares bcc tungsten pressure-volume relations as computed with MEAM, GGA-DFT, and measured through shock experimentsKinslow and Cable (1970). MEAM and GGA-DFT curves, obtained by static calculations of volumetric strain, are indistinguishable for pressures through 800 GPa and in excellent agreement with experimental results up to 300 GPa, indicating applicability of the fitted MEAM potential to high-pressure physics in tungsten.

Figure 1(c) compares linear thermal expansion predictions by MEAM to experimental resultsTouloukian et al. (1975) for temperatures between 300 K and the experimental melting point of 3695 K. Constant N-P-T MD simulations of 2000 atoms at P = 1 atm yield the thermal-expansion curve. Each MD simulation runs for 50,000 steps with a 1 fs timestep and the lattice constant for each temperature is determined by averaging over the last 5,000 simulation steps. MEAM shows excellent agreement with experiment up to 1,000 K and remains within 1 % of the experimental fit for all temperatures considered, indicating that the potential interpolates between temperatures included in the fitting database.

Table 1 shows the zero-pressure bcc elastic constants of the present MEAM and GGA-DFT results, compared to previous ab initio calculations and other interatomic potentials. The bulk modulus B and C11 predicted by MEAM are higher than experimental and GGA results but consistent with the LDA work of Einarsdotter et al.Einarsdotter et al. (1997). The pressure-dependence of bcc elastic constants is shown in Figure 2; MEAM does not predict a monotonic increase of Cij but remains within 21 % of the GGA-DFT values.

Figure 3 shows phonon dispersion of equilibrium bcc tungsten as computed with MEAM and DFT, compared to inelastic neutron scattering results of Chen and BrockhouseChen and Brockhouse (1964). Dispersions are calculated using the finite-displacement method in a 7$\times$7$\times$7 primitive bcc supercell. DFT dispersion agrees well with experiment but exhibits oscillations in the longitudinal (low-lying) branch near the $H$-point, a feature also found in density-functional perturbation theory results within LDA-DFTEinarsdotter et al. (1997). Overall MEAM tracks both DFT and experiment but underestimates the frequency along the L[$\xi\xi\xi$] branch, particularly near the $\omega$ mode at $\xi=\frac{2}{3}$.

III.2 Point and planar defects

Table 2 lists the energies of vacancies and self-interstitial atoms (SIAs) in bcc tungsten, essential quantities for the accurate modeling of plasticity. Present MEAM and DFT calculations use a 5$\times$5$\times$5 cubic supercell. Atomic positions are relaxed to 0.01 eV. Geometric details of bcc SIA calculations can be found in Xu and MoriartyXu and Moriarty (1996). GGA-DFT calculations of Becquart *et al.*Becquart and Domain (2007) and the present work indicate the $\langle$111$\rangle$-dumbell to be the most energetically-favorable self-interstitial, as do the present MEAM potential and F-S potentials of DerletDerlet et al. (2007) and AcklandAckland and Thetford (1987). ExperimentsDiCarlo et al. (1969); Okuda and Mizubayashi (1975) and previous MD studiesZhou et al. (2004) found the $\langle$011$\rangle$-dumbell to be the favored self-interstitial structure in tungsten, but recent workAmino et al. (2016) combining the object kinetic monte carlo (OKMC) method with dislocation loop measurements found OKMC simulations of $\langle$111$\rangle$ interstitials and 1D migration best match experiment. Vacancy formation and migration energies predicted by MEAM compare favorably with present ab-initio results and those of Becquart and DomainBecquart and Domain (2007) while existing F-S and EAM tungsten potentials are in closer agreement with GGA results of Nguyen-Manh *et al.*Nguyen-Manh et al. (2006).

Figure 4 presents unrelaxed vacancy migration pathways at five equally-spaced pressures between 0 and 100 GPa. Calculations are performed using a 127-atom 4$\times$4$\times$4 cubic bcc supercell with migration in the $\langle$111$\rangle$ direction. Overall MEAM tracks well with the DFT results; minor discrepancies are found when the vacancy lies near the lattice site and halfway between two lattice sites.

Figure 5 shows unrelaxed generalized stacking fault energies (GSFEs) at five pressures on the $\{112\}$ and $\{110\}$ planes as a function of relative displacement along $\langle$111$\rangle$ for MEAM and DFT. While bcc metals are less prone to stacking-fault formation than their close-packed counterparts, they have been observed in Fe, Nb, W and Mo-35%Re to exist on $\{112\}$ and $\{110\}$ planes, formed by the dissociation of \frac{1}{2}$$\langle111$\rangle$ dislocationsHirschhorn (1963). Relaxed GSFE curves, computed with MEAM at zero pressure, do not predict the presence of any metastable stacking fault configurations. At all pressures, MEAM agrees with DFT to within a few meV/Å2, and thus should be suitable for studying the effect pressure on \{112\}$$\langle111$\rangle$ and \{110\}$$\langle111$\rangle$ slip systems.

Table 3 shows energies and interplanar relaxations of low-index free surfaces. Present calculations employ 48-atom supercells, replicated along the surface normal with an equally-sized vacuum region and periodic boundary conditions. All results presented here predict the $\langle$011$\rangle$ surface to have the lowest energy, followed by $\langle$111$\rangle$. Finnis-Sinclair potentialsAckland and Thetford (1987); Wang et al. (2004a) tend to underestimate the surface energy with respect to GGA-DFT. The present MEAM potential compares favorably with present ab-initio results and those of Vitos *et al.*Vitos et al. (1998) but underestimates the inter-planar relaxation of the high-energy $\langle$100$\rangle$ surface by 50%. The origin of the discrepancy between GGA-DFT results of Moitra et al. and the others is unclear.

Figure 6 presents the unrelaxed ideal shear stresses and energy barriers for pressures up to 100 GPa. Ideal shear defines the upper limit of stress required to deform a perfect crystal and is fundamental to our current understanding of the strength of materials. Calculations are performed following the methodology of Paxton *et al.*Paxton et al. (1991), which uses a bcc primitive cell. MEAM accurately reproduces the GGA-DFT results for all pressures, with small discrepancies in shear stress around the extrema.

IV Transferability of the fitted potential

Transferability of the fitted potential is demonstrated by application to screw dislocation core structure, deformation twinning in a bicrystal nanorod, and the high-pressure bcc-to-fcc phase transformation.

IV.1 Dislocation core and deformation twinning

Core structure of the \frac{1}{2}$$\langle111$\rangle$ screw dislocation is determined using a cell with lattice directions $\left[1\bar{2}1\right]$, $\left[\bar{1}01\right]$, $\left[111\right]$ and periodic boundary conditions along the dislocation line. The first two lattice vectors are repeated to form a large cell containing 92,277 atoms which are displaced according to the appropriate elastic strain field. The core structure is then determined by relaxing a central region containing 54,396 atoms while the remaining atoms are fixed, ensuring that the correct boundary conditions are satisfied by the long-range anisotropic solution. This methodology is further explained in this group’s previous work on NbFellinger et al. (2010) and MoPark et al. (2012).

Figure 7 shows a non-degenerate symmetric core structure predicted by MEAM, presented as a differential displacement mapVitek et al. (1970), is in agreement with results from an existing bond-order potentialMrovec et al. (2007) and DFT-GGARomaner et al. (2010) calculation for tungsten. Existing F-S potentials predict an asymmetric coreTian and Woo (2004); Gröger et al. (2008). Our potential is also consistent with the criterion of Duesbery and VitekDuesbery and Vitek (1998), which is based on F-S calculations and states that the \frac{1}{2}$$\langle111$\rangle$ screw dislocation in bcc metals will have a symmetric core structure if $\gamma_{\{110\}}(b/3)>2\gamma_{\{110\}}(b/6)$, where $\gamma_{\{110\}}$ is the relaxed $\{110\}$ $\gamma$-surface and $b=a\sqrt{3}/2$ is the burgers vector magnitude. Relaxed values taken from Figure 5 for MEAM are $\gamma_{\{110\}}(b/3)$ = 100 meV/Å2 and $\gamma_{\{110\}}(b/6)$ = 39 meV/Å2.

While dislocation slip is fundamental to plastic deformation of bulk transition metals, twinning has been found to dominate deformation in nanocrystalline Mo, Ta and FeZhu et al. (2012). A recent studyWang et al. (2015) observed deformation twinning and detwinning during uniaxial loading and unloading of a bicrystal nanorod. The Finnis-Sinclair potential of Ackland and ThetfordAckland and Thetford (1987) was used to model this twinning and detwinning in good agreement with experiment. We simulate this deformation as a challenge for our fitted MEAM potential and to demonstrate transferability to non-equilibrium conditions and consistency with existing models.

Figure 8 displays cross-sections of a bicrystal tungsten nanorod under uniaxial stress at 300 K. The nanorod is 128 Å in diameter and 510 Å in length, with periodic boundary conditions parallel to the rod axis. A compressive strain of 10 % is applied from the top of the rod over 20 ps while atomic positions are updated using the Velocity VerletSwope et al. (1982) integrator and canonical ensemble with 1 fs timestep. The strain is then unloaded over an additional 20 ps. Multiple \{112\}$$\langle111$\rangle$ deformation twins can be seen in Figure 8(a) through (c) to nucleate at the grain boundary and grow with increasing stress. At full loading, strain is accommodated primarily by a single large deformation twin extending from the grain boundary to the rod surface. During unloading the accumulated strain is released by detwinning as can be seen in panels (d) through (f). This deformation behavior is nearly identical to the results of Wang *et al.*Wang et al. (2015), indicating the transferability of the fitted MEAM potential to modelling tungsten nanostructures and consistency with the successes of previously published potentialsAckland and Thetford (1987); Wang et al. (2015). Given that the F-S potential of Ackland and Thetford predicts an asymmetric core structure but accurately describes nanorod deformationWang et al. (2015), our current MEAM potential is well suited to study the interplay of deformation twinning and dislocation-induced plasticity in tungsten.

IV.2 High-pressure phase transformation

Finally, we investigate the transformation of tungsten to fcc. Theoretical studies have predicted that bcc tungsten becomes unstable with respect to close-packed fcc and hcp phases at extreme pressuresEinarsdotter et al. (1997) and under the conditions of strong electronic excitation during laser irradiationGiret et al. (2014), for which a $T_{e}$-dependent interatomic potential was developed to study the transitionMurphy et al. (2015). To the authors’ knowledge, fcc tungsten has only been observed in thin films formed by sputter deposition between 200 and 400 ${}^{\circ}C$ on glass, mica and rock-salt substratesChopra et al. (1967). The predicted zero-pressure lattice constants of fcc tungsten for MEAM and DFT are 4.049 Å and 4.044 Å, respectively, while Chopra *et al.*Chopra et al. (1967) found an fcc lattice constant of 4.13 Å in the aforementioned tungsten films. Here we consider the stabilization of fcc over bcc at high pressures, some accessible via diamond-anvil experiments.

Figure 9(a) compares MEAM phonon dispersions for bcc W at pressures of 30 to 1200 GPa with LDA-DFT results of Einarsdotter *et al.*Einarsdotter et al. (1997). MEAM force constants are computed using the small-displacement method, implemented in the Atomic Simulation EnvironmentBahn and Jacobsen (2002), on a 10$\times$10$\times$10 supercell with $\delta=a(P)/100$ where $a(P)$ is the cubic lattice constant at pressure $P$. LDA-DFT results employed the density-functional linear response method, norm-conserving pseudopotentials, and 5$s$5$p$5$d$6$s$6$p$ valence. As seen in Figure 3, MEAM predicts the L-$\frac{2}{3}[111]$ ($\omega$) phonon to have lower frequency compared with DFT and experiment. This mode softens with increasing pressure, albeit at a lower rate than predicted by LDA calculations. Otherwise MEAM accurately captures the other important features of bcc dispersion up to 1200 GPa. Low-pressure results (30-60 GPa) also compare favorably with the AMEAM results of Zhang and ChenZhang and Chen (2013).

Figure 9(b) compares the fcc phonon dispersion predicted by MEAM and LDA-DFT at pressures from 0 to 1200 GPa. At low pressure, where fcc is a highly unfavorable structure, MEAM does not compare well to ab initio results but correctly predicts unstable soft modes in the T[$\xi\xi 0$] and T[$\xi\xi\xi$] branches. However, the stabilization of these modes with increasing pressure is non-monatonic and particularly anomalous on the T[$\xi\xi\xi$] branch at intermediate pressures. By 1200 GPa, MEAM predicts fcc tungsten to be dynamically stable and shows excellent agreement with the LDA-DFT dispersion.

Figure 10(a) shows the elasic moduli C44 and C*′* between 400 and 500 GPa, where all Cij are positive definite. It can be seen that $\mathrm{C}^{\prime}=\frac{1}{2}\left(\mathrm{C}_{11}-\mathrm{C}_{12}\right)$ is negative for pressures below 455 GPa, reflecting the slope of the T${}_{\left[1\bar{1}0\right]}\left[\xi\xi 0\right]$ branch arbitrarily close to the $\Gamma$-point. Figure 10(b) shows this mode for pressures around 540 GPa, where long-wavelength modes are stable but the $\xi$ = 0.40 mode remains unstable. According to MEAM, this mode is the last unstable phonon in any of the considered high-symmetry lines in the Brillouin zone and stabilizes at 543 GPa. However Figure 10(c), which displays the enthalpy difference $\Delta H=H_{fcc}-H_{bcc}$ versus pressure, shows that the bcc phase remains energetically favorable until about 1134 GPa.

To summarize, present MEAM results are consistent with LDA-DFT predictions of Einarsdotter *et al.*Einarsdotter et al. (1997) in that for fcc $C_{44}$ is stable at relatively low pressures, $C^{\prime}$ stabilizes before fcc is thermodynamically favorable, the last phonon mode to become real is the $T_{[1\bar{1}0]}$[$\xi\xi 0$] mode at $\xi\approx 0.4$, and that bcc remains energetically favorable until about 1200 GPa. This indicates that the fitted potential is suitable for further study of this high-pressure transformation.

V Conclusion

We have developed and applied a novel semi-empirical interatomic potential for tungsten, based on the MEAM and SW formalisms, parameterized using bias-free quintic splines and force-matched to a large database of highly-converged DFT data using an evolutionary global optimization scheme. We have demonstrated accuracy of the fit by reproducing phonon frequencies, compression and thermal-expansion curves, formation energies of unfavorable crystal structures, self-interstitial defects, free surfaces, vacancies, stacking faults and ideal shear at multiple pressures. Transferability of the fitted potential has been demonstrated by description of the high-pressure bcc to fcc phase transformation, dislocation core structure and deformation twinning and detwinning of a tungsten nanorod. Given the accurate description of both deformation twinning and dislocation structure this potential is more suitable than previous models for studying their interplay. Accuracy of elastic and vibrational properties at high pressures will enable quality shock simulations, and the combination of accurate free-surfaces and non-equilibrium crystal structures should produce reliable descriptions of tungsten nanostructures.

ACKNOWLEDGMENTS

Many of the $ab$-$initio$ calculations in our fitting database were performed by Dr. Michael Fellinger. This work was supported by the US Department of Energy under Contract No. DE-FG02-99ER45795 and used computational resources provided by the Ohio Supercomputer Center.