In an [presentation][1] by Markus van Almsick, he gives an solution to visualize atomic orbitals using Image3D.

Radius wave function (hydrogen):

    R[n_Integer?Positive, l_Integer?NonNegative, r_] :=
     Block[{ρ = (2 r)/n}, 
       Sqrt[(2/n)^3 (n - l - 1)!/(2 n (n + l)!)] E^(-ρ/2) ρ^l LaguerreL[n - l - 1, 2 l + 1, ρ]] /; l < n

full wave function:

    ψ[n_, l_, m_, r_, ϑ_, φ_] := ψ[n, l, m, r, ϑ, φ] = 
      FullSimplify[R[n, l, r] SphericalHarmonicY[l, m, ϑ, φ], {r >= 0, ϑ ∈ Reals, φ ∈ Reals}]
    

    CompileWaveFunction = Compile[{{x, _Real}, {y, _Real}, {z, _Real}},
       Block[{ρ = x^2 + y^2, r, ϑ, φ},
        If[ρ > 0,
         r = Sqrt[ρ + z^2]; ϑ = ArcCos[z/r]; φ = ArcTan[x, y],
         r = Abs[z]; ϑ = π/2 Sign[z]; φ = 0];
         #
        ],
       CompilationTarget -> "C"
       ] &;

color function:

    colorFunction = (Blend[{
          {0., RGBColor[0.7, 0.8, 1., 0.]}, 
          {0.1, RGBColor[0., 0.7, 0.1, 0.012]}, 
          {0.4, RGBColor[1., 0.1, 0.03169, 0.06723]}, 
          {1., RGBColor[1., 0.95051, 0., 0.10963]}}, #] &)

plot 3p orbital

    Block[
     {nψ = 
       CompileWaveFunction[ψ[3, 1, 0, r, ϑ, φ]], data, vol},
     data = Table[Abs[nψ[x, y, z]]^2, {z, -20, 20, 0.25}, {y, -20, 20, 0.25}, {x, -20, 20, 0.25}];
     vol = RawArray["Byte", Round[(255/Max[data]) data]];
     Image3D[vol, "Byte", Background -> Black, 
      Method -> {"FastRendering" -> True, "InterpolateValues" -> True}, 
      ColorFunction -> colorFunction, BoxRatios -> 1]
     ]

![enter image description here][2]

it's quite smooth to rotate the Image3D object

![enter image description here][3]



----


We can also visualize the atomic orbital by plotting the isosurface:

    Block[{nψ = CompileWaveFunction[ψ[3, 2, 0, r, ϑ, φ]]},
     ContourPlot3D[Abs[nψ[x, y, z]]^2, {x, -20, 20}, {y, -20, 20}, {z, -20, 20}, 
      PlotPoints -> 15, Contours -> {0.00002}, 
      ColorFunction -> Function[{x, y, z}, ColorData["Rainbow"]@Rescale[Arg[nψ[x, y, z]], {-π, π}]], 
      ColorFunctionScaling -> False, Mesh -> None]
     ]

![enter image description here][4]


and make a animation showing different isosurfaces

    plots = ParallelTable[
       Block[{nψ = CompileWaveFunction[ψ[4, 2, 1, r, ϑ, φ]]}, 
         ContourPlot3D[
          Abs[nψ[x, y, z]]^2, {x, -20, 20}, {y, -20, 20}, {z, -20, 20}, PlotPoints -> 17, Contours -> {ct}, 
          ColorFunction -> Function[{x, y, z}, ColorData["Rainbow"]@Rescale[Arg[nψ[x, y, z]], {-π, π}]], 
          Boxed -> False, Axes -> False, ColorFunctionScaling -> False, Mesh -> None, 
          ViewPoint -> {0.98, -2.76, 1.7}, ViewVertical -> {-0.004, -0.117, 0.993}]],
       {ct, 0.00003, 0.000015, -0.0000005}
       ];
    ListAnimate[plots]

![enter image description here][5]


  [1]: http://www.wolfram.com/broadcast/video.php?channel=105/?fp=left&video=1554
  [2]: http://i.stack.imgur.com/b0GPf.png
  [3]: http://i.stack.imgur.com/WIx00.gif
  [4]: http://i.stack.imgur.com/KJIrd.png
  [5]: http://i.stack.imgur.com/TACP6.gif