Why are you limiting it to a particular potential? If you are using forces
to find the phonon dispersion curves, you might as well use a lattice
dynamics program such as GULP.
I have thought about this, and the MD way is to just use the atom motions
since phonons are just wave decompositions of aggregate motions.
In LAMMPS this could be generalized to any region described by a lattice,
somehow mapping a 'type' to the atoms for unit cells with basis of more than
one atom.
Also, I'm not sure what the benefit is to have LAMMPS do it since you have
to have a time series of the motions. So it would be best to do it as a
post-process.