9.1 Bulk crystalline silicon

Calculations were performed for a cubic simulation cell containing 216 silicon atoms and using pseudopotentials generated according to the method of Troullier and Martins [74]. The kinetic energy cut-off for the basis-set was set at 200 eV and the FFT grid contained $72 \times 72 \times 72$ points. The basis-set contained spherical-waves with angular momentum components up to $\ell=2$, and the support regions were chosen to be centred on the bonds, with one support function per region, so that no unoccupied bands were included in this calculation. The Brillouin zone was sampled using only the $\Gamma$-point. A value of 100 eV was used for the penalty functional prefactor $\alpha$.

Figure 9.1: Convergence of total energy with respect to support region radius.

Figure 9.2: Convergence of total energy with respect to density-kernel cut-off.

9.1.1 Convergence with density-matrix cut-off

In figure 9.1 we plot the total energy per atom against the support region radius $r_{\mathrm{reg}}$ for a density-kernel cut-off $r_K$ of 4.0 Å. In figure 9.2 we plot the total energy per atom against the density-kernel cut-off $r_K$ for a support region radius $r_{\mathrm{reg}}$ of 3.1 Å. In both cases, the energy converges to its limiting value from above, as expected, since the total energy is variational with respect to the density-matrix cut-off.

These results agree roughly with the calculations of Hernández et al. [127]. In their case, they used atom-centred support regions and included as many unoccupied bands as occupied bands, so that the convergence with respect to density-matrix cut-off should be more rapid in their case. They also used a local pseudopotential, which reduces the range of the Hamiltonian. Our calculations suggest a combined density-matrix cut-off of the order of 7.0 Å to obtain the same accuracy as they obtained with a cut-off of about 6.0 Å. We note that the band gap of silicon is relatively small (particularly within the LDA) so that the density-matrix decay is therefore slow. This makes silicon a difficult test case, and we can be confident that if we obtain reasonable results in this system, we shall be successful in others.

Figure 9.3: Electronic density of silicon in the (110) plane (units Å$^{-3}$).

Figure 9.4: Diamond structure of silicon, highlighting a {110} plane.

Figure 9.5: Electronic density of silicon in the (110) plane calculated by the CASTEP code (units Å$^{-3}$).

Figure 9.6: Difference between electronic densities of silicon calculated by the linear-scaling method and CASTEP (units Å$^{-3}$).

9.1.2 Electronic density

Since the minimising density-matrix is only approximately idempotent, the electronic density derived from it is not the exact ground-state density. However, in figure 9.3 we plot the electronic density in the (110) plane (containing the atoms highlighted in silver in figure 9.4) obtained from the minimising density-matrix with $r_{\mathrm{reg}} = 3.1$ Å and $r_K = 6.0$ Å, and observe that it is still qualitatively correct. Since electronic densities are used primarily for visualisation purposes, rather than for quantitative analysis, the information most commonly required can still be obtained from the minimising density-matrix.

In figure 9.5 we plot the electronic density obtained using the CASTEP plane-wave code using the same pseudopotential and energy cut-off, and equivalent Brillouin zone sampling (a $3 \times 3 \times 3$ Monkhorst-Pack mesh for an 8-atom unit cell). We observe that there is less density concentrated in the bonds for the CASTEP density compared with the linear-scaling density, and this is to be expected since in the linear-scaling calculation, the lowest energy (i.e. most bonding) orbital is over-occupied and the highest energy (i.e. least bonding) orbital is under-occupied. In figure 9.6 we plot the difference obtained by subtracting the CASTEP density from the linear-scaling density and confirm this observation.

9.1.3 Predictions of physical properties

The converged value of the total energy agrees with the CASTEP value (again with equivalent energy cut-off, Brillouin zone sampling and pseudopotential) only within 3%. However, if energy differences are caused by density-matrix variations which are short-ranged, then they will be much better converged than absolute energies. In figure 9.7 we plot the total energy against volume $V$. For a lattice parameter $a = 5.430$ Å, which was used for the calculations in section 9.1.1, we used values of $r_{\mathrm{reg}} = 3.10$ Å, $r_K = 6.00$ Å and an energy cut-off of 200 eV. As the volume $V$ was changed, the support region radius $r_{\mathrm{reg}}$, density-kernel cut-off $r_K$ and energy cut-off $E_{\mathrm{cut}}$ were changed proportionally. The parameters used are listed in table 9.1.

Table 9.1: Parameters used to calculate energy-volume curve.

$a$ / Å	$V$ / Å$^3$	$r_{\mathrm{reg}}$ / Å	$r_K$ / Å	$E_{\mathrm{cut}}$ / eV
5.31	149.72	3.03	5.87	209
5.35	153.13	3.05	5.91	206
5.39	156.59	3.08	5.96	203
5.43	160.10	3.10	6.00	200
5.47	163.67	3.12	6.04	197
5.51	167.28	3.15	6.09	194

Figure 9.7: Energy-volume curve for silicon. The line is the best fit to the Birch-Murnaghan equation of state.

Fitting the data to the Birch-Murnaghan equation of state allows values for the equilibrium volume $V$ and bulk modulus $B$ to be calculated, which are compared with the results obtained from a CASTEP calculation using the same pseudopotential and to experiment [174] in table 9.2. Generally we expect to obtain lattice parameters to within 2% and bulk moduli to within 10%. The bulk modulus is very sensitive to the data, so that whereas the prediction of the lattice parameter is in excellent agreement with both the CASTEP and experimental values, the value of the bulk modulus predicted by the linear-scaling method is about 8% too large. These results indicate that the linear-scaling calculation is not quite fully converged with the set of parameters used, but are very encouraging overall.

Table 9.2: Comparison of calculated and experimental data for silicon.

Calculation	Linear-scaling	CASTEP	Experiment
$a$ / Å	5.423	5.390	5.430
$V$ / Å$^3$	159.47	156.56	160.10
$B$ / GPa	108.8	101.7	100.0