version 1.0.0.0 (1.32 MB) by
Evan

mix, hybridize, and visualize hydrogen orbital wavefunctions

The script PlotHydrogenMolecularOrbital.m is a stand-alone sandbox for hydrogen electronic orbitals. Wavefunctions are created using Laguerre and Legendre polynomials, exponential functions, and the Bohr radius. The result is an easy way to plot orbitals with realistic sizes and density distributions.

To get started simply run any of the included scripts. Each script should contain all necessary functions. Each file uses PlotHydrogenMolecularOrbital.m to create and plot simple hybrid orbitals along with basic s, p, d, and f orbitals. These should suffice for most basic users simply looking for a visualization of orbitals.

PlotHydrogenMolecularOrbital.m by default plots a 3d_z^2 orbital.

For advanced users, any of the wavefunction outputs can be converted to DX-format files using mat2dx.m (check my file exchange), which can be used to overlay the orbitals with other atomic/orbital data. Don't be afraid to send an email should there be an issue with formats or use.

Evan (2020). plot Hydrogen Atom Molecular Orbital (https://www.mathworks.com/matlabcentral/fileexchange/44604-plot-hydrogen-atom-molecular-orbital), MATLAB Central File Exchange. Retrieved .

Created with
R2010a

Compatible with any release

Create scripts with code, output, and formatted text in a single executable document.

YCJoseph SmerdonMaybe it shouldn't be called 'molecular' orbital, since it isn't around a molecule.

Jaehong JeongThank you for your effort. But, be careful! There are quite many errors in the code.

1. cons.scale = 10^(-9.5); scaling length of calculation? It doesn't make sense.

2. cons.meter2ang: we don't need it. Better to play on nanometer or angstrom unit. You can easily achieve it just by defining Bohr radius a0 in the unit you want.

3. cons.a0 = 5.2820 * 10^(-11); the Bohr radius is ~0.52"9"2 ang or 0.052"9"2 nm.

4. in the radial function, you don't need scalFac2 if you define LaguerreP correctly.

5. the nuclear component should be (2*r/(a0*n)).^l .* exp(-r/(a0*n)). ".^l" is missing.

6. the angular function is wrong. The first part for m==2 is just a trick. normFac also doesn't make sense. It would be better to take references and use correct formula for the wave function.

Mark TVery nice - thank you for sharing this. One comment: it is often more appropriate to plot the surface(s) encompassing a specified fraction of the total *integral* of the probability density (also yields a more intuitive choice of the threshold to plot the surface).

I did this in a similar code:

[Ps,is]=sort(P(:),1,'descend');

Ps=cumsum(Ps);

Ps=1-Ps/Ps(end); % Normalised decreasing integral

Psi=zeros(size(P));

Psi(is)=Ps;

Then you can do an isosurface on Psi with say 0.05 tolerance to encompass 95% of the total probability. Maybe youcoudl add this as an option one day...

mohamed khoutoulvery good work

Zhang kewenwhy use isosurface rather than surf?

Fuad AlamiHey, thanks for the amazing work. I'm looking to use this as the base for code I am trying to implement in order to model Benzene. I was wondering where I could see the derivations of your formula inorder to adjust for the Z value.

Cheers!

David VargasThis is so cool!! I set the n,l,m numbers to inputs so I can generate new plots quickly. I really wish my Chem 1 professor had used this to show us orbitals with different quantum numbers. This blows my mind!

Try, n=5, l=3, m=3