LDA DFT simulations of an isolated silicon donor on the (110) surface of GaAs

The convergence of the band gap state of a single silicon dopant on the (110) surface of GaAs was investigated. By simulating different sized super-cells we were able to show that a 3x4 super-cell provides a well converged calculation for modelling an isolated dopant, with the total energy being converged to 1 part in 1000. The local density of the silicon band gap state was then checked against a number of more intensive calculations and was found to be well converged, with an eigenvalue accurate to within 3 meV.


Introduction
Recent STM experiments have shown interesting properties for silicon on the GaAs (110) surface, including positional bistability [1] and enhanced binding energy [2]. We have performed LDA DFT calculations for an isolated silicon donor on the (110) surface of GaAs, and by looking in detail at the localised band gap state of the dopant we can obtain a direct link with the STM results [1,2] and other theoretical work [3]. This is an important first step in understanding the properties of all group IV dopants on the GaAs (110) surface.
In order to be confident in any theoretical simulation a host of convergence checks must be employed. When simulating isolated donors it is particularly important to check the convergence with respect to the size of the super-cell. Larger cells provide better approximations to isolated dopant systems, but computational costs place upper limits on what can realistically be modelled. A common choice of system size for looking at donors on the (110) surface is a 3×2 lateral super-cell [3] and work has also been done on systems as small as 2×2 [4], but to our knowledge no systematic tests have been reported to show the validity of cells of these sizes.
For this system we are interested in the highest occupied molecular orbital (HOMO). We calculate the total energy and the HOMO energy for the system and converge these with respect to E cut , k-point grid and the number of atoms in the system. Through the simulation of more computationally costly systems we show that, in addition to the energies, the LDOS for the HOMO is also well converged.

Computational Details
We used the ABINT software package [5] to perform LDA DFT calculations of the (110) surface of GaAs. The psuedopotentials used for the GaAs slab and the silicon dopant were taken from the set developed by the Fritz-Haber Institute which are available through the ABINIT website. We investigated the dependence on super-cell size using twelve super-cells, all consisting of 7 atomic layers and 12Å of vacuum but varying in their in-plane dimensions. The super-cells were: 2×2 (56 atom), 2×3 and 3×2 (84 atoms), 2×4 and 4×2 (112 atoms), 3×3 (126 atoms), 4×3 and 3×4 (168 atoms), 5×2 and 2×5 (140 atoms) and 5×3 and 3×5 (210 atoms). For all cells up to 126 atoms the silicon atom replaced a gallium atom in the surface layer of the slab and the entire system was allowed to fully relax, with a tolerance on the relaxation of 1 × 10 −4 Hartrees per Bohr. The super-cells larger than 126 atoms were constructed from the final positions of the relaxed 126 atom super-cell with the additional atoms added using the relaxed geometry of the clean GaAs slab. To improve the speed of the relaxation procedure the initial position of the silicon was calculated using a simple geometrical model that treats the atoms as close-packed, hard-walled spheres with a radius equal to the experimentally measured covalent radius. For each super-cell we checked the convergence with respect to E cut and k-point grid and across all the different cells we looked at the behaviour of the total energy and the energy of the HOMO.

Results and Discussions
The results of the convergence tests can be seen in figure 1. As a result of these checks we developed a standard calculation of the silicon dopant which had a cut-off energy of 30 Ry and a 2×2×1 k-point grid. These values of E cut and k-grid provided a convergence on the total energy of the system of 1 part in 1000 and a convergence of the HOMO energy of 1 part in 5000 (1 meV) for each super-cell size. This high level of convergence on the HOMO removed any source of uncertainty when it came to looking at the state's energy variation as a result of the changing system size.  In figure 2 we show the variation of the HOMO energy with respect to system size. Panels (a) and (b) indicate that the variation in the energy of the dopant state is much more sensitive to changes in the x dimension of the super-cell compared to changes in y. Taking the 210 atom cell to be the converged case we find in (c) that the super-cells with a width of 2 in the x direction produce energies considerably below the converged level regardless of the length in y. This is because the silicon dopant is only 7.82Å away from its closest periodic image and so a sizeable silicon-silicon interaction is expected. Once the x separation has changed to 3 in (d) we see the energy level settle on the converged value quite rapidly. From these results we chose the 3×4 cell as the size for our standard calculations, giving us a HOMO energy that was converged to within 10 meV, with respect to system size, of its value in the 210 super-cell.  The different behaviour of the energy for variations in x and y is primarily caused by the non-equivalence in the lengths of the primitive (110) cell; the y repeating length is √ 2 times longer than x. This causes increases in y to separate the periodic dopants much more rapidly than for changes in x. Another factor is that the (110) slab has zig-zag chains in the plane of the surface in the x direction which may lead to greater defect-defect interactions than are possible in the y direction [6].
Once we had converged the energy of the system we were able to check the robustness of the calculation by simulating the LDOS for a range of different input parameters. Three of the additional simulations involved increasing numerical variables of the calculation while the fourth was run with a different choice of pseudopotential. The different pseudopotentials used were taken from the Hartwigsen-Goedecker-Hutter (HGH) library on the ABINIT website. Figure  3 (a) shows the silicon LDOS for the standard system, a 3×4 cell with the cut-off energy and k-point grid described at the start of section 2. We also present the LDOS for the four other systems used to check its robustness. We show the result of the simulation with the different pseudopotentials (b), increasing E cut to 34 Ry (c), the k-point grid being increased to a 3×3×1 grid (d), and a system where the number of atoms was increased to 210 (e).
In figure 4 we have taken a line scan of the LDOS through the slab at the location of the    silicon. Here we can clearly see the agreement between the five different simulations, indicating that our standard values of E cut , k-points, super-cell size and pseudopotential lead to a very well converged simulation, with the peak of the density agreeing to better than 1 part in 30 in all cases. The inset in figure 4 shows the decay of the density into the vacuum on a log scale, indicating the continuation of the agreement into the vacuum tail of the density.