The 2^nd derivatives of the electronic energy with respect to the nuclear coordinates (the hessian matrix) are the basic quantities required for calculating vibrational frequencies calculated at a stationary point, thus enabling the characterisation of that point. At a geometry that corresponds to a minimum of the energy all the vibrational frequencies will be real, whereas at a transition state there will be one negative eigenvalue, leading to a so called imaginary frequency. Finally, with vibrational frequencies, the vibrational modes are also obtained. These modes can be used as search directions to find molecular conformations that correspond to stable structures or transition states. In addition to this the 2^nd derivatives require algorithmic steps that provide elementary quantities for the calculation of infrared spectra and polarisabilities. Thus there are a number of key properties that become available once 2^nd derivatives of the energy can be calculated.

The initial implementation of the DFT hessian within the CCP1-developed DFT closely follows Handy et. al. [1] and Johnson et. al. [2] and was reported on previously [3]. However at that point no attempt had been made to reduce the costs of the calculations. This resulted in a formalism that scales as O(N⁴) where N is the number of atoms; this is clearly unacceptable as it prohibits calculations on any sizeable molecule. The performance can be greatly improved by eliminating negligibly small terms and by moving work out of inner-most loops. Examples of negligibly small terms are basis functions situated on atoms far from the current batch of grid points and contributions from small elements of the density matrix. We call the process of eliminating small terms “screening”. The most dramatic examples of moving work out of inner-most loops can be found comparing the MO-based and AO-based expressions in the CPHF equations. In the AO-basis many quantities involve matrices of dimension N_basis ², whereas in the MO-basis the dimension is only N_occ x N_unocc . Given that N_basis equals N_occ+ N_unocc, then the maximum number of elements to treat in the MO-basis is a quarter of that in the AO-basis. Thus a factor of 4 in performance can be achieved simply by choosing whether one transforms all quantities to the MO-basis within the DFT module or outside the DFT module.

Although the above example seems to suggest it is always worthwhile to evaluate all quantities in the MO-basis, in practice the situation is more complex. The reason is that MO’s are dispersed over the whole molecule, removing the potential for applying screening techniques on the values of the MO’s. Localising the MO’s is not an option because that would result in a non-diagonal Fock matrix and corresponding performance penalties. AO’s on the other hand are localised by definition and screening techniques can be applied straightforwardly and effectively. This means in the context of the example given above, that matrices in the MO-basis still scale like N_occ x N_unocc for large molecules, whereas in the AO-basis we need to treat only N_non-zero x N_non-zero elements, where N_non-zero is the number of non-zero AO’s and approaches a constant. Clearly in the limit of large molecules an AO-based algorithm will become more efficient.

The conclusion of the above discussion is that ideally we have to be able to perform the calculations in either AO or MO basis and select the most appropriate basis on a case-by-case basis. In the limit of small molecules, it is obvious that the MO basis algorithm should be used as screening techniques cannot be expected to achieve significant cost reduction. In the limit of very large molecules the obvious choice has to be to use the AO basis algorithm. For the intermediate sized molecules, which are probably the most widely studied systems, there are two options:

• First, examine the problem in hand and decide which of the AO basis or MO basis algorithm is the most efficient. For each step in the calculation the most appropriate algorithm is then used. This approach requires the calculation of statistics on which to base this decision, for which the integrator must be run at least once. The most appropriate statistic is the average number of non-zero AO’s per grid point. This statistic can be collected at little cost and used straightforwardly in expressions like those above to decide which approach is the most efficient.

• Second, as the actual number of non-zero AO’s differs for each batch of grid points, it is clear that decisions to treat either all batches in the AO or MO basis have to be sub-optimal compared to treating each batch in its most efficient way. Essentially this approach is easily implemented but suffers from the drawback that it is memory intensive as some data structures will have to be stored twice, in both the AO and MO basis. Therefore it may not always be possible to exploit this efficient scheme.

Having presented an outline of the various performance enhancement options that are currently being evaluated, it seems appropriate to present an example. In figure 1 the infrared spectrum of nitrobenzene is shown. The calculation was performed using the B3LYP functional as a typical example of a gradient corrected functional. The 6-31G* basis set was used, resulting in 145 basis functions of which 32 are occupied and 113 are unoccupied. In this case the MO basis calculation was 5.6 times faster than that in the AO basis. This is consistent with the fact that solving the CPHF equations is the step that dominates the cost and that the tight basis function cut-off of 10^-13 effectively eliminates all screening. Assuming that no screening is used, one would predict from theory that the calculation in the MO basis should be about 5.8 times faster than that in the AO basis.

Figure 1. The infrared spectrum of nitrobenzene calculated with GAMESS-UK using the 6-31G* basis set and the B3LYP functional.

references

[1] N.C. Handy, D.J. Tozer, G.J. Laming, C.W. Murray and R.D. Amos, Isr. J. Chem. 33 (1993) 331-344.

[2] B.G. Johnson and M.J. Fisch, J. Chem. Phys. 100 (1994) 7429-7442.

[3] H.J.J. van Dam, 2nd Derivatives of the Electronic Energy in Density Functional Theory, CCLRC Daresbury Laboratory, Warrington, 2001, DL-TR-01-002.

[4] A.D. Becke, J. Chem. Phys. 88 (1988) 2547-2553.

[5] M.F. Guest, J.H. van Lenthe, J. Kendrick, and P. Sherwood, GAMESS-UK - A package of ab initio programs, Computing for Science Ltd., Computational Science and Engineering Department, Daresbury Laboratory, Warrington, Cheshire, WA4 4AD, UK, http://www.cfs.dl.ac.uk.