This notebook demonstrates how to generate various slabs for a material using In addition, it will also demonstrate how to obtain the Wulff shape with a set of surface energies associated with Miller indices.


R. Tran, Z. Xu, B. Radhakrishnan, D. Winston, W. Sun, K. A. Persson, S. P. Ong, "Surface Energies of Elemental Crystals", Scientific Data, 2016, 3:160080, doi:10.1038/sdata.2016.80.

Sun, W.; Ceder, G. Efficient creation and convergence of surface slabs, Surface Science, 2013, 617, 53–59, doi:10.1016/j.susc.2013.05.016.

Calculating the surface energy

To do this, we actually need to calculate (from first principles) the total energy of two structures. The total energy of the bulk ($E_{bulk}$) (preferably the oriented unit cell corresponding to the slab) and the total energy of the slab($E_{slab}$) which is several layers of atoms interfaced with a vacuum. Using the following equation, we can get the surface energy $\gamma = \frac{E_{slab}-E_{bulk}}{2A}$ where both slab and bulk systems have the same number of atoms and species (you can multiply $E_{bulk}$ where n is the number of oriented unit cell layers contained in your slab system relative ot your bulk system). The factor of 1/2 accounts for the two surfaces in our slab model so $\gamma$ is the surface energy of one surface.

Some basic properties of the Wulff shape

Using the WulffShape class, there are a couple of properties we can get. Here we will explain what some of them are.


Surface Gibbs free energy for a crystal is given by $\Delta G=\sum\limits_{hkl}\gamma_{hkl}A_{hkl}$. Where $\gamma_{hkl}$ is the surface energy of facet (hkl) and $A_{hkl}$ is the surface area of that particular facet that occupies the Wulff shape. We can normalize this value with the total surface area of the Wulff shape to get the weighted (average) surface energy for a particular material $\bar{\gamma}=\frac{\Delta G}{\sum\limits_{hkl}A_{hkl}}$


Typically in the literature when discussing surface anisotropy, we would only look at the ratios of 2 surface energies when talking about anisotropy. eg. the ratio of a generic fcc (111) to (100) surface energy should be less than 1 as the (111) facet is the closest packed surface of an fcc structure and should have the lowest surface energy. However this method of determining surface anisotropy does not allow us to determine an overall anisotropy of a material, ie. how different are all the surface energies for a material. As such, we used the Coefficient of Variation from the weighted surface energy. For reference, an ideal sphere Wulff shape (eg. completely isotropic) has a anisotropy of 0.


An alternative to anisotropy. This is useful for determining the critical nucleus size. A large shape factor indicates great anisotropy. See Ballufi, R. W., Allen, S. M. & Carter, W. C. Kinetics of Materials. (John Wiley & Sons, 2005), p.461

When generating a slab of LiFePO4, we also want to be careful not break any of the P-O bonds. These bonds are strong enough that they will result in a significantly high surface energy when broken so its reasonable to say that any terminations with such broken bonds will not yield the lowest surface energy. To implement this, we add the bonds parameter to get_slabs, a dictionary where the key will be two atoms whose bonds we do not want to break and the element of that value would be their maximum bond length.

There are a couple of rules before we actually run calculations on some of these slab models. First off, we need to ensure that all slabs we will be calculating have the same surface on both sides. To do this, we need to ensure the slab model has Laue point group symmetry ie. contains inversion symmetry. We use the is_symmetric() property of our slab object to check this. It's important that both surfaces are the same as the above equation for surface energy is used to get the energy of one surface, hence the 1/2 factor in the equation. If the surfaces are different (the slab is not symmetric), we would be calculating the average surface energy of two different surfaces in our slab rather than the surface energy for one slab in our calculation. Secondly, for structures containing oxidation states, we need to ensure that our surfaces are nonpolar. A polar termination will lead to a very high surface energy, so we can skip those particular structures. We can check polarity using the is_polar() property of our slab object. Both these criterias (nonpolar and symmetric) should be satisfied before calculating a particular slab model.