Development of a second-nearest-neighbor modified embedded atom method potential for silicon–phosphorus binary system

Computational Materials Science 120 (2016) 1–12

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Development of a second-nearest-neighbor modified embedded atom method potential for silicon–phosphorus binary system Bin Liu a,b, Hao Zhang b, Junyong Tao a,⇑, Ziran Liu b,c, Xun Chen a, Yun’an Zhang a a

Science and Technology on Integrated Logistics Support Laboratory, College of Mechanics and Automation, National University of Defense Technology, Changsha 410073, China Department of Chemical and Materials Engineering, University of Alberta, Edmonton, AB T6G 2V4, Canada c College of Physics and Information Science, Hunan Normal University, Changsha 410073, China b

a r t i c l e

i n f o

Article history: Received 22 December 2015 Received in revised form 30 March 2016 Accepted 2 April 2016

Keywords: Second nearest-neighbor modified embedded atom method Genetic algorithm Silicon Phosphorus Mechanical properties

a b s t r a c t Phosphorus (P) is one of the most common impurities in silicon (Si). To investigate the effects of P on the mechanical properties of Si at nano or atomic scale, the second-nearest-neighbor Modified Embedded Atom Method model (2NN MEAM) potentials for pure P and Si–P binary system are developed using genetic algorithm (GA). The reference physical properties for the parameterization for pure P include cohesive energies, lattice constants, bulk moduli and first order pressure derivatives of experimentallyexisting phases. The van der Waals interactions in P are not considered since it’s beyond the description of MEAM and it’s less important in covalently bonded Si–P compounds. For Si–P binary system, cohesive or formation energies, lattice constants and elastic constants of typical Si–P structures, properties of point defects and P–vacancy pairs are considered. The robustness of the newly developed potentials is examined by testing the thermal stability and transformation from non-equilibrium to ordered phase. The effects of P on the tensile behaviors of bulk Si are demonstrated. Results show that phosphorus lowers the fracture threshold and stiffness of bulk Si. Ó 2016 Elsevier B.V. All rights reserved.

1. Introduction Si is a technologically important material in integrated circuit (IC), micro-/nano- electro-mechanical systems (MEMS/NEMS), solar cell, etc. As a representative covalent element, it has great significance in theoretical research and has drawn much attention during the past decades. In engineering, P is one of the most common dopants in Si and plays a critical role in its electrical properties. In addition, P significantly influences the mechanical [1] and thermal [2] properties of Si, about which many experimental reports [1,3–5] can be referred to. In order to reveal the underlying mechanisms behind the experimental phenomena, it’s essential to raise an atomistic-scale technique. Several techniques enable investigation of atomic interactions, among which quantum mechanics-based density functional theory (DFT) and empirical method molecular dynamics (MD) are most widely used. DFT provides the most reliable information about materials at atomistic-scale, which, however, cannot deal with a system with thousands of atoms despite of rapid development of computer technologies. Somehow, many physical phenomena, such as dislocation emission and expansion, can be only observed ⇑ Corresponding author. E-mail address: [email protected] (J. Tao). http://dx.doi.org/10.1016/j.commatsci.2016.04.002 0927-0256/Ó 2016 Elsevier B.V. All rights reserved.

in large systems. Therefore, empirical potentials are applied extensively in simulation experiments. Only a few empirical potentials for Si–P binary system can be referenced in the literatures [2,6]. Lee and Hwang [2] developed a SW potential for Si–P binaries to better describe the effects of P impurities on the thermal conductivity of Si. Wilson et al. [6] investigated the diffusion of P on the surface of Si crystal using Environment-Dependent Interatomic Potential (EDIP) model. In addition, some other P-based multi-element potentials can be referred to. For example, Ackland et al. [7] developed a Fe–P potential based on DFT calculations. This model, as the authors pointed out, is unable to describe covalent materials, thus is only suitable for low P concentration systems. Above-mentioned empirical models are applicable for only a few elements and have unsatisfactory extendibility for the other elements. Moreover, they are incapable of maintaining the brittle nature of Si. Although Hauch et al. [8] observed brittleness using the modified SW potential, some other important properties such as elastic constants and melting point become unacceptable. Mattila et al. [9,10] and Sheng et al. [11] developed a potential for Ni–P binary system using Embedded Atom Method (EAM) model which has been applied successfully for many elements and multi-element systems. However, it cannot perform well in describing amorphous Si. Mattila’s and Sheng’s work have enlightened us that EAM model can be used for

2

B. Liu et al. / Computational Materials Science 120 (2016) 1–12

P-based system, then the MEAM model developed from EAM can accordingly be applied in our work. A 2NN MEAM potential for pure P has been developed by Ko et al. [12]. Considering the critical role of reference structure in the thermal stability of Si crystal with high P concentration, we have to develop a new parameter set for P according to the specific requirements in this work. MEAM potential can maintain the brittleness of Si [13–15], and give good predictions for thermal [16], vacancy [17], phase stability [18], stacking faults [19] and amorphous [18] properties. It is thus clear that MEAM is a solid choice for our work. MEAM was constructed by Baskes et al. [20] by adding angular terms to EAM formalism [21], and has been extensively applied in a wide range of fields. Currently, 26 elements have been covered by MEAM [22], including metal materials, semiconductors, gas elements like nitrogen (N), oxygen (O), hydrogen (H), which shows the universality of MEAM formalism. Lee and Baskes developed a more generalized version, second-nearest-neighbor modified embedded atom method model (2NN MEAM) [23], which has been widely used for binary and ternary systems [24]. Up to now, Sibased multi-element MEAM potential include Si–Ni [25], Si–Li [26], Si–Al/Fe/Cu/Mg [27]. To lay the foundation of investigating the effects of P impurities on the mechanical properties of Si in our future work, this article reports the development of 2NN MEAM potentials for pure P and Si–P binary system using genetic algorithm (GA). The reference physical properties for parameterization are all from experiments or higher level calculations. The van der Waals interactions (existing in white (a) P, black P (A17) and a transformed phase (A7)) are omitted by only taking single molecule white (a) P (written P4), single layer A17 and A7 into account. Then the potential is extended to the Si–P binary system to describe various structural arrangements including typical diamond-like compounds, phases reported in phase diagram, and dopant P related defects. The robustness of the newly developed potentials is examined by testing the thermal stability and transformation of Si crystal with high P concentration from non-equilibrium to ordered arrangement. Last, the effects of P on the tensile behaviors of Si are tested to give explanations to experimental phenomena. 2. Methodology 2.1. 2NN MEAM interatomic potential The complete 2NN MEAM theory can be found in Ref. [24]. Here we introduce its primary framework and illustrate the differences from other related literatures. In the 2NN MEAM model, the total energy of a system is given as:

" # X 1X  E¼ F i ðqi Þ þ Sij /ðRij Þ 2 j–i i

ð1Þ

 i Þ is the embedding function, q  i is the background elecwhere F i ðq tron density at site i. Sij and /ðRij Þ are respectively the screening function and the pair interaction between atom i and j separated  consists of four by a distance Rij . The background electron density q terms (one spherically symmetric partial term qi

ð0Þ

and three angu-

q and q associated by the adjustable lar contributions q parameters t(h) (h = 1–3).  is calculated as follows: Background electron density q ð1Þ i ,

ð2Þ i

ð3Þ i )

q i ¼ qð0Þ i GðCÞ

ð2Þ

Gamma function GðCÞ can be expressed as:

GðCÞ ¼

2 1 þ expðCÞ

ð3:1Þ

( pffiffiffiffiffiffiffiffiffiffiffiffiffi C P 1 1 þ C; pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi GðCÞ ¼  j1 þ Cj; C < 1

ð3:2Þ

where C combines qi , qi and t(1–3), see Ref. [24]. (3.1) is used in this work, while some other literatures applied (3.2) [17,22,27]. The atomic electron density is given by: ð0Þ

ð1—3Þ

qaðhÞ ðRÞ ¼ q0 exp½bðhÞ ðR=re  1Þ (h)

ð4Þ

where b (h = 0–3) are adjustable parameters, re is the first neighbor distance of the reference structure. The scaling factor q0 has no effect on pure element systems while greatly affect multi-element systems. In MEAM formalism, pair interaction /ðRÞ doesn’t have any formula. Instead, it is evaluated from the values of the total energy per atom in the reference structure, where individual atoms are on the exact lattice point, and the embedding energy, as a function of the nearest-neighbor distance. The total energy per atom is estimated from the zero-temperature universal equation of state of Rose et al. [28]. Reference structure is crucial for a MEAM potential since it directly influence its reliability [29]. The reference structure applied in this work is diamond-structured, i.e. interlaced FCC for  alloy, with the Space Group F43m, see Fig. 1(a), which has been reported to experimentally exist [30]. We have tested B1-NaCl structure, the most commonly used reference structure [24,27,29], which can well reproduce the reference physical properties. However, the thermal stability of Si7P (Fig. 1(b)) is not acceptable. B2, l12 and C11 are also considered and passed for  the same reason. F43m gives overall better predictions on various  properties, especially lattice constant of F43m, elastic constants of Si7P, thermal stability and transformation from non-equilibrium to  ordered arrangement of Si7P. For F43m -SiP structure, the total energy per atom (1/2Si + 1/2P) is given as:

Phosphorus

Silicon

(a)  Fig. 1. Schematic view of (a) F43M and (b) Si7P.

(b)

B. Liu et al. / Computational Materials Science 120 (2016) 1–12

EuSiP ðRÞ ¼

1 1 Z Z  Si Þ þ F P ðq  P Þ þ 1 /SiP ðRÞ þ 2 ½SSi /SiSi ðaRÞ F Si ðq 2 2 2 4 þ SP /PP ðaRÞ

2.2. Application of GA in the parameters optimization of 2NN MEAM for Si–P binary system

ð5Þ

where Z1 and Z2 are respectively the number of first- and second nearest atoms in F43m structure, SSi and SP are the screening functions for the second nearest neighbor interactions between Si atoms and between P atoms respectively, and a is the ratio between the second- and first-nearest distances in the reference structure. The pair interaction between Si and P can now be obtained:

/SiP ðRÞ ¼

  1 Z  Si Þ  F P ðq  P Þ  2 ½SSi /SiSi ðaRÞ þ SP /PP ðaRÞ 2EuSiP ðRÞ  F Si ðq 2 Z1 ð6Þ

where the embedding function F Si and F P , the pair interaction /SiSi and /PP are calculated by single element 2NN MEAM potential for Si and P. The universal equation of state EuSiP can be obtained by Rose equation [28]:

EuSiP ðRÞ ¼ Ec ð1 þ a þ d  ða Þ3 Þea a ¼ aðR=re  1Þ a ¼ ð9BX=Ec Þ1=2 



d¼

ð7Þ ð8Þ ð9Þ



drepuls; a < 0 dattrac;



a > 0

3

ð10Þ

where B, X and Ec are the bulk modulus, equilibrium atomic volume and cohesive energy of the reference structure, respectively. drepuls and dattrac are adjustable parameters. To our knowledge, most of the publications set d = drepuls = dattrac. However, we found that drepuls and dattrac have conflicting effects on pair interactions. So they are optimized independently. In a 2NN MEAM potential for a binary system, the following parameters need to be determined: Ec, re, B, drepuls and dattrac in universal state of equation; scaling factor q0 in atomic electron density; two model parameters Cmin and Cmax that determine the extent of screening of an atom (k) from the interaction between two neighboring atoms (i and j). In current work, Ec, re and B for  F43m are calculated by DFT. However, in order to get the best overall agreement of properties produced by MEAM with reference values, they are tuned within a small range. For Cmin and Cmax in a Si–P binary system, there are four cases, namely Si–P–Si, Si–Si–P, P–Si–P and P–P–Si. The formula of Sijk expressed by Cmin and Cmax can be found in Ref. [24]. Besides the above mentioned adjustable parameters, rc (cut-off distance) and Dr (for designating the cut-off region) should be cautiously determined. Kang’s work [14] has shown that the brittleness of Si is closely related to rc. In addition, our recent study on ideal tensile strength of Si showed that when b(2) in the parameter set for Si is below 2.5, a sudden leap of the excess energy increment can be observed, leading to the discontinuity of fracture [19]. This artificially introduced stress barrier can be moved out of the typical opening distance (distance between two fractured surface) by increasing rc according to Ko et al. [31] (7.0 Å in Ref. [31]). However, increment in rc brings much more computational burden. Thus 6.0 Å (4.5 Å in Ref. [18]) is accordingly chosen for this work. According Ko’s recent research [31], enlarging Dr can improve the description of cleavage fracture as the artificial stress barrier is effectively relieved. However, larger Dr would lead to worse prediction for elastic constants, vacancy and surface properties of Si, while this effect on Si–P binary system is much less marked. Therefore, as a compromise, Dr is set 3.0 Å (4.5 Å in Ref. [31]).

Currently, many MEAM potentials for Si have been developed [16,17,22,26,27,32,33] according to different requirement of simulation. We have optimized the 2NN MEAM parameter set for Si to improve predictions for surface and stacking fault properties in our recent work [19]. However, it cannot yield as good amorphous properties as Lee’s potential [18]. Considering the mechanism that P impurities affect the mechanical behavior of amorphous Si will be investigated in our future work, Lee’s parameter set for Si is used in this work. Most of the adjustable parameters have effects on several reference physical properties, and each property is usually affected by more than one parameter. Thus tuning the parameters one by one manually is time consuming. Therefore simultaneous adjustment is worth trying, which can be realized by GA that a population of candidate solutions (called individuals, creatures, or phenotypes) to an optimization problem is evolved toward better solutions. GA would be efficient in global searching when given suitable fitness function. Another global optimization technique, particle swarm optimization (PSO), has been applied in a similar way [26,34]. The optimization or evolution starts from randomly generated individuals, namely multiple 2NN MEAM parameter sets in present work. Then the physical properties corresponding to each parameter set are calculated. The fitness function is used to measure the quality of each parameter set and determine its successive probability. Selection operation decides fitting parameter sets that can be used to breed new generation based on fitness-estimation. After selecting parental candidates, the next parameter sets are generated by crossover and mutation operations. Repeat these steps until the parameter set produces good physical properties compared with reference values. The whole technical process is completed by Matlab and LAMMPS [35]. The former executes GA algorithm, and the latter carries out Molecular Dynamics (MD) calculations and feeds the results back to Matlab. Assume that vector x ¼ ½x1 x2 x3    x13  consists of candidate parameters. The optimization procedure is supposed to minimize the overall prediction error that can be regarded as the objective function defined as:

f ðxÞ ¼

N X wi ½pi ðxÞ=Pi  12

ð11Þ

i¼1

N is the number of reference physical properties in the data base, w is normalized weight factor, p(x) and P are respectively a physical property calculated with parameter set x and the corresponding reference value. For Pure P, the weights are evenly allocated. For Si–P, smaller weights are allocated to elastic constants and structural  properties of F43M, Si7P and Si63P. Actually, the agreement between MEAM predicted elastic constants (C12 and C44) and reference values is unsatisfactory, even given much larger weight.  Besides, the structural properties of F43M, Si7P and Si63P can be well predicted with smaller weight. The 2NN MEAM potential parameters for pure Si and P are listed in Table 1. Table 2 lists the parameters for Si–P binary system. It’s worth noting that the reference structure for P is diamond in this work. We first examined body-centered cubic (BCC), an experimentally-existing phase under high pressure, however, when the P concentration reaches up to 12.5%, the crystal structure (Si7P) is unstable under 300 K and 1 bar. Since the diamond is a hypothetical structure for P, its cohesive energy Ec, bulk modulus B and nearest neighbor distance re are tuned independently. dattrac and drepuls have higher precision of 103 since a slight change of

4

B. Liu et al. / Computational Materials Science 120 (2016) 1–12

Table 1 2NN MEAM parameters for Si and P. rc, re and Dr in Å, Ec in eV, B in GPa.

Si P

rc

Dr

Ec

B

re

b(0)

b(1)

b(2)

b(3)

A

t(1)

t(2)

t(3)

Cmin

Cmax

q0

dattrac

drepul

6.0 6.0

3.0 3.0

4.63 2.79

99 44.3

2.35 2.36

3.55 2.07

2.50 4.93

0.00 0.57

7.50 0.39

0.58 1.23

1.80 0.56

5.25 0.01

2.61 2.06

1.41 0.52

2.80 4.36

1.00 0.89

0 0.094

0 0.196

Table 2 2NN MEAM parameters for Si–P binary system. rc, re and Dr in Å, Ec in eV, B in GPa. Ec

rc

Dr

B

re

q0 (q0P/q0Si)

Cmin121

Cmin112

Cmin122

Cmin212

Cmax121

Cmax112

Cmax122

Cmax212

dattrac

drepul

3.94

6.0

3.0

62.68

2.31

0.89

0.40

0.76

1.00

1.00

4.99

2.77

2.80

2.80

0

0.228

necessary to consider the van der Waals interactions to generally well predict the properties of pure P using an empirical formalism. However, when dissolved in Si, the van der Waals forces of P must be less important compared with covalent bonds. In this sense, the properties concerned with van der Waals forces are not taken into account in the parameterization process. In detail, only the single P4 molecule (white P), single layer A17 and A7 structures (black P), and transformed phases (SC, Imma, SH and BCC) are considered in the GA optimization, while the amorphous-like red P is not considered. This processing method is the same as Ko’s work [12]. A similar approach was applied in the development of a 2NN MEAM potential for pure carbon [46].

these two parameters results in big change of MD predictions, while the precision of the other parameters is 102. First-principle calculations in this work are based on density functional theory (DFT), using the projector augmented-wave (PAW) method [36] as implemented in the Vienna Ab-initio Simulation Package (VASP) [37]. Exchange–correlation effects are approximated by the generalized gradient approximation (GGA) as parameterized by Perdew et al. [38]. Integration over the irreducible Brillouinzone was performed using the Monkhorst–Pack scheme. In the cases of bulk A17 and A7, van der Waals corrections [39] are applied. Before DFT calculations, the rationality of pseudopotential, convergence condition and Brillouin zone sampling are verified. For Si, DFT predictions of elastic constants, cohesive energy and lattice constants are compared with experimental values, given in parentheses: C11 = 162.69 (168[40]), C12 = 65.89 (65 [40]), C44 = 76.07 (80[40]), B = 98.16 (99[40]), Ecohesive = 4.55 (4.63 [18]), a = 5.46 (5.431[18]). For P, DFT calculations of cohesive energy of A17 structure, and atomic volumes of various transformed phases are compared with experimental values: Ecohesive_A17 = 3.48 (3.43[41]), Vatom_A17 = 19.55 (19.0[42]), Vatom_A7 = 15.71 (16.5[42]), Vatom_SC = 14.53 (15.2[42], SC stands for simple cubic), Vatom_SH = 14.45 (14.2[42], SH stands for simple hexagonal). The good agreement gives us confidence for reliable DFT calculations.

3.1. Cohesive energy and lattice constants DFT and MEAM calculations of cohesive energies and lattice constants are listed in Table 3. One can see that the results from MEAM are in good agreement with DFT. It’s worth noting that DFT calculations show that the cohesive energy of A17 is the highest, proving that A17 is the most stable phase from the perspective of energy [41]. It should be mentioned that the DFT value of lattice constant for BCC is 2.36 Å in Ref. [47], smaller than that under 280 GPa, 2.42 Å, which is incorrect. It is worthwhile to note that the bulk phase structures of A17 and A7 predicted by MEAM potential are in accordance with experimental configurations. Hence the bulk structure information of A17 and A7 is given in Table 3. As for the cohesive energies of A17 and A7, MEAM can’t yield correct results as the distances between layers are smaller than cutoff distance. It means that neighbor layers are within the range of covalent interactions which should be van der Waals interactions in physical realities. Therefore, MEAM can only deal with the energy properties of single layer A17 and A7 (see Fig. 3). Actually, this process is reasonable because the DFT calculation of cohesive energies of single layer A17 and A7 are not much different from that of bulk A17 and A7, as shown in Tables 3 and 4. This indicates that van der Waals interactions have marginal effects on cohesive energies, meaning that the cohesive energies of single layer A17 and A7 can represent those of their bulk counterparts. The MEAM predictions of energies and structural information for

3. Calculation of physical properties for P A17, A7, SC, SH and BCC (see Fig. 2), with space group Cmca,  R3m, Pm3m, P6/mmm and Im-3M, respectively, are experimentally observed crystal phases. A17-P is known as black P, while the others are transformed phases. Besides, an experimentally unidentified intermediate phase between SC and SH has been confirmed to be Imma by Rajeev et al. [42] using first principle calculations. One difficulty in using MEAM to describe P is that though MEAM model can deal with covalent and metallic bonds, van der Waals interactions are beyond its description. Van der Waals forces exist in white P (a-P, bonded cubic aggregation of condensed P4 molecules) [43], red P (the most complex amorphous-like structure) [44], and layered A17 and A7 phases [45]. It is certainly

a

b

b

c

c

(a)

a

(b)

a (c)

a

60o

c (d)

a

(e)

Fig. 2. Schematic view of (a) bulk black P A17, and transformed phases (b) bulk A7, (c) SC, (d) Imma, (e) SH and (f) BCC.

a

(f)

5

B. Liu et al. / Computational Materials Science 120 (2016) 1–12

Table 3 MEAM predictions for cohesive energy, in eV, and lattice constants, in Å or °, for various phases of phosphorus, in comparison with experimental and DFT values. The meanings of the structural symbols are illustrated in Fig. 2. Here, the A17 and A7 are bulk structures whose configurations can be well predicted while the cohesive energies cannot because of the incapability of MEAM model in describing van der Waals interaction. Structure

Properties

DFT

Expt.

MEAM

A17

Cohesive energy Lattice constants

3.48a, 3.413d 3.34a 10.62a 4.41a 3.44a 3.62a 57.25a 3.38a, 3.335d 2.44a, 2.50d, 2.40d 3.29a 5.04a, 4.90i 3.71a, 4.41i 3.43a, 2.65i 3.04a, 2.949d 2.61a, 2.58i 2.45a, 2.45i 2.80a, 2.693d 3.05a, 3.08d, 3.01d

3.43b, 3.44c 3.31e 10.47e 4.37e – – – – 2.48g – – – – – – – – –

– 3.46f 10.65f 3.46f – 3.63f 57f 3.34f, 3.356h 2.50f, 2.48h 3.22f 5.10f 4.46f 2.58f 3.18f, 3.115h 2.58f 2.57f 2.79f, 2.720h 3.01f, 3.00h

A7

SC Imma

SH

BCC a b c d e f g h i

a b c

Cohesive energy Lattice constants

a h

Cohesive energy Lattice constants Cohesive energy Lattice constants

a a b c

Cohesive energy Lattice constants

a c

Cohesive energy Lattice constants

a

Present DFT calculation. Experimental data from Ref. [41]. Experimental data from Ref. [48]. First principle calculation from Ref. [12] and references therein, the sizes are converted from volumes. Experimental data from Ref. [49]. Present MEAM calculation. Experimental data from Ref. [42]. MEAM calculation from Ref. [12], the sizes are converted from volumes. First principle calculation from Ref. [42].

(a)

β

l2

(c)

(b)

α

α

α

l1

l

l

Fig. 3. Schematic view of (a) single layer A17, (b) single layer A7 and (c) single molecule P4.

Table 4 MEAM predictions of structural information for single layer A17 and A7 and single molecule P4, in comparison with experiments and DFT calculations, lengths in Å and angles in °. Note: the experimental information is extracted from bulk structures while the DFT and MEAM values are from single layer or molecule models. Structure Single layer A17

Properties Cohesive energy Structural information

Expt. l1 l2

a b

Single layer A7

Cohesive energy Structural information

l

a Single molecule P4

Cohesive energy Structural information

l

a a b c d e f

Experimental data from Ref. [49]. Present DFT calculation. Present MEAM calculation. MEAM calculation from Ref. [12]. Experimental data from Ref. [50]. Experimental data from Ref. [51].

DFT

MEAM

– 2.22a 2.28a 96a 102a

b

3.47 2.22b 2.26b 96b 104b

3.39c 2.45c, 2.44d 2.31c, 2.46d 90c, 90d 90c, 90d

– 2.13e 105e

3.39b 2.28b 95b

3.40c 2.22c, 2.32d 74c, 112d

– 2.21f 60f

3.34b 2.20b 60b

3.28c 2.19c, 2.43d 60c, 60d

6

B. Liu et al. / Computational Materials Science 120 (2016) 1–12

single layer A17 and A7, and single molecule P4 are listed in Table 4. Although the agreement of crystalline structure (bond lengths) of single layer A17 and A7 is not good, the cohesive energies and properties of single molecule P4 are well reproduced. While the cohesive energy of single layer A7 is slightly larger than that of single layer A17, this work takes A17 as reference state for defects-related calculations of Si–P binary systems as A17 is experimentally the most stable phase.

experiments. It’s worth noting that MEAM predictions for A17 are not given as its thermal stability is not good. Imma, an intermediate phase between SC and SH as reported by Ahuja et al. [42], exists in the range between 107 GPa and 136 GPa. However, when described by present MEAM potential, it remains Imma configuration only below 118 GPa. Hence, MEAM values for Imma are given from 107 GPa to 118 GPa. Likewise, MEAM predictions for SH are given from 137 GPa to 149 GPa.

3.2. Bulk modulus band first order pressure derivative B0

4. Calculation of physical properties for Si–P binary systems

B and B0 are fitted using Birch equation [52] defined as:

4.1. Cohesive energy and lattice constants of typical Si–P structures

where P is the pressure, V0 is the initial volume, and V is the volume under P. Thermal stability of the existing phases is an important criterion for the robustness of a potential [53]. All the phases are thermally stable except A17 whose B and B0 cannot be fitted using Eq. (12). A17, SC, Imma, SH and BCC are modeled as orthorhombic structures in this work. Orthorhombic SH is constructed using the methods in Ref. [54]. The bulk moduli of these five phases can also be calculated as (C11 + C22 + C33 + 2 ⁄ C12 + 2 ⁄ C13 + 2 ⁄ C23)/9 [55], where C11, C22, C33, C12, C13 and C23 are elastic constants. The bulk moduli of SC and BCC calculated by above mentioned two methods differ only a little, of which the results from the first method are listed in Table 5. The results are shown in Table 5. The MEAM predicted bulk moduli and pressure derivatives are overall in good agreement with reference values except A17 and A7 phases which differ by a factor of about 2, probably due to the absence of van der Waals interactions. It should be mentioned that A17 is not thermally stable so that its pressure derivative cannot be calculated using present MEAM function. 3.3. Pressure dependence of atomic volume of each phase The pressure dependence of atomic volume of each phase of black P is plotted in Fig. 4. It shows good agreement with

16 15

A7-Akahama A7-MEAM SC-Akahama SC-MEAM Imma-Ahuja Imma-MEAM SH-Akahama SH-MEAM BCC-Akahama BCC-MEAM

14 3

ð12Þ

The concentration of P, as a guest element dissolved in crystal Si, is usually not high. The solubility limit of P in Czochralski (CZ) Si is 1.3  1021 cm3 [3], about 2.6% compared with the number density of single crystalline silicon (SCS) 4.90  1022/cm3 [59]. It

Atomic volume (A)

"  "  #)  5=3 #( 7=3 2=3 3 V0 V0 3 V0 1 þ ðB0  4Þ P¼ B   1Þ 2 4 V V V

R

13 12 11 10 9 8 118GPa

7 6

149GPa

0

50

100

150

200

250

300

Pressure (GPa) Fig. 4. Pressure dependence of atomic volume of different phases of black P, in comparison with experiments (Akahama [47]) and DFT calculations (Ahuja [42]). MEAM predictions for A17 are not given because of its unacceptable thermal stability. Though Imma and SH exist within a wide range of pressures, present MEAM potential can only predict results up to 118 GPa and 149 GPa, respectively.

Table 5 Bulk modulus B, in GPa, and first-order pressure derivative B0 of various phases of phosphorus, in comparison with DFT and experimental values. Structure

Properties

Expt.

DFT

MEAM

A17

Bulk modulus B Pressure derivative B0

36.0 ± 2.0a 4.5 ± 0.5a

32b 4.5b

62.5c –

A7

Bulk modulus B Pressure derivative B’

46.0 ± 4.0a 3.0 ± 0.6a

62b 5.3 b

93.1c 3.4c

SC

Bulk modulus B Pressure derivative B0

70.7 ± 0.9a, 95 ± 5d 4.69 ± 0.1a, 2.1 ± 0.8d

132.0b, 127.0e, 104g, 113h 3.3b, 4.39e

110.5c, 88.2f 4.4c

Imma

Bulk modulus B Pressure derivative B’ Bulk modulus B Pressure derivative B’

– – 84.1 ± 1.0a 4.69a

90.0b 3.6b 99.0b, 87.5h 3.5b

93.9c 3.9c 93.8c, 99.7f 4.6c

Bulk modulus B Pressure derivative B’

– –

108.2i, 111.0b, 98.8g, 100h 3.84i, 3.0b

116.2c, 108.2f 4.3c

SH BCC a b c d e f g h i

Experimental value from Ref. [47]. DFT calculation from Ref. [42]. Present MEAM calculation. Experimental value from Ref. [56]. DFT calculation from Ref. [57]. MEAM calculation from Ref. [42]. DFT calculation from Ref. [12]. DFT calculation from Ref. [58]. DFT calculation from Ref. [47].

7

B. Liu et al. / Computational Materials Science 120 (2016) 1–12

is also growth orientation dependant. The maximum P concentration in (1 1 1)-grown CZ Si is reported to be 1.11  1020/cm3 [3], only 0.23%. However, according to Weeks’s recent experiments [60], the concentration of substitutional P in epitaxial Si can reach up to 10%. In view of the local high concentration which practically exists in P-doped Si structures, a Si–P structure with the concentration of 12.5% (one Si atom replaced by one P atom in the basic cell) is considered in this work, see Fig. 1(b). Cohesive energy is defined as the needed energy for an atom to transform from solid to free state, calculated as:

Ec ¼

NSi ESi þ N P EP  Ecrystal NSi þ NP

ð13Þ

where ESi and EP are respectively the free energy of single Si and P atom, NSi and NP are respectively the atom number of Si and P in the system, Ecrystal is the total energy of the crystal system. The heat of formation per atom is defined as:

Hf ¼

Ecrystal  NSi eSi  NP eP NSi þ NP

ð14Þ

where eSi and eP are the total energies per atom for Si and P in their ideal bulk structures. For Si, the ideal bulk structure is diamond, and that for P is white (a) P. Table 6 shows the formation heats, cohesive energies and lattice  constants of diamond-like structures F43m, Si7P and Si63P (a 2  2  2 unit cell with a substitutional P atom), and three relatively widely accepted to experimentally-existing phases, including CMC21, PBAM and PA3 (see Fig. 5) [63–65]. The reported data about the solid solubility of P in Si are rather conflicting as the sample fabrication is extremely restricted by temperature, pressure, proportion of P, etc. [66]. Probably the most convincing one is Kooi’s work [67], which showed that a stable phase of Si–P

solid was written as SiP. According to Wadsten et al. [63], the space  group of SiP is CMC21, while Osugi et al. [30] reported F43m. However, the former seems more acceptable [62]. This is supported by our DFT calculation which yields negative and positive formation  heat for CMC21 and F43m, respectively.Si2P [68] and SiP2 [65] have also been reported. However, the structure of Si2P is not verified till now, thus is not considered in this study. The SiP2 phase is reported to exist in two configurations, namely orthorhombic (PBAM) [64] and pyrite-type (PA3) [65]. Among the above Si phosphides, CMC21 is believed to be the most stable and the others are metastable [66]. The experimental and CALPHAD values for formation heat in Table 6 are based on diamond Si and white (a) P as reference states. As mentioned above, white (a) P cannot be used as reference state for MEAM calculation because of the absence of van der Waals interactions. Though the single molecule P4 of white P has been considered in MEAM parameter set optimization, the MEAM prediction of cohesive energy for bulk white (a) P is not acceptable. MEAM can only give formation heats based on black P as reference state. To make direct comparisons, the experimental or DFT values of the transformation heat between black and white P is necessary. A similar approximation was applied in Ko’s research [12]. The agreement of structural sizes, cohesive energies, especially formation heats of CMC21, PBAM and PA3 is not perfect. Actually, we fail to reproduce good results for both diamond-structured Si–P crystals and high P Si–P compounds, namely CMC21, PBAM and PA3 at the same time. Therefore, only the formation energy of CMC21 is involved in the process of parameter optimization since it is a stable phase in the phase diagram while the other two are metastable. Three points should be mentioned. First, the purposes of taking Si63P into account in this work are to investigate the

Table 6  Formation heats (in kJ/mol), cohesive energies (in eV), and lattice constants (in Å) of diamond-like crystals F43M, Si7P and Si63P, and high P phases CMC21, PBAM and PA3, in comparison with experiments and other calculations. The reference states are diamond Si and white (a) P. Compound

Properties

Expt.

CALPHADb

DFTa

MEAMc

 F43M

Formation heat, Ef Cohesive energy, Ec Lattice constants

a

– – 5.24d

– – –

10.11 3.87 5.33

6.72 3.94 5.33

Si7P

Formation heat, Ef Cohesive energy, Ec Lattice constants

a

– – 5.38e

– – –

0.17 4.38 5.42

0.96 4.45 5.42

Si63P

Formation heat, Ef Cohesive energy, Ec Lattice constants

a

– – 5.422f

– – –

1.58 4.54 5.459

0.08 4.61 5.429

CMC21

Formation heat, Ef Cohesive energy, Ec Lattice constants

30.95b

26.71 4.19 4.73 3.54 20.52

43.21 4.35 4.88 2.84 18.53

PBAM

PA3

a b c d e f g h i

Formation heat, Ef Cohesive energy, Ec Lattice constants

Formation heat, Ef Cohesive energy, Ec Lattice constants

a b c

4.67g 3.51g 20.49g

31.65 – – – –

a b c

– – 13.97h 3.44h 10.08h

26.43, 26.69 – – – –

 – – – –

38.40 3.26 13.86 2.95 10.07

a

– – 5.70i

– – –

25.67 3.97 5.74

55.68 3.23 5.63

Present DFT calculation. The heat of transformation between black P and white (a) P is 21 (2 kJ/mol[61]. Experimental data and CALPHAD calculations from Ref. [62]. Present MEAM calculation. Experimental data from Ref. [30], measured under 4 GPa. Measured using the sample with the concentration of 12.5% in Ref. [60] Measured using the sample with the concentration of 1.56% in Ref. [60]. Experimental data from Ref. [63]. Experimental data from Ref. [64]. Experimental data from Ref. [65].

8

B. Liu et al. / Computational Materials Science 120 (2016) 1–12

a

c

(c)

(b)

(a)

c

b

b

a

a

Fig. 5. Schematic view of (a) CMC21 (single chain), (b) PBAM and (c) PA3. The meanings of colors are in accordance with Fig. 1. Single chain of CMC21 is considered in this work to avoid van der Waals interactions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

changing trend of cohesive energy and lattice constant with P concentration, and to optimize the structure for calculating its elastic constants in the following part. In addition, the lattice constant of Si63P is divided by 2 for better comparing with other structures. Second, our DFT calculation cannot predict PBAM to be a stable structure. Third, only the single chain of CMC21 structure is considered for the chains interact with each other through van der Waals forces, so a for CMC21 in Table 6 means the width of the chain.

4.3. Point defects The P dopant mainly exists substitutionally with a small proportion in SCS. Another important form is interstitial. Hexagonal, Tetrahedral and Split-h1 1 0i interstitial (Fig. 6) are considered in present work. The formation energy of a substitutional point defect Esub is defined by:

Esub ¼ Etot ½ðN  1Þ  Si þ PSi   fEtot ½N  Si  eSi g  eP

4.2. Elastic constants The DFT calculation for elastic constants and bulk modulus of  F43M, Si7P and Si63Pis done with elastic energy theory [69]. The results are listed in Table 7. Both DFT and MEAM values indicate that the elastic constants and bulk modulus would change with the concentration of P. DFT calculations show that C11, C44 and B decrease as the concentration increases, while C12 increases, which is consistent with experimental phenomenon [1,70] and theoretical calculations [70–72]. However, MEAM values for C12 are on the contrary.

ð15Þ

where Etot ½ðN  1Þ  Si þ P Si  is the total energy of the Si system with N  1 Si atoms a neutral substitutional P atom (P Si ; Kröger–Vink notation [73]), Etot ½N  Si is the total energy of perfect Si system, eSi and eP are respectively the total energy per atom for Si and P in their ideal bulk structures. The formation energy of an interstitial point defect Eint is defined as:

Eint ¼ Etot ½N  Si þ Pi   Etot ½N  Si  eP

ð16Þ

Table 7  Elastic constants and bulk modulus of F43M, Si7P and Si63P, in GPa. Crystal structure

 F43M Si7P Si63P

C11

C12

C44

B

DFT

MEAM

DFT

MEAM

DFT

MEAM

DFT

MEAM

70.67 128.80 147.39

79.23 136.26 156.39

87.79 76.17 67.84

54.62 62.43 66.09

26.51 39.63 63.11

15.37 53.05 75.32

82.09 93.72 94.35

62.82 87.04 96.18

(a) Hexagonal

(b) Tetrahedral

(c) Split-<110>

Fig. 6. Schematic view of (a) Hexagonal, (b) Tetrahedral and (c) Split-h1 1 0i interstitial.

9

B. Liu et al. / Computational Materials Science 120 (2016) 1–12 Table 8 Formation energies (in eV) of substitutional and interstitial defects, in comparison with experimental and DFT values. The reference states are diamond-Si and A17-P. Point defect

64 (65) atoms Substitution Interstitial

a b c d

MEAMd

DFT

Hexagonal Tetrahedral Split-h1 1 0i

0.23a 3.23a, 3.5b 4.34a 3.40a

Table 9 MEAM predictions of formation energy for P–V, activation energy of migration of P, and binding energy of P and vacancy, all in eV, in comparison with experimental and DFT values. DFT

216 (217) atoms 0.18a, 0.15c 3.19a 4.19a 3.25a

Eact Ef(P–V)

0.17 3.05 4.26 3.21

Eb(P–V)

Present DFT calculations. Ref. [74]. Ref. [75]. Present MEAM calculations.

a b c d

where Etot ½N  Si þ P i  is the total energy of a system with N Si atoms and a neutral interstitial P atom (P i , Kröger–Vink notation [73]), Etot ½N  Si is the total energy of perfect Si system. DFT and MEAM calculations about point defects are listed in Table 8. The formation energy of substitution is given higher weight during the optimization process, thus a fairly good MEAM prediction is obtained. The MEAM predictions for interstitial formation energies are also in fairly good agreement with reference values. We can see that both DFT and MEAM calculation show that Hexagonal form is the most stable interstitial arrangement. 4.4. P–vacancy pair and migration of P P–vacancy pair (P–V), also known as E-centre, consists of a substitutional P atom and a neighboring Si vacancy. It is the most common observable defects introduced by electron irradiation of P-doped floating-zone Si. In the fabrication of Si-based components, P–V is of great technological significance as the lifetime of the charge carriers is largely controlled by the interactions between P impurities and vacancies [76]. Therefore, higher weights are given to P–V related properties. Currently, only data for firstnearest-neighbor (1NN) P–V are accessible in published literatures. In this work, second- and third-nearest-neighbor (2NN and 3NN) P–V pairs are considered. The formation energy of P–V is given by [77]:

Ef ðP—VÞ ¼ EðP—VÞ þ

1 Ebulk  EðPs Þ þ ðEV þ le ÞQ N

ð17Þ

where E(P–V) is the total energy of a system with a P–V pair, Ebulk is the total energy of perfect Si system with N atoms, E(Ps) is the total energy of a Si system containing a substitutional P atom, EV denotes the energy of the top of the valence band in perfect Si lattice, le is the chemical potential (Fermi level) measured relative to EV, and Q is the charge state of the defect, usually 0 or 1. Only 0 is considered for Q in present work for MEAM model cannot take charge into account. The diffusion activation energy of P–V Eact is defined as the sum of the formation energy Ef and the migration energy Em. MEAM prediction for Em is 0.967 eV, and Eact 3.12 eV. The binding energy of the P atom to the vacancy Eb(P–V) is calculated by:

Eb ðP—VÞ ¼ EðPs Þ  EðP—VÞ  Ebulk þ EðVÞ

ð18Þ

where E(V) is the total energy of a Si system with a vacancy, the other energies are the same as above. Above introduced three properties for P–V are calculated using MEAM, in comparison with DFT and experimental values, as shown in Table 9. One can see an excellent agreement between present calculation and DFT or experimental data. The MEAM predicted formation and binding energy of 1NN P–V pair are respectively smaller and larger than those of 2NN and 2NN P–V pairs, indicating

e f h

a

1NN 2NN 3NN 1NN 2NN 3NN

2.5–3.0 2.2a, 2.5b, 2.0c, 2.39d, 2.27e 3.00d, 2.97e 3.01d, 3.09e 0.6-1.4a, 1.0b, 1.8c, 1.21d, 1.28e 0.60d, 0.59e 0.59d, 0.47e

Expt.f

MEAMh

3.1 – – – P1.04 – –

3.12 2.15 3.14 3.26 1.10 0.09 0.01

DFT calculations from Ref. [75]. DFT calculations from Ref. [77]. DFT calculations from Ref. [78]. Present DFT calculations for 63-atom systems. Present DFT calculation for 215-atom systems. Experimental data from Ref. [76]. Present MEAM calculations.

that 1NN P–V pair is the most favored site, which agrees with our DFT results. 4.5. Transformation from non-equilibrium to ordered state and thermal stability of Si7P In addition to above mentioned important physical properties, another two significant criteria for the robustness of the newly developed 2NN MEAM potential for Si–P binary system are transformation from non-equilibrium to ordered structure [26,34,79,80] and the thermal stabilities [53] of the typical Si–P crystal compounds. As mentioned in Section 4.1, the concentration of substitutional P in SCS can reach up to 10%. A 5  5  5 Si7P model (substitutional concentration is 12.5%, see Fig. 1(b)) is used to verify the abovementioned two capacities of the newly developed MEAM potential. The parameter set directly given by GA might not render Si7P thermally stable. To address this problem, we have to manually find out the critical factors that affect the thermal stability, then modify them and rerun the GA to get the final parameter set. In the second running, the critical factors are fixed. We find that Cmax121 and the reference structure of 2NN MEAM potential for pure P are critical factors. Larger Cmax121 value and diamond structure are the solution. The non-equilibrium structures are modeled by shifting every atom D (0.8 Å) away from the original position in random direction. Fig. 7 shows the transition of Si7P from non-equilibrium to ordered state through two distinct equilibration processes. The first is energy minimization, and the second is MD relaxation on the NPT ensemble, which is conducted under specific 1300 K and 1 bar for 20 ps then cooled down to 0 K. The cohesive energy and lattice constants of the final crystalline structures are identical to the original ordered structures, indicating the good robustness of the newly developed potential. Above steps are repeated 30 times, all resulting in crystal state. The RDF curves in Fig. 7(a)–(c) are averaged from the repeated processes. It’s worth noting that in Fig. 7(c) the principal peaks are accompanied by accessory peaks, which is caused by the slight offset of Si atoms around P atoms, indicating that P doping introduces slight lattice distortion. 4.6. Effects of P on the uniaxial tensile properties of bulk Si P impurities have crucial effects not only on the electrical characteristics but also on the mechanical properties of Si, like elastic constants [70], dislocation mobility [81,82], indentation fracture behavior [1], fracture strength [83], fatigue failure in atmospheric environment [84] and so on. This work only demonstrates how

10

B. Liu et al. / Computational Materials Science 120 (2016) 1–12

1.5

8

(a)

120

(b)

(c)

100

0.5

RDF g(r)

RDF g(r)

RDF g(r)

6

1

4

80 60 40

2 20

0

0

0

1

2

3

4

5

6

0

1

2

3

4

5

0

6

0

1

2

R

R

Si-P distance (A)

(d)

3

4

5

6

R

Si-P distance (A)

Si-P distance (A)

(e)

(f)

Fig. 7. RDF(Si–P) curves of Si7P during transition from non-equilibrium to ordered arrangement (a) Before energy minimization, (b) After energy minimization of 102, (c) The crystal state after the relaxation controlled by NPT, and corresponding snapshots (d)–(f). Pink and blue respectively represent Si and P.

the P impurities affect the tensile properties of bulk Si, while the rest of the above-mentioned points will be investigated in our future work. Two aspects are concerned, namely concentration and distribution form. Regarding to the first one, three models with different concentrations, 0.0578%, 0.463% and 1.563%, are considered, in which P atoms are evenly distributed. Periodic boundary conditions (PBC) are prescribed in all three dimensions (x-[1 0 0], y-[0 1 0], z-[0 0 1]) of the simulation cell (65.17 Å  65.17 Å  32.59 Å, 6912 atoms). As for the second, three models with the same concentration (0.463%) but different distribution forms (or different regional concentrations) are applied in the simulations, as the schematic diagram inserted in Fig. 7(b). Constant temperature 300 K and pressure 1bars are controlled by NPT ensemble. Uniaxial tensile direction is along x, with loading rate 5.0  107 s1. The stresses in y- and z-dimension are monitored to fluctuate around 0, in order to make sure that the tension is uniaxial. Three rates, 5.0  107 s1, 1.0  108 s1 and 1.0  109 s1, are tested to verify the effects of strain rate. Results show that it doesn’t change the tensile stiffness but changes the fracture stress and corresponding strain, as shown in supplementary Figs. S1 and 2. Stress–strain curves are plotted in Fig. 8, showing that every tensile fracture is brittle. Three conclusions are drawn from Fig. 8(a). First, the fracture threshold decreases as the concentration increases, in accordance

40

30

(a) Stress (GPa)

20

(b) data1 l=16.29 Å l=10.86 Å data2 l= 5.43 Å data3

25

0% 0.058% 0.463% 1.563%

30

Stress (GPa)

with experimental phenomena [70]. Second, the stiffness decreases as the concentration increases, which also agrees with experiments [1,70] and can be reflected by the effects of P on the elastic constants. Third, when the concentration is relatively high, such as 0.463%, the incremental stress discontinuity can be observed, e.g. the sudden indenting at the strain of 0.19. A similar phenomenon has been reported in the nanoidentation behavior of Si [85], which is called ‘‘pop-in”. By observing the evolution of the configuration when the sudden indenting happens, we find that the P atoms are noticeably displaced from their lattice positions to release local high strain energy. This phenomenon exists in each of the above mentioned tests for strain rate effects, indicating that this is an intrinsic behavior and not caused by unreasonably high strain rate. From Fig. 8(b), one can see that regional high concentration decreases the fracture threshold while has no effect on the stiffness. The lattice distortion introduced by P atoms spreads out as strain increases. The radius of the largest distortion region (visible to naked eye) right before fracture is about 7 Å. When P atoms are close enough, the distortion regions overlap and interact, making it easier to fracture. This work has only studied the effect of P impurities on bulk Si. In our follow-up research, the effects on crystal and amorphous film Si will be presented. The bond strengths of P–P, Si–P and Si– Si calculated with current MEAM are, respectively, 413.87 kJ/mol, 371.09 kJ/mol, 302.99 kJ/mol, in good agreement with

Pop-in

20 15 10

10 5 l

0

0 0

0.05

0.1 0.15 Strain,

0.2

0.25

0

0.05

0.1 0.15 Strain,

0.2

0.25

Fig. 8. Uniaxial tensile test curves of (a) different concentrations and (b) different distribution forms (regional distribution concentrations), where l denotes the distance between P atoms.

B. Liu et al. / Computational Materials Science 120 (2016) 1–12

experimental values 489.5 ± 10.5 kJ/mol, 363.6 kJ/mol and 325 ± 7 kJ/mol [86]. This indicates that the newly developed potential can be applied to investigate the mechanism that P impurities affect the mechanical properties of Si, including fatigue (our preliminary results show that P impurities can effectively impede the crack initiation on surface and its propagation). Clearly, more studies are required to confirm the suitability of this potential, which will be discussed in our next paper. 5. Conclusions (1) The parameter sets of 2NN MEAM potentials for pure P and Si–P binary system are obtained by GA based on various fundamental physical properties. It is worth mentioning that the van der Waals interactions are omitted during the parameterization for pure P, and that special attention is paid to the thermal stability of Si7P when developing the parameters for Si–P system. The physical properties predicted by present MEAM potentials are overall in fairly good agreement with reference values. (2) Further discussions about the dopant dependent tensile properties are demonstrated, including the concentration and distribution form of P impurities. It is found that higher concentration leads to lower stiffness and fracture threshold, and higher regional concentration leads to lower tensile fracture threshold. This work has investigated the effects of distribution form preliminarily, which will be studied comprehensively in our future work. (3) In present work, different atomic arrangements such as bulk, point defects and P–V pair are considered to obtain a robust potential for Si–P binary system. The validated potential can be used for a wide range of modeling and simulation for Si–P compounds such as the mechanisms that the P impurities affect the cracking and fracture behavior, dislocation mobility, hardness, brittleness of Si crystals. Moreover, in consideration of the good description for surface and amorphous properties by Lee’s 2NN MEAM potential, the newly developed potential can be applied for studying the effects of P impurities on interface behaviors, surface adhesion, and the mechanical properties of Si-based nano wires and amorphous structures. In summary, present work has provided a useful tool for researches into the nano-structural mechanical properties of Si–P solids, which are not experimentally accessible.

Acknowledgements This work is financially supported by National Natural Science Foundation of China (No. 51175503), Doctoral Innovation Foundation of NUDT (No. B140305) and China Scholarship Counsil (No. 201403170367). Part of the computation is done by West Grid. We gratefully acknowledge the helpful instructions of GA techniques from Dr. Yubao Zhen, Department of Astronautics Science and Mechanics. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.commatsci.2016. 04.002. References [1] Z.D. Zeng, X.Y. Ma, J.H. Chen, Y.H. Zeng, D.R. Yang, Y.G. Liu, J. Appl. Phys. 107 (2010).

11

[2] Y. Lee, G.S. Hwang, Phys. Rev. B 86 (2012). [3] H.D. Chiou, J. Electrochem. Soc. 147 (2000) 345–349. [4] H. Siethoff, H.G. Brion, Mater. Sci. Eng. A – Struct. Mater. Prop. Microstruct. Process. 355 (2003) 311–314. [5] Y. Ohno, T. Shirakawa, T. Taishi, I. Yonenaga, Appl. Phys. Lett. 95 (2009). [6] H.F. Wilson, N.A. Marks, D.R. McKenzie, K.H. Lee, Nucl Instrum. Methods Phys. Res. Sec. B – Beam Interact. Mater. Atoms 215 (2004) 99–108. [7] G.J. Ackland, M.I. Mendelev, D.J. Srolovitz, S. Han, A.V. Barashev, J. Phys. – Condens. Matter 16 (2004) S2629–S2642. [8] J.A. Hauch, D. Holland, M.P. Marder, H.L. Swinney, Phys. Rev. Lett. 82 (1999) 3823–3826. [9] T. Mattila, R.M. Nieminen, M. Dzugutov, Nucl. Instrum. Methods Phys. Res. Sect. B – Beam Interact. Mater. Atoms 102 (1995) 119–124. [10] T. Mattila, R.M. Nieminen, M. Dzugutov, Phys. Rev. B 53 (1996) 192–200. [11] H.W. Sheng, E. Ma, M.J. Kramer, JOM 64 (2012) 856–881. [12] W.S. Ko, N.J. Kim, B.J. Lee, J. Phys. – Condens. Matter 24 (2012). [13] N.P. Bailey, J.P. Sethna, Phys. Rev. B 68 (2003). [14] K. Kang, W. Cai, Phil. Mag. 87 (2007) 2169–2189. [15] J.G. Swadener, M.I. Baskes, M. Nastasi, Phys. Rev. Lett. 89 (2002). [16] S. Ryu, C.R. Weinberger, M.I. Baskes, W. Cai, Modell. Simul. Mater. Sci. Eng. 17 (2009). [17] J.G. Swadener, M.I. Baskes, M. Nastasi, Phys. Rev. B 72 (2005). [18] B.J. Lee, CALPHAD – Comput. Coupling Phase Diagrams Thermochem. 31 (2007) 95–104. [19] B. Liu, H. Zhang, J. Tao, X. Chen, Y.a. Zhang, Comput. Mater. Sci. 109 (2015) 277–286. [20] M.I. Baskes, J.S. Nelson, A.F. Wright, Phys. Rev. B 40 (1989) 6085–6100. [21] M.S. Daw, M.I. Baskes, Phys. Rev. B 29 (1984) 6443–6453. [22] M.I. Baskes, Phys. Rev. B 46 (1992) 2727–2742. [23] B.J. Lee, M.I. Baskes, Phys. Rev. B 62 (2000) 8564–8567. [24] B.J. Lee, W.S. Ko, H.K. Kim, E.H. Kim, CALPHAD – Comput. Coupling Phase Diagrams Thermochem. 34 (2010) 510–522. [25] M.I. Baskes, J.E. Angelo, C.L. Bisson, Modell. Simul. Mater. Sci. Eng. 2 (1994) 505–518. [26] Z.W. Cui, F. Gao, Z.H. Cui, J.M. Qu, J. Power Sources 207 (2012) 150–159. [27] B. Jelinek, S. Groh, M.F. Horstemeyer, J. Houze, S.G. Kim, G.J. Wagner, A. Moitra, M.I. Baskes, Phys. Rev. B 85 (2012). [28] J.H. Rose, J.R. Smith, F. Guinea, J. Ferrante, Phys. Rev. B 29 (1984) 2963– 2969. [29] C.L. Kuo, P. Clancy, Modell. Simul. Mater. Sci. Eng. 13 (2005) 1309–1329. [30] J. Osugi, R. Namikawa, Y. Tanaka, Rev. Phys. Chem. Jpn. 36 (1966) 35–43. [31] W.S. Ko, B.J. Lee, Phil. Mag. 94 (2014) 1745–1753. [32] B.J. Lee, Acta Mater. 54 (2006) 701–711. [33] T.J. Lenosky, B. Sadigh, E. Alonso, V.V. Bulatov, T.D. de la Rubia, J. Kim, A.F. Voter, J.D. Kress, Modell. Simul. Mater. Sci. Eng. 8 (2000) 825–841. [34] Z.W. Cui, F. Gao, Z.H. Cui, J.M. Qu, Modell. Simul. Mater. Sci. Eng. 20 (2012). [35] S. Plimpton, J. Comput. Phys. 117 (1995) 1–19. [36] G. Kresse, D. Joubert, Phys. Rev. B 59 (1999) 1758–1775. [37] G. Kresse, J. Furthmuller, Phys. Rev. B 54 (1996) 11169–11186. [38] J.P. Perdew, J.A. Chevary, S.H. Vosko, K.A. Jackson, M.R. Pederson, D.J. Singh, C. Fiolhais, Phys. Rev. B 46 (1992) 6671–6687. [39] J. Park, B.D. Yu, S. Hong, J. Korean Phys. Soc. 59 (2011) 196–199. [40] H.J. Mcskimin, P. Andreatch, J. Appl. Phys. 35 (1964) 3312. [41] C. Kittel, Introduction to Solid State Physics, fourth ed., Wiley, New York, 1971. xv, 766p. [42] R. Ahuja, Phys. Status Solidi B – Basic Res. 235 (2003) 282–287. [43] A.T. Dinsdale, CALPHAD – Comput. Coupling Phase Diagrams Thermochem. 15 (1991) 317–425. [44] H. Thurn, H. Krebs, Acta Crystallogr. Sect. B – Struct. Crystallogr. Crystal Chem. B 25 (1969) 125–135. [45] D.K. Seo, R. Hoffmann, J. Solid State Chem. 147 (1999) 26–37. [46] B.J. Lee, J.W. Lee, CALPHAD – Comput. Coupling Phase Diagrams Thermochem. 29 (2005) 7–16. [47] Y. Akahama, H. Kawamura, S. Carlson, T. Le Bihan, D. Hausermann, Phys. Rev. B 61 (2000) 3139–3142. [48] D.R. Stull, H. Prophet, JANAF Thermochemical Tables, National Bureau of Standards, Washington, DC, 1971. [49] L. Cartz, S.R. Srinivasa, R.J. Riedner, J.D. Jorgensen, T.G. Worlton, J. Chem. Phys. 71 (1979) 1718–1721. [50] J.C. Jamieson, Science 139 (1963) 1291–1292. [51] E. Wiberg, N. Wiberg, Inorganic Chemistry, Academic, San Diego, CA, 2001. [52] L. Gerward, J. Phys. Chem. Solids 46 (1985) 925–927. [53] B.J. Lee, T.H. Lee, S.J. Kim, Acta Mater. 54 (2006) 4597–4607. [54] Z.X. Hui, Mech. Electron. Eng. Technol. (Icmeet 2014) 538 (2014) 32–35. [55] R. Hill, Proc. Phys. Soc. London Sect. B 65 (1952). 396-396. [56] T. Kikegawa, H. Iwasaki, Acta Crystallogr. Sect. B – Struct. Sci. 39 (1983) 158– 164. [57] T. Sasaki, K. Shindo, K. Niizeki, A. Morita, J. Phys. Soc. Jpn. 57 (1988) 978–987. [58] F.J.H. Ehlers, N.E. Christensen, Phys. Rev. B 69 (2004). [59] J.S. Custer, M.O. Thompson, D.C. Jacobson, J.M. Poate, S. Roorda, W.C. Sinke, F. Spaepen, Appl. Phys. Lett. 64 (1994) 437–439. [60] K.D. Weeks, S.G. Thomas, P. Dholabhai, J. Adams, Thin Solid Films 520 (2012) 3158–3162. [61] M.E. Schlesinger, Chem. Rev. 102 (2002) 4267–4301. [62] S.M. Liang, R. Schmid-Fetzer, J. Phase Equilib. Diffus. 35 (2014) 24–35. [63] T. Wadsten, Chem. Scr. 8 (1975) 63–69.

12

B. Liu et al. / Computational Materials Science 120 (2016) 1–12

[64] T. Wadsten, Acta Chem. Scand. 21 (1967) 593. [65] T. Wadsten, Acta Chem. Scand. 21 (1967) 1374. [66] R.W. Olesinski, N. Kanani, G.J. Abbaschian, Bull. Alloy Phase Diagrams 6 (1985) 130–133. [67] E. Kooi, J. Electrochem. Soc. 111 (1964) 1383–1387. [68] G. Fritz, H.O. Berkenhoff, Z. Anorg. Allg. Chem. 300 (1959) 205–209. [69] P. Ravindran, L. Fast, P.A. Korzhavyi, B. Johansson, J. Wills, O. Eriksson, J. Appl. Phys. 84 (1998) 4891–4904. [70] J.J. Hall, Phys. Rev. 161 (1967) 756. [71] R.N. Thurston, K. Brugger, Phys. Rev. 133 (1964) 1604. [72] K. Brugger, Phys. Rev. 133 (1964) 1611. [73] F.A. Kröger, H.J. Vink, Relations between the concentrations of imperfections in crystalline solids, in: S. Frederick, T. David (Eds.), Solid State Physics, Academic Press, 1956, pp. 307–435. [74] O. Sugino, A. Oshiyama, Phys. Rev. B 46 (1992) 12335–12341. [75] R. Virkkunen, R.M. Nieminen, Comput. Mater. Sci. 1 (1993) 7. [76] P.M. Fahey, P.B. Griffin, J.D. Plummer, Rev. Mod. Phys. 61 (1989) 289–384.

[77] C.S. Nichols, C.G. Van De Walle, S.T. Pantelides, Phys. Rev. B 40 (1989) 5484– 5496. [78] R. Car, P.J. Kelly, A. Oshiyama, S.T. Pantelides, Phys. Rev. Lett. 52 (1984) 1814– 1817. [79] S. Zhang, N.X. Chen, Phys. Rev. B 66 (2002). [80] C. Wang, S. Zhang, N.X. Chen, J. Alloys Compd. 388 (2005) 195–207. [81] J.R. Patel, L.R. Testardi, P.E. Freeland, Phys. Rev. B 13 (1976) 3548–3557. [82] I. Yonenaga, Mater. Sci. Eng. B – Solid State Mater. Adv. Technol. 124 (2005) 293–296. [83] B. Liu, J.Y. Tao, Y.A. Zhang, X. Chen, X.J. Wang, Micro Nano Lett. 8 (2013) 726– 730. [84] Y.A. Zhang, J.Y. Tao, Y.L. Wang, Z.Q. Ren, B. Liu, X. Chen, Microelectron. Reliab. 53 (2013) 1672–1675. [85] Y.B. Gerbig, S.J. Stranick, D.J. Morris, M.D. Vaudin, R.F. Cook, J. Mater. Res. 24 (2009) 1172–1183. [86] J.A. Kerr, CRC Handbook of Chemistry and Physics 1999–2000, 81st ed., CRC Press, Florida, 2000.