Our discussion so far has demonstrated how to localise the NGWFs in real-space in a manner that ensures that they are composed by a number of delta functions that is constant with system-size. On the other hand, each delta function is expanded in the plane-waves that can be supported by the regular grid of the simulation cell. The number of these plane-waves is proportional to the size of the system. As a result, the cost of a calculation still scales cubically with system-size as in the traditional plane-wave case.
In order to reduce the computational cost we must restrict the number of plane-waves that contribute to each delta function so that it is independent of system-size. This is not such a straightforward matter as the restriction of the NGWFs in real-space. Several factors have to be taken into account the most important of which is the Hermiticity of operators in integrals between NGWFs and also a representation of operators that is consistent when they are acting on different NGWFs. We have investigated these matters in detail in the context of the evaluation of kinetic energy integrals with NGWFs in a previous paper  where we have proposed the `FFT box' technique as an accurate and efficient solution. It is shown there that all the imposed conditions are satisfied if the plane-waves that are used to expand the delta functions are restricted to belong to a miniature simulation cell which we call the `FFT box'. The FFT box must be large enough to contain any possible orientation of overlapping NGWFs, and must be a parallelepiped with a shape commensurate with the simulation cell. It should contain a regular grid which is a subset of the simulation cell grid and thus the origin of the grid of the FFT box should coincide with a particular grid point of the simulation cell. FFT techniques using smaller boxes have been used in the past to take advantage of localised functions: Pasquarello et al. [44,45] use them to efficiently deal with the augmentation charges that arise when using ultrasoft pseudopotentials; Hutter et. al  use them to FFT localised Gaussian orbitals. However, to the best of our knowledge, this is the first time that a small FFT box has been used to define a systematic basis set for the entire total energy calculation.
Figure 1 demonstrates in two dimensions the FFT box inside the simulation cell as defined for two overlapping NGWFs and . In the same Figure the regular grid is also visible and we can observe that a portion of it is included in the FFT box. For the subset of grid points inside the FFT box we can define delta functions as we did for the simulation cell. We represent these FFT box delta functions by and in general we follow the convention of using lowercase letters to denote quantities associated with the FFT box and uppercase letters for quantities associated with the simulation cell.
Based on the knowledge gained through the use of the
FFT box to compute kinetic energy integrals,
we extend here its use to the calculation of all the terms
that constitute the total energy (section
2) and all the terms needed for
its subsequent optimisation with respect to
the density kernel and NGWF coefficients
(section 3). For this purpose we need to be
able to express or `project' a function from the simulation cell
to the FFT box and back. Assuming the function is expressed
entirely by delta functions
on grid points common to the simulation cell and the FFT box,
this task is straighforward
and all we have to do is to rewrite the function
as a linear combination of the FFT box delta functions
that lie on the same centres as the simulation cell
delta functions. This task is expressed formally
When this operator acts on a function in the
simulation cell, it maps it onto a function in the
FFT box. This is demonstrated for the function
of equation (7) by the following
expression which can also be considered as a definition
(a more explicit expression for
is given in Appendix B):
With this compact notation it is relatively straightforward to devise a way of calculating the total energy by using only the delta functions (and hence the plane-waves) periodic in the FFT box as the basis set for each NWGF. It is also equally straightforward to write formulae that represent this process in a concise way.
We start with the charge density of equation (2),
which should of course extend over the whole simulation cell, however
of NGWFs need not, and do not, when they are calculated
with the FFT box technique. Therefore the charge density is calculated by
replacing these pair contributions in (10)
by the following expression:
The matrix elements of one-electron operators, such as the kinetic
energy and non-local potential can be easily calculated in the
FFT box rather than in the simulation cell by simply transferring the
NGWFs in the FFT box before evaluating the integrals.
The notation for this set of operations simply involves
`sandwiching' the one-electon operator between the standard
grid projector operator and its adjoint. The kinetic energy
expression of equation (12) becomes
The matrix elements of the Hartree and local potentials can
be calculated with the same ease provided however we have
a way to transfer these potentials to the FFT box. This
is indeed possible as these local potentials are expressed in terms
of the fine grid delta functions and therefore we can
transfer them if we use the fine grid delta function
. Therefore the
Hartree energy of equation
Nothing changes in the evaluation of the exchange-correlation energy which is simply an integral over the fine grid of the whole simulation cell except of course for the fact that the charge density is now calculated by summing terms in equation (10) in place of the .
As far as the optimisation of the energy with respect to the density kernel is concerned, the results of section 3 are still valid provided the Hamiltonian matrix elements that constitute the density kernel gradient (23) are calculated with the FFT box method as we have described above.
In the same way, the calculation of the total energy gradient
with respect to the NGWF expansion coefficients of equation
(24) can be calculated in the FFT box provided
we use only the part of the Kohn-Sham Hamiltonian that exists
in the FFT box. Therefore all we have to do is to substitute
of equation (24) by
As far as implementation is concerned, there are many more issues and algorithmic details about the use of the FFT box that are beyond the scope of this paper. We will describe these in another paper .