Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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 hybrid orbital sample plot in x y plane

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?

share|improve this question
3  
You may be interested in heatmap density 3D –  Kuba Sep 15 '13 at 16:33
2  
Or how to visualize 4D function –  Kuba Sep 15 '13 at 16:34
2  
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
3  
@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

2 Answers 2

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

3d image

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", %];

3dGif

share|improve this answer
    
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.

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

it's quite smooth to rotate the Image3D object

enter image description here


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

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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.