The situation is a bit complex. The diagonal elements of the d/dk=r operator for Bloch wavefunctions can always be made zero thanks to a k-dependent phase factor (or a set of unitary transformations among degenerate wavefunctions), for every localized region around some k. However, this is not possible when one is interested in the computation of the static electric polarization, which must be done using the Berry phase approach, as the k-dependent phase factors cannot always be modified globally in the whole Brillouin Zone to have a vanishing diagonal component contribution, as one must respect the k-point periodicity. By contrast, for the linear optical reponse, it is clear that only the non-diagonal matrix elements contribute to the final value. For the non-linear optical response at finite frequency, this is not discussed in the Sharma et al papers. I have not followed the litterature recently, while I know some people have continued to work on this topics, for the MBPT treatment of non-linear optic. Anyhow, in the ABINIT implementation, one sticks to the formulation published by Sharma et al., where the diagonal matrix elements are neglected. Whether this is an approximation, or an exact result at that level of the formalism is not clear to me.