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.
Thank 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.
Very 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=1-Ps/Ps(end); % Normalised decreasing integral
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...
very good work
why use isosurface rather than surf?
Hey, 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.
This 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
Create scripts with code, output, and formatted text in a single executable document.