The FFT box technique

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 [43] 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* [46] 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[43],
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
by the
projection operator.
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
of
(a more explicit expression for
is given in Appendix B):

while the adjoint operator can act on a function in the FFT box and turn it into a function in the simulation cell in an analogous way.

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
its contributions
from pairs
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

(27) |

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
projector
. Therefore the
Hartree energy of equation
(14) becomes

(28) |

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

(29) |

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 [47].