Recently an old subroutine to calculate spherical harmonics, harmonicsr, was
replaced in the hybrid functionals code by the more commonly used ylm4
subroutine. It turns out that there is a case in which these subroutines yield
differing results:

The theta and phi angles in the routines are calculated from a vector that is
passed to them. But at some places this vector has a length of 0.0. In ylm4
this leads to the assumption that the vector points in z direction. This leads
to all m=0 Y_lm(theta,phi) having finite values. In the old subroutine a
vanishing vector length lead to an averaging over all angles. For l=0, m=0
this is the same as in ylm4. But for all other l,m combinations the
Y_lm(theta,phi) values vanish.

In most cases these different behaviors are irrelevant. But there seem to
be places where this leads to non-negligible differences in out-file values.