For the zone-center (Gamma point, i.e. in the long wavelength limit) modes in polar materials, it goes like this: First, determine the standard, regular dynamical matrix (which might be obtained by very different methods) that doesn't contain the polarization effects; its eigenvalues do not contain LO modes. Additionally, determine the Born effective charges and the dielectric constants. From these, the non-analytic part of the dynamical matrix is obtained. Diagonalization of the full dynamical matrix for a chosen wave vector direction yields the full spectrum of the phonon frequencies including those of the LO modes. Further details (and further references) you can find, e.g., in my Ph.D. dissertation (https://www.researchgate.net/publication/283271282), subsection 2.3.3.
Thesis Structure and Lattice Dynamics of GaN and AlN: Ab-Initio Inv...