3.2 Format of arrays containing charge density, potential, etc.

The index of arrays used to store functions defined on 3D meshes is actually a shorthand for three indices, following the FORTRAN convention ("leftmost index runs faster"). An example will explain this better. Suppose you have a 3D array psi(nr1x,nr2x,nr3x). FORTRAN compilers store this array sequentially in the computer RAM in the following way:

        psi(   1,   1,   1)
        psi(   2,   1,   1)
        ...
        psi(nr1x,   1,   1)
        psi(   1,   2,   1)
        psi(   2,   2,   1)
        ...
        psi(nr1x,   2,   1)
        ...
        ...
        psi(nr1x,nr2x,   1)
        ...
        psi(nr1x,nr2x,nr3x)
etc

Let ind be the position of the (i,j,k) element in the above list: the following relation

        ind = i + (j - 1) * nr1x + (k - 1) *  nr2x * nr1x

holds. This should clarify the relation between 1D and 3D indexing. In real space, the (i,j,k) point of the FFT grid with dimensions nr1 ($\le$nr1x), nr2 ($\le$nr2x), , nr3 ($\le$nr3x), is

r_ijk = $${\frac{{i-1}}{{nr1}}}$$$$\tau_{1}^{}$$ + $${\frac{{j-1}}{{nr2}}}$$$$\tau_{2}^{}$$ + $${\frac{{k-1}}{{nr3}}}$$$$\tau_{3}^{}$$

where the $\tau_{i}^{}$ are the basis vectors of the Bravais lattice. The latter are stored row-wise in the at array: $\tau_{1}^{}$ = at(:, 1), $\tau_{2}^{}$ = at(:, 2), $\tau_{3}^{}$ = at(:, 3).

The distinction between the dimensions of the FFT grid, (nr1,nr2,nr3) and the physical dimensions of the array, (nr1x,nr2x,nr3x) is done only because it is computationally convenient in some cases that the two sets are not the same. In particular, it is often convenient to have nrx1=nr1+1 to reduce memory conflicts.