orbitals.mws

Orbitals

(c)  David A. Harrington, Chemistry Department, University of Victoria, B.C., Canada,

http://surface.chem.uvic.ca, email:dharr@uvic.ca

The Orbitals  package evaluates, plots and calculates overlap integrals of hydrogenic or Slater-type orbitals.

Documentation and references are given in the source code file "orbitals.txt" .

>    restart: with(plots): with(plottools): interface(showassumed=0):

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

Set libname  to the directory where you have saved the Orbitals  package files maple.ind and maple.lib.

>    libname := "C:/mylib/Orbitals",libname:

>    with( Orbitals );

[HYD, HYDr, S, STO, STOr, Y, cartesian, fullcartesian, overlap, realY]

d[z^2]  orbital has an angular part which is the spherical harmonic Y(2,0,theta,phi) .

The function cartesian  converts this to the usual form in terms of x,y,z and r (use the function fullcartesian  to remove the last r)

>    dz2:=Y(2,0,theta,phi);

dz2 := 1/4*5^(1/2)*(3*cos(theta)^2-1)/Pi^(1/2)

>    cartesian(dz2,r,theta,phi,x,y,z);

-1/4*5^(1/2)*(-2*z^2+x^2+y^2)/r^2/Pi^(1/2)

The d[xz]  and d[yz] orbitals are real linear combinations of two (complex) spherical harmonics.

realY(l,m,theta,phi)  gives (Y(l,m,theta,phi)+Y(l,m,theta,phi))/sqrt(2)

realY(l,-m,theta,phi)  gives (Y(l,m,theta,phi)-Y(l,m,theta,phi))/(I*sqrt(2))

A complete list of all orbitals up to f is given at the end of this worksheet

>    dxz:=realY(2,1,theta,phi):
cartesian(dxz,r,theta,phi,x,y,z);

dyz:=realY(2,-1,theta,phi):
cartesian(dyz,r,theta,phi,x,y,z);

1/2*15^(1/2)/Pi^(1/2)*z/r^2*x

1/2*15^(1/2)/Pi^(1/2)*z/r^2*y

The plots of orbitals usually seen are just plots of the squares of their angular parts (for contour plots of the wavefunction with both radial and angular parts, see below). Remember that Maple's ( theta, phi ) is ( phi, theta ) in quantum chemistry.

>    orb:=plot3d(dz2^2,phi=0..2*Pi,theta=0..Pi,
coords=spherical,style=patchnogrid,scaling=constrained,grid=[50,50]):
display(orb, lightmodel=light1);

[Maple Plot]

Now animate it. Click on the plot and then on the play button in the toolbar.

>    viewarea:=[-0.4..0.4,-0.4..0.4,-0.4..0.4]:
orbframes:=[seq(display(rotate(orb,0,j*15/180*Pi,0),view=viewarea),j=0..11)]:
display(orbframes,insequence=true,view=viewarea, lightmodel=light1, orientation=[0,-100]);

[Maple Plot]

Now put in the radial part as well. The radial part of a 3p hydrogen orbital (atomic units) is given by the following.

(Note: arguments are n,l, orbital exponent = Zeff/n (= 1/n for hydrogen), radial co-ordinate)

>    threep:=HYDr(3,1,1/3,r);

threep := 1/81*2^(1/2)*3^(1/2)*r*exp(-1/3*r)*(4-2/3*r)

We can do Slater type orbitals (STO) as well

>    threepSTO:=STOr(3,1,1/3,r);

threepSTO := 2/1215*2^(1/2)*15^(1/2)*r^2*exp(-1/3*r)

>    plot([threep,threepSTO],r=0..20);

[Maple Plot]

If a[0]  is the Bohr radius, then in real units, the hydrogenic orbital is

>    HYDr(3,1,1/3/a[0],r);

1/81*2^(1/2)*3^(1/2)*(1/a[0])^(3/2)/a[0]*r*exp(-1/3*1/a[0]*r)*(4-2/3*1/a[0]*r)

Complete 3pz orbital with radial and angular components. Do it explicitly or use the built-in function with arguments n, l, m, exp, r, theta, phi

>    HYDr(3,1,1/3,r)*Y(1,0,theta,phi);
threepz:=HYD(3,1,0,1/3,r,theta,phi);

1/54*2^(1/2)*r*exp(-1/3*r)*(4-2/3*r)/Pi^(1/2)*cos(theta)

threepz := 1/54*2^(1/2)*r*exp(-1/3*r)*(4-2/3*r)/Pi^(1/2)*cos(theta)

Plot of a contour surface of the square of the complete 3pz wavefn, i.e. a probability surface,

though it is difficult to say what probability it encloses.

>    implicitplot3d(threepz^2=0.0001,r=0..15,phi=0..2*Pi,theta=0..Pi,
coords=spherical,style=patchnogrid,scaling=constrained,grid=[30,30,30],orientation=[0,-100], lightmodel=light1);

[Maple Plot]

Maybe it looks better done in cartesian co-ordinates?

(neither of them can be animated using rotate )

>    fullcartesian(threepz^2,r,theta,phi,'x','y','z'):
orb3pz:=implicitplot3d(%=0.0001,x=-15..15,y=-15..15,z=-15..15,style=patchnogrid,scaling=constrained,grid=[40,40,40],orientation=[0,-100], lightmodel=light1):
display(orb3pz);

[Maple Plot]

The general formula for the hydrogen orbitals in real units

>    psi(n,l,m,r,theta,phi)=HYD(n,l,m,1/n/a[0],r,theta,phi);
a[0]=epsilon[0]*h^2/Pi/m[e]/e^2;

psi(n,l,m,r,theta,phi) = ((n-l-1)!/n/(n+l)!^3)^(1/2)*(1/(n*a[0]))^(3/2)*(1/n/a[0]*r)^l*exp(-1/n/a[0]*r)*(n+l)!*Sum((-1)^i*(n+l)!/(n-l-1-i)!/(2*l+1+i)!/i!*(2/n/a[0]*r)^i,i = 0 .. n-l-1)*((2*l+1)*(l-abs(...
psi(n,l,m,r,theta,phi) = ((n-l-1)!/n/(n+l)!^3)^(1/2)*(1/(n*a[0]))^(3/2)*(1/n/a[0]*r)^l*exp(-1/n/a[0]*r)*(n+l)!*Sum((-1)^i*(n+l)!/(n-l-1-i)!/(2*l+1+i)!/i!*(2/n/a[0]*r)^i,i = 0 .. n-l-1)*((2*l+1)*(l-abs(...

a[0] = epsilon[0]*h^2/Pi/m[e]/e^2

Overlap integrals. First the general formula. STOr  means use the Slater Type Orbital radial function (could also use HYDr )

n1,l1 and zeta1 are the principal QN, ang. mmtm QN and orbital exponent of orbital 1; likewise for orbital 2;

m is the type of overlap (=0 for sigma, =1 for pi, =2 for delta) = m QN of the two orbitals, R is the distance between nuclei.

Atomic units are used. Prevent all the details of the STOr  being shown by using delayed evaluation of STOr

>    overlap(''STOr'',n1,l1,zeta1,n2,l2,zeta2,m,R);

Int(Int(1/4*('STOr')(n1,l1,rho*(tau+1)/R,1/2*(mu+nu)*R)*('STOr')(n2,l2,rho*(1-tau)/R,1/2*(mu-nu)*R)*((2*l1+1)*(l1-abs(m))!*(l1+abs(m))!)^(1/2)/(2^(l1+1))/l1!*(1-(mu*nu+1)^2/(mu+nu)^2)^(1/2*abs(m))*Sum(...
Int(Int(1/4*('STOr')(n1,l1,rho*(tau+1)/R,1/2*(mu+nu)*R)*('STOr')(n2,l2,rho*(1-tau)/R,1/2*(mu-nu)*R)*((2*l1+1)*(l1-abs(m))!*(l1+abs(m))!)^(1/2)/(2^(l1+1))/l1!*(1-(mu*nu+1)^2/(mu+nu)^2)^(1/2*abs(m))*Sum(...
Int(Int(1/4*('STOr')(n1,l1,rho*(tau+1)/R,1/2*(mu+nu)*R)*('STOr')(n2,l2,rho*(1-tau)/R,1/2*(mu-nu)*R)*((2*l1+1)*(l1-abs(m))!*(l1+abs(m))!)^(1/2)/(2^(l1+1))/l1!*(1-(mu*nu+1)^2/(mu+nu)^2)^(1/2*abs(m))*Sum(...

Try a numerical one. Two 2pz orbitals along the z axis 1 a.u. apart

>    overlap(STOr,2,1,1/2,2,1,1/2,0,1);

367/240*exp(-1/2)

>    overlap(STOr,2,1,0.5,2,1,1/2,0,1);

.9274864669

Here is the overlap as a function of distance. We must tell Maple the sign of R

>    assume( R > 0 );
olap:=overlap(STOr,2,1,1/2,2,1,1/2,0,R);

olap := -1/240*exp(-1/2*R)*(R^4-240+4*R^3-12*R^2-120*R)

>    plot(olap, R=0..20);

[Maple Plot]

Here are the angular parts of all the real orbitals up to the f orbitals.

Note that there are two common sets of real f orbitals in use, depending on which linear combinations of the Lz eigenfunctions are chosen. The other set (not shown) is more useful for situations with cubic symmetry (see F.A. Cotton, Chemical Applications of Group Theory, 3rd ed., Wiley, NY, 1990, Appendix IV).

>    s:=Y(0,0,theta,phi):'s'=cartesian(%,r,theta,phi,x,y,z);
pz:=Y(1,0,theta,phi):'p'[z]=cartesian(%,r,theta,phi,x,y,z);
px:=realY(1,1,theta,phi):'p'[x]=cartesian(%,r,theta,phi,x,y,z);
py:=realY(1,-1,theta,phi):'p'[y]=cartesian(%,r,theta,phi,x,y,z);
dz2:=Y(2,0,theta,phi):'d'[z^2]=cartesian(%,r,theta,phi,x,y,z);
dxz:=realY(2,1,theta,phi):'d'[xz]=cartesian(%,r,theta,phi,x,y,z);
dyz:=realY(2,-1,theta,phi):'d'[yz]=cartesian(%,r,theta,phi,x,y,z);
dxy:=realY(2,-2,theta,phi):'d'[xy]=cartesian(%,r,theta,phi,x,y,z);
dx2y2:=realY(2,2,theta,phi):'d'[x^2-y^2]=cartesian(%,r,theta,phi,x,y,z);
fz3:=Y(3,0,theta,phi):'f'[z^3]=cartesian(%,r,theta,phi,x,y,z);
fxz2:=realY(3,1,theta,phi):'f'[xz^2]=cartesian(%,r,theta,phi,x,y,z);
fyz2:=realY(3,-1,theta,phi):'f'[yz^2]=cartesian(%,r,theta,phi,x,y,z);
fxyz:=realY(3,-2,theta,phi):'f'[xyz]=cartesian(%,r,theta,phi,x,y,z);
fzx2y2:=realY(3,2,theta,phi):'f'[z*(x^2-y^2)]=cartesian(%,r,theta,phi,x,y,z);
fxx23y2:=realY(3,3,theta,phi):'f'[x*(x^2-3*y^2)]=cartesian(%,r,theta,phi,x,y,z);
fy3x2y2:=realY(3,-3,theta,phi):'f'[y*(3*x^2-y^2)]=cartesian(%,r,theta,phi,x,y,z);

s = 1/(2*Pi^(1/2))

p[z] = 1/2*3^(1/2)/Pi^(1/2)*z/r

p[x] = 1/2*3^(1/2)/Pi^(1/2)*x/r

p[y] = 1/2*3^(1/2)/Pi^(1/2)*y/r

d[z^2] = -1/4*5^(1/2)*(-2*z^2+x^2+y^2)/r^2/Pi^(1/2)

d[xz] = 1/2*15^(1/2)/Pi^(1/2)*z/r^2*x

d[yz] = 1/2*15^(1/2)/Pi^(1/2)*z/r^2*y

d[xy] = 1/2*15^(1/2)/Pi^(1/2)*y/r^2*x

d[x^2-y^2] = -1/4*15^(1/2)/r^2/Pi^(1/2)*(-x^2+y^2)

f[z^3] = -1/4*7^(1/2)/Pi^(1/2)*z/r^3*(-2*z^2+3*x^2+3*y^2)

f[xz^2] = -1/8*21^(1/2)*(-4*z^2+x^2+y^2)/r^3*x*2^(1/2)/Pi^(1/2)

f[yz^2] = -1/8*21^(1/2)*(-4*z^2+x^2+y^2)/r^3*y*2^(1/2)/Pi^(1/2)

f[xyz] = 1/2*105^(1/2)/Pi^(1/2)*z/r^3*y*x

f[z*(x^2-y^2)] = -1/4*105^(1/2)/r^3/Pi^(1/2)*z*(-x^2+y^2)

f[x*(x^2-3*y^2)] = -1/8*35^(1/2)/r^3*x*(-x^2+3*y^2)*2^(1/2)/Pi^(1/2)

f[y*(3*x^2-y^2)] = -1/8*35^(1/2)/r^3*y*(-3*x^2+y^2)*2^(1/2)/Pi^(1/2)

>   

>   

Disclaimer:  While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material.