Visualization of atomic orbitals

Share icon

Hello,

Did someone ever try to Surface Plot hydrogen atomic orbitals? Or simply put, how to quicker visualize the following function in igor:

(r*r*exp(-r/3))^2*magsqr(sphericalHarmonics(2,0,theta,phee))

where r, theta, and phee are the spherical coordinates of a position in 3D space. The above expression is called the d_z2 orbital in chemistry.

The surface plot of the above expression should be similar to the 2nd graph from top in this page generated by Mathematica:
http://mathematica.stackexchange.com/questions/32378/is-there-something…

Thank you.

zh
This is at least part of the solution:
Function HOrb(xx, yy, phee)
    Variable xx, yy, phee
 
    Variable rr = sqrt(xx^2 + yy^2)
    Variable theta = atan2(xx, yy)
    return (rr*rr*exp(-rr/3))^2*magsqr(sphericalHarmonics(2,0,theta,phee))
end

But if I'm not mistaken, the spherical harmonic is rotationally symmetric, so the phee input isn't really needed...
Using these commands:
make/n=(100,100) image
setscale/I x -10,10,image
setscale/I y -10,10,image
display;appendimage image
image = HOrb(x, y, 0)
ModifyImage image ctab= {*,*,BrownViolet,0}

results in the attached image.

That's not a "surface plot", but neither was the picture you linked. I presume you actually want an isosurface plot. I'll let the Gizmo wizard weigh in of that's what you're looking for.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
HOrbital.png (29.77 KB)

May 2, 2016 at 10:06 am - Permalink

BRDF wrote:
Hello,

Did someone ever try to Surface Plot hydrogen atomic orbitals? Or simply put, how to quicker visualize the following function in igor:

(r*r*exp(-r/3))^2*magsqr(sphericalHarmonics(2,0,theta,phee))

where r, theta, and phee are the spherical coordinates of a position in 3D space. The above expression is called the d_z2 orbital in chemistry.

The surface plot of the above expression should be similar to the 2nd graph from top in this page generated by Mathematica:
http://mathematica.stackexchange.com/questions/32378/is-there-something…

Thank you.

zh


I suggest you check out the SphericalHarmonicsDemo (File Menu->Example Experiments->Visualization...)

A.G.
WaveMetrics, Inc.

May 2, 2016 at 11:44 am - Permalink

John,

Right, what I wanted is an isosurface plot. The atomic wave function (actually its magnitude squared) I wanted to visualize is given by the first attached picture (taken from the middle of the following page):
https://en.wikipedia.org/wiki/Hydrogen_atom

If visualized efficiently, the n=3, L=2, m=0 case is given by the second attached picture (the link I gave in the original post);another example is (for n=4, l=3, m=1) seen here:
https://en.wikipedia.org/wiki/Hydrogen_atom#/media/File:Hydrogen_eigens…
I am sure Igor can do this as well.

Below is my code to calculate the magnitude squared wave function. Then what's the easiest way to visualize it as the above examples show?

Thank you.
zh

constant a0=1   //5.29e-11 meters,  Bohr radius
// Hydrogen atomic wavefunction
function WF(n, L, M, xx, yy, zz) // n, L, M are the principal, angular and magnetic quantum numbers, respectively
variable n, L, M, xx, yy, zz            // transform from Cartesian to spherical coordinate
variable rr=sqrt(xx^2+yy^2+zz^2)  // radial distance of the electron from the nucleaus
variable theta= atan2(zz, rr)         // zenith angle
variable phee=atan2(xx, yy)      // azimuth angle 
variable rou=2*rr/(n*a0)     // rou is basically r, see (middle of the page) en.wikipedia.org/wiki/Hydrogen_atom
variable c= factorial(n-L-1)/(2*n)/factorial(n+L)   // a bunch of constant values
variable radial= sqrt((2/n*a0)^3*c)*exp(-rou/2)*rou^L*LaguerreA((n-L-1), (2*L+1), rou)
variable angular_squared= magsqr(SphericalHarmonics(L, M, theta, phee))
return radial*radial*angular_squared 
End
hydrogen_wave_fct.png (8.27 KB) n3_L2_m0.png (170.76 KB)

May 4, 2016 at 05:10 am - Permalink

Quote:

I suggest you check out the SphericalHarmonicsDemo (File Menu->Example Experiments->Visualization...)
A.G.
WaveMetrics, Inc.


Thanks A.G., I did but couldn't figure out how to do it when radial coordinate kicks in. Please see my above post.

I wonder why the command window history is empty for many example experiments such as SphericalHarmonicsDemo?
zh

May 4, 2016 at 05:19 am - Permalink

BRDF][quote wrote:

Thanks A.G., I did but couldn't figure out how to do it when radial coordinate kicks in.

You need to take a look at the procedure that is executed by the panel. That would lead you to the calcParametric() function. 
Function calcParametric(L,M,size)
    Variable L,M,size
    
    Make/O/N=(size+1,size+1,3) M_Parametric
    
    Variable i,j,dt,df,rr,theta,phi
    dt=pi/size
    df=2*dt
    
    SetScale/P x 0,1,"", M_Parametric
    SetScale/P y 0,1,"", M_Parametric
    
    for(i=0;i<=size;i+=1)
        theta=i*dt
        for(j=0;j<=size;j+=1)
            phi=j*df
            rr=sqrt(magsqr(sphericalHarmonics(L,m,theta,phi)))
            M_Parametric[i][j][0]=rr*sin(theta)*cos(phi)
            M_Parametric[i][j][1]=rr*sin(theta)*sin(phi)
            M_Parametric[i][j][2]=rr*cos(theta)
        endfor
    endfor
End


As you can see, I'm generating a parametric surface which is appropriate for displaying the magnitude of the spherical harmonics. I expect that your expression is for probability density and as such you should generate the data for some range of radius values and then plot an isosurface for a desired constant or volume slices in Gizmo. To do that I recommend starting from a sufficiently dense cube, e.g.,
 Make/N=(301,301,301) ddd
(this is purposely odd so that you can have symmetric range about an origin.)
Then write a function to fill up the values. Here is a rough draft which I have not tested:
Function fillProb()
 
    Make/O/N=(301,301,301) ddd  
    Variable center=150
    Variable L=1,M=1            // set this from your panel
    Variable i,j,k,xx,yy,zz,rr,rho,phi,theta
    for(k=0;k<301;k+=1)
        zz=(k-center)/center  // normalized +-1 range
        for(j=0;j<301;j+=1)
            yy=(j-center)/center
            for(i=0;i<301;i+=1)
                xx=(i-center)/center
                // compute theta phi and rr:
                rho=xx^2+yy^2
                rr=sqrt(rho+zz^2)
                phi=atan2(yy,xx)
                theta=atan2(zz,sqrt(rho))
                // the expression you posted takes the form:
                ddd[i][j][k]=(rr^2*exp(rr/3))^2*magsqr(sphericalHarmonics(L,m,theta,phi))
            endfor
        endfor
    endfor
End


I hope this helps.


Quote:
I wonder why the command window history is empty for many example experiments such as SphericalHarmonicsDemo?
zh

That's because we want to clean up the clutter left by all the attempts to tweak the display.

A.G.
WaveMetrics, Inc.

May 4, 2016 at 11:31 am - Permalink

Thanks A. G. Now I can generate the desired patterns. The attached picture is for n=3, L=2 and m=0, or the d_z^2 orbital.

To be a little greedy, would it be likely something like the DensityPlot3D be implemented in future igor versions?
https://reference.wolfram.com/language/ref/DensityPlot3D.html
(in Applications->Potentials and Wave Functions/, then a one command takes all the jobs)

Cheers,
zh

Function fillProb()   
      variable length=41, height=61
       variable length_center=floor(length/2), height_center=floor(height/2)
    Make/O/N=((length), (length), (height)) ddd 
    Variable n=3, L=2,M=0           // set this from your panel
    Variable i,j,k,xx,yy,zz,rr,r_xy,phi,theta
    for(k=0;k<(height);k+=1)
        zz=(k-(height_center))//center  // normalized +-1 range
        for(j=0;j<(length);j+=1)
            yy=(j-length_center)//center
            for(i=0;i<(length);i+=1)
                xx=(i-length_center)//center
                // compute theta phi and rr:
                r_xy=sqrt(xx^2+yy^2)
                rr=sqrt(r_xy^2+zz^2)
                phi=atan2(xx, yy)
                theta=atan2(r_xy, zz)
                ddd[i][j][k]=WF(n, L, m, rr, theta, phi)
              endfor
        endfor
    endfor
End
 
constant a0=1  //5.29e-11 //meters,  Bohr radius
// Hydrogen atomic wavefunction
function WF(n, L, M, rr, theta, phee) // n, L, M are the principal, angular and magnetic quantum numbers, respectively
variable n, L, M, rr, theta, phee           // transform from Cartesian to spherical coordinate
variable rou=2*rr/(n*a0)     // rou is basically r, see (middle of the page) en.wikipedia.org/wiki/Hydrogen_atom
variable c= factorial(n-L-1)/(2*n)/factorial(n+L)   // a bunch of constant values
variable radial= sqrt((2/n*a0)^3*c)*exp(-rou/2)*rou^L*LaguerreA((n-L-1), (2*L+1), rou)
variable angular_squared= magsqr(SphericalHarmonics(L, M, theta, phee))
return radial*radial*angular_squared 
End
orbital_d_z2.png (7.66 KB)

May 5, 2016 at 03:12 am - Permalink

BRDF wrote:
Thanks A. G. Now I can generate the desired patterns. The attached picture is for n=3, L=2 and m=0, or the d_z^2 orbital.

I'm glad to see that you are now generating the desired patterns.
[quote]To be a little greedy, would it be likely something like the DensityPlot3D be implemented in future igor versions?
https://reference.wolfram.com/language/ref/DensityPlot3D.html[/quote]
I'm not exactly sure what I am looking at (besides an off-putting non orthographic projection). IMO, the problems with graphs like that is that they do not really give you sufficient information as dense data closer to the viewer obscures any structure behind it. I agree that an isosurface is not always ideal so it really depends if your display is stationary (i.e., single frame) or dynamic. My own preference is closer to the attached example which I generated from your code above (your resolution appears to be low).

If you would like the generating experiment and you have IP7 contact me at support@wavemetrics.com.

A.G.
WaveMetrics, Inc.
HProb.png (173.37 KB)

May 5, 2016 at 10:45 am - Permalink

Quote:
I'm not exactly sure what I am looking at (besides an off-putting non orthographic projection). IMO, the problems with graphs like that is that they do not really give you sufficient information as dense data closer to the viewer obscures any structure behind it. I agree that an isosurface is not always ideal so it really depends if your display is stationary (i.e., single frame) or dynamic. My own preference is closer to the attached example which I generated from your code above (your resolution appears to be low).

If you would like the generating experiment and you have IP7 contact me at support@wavemetrics.com.

A.G.
WaveMetrics, Inc.


I see.
I do like your experiment but...I'm still with IP6.37 and OSX 10.6.8...I will contact you when I am ready.
Thank you!

zh

May 5, 2016 at 05:45 pm - Permalink

Quote:
I do like your experiment but...I'm still with IP6.37 and OSX 10.6.8...I will contact you when I am ready.
Here is what you can do in IP6.37 (using Win 7) to improve your graphic. The attached Gizmo figure used a triangle-wave surface (converted from an isoSurface) to enhance the appearance and use various lighting effects. See the Gizmo Help and Gizmo Reference files (as well as various Visualization Example pxps) for details.
Orbital.PNG (16.99 KB)

May 10, 2016 at 09:46 am - Permalink

Thanks...I haven't come back for a while, just saw your reply.

s.r.chinn wrote:
Quote:
I do like your experiment but...I'm still with IP6.37 and OSX 10.6.8...I will contact you when I am ready.
Here is what you can do in IP6.37 (using Win 7) to improve your graphic. The attached Gizmo figure used a triangle-wave surface (converted from an isoSurface) to enhance the appearance and use various lighting effects. See the Gizmo Help and Gizmo Reference files (as well as various Visualization Example pxps) for details.

July 10, 2016 at 04:24 am - Permalink