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

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.


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 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.

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: where U 0 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 11 , C 12 and C 44 , while hexagonal crystals have 5 independent elastic constants (take axis-3 as the axis of symmetry): C 11 , C 12 , C 13 , C 33 and C 66 .
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: where E i is the energy of the i-th atom and V 0 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.

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]: where Γ (0) is the surface energy density at the relaxed and unstrained state of the surface. Γ (1) is the surface stress with Γ βα that exists when the surface strain is absent. Γ (2) is the surface stiffness tensor with Γ γδαβ . 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, Γ (0) , Γ (1) and Γ (2) 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 Γ

2
where Γ (0) is the surface energy density at the relaxed and unstrained state of the surface. Γ (1) is the surface stress with Γ ( ) = Γ ( ) that exists when the surface strain is absent. Γ (2) is the surface stiffness tensor with Γ ( ) = Γ ( ) = Γ ( ) = Γ ( ) . 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, Γ (0) , Γ (1) and Γ (2) 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 Γ 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]: 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]: where E i is the energy of the i-th surface atom, E 0 is the energy of any bulk atom deep in the interior of a large crystal, and A 0 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 where the constant surface energy density is obtained as 1.2879 J/m 2 , the surface elastic constant is obtained as −0.2831 J/m 2 , 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 × 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.

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. 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.

Results
In this work, 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).

Bulk Elastic Constants
The bulk elastic constants of cubic and hexagonal crystals are calculated and listed in Tables 1 and 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 Tables 1  and 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.

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 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, Γ (2) 1212 , in the surface stiffness of Cu(111) and Ni(111) surfaces, and the constant, Γ (2) 2222 , 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 and those of hexagonal crystals are listed in Table 5, 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 Tables 4 and 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. Table 3. Surface properties of face-centered-cubic (FCC) crystals (unit: J/m 2 ).

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 (1120) > (1010) > (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 (1010) and (1120) surfaces are different along the two in-plane directions with Γ (1) 11 Γ (1) 22 , 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) > (1120) > (1010), but different ordering of (0001) > (1010) > (1120) 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 2 , but the results of the surface stiffness could reach tens of J/m 2 . 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.

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 2 , but the results of the surface stiffness could reach tens of J/m 2 . Moreover, comparative studies of surface energy density, surface stress and surface stiffness are carried out, to uncover trends in their behaviors.