4. Hydrogen orbitals calculations.¶
4.1 Importing modules.¶
In [1]:
import numpy
import scipy.special
import scipy.misc
import npy2cube
import os
import IPython
import math
import matplotlib.pyplot as plt
4.2 Defining wave function.¶
Below are all the failed attempts to return the correct function values.¶
Returning the probability density that the particle is at particular point does NOT work.¶
After several painful, bloody hours it turned out that displaying real and imaginary parts of WF should work.¶
BUT¶
Returning the square root of the sum of WF real and imaginary parts squared does NOT work.¶
Returning the sum of real and imaginary parts of WF (parts are defined as lambda functions) does NOT work.¶
Finally¶
Returning the sum of real and imaginary parts of WF (parts are defined as integers) WORKS.¶
Have no clue why does it work.¶
In [2]:
def w(n, l, m, d):
# returns a dense multi-dimensional "meshgrid" = cube
# note: 30j -- complex number, where real part (30) is the
# number of points to create between the start (-d) and stop (d) values
# so, the size of the cube is 30*30*30
x, y, z = numpy.mgrid[-d:d:30j, -d:d:30j, -d:d:30j]
# The solutions of Schrodinger equation can be expressed in spherical coordinates (r, theta, phi)
# r -- distance from the center of the nucleus
# theta -- represents the angle to the positive z-axis
# phi -- represents the angle to the positive x-axis in the xy-plan
# arctan2 quantifies arctan in all 4 quadrants
r = numpy.sqrt(x**2 + y**2 + z**2)
theta = numpy.arccos(z / r)
phi = numpy.arctan2(y, x)
# wave_function(r, theta, phi) = R(r)*Y(theta, phi), where R() -- radial variable, Y() -- angular variable
# R(r) -- depends on n and l, Y(theta, phi) -- depends on l, m
# a0 -- Bohr radius = the most probable distance between the nucleus
# and the electron in a hydrogen atom
a0 = 1.
norm_coef = math.sqrt((2/n/a0)**3 * math.factorial(n-l-1)/(2*n*math.factorial(n+l)))
# R(r) = N*p(r)*e^(-kr), where N -- positive normalizing constant, p(r) -- polynomial in r
# numpy.exp() is a e^(-kr), scipy..genlaguerre is a polynomial in r
R = (2*r/n/a0)**l * numpy.exp(-r/n/a0) * scipy.special.genlaguerre(n-l-1, 2*l+1)(2*r/n/a0)
# WF -- wave_function(r, theta, phi) = norm_coef * R(r)*Y(theta, phi), where scipy..sph_harm is a Y(theta, phi)
WF = norm_coef * R * scipy.special.sph_harm(m, l, phi, theta)
# absWF = WF**2 quantifies the probability of the electron being at a particular point
absWF = numpy.absolute(WF)**2
return numpy.real(WF) + numpy.imag(WF)
4.3 Creating .cube files for our defined WF.¶
In [3]:
where_am_i = os.getcwd()
for n in range(1,4):
d = 6.5*n
step = 2*d/30
for l in range(0,n):
for m in range(-l,l+1):
grid = w(n,l,m,d)
name = '{}_{}_{}'.format(n, l, m)
file_path = os.path.join(where_am_i, 'wf_cube', name + '.cube')
npy2cube.npy2cube(grid, (-d,-d,-d), (step,step,step), file_path)
4.4 Visualising our .cube files in pymol.¶
In [2]:
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
import pymol
pymol.finish_launching()
from pymol import cmd
In [37]:
twist_90 = ['3_1_0', '3_1_1', '3_1_-1', '3_2_2', '3_2_-2']
cmd.reinitialize()
cmd.bg_color('white')
cmd.volume_ramp_new('my_ramp', [-0.015, 1.00, 1.00, 0.00, 0.50,
-0.01, 0.20, 0.80, 0.00, 0.20,
-0.005, 0.60, 0.60, 0.00, 0.00,
0.005, 0.80, 0.40, 0.00, 0.00,
0.01, 0.80, 0.20, 1.00, 0.20,
0.015, 1.00, 0.00, 1.00, 0.50])
for n in range(1,4):
for l in range(0,n):
for m in range(-l,l+1):
name='{}_{}_{}'.format(n, l, m)
fname = 'wf_cube/' + name + '.cube'
vname = name + '_vol'
cmd.load(fname, name)
cmd.volume(vname, name, ramp='my_ramp')
if twist_90[m] in vname:
cmd.rotate('y', 90)
cmd.hide('all')
cmd.show('volume', vname)
cmd.draw(720, 720, antialias=2)
cmd.png('wf_cube/pics/' + name)
4.5 Displaying created pics.¶
In [6]:
f, ax = plt.subplots(5, 3, figsize=(10, 10))
row, col = 0, 0
for n in range(1, 4):
for l in range(0, n):
for m in range(-l, l+1):
name = '{}_{}_{}'.format(n, l, m)
pic_path = os.path.join('wf_cube', 'pics', name)
f = plt.imread('{}.png'.format(pic_path))
ax[row, col].imshow(f)
ax[row, col].axis('off')
ax[row, col].set_title('{}'.format(name))
col += 1
if col == 3: row += 1; col = 0
ax[4,2].set_visible(False)
plt.tight_layout()
4.6 Creating .cube files using orca.¶
In [5]:
# MO(name_of_the_file.cube, molecular orbital, operator)
# operator:
# electrons can have one of two spins alpha (ms = +1/2) or beta (ms = -1/2)
# in our record alpha spin -- last 0, beta spin -- last 1
# Total spin -- S = (# of alpha electrons)(+1/2) + (# of beta electrons)(-1/2)
# Molecules that contain only paired electrons would have a total spin of 0
# (all the alpha electrons cancel all of the beta electrons), and a multiplicity of one
# multiplicity is a property that describes an atom’s or molecule’s electronic structure
# Multiplicity = |2S| + 1
# a multiplicity of 1 is referred to as a singlet, two is a doublet, etc.
with open('orca_file.inp', 'w') as orca_file:
orca_file.write('''! UHF SVP XYZFile
%plots Format Cube
MO("H_1.cube",1,0);
MO("H_2.cube",2,0);
MO("H_3.cube",3,0);
MO("H_4.cube",4,0);
MO("H_5.cube",1,1);
MO("H_6.cube",2,1);
MO("H_7.cube",3,1);
MO("H_8.cube",4,1);
end
* xyz 0 2
H 0 0 0
*''')
In [3]:
import subprocess
subprocess.run('echo export PATH=$HOME/orca:$PATH >> ~/.bash_profile', shell=True)
subprocess.run('. ~/.bash_profile', shell=True)
subprocess.run('orca orca_file.inp > orca_file.out', shell=True)
4.7 Visualising orca .cube files in pymol.¶
In [7]:
cmd.reinitialize()
cmd.bg_color('white')
cmd.volume_ramp_new('my_ramp', [-0.015, 1.00, 1.00, 0.00, 0.50,
-0.01, 0.20, 0.80, 0.00, 0.20,
-0.005, 0.60, 0.60, 0.00, 0.00,
0.005, 0.80, 0.40, 0.00, 0.00,
0.01, 0.80, 0.20, 1.00, 0.20,
0.015, 1.00, 0.00, 1.00, 0.50])
for n in range(1,9):
name='H_{}'.format(n)
fname = 'orca_cube/' + name + '.cube'
vname = name + '_vol'
cmd.load(fname, name)
cmd.volume(vname, name, ramp='my_ramp')
cmd.hide('all')
cmd.show('volume', vname)
cmd.draw(720, 720, antialias=2)
cmd.png('orca_cube/pics/' + name)
4.8 Displaying created pics.¶
In [8]:
f, ax = plt.subplots(2, 4, figsize=(10, 3))
row, col = 0, 0
for i in range(1,9):
name = 'H_{}'.format(i)
pic_path = os.path.join('orca_cube', 'pics', name)
f = plt.imread('{}.png'.format(pic_path))
ax[row, col].imshow(f)
ax[row, col].axis('off')
ax[row, col].set_title('{}'.format(name))
col += 1
if col == 4: row += 1; col = 0
plt.tight_layout()
