HYDROGEN WAVE FUNCTION PROBABILITY DENSITIES by Corinne Manogue, Kerry Browne and Steve Sahyun Copyright 2010 Corinne Manogue This worksheet is intended to be used by the teacher for display purposes. Three useful plotting functions are defined int eh section below. In order to begin using this you must first open the section by clicking on the triangle on the left hand side below. Then you must run the command inside in order to load the functions. These functions seem to work fine for linear combinations of hydrogen orbitals. The scaling factors are difficult to calculate so I have opted for a simplified scheme. As a result you may need to play a bit with the scaling factors to get things to work. If you run into any problems or need help with this, please feel free to e-mail me at kerry.p.browne@gmail.com.
<Text-field style="Heading 1" layout="Heading 1"><Font size="12">The commands in this section setup the worksheet and define the function hprobplot. This function can be used to plot a linear combination of hydrogen atom orbitals. Details are listed inside this section: </Font><Font underline="true" italic="true" size="16">These commands will automatically execute when you open this file. If you experience problems with Maple running slowly you may need to either close and reopen the file or reexecute these commands manually.</Font></Text-field> hprobplot(Psiex, nmax, options) is a command for plotting hydrogen atom orbitals in a variety of ways. Psiex and nmax must be defined in order for the function to execute properly. Psiex is the input expression, usually a linear combination of hydrogen atom eigenstates. nmax is the largest value of n included in the linear combination. This value is used primarily to set the spatial range. Options: Options are defined by including optionname=option_value in the paratmeter list. For example including plottype=timeanim sets the plottype to timeanim. plottype determines what type of plot will be displayed. Plottype can take on several different values including: static1: a static plot of a single slice through the xz-plane. This type plots very quickly. static2: a static plot including two slices; one though the xz-plane and one parallel to the xy-plane with the z value set by the optional parameter zvalue. This type plots very quickly staticanim: this plot type generates an animation with a static slice through the xz-plane and a moving slice parallel to the xy-plane. The animation varies the value of the z coordinate of the slice. By dragging the slider in the animation, it is easier to get a 3-D picture of the probability density. timeanim1: a plot of a single slice through the xz-plane is generated for each value of time. These plots are then displayed as an animation. timeanim2: a plot including two slices; one though the xz-plane and one parallel to the xy-plane with the z value set by the optional parameter zvalue is generated for each value of time. These plots are then displayed as an animation in time. colorscale is a factor used to adjust the color display. If there is too much purple, decrease this number, if there is too much red, increase this value. The default value is 1 if this value is not set. rsccale is a factor used to adjust the radial range. If you want to zoom out, increase rsc. If you want to zoom in closer to the origin, decrease rsc. timevalue is the time value used for calculating the plot.The default value is zero if this value is not set. numframes sets the number of frames calculated for the animation. A higher value for numframes makes the animation smoother, but the calculation slower. The default value is 25 if this parameter is not set. maxtime adjusts the amount of time in the animation. By increasing maxtime, the animation will show a longer period of time without adjusting the number of frames. zvalue sets the vertical position of the slice parallel to the xy-plane when it is displayed. restart:with(plots):with(orthopoly):with(StringTools): numpts:=4096: setoptions3d(axes=BOXED, style=patchnogrid, numpoints=100000000, scaling=constrained): assume(r>0,theta>0,theta<Pi,phi>0,phi<2*Pi, x,real,y,real,z,real, t, real); Zq:=1: a0:=1: Laguerre:=(q)-> exp(rho)*diff(rho^q*exp(-rho),[rho\$q]): ALaguerre := (q,p) -> diff(Laguerre(q, rho), [rho\$p]): Nr:=(n,el) -> sqrt((2*Zq/(a0*n))^3*(n-el-1)!/(2*n*((n+el)!)^3))*(2*Zq/(n*a0))^el: R:=(n,el) -> Nr(n,el)*r^el*exp(-Zq*r/(n*a0))*subs(rho=2*Zq*r/(n*a0),ALaguerre(n+el, 2*el+1)): Ny:=(el,m)->(-1)^((m+abs(m))/2)*sqrt((2*el+1)*(el-abs(m))!/(4*Pi*(el+abs(m))!)): ALegendre:=(el,m)->(1-z^2)^(abs(m)/2)*diff(P(el,z),[z\$abs(m)]): Y:=(el,m)->Ny(el,m)*subs(z=cos(theta),ALegendre(el,m))*exp(I*m*phi): PsiH:=(n,l,m)-> R(n,l)*Y(l,m)*exp(-I*t/n^2): hprobplot:= proc(Psiex,nmax,{plottype:=static1,colorscale:=1,rscale:=1,timevalue:=0,numframes:=25, maxtime:=10, zvalue:=0, titlevalue:=" "}) option autocompile: local Psiexsimp, Psiexsq, fmin, fmax, scale, CPsiexsq, CPsiexsqxyz, maxr, plot1, plot2, maxfvalue, maxt, offset: if (plottype=timeanim1 or plottype=timeanim2) then Psiexsimp:=simplify(subs(t=t, Psiex)) else Psiexsimp:=simplify(subs(t=timevalue, Psiex)) fi: Psiexsq:=simplify(expand(abs(Psiexsimp)^2)): maxr:=rscale*2*nmax^2: offset:=0.005*maxr: fmin:=0: maxfvalue := [3., 300., 3000., 20000., 50000., 1.*10^5, 3.*10^5, 5.*10^5, 1.*10^6, 2.*10^6, 3.*10^6, 5.*10^6, 1.*10^7, 1.*10^7, 1.*10^7, 2.*10^7, 2.*10^7, 3.*10^7, 5.*10^7, 5.*10^7, 1.*10^8, 1.*10^8, 1.*10^8, 1.*10^8, 2.*10^8, 2.*10^8, 2.*10^8, 3.*10^8, 3.*10^8, 5.*10^8, 5.*10^8, 5.*10^8, 1.*10^9, 1.*10^9, 1.*10^9, 1.*10^9, 1.*10^9, 1.*10^9, 2.*10^9, 2.*10^9, 2.*10^9, 2.*10^9, 2.*10^9, 3.*10^9, 3.*10^9, 3.*10^9, 5.*10^9, 5.*10^9, 5.*10^9, 5.*10^9]: scale:=colorscale*maxfvalue[nmax]/2: CPsiexsq:=scale*subs(theta=arccos(z/r), r=sqrt(x^2+y^2+z^2), phi=piecewise(x>0, arctan(y/x),Pi+arctan(y/x)), Psiexsq): CPsiexsqxyz:=piecewise(CPsiexsq>0.9,0.9,CPsiexsq): if (plottype=static2) then plot1 := plot3d([x, offset, z], x = -maxr+offset .. maxr, z = -maxr+offset .. maxr, color = subs(y = offset, CPsiexsqxyz), numpoints = numpts, view = [-maxr .. maxr, -maxr .. maxr, -maxr .. maxr]): plot2 := plot3d([x, y, zvalue+offset], x = -maxr+offset .. maxr, y = -maxr+offset .. maxr, color = subs(z = zvalue+offset, CPsiexsqxyz), numpoints = numpts, view = [-maxr .. maxr, -maxr .. maxr, -maxr .. maxr], transparency = 0.3): display([plot1,plot2], orientation = [80, 65], title=titlevalue): elif (plottype=static1) then plot3d([x, offset, z], x = -maxr+offset .. maxr, z = -maxr+offset .. maxr, color = subs(y = offset, CPsiexsqxyz), numpoints = numpts, view = [-maxr .. maxr, -maxr .. maxr, -maxr .. maxr], orientation = [90,90], title=titlevalue): elif (plottype=staticanim) then plot1 := animate(plot3d, [[x, y, zval], x = -maxr+offset .. maxr, y = -maxr+offset .. maxr, color = subs(z = zval, CPsiexsqxyz), numpoints = numpts, view = [-maxr .. maxr, -maxr .. maxr, -maxr .. maxr]], zval = -maxr+offset .. maxr, frames = numframes): plot2 := plot3d([x, offset, z], x = -maxr+offset .. maxr, z = -maxr+offset .. maxr, color = subs(y = offset, CPsiexsqxyz), numpoints = numpts, view = [-maxr .. maxr, -maxr .. maxr, -maxr .. maxr]): display({plot1, plot2}, orientation = [80, 65], title=titlevalue): elif (plottype=timeanim1) then animate(plot3d,[[x, 0.1e-1, z], x = -maxr+offset .. maxr, z = -maxr+offset .. maxr, color = subs(y = offset, CPsiexsqxyz), numpoints = numpts, view = [-maxr .. maxr, -maxr .. maxr, -maxr .. maxr]], t=0..maxtime, frames=numframes, orientation=[90,90], title=titlevalue): elif (plottype=timeanim2) then plot1 := animate(plot3d,[[x, offset, z], x = -maxr+offset .. maxr, z = -maxr+offset .. maxr, color = subs(y = offset, CPsiexsqxyz), numpoints = numpts, view = [-maxr .. maxr, -maxr .. maxr, -maxr .. maxr]], t=0..maxtime, frames=numframes): plot2 := animate(plot3d,[[x, y, zvalue+offset], x = -maxr+offset .. maxr, y = -maxr+offset .. maxr, color = subs(z = zvalue+offset, CPsiexsqxyz), numpoints = numpts, view = [-maxr .. maxr, -maxr .. maxr, -maxr .. maxr], transparency = 0.3], t=0..maxtime, frames=numframes): display([plot1, plot2], orientation = [80, 65], title=titlevalue): else print("Incorrect value for plottype") fi: end proc:
Psiex:=1/sqrt(2)*(simplify(PsiH(3,0,0))+simplify(PsiH(2,1,0))); hprobplot(Psiex,3, plottype=static1,colorscale=1, rscale=1, zvalue=0); hprobplot(Psiex,3, plottype=static2,colorscale=1, rscale=0.5, zvalue=-4, maxtime=48); hprobplot(Psiex,3, plottype=staticanim,colorscale=1, rscale=1); hprobplot(Psiex,3, plottype=timeanim1,colorscale=1, rscale=1, maxtime=48); hprobplot(Psiex,3, plottype=timeanim2,colorscale=1, rscale=0.75, zvalue=-4, maxtime=48);