Your unmodified run.in does not work with fixed boundaries (boundary f f f) instead of shrink-wrapped boundaries (boundary s s s), a classic subtle hint of numerical instability. Looking at the following settings:
region nucleus sphere 500.0 500.0 500.0 15.0 side in
fix wall1 all wall/region nucleus lj126 1.0 1.0 1.12
fix wall2 te spring tether 50.0 500.0 500.0 500.0 14.0
fix wall3 ce spring tether 50.0 500.0 500.0 500.0 14.0
a particle always experiences very intense forces at a radial distance between 14 and 15 units away from the sphere center, resulting in numerical instability.
My guess is that the instability causes particles to “freeze” on the sphere (since any time they try to move transverse to the sphere they end up experiencing very high forces and strong damping from the thermostat). They don’t so much form clusters as retain their initial clustering that they can’t diffuse away from.
More importantly, removing the wall/region fix (and retaining the tether fixes) results in a simulation that runs with fixed boundaries and shows particles diffusing along the surface of the sphere as desired.