import numpy
import scipy.special
import scipy.misc
import math
import npy2cube
import time
def w(n,l,m,d):
#На вход подаем квантовые числа n, l, m
x,y,z = numpy.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
#Переходим из декартовых в сферические координаты
r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2)
theta = lambda x,y,z: -numpy.arcsin(z/r(x,y,z)) + math.pi/2
phi = lambda x,y,z: numpy.arctan2(y, x) + math.pi
#Задаем радиус Бора
a0 = 1
#Задаем радиальную часть волновой функции
R = lambda r,n,l: (2*r/n/a0)**l * numpy.exp(-r/n/a0) * scipy.special.genlaguerre(n-l-1,2*l+1)(2*r/n/a0)
#Нормировочный коэффициент
b = lambda n,l: numpy.sqrt((math.factorial(n-l-1)/(2*n*math.factorial(n+l)))*(2/n/a0)**3)
#Волновая функция
WF = lambda r,theta,phi,n,l,m: b(n, l) * R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
#Плотность вероятности волновой функции
absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2
out = WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
return numpy.real(out) + numpy.imag(out)
#Цикл, перебирающий квантовые числа
for n in range(0,4):
d = 9 * n
step = d / 17
for l in range(0,n):
for m in range(-l, l+1):
grid = w(n,l,m,d)
name = f'{n}-{l}-{m}'
#задаем координаты старта grid и шаг по каждому направлению
npy2cube.npy2cube(grid,(-d,-d,-d),(step,step,step),name+'.cube')
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
import pymol
pymol.finish_launching()
#Делаем скрипт для визуализации в PyMOL:
from pymol import cmd,stored
cmd.reinitialize()
# Эти цифры отвечают за цвет орбиталей
cmd.volume_ramp_new('ramp007', [-0.015, 1.00, 0.00, 0.00, 0.00,
-0.01, 1.00, 1.00, 0.00, 0.20,
-0.005, 0.00, 0.00, 1.00, 0.00,
0.005, 0.00, 0.00, 1.00, 0.00,
0.01, 0.00, 1.00, 1.00, 0.20,
0.015, 0.00, 0.00, 1.00, 0.00])
for n in range(1, 4):
for l in range(0, n):
for m in range(-l, l+1):
name = f'{n}-{l}-{m}'
cmd.load(f'{name}.cube')
cmd.volume(f'{name}-vol', name, ramp='ramp007')
#Сохраняем картинки орбиталей
slab = ['1-0-0', '2-0-0', '2-1-0', '3-0-0', '3-1-0', '3-2-0']
for n in range(1, 4):
for l in range(0, n):
for m in range(-l, l+1):
name = f'{n}-{l}-{m}'
cmd.hide()
cmd.show('volume', f'{name}-vol')
cmd.reset()
cmd.turn('x', 90)
if name in slab:
cmd.clip('slab', 0)
time.sleep(0.1)
cmd.do(f'png {name}, width={480}, height={270}')
import matplotlib.pyplot as plt
orbitales = {0: 's', 1: 'p', 2: 'd'}
_, axarr = plt.subplots(7, 2, figsize=(10, 21.33)) #будет 7 строк по 2 столбца с изображениями
row, col = 0, 0
for n in range(1, 4):
for l in range(0, n):
for m in range(-l, l+1):
name = f'{n}-{l}-{m}'
f = plt.imread(f'{name}.png')
axarr[row, col].axis('off')
axarr[row, col].imshow(f)
axarr[row, col].set_title(f'{n}{orbitales[l]}',
{'fontsize': 35}, pad=10)
col += 1
if col == 2:
row += 1
col = 0
plt.tight_layout()
#Рассчитаем орбитали в программе Orca. Создадим текстовый файл h.inp:
with open('hydrogen.inp', 'w') as f:
f.write('''! UHF QZVPP XYZFile
%plots Format Cube
MO("H-0.cube",0,0);
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",5,0);
MO("H-6.cube",6,0);
MO("H-7.cube",7,0);
MO("H-8.cube",8,0);
MO("H-9.cube",9,0);
MO("H-10.cube",10,0);
MO("H-11.cube",11,0);
MO("H-12.cube",12,0);
MO("H-13.cube",13,0);
end
* xyz 0 2
H 0 0 0
* ''')
cmd.reinitialize()
cmd.volume_ramp_new('ramp007', [-0.015, 0.00, 1.00, 1.00, 0.50,
-0.01, 0.20, 0.80, 1.00, 0.20,
-0.005, 0.40, 0.60, 1.00, 0.00,
0.005, 0.60, 0.40, 1.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(14):
name = f'H-{n}'
cmd.load(f'{name}.cube')
cmd.volume(f'{name}-vol', name, ramp='ramp007')
slab = ['H-0', 'H-1', 'H-2', 'H-5', 'H-13', 'H-14']
for n in range(14):
name = f'H-{n}'
cmd.hide()
cmd.show('volume', f'{name}-vol')
cmd.reset()
if name in slab:
cmd.clip('slab', 0)
time.sleep(0.1)
cmd.do(f'png {name}')
orbitales = ['1s', '2s', '2p', '2p', '2p', '3s', '3d', '3d', '3d', '3d', '3d', '3p', '3p', '3p']
orbitales = (n for n in orbitales)
_, axarr = plt.subplots(7, 2, figsize=(10, 21.33))
row, col = 0, 0
for n in range(14):
name = f'H-{n}'
f = plt.imread(f'{name}.png')
axarr[row, col].axis('off')
axarr[row, col].imshow(f)
axarr[row, col].set_title(f'{next(orbitales)}', {'fontsize': 35}, pad=10)
col += 1
if col == 2:
row += 1
col = 0
plt.tight_layout()