Calculation of Surface Properties of Cubic and Hexagonal Crystals through Molecular Statics Simulations

Abstract

Surface property is an important factor that is widely considered in crystal growth and design. It is also found to play a critical role in changing the constitutive law seen in the classical elasticity theory for nanomaterials. Through molecular static simulations, this work presents the calculation of surface properties (surface energy density, surface stress and surface stiffness) of some typical cubic and hexagonal crystals: face-centered-cubic (FCC) pure metals (Cu, Ni, Pd and Ag), body-centered-cubic (BCC) pure metals (Mo and W), diamond Si, zincblende GaAs and GaN, hexagonal-close-packed (HCP) pure metals (Mg, Zr and Ti), and wurzite GaN. Sound agreements of the bulk and surface properties between this work and the literature are found. New results are first reported for the surface stiffness of BCC pure metals, surface stress and surface stiffness of HCP pure metals, Si, GaAs and GaN. Comparative studies of the surface properties are carried out to uncover trends in their behaviors. The results in this work could be helpful to the investigation of material properties and structure performances of crystals.

1. Introduction

In crystal growth and design, surface property could affect surface morphology [1], surface reconstruction [2], growth rate [3], adatom absorption [4], etc. For nanomaterials, it is also found to play a critical role in changing conventional bulk properties, such as elastic modulus [5], piezoelectricity [6], thermal conductivity [7], etc. Thus, surface property is important for the determination of material properties and structure performances [8,9,10,11]. The concept of surface free energy (called surface energy for short) was first introduced in fluid mechanics for studying surface phenomena. It is equal to the reversible work per unit area needed to create a surface, which is also widely referred to as surface tension in fluid mechanics. For solids, a new area can be created by a process such as cleavage, but another surface quantity, surface stress, could also critically affect the behavior of the surface, which is defined as the reversible work per unit area needed to elastically stretch a pre-existing surface [12]. In this case, the change in surface energy should be equal to the work done by the surface stress, as it deforms the surface area in the Lagrangian description of the well-known Shuttleworth relation [13]. Due to the difficulties of obtaining surface properties by experiments, they were widely investigated by atomistic simulations, such as first principle (FP) calculations or molecular statics (MS) calculations with empirical potentials. Although MS calculations are usually less accurate than FP calculations, they remain useful in understanding trends and developing concepts. Earlier, using embedded-atom-method (EAM) potentials [14], Foiles et al. [15] presented surface energies of low-index surfaces of some pure face-centered-cubic (FCC) metals, such as Cu, Ag, Au, Ni, Pd and Pt. Finnis and Sinclair [16] proposed simple empirical N-body Finnis–Sinclair (FS) potentials and calculated surface energies of some pure body-centered-cubic (BCC) metals, such as V, Nb, Ta, Cr, Mo, W and Fe. Afterwards, surface stresses of pure metals were also investigated by different researchers. Typically, Gumbsch and Daw [17] calculated the isotropic surface stresses of FCC(100) and FCC(111) surfaces of Cu, Ag, Au, Ni and Pt with EAM potentials. Feibelman [18] investigated the anisotropic surface stresses of FCC(110) surfaces of Pd and Pt through FP calculations. On the other hand, with FS potentials, Ackland and coworkers calculated the surface stresses of BCC(100), BCC(110) and BCC(111) surfaces of V, Nb, Ta, Mo and W [19,20], as well as the basal (0001) and prismatic $(10\overline{1}0)$ and $(11\overline{2}0)$ surfaces of hexagonal-close-packed (HCP) crystals like Ti [21] and Zr [22]. Moreover, surface energy and surface stress were also investigated for semiconductor materials. For diamond-structure crystals, surface energies of Si(100) and Si(111) surfaces were presented by Northrup and coworkers [23,24] through FP calculations, while the surface stresses of Si(111) and Ge(111) surfaces were given by Meade and Vanderbilt [25,26] through FP calculations. However, for compound semiconductor materials, FP calculations were mostly performed to obtain the surface energy only, such as (100), (110) and (111) surfaces of zincblende GaAs [27,28,29], and the basal (0001) and prismatic $(10\overline{1}0)$ and $(11\overline{2}0)$ surfaces of wurzite GaN [3,30].

Nevertheless, surface stress was treated as a constant quantity in aforementioned works. Referring back to the concept of surface energy and surface stress, surface energy could be generalized as a deformation-dependent quantity, so it could be expanded as Taylor series in terms of surface strain, where the constant surface stress is just the first-order coefficient, while the second-order coefficient could be defined as surface stiffness tensor similar to the bulk stiffness tensor. In fact, surface stiffness could play an important role in changing the constitutive law seen in the classical elasticity theory [31,32,33,34,35,36]. Miller and Shenoy [37] investigated the size-dependent elastic properties of nano-sized structural elements, by incorporating surface stiffness. Dingreville et al. [38] proposed to formulate the surface energy density as a quadratic function of surface strain, where surface properties (surface energy density, surface stress and surface stiffness) could be defined, and surface stiffness was demonstrated to affect the elastic behavior of nano-sized particles, wires and films. Correspondingly, Shenoy [39] presented MS calculations of surface properties of some FCC crystals in the Eulerian frame. Dingreville and Qu [40] presented a semi-analytical method to compute surface properties of some FCC crystals in the Lagrangian frame, which were given analytically in terms of the interatomic potentials, with the relaxed positions of the atoms near the free surface obtained through one single MS simulation.

Practically, the surface energy of different planes in a crystal could be used to determine the equilibrium morphology and equilibrium growth rate in the direction normal to each plane [3], while surface stress is usually invoked in crystal growth mechanism to explain and understand adatom absorption and surface reconstruction behaviors [18]. For nanomaterials, surface stiffness has been found to change the constitutive law seen in the classical elasticity theory, and the conventional bulk properties become size-dependent [5,34]. Therefore, the results of surface properties are very helpful to the investigation of material properties and structure performances of crystals. This work investigates surface properties of some typical cubic and hexagonal crystals through MS calculations. FCC pure metals (Cu, Ni, Pd and Ag), BCC pure metals (Mo and W), diamond Si, zincblende GaAs and GaN are taken as typical examples of cubic crystals, while HCP pure metals (Mg, Zr and Ti), and wurzite GaN are taken as typical examples of hexagonal crystals. Since surface properties can be regarded as the coefficients in the Taylor series of surface energy density, with respect to surface strain similar to the relationship between bulk energy density and bulk strain, elastic constants of the bulk stiffness tensor are first obtained for all crystals in this work, to validate the methodology by comparing with available results in literature. Afterwards, certain surface properties of some crystals in this work are also compared with available results in literature, to guarantee the correctness of the calculation. Sound agreements of the bulk and surface properties between this work and the literature are found. New results are presented for the surface stiffness of BCC pure metals, surface stress and surface stiffness of HCP pure metals, Si, GaAs and GaN. In this work, the surface properties are calculated for ideally cleaved surfaces, so any surface defects (adatoms, vacancies, etc.) are not considered. It should be mentioned that the methodology is in general, but the results from MS calculations rely on the choice of interatomic potentials. The readers could follow the same methodology to obtain the property of other materials according to their own interest and with their own choice of interatomic potentials.

2. Methodology

2.1. Bulk Elastic Constants

In the classical elasticity theory, the bulk energy density could be related to the bulk strain of the material as Taylor series:

$$U=U_0+L_{ij}{\epsilon }_{ij}+\frac{1}{2}{C}_{ijkl}{\epsilon }_{ij}{\epsilon }_{kl}+\dots ,$$

(1)

where U₀ is the cohesive energy density and it is not considered to contribute to the elastic deformation. L_ij is the residual stress to be zero in this work. C_ijkl is the fourth order bulk stiffness tensor of the material. Usually, the bulk stiffness tensor is adequate to describe the elastic behavior of the bulk material.

For anisotropic materials, the bulk stiffness tensor has symmetric properties as C_ijkl = C_ijlk = C_jikl = C_klij, so it has only 21 independent elastic constants, which could be expressed as a 6x6 stiffness matrix by Voigt’s compact notation. With planes of symmetry and axis of symmetry, cubic crystals have only 3 independent elastic constants: C₁₁, C₁₂ and C₄₄, while hexagonal crystals have 5 independent elastic constants (take axis-3 as the axis of symmetry): C₁₁, C₁₂, C₁₃, C₃₃ and C₆₆.

Through MS calculations, the relationship in Equation (1) between the bulk energy density and the bulk strain is fitted to obtain the bulk elastic constants. For a given crystal composed of N atoms with periodic boundary conditions in all directions, the bulk energy density is simply calculated as:

$$U=\frac{1}{V_0}{\displaystyle \sum _{i=1}^N{E}_i,}$$

(2)

where E_i is the energy of the i-th atom and V₀ is the volume of the crystal in the unstrained reference state.

The methodology is to yield a mesh of the bulk energy density in the bulk strain space. First, the initial assembly of atoms is created using the given material properties, and a self-equilibrium state of the assembly is obtained as the unstrained reference state. Afterwards, a small strain is applied to the assembly and re-equilibrate. The bulk energy density can be computed for this corresponding strain. Repeating this procedure with different strains (−1~1%), we obtain a mesh of the bulk energy density as a function of the bulk strain, and the bulk elastic constants are fitted in the end.

2.2. Surface Properties

Similar to the relationship in Equation (1) between the bulk energy density and the bulk strain, assuming the surface energy density is a smooth function of the surface strain, one may expand the surface energy density as Taylor series of surface strain [38,40]:

$$\mathsf{\Gamma }={\mathsf{\Gamma }}^{(0)}+{\mathsf{\Gamma }}_{\alpha \beta }^{(1)}{\epsilon }_{\alpha \beta }^s+\frac{1}{2}{\mathsf{\Gamma }}_{\alpha \beta \gamma \delta }^{(2)}{\epsilon }_{\alpha \beta }^s{\epsilon }_{\gamma \delta }^s+\dots$$

(3)

where Γ⁽⁰⁾ is the surface energy density at the relaxed and unstrained state of the surface. Γ⁽¹⁾ is the surface stress with ${\mathsf{\Gamma }}_{\alpha \beta }^{\left(1\right)}={\mathsf{\Gamma }}_{\beta \alpha }^{\left(1\right)}$ that exists when the surface strain is absent. Γ⁽²⁾ is the surface stiffness tensor with ${\mathsf{\Gamma }}_{\alpha \beta \gamma \delta }^{\left(2\right)}={\mathsf{\Gamma }}_{\beta \alpha \gamma \delta }^{\left(2\right)}={\mathsf{\Gamma }}_{\alpha \beta \delta \gamma }^{\left(2\right)}={\mathsf{\Gamma }}_{\gamma \delta \alpha \beta }^{\left(2\right)}$. It should be noted that the surface is a two-dimensional object and the bulk is three-dimensional. Thus, the Greek indices take the value of 1 or 2 in Equation (3), while Latin subscripts adopt values from 1 to 3 in Equation (1). Together, Γ⁽⁰⁾, Γ⁽¹⁾ and Γ⁽²⁾ are the surface properties that fully characterize the linear elastic behavior of the surface.

Due to the symmetrically atomic arrangements along different in-plane directions of the surface, the numbers of independent surface elastic constants could be reduced. As shown in Figure 1, for cubic (001), cubic (111) and hexagonal (0001) surfaces, the atomic arrangements are symmetric along the two in-plane directions, so ${\mathsf{\Gamma }}_{11}^{\left(1\right)}={\mathsf{\Gamma }}_{22}^{\left(1\right)}$ and ${\mathsf{\Gamma }}_{1111}^{\left(2\right)}={\mathsf{\Gamma }}_{2222}^{\left(2\right)}$, while the atomic arrangements in cubic (110), hexagonal $(10\overline{1}0)$ and $(11\overline{2}0)$ surfaces are distinct along the two in-plane directions, so ${\mathsf{\Gamma }}_{11}^{\left(1\right)}\ne {\mathsf{\Gamma }}_{22}^{\left(1\right)}$ and ${\mathsf{\Gamma }}_{1111}^{\left(2\right)}\ne {\mathsf{\Gamma }}_{2222}^{\left(2\right)}$.

Through MS simulations, the relationship in Equation (3) between the surface energy density and the strain is fitted to obtain the surface properties. It should be noted that surface thickness is neglected in Equation (3), based on the concept of a dividing surface, but the surface energy is actually contributed from all atoms within a few atomic layers near the surface. If there are N atoms within these atomic layers associated with the surface, the surface energy density is defined as the summation of individual excess energy between a surface atom and the bulk atom deep in the interior of a large crystal [40]:

$$\mathsf{\Gamma }=\frac{1}{A_0}{\displaystyle \sum _{i=1}^N(E_i-E_0),}$$

(4)

where E_i is the energy of the i-th surface atom, E₀ is the energy of any bulk atom deep in the interior of a large crystal, and A₀ is the surface area in the relaxed and unstrained reference state.

The methodology to obtain surface properties is similar to that of bulk elastic constants, with only a few differences. Periodic boundary conditions are only applied to the two in-plane surface directions of a given crystal, while a fixed boundary condition is imposed in the third direction where the top and bottom surfaces locate. There should be enough layers of atoms in the vertical direction to avoid the interaction between the two surfaces. First, the initial assembly of atoms is created using the given material properties, and a self-equilibrium relaxed state of the assembly is obtained as the reference state. Afterwards, a small strain along in-plane surface directions is applied to the assembly and re-equilibrated for relaxation. The surface energy density can be computed for this corresponding strain. Repeating this procedure with different strains (−1~1%), we obtain a mesh of surface energy density as a function of the strain, and the surface properties are fitted in the end.

As an example, Figure 2 shows the fitted result between the surface energy and shear strain for the Cu(100) surface. In the MS simulation, there are 4 lattice units in each of the two in-plane surface directions and 20 lattice units in the direction vertical to the surface. The shear strain varies from −1~1%, with 40 incremental steps. The fitted curve is obtained as a quadratic polynomial as:

$$\mathsf{\Gamma }=1.2879+1.3965\times {10}^{-12}{\epsilon }_{12}-0.2831{\epsilon }_{12}^2,$$

(5)

where the constant surface energy density is obtained as 1.2879 J/m², the surface elastic constant is obtained as −0.2831 J/m², and the surface stress is very small that can be regarded as zero. Furthermore, the fitting R-square statistic is 0.9999991 and the estimate of error variance is $1.86\times {10}^{-17}$. For details, the readers could refer to the LAMMPS input file and MATLAB post-process file included in the Supplementary Materials for reproducibility.

In this methodology, surface relaxation is fully taken into consideration, which is corresponded to the correct physical situation [39], compared with unrelaxed surfaces. However, surface reconstruction is out of the scope, due to the complexity when surface atoms rearrange to form different structural configurations. Furthermore, there could be two types of the surface with the same orientation (e.g., the wurzite GaN(0001) surface could be terminated by Ga or N atoms), however, only the result of the surface with the lower surface energy is reported, since it is the more stable and common case in reality.

3. Results

In this work, FCC pure metals (Cu, Ni, Pd and Ag), BCC pure metals (Mo and W), diamond Si, zincblende GaAs and GaN are taken as typical examples of cubic crystals, while HCP pure metals (Mg, Zr and Ti), and wurzite GaN are taken as typical examples of hexagonal crystals. Both bulk and surface properties of cubic and hexagonal crystals are calculated through MS simulations in LAMMPS [41] (https://lammps.sandia.gov/). All interatomic potentials are retrieved from NIST Interatomic Potentials Repository (https://www.ctcms.nist.gov/potentials/testing/). It should be mentioned that the methodology in Section 2 is in general, but the results from MS calculations could vary on the choice of different interatomic potentials. EAM and FS potentials are adopted for pure metals, since they are widely used in atomistic simulations of metallic materials and have produced many reasonable surface-related properties of metals [42], such as stacking fault energies, vacancy formation and migration energies, surface energies, surface relaxations and reconstructions. On the other hand, interatomic potentials for semiconductor materials are less mature and they have still been evolving in recent years, and continue to evolve. Here, we choose the Stillinger–Weber potential for Si and Tersoff potential for GaAs and GaN, as they are the most widely used in literature nowadays. In case the readers would like to reproduce any result in this work with the same interatomic potential, each result from this work is listed under the potential file name, respectively, which is identical to that from NIST Interatomic Potentials Repository. In the following, as the prime concern, each result from this work is compared with that calculated using the same potential from the original paper, otherwise the referenced result is taken from different sources, e.g., FP calculations. However, the methodology to obtain bulk and surface properties would be different from Equations (2) and (4), because the energy cannot be decomposed into atomic contributions in FP calculations, where only the total energy of the assembly of atoms is available. Specifically, Equations (2) and (4) are limited to classical molecular simulations, but the original definitions of bulk and surface properties in Equations (1) and (3) are still valid, regardless of the type of atomistic simulations. In that case, it is still possible to obtain the properties through FP calculations, and the readers could refer to the literature to first compute the energies through FP calculations and then obtain the properties by considering Equations (1) and (3).

3.1. Bulk Elastic Constants

The bulk elastic constants of cubic and hexagonal crystals are calculated and listed in

Table 1. Bulk elastic constants of cubic crystals (unit: GPa).

Crystals	Potentials	C11	C12	C44
Cu	Ref. [15]	167	124	76
	Cu_u3.eam	167.3	124.2	76.4
Ni	Ref. [15]	233	154	128
	Ni_u3.eam	233.3	154.3	127.6
Pd	Ref. [15]	218	184	65
	Pd_u3.eam	218.2	184.5	64.9
Ag	Ref. [15]	129	91	57
	Ag_u3.eam	129.1	91.0	56.8
Mo	Ref. [16]	464.7	161.5	108.9
	mo.fs.eam.alloy	465.9	161.7	108.9
W	Ref. [16]	522.4	204.4	160.6
	W.eam.fs	534.3	210.2	164.9
Si	Ref. [43]	151.4	76.4	56.4
	Si.sw	151.6	76.5	56.5
GaAs	Ref. [44]	123.6	48.2	39.4
	GaAs.tersoff	123.9	48.3	39.2
GaN(zincblende)	Ref. [45]	287	169	128
	GaN.tersoff	287.2	169.1	128.5

and

Table 2. Bulk elastic constants of hexagonal crystals (unit: GPa).

Crystals	Potentials	C11	C33	C12	C13	C66
Mg	Ref. [46]	69.6	69.5	25.3	16	12.8
	Mg_mm.eam.fs	69.1	69.4	26.0	16.1	12.8
Zr	Ref. [22]	150	175	85	67	36.0
	Zr.eam.fs	150.9	177.0	84.5	70.3	35.8
Ti	Ref. [21]	180.0	217.1	87.4	76.6	51.4
	Ti_v2.eam.fs	179.0	216.8	85.5	74.7	51.3
GaN(wurzite)	Ref. [45]	347	381	154	123	81
	GaN.tersoff	345.4	380.4	158.1	123.2	82.2

, respectively. If no result calculated using the same potential is found in the original paper, the crystal symbol will be denoted by an underline, to indicate that the referenced result is taken from different sources. Specifically, in Table 1, the referenced results of BCC pure metals of Mo and W are taken from an earlier publication [16] of the same group of authors, because the original paper [20] with the potentials did not present any results of bulk elastic constants. The referenced result of Si is taken from a third-party paper [43], using the same potential which was obtained through lattice dynamics method instead of MS calculations. As expected, excellent agreements between this work and the literature are found for the bulk elastic constants of all crystals. In fact, the results in Table 1 and Table 2 are also very close to those measured by experiments, because bulk elastic constants are usually important parameters that are fitted to obtain the potentials.

3.2. Surface Properties

Through a semi-analytical method, Dingreville and Qu [40] calculated all surface properties of some FCC crystals in the Lagrangian frame with the same potential in this work, so their results and those from this work are first compared in

Table 3. Surface properties of face-centered-cubic (FCC) crystals (unit: J/m²).

	Cu			Ni			Pd			Ag
	Ref. [15]	Ref. [40]	Cu_u3.eam	Ref. [15]	Ref. [40]	Ni_u3.eam	Ref. [15]	Ref. [40]	Pd_u3.eam	Ref. [15]	Ref. [40]	Ag_u3.eam
(100)
${\mathsf{\Gamma }}^{\left(0\right)}$	1.28	1.288	1.288	1.58	1.572	1.570	1.37	1.377	1.364	0.705	0.703	0.702
${\mathsf{\Gamma }}_{11}^{\left(1\right)}$		1.396	1.396		1.321	1.319		1.981	2.000		0.816	0.815
${\mathsf{\Gamma }}_{1111}^{\left(2\right)}$		−0.712	−0.748		−0.865	−0.771		−2.36	−2.706		−1.245	−1.353
${\mathsf{\Gamma }}_{1122}^{\left(2\right)}$		5.914	5.885		10.722	10.600		2.611	2.089		3.343	3.331
${\mathsf{\Gamma }}_{1212}^{\left(2\right)}$		−0.992	−0.283		−0.927	−0.282		−3.25	−2.366		−1.666	−1.260
(110)
${\mathsf{\Gamma }}^{\left(0\right)}$	1.4	1.413	1.412	1.73	1.721	1.719	1.49	1.482	1.478	0.77	0.768	0.764
${\mathsf{\Gamma }}_{11}^{\left(1\right)}$		1.126	1.125		1.054	1.062		1.23	1.337		0.492	0.503
${\mathsf{\Gamma }}_{22}^{\left(1\right)}$		0.993	0.992		0.706	0.710		1.656	1.771		0.684	0.690
${\mathsf{\Gamma }}_{1111}^{\left(2\right)}$		−7.798	−10.656		−13.031	−17.581		−4.775	−7.561		−5.51	−7.741
${\mathsf{\Gamma }}_{2222}^{\left(2\right)}$		−2.263	−3.933		0.95	−2.288		−6.654	−8.697		−2.246	−3.923
${\mathsf{\Gamma }}_{1122}^{\left(2\right)}$		−3.6	−5.694		−5.045	−8.872		−2.086	−4.461		−2.332	−4.226
${\mathsf{\Gamma }}_{1212}^{\left(2\right)}$		−4.436	−3.980		−7.827	−7.604		−3.378	−2.552		−3.296	−2.965
(111)
${\mathsf{\Gamma }}^{\left(0\right)}$	1.17	1.181	1.181	1.45	1.436	1.436	1.22	1.224	1.215	0.62	0.62	0.617
${\mathsf{\Gamma }}_{11}^{\left(1\right)}$		0.866	0.863		0.457	0.451		1.848	1.899		0.636	0.633
${\mathsf{\Gamma }}_{1111}^{\left(2\right)}$		2.054	1.922		6.526	6.281		−2.914	−3.531		0.888	0.749
${\mathsf{\Gamma }}_{1122}^{\left(2\right)}$		1.086	1.133		3.986	3.857		−1.014	−0.915		1.194	1.116
${\mathsf{\Gamma }}_{1212}^{\left(2\right)}$		−1.071	0.798		−1.188	1.450		−2.354	−0.367		1.173	0.119

to validate the methodology. The results of surface energy density from the original potential paper [15] are also listed for reference. It can be seen from Table 3 that all the results of surface energy density are in excellent agreements. Furthermore, all the results of surface stresses from this work and Ref. [40] also agree very well with each other. On the contrary, some discrepancies can be seen in the results of surface stiffness from this work and Ref. [40], although most of them are still in fair agreement. Particularly, the constant, ${\mathsf{\Gamma }}_{1212}^{\left(2\right)}$, in the surface stiffness of Cu(111) and Ni(111) surfaces, and the constant, ${\mathsf{\Gamma }}_{2222}^{\left(2\right)}$, in surface stiffness of Ni(110) surface are calculated to have different signs between this work and Ref. [40]. The sound agreements of the results of all surface energy density, all surface stresses and many of the surface stiffness values guarantee the correctness of the calculation. In case that the situation demands rough surface stiffness values, one could use the data in Table 3 for reference. Otherwise, if one requires accurate surface stiffness values and finds conflicting values in Table 3, then one is encouraged to recalculate it through FP calculations.

The surface properties of BCC, diamond and zincblende crystals are listed in

Table 4. Surface properties of body-centered-cubic (BCC), diamond and zincblende crystals (unit: J/m²).

	Mo		W		Si		GaAs			GaN(zincblende)
	Ref. [19]	mo.fs.eam.alloy	Ref. [19]	W.eam.fs	Ref. [47]	Si.sw	Ref. [47]	Ref. [28]	GaAs.tersoff	GaN.tersoff
(100)
${\mathsf{\Gamma }}^{\left(0\right)}$	2.1	2.100	2.924	2.975	2.53	2.356	2.2		1.429	3.380
${\mathsf{\Gamma }}_{11}^{\left(1\right)}$	2.241	2.240	3.032	3.393		0.002			0.524	−0.636
${\mathsf{\Gamma }}_{1111}^{\left(2\right)}$		−18.862		−13.794		−11.423			−11.465	−13.906
${\mathsf{\Gamma }}_{1122}^{\left(2\right)}$		−16.093		−14.505		−1.242			−1.416	−0.592
${\mathsf{\Gamma }}_{1212}^{\left(2\right)}$		−2.168		−6.993		−3.832			−2.153	−7.436
(110)
${\mathsf{\Gamma }}^{\left(0\right)}$	1.829	1.829	2.575	2.608	1.78	1.666	1.5	0.91	0.924	2.307
${\mathsf{\Gamma }}_{11}^{\left(1\right)}$	2.019	2.019	2.385	2.489		0.005			0.521	−0.852
${\mathsf{\Gamma }}_{22}^{\left(1\right)}$	0.775	0.774	0.271	0.258		0.003			0.988	−1.527
${\mathsf{\Gamma }}_{1111}^{\left(2\right)}$		−3.568		−0.257		−16.350			−15.639	−22.743
${\mathsf{\Gamma }}_{2222}^{\left(2\right)}$		−9.262		−9.778		−6.163			−4.037	−6.115
${\mathsf{\Gamma }}_{1122}^{\left(2\right)}$		7.617		9.342		−8.444			−7.475	−13.152
${\mathsf{\Gamma }}_{1212}^{\left(2\right)}$		−0.707		0.476		−4.699			−2.067	−10.133
(111)
${\mathsf{\Gamma }}^{\left(0\right)}$		2.356		3.337	1.46	1.361	1.3		0.751	1.864
${\mathsf{\Gamma }}_{11}^{\left(1\right)}$		0.801		0.392		0.001	0.5		0.988	−1.502
${\mathsf{\Gamma }}_{1111}^{\left(2\right)}$		−20.118		−25.246		−3.829			−3.058	1.139
${\mathsf{\Gamma }}_{1122}^{\left(2\right)}$		6.041		3.730		−3.836			−5.056	−1.472
${\mathsf{\Gamma }}_{1212}^{\left(2\right)}$		−12.560		−14.568		−7.574			1.486	0.572

and those of hexagonal crystals are listed in

Table 5. Surface properties of hexagonal crystals (unit: J/m²).

	Mg		Zr		Ti		GaN(wurzite)
	Ref. [48]	Mg_mm.eam.fs	Ref. [22]	Zr mm.eam.fs	Ref. [21]	Ti_v2.eam.fs	Ref. [30]	GaN.tersoff
(0001)
${\mathsf{\Gamma }}^{\left(0\right)}$	0.238	0.238	1.022	1.022	0.993	0.998	2.0027	1.865
${\mathsf{\Gamma }}_{11}^{\left(1\right)}$		0.289	1.161	1.140	0.96	0.959		−1.502
${\mathsf{\Gamma }}_{1111}^{\left(2\right)}$		1.041		−0.347		−0.191		0.791
${\mathsf{\Gamma }}_{1122}^{\left(2\right)}$		0.947		0.502		0.616		−1.094
${\mathsf{\Gamma }}_{1212}^{\left(2\right)}$		0.148		0.231		0.138		0.247
$(10\overline{1}0)$
${\mathsf{\Gamma }}^{\left(0\right)}$		0.263	1.086	1.086	1.061	1.061	1.7624	1.995
${\mathsf{\Gamma }}_{11}^{\left(1\right)}$		0.063	0.655	0.504	0.481	0.477		−1.348
${\mathsf{\Gamma }}_{22}^{\left(1\right)}$		0.168	0.476	0.666	0.68	0.648		−1.354
${\mathsf{\Gamma }}_{1111}^{\left(2\right)}$		−0.359		−9.039		−2.991		−1.003
${\mathsf{\Gamma }}_{2222}^{\left(2\right)}$		−2.158		−8.234		−10.369		−18.547
${\mathsf{\Gamma }}_{1122}^{\left(2\right)}$		0.791		−5.305		−2.913		−4.437
${\mathsf{\Gamma }}_{1212}^{\left(2\right)}$		−3.594		1.214		−0.055		−12.293
$(11\overline{2}0)$
${\mathsf{\Gamma }}^{\left(0\right)}$		0.343	1.23	1.232	1.187	1.188	1.9707	2.315
${\mathsf{\Gamma }}_{11}^{\left(1\right)}$		0.123	0.844	0.905	0.915	0.916		−0.906
${\mathsf{\Gamma }}_{22}^{\left(1\right)}$		0.317	0.989	1.010	1.019	1.018		−1.303
${\mathsf{\Gamma }}_{1111}^{\left(2\right)}$		−7.452		−16.977		−16.299		−27.498
${\mathsf{\Gamma }}_{2222}^{\left(2\right)}$		−3.681		−9.752		−12.866		−22.276
${\mathsf{\Gamma }}_{1122}^{\left(2\right)}$		−1.649		−7.660		−7.211		−10.043
${\mathsf{\Gamma }}_{1212}^{\left(2\right)}$		−0.194		−2.688		−3.642		−5.771

, respectively. By comparing these with the available results in literature, it can be seen that the surface energy density and surface stress of BCC pure metals of Mo and W in Table 4 and the HCP pure metals of Mg, Zr and Ti in Table 5, are in excellent agreement, since they are obtained using the same potential through MS calculations. For semiconductor materials, the surface energy density of the three wurzite GaN surfaces are taken from the FP calculations [30] in Table 5. There is no report about those of zincblende GaN surfaces, because GaN usually crystallizes in the wurtzite lattice, while the much less common zincblende GaN could only be grown under certain conditions. In Table 4, the surface energy density of the three Si and GaAs surfaces, as well as the surface stress of the GaAs(111) surface, are taken from Ref. [47], which were obtained through analytical calculations with force constants and surface bonding model. There are some results from FP calculations in literature [23,24,25,26,27,29] about the surface energy density of Si and GaAs surfaces, but they were given in units of eV per cell or per surface atom, and only the surface energy density of the GaAs(110) surface [28] was mentioned in unit of eV per unit area for our comparison. However, although the referenced results of semiconductor materials are taken from different resources, good agreements with those of this work are found in Table 4 and Table 5, which demonstrate the validity of the methodology in this work. In this case, the readers could follow the same methodology to obtain the property of other materials, according to their own interest and with their own choice of interatomic potentials.

4. Discussions

It is of interest to make comparative studies of the surface properties to uncover trends in their behaviors. For the surface energy density, the results of all the crystals in this work are positive. In FCC pure metals, the surface energy density is ordered as (110) > (100) > (111). This is as expected, since the (111) plane is the most densely packed plane for FCC crystals, while the (110) plane is the least densely packed. Similar trends can be found for other crystals following the packing density of each plane. In BCC pure metals, it is ordered as (111) > (100) > (110). Moreover, the two different elements in the zincblende structure are arranged as two diamond structures individually, with some offset; therefore, in diamond Si, zincblende GaAs and GaN, all the surface energy density is ordered as (100) > (110) > (111). Similarly, the two different elements in the wurzite structure are arranged as two HCP structures individually with some offset; therefore, in HCP pure metals and wurzite GaN, it is ordered as $(11\overline{2}0)$ > $(10\overline{1}0)$ > (0001). Furthermore, the surface energy of different planes in a crystal could be used to determine the equilibrium morphology and equilibrium growth rate in the direction normal to each plane, e.g., the surface energy data could be fed into Wulff construction, to determine the shape and exposed surfaces of nanoparticles [49].

For the surface stress, the results of almost all crystals in this work are positive, but those of zincblende and wurzite GaN are all negative. In this case, the absolute values of GaN are used for comparison in the following. Since the surface stresses in cubic (110), hexagonal $(10\overline{1}0)$ and $(11\overline{2}0)$ surfaces are different along the two in-plane directions with ${\mathsf{\Gamma }}_{11}^{\left(1\right)}\ne {\mathsf{\Gamma }}_{22}^{\left(1\right)}$, averaged surface stress is used for comparison. In BCC pure metals of Mo and W, it is ordered as (100) > (110) > (111), and the same trend is found in other BCC pure metals of V, Nb, Ta [19]. On the contrary, in FCC pure metals of Cu and Ni, it is ordered as (100) > (110) > (111), but different ordering of (100) > (111) > (110) is found in Pd and Ag. In diamond Si, zincblende GaAs and GaN, it is ordered as (111) > (110) > (100). On the other hand, the surface stress in HCP pure metals is ordered as (0001) > $(11\overline{2}0)$ > $(10\overline{1}0)$, but different ordering of (0001) > $(10\overline{1}0)$ > $(11\overline{2}0)$ is found in wurzite GaN.

For the surface stiffness, the results are found to vary with different signs in different surfaces for all crystals, so no clear trend can be drawn for the surface stiffness. As pointed out by Shenoy [39], unlike the positive definiteness of the bulk stiffness tensor to guarantee the stability of the crystal, the surface stiffness tensor does not need to be positive definite, because a surface cannot exist independent of the bulk, so the stability of the crystal is maintained as long as the total energy (bulk + surface) satisfies the positive definiteness condition. Furthermore, it should be noted that all the results of surface energy density and surface stress are on the order of several J/m², but the results of the surface stiffness could reach tens of J/m². This indicates that surface stiffness could play an important role in changing the constitutive law seen in the classical elasticity theory, and it should not be neglected for the investigation of nano-sized structural elements.

5. Conclusions

This work investigates surface properties of some typical cubic and hexagonal crystals through MS calculations. FCC pure metals (Cu, Ni, Pd and Ag), BCC pure metals (Mo and W), diamond Si, zincblende GaAs and GaN are taken as typical examples of cubic crystals, while HCP pure metals (Mg, Zr and Ti), and wurzite GaN are taken as typical examples of hexagonal crystals. Bulk elastic constants are first obtained for all crystals to validate the methodology. Certain surface properties of some crystals in this work are also compared with available results in literature, to guarantee the correctness of the methodology. Sound agreements of the bulk and surface properties between this work and the literature are found. New results are presented for the surface stiffness of BCC pure metals, surface stress and surface stiffness of HCP pure metals, Si and compound crystals. The surface energy densities are calculated to be positive for all crystals, and the surface stresses are positive for almost all crystals, but those of zincblende and wurzite GaN are negative. On the contrary, the surface stiffnesses are found to vary, with different signs in different surfaces for all crystals. Furthermore, the results of surface energy density and surface stress are on the order of several J/m², but the results of the surface stiffness could reach tens of J/m². Moreover, comparative studies of surface energy density, surface stress and surface stiffness are carried out, to uncover trends in their behaviors.