FullSimplify completely changes result
Floats cause that. The Rationalize
command helps.
FullSimplify[Rationalize[E^(-33.3333 (0.8 - (x^2 + y^2 + z^2)^0.26)^2) +
E^(-33.3333 (0.8 + (x^2 + y^2 + z^2)^0.26)^2), 10^-40]]
E^(-((333333 (4 - 5 (x^2 + y^2 + z^2)^(13/50))^2)/250000)) + E^(-(( 333333 (4/5 + (x^2 + y^2 + z^2)^(13/50))^2)/10000))
Addition. Maybe, the result of
FullSimplify[E^Expand@Rationalize[(-33.3333 (0.8 - (x^2 + y^2 + z^2)^0.26)^2),
10^-40] + E^Expand@Rationalize[(-33.3333 (0.8 + (x^2 + y^2 + z^2)^0.26)^2),
10^-40], Assumptions -> {x, y, z} \[Element] Reals]
E^(-((333333 (4 + 5 (x^2 + y^2 + z^2)^(13/50))^2)/ 250000)) (1 + E^((333333 (x^2 + y^2 + z^2)^(13/50))/3125))
is more preferable.
Try first to simplify, then to substitute:
Simplify[func, {\[Rho] > 0, p > 0, \[Sigma] > 0}] /. {\[Sigma] -> 3./
100, p -> 0.26, \[Rho] -> 0.8}
(* E^(-33.3333 (-0.8 + (x^2 + y^2 +
z^2)^0.26)^2) + E^(-33.3333 (0.8 + (x^2 + y^2 + z^2)^0.26)^2) *)
Have fun!