Surface energies of elemental crystals

1The surface properties of a crystal are crucial to the understanding and design of materials for many applications. For instance, technologies such as fuel cells and industrial chemical manufacturing require the use of catalysts to accelerate chemical reactions, which is fundamentally a surface-driven process. Surface effects are especially important in nanomaterials, where relatively large surface area to volume ratios lead to properties that differ significantly from the bulk material. For example, the nanoscale stability of metastable polymorphs is determined from the competition between surface and bulk energy of the nanoparticle.

2The stability of a surface is described by its surface energy γ, a measure of the excess energy of surface atoms due to a variety of factors, such as the broken bonds yielding undercoordinated atoms. This fundamental quantity is important in understanding surface structure, reconstruction, roughening and the crystal’s equilibrium shape. Despite its importance, experimental determination of surface energies, especially for specific facets, is difficult and rare. Furthermore, experimentally observed Wulff shapes are often inconsistent across studies due to the sensitivity of high energy facets to temperature and impurities. References and have accumulated a large set of metallic elemental surface energy data by extrapolating surface tension of liquid phases for solid surfaces. Reviews of such surface tension techniques have been compiled by Mills et al. and Keene.

3First principles computations such as those based on density functional theory (DFT) methods are important complementary tools to experimental techniques in characterizing surface properties of a material. Computational techniques provide the means to precisely control the surface structure and composition. Indeed, fundamental and application-driven computational studies of surfaces in the literature are extensive,. On the broader scale, Vitos et al. have previously compiled a database of surface energies for all metals up to Pu using the full charge density (FCD) DFT method, a technique based on the coupling of the linear muffin-tin orbital method and the atomic-sphere approximation. However, this database is limited to surfaces of ground state crystals up to a maximum Miller index (MMI) of 1 only.

4The challenges for DFT determination of surface properties are three-fold. First, the choice of the exchange correlation functional as well as other parameters such as pseudopotentials, integration grid, etc. can significantly affect the accuracy and convergence of surface energies, which in turn leads to values that are generally difficult to compare across different works. Second, the typical ‘slab’ approach for performing surface energy calculations requires the use of large supercells with the introduction of a large vacuum region, which makes such calculations computationally intensive, especially for low symmetry materials and high Miller indices. Finally, some surfaces undergo significant reconstruction to reduce surface energy and increase stability, the most notorious example being the 7×7 reconstruction of the Si (111) surface. In order for the DFT geometry optimization algorithms to identify the global energy minimum surface structure, a reasonable initial guess for the reconstructed surface is needed to avoid relaxation to a local minimum, and such initial guesses are difficult to determine a priori without experimental input.

5In this work, we present a large database of surface energies and Wulff shapes of more than 100 polymorphs of about 70 elements constructed using high-throughput DFT calculations. The elemental crystals studied include all polymorphs available in the Materials Project (https://www.materialsproject.org) and the chemistries covered include both metals and non-metals. All surfaces up to a maximum absolute Miller index (MMI) of 2 and 3 were computed for non-cubic and cubic crystals respectively. To address the first and second challenges outlined above, a robust computational workflow based on the efficient convergence scheme proposed by Sun et al. was developed to automate the large number of multi-step calculations. To partially address the third challenge, we include well-known reconstructions such as the missing row (110) and hexagonal (111) reconstructions for the face-centered-cubic surfaces, in our data set. All data have been made publicly available via the Crystalium web application at the Materials Virtual Lab website (http://crystalium.materialsvirtuallab.org), as well as the Materials Project’s graphical user interface and RESTful Application Programming Interface (API),. The full dataset is also available as a JSON file at the Dryad-repository (Data Citation 1). We show that the calculated surface energies are in excellent agreement with available experimental values. We will also provide an assessment of the effect of surface reconstructions on the accuracy of the data.

6The approach used to compute surface energies in this work is the typical slab model, wherein a supercell of the crystal oriented to expose the relevant surface of interest is generated, and atoms are removed from a portion of the supercell to create a vacuum. An example is given in Fig. 1. Starting from the conventional unit cell, a series of lattice vector transformations is performed to create an ‘oriented’ unit cell (OUC) where the a and b lattice vectors are parallel to the plane with Miller indices (hkl) and the c lattice vector is not parallel to the plane. It should be noted that c vector is not necessarily perpendicular to the plane, though the ‘most orthogonal’ vector that can be obtained within a reasonable cell size is used. The OUC is then extended by multiples of the c lattice vector and atoms are removed to generate the slab model with a desired slab and vacuum thickness. To obtain all symmetrically distinct terminations for a given plane, atoms in the OUC are shifted in the c direction prior to extension to generate the slab model. All slabs generated are constrained to have symmetric top and bottom surfaces. The model generation algorithm has been implemented in the Python Materials Genomics (pymatgen) materials analysis library, with a comprehensive set of unit tests to ensure robust functioning of the code.

7For a given slab model of a facet with Miller index (hkl), the surface energy is given by the following expression: where, is the total energy of the slab model with termination σ, is the energy per atom of the bulk OUC, nslab is the total number of atoms in the slab structure, Aslab is the surface area of the slab structure, and the factor of 2 in the denominator accounts for the two surfaces in the slab model. To efficiently converge the surface energy, both the bulk and slab energies were calculated using the OUC, which allows consistent reciprocal integration grids to be used for the bulk and slab calculations,. The OUC was fully relaxed in both cell volume and atomic positions, while the slabs were relaxed in the atomic positions only.

8All DFT, energy calculations were performed using the Vienna Ab initio Simulation Package (VASP) within the projector augmented wave (PAW) approach. The exchange-correlation effects were modeled using the Perdew-Berke-Ernzerhof (PBE) generalized gradient approximation (GGA) functional, and all calculations were spin-polarized with a plane wave cutoff energy of 400 eV. The pseudopotentials used are similar to those used in the Materials Project. The energies and atomic forces of all calculations were converged to within 10−6 eV and 0.01 eV Å−1 respectively. The Methfessel-Paxton method was chosen as the smearing algorithm, the blocked Davidson iteration scheme was chosen as the electron minimization algorithm, and ions were updated using the conjugated gradient algorithm. Γ-centered k-point meshes of and were used for OUC and slab calculations respectively with non-integer values rounded up to the nearest integer. Through a series of comprehensive convergence tests, it was determined that vacuum and slab thicknesses of at least 10 Å were sufficient to ensure that the surface energies were converged to within 0.02 Jm−2.

9The Wulff construction gives the crystal shape under equilibrium conditions,,. In this construction, the distance of a facet from the crystal center is proportional to the surface energy of that facet, and the inner convex hull of all facets form the Wulff shape. Here, we use the methodology developed by Miracle-Sole to construct the Wulff shape. An example of the Wulff shape of Fe is given in Fig. 2a, which contains all surfaces up to a MMI of 3. To illustrate the importance of sampling higher Miller index surfaces in obtaining an accurate Wulff shape, the Wulff shape of Fe up to MMI of 1 is shown in Fig. 2b. It is evident that inclusion of higher-index surfaces is necessary to produce an accurate description of the equilibrium crystal shape.

10The area fraction of each symmetrically distinct facet under equilibrium conditions can be determined from the Wulff shape. We define the weighted surface energy using this fraction as given by the following equation: where γhkl is the surface energy of a unique facet existing in the Wulff shape, Ahkl is the total area of all facets in the {hkl} family in the Wulff shape, and is the area fraction of the {hkl} family in the Wulff shape. In this work, is used as a basis for comparison to experimentally determined surface energies, for which usually only one value is reported.

11Several measures of surface anisotropy have been proposed in the literature. For instance, the ratio γ111/γ110 is frequently used for fcc and bcc metals, but such surface-specific measures lack general applicability across crystals of different symmetry. The most commonly used general measure of anisotropy is the shape factor η, which is given by the following equation: where A and V are the surface area and volume of the Wulff shape, respectively. The shape factor is a useful quantity in determining the critical nucleus size. Typically, the shape factor is compared against that of an ideal sphere, and a larger η indicates greater anisotropy.

12Though generally applicable, η does not account for variation in surface energies except indirectly through its impact on the Wulff shape. In the database, we have provided an alternative definition of surface energy anisotropy αγ given by the following equation:

13αγ can effectively be viewed as a coefficient of variation of surface energies that is normalized for comparison across crystals with different weighted surface energies. A perfectly isotropic crystal would have αγ=0. αγ is comparable across different crystal systems and accounts for all surfaces based on their relative importance (in terms of contribution to the Wulff shape).

14The current release of the database contains surface energies and Wulff shapes of more than 100 polymorphs of about 70 elements. This database is far larger in scope than previous compilations of surface energies in several important ways:

15Both metals and non-metals are included.All polymorphs for each element available in the Materials Project database are included. This is in contrast with previous works where only the most stable ground state crystal structures were studied.All surfaces up to a maximum Miller index (MMI) of two and three were calculated for all non-cubic and cubic crystals respectively. The most common spacegroups among the crystals sampled are,, P63/mmc, and, which contains 13, 13, 12, and 13 symmetrically distinct surfaces respectively up to the MMI used in this work. The only exceptions are the ground state crystals of Mn, S, P, Se and B for which only the terminations and surfaces containing the least number of broken bonds were calculated due to the low symmetry and/or large unit cell sizes.Well-known reconstruction schemes are incorporated. At the time of writing, the schemes considered are the (110) 2×1 missing row and (111) 1×1 hexagonal and (111) 2×1, (110) 1×1 and (100)p 2×1 reconstructions.

16Both metals and non-metals are included.

17All polymorphs for each element available in the Materials Project database are included. This is in contrast with previous works where only the most stable ground state crystal structures were studied.

18All surfaces up to a maximum Miller index (MMI) of two and three were calculated for all non-cubic and cubic crystals respectively. The most common spacegroups among the crystals sampled are,, P63/mmc, and, which contains 13, 13, 12, and 13 symmetrically distinct surfaces respectively up to the MMI used in this work. The only exceptions are the ground state crystals of Mn, S, P, Se and B for which only the terminations and surfaces containing the least number of broken bonds were calculated due to the low symmetry and/or large unit cell sizes.

19Well-known reconstruction schemes are incorporated. At the time of writing, the schemes considered are the (110) 2×1 missing row and (111) 1×1 hexagonal and (111) 2×1, (110) 1×1 and (100)p 2×1 reconstructions.

20This database will be systematically improved through continuous updates. Future updates will include higher MMIs and other surface reconstruction schemes such as the well-known fcc(100)-hex surfaces. Further expansion to non-elements, e.g., binaries, ternaries, etc. is also planned at a later stage.

21Figure 3 shows an outline of the high-throughput workflow used in this work implemented using the FireWorks software package. The initial bulk crystals were obtained by querying for all elemental crystals from the Materials Project via the Materials API. All OUCs up to the MMI for each crystal were then automatically generated. For each OUC (corresponding to a specific Miller index (hkl)), a full relaxation was then carried out, and the slab models for distinct terminations were then generated from the fully relaxed OUC and calculated. To handle errors that may arise during calculations, the custodian software package was used as a wrapper around VASP together with a set of robust error handling rules. The results from successful calculations were then automatically inserted into a MongoDB document-based database. The metadata of the DFT calculations and the surface properties extracted (see Tables 1 and 2) were subsequently inserted into the Materials Project database and the Dryad repository.

22Consistency checks were devised as part of the workflow to detect possible errors and anomalous behavior in the surface calculations. Calculations containing such anomalies are tagged with a warning containing a list of their brief descriptions. The nature of these warnings are detailed in Table 3. Entries tagged with such warnings do not automatically indicate that the calculations are invalid. For example, although surface atom displacements on relaxation are typically expected to be less than 5% (ref. ), the relaxation of the (110) surface of diamond Si exceeds this value due to its tendency to reconstruct to achieve stability.

23Pymatgen is the primary materials analysis code used in this work, and the surface construction and input file generation algorithms are implemented in this package. Both pymatgen and the custodian error recovery library are open-source software under the MIT (Massachusetts Institute of Technology) License. The high-throughput workflow was implemented using the FireWorks library, which is freely accessible under a modified GPL (GNU General Public License). All implemented algorithms, including the model generation and Wulff shape calculation, come with a comprehensive set of unit tests. The open-source software codes are also continuously tested via a continuous integration service. The VASP DFT code used is copyrighted by the University of Vienna, Austria and is accessible under a paid license.

24A user-friendly web application, Crystalium, has been developed to allow users to efficiently explore the Wulff shapes, surface energies and slab structures of the polymorphs investigated. This web application is at http://crystalium.materialsvirtuallab.org (a screenshot is given in Fig. 4). In addition, the surface energies and Wulff shapes are also available on the Materials Project website (https://www.materialsproject.org) on the detailed data pages of specific crystals. Two JSON data files are also available in a Dryad repository (Data Citation 1). The first data file contains the complete set of metadata and surface properties for all materials studied so far while the second contains the key VASP input parameters and output data for all slab calculations performed in this study.

25The surface properties for each material is stored as an individual JSON document (Data Citation 1). The material is described by its metadata which contains information such as the spacegroup, formula, a unique Materials Project identifier (mp-id), the energy per atom above the hull (a measure of the relative stability of the bulk crystal) and the polymorph number. A description of the metadata can be found in Table 1. A similar JSON document provides the key VASP input parameters and output data for each slab calculation, which is also uniquely identified using the mp-id of the crystal.

26The JSON document for each entry contains an organized list of sub-entries that describes the properties of each surface in detail. Each sub-entry contains information such as the Miller index, surface energy and the fraction of the Wulff shape’s surface area occupied by this facet. For each Miller index, the lowest surface energy termination, including among different reconstructions investigated where applicable, is provided in each sub-entry. The slab structure used to model the surface is available as a string in the JSON document in the format of a Crystallographic Information File (cif), which can also be downloaded via the Materials Project website and Crystalium web application. In addition, the weighted surface energy (equation (2)), shape factor (equation (3)), and surface anisotropy (equation (4)) are given. Table 2 provides a full description of all properties available in each entry as well as their corresponding JSON key.

27The data was validated through an extensive comparison with surface energies from experiments and other DFT studies in the literature. Due to limitations in the available literature, only the data on ground state phases were compared.

28Experimental determination of surface energy typically involves measuring the liquid surface tension and solid-liquid interfacial energy of the material to estimate the solid surface energy at the melting temperature, which is then extrapolated to 0 K under isotropic approximations. Surface energies for individual crystal facets are rarely available experimentally. Figure 5 compares the weighted surface energies of all crystals (equation (2)) to experimental values in the literature,,. It should be noted that we have adopted the latest experimental values available for comparison, i.e., values were obtained from the 2016 review by Mills et al., followed by Keene, and finally Niessen et al. and Miller and Tyson. A one-factor linear regression line was fitted for the data points. The choice of the one factor fit is motivated by the fact that standard broken bond models show that there is a direct relationship between surface energies and cohesive energies, and previous studies have found no evidence that DFT errors in the cohesive energy scale with the magnitude of the cohesive energy itself.

29We find that the DFT weighted surface energies are in excellent agreement with experimental values, with an average underestimation of only 0.01 J m−2 and a standard error of the estimate (SEE) of 0.27 J m−2. The Pearson correlation coefficient r is 0.966. Crystals with surfaces that are well-known to undergo significant reconstruction tend to have errors in weighted surface energies that are larger than the SEE.

30The differences between the calculated and experimental surface energies can be attributed to three main factors. First, there are uncertainties in the experimental surface energies. The experimental values derived by Miller and Tyson are extrapolations from extreme temperatures beyond the melting point. The surface energy of Ge, Si, Te, and Se were determined at 77, 77, 432 and 313 K respectively while the energy of the (0001) surface of pyrolytic graphite was determined using the sessile drop technique under high temperatures. In addition, the possibility of contamination by surface active elements such as oxygen can lead to lower surface tension values. As a result, the higher value of surface tension measured in experiments is often argued to be the more accurate value,.

31Second, the limitations of the exchange-correlation functional used can also cause discrepancies,. Though the average difference between the computed PBE surface energies in this work and experimental surface energies is very small (0.01 J m−2), there is nevertheless non-negligible differences for specific elements. In addition, this study does not take into consideration the effects of Van der Waals forces in materials such as graphite where it is the dominant interaction between graphene layers in the (0001) direction. Of the 1,000 different surfaces studied, only the (0001) surfaces of the two graphite polymorphs have unphysical negative surface energies, which was also previously observed by Ooi et al..

32Finally, surface reconstructions could also contribute to differences between the computed and experimental values. To understand the effect of reconstructions on surface energies, we compared the surface energy of the relaxed (110) fcc metal surfaces from the database to the reconstructed configuration described by Zhang et al., as shown in Fig. 6. We find that the missing row reconstruction is predicted to be favored (lower in energy) over the unreconstructed surface for Pt, Au and Ir only, in agreement with previous experimental and computational results,. Even for the surfaces that undergo reconstruction, the differences between the reconstructed and unreconstructed surface energies are relatively small (<0.2 J m−2) in metals.

33Exceptional cases of reconstruction with large differences between the relaxed and reconstructed surface energies do exist. Semiconductors such as Si and Ge are known to have significantly lower surface energies in their reconstructed state. The experimental, reconstructed and relaxed surface energies for Si are shown in Table 4. The reconstructed surface energies are much closer to those found experimentally than the unreconstructed surface energies. The reconstructed (111) surface in particular shows the largest energy difference (0.85 J m−2). It should be noted that this HT work did not explore the well-known 7×7 reconstruction, of the Si(111) surface due to the large supercell required.

34Table 5 presents a comparison of the calculated surface energies in this work with those computed in previous works. Unlike the comparison with experimental data, surface-specific comparisons are presented to demonstrate the accuracy of the data. An ‘averaged’ experimental value is provided where available for reference. Despite the fact that the literature values come from studies with widely different parameters, including choice of exchange-correlation functional, pseudopotential choice/full electron calculations, etc., we find that the calculated values in this work differ from these literature values by only 2–13% (refs,). The (111) surface of fcc Pt from Da Silva et al. has the largest (13%) difference. This is largely due to use of pseudopotentials in this study while Da Silva et al. used an all-electron method.

35In general, we find that the calculated surface energies in this work are slightly lower than that in the work of Vitos et al. This observation may be attributed to the atomic sphere approximation used in the FCD calculations by Vitos et al., which is known to hinder relaxation at the surface, thus leading to higher surface energy values. Also, our calculated surface energies are significantly lower than those computed using the local density approximation (LDA) functional in the literature, which may be attributed to the LDA functional’s propensity to overbind compared to the GGA functional used in this work.

36The planes with the highest atomic density per unit area of fcc, bcc and hcp materials are the {111}, {110} and {0001} planes respectively. According to the classic broken bond model, the minimization of broken bonds in these surfaces leads to these surfaces having the lowest γ,. Our results support this empirical model, with a few notable exceptions. For example, the (0001) surfaces for hcp Sc and Y and the (110) surface for Li are not the lowest energy surfaces among the facets investigated. These exceptions are also observed in Vitos et al. for various materials and in other previous first-principle studies for the Li (110) surface,,.

37The database in this study is the most extensive collection of calculated surface energies for elemental crystalline solids to-date. When used with data mining and machine learning techniques, the database can be used to reveal trends in surface phenomena and guide the screening of potentially interesting materials for target surface properties. For instance, the relaxed surface structures and energies of many metals, particularly the noble metals, will be highly useful in the study of surface absorption of molecules, which is of great fundamental relevance in catalysis. Consideration of surfaces is especially crucial in nanomaterials design, where surface effects tend to dominate the overall performance and properties. We also anticipate the data presented to be a useful starting point for the study of the interfaces between materials, either within the same material (e.g., grain boundaries) or between different materials. In the near future, we will provide the facility for Materials Project users to upload experimentally proposed reconstructions to further improve the completeness of the database. Furthermore, future enhancements would include the surface properties of not just elements, but also multicomponent compounds (including alloys).

38How to cite this article: Tran, R. et al. Surface energies of elemental crystals. Sci. Data 3:160080 doi: 10.1038/sdata.2016.80 (2016).