No dice, that's still pretty "scalar" looking and messy. Is this as good as it gets maybe? It's a bit surprising to me that it doesn't do a better job here (or alternatively, that the "correct" way of doing this is so hard to figure out).

That make some sense, but I can't just get rid of the square root - it's part of the actual problem I'm solving (this feeds into a gauss newton optimzer).

Is there an easy way to tell Mathematica to just treat things as real? Ive tried using Refine and Simplify (both in the dist function and wrapped "around" the Grad call) with assumptions but it makes no difference (e.g. x [Element] Reals). I also tried putting the funky "/. Abs -> Simplify" just to get rid of the abs, but without squaring the norm and that at least got rid of all the redundant "Abs" stuff, but presumably the simplifier still doesn't know that this stuff is Real so the output is still kind of horrid:

dist[p_] :=  Norm[d d.p - p] /. Abs -> Identity
FullSimplify[Grad[dist[{x, y, z}], {x, y, z}]]
Out={(x - d d.{x, y, z} - 
  d (x + y + z - 
     3 d d.{x, y, z}) d.{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}})/Sqrt[
 x^2 + y^2 + z^2 + 
  d d.{x, y, z} (-2 (x + y + z) + 3 d d.{x, y, z})], (
 y - d d.{x, y, z} - 
  d (x + y + z - 
     3 d d.{x, y, z}) d.{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}})/Sqrt[
 x^2 + y^2 + z^2 + 
  d d.{x, y, z} (-2 (x + y + z) + 3 d d.{x, y, z})], (
 z - d d.{x, y, z} - 
  d (x + y + z - 
     3 d d.{x, y, z}) d.{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}})/Sqrt[
 x^2 + y^2 + z^2 + d d.{x, y, z} (-2 (x + y + z) + 3 d d.{x, y, z})]}

Any ideas? It just seems like I'm missing something.