Development and application of a Ni-Ti interatomic potential with high predictive accuracy of the martensitic phase transition

Phase transitions in nickel-titanium shape-memory alloys are investigated by means of atomistic simulations. A second nearest-neighbor modified embedded-atom method interatomic potential for the binary nickel-titanium system is determined by improving the unary descriptions of pure nickel and pure titanium, especially regarding the physical properties at finite temperatures. The resulting potential reproduces accurately the hexagonal-close-packed to body-centered-cubic phase transition in Ti and the martensitic B2-B19′ transformation in equiatomic NiTi. Subsequent large-scale molecular-dynamics simulations validate that the developed potential can be successfully applied for studies on temperatureand stress-induced martensitic phase transitions related to core applications of shape-memory alloys. A simulation of the temperature-induced phase transition provides insights into the effect of sizes and constraints on the formation of nanotwinned martensite structures with multiple domains. A simulation of the stress-induced phase transition of a nanosized pillar indicates a full recovery of the initial structure after the loading and unloading processes, illustrating a superelastic behavior of the target system.


I. INTRODUCTION
Shape-memory alloys are a class of materials with the property of recovering their original shape upon heat treatment (shape-memory effect) and of sustaining large elastic strains (superelasticity). These unique properties make shapememory alloys into widely used functional materials in many applications [1]. Among the shape-memory alloys discovered so far, nickel-titanium (NiTi) shape-memory alloys with equiatomic or nearly equiatomic compositions have received great attention owing to their excellent mechanical properties, corrosion resistance, biocompatibility, and their ability to transform close to room temperature [2]. In NiTi alloys, the shape-memory effect and superelasticity result from the reversible temperature-or stress-induced martensitic phase transition between cubic B2 (austenite) and monoclinic B19 (martensite), respectively [2].
Recently, the application of NiTi shape-memory alloys has been extended to micro-and nanoelectromechanical systems (MEMS/NEMS) [3,4]. To this end, the properties of miniaturized NiTi shape-memory alloys such as nanosized wires, pillars, and particles are of great interest as they can differ significantly from their bulk counterparts. Experimental studies focusing on the exceptional characteristics of phase transitions in miniaturized NiTi alloys are therefore ongoing [3,4]. To supplement experiments, the theoretical study of phase transitions by means of atomistic simulations such as molecular dynamics (MD) is highly desirable to provide a detailed understanding of the underlying mechanisms.
To enable large-scale MD simulations of phase transitions in NiTi shape-memory alloys, the availability of a reliable interatomic potential is of crucial importance. Surprisingly, presently available interatomic potentials cannot sufficiently well reproduce the phase transitions in the NiTi system. Interatomic potentials were developed by Farkas et al. [5] based on the embedded-atom method (EAM) [6] and by Lai and Liu [7] based on the Finnis-Sinclair model [8]. These potentials are not able to reproduce the observed phase transitions [9], as they were developed by focusing on the properties of the ternary Ni-Ti-Al alloy [5] and amorphous alloys [7]. Two further potentials were developed by Saitoh et al. [10] and Ishida and Hiwatari [11] based on the modified EAM (MEAM) [12] with a focus on the phase transitions in NiTi. However, these potentials were not found to satisfactorily reproduce the reversible temperature-and stress-induced phase transitions and crystallography of related phases [9].
Recently, the Finnis-Sinclair potential by Lai and Liu [7] was independently modified by Mutter and Nielaba [9] and Zhong et al. [13]. The authors [9,13] clarified the occurrence of a reversible temperature-induced phase transition; however, the motif of the low temperature martensite structure was clearly different from the experimentally reported B19 structure. As a direct consequence of the wrongly predicted martensite structure, twinning-being the most important deformation mechanism in shape-memory alloys [2,3]-cannot be properly reproduced. For example, one potential [9] predicts the martensitic transition without the occurrence of any twinning, and the other potential [13] predicts an unphysical twinning behavior with a negative twin boundary energy.
In the present paper, we have developed a potential based on the second nearest-neighbor (2NN) MEAM model [14][15][16] and applied it to study the phase transitions of NiTi shape-memory alloys. The 2NN MEAM potential parameterization has been selected here because of the notorious difficulties experienced in previous studies to properly describe the complex low temperature martensite structure (B19 ) with simple potential parameterizations. The characteristics of the B19 structure, a monoclinic angle of ≈98°and a shuffle of Ni and Ti atoms at the faces of the unit cell [2], lead to a strong directionality of the atomic bonds in the interatomic potentials. This feature is difficult to capture with simple models, whereas the 2NN MEAM model implicitly provides angle-dependent terms to reflect the directionality of atomic bonds.
In addition to the proper selection of the potential model, the selection of an appropriate optimization method of the potential parameters is necessary considering difficulties in previous studies based on the MEAM model [10,11]. Previous MEAM potentials were developed mostly focusing on energies of a few configurations at 0 K, and a transferability to properties at finite temperatures was never tested. In the present paper, the force-matching method proposed by Ercolessi and Adams [17] is utilized for the optimization. This method considers forces and energies related to various atomic configurations, including configurations at finite temperatures from density functional theory (DFT). It has been previously reported that this method provides robust potentials for various applications at finite temperatures [18,19].
To develop a potential for the binary Ni-Ti system, potentials for the pure Ni and Ti systems are necessary because the 2NN MEAM description of a binary system is based on the constituent unary potentials. 2NN MEAM potentials for the pure Ni and Ti systems have been developed previously by Lee et al. [20] and Kim et al. [21], respectively, but these potentials are not appropriate for investigating phase transitions at finite temperatures. The most serious problem of the previously developed Ti potential [21] is that it cannot reproduce the phase transition between body-centered-cubic (bcc) Ti (β-austenite phase) and hexagonal-close-packed (hcp) Ti (α-martensite phase), which is closely related to the transition mechanism in the NiTi shape-memory alloy. Both systems (pure Ti and NiTi) have very similar structures of entropically stabilized austenite phases (bcc in Ti and B2 in NiTi) and similar transition paths related to the imaginary phonon modes of the austenite phases [22][23][24][25]. In addition, the previous potentials [20,21] less accurately reproduce physical properties at finite temperatures than those at 0 K. To overcome these deficiencies, the development of potentials in the present paper proceeds in a systematic manner: Accurate potentials for the pure Ni and Ti systems are developed first, and the development of a binary potential for the Ni-Ti system is subsequently addressed.

A. Construction of a DFT database
For the potential development based on the force-matching method, a DFT database of atomic forces and energies related to various atomic configurations needs to be prepared. To ensure the transferability of the potential to a wide variety of atomistic situations, configurations resulting from various temperature and strain conditions as well as various defect configurations are included in our DFT database, as listed in Table I. Corresponding physical target properties are compiled in Table II.
The DFT calculations were performed using the VASP code [26][27][28] and the projector-augmented wave (PAW) method [29] within the Perdew-Burke-Ernzerhof generalized gradient approximation (GGA) [30,31] for the exchangecorrelation functional. In the PAW potential for Ti, 3p electrons were treated as part of the valence. A cutoff energy of 340 eV for the plane-wave basis set and the Methfessel-Paxton smearing method with a width of 0.1 eV were used. A k-point mesh of 19×19×19 was selected for the face-centered cubic (fcc) primitive unit cell, and the corresponding k-point density was employed for other unary structures, binary compound structures, and supercells. Magnetism was included by considering spin-polarized calculations for pure Ni, binary compounds, and Ni-rich supercells.
To obtain the equilibrium lattice constant and bulk modulus, we employed the Birch-Murnaghan equation of state [32][33][34] fitted to a volume range from 0.95V 0 to 1.05V 0 (V 0 = equilibrium volume). The elastic constants were calculated by applying strains from −1 to +1% to a cell in the equilibrium structure. In all defect calculations, atomic positions were relaxed at a constant volume and cell shape with the convergence criteria for energy and forces set to 10 −6 eV and 10 −2 eV/Å, respectively. For the calculation of the vacancy migration energy of the pure metals and the solute migration energy of the binary solid solutions, the nudged elastic band (NEB) method with the climbing-image extension [35,36] was used to find a suitable saddle-point configuration. Phonon calculations were performed using the supercells listed in Table I (108 atoms for fcc Ni and hcp Ti and 128 atoms for bcc Ti and B2 NiTi). The calculations were based on the direct force constant approach [37], as implemented in the Phonopy code [38,39] with the convergence criteria for energy and forces set to 10 −8 eV and 10 −4 eV/Å, respectively.
To obtain configurations at finite temperatures, ab initio MD simulations [26] were conducted. Following the concept of the upsampled thermodynamic integration using Langevin dynamics (UP-TILD) method [40], in a first step these calculations were performed based on relatively low DFT convergence parameters to provide an efficient sampling of the configuration space. The cutoff energy was set to the default value of the PAW potential from the VASP library and a single k-point was used. The MD simulations were performed for a total of 1000 steps with a time step of 1.5 fs. In a following step, uncorrelated MD snapshots were extracted from the MD trajectories and recalculated with a higher cutoff energy and denser k-point mesh to determine accurate forces and energies for the fitting process.

B. 2NN MEAM potential
The MEAM interatomic potential was proposed by Baskes [12] as an extension of the EAM potential by introducing additionally angular dependent terms. MEAM potentials are well suited for simulations of multicomponent systems composed of elements with different ground states because they can describe a wide range of phases (fcc, bcc, hcp, diamond-structured, and even gases) using a common mathematical formalism. The MEAM description was improved by Lee and Baskes (2NN MEAM) [14] to partially consider 2NN interactions, thereby overcoming some critical shortcomings of the original MEAM approach.
Within the MEAM approach, the total energy of a system is approximated as where F i is the embedding energy as a function of background electron densityρ i . Further, S ij and φ ij (R ij ) are the screening function and the pair interaction between atoms i and j TABLE I. Atomic configurations entering the DFT database used for fitting (marked by the asterisk in the Fitting column) and testing (unmarked in the Fitting column) the present MEAM potential. In the Stability column, "stable" indicates that the structure is reported in an equilibrium phase diagram, and "metastable" indicates that the structure has been reported as a metastable phase. Other phases are labeled as "hypothetical." The Strain column indicates the strain applied to the supercells, where H , O, and M denote hydrostatic, orthorhombic, and monoclinic strains, respectively.
where the atomic electron density scaling factor ρ 0 and the decay lengths β (h) are adjustable parameters and r e is the nearest-neighbor distance in the equilibrium reference structure.
To compute the total energy in Eq. (1), a functional form of the pair interaction φ ij (R ij ) is also necessary. In the MEAM approach, the value of the pair interaction φ ij (R ij ) is computed not from a specific functional form but from a known value of the total energy and the embedding function for an atom in an equilibrium reference structure. Here, the equilibrium reference structure is defined as a structure where individual atoms are sitting on exact lattice points. The total energy per atom for the equilibrium reference structure is obtained from the zero-temperature universal equation of state by Rose TABLE II. Physical properties of binary NiTi and constituent pure Ni and Ti systems considered in the present paper. Properties closely related to atomic configurations from the DFT database explicitly used in the parameter optimization process (marked by an asterisk in Table I) are in italics. The cohesive energy of fcc Ni and hcp Ti is written in bold italics because the rescaling procedure was used to adjust it to the experimental value (see Sec. II C). Properties used only for testing the transferability of the potential are written in regular, upright font, with bold upright indicating that only the experiment is presently available for comparison.
where d is an adjustable parameter, and Here, r e is the equilibrium nearest-neighbor distance, E c is the cohesive energy, B is the bulk modulus, and is the equilibrium atomic volume of the reference structure. The value of the pair interaction is evaluated from the known values of the total energy per atom and the embedding energy as a function of the nearest-neighbor distance. While the original MEAM considers only first nearest-neighbor interactions by using a strong many-body screening function [42], the 2NN MEAM partially considers also 2NN interactions by adjusting the screening parameters (C min , C max ) so that the many-body screening becomes less severe. A detailed formulation of the unary 2NN MEAM formalism is available in the literature [14][15][16]20,21].
To describe an alloy system, pair interactions between different elements need to be determined. For this purpose, a similar technique is employed as for the pair interactions of TABLE IV. Optimized 2NN MEAM potential parameter set for the binary Ni-Ti system. The following quantities are dimensionful: the enthalpy of formation of the reference structure (B2 NiTi) E f (eV/atom), the equilibrium nearest-neighbor distance r e (Å), and the bulk modulus B (10 12 dyne/cm 2 ).

Parameter
Selected value pure elements. A binary reference structure, where one type of atom has only the same type of atoms as second-nearest neighbors, is chosen. The total energy per atom of the reference structure is computed using the universal equation of state. Then, the pair interactions between different elements are obtained from the known values of the total energy per atom and the embedding energy of the reference structure. A detailed formulation of the binary 2NN MEAM formalism is available in the literature [43].

C. Optimization of potential parameters
To describe a pure element using the 2NN MEAM potential formalism, 14 independent parameters are necessary. Four of these parameters [the cohesive energy (E c ), the equilibrium nearest-neighbor distance (r e ), the bulk modulus (B) of the reference structure, and the adjustable parameter d] are related to the universal equation of state. Seven further parameters [the decay lengths (β (0) , β (1) , β (2) , β (3) ) and the weighting factors (t (1) , t (2) , t (3) )] are related to the electron density. The parameter A belongs to the embedding function, and the parameters C min and C max are responsible for the many-body screening. Details are given in Refs. [14][15][16]20,21]. For describing a binary alloy, 13 independent parameters are necessary, in addition to the constituent unary parameters. Four of these parameters [E c , r e , B, and d] are related to the universal equation of state. The atomic electron density scaling factor ρ 0 belongs to the electron density, and the remaining eight parameters (four C min and four C max ) are responsible for the many-body screening. Details are given in Refs. [16,43].
We determined the unary potential parameters by fitting to the DFT database of energies and forces introduced in Sec. II A. A genetic algorithm that allows for an efficient parameter search in high-dimensional spaces was used as an optimizer. The optimization proceeded iteratively. First, a reference structure, a radial cutoff distance, and fitting weights were specified. The optimization algorithm then adjusted the parameters to minimize the weighted error between the DFT database and the corresponding values produced by the parameters. If the fitting error was too large and the potential failed to reproduce physical properties satisfactorily, the parameters were refitted based on different reference structures, radial cutoff values, fitting weights, and/or DFT databases with added or removed configurations. This optimization process was repeated until reliable potentials were obtained.
The optimization of the pure Ni potential was conducted by fitting forces and energies related to various configurations of the fcc, bcc, hcp, and liquid phases, as listed in Table I. During the optimization, a rescaling of the structural energies was performed in order to match the experimental cohesive energy of the fcc phase (4.450 eV [20]), which we consider TABLE V. Calculated bulk and defect properties of pure Ni using the present 2NN MEAM potential, in comparison with experimental data, DFT data, and previous MEAM calculations by Lee et al. [20]. The following quantities are listed: the cohesive energy E c (eV/atom); the lattice constant a (Å); the bulk modulus B and the elastic constants C 11 , C 12 , and C 44 (10 12 dyne/cm 2 ); the structural energy differences E (eV/atom); the vacancy formation energy E vac f (eV); the vacancy migration energy E vac m (eV); the activation energy of vacancy diffusion Q vac (eV); the divacancy formation energy E div f (eV); and the surface energies E surf (erg/cm 2 ) for the orientations indicated by the superscript. more reliable than the corresponding DFT value (4.842 eV). The use of this rescaling procedure has only a minor effect on the overall fitting since it is related to an arbitrary choice of the reference state. The optimization of the pure Ti potential was started with forces and energies related to configurations of the hcp (α), bcc (β), and hexagonal (ω) phases, which are stable at ambient and high temperature and high pressure conditions, respectively [44]. The rescaling of the structural energies to obtain the experimental cohesive energy of hcp (α) Ti (4.870 eV [21]) was also performed. During the optimization process, we found that a stabilization of the hexagonal (ω) phase, a ground state phase predicted by DFT [24], over the hcp (α) phase always resulted in a significant worsening of the properties of the hcp (α) phase in particular in a negative vacancy formation energy. Considering the importance of the hcp phase and of the hcp-bcc (α-β) phase transition for the present paper, a final optimization was conducted using a reduced fitting weight for the hexagonal (ω) phase. Table III presents the finally determined potential parameter sets for pure Ni and Ti, which were the basis for the following optimization of the binary Ni-Ti potential. The binary Ni-Ti system shows several stable intermetallic compound phases at Ni 3 Ti, NiTi (equiatomic), and NiTi 2 compositions [45] and metastable compound phases at Ni 3 Ti 2 [46] and Ni 4 Ti 3 [47] compositions. During the optimization, we found that the binary parameters introduced in the 2NN MEAM model do not provide the necessary flexibility to accurately describe all these phases by the force-matching process. In particular, the small structural energy differences at the equiatomic composition predicted by DFT (on the order of a few meV/atom) could not be reproduced by any of our fitted potentials, which instead predicted significant energy differences after atomic relaxation. We could trace back this problem to the pair interaction description, which, as explained above, is determined not from a specific functional form but from the universal equation of state. The universal equation of state was originally introduced into the MEAM formalism [12] in order to allow a straightforward extension to multicomponent systems avoiding in particular overfitting. However, it also inevitably restricts the flexibility of the overall fitting procedure resulting in the observed discrepancies. A long-term route to overcome these restrictions would be a search for more flexible functional dependencies of the pair interactions. In the present paper, we have followed a more pragmatic route, as described in the following.
The binary potential developed in the present paper is mainly intended to properly describe the martensitic transformations close to the equiatomic composition, and to this end an accurate description of the corresponding small structural energy differences is essential. We have therefore included only configurations of these equiatomic compounds into the fitting procedure (cf . Tables I and II). Further, we excluded forces from the fitting procedure, and instead we focused on the following target properties: lattice constants, bulk moduli, TABLE VI. Calculated bulk and defect properties of pure Ti using the present 2NN MEAM potential, in comparison with experimental data, DFT data, and previous MEAM calculations by Hennig et al. [24] and by Kim et al. [21]. The following quantities are listed: the cohesive energy E c (eV/atom); the lattice constants a and c (Å); the bulk modulus B and elastic constants C 11 , C 12 , C 13 , C 33 , and C 44 (10 12 dyne/cm 2 ); the structural energy differences E (eV/atom); the vacancy formation energy E vac f (eV); the vacancy migration energy E vac m (eV); the activation energy of vacancy diffusion Q vac (eV); the divacancy formation energy E div f (eV); and the surface energies E surf (erg/cm 2 ) for the orientations indicated by the superscript. In-basal plane and out-of-basal plane vacancy formation and migration energies are designated by "in" and "out," respectively. enthalpies of formation, and the monoclinic angle, which were all obtained by fully including atomic relaxations for each configuration. More specifically, we employed the following stepwise fitting procedure in order to capture the important characteristics of the phase transition. Initially, a large number of candidate parameter sets was prepared. Then, using all these parameter sets, one of the target properties was computed, and only parameter sets yielding the target value were selected for further consideration and calculation of the other target properties. In this manner, the candidate parameter sets were sequentially reduced. The sequence of target properties was arranged as follows: (1) enthalpy of formation; (2) monoclinic angle; (3) lattice constants; (4) bulk moduli [with (1) to (4) corresponding to 0 K]; and, finally, (5) phase stabilities at finite temperatures. Table IV presents the determined potential parameter set for the binary Ni-Ti system. After the optimization of the binary potential, a radial cutoff distance of 5.0Å was confirmed to be sufficiently large to reproduce various physical properties as well as the phase transitions of pure metals and binary alloys. Therefore, the simulations in the following sections were performed based on this radial cutoff distance.

III. ACCURACY AND TRANSFERABILITY OF THE DEVELOPED 2NN MEAM POTENTIAL
In this section, the accuracy and transferability of the developed potential are examined by comparing them to a 134107-7 wide variety of physical properties from DFT and experiment. The corresponding 2NN MEAM calculations were performed using the LAMMPS code [48]. If not specifically designated as MD simulations, all calculated values represent results of molecular statics simulations, and the number of atoms in the supercells is at least 4000. Cell dimensions and individual atomic positions were allowed to fully relax. MD simulations were performed using a time step of 2 fs, the Nosé-Hoover thermostat [49,50], and the Parrinello-Rahman barostat [51] for controlling temperature and pressure, respectively.
Generally, the physical properties calculated in this section can be divided into two groups. The first group comprises properties that are closely related to the atomic configurations used in the parameter optimization (see configurations marked with an asterisk in Table I and properties written in italics in Table II). The comparison of these properties indicates the accuracy of the fitting. The other group comprises properties that were not used directly in the parameter optimization process (unmarked configurations in Table I and properties in regular, upright font in Table II). The comparison of these properties indicates the transferability of the developed potential. For the properties in the first group, DFT values are available and can be directly used for the comparison. Many of the properties in the second group correspond to finite temperature conditions and cannot be easily obtained by DFT calculations. For these properties, we use experimental data for comparison (bold properties in Table II).

Bulk and defect properties at T = 0 K
We focus first on bulk and defect properties, most of which were included in the fitting optimization (cf. Table II). Corresponding results are compared with DFT, experimental data, and previous potentials in Tables V and VI. The experimental cohesive energy is exactly reproduced for both Ni and Ti due to the employed rescaling procedure (Sec. II C). The other bulk properties should be compared with the DFT values, and we observe a satisfactory agreement, except for the structural energy differences of pure Ti. In particular, the present Ti potential cannot reproduce the stability of the hexagonal (ω) Ti phase because the optimization was TABLE VII. RMS errors for energies and forces of pure Ni and Ti with respect to the DFT database (cf. Fig. 1). Values using the present potential are compared to those using the previous potentials by Lee et al. [20] and Kim et al. [21]. conducted with a reduced fitting weight of this phase to avoid the worsening of the properties of the hcp (α) phase, as discussed in Sec. II C. This means that phase transitions related to the hexagonal (ω) phase cannot be studied by the present Ti potential. The calculated vacancy formation and activation energies by the present potentials agree well with their DFT counterparts except for the divacancy formation energy of pure Ti, which shows an underestimation of about 25%. Similar deviations can be observed also for the surface energies where the largest discrepancy in Ni is 19% for the [111] surface and in Ti 36% for the [1120] surface. These deviations reflect the fact that the surface energies were not explicitly included in the fitting process; this should be considered in future studies when applying our potentials to corresponding simulations.

RMS error
Overall, the accuracy of our Ni and Ti potentials in describing bulk and defect properties can be considered as satisfactory. However, comparing to the accuracy of the previous potentials [20,21] (see Tables V and VI), we observe no significant improvement; thus the refitting of the pure potentials is unjustified at this stage. The justification will be provided in the following three sections (Sec. III A 2-III A 4) where we stepwise extend the performance evaluation.

Statistical quality of energies and forces
We first extend the evaluation to a large set of configurations that were obtained from MD simulations at finite temperatures. To this end, we investigate the statistical correlation between the energies and forces from the present/previous MEAM potentials and the DFT database. The MEAM calculations were performed using atomic configurations, as obtained from DFT without further atomic relaxations. Figure 1 shows scatter plots for the statistical correlation, and Table VII shows   error for pure Ti is generally higher than that for pure Ni, reflecting the difficulty in describing configurations related to the hexagonal (ω) phase. (2) The correlation with the DFT values is considerably better for the present potential than for the previous ones. In particular, the RMS error for both the energies and forces is reduced by a factor of two, as evidenced by Table VII. This gives a first indication of the improved performance of our potential at finite temperatures. However, the configurations entering the statistical analysis here have been used explicitly in the optimization of the potential, and we thus have to further extend the evaluation.

Quality of perturbative forces: Phonon-dispersion relations
We now come to the group of properties that were not explicitly included in the fitting process but are especially important for the application of the developed potentials over a wide range of temperature conditions. In this section, we investigate phonon dispersion relations, which were obtained by perturbative calculations from small displacements around the equilibrium structure. This provides a very sensitive measure of the fitted atomic interactions. Figure 2 shows phonon spectra along high-symmetry directions in the Brillouin zone of fcc Ni, hcp Ti (α), and bcc (β) Ti for the MEAM potentials, DFT, and the experiment [53,64,65]. For fcc Ni [ Fig. 2(a)], the present potential and the previous one by Lee et al. [20] closely reproduce the overall dependence of the phonon branches with a small overestimation of the DFT/experimental frequencies at the high symmetry points.
For the hcp (α) Ti phase [ Fig. 2(b)], the acoustic branches by the present potential are in better agreement with experiment and DFT than those by the previous potential by Kim et al. [21], in particular along the [ξ 00] direction. The optical branches are likewise better described by the present potential, but the full complexity of the DFT/experimental dispersion cannot be captured as observed for example for the highest frequency band around .
For the bcc (β) Ti phase, imaginary phonon branches predicted by DFT reflect the instability of this phase at low temperatures. It was previously reported that the imaginary phonon branch along the [ξξ0] direction corresponds to the Burgers mechanism of the hcp-bcc (α−β) transition [24,66] and the imaginary phonon branch along the [ξξξ] is responsible for the (111) plane collapse mechanism of the bcc-hexagonal (β-ω) transition [24]. As shown in Fig. 2(c), the previous potential by Kim et al. [21] cannot reproduce any of the imaginary phonon branches, which is likely related to the failure of this potential [21] to reproduce the corresponding phase transitions. In contrast, the present potential can be expected to reproduce the hcp-bcc (α-β) transition correctly because it can successfully reproduce the imaginary phonon branch along the [ξξ0] direction. The difference in the long wavelength limit, i.e., that our potential predicts mechanical stability, whereas DFT shows a mechanical instability, is an artefact of the employed supercell size. For larger supercell sizes, DFT likewise predicts a mechanically stable system, i.e., positive frequencies in the long wavelength limit as explicitly investigated in Ref. [24]. The reproducibility of the hcp-bcc (α-β) phase transition by our potential will be explicitly tested in Sec. IV B below. In contrast, the bcc-hexagonal (β-ω) transition is not feasible by both potentials because the stability of the hexagonal (ω) phase is not properly reproduced, as discussed above.

Quality of thermodynamics: Expansion, heat capacity, and melting
We finally evaluate the performance of the developed pure Ti and Ni potentials in describing thermal properties such as the thermal expansion coefficient, the specific heat, the melting point, the enthalpy of melting, and the volume change upon melting. The thermal expansion coefficient and the specific heat were calculated based on the quasiharmonic approximation and using full MD simulations. The latter were performed using an isobaric-isothermal ensemble (NPT) at zero pressure and 300 K. The melting temperature was calculated using the interface method, which utilizes a simulation cell consisting of solid and liquid phases in contact with each other. The enthalpy of melting and the volume change upon melting were calculated using an NPT ensemble at zero pressure and at the obtained melting temperature.
The calculated thermal properties are compared with experimental data in Fig. 3 and Table VIII. For the thermal expansion coefficient and the specific heat of pure Ni, the present potential and the previous potential by Lee et al. [20] show no clear difference in the reproducibility of experiment. For the melting temperature of pure Ni, the result by the present potential shows better agreement with the experimental value. It is about 10% higher than the experimental value, while that by the previous potential by Lee et al. [20] shows about 16% overestimation. The calculated enthalpy of melting and the volume change upon melting of pure Ni by the present potential are also in better agreement with experimental values than those by the previous potential by Lee et al. [20].
In the case of the thermal properties of pure Ti, the present potential shows a considerably improved high temperature phase stability. The previous potential by Kim et al. [21] does not properly stabilize the bcc (β) phase at higher temperatures, as we confirmed by our own calculations. In fact, this potential [21] stabilizes the hcp (α) phase until the melting temperature of 1706 K; thus, it cannot be employed for finite temperature simulations of the bcc phase and in particular of the bcc to liquid transition. This incorrect behavior at finite temperatures is related to the failure of describing the imaginary phonon branch along the [ξξ0] direction, as discussed above. In contrast, our potential correctly predicts the transition to the bcc (β) phase (discussed in detail in Sec. IV B), and it shows a reasonable melting temperature of 1651 K, which is an underestimation of the experimental melting temperature by only 15%. The corresponding melting enthalpy likewise shows an acceptable agreement with the experiment (17% underestimation).

Summary for pure Ni and Ti potentials
Overall, we have shown that the present potentials for the pure Ni and Ti systems reproduce a wide range of fundamental physical properties. In particular, the good transferability of the potentials to properties relating to lattice vibrations and thermodynamic properties is crucial for finite temperature applications. A striking improvement over the previously available potential [21] is found for Ti, where the phase transition sequence hcp to bcc to liquid is now properly described. Based on these good pure-element descriptions, the potentials can be confidently extended towards the binary Ni-Ti system.

Accuracy and transferability for equiatomic NiTi compounds
At the equiatomic NiTi composition, several compounds were proposed in previous experiments [2,[71][72][73] and DFT calculations [74][75][76] (see Fig. 4). The cubic B2 and the monoclinic B19 are the structures of the experimentally well investigated austenite and martensite phases, respectively [2]. The B19 structure differs from the B2 structure by a monoclinic distortion (β = 97.9 • [71]) and additionally by a shuffling of Ni and Ti atoms on the (110) B2 plane. The orthorhombic B19 structure can be regarded as an intermediate structure between the B2 and B19 structures as it shows the shuffles similar to the B19 structure but without the monoclinic distortion. The B19 structure is experimentally observed when Cu is alloyed to the NiTi alloy [73]. The orthorhombic B33 structure, which also shows atomic shuffling albeit with 134107- 11  FIG. 4. (Color online) Atomic structures of the cubic B2, orthorhombic B19, monoclinic B19 , and orthorhombic B33 NiTi phases with the monoclinic angle (β) and the lattice constants (a and c) indicated. Ni and Ti atoms are represented by blue (dark gray) and orange (gray) balls, respectively. a higher monoclinic angle (β = 107 • ) [74][75][76], is a T = 0 K ground state structure predicted by DFT. So far, this is only a hypothetical phase not observed by any experiment, possibly because of kinetic constraints, as it is not geometrically related to the shape-memory transformation.
In Table IX, the 0 K properties of the equiatomic compounds calculated by the present binary potential are compared with DFT and experimental data. These properties were used in the parameter optimization process (cf . Tables I and II), and the comparison reflects therefore the accuracy of the potential fitting. We observe that our potential reproduces well the lattice constants, the atomic volume, the enthalpy of formation, and the bulk modulus of the compounds. It also successfully reproduces the DFT trend in the formation enthalpies, including a slightly more negative formation enthalpy of the B33 structure than that of the B19 structure (6 meV/atom). To evaluate the transferability of our binary potential to finite temperatures at the equiatomic composition, we focus on the temperature dependence of the lattice parameters of the B19 structure, which were not included into the fitting. The results  (blue circles in Fig. 5) indicate a nearly constant a lattice constant, an increase in the b lattice constant, and a decrease in the c lattice constant with increasing temperature. All these dependencies are in excellent agreement with the experimental data near room temperature (black squares in Fig. 5), indicating a good performance of the present potential at the equiatomic composition.
The phonon dispersion of the B2 structure, which contains critical information about the martensitic transition, is shown in Fig. 6 (blue lines). As for the pure elements, the B2 phonons were obtained by perturbative calculations from small displacements around the equilibrium structure, thereby providing a very sensitive measure. Overall, we observe a reasonable qualitative agreement with DFT (orange lines), although some of the complex dependencies of the optical branches cannot be well reproduced. However, the most important feature for the present purpose is the imaginary, acoustic phonon branch along the [ξξ0] direction, which is critical for the B2-B19 transition mechanism, similarly as the imaginary phonon branch along the [ξξ0] direction in pure bcc Ti is responsible for the hcp-bcc (α−β) transition (see Sec. III A 3). The present binary potential successfully reproduces this imaginary phonon branch and thus can be expected to properly capture the martensitic transformation B2-B19 . The details of this transformation will be investigated in Sec. IV C, and it will be shown that the transformation is indeed well described.
While the imaginary phonon branch along the [ξξ0] direction is well described, there is a disagreement in the branch along the [ξξξ] direction. DFT predicts an instability, whereas our potential predicts a stable phonon branch. Any attempts to avoid this disagreement without modifying the unary parameters to better describe the hexagonal (ω) Ti phase have failed. It therefore seems very plausible that the deficiency of the binary potential is closely related to the wrong description of the hexagonal (ω) Ti phase, which itself originates in the wrong description of the imaginary phonon branch along the [ξξξ] direction of the bcc (β) Ti phase [ Fig. 2(c)] by the pure Ti potential. An improvement in the properties of the hexagonal (ω) Ti phase always resulted, however, in a worsening of the properties of the hcp (α) phase and in a wrong description of the imaginary phonon branch along the [ξξ0] direction of the bcc (β) Ti and the B2 NiTi phases. Considering the intention of the present paper to develop a binary NiTi potential especially for the martensitic transformation, an accurate description of the hexagonal (ω) Ti phase and of the imaginary phonon branch along the [ξξξ] direction of both the bcc (β) Ti and the B2 NiTi phases was sacrificed.

Transferability to other compounds, solid solutions, and the liquid phase
We finally evaluate the performance of the present binary MEAM potential in describing stable and metastable compounds at various compositions other than the equiatomic composition. None of these compounds have been included into the optimization process (see the discussion at the end of Sec. II C); therefore, this evaluation is a stringent test of transferability. Table X reveals that except for the c lattice constant of the metastable Al 3 Os 2 type structure (33% overestimation), the lattice parameters, the enthalpy of formation, and the bulk modulus are reasonably well reproduced for the various structures. From these results, it cannot be yet concluded whether the experimentally reported stable phases (DO 24 and E9 3 phases) [45] are indeed correctly predicted as stable phases by the present potential due to a possible decomposition. We have therefore extended the test and calculated the stability of several additional, hypothetical compound phases with our potential and with DFT to determine the ground-state convex hull. Figure 7 shows the results, including the formation energies of the equiatomic compounds. The DFT convex hull is generally well reproduced by our potential, and the ground state structures match the ones in DFT. There are, however, some quantitative differences. At the Ni 3 Ti composition, DFT predicts that the enthalpy of formation of the L1 2 structure is more positive by 17.2 meV/atom than that of the experimentally reported DO 24 structure. The present potential reproduces this trend, but the difference is tiny (0.05 meV/atom); it is thus not clear whether the DO 24 structure is indeed stable at finite temperatures.

134107-13
We have therefore extended the test to finite temperatures to investigate whether unwanted phase transitions, e.g., to hypothetical compounds, occur. For that purpose, NPT ensemble MD simulations with gradually increasing temperatures were performed using initially the experimentally reported stable structures (DO 24 and E9 3 ) of the Ni 3 Ti and NiTi 2 phases. We could confirm that the Ni 3 Ti phase correctly maintains its  , and the solute migration energy E sol m . The reference states for the dilute heat of solution are fcc Ni and hcp Ti. In the Ni-rich fcc phase, the energies for first and second nearest-neighbor bindings are designated by "1NN" and "2NN" respectively. In the Ti-rich hcp phase, in-basal plane and out-of-basal plane binding and migration energies are designated by "in" and "out" respectively. The defect binding energies are defined such that a positive sign indicates attractive interaction.

Structure
Property original crystal structure (DO 24 ) upon heating to a corresponding overheated melting temperature. However, with the present potential, the NiTi 2 phase does not maintain its original crystal structure (E9 3 ) at finite temperatures, and it decomposes to a partly disordered phase. We attribute the wrong prediction of the NiTi 2 (E9 3 ) phase to the incompleteness of the present potential at the Ti-rich side. This should be kept in mind in future applications of the present potential. Table XI lists calculated physical properties of the two solid solutions (the dilute heat of solution, the vacancy-solute binding energy, the solute-solute binding energy, and the solute migration energy). The DFT calculation predicts negative dilute heats on both sides of the phase diagram, indicating strong bonding between Ni and Ti. This result is consistent with the experimental fact that the Ni-Ti system shows several intermetallic compound phases. The present potential correctly reproduces these trends qualitatively, but there is a rather large quantitative discrepancy in the dilute heat of the Ti-rich solid solution with a strong underestimation of the Ni formation energy (MEAM: −0.786 eV versus DFT: −0.144 eV). In the Ni-rich solid solution, the present potential and the DFT calculation consistently predict a small binding (positive sign) or repulsive (negative sign) interaction between a vacancy and a Ti atom and a strong repulsive interaction between Ti atoms. There is a general agreement with DFT also for the migration energy of a Ti atom toward the adjacent vacancy in the Ni-rich solid solution. In the Ti-rich solid solution, however, the present potential cannot reproduce trends in the binding and the migration properties. We therefore stress again that special attention is necessary in future simulations of the Ti-rich solid solution using the present potential. In contrast, the performance of the present potential on the Ni-rich side is significantly better. To support this statement further, we have investigated also the lattice constant composition dependence in the Ni-rich solid solution, as shown in Fig. 8. The experimental data [82] show an increase in the lattice constant as the concentration of Ti increases, and we observe that the calculated values using the present potential closely follow this trend.
As a last transferability test, we have investigated the composition dependence of the liquid phase. Figure 9 shows the enthalpy of mixing of the liquid phase at 2000 K calculated by the present potential compared to available experimental data [82]. The calculation was performed using NP T ensemble MD simulations with a random distribution of each element at a certain composition. The results by both the present calculation and the experimental data indicate a negative deviation from the ideal mixing with a minimum point at the composition of around 40% Ti. The present potential thus successfully reproduces the liquid phase over a large composition range.

C. Summary of the performance of the present potential
To summarize, the developed potential performs very well in describing equiatomic Ni-Ti compounds, Ni-rich compounds, Ni-rich solid solutions, the binary liquid phase, phonons related to the α-β phase transition in pure Ti, and phonons related to the B2-B19 transition in equiatomic Ni-Ti. It less suited for describing Ti-rich compounds and Ti-rich solid solutions, and the following properties are not captured: dynamical instability along the [ξξξ] direction in pure bcc Ti, the resulting β-ω transition in pure Ti, and the dynamical instability along the [ξξξ] direction in B2 NiTi.

IV. APPLICATIONS OF THE DEVELOPED POTENTIAL
We have shown in Sec. III that the present 2NN MEAM binary potential well describes various fundamental physical properties of the Ni-Ti system as well as those of the pure Ni and Ti systems. In this section, we present several core applications for which our potential is particularly well suited. We investigate the temperature dependence of the monoclinic angle of the B19 phase allowing us to resolve a previous discrepancy between DFT and experiment. Further, we study temperature-and stress-induced phase transitions of pure Ti and of NiTi shape-memory alloys. Corresponding simulations are based on MD within an NP T ensemble using the same conditions for the time step, thermostat, and barostat as in the previous section (Sec. III).

A. Temperature dependence of the monoclinic angle of B19
A closer inspection of Table IX [71]). This discrepancy led the authors of previous studies [75,76] to regard the DFT structure with the higher monoclinic angle as a new, separate phase and to label it B19 . However, the conclusion drawn in the previous studies [75,76] has been based on a comparison of DFT values computed at 0 K and experimental data obtained at room temperature, and it is not clear whether temperature-induced changes might influence the comparison.
An extension of the DFT calculations to finite temperatures could solve this issue but is computationally demanding. Instead, we can utilize for this purpose our potential, which describes very well the finite temperature properties of the B19 phase (Fig. 5). The result for the temperature dependence of the monoclinic angle is shown in Fig. 10 (blue solid curve) in comparison to experiment (black squares). We observe a good agreement in the temperature dependence (cf. dashed curve) with a small constant shift. At room temperature, the present potential predicts an angle of 97.0 • , which is close to the experimental value of 97.9 • [71]. Following the temperature dependence of the shifted curve down to 0 K, we obtain a  [71]) and T = 0K DFT results from Refs. [75] and [76] and from the present paper. The dashed MEAM curve has been shifted on top of the room temperature experiment to emphasize the similar temperature dependence. value of 100.4 • , which falls within the shown DFT values. The differences in the three DFT results and a possible impact of the exchange-correlation functional would require a further, detailed DFT investigation. Nevertheless, our result suggests that the previously reported discrepancy between DFT and experiment is mainly due the different temperature conditions for DFT (0 K) and experiment (room temperature) and that there is no need to introduce a new phase such as B19 .
The result of the present MD simulation might also explain why the ground state structure (B33) predicted by DFT is not observable by experiment. In the DFT study [74] that first proposed the B33 structure as a ground state structure of equiatomic NiTi, the monoclinic angle was artificially constrained to the experimentally reported value of 97.8 • , i.e., the room temperature value. The corresponding DFT calculation, however, was performed at T = 0 K, and this introduces a strain energy penalty because the structure has been deformed from its ground state monoclinic angle. Due to this strain energy, the B19 structure was energetically destabilized over the B33 structure by 8 meV/atom. However, if the energy difference between the B33 and the B19 structures is correctly obtained using a fully relaxed B19 structure with a higher monoclinic angle (100.0 • ) [76], the energy difference is very small (less than 0.5 meV/atom [76]), which is close to the resolution limit of DFT.

B. Temperature-induced phase transition of pure Ti
The temperature-induced phase transition of pure Ti was analyzed because of its similarity to the transition mechanism of the NiTi shape-memory alloy. The corresponding MD simulations were performed starting with a single-crystal supercell of pure Ti. Initially, the bcc structure was equilibrated at 1950 K. The temperature was then gradually decreased to 750 K and increased again to 2250 K with cooling and heating rates of ±0.5 K/ps. Periodic boundary conditions were applied along all three dimensions, and cell dimensions and angles were allowed to relax. Temperature-induced changes in the atomic volume were recorded to observe the occurrence of phase transitions. To analyze finite-size effects, several independent simulations were performed using supercells ranging in size from 462 to 250 000 atoms. Figure 11(a) shows a representative temperature dependence of the atomic volume during cooling and reheating for a supercell with 1848 atoms. As the bcc (austenite) phase is cooled down, it transforms into the hcp (martensite) phase. The discontinuous jump in the volume curve represents the phase transition event, and the corresponding temperature (900 K) is recorded as a martensite start temperature (M s ). The atomic structure of the transformed hcp phase is compared with the initial structure of the bcc phase in the inset of Fig. 11(a). The hcp structure corresponds to a perfect crystal without any sign of defects such as twins. The orientation relation between the initial bcc and the final hcp structure clearly indicates that the phase transition occurs via the Burgers mechanism [66], which has the orientation relation of (110) bcc || (0002) hcp and [111] bcc || [1120] hcp . Upon reheating, the hcp phase transforms back into the bcc phase at a much higher temperature (1880 K) than the M s temperature. This temperature is recorded as the austenite finish temperature (A f ).
The range of calculated transition temperatures (900 K and 1880 K) covers the experimentally reported transition temperature between the bcc and the hcp phases (1155 K [44]). However, a significant amount of thermal hysteresis (≈1000 K) caused by an overheating and undercooling is observed in the simulations, whereas the experimental measurements indicate a very narrow hysteresis (≈10 K [83]). We attribute the large theoretical hysteresis to (1) the absence of heterogeneous nucleation sites such as defects and (2) the prevention of active habit planes as the size of the simulation cell is tiny compared to that of experimental samples. The importance of heterogeneous nucleation sites was reported in previous MD simulations on the phase transition of a NiAl shapememory alloy [84][85][86][87][88]. These works also resulted in a large thermal hysteresis (≈1000 K) if simulations were started using defect-free crystals, in accordance with our present finding. Further simulations considering various kinds of defects, such as a dislocation [85], a grain boundary [84,88], an antiphase boundary [85], and a free surface [85,86], confirmed that such defects can significantly reduce the thermal hysteresis window. Figure 11(a) also illustrates the occurrence of a solid-liquid phase transition upon continued heating of the solid phase to higher temperatures. The overheated solid phase transforms into the liquid phase at a temperature of 2010 K. The amount of overheating is around 360 K, which we obtained by comparing the melting temperature of the overheated solid to the equilibrium melting temperature of 1651 K calculated by the interface method in Sec. III A 4. Figure 11(b) shows transition temperatures calculated using various supercell sizes. As the size of the system increases, the M s temperature decreases and the A f temperature increases, indicating larger undercooling and overheating windows. The reduced amount of undercooling and overheating in smaller sized supercells is possibly caused by an increased correlation in the thermal and stress fluctuations, providing more opportunities for the system to transform. If the number of atoms in a system is larger than 10000, the M s temperature converges to a value of 820 K. The converged value of the A f temperature cannot be accurately established because it exceeds the melting temperature of the overheated hcp phase.

C. Temperature-induced phase transition of a NiTi shape-memory alloy
We now turn to the temperature-induced phase transition of a NiTi shape-memory alloy. Corresponding MD simulations were performed starting with a single-crystal supercell of equiatomic NiTi. Initially, the B2 structure was equilibrated at 550 K. The temperature was then gradually decreased to 10 K and increased again to 550 K with cooling and heating rates of ±0.5 K/ps. Periodic boundary conditions were applied along all three dimensions, and cell dimensions were allowed to relax. To investigate the effect of mechanical constraints on 134107-17 the martensitic transformation, two sets of simulations, with and without the relaxation of the cell angles (angles fixed at 90 • for the latter), were independently performed. To analyze finite size effects, several independent simulations were performed using supercells ranging in size from 384 to 978 432 atoms. Figure 12(a) shows a representative temperature dependence of the atomic volume for a cell with 3072 atoms. The high-temperature B2 austenite phase transforms into the low-temperature B19 martensite phase during cooling, and the martensite phase transforms back into the austenite phase during reheating indicating M s and A f temperatures, respectively. The volume change due to the martensitic transition is positive, as also expected from experimental work [71]. The amount of the room-temperature volume change ((V B19 − V B2 )/V B2 ) calculated using the present potential (0.66%) compares well with experimental data (0.69%) [71]. Figure 12(b) shows the atomic configurations of a simulation cell with 3072 atoms before and after the martensitic transition. In contrast to the transition of pure Ti bulk, the transformed structure does not maintain a perfect single crystal but forms a twinned structure with finely dispersed (001) compound twin boundaries. This twinned structure has been frequently observed also in experiment for phase transforming NiTi nanocrystals [3,[89][90][91]. By performing several independent simulations, we confirmed that the distance between each twin boundary is not a fixed value but varies with cell sizes, boundary conditions, and thermal histories. Despite these small variations, a general observation is that the magnitude of the twin distance is in the nanometer regime. Our result agrees well with experimental studies [3,[89][90][91], which show an average twin width of a few nanometers in NiTi nanocrystals. Consequently, we can deduce that the corresponding twin boundary formation energy must be rather low. To quantify this, a DFT calculation and a molecular statics simulation based on the present potential were performed using a bicrystal cell, as shown in Fig. 13. The calculated twin boundary energies by DFT and the present molecular statics simulation are 9.6 and 5.3 mJ/m 2 , respectively. These values are very low as compared, e.g., to twin boundary energies of fcc metals (8 . . . 161 mJ/m 2 ) [92], indicating the special mechanical properties of NiTi.
Compared to the phase transition of pure Ti shown in Fig. 11, the phase transition of the NiTi shape-memory alloy shows a considerably smaller thermal hysteresis (≈1000 K for pure Ti and ≈200 K for NiTi). The reason for this difference is directly related to the occurrence of the nanometer-sized twins in the NiTi alloy, which were absent in pure Ti. As twinning partially relieves the transition strain, a barrier for the phase transition can be significantly reduced. The theoretical window of the thermal hysteresis (≈200 K) is, however, still too large when compared to the experimental window (41 K) [72]. This overestimation is due to the absence of heterogeneous nucleation sites in the MD simulation, as confirmed by our results for the nanopillar to be discussed in the following section (Sec. IV D). The same argument explains also why the theoretical M s temperature of the defect-free NiTi bulk calculated in this section (270 K) is underestimating the experimental one (339 K) [72]. Figure 14(a) shows transition temperatures calculated using various supercell sizes, and Fig. 14(b) shows corresponding snapshots of the low-temperature martensite structures. Similar to the case of pure Ti, we observe that small-sized supercells result in a reduced undercooling, i.e., in a too high M s temperature. For example, for the 384 atom supercell, the M s temperature is about 350 K. Increasing the supercell size, the undercooling increases, and the M s temperature converges to 230 K. For the A f temperature, a more complex convergence behavior is observed. For supercell sizes below about 30 000 atoms, it seemingly converges to a value of about 500 K. However, for larger supercell sizes, a sudden drop occurs, and subsequently the A f temperature shows an oscillatory behavior with an amplitude of 100 K. Moreover, the results of the two independent simulation sets, with relaxed cell angles ("unconstrained") and with constrained cell angles ("constrained"), illustrate that the A f temperatures of constrained supercells are generally lower than those of unconstrained supercells.
The drop and the oscillations in the A f temperature for large-sized supercells can be explained by the occurrence of multidomain martensite structures. As shown in Fig. 14(b), the martensite structures of unconstrained supercells with more than 220 320 atoms and constrained supercells with more than 73 728 atoms exhibit multiple domains that consist of twinned B19 martensite variants with finely dispersed (001) compound twin boundaries. A similar martensite structure with multiple domains was reported in previous experiments, and it was referred to as a "herringbone" structure [3,90]. During reheating, the domain boundaries of the herringbone structure act as heterogeneous nucleation sites for the phase transition, and consequently the A f temperatures of large-sized supercells are lower than the ones of small-sized supercells without domain boundaries. This argument is further supported by performing a separate simulation initiated with defect-free B19 single crystals [red filled triangles in Fig. 14(a)]. The calculated A f temperatures from these simulations show consistently a higher value of 520 K since no heterogeneous nucleation sites are available. For the same reason, the A f temperature of the defect-free B19 single crystals shows an excellent (i.e., non-oscillatory) convergence behavior. The contrastingly poor convergence of the A f temperature of the reheated samples [open and crossed triangles in Fig. 14(a)] can be explained by the different morphology of each herringbone structure resulting from different cell sizes, boundary conditions, and thermal histories. The lower A f temperatures of constrained supercells with respect to unconstrained supercells also can be explained by the difference in the morphology of domain structures. The fixed cell angles during the simulation provide less flexibility to form a single domain. Multiple domains and extended boundary areas are therefore energetically favored, as illustrated in Fig. 14(b).

D. Stress-induced phase transition of a NiTi shape-memory alloy
As a final application of the present paper, the developed potential was applied to a stress-induced phase transition in a NiTi shape-memory alloy. We focused in particular on the transition of a nanosized pillar under a compressive stress. The pillar was prepared as a single crystal consisting of the B2 structure with dimensions of 10.7×10.7×21.1 nm and 175 000 atoms. The sides of the pillar were {110}-type surfaces, and the longitudinal direction was aligned to the [001] direction of the B2 structure. Periodic boundary conditions were applied along the longitudinal direction.
Prior to the loading simulation, transition temperatures of the pillar under a zero-stress state were obtained using the same thermal loading process, as described in the previous section (Sec. IV C). The calculated M s and A f temperatures of the pillar are 320 K and 400 K, respectively, indicating a significantly smaller thermal hysteresis compared to the transition between the defect-free austenite and martensite phases in the previous section (Sec. IV C). This can be interpreted by the fact that the free surfaces provide heterogeneous nucleation sites during the phase transition, which were not available in the defect-free simulations. Incidentally,  [94] and [95]). In each snapshot, sky blue (gray) atoms correspond to the B2 structure and brown (dark gray) atoms to the B19 structure. The intermediate B19 structure cannot be distinguished by the CNA pattern, and respective atoms can fall into both regions, blue or brown. A few of the outmost surface layers (thickness of 0.5 nm) are not visualized for clarity. the theoretical transition temperatures and thermal hysteresis are now very close to experimental values (M s = 339 K, A f = 380 K [72]). Knowing the transition temperatures, the pillar was heated to 450 K, i.e., to a temperature higher than the obtained A f temperature, in order to maintain the B2 austenite phase under a zero-stress state. A stress-controlled uniaxial loading was then applied by adjusting the stress along the longitudinal direction of the pillar. The compressive stress was gradually increased to 2.1 GPa and decreased to 0 GPa with loading and unloading rates of ±3.5 MPa/ps. Figure 15(a) shows a resultant stress-strain response of the pillar, and corresponding snapshots are shown in Fig. 15(b). The discontinuous jumps in the curve represent the occurrence of phase transitions. While the temperature-induced phase transition of the bulk NiTi alloy (Sec. IV C) showed a direct transition between the B2 austenite and the B19 martensite phases, the stress-induced phase transition of the NiTi nanopillar shows a transition via an intermediate phase. During the loading process, the initial B2 phase transforms first into an intermediate B19 phase (B → C), and during further loading, the B19 phase transforms into the B19 phase (D → E). During the unloading process, the B19 phase transforms directly back into the B2 phase (G → H) without the occurrence of an intermediate transition.
Our results clearly indicate that the present potential can be successfully applied for studies of stress-induced as well as temperature-induced phase transitions of NiTi shape-memory alloys. In fact, the full recovery of the initial B2 austenite phase and of the cell dimensions after the loading and unloading processes (A → H) indicates that the nanopillar considered here possesses the important property of superelasticity, i.e., the capability to sustain large elastic strains. The detailed response of nanopillars to an applied stress depends on the size, orientation, and shape of the pillar. The analysis of these conditions is beyond the present scope and is therefore left to future studies.

V. CONCLUSION
An interatomic potential for the binary Ni-Ti system is now available based on the 2NN MEAM formalism. The potential has been developed by improving the unary descriptions of pure Ni and Ti utilizing the force-matching method and by extending it to the binary system. The resulting unary Ti-description can now successfully reproduce the hcp-bcc (α-β) transition that is closely related to the shape-memory transition in the NiTi alloy. The good performance extends to the equiatomic composition where the small structural energy differences and other properties of the various complex structures are faithfully reproduced, resulting in an accurate description of the martensitic shape-memory alloy transformation (B2-B19 ). For achieving the high quality description of the hcp-bcc and B2-B19 transitions, a compromise had to be made by reducing the fitting weight of the hexagonal (ω) phase and by neglecting nonequiatomic compounds in the fitting process. As a consequence, the Ti potential cannot reproduce the hcp-hexagonal (α-ω) phase transition, and the binary potential has deficiencies in describing Ti-rich alloys.
Several applications have been presented for which our potential is particularly well suited. The temperature dependence of the monoclinic angle of the B19 phase has been computed, and the result successfully resolves a previous discrepancy between DFT and experiment. Large-scale MD simulations have been performed to investigate both temperature-and stress-induced phase transitions of the equiatomic alloy related to core applications of the NiTi shape-memory alloy. The MD simulation of the temperature-induced phase transition indicates the occurrence of a nanotwinned martensite structure with multiple domains under varying sizes and constraints in accordance with experiments. The MD simulation of the stress-induced phase transition of the nanopillar indicates a full recovery of the initial structure after the loading and unloading processes correctly reproducing the superelastic behavior of a shape-memory alloy above the critical temperature. The proposed 2NN MEAM potential can be utilized in further studies to accurately describe the structural and mechanical response of shape-memory alloys on the nanometer length scale. The strategies outlined here to construct realistic empirical potential for systems that show a complex phase diagram with many competing phases and defect structures can be straightforwardly applied to other material systems.