Yes, there is a slightly hidden function for this:
calculateDipoleTransitionMatrixElements(configuration, kpoint=None, number_of_states=None)
Calculates the dipole transition matrix elements for a bulk configuration at a specified k-point.
The x-component of the dipole matrix elements is defined as
x_nm(k) = <mk | x | nk>,
where |nk> is a Bloch state in band 'n' at the k-point 'k', and likewise for the y- and z-components. When calculating the dipole matrix element, we use the commutator relation [H, x] = -i*hbar*p_x / m, where p_x is the momentum matrix element. Note, however, that this commutator relation is only approximately true for a normconserving DFT Hamiltonian due to the non-local pseudopotential operator. Using the commutator relation, the relation between the dipole transition matrix element and the momentum matrix element is x_nm = -i * hbar * p_nm^x / (m * E_nm), where 'm' is the electron mass and E_nm is the energy difference between the band energies E_n(k) - E_m(k).
configuration = The configuration for which to calculate the gradients
kpoint = The kpoint (as three floats representing fractional reciprocal space coordinates)
param number_of_states = The states in the index interval [0, number_of_states - 1] are used to construct the momentum operator matrix elements
Returns a matrix with shape (3, N, N), where N is the number of bands included (unpolarized, noncolinear and spin-orbit), or (2, 3, N, N) with the first index corresponding to Spin.Up and Spin.Down (for (spin-polarized)