A Modified Stillinger-Weber Potential for TlBr , and Its Polymorphic Extension

TlBr is promising for γand xradiation detection, but suffers from rapid performance degradation under the operating external electric fields. To enable molecular dynamics (MD) studies of this degradation, we have developed a Stillinger-Weber type of TlBr interatomic potential. During this process, we have also addressed two problems of wider interests. First, the conventional Stillinger-Weber potential format is only applicable for tetrahedral structures (e.g., diamond-cubic, zinc-blende, or wurtzite). Here we have modified the analytical functions of the Stillinger-Weber potential so that it can now be used for other crystal structures. Second, past modifications of interatomic potentials cannot always be applied by a broad community because any new analytical functions of the potential would require corresponding changes in the molecular dynamics codes. Here we have developed a polymorphic potential model that simultaneously incorporates Stillinger-Weber, Tersoff, embedded-atom method, and any variations (i.e., modified functions) of these potentials. We have implemented this polymorphic model in MD code LAMMPS, and demonstrated that our TlBr potential enables stable MD simulations under external electric fields.


Introduction
Thallium bromide (TlBr) has emerged to be one of the most promising semiconductors for γ-and x-radiation detection in recent years, achieving resolution as high as 1% at 662 keV (Shorohov et al., 2009;Hitomi, Matsumoto, Muroi, Shoji, & Hiratate, 2002).To collect the charges created during the radiation events, the material must be subject to an external electric field.Unfortunately, the performance of TlBr degrades under external electric fields after operation times as short as a few hours to a few weeks (Hitomi, Kikuchi, Shoji, & Ishii, 2009).It is unclear if this degradation comes from the structure evolution due to the ionic conduction under the electric fields.As ionic conduction may be affected by defective structures such as the open channels of edge dislocations, molecular dynamics (MD) that allows extended defects to be included in simulated crystals becomes a useful method to study the ionic conduction induced structure evolution of TlBr.Such MD simulations are not yet possible due to the lack of an interatomic potential for the Tl-Br system that has a CsCl type of crystal structure.
Past MD simulations (Zheng-Johansson & McGreevy, 1996;Zheng-Johansson, Ebbsjö, & McGreevy, 1995;Huang, Wei, Chen, & Chen, 2011;Tung, Chang, Hsiung, Chiang, & Li, 2010;Yamamura, Kawasaki, & Sakai, 1999;Chung, & De Leeuw, 2004;Beckers, Van Der Bent, & De Leeuw, 2000;Lindan & Gillan, 1993) mainly focused on ionic conductivity, but not the structural evolution under external electric fields.As a result, these studies typically do not apply external electric fields.Instead, the system is simply annealed at a sufficiently high temperature to cause thermally activated diffusion of ions.The trajectories of ions are then used to calculate diffusion coefficients of ions.These diffusion coefficients are in turn used to calculate ionic conductivity through Nernst-Einstein relation.Our interest is to understand the structural evolution under external electric fields.Hence, we will apply an external electric field to accelerate the ionic motion, therefore leading to sufficient structural changes needed for the analysis within the short time scale of MD simulations.This, however, requires a robust interatomic potential that can at least maintain the crystal structure at a large external field.
While relatively complex many-body potential (Marrocchelli, Madden, Norberg, & Hull, 2011) can be applied, it is not uncommon to model ionic materials with simple force fields such as pair potentials (Zheng-Johansson, Ebbsjö, & McGreevy, 1995;Huang, Wei, Chen, & Chen, 2011;Tung, Chang, Hsiung, Chiang, & Li, 2010;Yamamura, Kawasaki, & Sakai, 1999), or pair potentials plus angular energy penalty interactions (Chung, & De Leeuw, 2004;Beckers, Van Der Bent, & De Leeuw, 2000).The angular energy penalty approach is very similar to the widely-used Stillinger-Weber (SW) potential (Stillinger & Weber, 1985), and has the advantage to stabilize relatively complex crystal structures that the pair potentials cannot.For example, SW potentials use a parabolic energy penalty term to penalize non-tetrahedral bond angles.As a result, SW potentials have been successfully applied to tetrahedral structures such as diamond-cubic, zinc-blende, and wurtzite Note 1 .Depending on parameterization, SW potentials can also be used for fcc elements (Zhou et al., 2013;Zhou, Jones, Kimmer, Duda, & Hopkins, 2013;Zhou, Jones, Duda, & Hopkins, 2013).However, the conventional SW potentials have two limitations (Zhou et al., 2013;Zhou, Jones, Kimmer, Duda, & Hopkins, 2013;Zhou, Jones, Duda, & Hopkins, 2013): (1) they significantly overestimate the elastic constants of closely packed (e.g., fcc) elements; and (2) they usually cannot be used for other (non-tetrahedral) structures such as sc elements, and NaCl and CsCl compounds.Note that although most elements do not exhibit the lowest energy for the sc crystal, the sc crystal can have a near-lowest energy that cannot be captured by SW potentials.Potentials capable of prescribing a low energy for the sc structure, therefore, improve upon SW potentials in terms of the general energy trends when a variety of configurations are considered (Zhou, Foster, Van Swol, Martin, & Wong, 2014).
Many literature potentials are constructed using particular analytical functions.These potentials can be easily improved if alternative functions are used.For example, Rockett (Wang & Rockett, 1991) used alternative functions in the Tersoff (Tersoff, 1989;Albe, Nordlund, Nord, & Kuronen, 2002) potential format to better treat short-and long-range interactions in covalent systems, and we used alternative functions in the Tersoff potential to improve its prediction on thermal conductivity (Zhou & Jones, 2011).Although alternative functions do not change the potential format and should, therefore, be easily applied, such an approach has not been widely used because it does require new molecular dynamics codes to be developed for each modified function.
With the recognition that the problems described above can limit the atomistic level studies of emerging important materials, the objective of the present work is fourfold: (1) modify the SW potential so that it can better describe elastic constants of elements and be applicable to a wide range of crystal structures including sc, NaCl, and CsCl; (2) develop a polymorphic potential model that incorporates simultaneously SW potential, modified SW potential, Tersoff potential (Tersoff, 1989;Albe, Nordlund, Nord, & Kuronen, 2002), modified Tersoff potential (Wang & Rockett, 1991), and embedded-atom method (EAM) potential (Daw & Baskes, 1984); (3) implement this polymorphic potential model in the public MD code LAMMPS (Plimpton, 1995) so that future development of any alternative functions for these types of potentials no longer requires modification of the molecular dynamics codes; and (4) parameterize the modified SW potential for TlBr system and demonstrate the utility of the resulting potential under external electric fields.

Modified SW Potentials
In SW potential (Stillinger & Weber, 1985), the total energy of a system of N atoms is expressed as Comment: In Eqs.(1), ( 2), (3), ( 9), (15), ϕ should all be changed to φ where i 1 , i 2 , …, i N is a list of neighbors of atom i, θ jik is the bond angle formed by atoms j and k at the site of atom i, φ R,IJ (r ij ) and φ A,IJ (r ij ) are, respectively, pairwise repulsive and attractive functions, u IJ (r ij ) is another pair function for the three-body term, g JIK (cosθ jik ) is an angular energy penalty function, and subscripts i,j,k and I,J,K indicate, respectively, the atoms and the species of the atoms (note that three bodies JIK and KIJ are equivalent).The original SW potentials significantly overestimate the elastic constants of elements because the φ R,IJ (r) and φ A,IJ (r) functions used in these potentials do not allow independent adjustment of bond energy and its second derivative (Zhou et al., 2013).Here we propose to use modified Morse's functions capable of independent change of bond energy and its second derivative (Zhou et al., 2013) to represent φ R,IJ (r) and φ A,IJ (r): where E b,IJ , r 0,IJ , α IJ , and β IJ are four pair dependent parameters, and f c,IJ (r ij ) is a cutoff function.Note that the parameters introduced here have physical meanings: E b,IJ and r 0,IJ correspond respectively to the equilibrium bond energy and bond length, and α IJ >> β IJ control the curvature of the bond energy at the equilibrium bond length.The cutoff function f c,IJ (r ij ) is expressed as: .Note that the cutoff function approximately equals one at r ≤ r s,IJ and equals zero at r = r c,IJ (i.e., r c,IJ is the cutoff distance).Hence, multiplying any potential function with this cutoff function does not affect significantly the potential function at small distances but allows the potential function to be smoothly cut off at r c,IJ .Such a cutoff method is superior to the spline approach used by the Tersoff potential as the latter does not have continuous second and higher order derivatives.A general exponential decay function is used to represent the u IJ (r) function: where γ IJ is a pair parameter.
For the angular function, SW potentials use a parabolic energy penalty to the non-tetrahedral angle g JIK (cosθ jik ) = (cosθ jik -cosθ 0,JIK ) 2 where the parameter cosθ 0,JIK is fixed at the tetrahedral bond angle cosθ 0,JIK = -1/3.This imposes two constraints: not only the favorable bond angle is fixed at the tetrahedral angle, but also the favorable bond angle does not depend on species I, J, and K.The g JIK (cosθ jik ) function can be made more general by treating cosθ 0,JIK as fitting parameters that depend on I, J, and K.Even so, the resulting function is still not fully flexible, for instance, it does not have a scaling factor, and its second derivative with respect to cosθ jik is fixed at 2. In addition, the function is symmetric at cosθ jik = cosθ 0,JIK , and the energy penalty does not saturate (i.e., the absolute slope increases when the angle deviates from cosθ 0,JIK ).These cause difficulties for capturing the angular function derived from quantum mechanical theories (Pettifor et al., 2004).Here we consider a new angular function g JIK (cosθ jik ) = λ JIK •{1-exp[-ξ 0,JIK •(cosθ jik -cosθ 0,JIK ) 2 ]}, where λ JIK is three-body dependent scaling factor and ξ 0,JIK is another three-body dependent parameter.Note that for λ JIK = 1, ξ 0,JIK = 1 and cosθ 0,JIK = -1/3, the modified energy penalty is equivalent to the parabolic function near cosθ jik = -1/3 because the first term of the Taylor series of the modified function expanded at cosθ jik = cosθ 0,JIK is in fact the parabolic function.However, replacing the leading term with a full series does allow the energy penalty curvature to be adjusted through an added parameter ξ 0,JIK and the value of function to be saturated (while retaining the "penalty" effect, i.e., the function minimizes at cosθ jik = cosθ 0,JIK and monotonically increases when cosθ jik deviates from cosθ 0,JIK ).An even further flexible function will be to penalize the energy when the bond angle deviates from three independent values, which also results in asymmetric minimums.Based on this consideration, we propose a general angular function as: where λ JIK , ξ 0,JIK , ξ 1,JIK , ξ 2,JIK , cosθ 0,JIK , cosθ 1,JIK , cosθ 2,JIK , κ 1,JIK , and κ 2,JIK are all three-body dependent parameters.Note that when the parameters are given, the denominator in Equation ( 6) is essentially a normalization constant so that g JIK (cosθ jik =1) = λ JIK .It can be seen that when κ 1,JIK = κ 2,JIK = 0, Equation (6) penalizes the energy when the bond angle deviates from a single value cosθ jik = cosθ 0,JIK as in the conventional SW potential.Otherwise Equation ( 6) can penalize the energy when the bond angle deviates from three values cosθ jik = cosθ 0,JIK , cosθ jik = cosθ 1,JIK , cosθ jik = cosθ 2,JIK .Equations ( 1) -( 6) fully define our modified SW (MSW) potential.

Polymorphic Potential Model
Any improved interatomic potentials will not be applied unless molecular dynamics codes are available to run them.To provide potential developers with a great flexibility for modifying the interatomic potentials without worrying about MD codes, we have constructed a polymorphic potential model.In this model, the energy of the system is expressed as where δ ij is Kronecker delta (i.e., δ ij = 1 when i = j and δ ij = 0 when i ≠ j), η ij is an indicator of the potential type that can be set to either ) and V IJ (r ij ) are two pair functions, and F IJ (X ij ) is a function of a local variable X ij that will be discussed below.It can be seen that when η ij = δ ij , the summation in equation ( 7) excludes the self-interaction term i = j and is therefore over all pairs of different atoms.When


, which can be used to incorporate the embedding energy of the embedded-atom method as will be clear later in this section.The variable X ij essentially accounts for the environment surrounding the ij bond, and is defined as where P IJ (Δr jik ) is a function of weighted difference between atomic spacing r ij and r ik , which is written as Δr jik = r ij ξ IJ •r ik with the weighting factor ξ IJ being either 0 or 1 to include or exclude r ik , W IK (r ik ) is another pair function, and G JIK (θ jik ) is a three-body function of bond angle θ jik .It can be seen that this polymorphic interatomic potential model is fully defined when the indicators η ij and ξ IJ , and the six functions U IJ (r), V IJ (r), P IJ (Δr), W IJ (r), F IJ (X), and G JIK (θ) (for all the species I, J, K = 1, 2, …) are given.Note that these six functions can all be supplied as one-dimensional tables and can therefore be implemented in MD codes using cubic spline interpolation and/or extrapolation.As a result, users can easily perform simulations using different potentials by tabulating these functions (in a MD read-in table file) accordingly.For instance, the polymorphic potential reduces to our MSW potential if we tabulate the functions according to: where φ R,IJ (r), φ A,IJ (r), u IJ (r), and g JIK (cosθ) are defined by Equations ( 2), (3), ( 5) and ( 6).The polymorphic potential reduces to a conventional SW (Stillinger & Weber, 1985) potential if we tabulate the functions according to: where A IJ , B IJ , ε IJ , σ IJ , λ IJ , γ IJ , a IJ , p, and q are the normal parameters for the SW potential as described above.The polymorphic model represents Tersoff types of potential (Tersoff, 1989;Albe, Nordlund, Nord, & Kuronen, 2002) if we set where f c,IJ (r ij ) is a cutoff function defined as , , , , and D e,IJ , S IJ , r e,IJ , β IJ , μ IJ , γ IJ , c IJ , d IJ , h IJ , r s,IJ and r c,IJ are all pairwise parameters.The polymorphic potential can also represent the Rockett-Tersoff potential (Wang & Rockett, 1991) where f c,1,IJ (r) is a cutoff function similar to equation ( 12) but operates at a different cutoff range: and A IJ , B IJ , λ 1,IJ , λ 2,IJ , λ 3,IJ , β IJ , n IJ , c IJ , d IJ , h IJ , r s,1,IJ , and r c,1,IJ are all pairwise parameters.To use the polymorphic model for the embedded-atom method potential (Daw & Baskes, 1984), we can simply set: where φ IJ (r) is a pair function, f J (r) is an atomic electron density function, F I (X) is the embedding energy function, and X is used to represent electron density (X = ρ).We have implemented this polymorphic potential model in the public parallel MD code, LAMMPS (Plimpton, 1995).We plan to release the code in LAMMPS package after the present paper is published.Interested readers can contact us to get the code prior to the formal release of the code.

Parameterization
To enable structure evolution of ionic materials to be studied within the short time scales of MD simulations, we attempt to accelerate ionic migration by external electric fields.Such electric fields are approximated by applying opposite external forces to cations and anions.To ensure that the interatomic potential is robust enough to maintain the equilibrium crystal structure of TlBr when Tl+ and Br -ions are subject to large external forces, we require that the potential not only best captures the key experimental properties of the observed Tl, Br, and TlBr phases (e.g., Tl-hcp, TlBr-CsCl, etc.), but also predicts the crystalline growth of the ground state structures during MD simulations of growth (e.g., vapor deposition).Note that a potential is said to be incapable of crystalline growth simulations only when it predicts amorphous growth at all combinations of temperatures and growth rates possible with the MD simulations.At low temperatures and high growth rates where earlier adatoms are buried by later adatoms before the surface reaches a low energy configuration, amorphous growth is the correct prediction.The growth simulation tests are important because they sample a variety of configurations (at the growth surface) not considered a prior.If any of the random nuclei formed on the growth surface has a lower energy than the growth crystal, the simulations is likely to always give an amorphous growth regardless of temperature and growth rate.Hence, crystalline growth provides strong validation that the growth crystal has the lowest energy compared to any other configurations.When the equilibrium TlBr crystal has the lowest energy, large external forces can be applied to Tl+ and Br -ions without causing phase transformation.Hence, the growth simulation capability is essential for our applications.
We proceed by parameterizing first the Tl and Br potentials, and then the TlBr potential at the known Tl and Br parameters.The observed room-temperature equilibrium phases are hcp for Tl (Donnay & Ondik, 1973), diatomic (Br 2 ) liquid for Br, and CsCl for TlBr (Donnay & Ondik, 1973).Note that solid Br has an orthorhombic crystal structure; however, Br 2 liquid is the stable room-temperature phase.Like many other potentials, our MSW model is not intended to capture the Br 2 molecules (it is possible to capture the Br 2 molecules, but this does not necessarily result in a better potential).On the other hand, our density function theory (DFT) calculations, employing the optB86b-vdW functional (see Appendix for computational details), indicated that the sc Br phase has a lower cohesive energy than dc, bcc, fcc, and hcp phases.Hence, we target sc as the lowest energy lattice phase for Br (at 0 K) while at the same time ensure that the stable Br phase at room temperature is liquid.Note that if a correct negative heat of formation of the TlBr compound is captured, elemental phases do not form in MD simulations under stoichiometric conditions.Because they do not form, particular elemental structures are not important for studying stoichiometric compounds.
Lattice constants, cohesive energies, and elastic constants for the model lowest energy lattices (Tl-hcp, Br-sc, and TlBr-CsCl) are fitted under the constraints that the energies of all the other phases (e.g., dc, bcc, fcc, NaCl, wz, etc.) are higher than the targeted lowest energy phases.For the experimentally observed structures such as Tl-hcp and TlBr-CsCl, the experimental lattice constants (Donnay & Ondik, 1973), cohesive energies (Barin, 1993), and elastic constants (Simmons, 1965;Morse & Lawson, 1967) are directly used as the target values for the fitting.For the phases that are not observed (e.g., Br-sc), the available experimental properties of other phases and DFT results are used to guide the selection of the target values.In particular, DFT results may be quite different from the experimental values.For example, the hcp Tl cohesive energy obtained from experiments and DFT calculations is -1.85 and -2.40 eV/atom respectively.As a result, we scale the DFT results so that for the observed structures, the scaled DFT results match the experiments.
The software package Mathematica (Wolfram, 2004) is used to perform the parameterization.To promote global optimization, four different numerical optimization routines, namely a conjugate gradient method (Hestenes & Stiefel, 1952), the downhill simplex method of Nelder and Mead (Olsson & Nelson, 1975), a genetic algorithm (Storn & Price, 1997), and a biased random walk (simulated annealing) (Kirkpatrick, Gelatt, & Vecchi, 1983), are all used to determine the parameters that minimize the weighted mean-square deviation between the target and predicted properties.The goal to capture the crystalline growth is more challenging, requiring a highly iterative parameterization process.After each fitting iteration, the four sets of parameters from the four optimization routines are tested for vapor deposition simulations.If the potential does not pass these tests, the entire process is repeated with an appropriate adjustment of parameters bounds, target structures and target properties.The iterations continue until one of the four optimization routines results in a satisfactory set of potential parameters.The MSW potential thus determined are listed in Tables I and II for two-body, and three-body parameters respectively.
Comment: In Tables 3 and A-I, use horizontal lines to separate Tl, Br, and TlBr data.In the present version, it is not clear which structures are for Tl, Br, and TlBr.
Table III indicates that of the lattice structures explored, our MSW potential captures the Tl-hcp and TlBr-CsCl lattice crystals as the lowest energy phases, in agreement with experiments that Tl-hcp and TlBr-CsCl phases are observed at room temperature.Table III also indicates that our MSW potential captures the lowest energy for the Br-sc crystal, which is really the target of the MSW model.We cannot include the Br 2 liquid observed in experiments in Table III.Instead, we will show below through growth simulations that our potential gives Br liquid as the most stable phase at room temperature.Here we can see that in addition to capturing appropriate lowest energy phases, our MSW model also reproduces the experimental cohesive energies of Tl-hcp and TlBr-CsCl.Our energy of -0.827 eV/atom for Br-sc is not unreasonable compared with the experimental cohesive energy of -1.134 eV/atom for a more stable Br 2 -liquid, especially considering that this difference between MSW and experiments is not much more significant than that between DFT and experiments, not to mention that the elemental Br phase does not occur for our intended problem.Note that MSW model also accurately captures the experimental lattice constant of TlBr-CsCl.The lattice constants of Br-sc and Tl-hcp as determined by our MSW is not far off from the experimental or DFT values either (in case of Tl-hcp, for example, while the lattice constant a is slightly under-estimated, the lattice constant c is slightly over-estimated resulting in a good description of atomic volume).
Capturing non-tetrahedral crystals (such as hcp, sc, CsCl) as the lowest energy phases proves that the MSW potential is more flexible than the conventional SW potentials.SW potentials are not designed for capturing property trends of a variety of metastable structures.The purpose of listing some selected metastable phases in Table III is to show that the equilibrium phases have lower energies than these metastable phases.Interestingly, however, our MSW potential reproduces the DFT order of Tl structures with increasing (more negative) cohesive energies as in dc → sc → bcc → fcc.Furthermore, our MSW model correctly captures the lowest energy for Tl-hcp whereas the DFT does not.This observation suggests that MSW potential can possess good transferability to a variety of other phases.Usually such transferability is only possible with more sophisticated potentials.Since our MSW is designed for studying the TlBr bulk, the transferability to many metastable phases is not relevant and so we do not exploit this further for the Br and TlBr phases which would necessarily sacrifice the properties of the equilibrium phases.
Finally, we point out that Table III only explores limited number of phases and therefore does not prove that the predicted lowest energy phases indeed have the lowest energies as compared to any other configurations.As will be discussed below, we will prove this rigorously through vapor deposition simulations.

Elastic Constants and Melting Temperature of the Observed Tl and TlBr Phases
Elastic constants and melting temperature of the experimentally observed Tl-hcp and TlBr-CsCl phases are calculated, and the results are summarized in Table IV along with the experimental values.Here the melting temperature calculations follow the same method described previously (Zhou et al., 2013;Zhou, Foster, Van Swol, Martin, & Wong, 2014;Ward, Zhou, Wong, Doty, & Zimmerman, 2012;Ward, Zhou, Wong, Doty, & Zimmerman, 2012;Ward, Zhou, Wong, & Doty, 2013).It can be seen that our MSW model reproduces the experimental elastic constants for TlBr-CsCl.The model elastic constants Tl-hcp are generally lower than the experimental values.Note that we could have fit exactly the experimental elastic constants for Tl-hcp.However, this would result in a melting temperature that is significantly above the experimental value.On the other hand, the conventional SW potentials tend to significantly overestimate elastic constants for closely packed elements and hence the ability to prescribe low elastic constants is significant for the MSW model.We get the melting temperature of Tl-hcp reasonably close to the experimental value.Despite a prolonged effort, our current parameterization of the MSW potential still significantly overestimates the melting temperature of TlBr-CsCl.We can better capture the melting temperature by reducing the elastic constants.Considering that our objective is to study the behavior of the TlBr-CsCl crystal at ambient temperature, we chose to capture exactly the elastic constants at the cost of overestimating the melting temperature.We feel that this problem cannot be resolved by parameterization alone, and further modifications of the MSW format are needed in future efforts in order to capture both elastic constants and melting temperature.

D. Point Defects
Native point defects in TlBr-CsCl, including Tl vacancy V Tl , Br vacancy V Br , Tl at Br antisite Tl Br , Br at Tl antisite Br Tl , Tl interstitial between Tl sites Tl i,1 , Tl interstitial between Br sites Tl i,2 , Br interstitial between Tl sites Br i,1 , and Br interstitial between Br sites Br i,2 , are all studied.In particular, intrinsic defect energies ' D E are calculated as (Zhang & Northrup, 1991;Northrup & Zhang, 1993) where E, N Tl , and N Br are total energy, number of Tl atoms, and number of Br atoms of the system containing the defect, and E TlBr , E Tl , and E Br are cohesive energies (per atom unit) for the lowest energy phases of TlBr, Tl, and Br respectively.The results obtained from MSW and DFT calculations are summarized in Table V.It can be seen that MSW potential predicts a lower energy for Tl vacancy than for Br vacancy, a low energy for Br at Tl antisite than for Tl at Br antisite, and a lower energy for Br i,2 interstitial than for both Tl interstitials.All these are in good agreement with the DFT calculations.Our MSW potential indicates that the Br i,1 interstitial has a lower energy than both Tl interstitials, which differs from the DFT result that the Br i,1 interstitial energy is only lower than the Tl i,2 interstitial energy.Again here we only compare the trends, and do not make conclusions on the absolute values due to the lack of experimental data.

Molecular Dynamics Vapor Deposition Verifications
As mentioned above, only when a potential captures the crystalline growth during MD vapor deposition simulations will it capture the lowest energy phase and be robust enough to allow applications of high external electric fields.Here we perform vapor deposition simulations to validate that our TlBr MSW potential captures the crystalline growth of the lowest energy Tl-hcp and TlBr-CsCl crystals.In addition, we will demonstrate that our potential gives the crystalline growth of Br-sc at a low temperature (this is the lowest energy "model" phase as opposed to the orthorhombic experimental lowest energy phase) but gives a liquid Br structure as the stable phase at room temperature.

Tl-hcp Growth
For Tl-hcp growth, an initial substrate of an hcp crystal containing 1008 Tl atoms with 28 ( ) 2 1 10 layers in the x direction, 9 ( ) 0002 layers in the y direction, and 8 ( ) 0 110 layers in the z direction is used.Here layers refer to crystallographic planes so that one (0001) layer is equivalent n (000n) layers etc.The substrate temperature is set at T = 300 K by assigning velocities to atoms according to the Boltzmann distribution.During simulations, the bottom (-y) 2 (0002) layers are held fixed to prevent crystal shift upon adatom impact on the top surface.The next 3 (0002) layers are isothermally controlled at the substrate temperature.This leaves the top 4 layers free where the motion of atoms is solely determined by Newton's law.Injection of Tl adatoms from random locations far above the surface simulates the growth.All adatoms have an initial far-field incident kinetic energy E i = 0.05 eV and an incident angle θ = 0 o (i.e., the moving direction is perpendicular to the surface).The adatom injection frequency is chosen to give a deposition rate of R = 2.5 nm/ns.To approximately maintain a constant thickness of the free surface region, the isothermal region expands upward during simulations.Since surface roughness might develop, the isothermal region expands at about 80% of the surface growth rate to ensure that the upper boundary of the isothermal region never exceeds the surface valley locations.Fig. 1 shows the resulting configuration obtained after 0.66 ns deposition, where the original substrate is shaded in purple.It can be seen that the MSW potential correctly captures the crystalline growth of the Tl-hcp phase.This strongly validates that Tl-hcp has the lowest energy at room temperature as compared to any other configurations.

Br-sc Growth
MD simulations are also performed to grow Br on an Br-sc substrate.A sc crystal containing 1008 Br atoms with 14 ( ) 100 layers in the x direction, 9 ( ) 010 layers in the y direction, and 8 ( ) 001 layers in the z direction is used as the initial substrate.During simulations, the bottom 2 (010) layers are held fixed.The next 4 (010) layers are controlled at the desired growth temperature.Using the same approach as described above, the growth simulation is performed at two substrate temperatures T = 150 K and T = 300 K, an incident energy E i = 0.05 eV, an incident angle θ = 0 o , and a deposition rate R = 2.5 nm/ns.The resulting configurations obtained after 1.20 ns deposition are shown in Figs.2(a) and 2(b) respectively for the 150 K and 300 K temperatures.It can be seen from Fig. 2(a) that the sc-crystalline growth is achieved with our potential, strongly validating that the Br-sc crystal has the lowest (free) energy at 150 K as compared to any other configurations.Fig. 2(b), on the other hand, shows that amorphous growth is achieved at 300 K.This strongly validates that Br exhibits a liquid type phase at room temperature.020 layers are used to control the growth temperature.Following the same approach as described above, the growth simulations are performed at a substrate temperature T = 700 K, an incident energy E i = 0.05 eV, an incident angle θ = 0 o , a deposition rate R = 2.75 nm/ns, and a stoichiometric vapor flux ratio Tl:Br = 1:1.Fig. 3 shows the system configuration obtained at 0.78 ns deposition time.It is seen again that our MSW potential correctly captures the crystalline growth of the equilibrium (CsCl) phase of TlBr.In particular, the randomly injected Tl and Br atoms are reconstructed correctly to their corresponding sublattices.Because the potential captures the crystallization from a rather stochastic vapor phase, it enables robust simulations of structural evolution of TlBr-CsCl under conditions where the presence of both dislocations and external electric fields may induce configuration disorders if the potential has any deficiencies in capturing the lowest energy phase.

Molecular Dynamics Structural Evolution Under Fields
Our objective is to allow direct MD simulations of structural evolution of TlBr-CsCl crystal under an external electric field.The external electric field can be simulated by applying opposite biased forces to Tl and Br atoms.Note that in our model, we do not directly address charges.This is a reasonable approximation because charges only give two forces: the Coulomb forces between atoms, and biased forces under external fields.The Coulomb forces between atoms are digested into the interatomic potential in our model.This allows us to use biased forces to independently simulate the external electric field.
The simulations may become challenging at large external electric fields because the biased forces may induce phase transformation when the potential does not capture the lowest energy for the equilibrium phase.Here we demonstrate two cases to demonstrate that our potential allows simulations to be performed when atoms are subject to large forces of ±0.4 eV/Å.These forces correspond to a high electric field of 4×10 6 V/mm assuming Tl and Br atoms adopt full charges of ±1 e. Again note that atoms are not strictly point charges and atoms in perfect bulk may not be subject to big forces from the external field.Nonetheless, the model becomes robust if we pass this test.
A TlBr-CsCl crystal containing 16128 Tl atoms and 16128 Br atoms with 84 ( ) 110 layers in the x direction, 24 ( ) 110 layers in the y direction, and 32 ( ) 002 layers in the z direction is used.To prevent system from shifting, the bottom region of about 10 Å wide is fixed.To remove the boundary effects, periodic boundary conditions are used in all three coordinate directions.An MD simulation is then performed at 1200 K (a homologous temperature 0.832 T m ) with a biased force of ±0.4 eV/Å using the NVT ensemble (i.e., number of atoms, volume, and temperature are all constant).
In the first case, we assume that the system contains a pair of Tl and Br vacancies by removing a Tl and a Br atom far away from the fixed region.Comparison of atomic configurations between a time span of 0.06 ns is shown in Fig. 4. Fig. 4 verifies that our potential allows for stable TlBr-CsCl MD simulations to be performed with a high biased force.Interestingly, we found that Tl and Br vacancies are not very mobile even at the high biased force and temperature.In fact, during the 0.06 ns span, Tl vacancy jumped by one lattice spacing whereas the Br vacancy did not jump.In the second case, we assume that the system contains a pair of Tl and Br interstitials by inserting a Tl atom and a Br atom at locations far away from the fixed region.Comparison of atomic configurations between a time span of 0.12 ns is shown in Fig. 5. Fig. 5 again verifies that our potential allows stable TlBr-CsCl MD simulations to be performed with a high biased force.Unlike the vacancy case, we found that Tl and Br interstitials are very mobile, and both interstitials moved a significant distance during the 0.12 ns time span.The cases shown in Figs. 4 and 5 are designed to validate the robustness of our potential, but not to launch a thorough study of structure deterioration under external fields.We are currently using our model to study how defects (dislocations, vacancies, interstitials) evolve under external electric fields, and how diffusion of various species contribute to mass transport and changes of stoichiometry of the material.These studies are beyond the scope of this paper and will therefore be published separately.

Conclusions
Herein we have developed a new modified Stillinger-Weber potential.Unlike the conventional SW potential that significantly overestimates the elastic constants of closely packed elements and are limited mainly to tetrahedral structures, our modified potential can capture very low elastic constants for closely packed elements and can be used for many non-tetrahedral crystal structures.We have parameterized the modified SW potential for TlBr.Through rigorous vapor deposition simulation tests, we have demonstrated that our potential captures the experimental properties of the observed Tl and TlBr phases, and predicts the Br liquid to be the most stable phase at room temperature.Moreover, we have demonstrated that our potential is robust enough for challenging simulations of structure evolution of TlBr crystals under very high external electric fields.Test simulations indicate that interstitials migrate much faster than vacancies.We have also developed a polymorphic potential model and implemented it in the public MD code, LAMMPS.This essentially enables future material research to be performed at a higher fidelity level because improved potentials no longer require modification of the MD codes and therefore can immediately be utilized by a broader materials and physics community.
to the inclusion of exact Fock exchange.The results obtained using the optB86b-vdW functional are in relative good agreement with those obtained using the HSE06+D method.Moreover, the lattice constants obtained for TlBr-CsCl and Tl-hcp agree very well with experiment (~1%), and even slightly better than those obtained using the HSE06 and HSE06+D methods.Note, the PBE functional was also considered but the predicted lattice constants were in greater disagreement with experiment.Regarding the prediction of the cohesive energies, the HSE06 functional performs the best.This is attributed to a decrease in self-interaction error as a result of using a functional that incorporates exact Fock exchange.We have decided to included only the HSE06 values in the main text for internal consistency and simplicity; however, we have included the results obtained using the different methods for completeness and for interested readers.
4)where r s,IJ and r c,IJ (r s,IJ << r c,IJ ) are two independent pair parameters, and ζ IJ and υ IJ are two dependent pair parameters

Figure 1 .
Figure 1.Vapor deposited hcp Tl film obtained from MD simulations

Figure 2 .
Figure 2. MD simulations of Br growth on a Br-sc substrate at a temperature of (a) 150 K and (b) 300 K. Our potential prescribes a liquid Br as the most stable phase at room temperature

Figure 3 .
Figure 3. Vapor deposited CsCl phase of a TlBr film obtained from MD simulations

Figure 4 .
Figure 4. Structure evolution of a TlBr-CsCl crystal containing vacancies at a temperature of 0.832 T m and a biased force of ±0.4 eV/Å: (a) Starting time of observation; and (b) 0.06 ns later.The TlBr-CsCl crystal remains intact at the large electric field.During the 0.06 ns span, Tl vacancy jumped by one lattice spacing whereas the Br vacancy did not jump

Figure 5 .
Figure 5. Structure evolution of a TlBr-CsCl crystal at a temperature of 0.832 T m and a biased force of ±0.4 eV/Å: (a) Starting time of observation; and (b) 0.12 ns later.The TlBr-CsCl crystal remains intact at the large electric field.During the 0.12 ns span, Tl and Br interstitials migrated significant distances

Table 1 .
Two-body parameters of MSW potential (length in Å and energy in eV)

Table 3 .
Lattice constants a (for dimer, a refers to the dimer bond length) and c (Å), and cohesive energy E c (eV/atom), obtained from the MSW potential, DFT calculations, and experiments for selected material structures

Table 4 .
Elastic constants C 11 , C 12 , C 13 , C 33 , C 44 , bulk modulus B (eV/Å 3 ), and melting temperature T m (K) of the observed Tl and TlBr crystalline phases, obtained from the MSW potential and experiments

Table A -
I. Lattice constants a (for dimer, a refers to the dimer bond length) and c (Å), and cohesive energy E c (eV/atom) determined using optB88b-vdw, HSE06, and HSE06+D functionals Br Tl Br Br Tl Tl i,1 Tl i,2 Br i,1 Br i,2