Atomic Models

(* Rutherford model *)

es[t_] := {{Cos[t + Pi], 2.5*Sin[t + Pi], 2*Sin[t + Pi]}, 
           {Cos[t + (4*Pi)/5], 2*Sin[t + (4*Pi)/5], -1.5*Sin[t + (4*Pi)/5]}, 
           {2*Sin[t + (3*Pi)/5], Cos[t + (3*Pi)/5], 2*Sin[t + (3*Pi)/5]}, 
           {2.5*Sin[t + (2*Pi)/5], Cos[t + (2*Pi)/5], -1.5*Sin[t + (2*Pi)/5]}}
Manipulate[Show[
    ParametricPlot3D[Evaluate[es[2*u]], {u, 0, 2*Pi}, PlotStyle -> Directive[Thick, Dotted]], 
    Graphics3D[
           {Specularity[White, 200], 
            Flatten[Table[{ColorData[97, "ColorList"][[i]], Sphere[es[2*t][[i]], 0.2]}, {i, 1, 4}]], 
            Orange, Sphere[(1/5)*PolyhedronData["Icosahedron", "VertexCoordinates"], 1/6]}
    ]], {t, 0, 2*Pi, 0.001}]


(* QM model *)

a = 1; 
wf[n_, \[ScriptL]_, m_, r_, \[Theta]_, \[Phi]_] := (Sqrt[(2/(n*a))^3*((n - \[ScriptL] - 1)!/(2*n*(n + \[ScriptL])!^3))]*
    ((2*r)/(n*a))^\[ScriptL]*LaguerreL[n - \[ScriptL] - 1, 2*\[ScriptL] + 1, (2*r)/(n*a)]*
    SphericalHarmonicY[\[ScriptL], m, \[Theta], \[Phi]])/E^(r/(n*a))
pd[n_, \[ScriptL]_, m_, r_, \[Theta]_, \[Phi]_] := FullSimplify[Conjugate[wf[n, \[ScriptL], m, r, \[Theta], \[Phi]]]*wf[n, \[ScriptL], m, r, \[Theta], \[Phi]]]

SphericalPlot3D[
    pd[5, 3, 0, r, \[Theta], \[Phi]] /. r -> a,
    {\[Theta], 0, Pi}, {\[Phi], 0, 2*Pi}, 
    PlotPoints -> 20, PlotRange -> All]