Optimize fft_grid computation

[antoine-levitt]

On a test (4x1x1 supercell of silicon), our prod(fft_size) is about 25% higher than abinit's. Abinit uses a more clever iterative algorithm (see getng function)

[antoine-levitt]

This is not fixed by using an explicit search for the rectangular box containing the sphere of radius 2Gmax. I'm not completely sure how to do it without a quadratic algorithm.

[antoine-levitt]

For reference here's the code for doing that

# returns the lengths of the bounding rectangle in reciprocal space
# that encloses the sphere of radius Gmax
function bounding_rectangle(lattice::AbstractMatrix{T}, Gmax; tol=1e-8, alg=:crude) where {T}
    Glims = [0, 0, 0]
    if alg == :crude || det(lattice) == 0
        Glims_ = [norm(lattice[:, i]) / 2T(π) * Gmax for i in 1:3]
        # Round up, unless exactly zero (in which case keep it zero in
        # order to just have one G vector for 1D or 2D systems)
        for i = 1:3
            if Glims_[i] != 0
                Glims[i] = ceil(Int, Glims_[i] .- tol)
            end
        end
    else
        Glims_crude = bounding_rectangle(lattice, Gmax; tol=tol, alg=:crude)
        recip_lattice = 2T(π)*inv(lattice')
        for G1 = -Glims_crude[1]:Glims_crude[1]
            for G2 = -Glims_crude[2]:Glims_crude[2]
                for G3 = -Glims_crude[3]:Glims_crude[3]
                    G = recip_lattice * [G1; G2; G3]
                    if norm(G) ≤ Gmax
                        Glims .= max.(Glims, abs.((G1, G2, G3)))
                    end
                end
            end
        end
    end
    Glims
end

It doesn't decrease significantly the size of the grid

[antoine-levitt]

OK, I've implemented the exhaustive search algorithm. Eg for the 4x1x1 silicon with gamma-only sampling at Ecut=30, I get (144, 35, 35) (old algorithm) (147, 40, 40) (new algorithm) Too bad I forgot to record abinit's values while I was at it the other day... @mfherbst do you happen to have any abinit output file I could compare to?

[mfherbst]

No, but I'll run it quickly.

[mfherbst]

So on my experiments (current master) I got:

DFTK:  [150, 40, 40]
ABINIT: 144      36      36

[antoine-levitt]

Without the optimize flag right?

[mfherbst]

current master PR not merged.

[antoine-levitt]

Ah, right. So we match abinit, good

[antoine-levitt]

The status of that is that the optimized algorithm is implemented but has a significant overhead, so it's turned off by default