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 cutoff for the basisset was set at 200 eV and the FFT grid contained points. The basisset contained sphericalwaves 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 densitykernel cutoff of 4.0 Å. In figure 9.2 we plot the total energy per atom against the densitykernel cutoff 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 densitymatrix cutoff.
These results agree roughly with the calculations of Hernández et al. [127]. In their case, they used atomcentred support regions and included as many unoccupied bands as occupied bands, so that the convergence with respect to densitymatrix cutoff 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 densitymatrix cutoff of the order of 7.0 Å to obtain the same accuracy as they obtained with a cutoff of about 6.0 Å. We note that the band gap of silicon is relatively small (particularly within the LDA) so that the densitymatrix 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 densitymatrix is only approximately idempotent, the electronic density derived from it is not the exact groundstate 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 densitymatrix 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 densitymatrix.
In figure 9.5 we plot the electronic density obtained using the CASTEP planewave code using the same pseudopotential and energy cutoff, and equivalent Brillouin zone sampling (a MonkhorstPack mesh for an 8atom unit cell). We observe that there is less density concentrated in the bonds for the CASTEP density compared with the linearscaling density, and this is to be expected since in the linearscaling calculation, the lowest energy (i.e. most bonding) orbital is overoccupied and the highest energy (i.e. least bonding) orbital is underoccupied. In figure 9.6 we plot the difference obtained by subtracting the CASTEP density from the linearscaling density and confirm this observation.
The converged value of the total energy agrees with the CASTEP value
(again with equivalent energy cutoff, Brillouin zone sampling and
pseudopotential) only within 3%. However, if energy differences are
caused by densitymatrix variations which are shortranged, 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 cutoff of 200 eV.
As the volume was changed, the support region radius
,
densitykernel cutoff
and energy cutoff
were changed proportionally.
The parameters used are listed in table 9.1.

Fitting the data to the BirchMurnaghan 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 [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 linearscaling method is
about 8% too large. These results indicate that the linearscaling
calculation is not quite fully converged with the set of parameters used,
but are very encouraging overall.
