import numpy
import time
import scipy.special
# from dungeon import slaves
import math
import scipy.misc
import matplotlib.pyplot as plt
# Выданную функцию npy2cube пришлось немного изменить, поскольку не получалось брать элементы map() по индексу
# выдавалась ошибка, нужно заменить их на списки
def npy2cube(grid, start, step, cube_f):
bohr_to_angs = 0.529177210859 # const
start = [x / bohr_to_angs for x in start]
step = [x / bohr_to_angs for x in step]
with open(cube_f, 'w') as cube_file:
cube_file.write('CPMD CUBE FILE.\n'
'OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z\n')
cube_file.write(" 1 %f %f %f\n" % (start[0], start[1], start[2]))
cube_file.write(" %i %f 0.000000 0.000000\n" %
(grid.shape[0], step[0]))
cube_file.write(" %i 0.000000 %f 0.000000\n" %
(grid.shape[1], step[1]))
cube_file.write(" %i 0.000000 0.000000 %f\n" %
(grid.shape[2], step[2]))
cube_file.write(" 1 0.000000 %f %f %f\n" %
(start[0], start[1], start[2]))
i = 0
for grid_item in numpy.nditer(grid, order='C'):
if i < 5:
cube_file.write('%f ' % (float(grid_item)))
i += 1
elif i == 5:
cube_file.write('%f\n' % (float(grid_item)))
i = 0
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.arccos(z/r(x,y,z))
phi = lambda x,y,z: numpy.arctan(y/x) + math.pi
# задавания расстояния от ядра до ближайшей орбиты электрона (Боровского радиуса). В данном случае он равен 1
a0 = 1
# Радиальная функция R, которая зависит от r, n, l
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)
# Не хватает нормировочного коэффициента. Его нужно ввести
normal = lambda n,l: numpy.sqrt((2/n/a0)**3 * math.factorial(n-l-1)/(2*n*math.factorial(n+l)))
# К итоговой волновой функции нкжно добавить ещё и нормировочный коэффициент. Здесь sph_harm - подсчёт гармоники
WF = lambda r,theta,phi,n,l,m: R(r,n,l)*normal(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
# Итоговую плотность вероятности считаем как сумму действительной и мнимой частей WF
#(если её выводить как просто квадрат модуля, то выскакивают ошибки)
itog = result = WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
vivod = numpy.real(itog) + numpy.imag(itog)
return vivod
# Используем выданный скрипт
# Нужно задать цикл по перебору квантовых чисел,
# т.е. для каждого числа n от 0 до 3 нужно перебрать орбитальные числа l от 0 до n-1,
# и для них уже перебрать магнитные числа m как 0, 1 и l
# Скрипт применить для каждой тройки таких чисел
# Также определяем нужный шаг grid
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):
grid = w(n, l, m, d)
imago = '%s_%s_%s' % (n,l,m)
npy2cube(grid, (-d,-d,-d), (step,step,step), 'dungeon_master'+imago+'.cube')
# Теперь откроем Pymol
conda install -c schrodinger pymol
Collecting package metadata (current_repodata.json): ...working... done Note: you may need to restart the kernel to use updated packages. Solving environment: ...working... done ## Package Plan ## environment location: C:\Users\msk\anaconda3 added / updated specs: - pymol The following packages will be downloaded: package | build ---------------------------|----------------- _tkinter_pymol-1 | py38hfa6e2cd_0 29 KB schrodinger apbs-1.5 | 1 266 KB schrodinger conda-4.10.0 | py38haa95532_0 2.9 MB freemol-1.158 | py_2 6 KB schrodinger glew-2.0.0 | 0 682 KB schrodinger hdf4-4.2.13 | h712560f_2 1.3 MB libholoplaycore-0.1.0_rc4 | 0 298 KB schrodinger libnetcdf-4.7.3 | h1302dcc_0 516 KB mengine-1 | 0 422 KB schrodinger mpeg_encode-1 | 0 78 KB schrodinger mtz2ccp4_px-1.0 | 1 811 KB schrodinger openvr-1.0.17 | 0 387 KB schrodinger pdb2pqr-2.1.2+pymol | py_0 236 KB schrodinger pmw-2.0.1+3 | py_3 60 KB schrodinger pymol-2.4.1 | py38h43f30b1_0 8.3 MB schrodinger rigimol-1.3 | 2 561 KB schrodinger vs2013_runtime-12.0.21005 | 1 669 KB
DEBUG menuinst_win32:__init__(198): Menu: name: 'PyMOL (Anaconda${PY_VER} ${PLATFORM})', prefix: 'C:\Users\msk\anaconda3', env_name: 'None', mode: 'user', used_mode: 'user' DEBUG menuinst_win32:create(323): Shortcut cmd is C:\Users\msk\anaconda3\PyMOLWin.exe, args are [] DEBUG menuinst_win32:create(323): Shortcut cmd is C:\Users\msk\anaconda3\PyMOLWin.exe, args are ['+2'] DEBUG menuinst_win32:create(323): Shortcut cmd is C:\Users\msk\anaconda3\PyMOLWin.exe, args are ['-S'] DEBUG menuinst_win32:create(323): Shortcut cmd is C:\Users\msk\anaconda3\PyMOLWin.exe, args are ['-St', '6']
------------------------------------------------------------ Total: 17.4 MB The following NEW packages will be INSTALLED: _tkinter_pymol schrodinger/win-64::_tkinter_pymol-1-py38hfa6e2cd_0 apbs schrodinger/win-64::apbs-1.5-1 freemol schrodinger/noarch::freemol-1.158-py_2 glew schrodinger/win-64::glew-2.0.0-0 hdf4 pkgs/main/win-64::hdf4-4.2.13-h712560f_2 libholoplaycore schrodinger/win-64::libholoplaycore-0.1.0_rc4-0 libnetcdf pkgs/main/win-64::libnetcdf-4.7.3-h1302dcc_0 mengine schrodinger/win-64::mengine-1-0 mpeg_encode schrodinger/win-64::mpeg_encode-1-0 mtz2ccp4_px schrodinger/win-64::mtz2ccp4_px-1.0-1 openvr schrodinger/win-64::openvr-1.0.17-0 pdb2pqr schrodinger/noarch::pdb2pqr-2.1.2+pymol-py_0 pmw schrodinger/noarch::pmw-2.0.1+3-py_3 pymol schrodinger/win-64::pymol-2.4.1-py38h43f30b1_0 rigimol schrodinger/win-64::rigimol-1.3-2 vs2013_runtime pkgs/main/win-64::vs2013_runtime-12.0.21005-1 The following packages will be UPDATED: conda 4.9.2-py38haa95532_0 --> 4.10.0-py38haa95532_0 Downloading and Extracting Packages _tkinter_pymol-1 | 29 KB | | 0% _tkinter_pymol-1 | 29 KB | ########## | 100% _tkinter_pymol-1 | 29 KB | ########## | 100% pdb2pqr-2.1.2+pymol | 236 KB | | 0% pdb2pqr-2.1.2+pymol | 236 KB | ########## | 100% pdb2pqr-2.1.2+pymol | 236 KB | ########## | 100% glew-2.0.0 | 682 KB | | 0% glew-2.0.0 | 682 KB | ########## | 100% glew-2.0.0 | 682 KB | ########## | 100% libnetcdf-4.7.3 | 516 KB | | 0% libnetcdf-4.7.3 | 516 KB | ####3 | 43% libnetcdf-4.7.3 | 516 KB | ########## | 100% libnetcdf-4.7.3 | 516 KB | ########## | 100% openvr-1.0.17 | 387 KB | | 0% openvr-1.0.17 | 387 KB | ########## | 100% openvr-1.0.17 | 387 KB | ########## | 100% libholoplaycore-0.1. | 298 KB | | 0% libholoplaycore-0.1. | 298 KB | ########## | 100% libholoplaycore-0.1. | 298 KB | ########## | 100% conda-4.10.0 | 2.9 MB | | 0% conda-4.10.0 | 2.9 MB | #6 | 16% conda-4.10.0 | 2.9 MB | ###7 | 37% conda-4.10.0 | 2.9 MB | #####8 | 58% conda-4.10.0 | 2.9 MB | #######8 | 79% conda-4.10.0 | 2.9 MB | ########## | 100% conda-4.10.0 | 2.9 MB | ########## | 100% mengine-1 | 422 KB | | 0% mengine-1 | 422 KB | ########## | 100% mengine-1 | 422 KB | ########## | 100% mtz2ccp4_px-1.0 | 811 KB | | 0% mtz2ccp4_px-1.0 | 811 KB | ########## | 100% mtz2ccp4_px-1.0 | 811 KB | ########## | 100% rigimol-1.3 | 561 KB | | 0% rigimol-1.3 | 561 KB | ########## | 100% rigimol-1.3 | 561 KB | ########## | 100% pmw-2.0.1+3 | 60 KB | | 0% pmw-2.0.1+3 | 60 KB | ########## | 100% pmw-2.0.1+3 | 60 KB | ########## | 100% apbs-1.5 | 266 KB | | 0% apbs-1.5 | 266 KB | ########## | 100% apbs-1.5 | 266 KB | ########## | 100% vs2013_runtime-12.0. | 669 KB | | 0% vs2013_runtime-12.0. | 669 KB | 2 | 2% vs2013_runtime-12.0. | 669 KB | ##6 | 26% vs2013_runtime-12.0. | 669 KB | ########3 | 84% vs2013_runtime-12.0. | 669 KB | ########## | 100% mpeg_encode-1 | 78 KB | | 0% mpeg_encode-1 | 78 KB | ########## | 100% mpeg_encode-1 | 78 KB | ########## | 100% pymol-2.4.1 | 8.3 MB | | 0% pymol-2.4.1 | 8.3 MB | ########## | 100% pymol-2.4.1 | 8.3 MB | ########## | 100% freemol-1.158 | 6 KB | | 0% freemol-1.158 | 6 KB | ########## | 100% freemol-1.158 | 6 KB | ########## | 100% hdf4-4.2.13 | 1.3 MB | | 0% hdf4-4.2.13 | 1.3 MB | ###5 | 36% hdf4-4.2.13 | 1.3 MB | ########2 | 82% hdf4-4.2.13 | 1.3 MB | ########## | 100% Preparing transaction: ...working... done Verifying transaction: ...working... done Executing transaction: ...working... done
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
import pymol
pymol.finish_launching()
from pymol import cmd,stored
cmd.reinitialize()
# Нужно задать цвет орбиталей
cmd.volume_ramp_new('grib', [-0.015, 0.00, 0.00, 0.00, 0.00,
-0.01, 1.00, 0.00, 1.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, 0.00, 0.00])
# Для каждого главного числа n от 1 до 3 нужно будет перебрать орбитальные числа l от 0 до n-1, и для них перебрать
# магнитные m от -l до l, и уже для каждой тройки таких чисел загрузить соответствующий файл cube и визуализировать его
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'dungeon_master/{name}.cube')
cmd.volume(f'{name}-v', name, ramp='grib')
cmd.hide()
cmd.show('volume', f'{name}-v')
cmd.reset()
cmd.turn('x', 45)
cmd.turn('y', 45)
cmd.turn('z', 45)
time.sleep(0.3)
cmd.do(f'png {name}, width={480}, height={370}')
# Посмотрим на изображения орбиталей
import glob
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
numbers = ['1_0_0', '2_0_0', '2_1_0', '2_1_1', '2_1_-1', '3_0_0', '3_1_0', '3_1_1', '3_1_-1', '3_2_0', '3_2_1', '3_2_-1', '3_2_2', '3_2_-2']
img_paths = [f'{i}.png' for i in numbers]
zaborchik = []
for img_path in img_paths:
zaborchik.append(mpimg.imread(img_path))
orbitals = ['s1','s2','p2','p2','p2','s3','p3','p3','p3','d3','d3','d3','d3','d3']
f = plt.figure(figsize=(30,20))
columns = 3
for i, image in enumerate(zaborchik):
f.add_subplot(len(zaborchik) / columns + 1, columns, i + 1)
plt.imshow(image)
plt.xticks([])
plt.yticks([])
plt.title(f'{orbitals[i]}')
plt.show()
# Создаём файл для Орки
# Далее этот файл будет обработан на сервере кодомо программой orca
with open('orcanavt.inp', 'w') as bullet:
bullet.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 4
H 0 0 0
* '''
)
# Теперь с этим проводим всё то же самое, что делали с прошлой выдачей
cmd.reinitialize()
cmd.volume_ramp_new('grib', [-0.015, 0.00, 0.00, 0.00, 0.00,
-0.01, 1.00, 0.00, 1.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, 0.00, 0.00])
for kl in range(0, 14):
gorgonopsis = f'H-{kl}'
cmd.load(f'dungeon_master/{gorgonopsis}.cube')
cmd.volume(f'{gorgonopsis}-v', gorgonopsis, ramp='grib')
cmd.hide()
cmd.show('volume', f'{gorgonopsis}-v')
cmd.reset()
cmd.turn('x', 45)
cmd.turn('y', 45)
cmd.turn('z', 45)
time.sleep(0.3)
cmd.do(f'png boba_{gorgonopsis}, width={480}, height={370}')
# Теперь загрузим полученные картинки
import glob
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
img_paths = [f'boba_H-{i}.png' for i in range(14)]
kartinki = []
for img_path in img_paths:
kartinki.append(mpimg.imread(img_path))
orbitals = ['s1','s2','p2','p2','p2','s3','p3','p3','p3','d3','d3','d3','d3','d3']
f = plt.figure(figsize=(30,20))
columns = 3
for i, image in enumerate(kartinki):
f.add_subplot(len(kartinki) / columns + 1, columns, i + 1)
plt.imshow(image)
plt.xticks([])
plt.yticks([])
plt.title(f'{orbitals[i]}')
plt.show()
# Полученные с помощью программы orca изображения орбиталей практически не отличаются от рассчётных (как будто стали
# более распухшими)