Running HF with anisotropic contracted basis set

Problem statement:
For orbitals of angular momentum l > 0, standard quantum chem software takes the contraction for a given orbital and creates 2l+1 identical copies. For example, for p- or d-orbitals, it creates 3 and 5 copies of a given p- and d-contraction specified in the basis set input, corresponding to the different spherical harmonics (px, py, pz, or dxz, dxy, dyz, dz2, dx2-y2). Each of these copies has the same set of contraction coefficients, e.g. the contraction coefficients for px are the same as for py and pz.

My question:
Would it be possible to use or modify Psi4 (and if so, how?) to do the following:
Specify a different contraction for every spherical harmonic, i.e. different coefficients for px than for py and pz (but same basis functions i.e. same exponents). Then, when running the calculation (say, a Hartree-Fock calculation), the program should not triplicate the contraction for a p-orbital, but instead, build molecular orbitals with each contraction (px, py, pz) appearing only once in the linear combination of atomic orbitals.

I don’t believe this is possible without rewriting much of our integral and AO code.