10.2 Further work

Starting from the underlying quantum-mechanical theory, in this dissertation an original scheme to perform linear-scaling total energy calculations based upon the density-matrix and using a new basis-set has been outlined and its computational implementation discussed. The results obtained are very promising, but there is still much work to be done in developing and optimising the scheme e.g. the occupation number preconditioning outlined in section 7.5. One of the most computationally expensive steps is currently the evaluation of the electronic density on the real-space grid. By choosing to minimise the non-interacting Kohn-Sham energy non-self-consistently rather than the interacting energy self-consistently the number of times this evaluation has to be performed can be greatly reduced. Self-consistency can then be obtained by density mixing [179,180] . The problem of ``charge-sloshing'' which can arise in this case in traditional methods appears not to be present when localised functions are used, and can anyway be eliminated[181,83].

Once an efficient electronic minimisation scheme is in place, the next step is to calculate ionic forces in order to perform ionic relaxation or even molecular dynamics. First it is worth noting that there are of ${\cal O}(N)$ ions in the system, and the computational effort to calculate the force on each (from the local pseudopotential) is also of ${\cal O}(N)$ so that a computational effort of ${\cal O}(N^2)$ is necessary to calculate all the forces. Secondly, the typical time-scale of molecular dynamics required by a system of volume $V$ is $\tau \propto V^{2 \over 3}$ so that this introduces an extra factor of ${\cal O}(N^{2 \over 3})$ in the computational effort. The Hellmann-Feynman forces [182,183] resulting from the derivative of the Hamiltonian with respect to the ionic positions have a contribution from the local pseudopotential (which is calculated using the reciprocal-space grid) and the non-local pseudopotential (which is calculated by taking analytic derivatives of the results for the non-local pseudopotential matrix elements). These forces have in fact been calculated and tested using the spherical-wave basis-set. However, because the support regions move with the ions, other contributions known as Pulay forces [184] arise which must be calculated. Since the correction to the energy derived from the penalty functional can be expressed analytically, it should be possible to calculate accurate forces which are consistent with the corrected energy.

Given the expected improvement in efficiency of the method when applied to molecular rather than crystalline systems (section 9.2.2), serious consideration should be given to converting the code to perform calculations on clusters, rather than using the supercell approximation. This would have the added advantage of eliminating the ${\cal O}(N^2)$ step to calculate the structure factor and the calculation of the ion-ion interaction energy (although still an ${\cal O}(N^2)$ step) would also become much cheaper.

In the long term it is essential that advantage be taken of the natural way in which linear-scaling methods based in real-space can be separated into simultaneous calculation of interacting fragments. This property means that the problem lends itself to implementation on massively parallel computers [185]. Even with the serial code on a single workstation it has been possible to model 512 silicon atoms, and there exist parallel computers consisting of a few hundred nodes each with the same power, so that in principle calculations of 100000 atoms are feasible. The development of schemes which will exploit the advantages of parallel processing is therefore essential to reap the full benefit of linear-scaling methods.