4. Tests of the algorithm

In this section, we present the results obtained from the calculations based on the block update procedure. Test cases are taken from the molecular chlorine and bulk crystalline silicon systems. The localized spherical-wave basis set [15] is used, where the basis functions are chosen to be centered on the atoms. We have used norm-conserving Troullier-Martins pseudopotentials [22] in the Kleinman-Bylander form [23], with angular momentum components up to $l=2$. We use an LDA [24] for the exchange and correlation term. Periodicity of the supercell is assumed and the $\Gamma$ point is used for the $k$-point sampling.

Figure 1: Convergence of the sum of eigenvalues $\Omega$ (eV) as a function of iteration number for the chlorine molecule calculations using a cutoff energy of 640 eV. $\Omega _0$ is the ``exact'' value from the direct matrix diagonalization. The curve labelled by empty diamonds ($\diamond$) corresponds to the calculation where $\tau$ is updated according to the highest kinetic energy of all approximate eigenvectors.

A chlorine molecule of bond length 2.0 Å is placed in a cubic box of side 10 Å. With a cutoff energy of 640 eV and the basis-function radius $R$ of 4.0 Å, a total of $2 \times 139 = 278$ basis functions are used. In Fig. 1 we display the convergence of the sum of Kohn-Sham eigenvalues toward the ``exact'' value obtained from direct matrix diagonalization, as a function of the iteration number. The convergence of solution is seen to be linear when the number of iterations is smaller than the number of basis functions. To investigate the effect of preconditioning on the convergence of the solution, we have used a number of fixed $\tau$ values. It is seen that the performance of the method improves with moderate preconditioning. Fig. 1 shows that $\tau$ should be about 10 eV for good convergence. We have performed another calculation with $\tau$ updated according to the highest kinetic energy of all approximate eigenvectors, which converges to 24 eV. This is the natural choice for $\tau$ used in other preconditioning schemes, and the performance of this calculation (the curve labelled by open diamonds $\diamond$) is seen to be rather similar to that of the `optimal' case with $\tau = 10$ eV. This method therefore allows $\tau$ to be chosen automatically, and optimized during the calculation, rather than being another parameter which the user must specify.

Figure 2: Convergence of the sum of eigenvalues $\Omega$ (eV) as a function of iteration number for the chlorine molecule calculation using a cutoff energy of 640 eV. The curves show the difference between the tensor-nature-preserving (TNP) calculations and non-tensor-nature-preserving (NTNP) calculations.

To investigate the importance of preserving the tensor nature of the search direction, we have performed calculations with the same cutoff energy of 640 eV on the molecular chlorine system, but this time with $S$ set to the identity matrix (this corresponds to the case where tensor nature of the search direction is not preserved) and the off-diagonal elements of $T$ set to zero (this corresponds to the diagonal approximation used in Ref. [14]) when we solve Eq. 22. The results of the calculations are presented in Fig. 2 where we have included the tensor-nature-preserving (TNP) curves for comparison. It is found that that the non-tensor-property-preserving (NTNP) cases fail to converge to the right solution. We conclude that it is essential to take tensor properties into account when one is dealing with a non-orthogonal basis set.

Figure 3: Convergence of the sum of eigenvalues $\Omega$ (eV) as a function of iteration number for the chlorine molecule calculations using a cutoff energy of 4800 eV. $\Omega _0$ is the ``exact'' value from the direct matrix diagonalization. The curve labelled by empty diamonds ($\diamond$) corresponds to the calculation where $\tau$ is updated according to the highest kinetic energy of all approximate eigenvectors.

With a cutoff energy as high as 4800 eV (a total of $2 \times 392 = 784$ basis functions are used in these calculations), Fig. 3 clearly indicates that it is crucial to use the preconditioning scheme. A comparison between Figs. 1 and 3 reveals that when the optimal value of $\tau$ is used, the number of iterations to achieve the same accuracy remains roughly the same, even though the number of basis functions has more than doubled, which shows that the preconditioning scheme is indeed working.

Figure 4: Convergence of the sum of eigenvalues $\Omega$ (eV) as a function of iteration number for the 64-atom Si crystal calculations using a cutoff energy of 200 eV. $\Omega _0$ is the ``exact'' value from the direct matrix diagonalization. The curve labelled by empty diamonds ($\diamond$) corresponds to the calculation where $\tau$ is updated according to the highest kinetic energy of all approximate eigenvectors.

Similar tests are performed on the bulk crystalline silicon system. The calculations on a 64-atom silicon unit cell are performed at the equilibrium lattice parameter of 5.43 Å with an energy cutoff of 200 eV. We have chosen $R$ to be 3.1 Å which is sufficient for this purpose. These settings result in a total of $64 \times 55 = 3520$ basis functions for the calculations. In Fig. 4 we we note that our `best' $\tau \sim 1$ eV is comparable with the value of 3.8 eV used by Bowler and Gillan [14]. We have performed another calculation with $\tau$ updated according to the highest kinetic energy of all approximate eigenvectors, which converges to 12 eV. The performance of this calculation is seen to be rather similar to that of the optimal cases with $\tau = 1$ or 10 eV.