We consider the classical problem to find optimal distributions of interacting particles on a sphere by solving an evolution problem for a particle density. The higher order surface partial differential equation is an approximation of a surface dynamic density functional theory. We motivate the approach phenomenologically and sketch a derivation of the model starting from an interatomic potential. Different numerical approaches are discussed to solve the evolution problem: (a) an implicit approach to describe the surface using a phase-field description, (b) a parametric finite element approach, and (c) a spectral method based on nonequispaced fast Fourier transforms on the sphere. Results for computed minimal energy configurations are discussed for various particle numbers and are compared with known rigorous asymptotic results. Furthermore extensions to other more complex and evolving surfaces are mentioned.