$$\mathbf{p}=\frac{\rho c^2 \mathbf{Q}}{V}\sum_n \frac{\epsilon_{n_x}\epsilon_{n_y}\epsilon_{n_z}\Psi(S)\Psi(R)}{2\omega_n k_n/\omega +i(\omega_n^2/\omega-\omega)}$$