HEALpy - Positive AND negative m spherical harmonics as pixel maps

Question

I am interested in storing spherical harmonic functions Y_lm as maps in pixels, which I can then plot with Healpy. This question has already been asked and resolved here: HEALpy - Getting spherical harmonics from as function of pixel, l and m, but ONLY for positive m (Healpy convention). However, I need the complete set of spherical harmonics, but Healpy does not provide the spherical harmonics for negative m. I could really use some help on how to get these - for example, once I have the positive m spherical harmonics (as explained how to get in the link above), how can I get the negative m ones from these as well. Thanks.

I tried using scipy.special.sph_harm:

def spher_harm_test(l,m,nside):
    npix = hp.nside2npix(nside)
    sky_dec, sky_ra = calc_skymap_in_radec(nside)
    theta = np.pi/2 - sky_dec
    phi = sky_ra
    s_map=spec.sph_harm(m,l,phi,theta)
    return s_map

This seems to yield the right answer for positive m, but not for negative m (which I can see from healpy.mollview(spher_harm_test(l,m,nside)), for some specific choice of parameters - for example, the maps for l=4, m=4 and l=4, m=-4 look identical). But even if the spherical harmonics seemed right in mollweide projection for both positive and negative m, I am a bit sceptical about this method because I do not know if Healpy uses the exact same definition and computation for spherical harmonics as scipy...

Another thing I tried was to define a function spher_harm_posm for the positive m spherical harmonics, and another function spher_harm_negm that gets the negative m spherical harmonics from the positive m ones after applying a rotation and an interpolation. However, this only seems to work well for very high nside; for very low nside, it fails quite miserably as can be seen when plotting with healpy.mollview... See below:

def spher_harm(l_max, nside):
    npix = hp.nside2npix(nside)
    alm = np.zeros(hp.Alm.getidx(l_max, l_max, l_max) + 1, dtype=np.complex128)
    sph_harm_lmax_real=np.zeros((len(alm),npix))
    sph_harm_lmax_imag=np.zeros((len(alm),npix), dtype='complex')
    for count in range(len(alm)):
        alm[count]=0.5
        sph_harm_lmax_real[count,:]=hp.alm2map(alm, nside=nside)
        alm[count]=-0.5j
        sph_harm_lmax_imag[count,:]=hp.alm2map(alm, nside=nside)*1j
        alm[count]=0
    return sph_harm_lmax_real + sph_harm_lmax_imag   #(i,p) array, where i gives the (l,m) pair

def spher_harm_negm(l_max, nside):
    posm=spher_harm(l_max, nside)
    dim1,dim2=posm.shape
    rot_map=posm
    for i in range(dim1):
        l, m=hp.Alm.getlm(l_max,i)
        if m!=0:
            r = hp.Rotator(rot=[90/m,0])
            t,p = hp.pix2ang(nside, np.arange(hp.nside2npix(nside)))
            trot, prot = r(t,p)
            rot_map[i,:]=hp.get_interp_val(posm[i,:], trot, prot)
    negm=rot_map
    for i in range(l_max+1):
        negm=np.delete(negm,0,0)
    dim1,dim2=negm.shape
    for i in range(dim2):
        negm[:,i]=np.flip(negm[:,i])
    return negm

So both of these attempts failed, I am open to suggestions. Thank you in advance.