Introduction

This worksheet gives some examples of using the Orbitals Package. The Orbitals package evaluates, plots and calculates orbitals, overlap integrals, and atomic four-electron interals for hydrogenic or Slater-type orbitals. To set up the package, copy the three files Orbitals.hdb, Orbitals.lib and Orbitals.ind to a directory of your choice, such as "c:\Maplelibs". Further details are given in the file "Orbitals ReadMe.txt". Documentation and references are given in the help pages, which are accessible after the command "with(Orbitals)" in the initialization section of this worksheet. These pages are accessed using the command ?Orbitals, or are found under "Science and Engineering/Orbitals" in the help browser.

Initialization

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

Warning, the name changecoords has been redefined

Warning, the assigned name arrow now has a global binding

Set libname to include the directory where you have saved the Orbitals package files Orbitals.ind, Orbitals.lib and Orbitals.hdb to, e.g., "c:\Maplelibs".

>	libname := "C:/MapleLibs",libname: with(Orbitals);

Working with and plotting orbitals

orbital has an angular part which is the spherical harmonic .

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);

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

The and orbitals are real linear combinations of two (complex) spherical harmonics.

gives

gives

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);

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 ( ) is ( ) 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);

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

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);

We can do Slater type orbitals (STO) as well

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

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

If is the Bohr radius, then in SI units, the hydrogenic orbital is

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

Complete orbital with radial and angular components. Do it explicitly or use the built-in function with arguments

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

Plot of a contour surface of the square of the complete 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);

The general formula for the hydrogen orbitals in real units

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

Plots of rigid rotor wavefunctions

The spherical harmonics are used also in the rigid rotor wavefunctions:

>	psi:=(l,m)->exp(-I*l*(l+1)*T)*Y(l,m,theta,phi);

Plot the phase on the surface of the unit sphere, using a color code for the phase. The nodes around two lines of latitude are easily seen.
Here is the (l,m)=(3,1) case. The factor of 4 in front of psi is for convenience in colouring.

>	plot3d(1,phi=0..2*Pi,theta=0..Pi,color=argument(subs(T=0,4*psi(3,1)))/2/Pi, coords=spherical,scaling=constrained,grid=[50,50],style=patchnogrid,orientation=[0,60]);

Now do the probability density on a grayscale plot for the same wavefunction. Darker is higher probability density.

>	pdplot:=proc(l,m) local col;    col:=1-evalc(abs(4*psi(l,m)))^2: # colour    plot3d(1,phi=0..2*Pi,theta=0..Pi,           coords=spherical,scaling=constrained,           grid=[50,50],color=[col,col,col],           style=patchnogrid,orientation=[0,60]);    end proc:

>	pdplot(3,1);

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) and is just the m QN of the two orbitals, R is the distance between nuclei.

Atomic units are used. For the general formula, we can just name the radial function rho; later it can be STOr or HYDr

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

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

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

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

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);

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

Variational calculation for Helium

In the simplest variational treatment of Helium, we use a configuration with the orbital exponent zeta ( ) as a variational parameter. The Hamiltonian is , and can be rearranged to . The first term in brackets is two hydrogenlike atoms with nuclear charge and gives energy twice .

>	epsilon[`1s`]:=-zeta^2/2;

The next term has energy and the integral can be evaluated straightforwardly:

>

`1s(1)`:=HYD(1,0,0,zeta,r1,theta1,phi1); `1s(2)`:=HYD(1,0,0,zeta,r2,theta2,phi2); Int(Int(Int(Int(Int(Int(`1s(1)`*`1s(2)`*(zeta-2)/r1*`1s(1)`*`1s(2)`*r1^2*sin(theta1)*r2^2*sin(theta2),phi1=0..2*Pi),theta1=0..Pi),r1=0..infinity),phi2=0..2*Pi),theta2=0..Pi),r2=0..infinity); E2:=value(%) assuming positive;

The next term has the same value and the last term is the coulomb integral J1s1s

>	J1s1s:=atomJ(HYDr,1,0,0,zeta,1,0,0,zeta);

We find the total energy and minimize it with respect to zeta. We find zeta = 27/16 and the energy is -(27/16)^2= -729/256 = -2.847656250 a.u.; experimental energy is -2.9033

>	E:=2*epsilon[`1s`]+2*E2+J1s1s; diff(E,zeta); solve(%,zeta); subs(zeta=%,E); evalf(%);

List of angular parts of orbitals

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);

>

References

References used in developing the package functions are given in the help page ?Orbitals,references.

The variational calculation for Helium is given in many quantum texts, e.g., D.A. McQuarrie, Quantum Chemistry, University Science Books, 1983, or I.R. Levine, Quantum Chemistry, 2nd Ed., Allyn and Bacon, 1974.

Conclusions

The Orbitals package is useful for plotting and manipulating atomic orbitals. Numerical results are available by specifying all parameters, but symbolic arguments gnerally lead to general formulas. If one or two parameters are left symbolic, plots readily show the variation with the parameter(s), e.g., overlap integrals as a function of internuclear distance.