Parameters setup for first-principles on well-formedness and thermodynamic properties of bulk and habit surfaces of Mg2Si

The well-formedness, thermodynamic properties and mechanical properties of Mg2Si bulk and Mg2Si (1 0 0), (1 1 0) and (1 1 1) surfaces were calculated by using first-principles based on density functional theory (DFT). Some accurate and efficient parameters were programed by computing massively and repeatedly. The accurate and efficient parameters used in Mg2Si bulk is Energy cutoff (500 eV), k-point (5 × 5 × 5) and exchange-correlation interaction (PBEsol). The relaxation of surficial atoms shows that the rule or rate of expansion and shrink are similar from the second column. The rate of the innermost layer decreases gradually with increasing slab thickness. The rate of the innermost layer in Mg-terminated or Si-terminated Mg2Si (1 0 0) and (1 1 0) is no more than 0.1% when convergent layers are no less than 11. The convergent layers of Mg2-terminated and Si-terminated Mg2Si (1 1 1) are 14 and 16, respectively. The Mg1-terminated Mg2Si (1 1 1) surfaces hardly converge even if the layers reach to 18. Si-terminated surfaces are harder to keep stability than Mg-terminated surfaces in Mg2Si (1 0 0) and (1 1 1) surfaces because of a higher surface energy. Mg2Si (1 1 0) surface energy is a constant because of the nonpolar and stoichiometric Mg-Si terminated surface. Si-terminated and Mg-terminated Mg2Si (1 0 0) surfaces are more stable than Si-terminated and Mg2-terminated Mg2Si (1 1 1) surfaces, respectively. Mg1-terminated Mg2Si (1 1 1) surface is the most stable surface over the entire surfaces in Mg2Si.


Introduction
Mg 2 Si attracted a considerable attention in recent years because of its narrow band gap as n-type semiconductors, which has been applied to eco-friendly thermoelectric material at mid-temperature (600-900 K) [1][2][3]. The electron mobility of Mg 2 Si is significantly higher than vacancy mobility, that will become a potential function material exhibiting excellent electrical performance [4,5]. In addition, because of their inherent high-melting, low density, high hardness, low thermal expansion coefficient (TEC) and high elastic modulus, Mg 2 Si is used as reinforced particles in Mg-or Al-alloys [6]. Therefore, some kinds of Mg-or Al-metal matrix composites (MMCs), which exhibit high ductility and high strength and hardness at the same time. However, brittle phase Mg 2 Si always has lots of sharp corners, which will dissever matrix in its application [7,8]. Therefore, modifying the morphology of reinforced phase Mg 2 Si will improve strength and toughness observably. The development of reinforced MMCs will be no limit after eliminating the danger caused by the morphology of Mg 2 Si.
The morphology of crystal covers micro-morphology and macro-morphology. The macro-morphology is determined by the growing environment of crystal, but the intrinsic characteristics define the micromorphology. Therefore, it is significant to study on the intrinsic characteristics of Mg 2 Si, such as mode of nucleation, stack method of atoms, surface properties of crystals. Especially, according to Gibbs-Wulff theorem [9], crystal will adjust its morphology to a status that exhibits the minimum surface free energy in thermodynamic equilibrium. However, the intrinsic characteristics of crystal are hardly characterized by Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. experimental means. Using first-principles, which does not depend on any experience parameters, is more efficient and accurate to reveal the configuration and surface properties of Mg 2 Si surfaces [10]. Lots of strong evidences will be provided on study of the initial micro-morphology of Mg 2 Si emerging in solution.
The configuration and surface properties of Mg 2 Si (1 0 0) surface has been performed by Liao [11] and Niranjan [12], respectively, up to now. Some other habit surfaces with low Miller indices has not been performed widely. In present investigation, some exploration has been done to reveal the well-formedness and thermodynamic properties of Mg 2 Si bulk and surfaces. In addition, some simulated results have great error, but some consume significant computer resources because of aspiring to over-accurate parameters setup. Imai et al [3] simulated by using an energy cutoff of 380 eV and a grid of 3×3×3, so that the lattice didn't converge better; and the error of lattice and formation enthalpy is greater than most work. Wang et al [13] suggested that calculating the lattice and bulk modulus used an energy cutoff of 1360 eV and a grid of 8×8×8, which will consume significant computer resources. However, the error didn't be reduced because of using LDA functional. Therefore, the current work compute massively and repeatedly in order to provide a suit of accurate and efficient parameters. Researchers can further optimized referring to these parameters setup.

Calculation method and details
The calculations in this work were all performed by density functional theory (DFT) using the Cambridge Serial Total Energy Package (CASTEP) Code [14,15]. The interaction between the ionic core and the valence electrons was described by the plane-wave ultra-soft pseudopotential method [16]. Furthermore, the 3s 2 for Mg and 3s 2 3p 2 for Si were chosen as the valence electrons of the atoms. The PW91 (Perdew-Wang) [17], PBE (Perdew-Burke-Enzerhof) [18] and PBEsol (PBE solid) [19] functional of GGA (generalized gradient approximation) and LDA (local-density approximation) [20,21] were utilized to treat the exchange-correlation interactions. The Brillouin zone was sampled by means of the Monkhorst-Pack k-point grid using a grid from 2×2×2 to 11×11×11 for Mg 2 Si bulk, and also a suitable grid for different supercell slabs of Mg 2 Si surfaces. A kinetic energy cutoff value, which is dependent on elements Mg and Si and independent on cell structure, was filtrated from 300 eV to 800 eV for both Mg 2 Si bulk and surfaces. For obtaining well-formedness of bulk and surfaces, the BFGS (Broyden-Fletcher-Goldfarb-Shanno) algorithm [22] was applied to geometry optimization to reach the ground state. In addition, the vacuum layer of Mg 2 Si slabs were selected from 5 Å to 15 Å for each surface structure to eliminate the interactions between the surface atoms. The following convergence parameters were set in the calculations: the total energy tolerance was 5.0×10 −6 eV/atom, the maximum force tolerance was 0.01 eV Å, and maximal displacement was 5.0×10 −4 nm.

Results and discussions
3.1. Parameters setup, well-formedness and thermodynamic properties of Mg 2 Si bulk To guarantee the rationality and accuracy of the calculation, a suit of accurate and efficient parameters, especially energy cutoff, k-point grid and functional of exchange-correlation, for Mg 2 Si bulk calculation was designed. Therefore, convergence test should be performed by using above parameters; and then the lattice constant, modulus (elastic modulus, shear modulus, bulk modulus, Poisson's ratio) and formation enthalpy for Mg 2 Si were deduced.
Generally, the denser the k-point grid is and the greater the energy cutoff is, the more accurate the simulation of total energy will be. The results of convergence test are plotted in figure 1, which indicates that the lattice parameters of Mg 2 Si will tend to converge (using a grid of 5×5×5) when energy cutoff is above 450 eV. However, the efficiency of simulation will slump with increasing energy cutoff. It also suggests that a grid of 5×5×5 k-point will be more efficient when energy cutoff is programmed to 500 eV. Therefore, the grid smaller than 5×5×5 and the energy cutoff lower than 450 eV are considered that the precision is not enough. Imai et al indicated that a lattice parameter 7.99 Å, which show the maximum error than experiment and other calculation, was simulated using an energy cutoff of 380 eV and a grid of 3×3×3 [3]. Therefore, accurate parameters setup is extremely important.
On the process of geometry optimization, the most common functional of exchange-correlation are LDA and GGA. Generally, the PW91, PBE and PBEsol are always regarded as accurate approximation among GGA. The lattice parameter of Mg 2 Si was simulated by above approximation method. The lattice parameters of Mg 2 Si reached to 6.259 Å, when LDA was performed on CASTEP, because of an over-binding effect. Comparing to lots of calculations [3,13,23] and experiment [24,25], PBEsol exhibited considerable accuracy (refer to table 1). Because the equilibrium properties of solids and their surfaces, comparing with PBE functional, are improved further. In addition, comparing to the time of simulation, the test of PBE and PW91 took more than twice.
According to above parameters setup, the reasonable configuration and thermodynamic and mechanical properties are determined. The mechanical properties were deduced by elastic constants obtained through DFT directly; and the formation enthalpy of Mg 2 Si bulk was calculated as follows [26,27]:    [28]. With increasing the atomic layers of slab, more accurate results could be obtained, yet more hardware resources and longer times are needed. Therefore, initial convergence tests of the habit surfaces with low Miller indices, such as (1 0 0), (1 1 0) and (1 1 1) surfaces of Mg 2 Si, were performed to determine the minimum atomic layers.
In the current work, Mg 2 Si surfaces were built by cleaving the optimized Mg 2 Si bulk supercell and structural terminations consisted of Mg atom, Si atom or Mg-Si atom. The periodic boundary condition was used in the surface supercell, including a slab model and a vacuum region [29]. As shown in figure 2, the Mg 2 Si (1 0 0) plane exhibit two different terminations, such as the Mg-or Si-terminated surface. However, there are three terminations, such as Mg 1 -, Mg 2 -or Si-termination, emerging in Mg 2 Si (1 1 1) surface. The kind of terminations is determined by the order of layers. Mg 2 Si (1 1 0) has only one considerable termination, because different layers of Mg 2 Si (1 1 0) surface are Mg-Si, which exhibit same texture and properties.
No matter which kind of termination in a surface, there is the same order of atomic layer at the two ends of a slab to reduce the dipole effect. Furthermore, it is significant to understand that the Mg 2 Si slab thickness of the two polar surfaces is defined by the layer number of surface relaxation. The percentage of increase or decrease between different layers are calculated as follows: whered i j is the layer distance between i-layer and j-layer of the Mg 2 Si slab after surface relaxation, and d 0 is the distance of the structure which has been not relaxed. Therefore, Dd i j is less than Zero, when the distance of the layers decreased.
In addition, in order to prevent other atomic interaction between layers, surface models should include a vacuum region with the thickness of 10 Å. The vacuum regions were tested from 5 Å to 15 Å and the results showed that 10 Å is enough for all Mg 2 Si surfaces. Therefore, vacuum region (10 Å), energy cutoff (500 eV) and exchange-correlation interaction (PBEsol) were programed during the geometry optimization, which can relax all atoms fully with peak efficiency.
The surface relaxation of Mg 2 Si (1 0 0), (1 1 0) and (1 1 1) as a function of terminations and slab thickness are listed in table 2, where n is the number of layers. According to Gibbs-Wulff theorem [9], crystal will adjust its morphology to a status that exhibits the minimum surface free energy in thermodynamic equilibrium. In order to keep a steady state, the surface atoms of Mg 2 Si have to shift the position toward low energy. There are lots of conclusions can be drawn related to slab thickness, terminations and Mg 2 Si surfaces. Firstly, the first column of the data in table 2 has an obvious lattice distortion and irregularity, so that the fewest layers exhibit maximum dimension changes. Secondly, the rule and rate of expansion and shrink are similar from second column to the last. Therefore, the movement of atoms around Mg 2 Si surfaces are according to the numbers of layer. The last and the most, the rate of the innermost layer decreases gradually with increasing slab thickness. The relaxation occurs on different layers to the different surfaces of Mg 2 Si with different terminations. For Mg 2 Si (1 0 0) and (1 1 0) with different terminations, Δd ij of the innermost layer is no more than 0.1% when n11, which indicats that 11-layered Mg 2 Si (1 0 0) and (1 1 0) exhibit a better convergence. However, the convergence has great distinction on three terminations of Mg 2 Si (1 1 1) surfaces. The number of convergent layers are 14-layers and 16-layers in Mg 2 -terminated and Si-terminated Mg 2 Si (1 1 1) surfaces, respectively. The Mg 1 -terminated Mg 2 Si (1 1 1) surfaces still hardly reach to an ideal value, even if the layers reach to 18. Above all, the parameters setup of the atomic layers was defined according to the different Mg 2 Si surfaces with different terminations.

Surface energy of Mg 2 Si surfaces
According to the lowest energy principle, crystal will tend to a state that exhibit the lowest surface free energy. Therefore, the stability of different terminated atoms and different surfaces can also be indicated by the surface energy of the compounds.  Si (1 0 0), (1 1 0) and (1 1 1)  The chemical potentials of Mg and Si atoms should be taken into account during the calculation of surface energy. The surface energy of Mg 2 Si (1 0 0), Mg 2 Si (1 1 0) and Mg 2 Si (1 1 1) with different terminations is determined from the following equation [30]: The chemical potential of each element is smaller than that of the corresponding bulk substances. Even more, combining the description as equation (6), the relationship between chemical potentials can be described as follow: The surface energy of all above surfaces with different terminations were calculated according to equation (5), and One more accurate formation enthalpy (−20.15 kJ mol −1 ) obtained by GGA-PBEsol was selected for discussing of surface energy. Finally, the relationship between surface energy and chemical potential were plotted in figure 3. The surface energy is not a constant in the polar surface, such as Si-terminated and Mgterminated Mg 2 Si (1 0 0) and Si-terminated and Mg 2 -terminated Mg 2 Si (1 1 1). The surface energy will shift with the changes of chemical composition in a system. The upper limit of Mg chemical potential, where As a conclusion, Si-terminated surfaces always exhibit a higher surface energy, so that Si-terminated surfaces are hard to keep stability in Mg 2 Si (1 0 0) and (1 1 1) surfaces. In addition, Si-terminated and Mg-terminated Mg 2 Si (1 0 0) surfaces are more stable than Si-terminated and Mg 2 -terminated Mg 2 Si (1 1 1) surfaces. Mg 2 Si (1 1 0) surface energy exhibits a constant value, which is smaller than above four terminated surfaces, because of the nonpolar and stoichiometric Mg-Si terminated surface. The Mg 1 -terminated Mg 2 Si (1 1 1) surface is the most stable surface over the entire surfaces of Mg 2 Si. Although Mg 1 -terminated Mg 2 Si (1 1 1) surface may form at the initial stage of micro-morphology emerging in solution, the Mg 2 Si (1 1 1) surface is still hardly piled up because of the highest energy of Si-terminated and Mg 2 -terminated Mg 2 Si (1 1 1) surfaces. Therefore, it is great significant to expose the stack sequence of atoms by researching and discussing the surface properties of Mg 2 Si in depth.