Force-matched embedded-atom method potential for niobium

Large-scale simulations of plastic deformation and phase transformations in alloys require reliable classical interatomic potentials. We construct an embedded-atom method potential for niobium as the first step in alloy potential development. Optimization of the potential parameters to a well-converged set of density-functional theory (DFT) forces, energies, and stresses produces a reliable and transferable potential for molecular-dynamics simulations. The potential accurately describes properties related to the fitting data and also produces excellent results for quantities outside the fitting range. Structural and elastic properties, defect energetics, and thermal behavior compare well with DFT results and experimental data, e.g., DFT surface energies are reproduced with less than 4% error, generalized stacking-fault energies differ from DFT values by less than 15%, and the melting temperature is within 2% of the experimental value.

Figure 1

The three cubic splines of the EAM potential. The points are optimized spline knots and the solid lines are cubic polynomials that interpolate between the knots. (a) The pair potential $\varphi$ and (b) the electronic density $\rho$ are functions of the distance $r$ between pairs of atoms. Both of these functions have 17 optimized spline knots and a cutoff radius of $4.75\text{Å}$, which includes first-, second-, and third-nearest-neighbor interactions in bcc Nb. (c) The embedding function $F$ depends on the local electronic density $n$. Eight optimized spline knots parameterize $F$.

Figure 2

Errors between the 1895 DFT forces in the fitting database and the corresponding EAM forces. The relative errors in the forces tend to decrease as the force magnitudes increase. (a) The percent relative deviation of the EAM force magnitudes from the DFT values, versus the magnitudes of the DFT forces. The weighted relative rms deviation of all the forces is 25%. The inset shows the absolute errors in the forces. (b) The angular deviations of the EAM force directions from the corresponding DFT force directions. The weighted average of all the angular deviations is $15.1°$. (c) Histogram of the angular deviations: 81% of the deviations are within the range $0<\theta <30°$, 96% of the deviations are within the range $0<\theta <60°$, and 98% of the deviations are within the range $0<\theta <90°$.

Figure 3

Schematic illustration of the (a) $⟨100⟩$ dumbbell, (b) $⟨110⟩$ dumbbell, (c) $⟨111⟩$ dumbbell, (d) $⟨111⟩$ crowdion, (e) octahedral, and (f) tetrahedral interstitials.

Figure 4

Phonon-dispersion curves along high-symmetry directions. Our EAM results agree well with experiment (Ref. 59) for small wave vectors but the classical EAM potential is unable to fully capture the spectrum. However, the potential improves upon the EAM results of Guellil and Adams (Ref. 8) and Hu et al. (Ref. 12). The DFT phonon results closely match the experimental data over much of the Brillouin zone. The transverse modes in the $\left[\xi 00\right]$ direction show a plateau around $\xi =0.25$ which is an artifact of the interpolation scheme used to generate the curves. The DFT phonon results were provided by Hennig.

Figure 5

(a) Thermal expansion curve. The thermal expansion of our EAM potential agrees closely with experiment (Ref. 60) from 0 K to $T_{\text{melt}}^{\text{exp}}=2742\text{K}$. Our EAM potential slightly overestimates the thermal expansion while the EAM potential of Guellil and Adams (Ref. 8) underestimates it. (b) Pressure versus volume curve. The experimental data are from shock experiments (Ref. 61). For pressures to 75 GPa, our EAM potential agrees well with experimental data and our DFT calculations.

Figure 6

$\gamma$-surface sections in the $⟨111⟩$ direction. The absence of minima indicate that there are no stable stacking faults in the {112} and {110} planes along this direction. The EAM and DFT results agree well, even though no data from configurations with stacking faults is used to construct the potential.

Figure 7

Differential-displacement maps of the (a) degenerate and (b) nondegenerate (or symmetric) core structures found for (1/2)[111] screw dislocations in bcc metals. For the degenerate core, the structures on the left and right have the same energy. (c) In all cases, the core spreads into three {110} planes of the [111] zone.

Figure 8

Differential-displacement maps for the (1/2)[111] screw dislocation core structure produced by our EAM potential. (a) The degenerate core-structure results when the atoms are relaxed. (b) Shear stress applied along [111] produces a net displacement of the dislocation in the $\left(\overline{1}\overline{1}2\right)$ plane. The figure shows the core structure when the shear stress is just above the critical-resolved shear stress for dislocation motion.

Figure 9

Orientation of the maximum resolved shear stress plane (MRSSP) with respect to the $\left(\overline{1}01\right)$ plane. (a) The angle between the MRSSP and the $\left(\overline{1}01\right)$ plane is $\chi$. (b) Due to symmetry, we only need to consider $-30°<\chi <+30°$ (the shaded interval). This range of angles includes the $\left(\overline{2}11\right)$, $\left(\overline{1}01\right)$, and $\left(\overline{1}\overline{1}2\right)$ planes.

Figure 10

The critical resolved shear stress (CRSS) for dislocation motion as a function of MRSSP orientation. The figure shows a deviation from the Schmid law when the angle between the MRSSP and the $\left(\overline{1}01\right)$ plane is greater than about $15°$.

Figure 11

Two-phase melting simulations determine the melting temperature of our EAM potential. Initially, half the simulation cell contains liquid Nb and the other half contains bcc Nb. The liquid region of the simulation cell solidifies below the melting temperature and the solid region melts above the melting temperature. The figure shows the equilibrium volume for each simulation temperature. There is a sharp jump in volume upon melting.

Figure 12

The melting curve of Nb. The points are EAM results and the solid line fits the values. The melting temperature is computed for pressures to 2.5 GPa with an error of $\pm 5\text{K}$ for each temperature. The EAM potential produces a melting temperature of 2685 K at $P=1\text{atm}$. The experimental melting temperature at this pressure is 2742 K. The error between EAM and experiment is only 2%. Melting temperatures have not been measured for higher pressures.

Figure 13

Radial distribution functions (RDF). (a) Nb is bcc at 293 K. (b) Liquid Nb at 2750 K. No experimental data is available, so we compare the liquid result to the RDF from an EAM potential intended for simulating liquid Nb (Ref. 15). The two potentials produce similar results. The first and second neighbor shells in the solid merge into a single peak in the liquid with similar behavior for higher-order neighbor shells.