# Is there something like DensityPlot3D to visualize atomic orbitals?

I'm visualizing some hydrogen like atomic orbitals. For looking at plane slices of the probability density, the DensityPlot function works well, and with something like:

Manipulate[
DensityPlot[ psi1XYsq[u, v, z], {u, -w, w}, {v, -w, w} ,
Mesh -> False, Frame -> False, PlotPoints -> 45,
ColorFunctionScaling -> True, ColorFunction -> "SunsetColors"]
, {{w, 10}, 1, 20}
, {z, 1, 20, 1}
]


I can get a nice plot

I was hoping that there was something like a DensityPlot3D so that I could visualize these in 3D, but I don't see such a function. I was wondering how DensityPlot be simulated using other plot functions, so that the same idea could be applied to a 3D plot to construct a DensityPlot3D like function?

-
You may be interested in heatmap density 3D –  Kuba Sep 15 '13 at 16:33
–  Kuba Sep 15 '13 at 16:34
Image3D works, apart from interpolation. A question arises, though: how you would be able to interpret the results? In general, isosurfaces are much more practical than images which look like fuzzy, vague clouds at best, opaque mass with hidden internal structure at worst. –  kirma Sep 15 '13 at 16:47
@kirma "fuzzy, vague clouds" - sounds like an orbital ;) –  Kuba Sep 15 '13 at 16:49
@Kuba, yes, that 4D visualization answer by halmir looks like what I was looking for. –  Peeter Joot Sep 18 '13 at 13:19

My preferred method for this kind of thing is projecting each dimension onto a plane and then combining them together. I think MATLAB has similar functionality. Mind you, the answers and comments on my question about projecting are right in pointing out that this will become inefficient for high polygon counts (essentially more PlotPoints) so if you want to Manipulate in a smooth way, you may want to use Texture. See this relevant question for details on that.

Now, I haven't got your function psi1XYsq so I will pick a double gaussian:

E^(-(x^2 + (y - 2)^2 + (z - 3)^2)/10) + E^(-((x + 1)^2 + (y + 2)^2 + (z + 2)^2)/25)


The idea is to use something similar to the function @Jens is using in his answer to How to make a drop-shadow for Graphics3D objects.

Block[{d1, d2, d3, function, options, opacity, x0, y0, z0, min, max},
{min, max} = {-9, 9};
{x0, y0, z0} = {0., 0., 0.};
opacity = 0.9;
function[x_, y_, z_] :=
E^(-(x^2 + (y - 2)^2 + (z - 3)^2)/10) +
E^(-((x + 1)^2 + (y + 2)^2 + (z + 2)^2)/25);
options = Sequence @@ {PlotPoints -> 45, Mesh -> None,
ColorFunctionScaling -> False, ColorFunction -> "SunsetColors"};
d1 = First@DensityPlot[function[x, y, z0], {x, min, max}, {y, min, max},
Evaluate@options] /. {x_?AtomQ, y_?AtomQ} -> {x, y, z0};
d2 = First@DensityPlot[function[x, y0, z], {x, min, max}, {z, min, max},
Evaluate@options] /. {x_?AtomQ, z_?AtomQ} -> {x, y0, z};
d3 = First@DensityPlot[function[x0, y, z], {y, min, max}, {z, min, max},
Evaluate@options] /. {y_?AtomQ, z_?AtomQ} -> {x0, y, z};
Show[Graphics3D[{d1, d2, d3}, Lighting -> "Neutral"] /.
GraphicsComplex[xx__] -> {Opacity[opacity], GraphicsComplex[xx]}]
]


The code is straightforward: min, max define the range for each variable, {x0, y0, z0} define the projection planes, and opacity the Opacity. You will notice I have turned off ColorFunctionScaling so that each slice is bright according to an absolute value and they merge together nicely. If your function is not normalised you may want to normalise it before doing that.

If you can afford lowering the PlotPoints, Manipulate isn't too bad, and you can make animations that look like volumetric rendering (apologies for the 300K gif):

Table[Block[{d1, d2, d3, function, options, opacity, x0, y0, z0, min,
max},
{min, max} = {-9, 9};
{x0, y0, z0} = {0., t, 0.};
opacity = 0.9;
function[x_, y_, z_] :=
E^(-(x^2 + (y - 2)^2 + (z - 3)^2)/10) +
E^(-((x + 1)^2 + (y + 2)^2 + (z + 2)^2)/45);
options =
Sequence @@ {PlotPoints -> 25, Mesh -> None,
ColorFunctionScaling -> False,
ColorFunction -> "SunsetColors"};
d1 = First@
DensityPlot[function[x, y, z0], {x, min, max}, {y, min, max},
Evaluate@options] /. {x_?AtomQ, y_?AtomQ} -> {x, y, z0};
d2 = First@
DensityPlot[function[x, y0, z], {x, min, max}, {z, min, max},
Evaluate@options] /. {x_?AtomQ, z_?AtomQ} -> {x, y0, z};
d3 = First@
DensityPlot[function[x0, y, z], {y, min, max}, {z, min, max},
Evaluate@options] /. {y_?AtomQ, z_?AtomQ} -> {x0, y, z};
Show[Graphics3D[{d1, d2, d3}, Lighting -> "Neutral"] /.
GraphicsComplex[xx__] -> {Opacity[opacity],
GraphicsComplex[xx]}]
], {t, 9, -9, -1.5}];
Export["3d.gif", %];


-
What does your {x_?AtomQ, y_?AtomQ} -> {x, y, z0} replacement in this code do? –  Peeter Joot Sep 17 '13 at 22:17
I adds a third dimension to every pair of 2d coordinates of the density plot. –  gpap Sep 17 '13 at 23:46

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

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]
]


it's quite smooth to rotate the Image3D object

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]
]


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]


-