[Wannier] How to write WF basis set as a combination of natural pure orbitals (s, p, d, etc.)

[Shima Sharifi]

Dear Wannier90 users,
I am using Wannier90 code to do some calculations, including spin-orbit coupling. I noticed that after Wannierisation process, the code only prints out the final spread of each Wannier function and their centers. 

I also have extracted _hr.dat file for the metric elements.
In addition, I used "wannier_plot_spinor_mode= true" and "wannier_plot_spinor_phase= up/down" and plotted MLWFs in Xcrysden format. They look like p orbitals with some arbitrary orientation in the unit cell (my initial guesses were p orbitals centered on different atoms, as well).

As the next step, I am interested to know what is the exact orbital character of my WFs in terms of natural p orbitals.
I mean, I'd like to have them like this:
|WF> = a |pz↑> + b |pz↓> + c |px↑> +d |px↓> + e |py↑> + f |py↓>   (where px, py and pz are typical p orbitals. ↑ and ↓ are spin up and down. a,b, c , ...,f are the coefficients which show their contribution in constructing that WF). 
In other words, I am looking for coefficients a, ..., f.

I'd appreciate it if you could guide me if there is any way to extract exact final WF basis set, like above mentioned one.

Thank you in advance!


Regards,
Shima