Calculations were performed for a cubic simulation cell containing 216 silicon atoms and using pseudopotentials generated according to the method of Troullier and Martins . The kinetic energy cut-off for the basis-set was set at 200 eV and the FFT grid contained points. The basis-set contained spherical-waves with angular momentum components up to , 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 -point. A value of 100 eV was used for the penalty functional prefactor .
In figure 9.1 we plot the total energy per atom against the support region radius for a density-kernel cut-off of 4.0 Å. In figure 9.2 we plot the total energy per atom against the density-kernel cut-off for a support region radius 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. . 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.
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 Å and Å, 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 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.
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 . For a lattice parameter
Å, which was used for the calculations in section
9.1.1, we used values of
Å and an energy cut-off of 200 eV.
As the volume was changed, the support region radius
and energy cut-off
were changed proportionally.
The parameters used are listed in table 9.1.
Fitting the data to the Birch-Murnaghan equation of state allows values
for the equilibrium volume and bulk modulus to be calculated, which are
compared with the results obtained from a CASTEP calculation using the
same pseudopotential and to experiment  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.