Commit 8f56cc05 authored by Ivo Filot's avatar Ivo Filot
Browse files

Initial commit

parents
# -*- coding: utf-8 -*-
#
# Create contour plots of the atomic orbitals
#
# MIT License
#
# Copyright (c) 2020 Ivo Filot <i.a.w.filot@tue.nl>
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
import numpy as np
from math import factorial
from scipy.special import assoc_laguerre
from scipy.special import lpmv
import matplotlib.pyplot as plt
from scipy.special import sph_harm
def main():
# 1s orbital
plot_orbital(1, 0, 0, 15)
# 2pz orbital
plot_orbital(2, 1, 0, 15, plane = 'xz')
# 2px orbital
plot_orbital(2, 1, 1, 15)
#
# No need to change the contents below
# (or do so at your own risk)
#
def plot_orbital(n, l, m, sz, plane='xy'):
"""
Plot orbital
@param n : principle quantum number
@param l : angular quantum number
@param m : magnetic quantum number
@param sz : size of the domain (in Bohr units) for the contour plot
"""
x = np.linspace(-sz,sz,100)
y = np.linspace(-sz,sz,100)
xx, yy = np.meshgrid(x, y)
# build containers
r = []
theta = []
phi = []
if plane == 'xy':
for x,y in zip(xx,yy):
r.append(np.sqrt(x**2 + y**2))
theta.append(np.arctan2(y,x))
phi.append(np.pi/2)
r = np.array(r)
theta = np.array(theta)
elif plane == 'xz':
for x,y in zip(xx,yy):
r.append(np.sqrt(x**2 + y**2))
theta.append(0)
phi.append(np.arccos(y/r[-1]))
elif plane == 'yz':
for x,y in zip(xx,yy):
r.append(np.sqrt(x**2 + y**2))
theta.append(np.pi / 2)
phi.append(np.arcsin(y/r[-1]))
else:
raise "Invalid plane set: %s" % plane
# cast to arrays
r = np.array(r)
theta = np.array(theta)
phi = np.array(phi)
# note that the longitudinal angle between pole and xy-plane is pi/2
psi = wf(n, l, m, r, theta, phi)
fig, ax = plt.subplots(dpi=72, figsize=(6,5))
ax.contour(xx, yy, psi, colors='black')
c = ax.contourf(xx, yy, psi, cmap=plt.get_cmap('PiYG'))
plt.grid(linestyle='--', color='black', alpha=0.5)
plt.xlim(-sz,sz)
plt.ylim(-sz,sz)
# set axis labels
if plane == 'xy':
plt.xlabel("$x$ in Bohr")
plt.ylabel("$y$ in Bohr")
elif plane == 'xz':
plt.xlabel("$x$ in Bohr")
plt.ylabel("$z$ in Bohr")
elif plane == 'yz':
plt.xlabel("$y$ in Bohr")
plt.ylabel("$z$ in Bohr")
cmap = plt.colorbar(c)
plt.tight_layout()
plt.show()
def wf(n,l,m,r,theta,phi):
"""
Construct the wave function
"""
return radial(n,l,r) * angular(l,m,theta,phi) * r
def angular(l,m,theta,phi):
"""
Construct the angular part of the wave function
"""
# see: https://en.wikipedia.org/wiki/Spherical_harmonics#Real_form
#
# this create so-called Tesseral spherical harmonics
#
if m == 0:
return np.real(sph_harm(m,l,theta,phi))
elif m > 0:
return np.real(sph_harm(-m,l,theta,phi) + (-1)**m * sph_harm(m,l,theta,phi)) / np.sqrt(2.0)
elif m < 0:
return np.imag(sph_harm(m,l,theta,phi) - (-1)**m * sph_harm(-m,l,theta,phi)) / np.sqrt(2.0)
def radial(n,l,r):
"""
This is the formulation for the radial wave function as encountered in
Griffiths "Introduction to Quantum Mechanics 3rd edition"
"""
n = float(n)
l = float(l)
a = 1.0
rho = 2.0 * r / (n * a)
val = np.sqrt((2.0 / (n * a))**3) * \
np.sqrt(factorial(n - l - 1) / (2 * n * factorial(n + l))) * \
np.exp(-0.5 * rho) * \
rho**l * \
assoc_laguerre(rho, n-l-1, 2*l+1)
return val
def polar(l,m,theta,phi):
"""
Construct polar part of the angular wave function
"""
pre = np.sqrt((2 * l + 1) * factorial(l - np.abs(m)) / factorial(l + np.abs(m)))
return pre * lpmv(np.abs(m), l, np.cos(theta))
def azimuthal(m, phi):
"""
Construct azimuthal part of the angular wave function
"""
pre = 1.0 / np.sqrt(4.0 * np.pi)
if m == 0:
return pre;
if(m > 0):
return pre * np.cos(m * phi)
else:
return pre * np.sin(-m * phi)
if __name__ == '__main__':
main()
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment