In this test case, let us consider a body-centered cubic (bcc) crystal of Tantalum (Ta), described by the Finnis-Sinclair (FS) potential.31 The calculations are performed using the MD++ program. The source code and the input files for this and subsequent test cases in this chapter can be downloaded from http://micro. stanford. edu/wiki/Comprehensive_ Nuclear_Materials_MD_Case_Studies.

The cut-off radius of the FS potential for Ta is 4.20 A. To avoid interaction between an atom with its own periodic images, we consider a cubic simulation cell whose size is much larger than the cut-off radius. The cell dimensions are 5[100], 5[010], and 5[001] along x, y, and z directions, and the cell contains N = 250 atoms (because each unit cell of a bcc crystal contains two atoms). PBC are applied in all three directions. The experimental value of the equi­librium lattice constant ofTa is 3.3058 A. Therefore, to compute the equilibrium lattice constant of this potential model, we vary the lattice constant a from 3.296 to 3.316 iA, in steps of 0.001 JA. The potential energy per atom E as a function of a is plotted in Figure 5. The data can be fitted to a parabola. The

ao (A)

location of the minimum is the equilibrium lattice constant, a0 = 3.3058 /A. This exactly matches the experimental data because a0 is one of the fitted parameters of the potential. The energy per atom at a0 is the cohesive energy, Ecoh = —8.100 eV, which is another fitted parameter. The curvature of parabolic curve at a0 gives an estimate of the bulk modulus, B = 197.2 GPa. However, this is not a very accurate estimate of the bulk modulus because the range of a is still too large. For a more accurate determination of the bulk modulus, we need to compute the E(a) curve again in the range of |a — a01 < 10 4 A. The curvature of the E(a) curve at a0 evaluated in the second calcu­lation gives B = 196.1 GPa, which is the fitted bulk modulus value of this potential model.31

When the crystal has several competing phases (such as bcc, face-centered cubic, and hexagonal — closed-packed), plotting the energy versus volume (per atom) curves for all the phases on the same graph allows us to determine the most stable phase at zero temperature and zero pressure. It also allows us to predict whether the crystal will undergo a phase transition under pressure.32

Other elastic constants besides B can be computed using similar approaches, that is, by imposing a strain on the crystal and monitoring the changes in potential energy. In practice, it is more convenient to extract the elastic constant information from the stress-strain relationship. For cubic crystals, such as Ta considered here, there are only three independent elastic constants, Cn, C12, and C44. C11 and C12 can be obtained by elongating the simulation cell in the x-direction, that is, by changing the cell length into L =(1 + exx) • L0, where L0 = 5a0 in this test case. This leads to nonzero stress components sxx, ayy szz, as computed from the Virial stress formula [19], as shown in Figure 6 (the atomic velocities are zero because this calculation is quasistatic). The slope of these curves gives two of the elastic con­stants C11 = 266.0 GPa and C12 = 161.2 GPa. These results can be compared with the bulk modulus obtained from potential energy, due to the relation B = (Cn + 2 C12) / 3 = 196.1GPa.

C44 can be obtained by computing the shear stress axy caused by a shear strain £xy Shear strain exy can be applied by adding an off-diagonal element in matrix H that relates scaled and real coordinates of atoms.

The slope of the shear stress-strain curve gives the elastic constant C44 = 82.4 GPa.

In this test case, all atoms are displaced according to a uniform strain, that is, the scaled coordinates of all atoms remain unchanged. This is correct for simple crystal structures where the basis contains only one atom. For complex crystal structures with more than one basis atom (such as the diamond-cubic structure of silicon), the relative positions of atoms in the basis set will undergo additional adjustments when the crystal is subjected to a macroscopically uniform strain. This effect can be captured by performing energy minimization at each value of the strain before recording the potential energy or the Virial stress values. The resulting ‘relaxed’ elastic constants corre­spond well with the experimentally measured values, whereas the ‘unrelaxed’ elastic constants usually overestimate the experimental values.