This came up now with me computing DFPT for various structures (metals, semiconductors, insulators) with temperature=1e-3 irrespective of the structure. It worked fine on Al and Si, but NaCl (salt) failed once it got to the response solver. Here's a cooked down self-contained reproducer:

using DFTK
using PseudoPotentialData
using ForwardDiff
using LinearAlgebra
DFTK.setup_threading()

pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_5.standard.upf")
atoms = [
    ElementPsp(:Na, pseudopotentials),
    ElementPsp(:Cl, pseudopotentials),
]
positions = [[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]]

a = 10.92
lattice = (a/2) * [0 1 1; 1 0 1; 1 1 0]
Ecut = 44

function compute_stress(strain)
    new_lattice = (1 + strain) * lattice
    model = model_DFT(new_lattice, atoms, positions;
                      functionals=PBE(), kinetic_blowup=BlowupCHV(),
                      temperature=1e-3, smearing=Smearing.Gaussian())
    basis = PlaneWaveBasis(model; Ecut, kgrid=(4, 4, 4))
    scfres = self_consistent_field(basis; tol=1e-9)
    stress = compute_stresses_cart(scfres)
    tr(stress)
end

@show compute_stress(0.)  # This works fine
@show ForwardDiff.derivative(compute_stress, 0.)  # This gives NaN in response solver

Tracking this NaN down, I think I found the issue in chi0.jl:

For an insulator with a large band gap (such as NaCl with ~5.0 eV at PBE), we get numerical underflow in all the occupation derivatives, which in this case indeed leads to exact zeros δocc_tot = 0.0 and D = 0.0 thus δocc_tot / D = NaN .

It seems the quickest fix is just to check iszero(D) before that division. What I am not sure about is whether returning zero as the Fermi level derivative in this case is even valid.

With fixed temperature, but increasing band gap the Fermi level should be converging to the center of the gap, right? That would mean that it should asymptotically observe half the derivative of the gap. Does this make sense?