Функция npy2cube была загружена отсюда.
import numpy as np
import math
import scipy.special
import scipy.misc
from npy2cube import npy2cube
Зададим волновую функцию
def w(n,l,m,d):
# Построение сетки
x,y,z = np.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
# Переход в сферические координаты
r = lambda x,y,z: np.sqrt(x**2+y**2+z**2)
theta = lambda x,y,z: np.arccos(z/r(x,y,z))
phi = lambda x,y,z: np.arctan(y/x)
# Боровский радиус
a0 = 1.
# Радиальная функция, n - главное квантовое число, l - орбитальное квантовое число
R = lambda r,n,l: (2*r/n/a0)**l * np.exp(-r/n/a0) * scipy.special.genlaguerre(n-l-1,2*l+1)(2*r/n/a0)
# Нормировочный коэффициент
base = lambda n,l: np.sqrt((2/n/a0)**3 * math.factorial(n-l-1)/(2*n*math.factorial(n+l)))
# Волновая функция, m - магнитное квантовое число
WF = lambda r,theta,phi,n,l,m: base(n, l) * R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
# Квадрат волновой функции (действительное число)
absWF = lambda r,theta,phi,n,l,m: np.absolute(WF(r,theta,phi,n,l,m))**2
# Из опыта коллег стало известно, что орибтали получаются более похожими на истинные,
# если возвращать сумму истинной и мнимой частей.
res = WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
return np.real(res) + np.imag(res)
# Зададим цикл по перебору квантовых чисел
for n in range(1, 4):
d = 10 * n
step = d / 10
for l in range(0,n):
for m in range(-l,l+1,1):
grid= w(n, l, m, d)
name='%s_%s_%s' % (n,l,m)
# для сохранения нужно задать координаты старта grid и шаг по каждому направлению
npy2cube(grid, (-d,-d,-d), (step,step,step), f'{name}.cube')
Теперь визуализиурем полученные файлы в Pymol. Очень хотелось сделать всё автоматически и сразу, но почему-то не вышло: команды по отдельности в pymol работают, а скриптом нет. Pymol зависает и не отвечает на запросы. Вот такой вот скрипт использовался.
import __main__
__main__.pymol_argv = [ 'pymol', '-cp' ]
import pymol
pymol.finish_launching()
from pymol import cmd,stored
cmd.reinitialize()
cmd.volume_ramp_new('colors', [-0.015, 0.00, 0.00, 0.00, 0.00,
-0.01, 1.00, 0.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, 0.00, 0.20,
0.015, 0.00, 0.00, 0.00, 0.00])
names = []
for n in range(1, 4):
for l in range(0, n):
for m in range(-l, l+1):
def_name = f'{n}_{l}_{m}'
names.append(def_name)
cmd.load(f'{def_name}.cube')
cmd.volume(f'{def_name}-volume', def_name, ramp='colors')
cmd.hide()
cmd.show('volume', f'{def_name}-volume')
cmd.reset()
cmd.turn('x', 45)
cmd.turn('y', 45)
cmd.turn('z', 45)
cmd.do(f'png {name}, width=600, height=400')
Поскольку автоматически сделать картинки не удалось, а самому их все вручную сохранять не хочется, посмотрим по одному представителю s, p и d орбитали.
Image('s2.png')
Image('p3.png')
Image('d3.png')
Входной файл для программы Orca получили вот так:
with open('orca.inp', 'w') as f:
f.write('''! UHF QZVPP XYZFile
%plots Format Cube
''')
for i in range(0,14):
f.write(f'MO("H-{i}.cube",{i},0);\n')
f.write('''end
* xyz 0 4
H 0 0 0
* ''')
Затем запустили вот эти команды в терминале и пролучили модели орбиталей
export PATH=${PATH}:/srv/databases/orca
orca h.inp
Скриптом ниже тоже не получилось автоматически визуализировать по той же причине.
import __main__
__main__.pymol_argv = [ 'pymol', '-cp' ]
import pymol
pymol.finish_launching()
from pymol import cmd,stored
cmd.reinitialize()
cmd.volume_ramp_new('colors', [-0.015, 0.00, 0.00, 0.00, 0.00,
-0.01, 1.00, 0.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, 0.00, 0.20,
0.015, 0.00, 0.00, 0.00, 0.00])
for n in range(14):
cmd.load(f'H-{n}.cube')
time.sleep(0.4)
cmd.volume(f'H-{n}-volume', f'H-{n}', ramp='colors')
time.sleep(0.4)
cmd.hide()
cmd.show('volume', f'H-{n}-volume')
cmd.reset()
cmd.turn('x', 45)
cmd.turn('y', 45)
cmd.turn('z', 45)
time.sleep(0.4)
cmd.do(f'png orca_{n}, width=500, height=500')
Попросил товарища Антона Иззи визуализировать орбитали, полученные оркой, на его компьютере. У него запустилось. Так что можно посмотреть на все орбитали.
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
orcas = [f'orcas/biba_H-{i}.png' for i in range(0,14)]
pics = []
for img_path in img_paths:
pics.append(mpimg.imread(img_path))
orbs = ['s1','s2','p2','p2','p2','s3','p3','p3','p3','d3','d3','d3','d3','d3']
f = plt.figure(figsize=(17,10))
columns = 3
for i, image in enumerate(images):
f.add_subplot(len(images) / columns + 1, columns, i + 1)
plt.imshow(image)
plt.xticks([])
plt.yticks([])
plt.title(f'{orbs[i]}')
plt.show()
К сожалению, из-за каких-то неполадок с pymol не удалось получить все орбитали первым методом. Те орбитали, которые сравнил в pymol вручную, оказались очень схожими с теми, которые получены программой Orca.