If $H|n\rangle = E_n|n\rangle$ for each $R$ (with the dependence on $R$ suppressed for readability), then
$$\nabla_R \langle m|H|n\rangle = \nabla_RE_n\delta_{nm} = E_n(R)\bigg[\nabla_R \langle m|\bigg]|n\rangle + E_m\langle m|\bigg[\nabla_R|n\rangle\bigg] + \langle m|\big(\nabla_RH\big)| n\rangle$$
Noting that $\nabla_R\langle m|n\rangle = 0 \implies \bigg[\nabla_R\langle m|\bigg]|n\rangle = - \langle m|\bigg[\nabla_R |n\rangle\bigg]$, it follows that for $E_n\neq E_m$,
$$\langle m|\nabla_R|n\rangle = \frac{\langle m|\big(\nabla_RH\big)|n\rangle}{E_m-E_n}$$
I thought it was a strict equality for any energy scale, not a perturbative one [...]
This is not an approximate result. It can be obtained by expanding $|n\rangle$ in a perturbation series about $R=R_0$ and then reading off the first term, but that's the same as finding the derivative of a function by computing its power series and reading off the coefficient of the first term:
$$\frac{1}{1-x} = 1 + \color{red}{1}x + x^2 + \ldots \implies \frac{d}{dx}\left[\frac{1}{1-x}\right]\bigg|_{x=0} = \color{red}{1}$$
The power series may be an approximation, but the derivative is exact.