GULP

Home

Downloads

Help

   FAQ

   Examples

   Running GULP

   GULP manuals

Overview of GULP capabilities

Potential Models

GULP News

GULP Bugs

Calculation of surface properties

The properties of the surfaces of materials are every bit as important as the bulk properties, since they control the interaction between the substance and the external environment. At the most obvious level, the very shape of the particles or crystallites formed is determined by the properties of the surface relative to the bulk, while catalysis and reactions of the material also predominantly occur at the surface. From the earlier discussion of electrostatics, it is apparent that the surface structure is a key factor in determining the bulk polarisation and net dipole of the material, which has consequences even for the bulk region.

Due to the interest in surface phenomena, the history of modelling surfaces using atomistic simulation is also a long one spanning a number of different computer codes. Of the recent era, the dominant computer code for surface simulation of ionic materials has been MARVIN [90], which has been developed alongside GULP for many years. However, to ensure the maximum integration of functionality between bulk and surface simulations, both in terms of forcefield models, as well as accessible properties, it was decided to incorporate much of the functionality of MARVIN into GULP to produce a single code capable of simulating systems with any type of boundary condition from 0-D, through to 3-D. Here we describe the methodology employed for surface calculations.

There are two phases to any surface calculation, namely the creation of the surface from the bulk material and the subsequent calculation of its optimised structure and properties. Each surface is specified by at least two pieces of information. Firstly, there are the Miller indices $\left(hkl\right)$ of the plane that defines the orientation of the bulk cleavage. Secondly, there is the so-called shift - i.e. the displacement of the plane relative to the unit cell origin. For simple cases, such as the $\left(001\right)$ surface of a rock salt structure there is only one unique choice of shift. However, for more complex cases there may be several shifts for a given plane that lead to distinct surfaces. When cleaving surfaces there are also other important considerations to take into account, in particular the type of surface. As illustrated in Figure 1.1, there are three basic types of surface. In type 1, the atomic structure consists of charge neutral sheets of ions parallel to the surface plane and thus all shifts are guaranteed to yield a surface with no dipole moment in the direction of the surface normal. For type 2 surfaces, there are combinations of layers of cations and anions that possess zero net dipole in the appropriate direction. Hence for well chosen values of the shift a non-dipolar surface can be obtained. Finally, for type 3, all cleavage planes result in a dipolar surface, which is therefore likely to be less stable. While dipolar surfaces can exist, this normally leads to twinning of crystals or strong solvation in order to anul the dipole. Alternatively, the surfaces can reconstruct in order to remove the dipole. This typically involves the creation of cation or anion vacancies at the surface or chemical modification. When an ion is conceptually removed, in practice it is moved to the bottom of the surface slab in a calculation in order to maintain charge neutrality [91]. Real examples of the three types of surface are presented for polymorphs of calcium carbonate in Figure 1.2.

Figure 1.1: The three types of surface as categorised by Tasker (a) type I, where each layer consists of coplanar anions and cations, and all planar cuts lead to a nonpolar surface; (b) type IIa, where the anions and cations comprising the layers are not coplanar, but which allows for some surface cuts to split the layers in such a way as to produce no dipole; (c) type IIb, which is as per type IIa, except that some ions must be moved from the surface to the bottom of region 2 in order to achieve zero dipole; and (d) type III, where there are alternating layers of cations and anions, and all possible planar cuts result in a surface with a dipole moment.

(a)	(b)
(c)	(d)

Figure 1.2: The examples of the three types of surface illustrated using polymorphs of calcium carbonate; (a) type I (b) type II and (c) type III. The type I and III surfaces illustrated are from the crystal structure of calcite, which due to its high symmetry has no type II surfaces present. The type II example is from aragonite.

[Type I] [Type II][Type III]

From the above, it can be seen that often the creation of the surface is a significant task in its own right and can be a complex process. In previous programs for surface calculations, the structure manipulation has usually been performed via the input deck of the code. Clearly, when reconstructions are involved this can become rather unwieldy. As a result, a different strategy has been adopted for use with GULP. All construction of the surface and structural manipulation is performed independently from the main forcefield engine by graphical means. This is achieved through an interface to the freely available program GDIS developed by Dr Sean Fleming (http://gdis.seul.org/). This interface allows surfaces to be specified by their Miller indices, valid shifts to be searched for, and the geometries then to be manipulated, if necessary. Once the desired surface structure has been generated, then the necessary GULP calculation can be performed.

Surface energy

The thermodynamic penalty for cleaving a surface from a bulk material is measured according to the surface energy. Given a bulk energy of $U_{bulk}$ and an energy for the same system with a surface created of $U_{surface}$, then the surface energy, $\Delta U_{SE}$, is defined as an intensive quantity according to:

$$\Delta U_{SE}=\frac{\left(U_{surface}-U_{bulk}\right)}{A}$$

where $A$ is the resulting surface area. By definition, for any stable material the surface energy will be endothermic. A calculated negative surface energy implies that a material should dissociate, i.e. the crystal should disperse into the surrounding medium.

There are two practical approaches that are widely used to determine the surface energy by computational means. In the first, a two-dimensional slab of material is created from the bulk, thus creating two surfaces overall. This method has the advantage that it can be used within programs that only allow for three-dimensional boundary conditions through the introduction of a sufficently large vacuum gap between the images, such that the surfaces don't interact across the vacuum. In addition, it becomes necessary to assess whether the slab is also thick enough since the properties must converge to those of the bulk at the centre of the slab. In the second method, a single surface is created by employing a two region strategy, as shown in Figure 1.3. Here the solid is divided into region 1, which contains the surface and all layers of atoms below it that exhibit a significant atomic relaxation, and region 2, which contains the rest of the bulk material where it is assumed that no displacements from the three-dimensional crystal structure are induced. In practice, only the atoms of region 2 that have an interaction with region 1 need be explicitly considered, and so the depth of region 2 is controlled by the cut-offs of the forcefield and the Parry sum used for the electrostatics. This second method is the most efficient and precise for atomistic techniques. However, it is considerably harder to extend to quantum mechanical methods since the electronic perturbations may extend further into the bulk and embedding, typically via Green's functions, is required to determine the influence of the electronic structure of region 2 on region 1. Through the GDIS interface, it is possible to automatically estimate the required region 1 and 2 sizes needed to converge the surface energy (with a default tolerance of $0.001J/m^{2}$) based on the unrelaxed surface energies. However, if there are strong relaxations of the surface, it may be necessary to further verify that the relaxed surface energy is sufficiently converged.

Figure 1.3: The two region surface simulation cell viewed at right angles to the surface normal, where solid vertical lines denote the boundaries between two-dimensional periodic images of the cell and the dash line indicates the boundary between region 1 and the frozen region 2.

The total energy for a surface calculation comprising two regions can be written in terms of the interaction energies within, and between, the different regions:

$$U_{tot}=U_{11}+U_{12}+U_{22}$$

The energy of region 2 with itself, $U_{22}$, is not particularly meaningful, since the region 2 is just a partial representation of the effectively infinite bulk material below, and any given particle within this region will not experience the full set of interactions that it should. However, this term is just an additive constant that is unaltered on energy minimisation, or any other displacement of region 1. Consequently it can be ignored in calculations. In the two region model, the energy required to determine the surface energy is given by:

$$U_{surface}=U_{11}+\frac{1}{2}U_{12}$$

Note that while the above energy only includes half of the region 1-region 2 interaction energy, the objective quantity for energy minimisation is the total energy, which includes the full value of $U_{12}$. This is because the energy change of region 2 must be allowed for when optimising.

Attachment energy

The surface energy, described above, provides a measure of the thermodynamic stability of a cleavage plane. However, there is a widely used alternative criterion which is the attachment energy. This is the energy associated with the addition of a stoichiometric layer of material on to the surface cut:

$$U_{attachment}=U_{tot}^{n+1}-U_{tot}^{n}-U_{tot}^{1}$$

where $U_{tot}^{n}$ represents the total internal energy of a surface model consisting of $n$ growth layers, and $U_{tot}^{1}$ is the energy of the growth layer alone. For any stable material, this implies that the attachment energy must be an exothermic quantity. In practice, the calculation of the attachment energy is obtained from the energy of interaction of the growth layer at the surface with the rest of the underlying material. This benefits from the fact that the attachment energy can be obtained from a single calculation, just as is the case for the surface energy, rather than by performing the actual addition of a layer as part of a multistage process.

Although the attachment energy is also a strictly thermodynamic quantity, it is often regarded as representing the kinetics of crystal growth because of the conceptual link between the ease of the addition of a growth layer and the rate at which a surface is added to. Consequently, those faces where the attachment energy is most exothermic will tend to grow most rapidly.

Morphology

The morphology of a crystal is the macroscopic shape that it adopts. Because this can be readily observed for nearly all materials, either under an electron microscope or, in the case of many naturally occuring minerals, by visual inspection with the naked eye, it should provide a ready means to test the validity of a simulation model. Of course, the reality is rather more complex, since the morphology is sensitive to the presence of impurities, the nature of the solvent used as the growth medium, and many other factors relating to the sample preparation. Consequently there are materials where many different morphologies can be observed for the same compound. A classic example, is that of calcite ($CaCO_{3}$), where there are both several different polymorphs of the bulk material and several hundred different known morphologies. Many of these variations result from biomineralisation by different species. Despite this, for many pure inorganic materials morphological prediction using atomistic techniques is surprisingly successful.

Crystal morphologies can be calculated based on either the surface energy or attachment energy, which are typically taken to represent growth under conditions of thermodynamic and kinetic control, respectively. In order to do this, it is first necessary to determine the objective quantity for all significant faces. Given that dipolar faces are usually unstable, the number of likely cleavage planes for most materials is actually considerably smaller than initially might be conceived of based on permutations of the Miller indices. Furthermore, where there is space group symmetry present for the bulk, many surface planes are equivalent, thus reducing the number of unique faces. Finally, only those faces with the largest interplanar spacings are likely to appear in the morphology [92]. The actual morphology is generated as a three-dimensional Wulff plot [93]. Here the ratio of the surface normal distances of all planes from the centre of the polyhedron are determined according to the either the surface or attachment energy. The final shape of the polyhedron is then determined by the intersection of the cleavage planes. Unstable surfaces lie outside the polyhedron and never intersect. Morphological plots can also be produced through the GDIS interface to GULP.

In the equilibrium morphology approach, the contribution of a given plane to the total surface area is inversely proportional to its surface energy [94]. For the growth morphology methodology [95], the surface area contribution is again inversely proportional, but now to the negative of the attachment energy. This is because surfaces with a highly exothermic attachment energy will rapidly grow out of the morphology to leave the slow growing bounding faces.

Surface phonons

The calculation of the phonons of a 2-D slab is exactly analogous to that for the three dimensional case, except that the Brillouin zone is also now two dimensional leading to dispersion only in the $xy$-plane (taking $z$ to be the direction of the surface normal).

Turning to consider the standard two region surface model, there are some important issues to consider. Because region 2 is quasi-infinite, it is only possible to determine the phonons for the particles in region 1. Therefore, the dynamical matrix is constructed based on region 1 alone, but with contributions to the on-diagonal matrix blocks from the potential experienced due to a rigid, non-vibrating region 2. Consequently it is assumed that the vibrations of the two regions are completely decoupled. Since this is an approximation, the frequencies near the interface of the regions will be slightly in error, particularly in the low frequency regime where coupling is generally strongest. However, the surface modes, which are usually the ones of primary interest, will be less affected. As always, it is essential to monitor the convergence of all quantities with respect to increasing the region 1 depth. Finally, because region 1 is vibrating under the influence of an external potential, the first three frequencies at the $\Gamma$-point will no longer be zero, though they typically will be small.