HW4. Вычисление атомных орбиталей водорода

import numpy
import scipy.special
import scipy.misc

import npy2cube

import os
os.getcwd()

'/Users/Nataliyadi/Downloads/Golovin/4hw'

Волновая функция

$$ \psi_{nlm}(r,\theta, \varphi) = \sqrt{\left(\frac{2}{n a_{0}}\right) ^3 \frac{(n-l-1)!}{2n(n+l)!}}e^{\frac{-r}{na_0}}(\frac{2r}{na_0})^{l} L^{2l+1}_{n-l-1}(\frac{2r}{na_0})Y^{m}_{l}(\varphi,\theta)$$

n,l,m - основные квантовые числа

• n- основное число (1,2,3..)
• l – орбитальное число (0,1,2.. n-1)
• m – магнитное число (-l..+l)

def w(n,l,m,d):
    
    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)
    
    a0 = 1.0
    
    # Радиальная часть волновой функции:
    R = lambda r,n,l: (2.0*r/n/a0)**l * numpy.exp(-r/n/a0) * scipy.special.genlaguerre(n-l-1,2*l+1)(2.0*r/n/a0)
    # Умножение радиальной части на угловую:
    WF = lambda r,theta,phi,n,l,m: 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
    
    return WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)

Рассчитаем значения для первых трех уровней:

d = 30
step = float(2.*d/27)

# Зададим цикл по перебору квантовых чисел

for n in range(1,4):
    for l in range(0,n):
        for m in range(0,l+1):
            grid = w(n,l,m,d)
            name='%s-%s-%s' % (n,l,m)
            npy2cube.npy2cube(grid,(-d, -d, -d),(step, step, step),name+'.cube')

name

'3-2-2'

Скрипт для автоматической отрисовки в PyMol всех файлов .cube:

f = open ('volume_color.pml', 'w')

f.write('cmd.volume_ramp_new(\'ramp007\', [\
      -0.035, 0.00, 0.00, 1.00, 0.00, \
      -0.03,  0.00, 1.00, 1.00, 0.20, \
      -0.025, 0.00, 1.00, 0.00, 0.00, \
      0.025, 0.00, 1.00, 0.00, 0.00, \
      0.03,  0.00, 1.00, 1.00, 0.20, \
      0.035, 0.00, 0.00, 1.00, 0.00, \
    ])\n')

for n in range(1,4):
    for l in range(0,n):
        for m in range(0,l+1):
            name='%s-%s-%s' % (n,l,m)
            full_name = name+'.cube'
            vol_name = name+'_v'
            f.write('load /Users/Nataliyadi/Downloads/Golovin/4hw/first/%s, %s\n' % (full_name, name)) 
            f.write('volume %s, %s\n' % (vol_name, name))
            f.write('volume_color %s, ramp007\n' % (vol_name)) 
f.close()

Хотелось бы еще автоматическое сохранение, но почему-то через скрипт сохраняет черный экран.

from IPython.display import display,Image

image_set=[]
for n in range(1,4):
    for l in range(n):
        for m in range(l+1):
            name='%s%s%s' % (n,l,m)
            full_name = '/Users/Nataliyadi/Downloads/Golovin/4hw/png/'+name+'.png'
            image_set.append(full_name)

for i in image_set:
    display(Image(i))

%%writefile h.in
! UHF SVP 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);
end
* xyz 0 4
 H 0 0 0
*

Writing h.in

Далее работа частично шла на кластере. Запустить orca с локального компьютера, к сожалению, не получилось.

• код на bash:

#ssh -p 22025 shad@vsb.fbb.msu.ru
#cd ./hse
#mkdir didkovskaya
#cp -r /home/shad/progs/bin/orca ./didkovskaya
#orca h.inp 
#Далее с локального компьютера запущено копирование:
#scp -r -P 22025 shad@vsb.fbb.msu.ru:/home/shad/hse/didkovskaya/*.cube /Users/Nataliyadi/Downloads/Golovin/4hw

Скрипт для автоматической отрисовки файлов .cube от orca

f2 = open ('vol_color_orca.pml', 'w')

f2.write('cmd.volume_ramp_new(\'ramp008\', [\
      -0.035, 0.00, 0.00, 1.00, 0.00, \
      -0.03,  0.00, 1.00, 1.00, 0.20, \
      -0.025, 0.00, 1.00, 0.00, 0.00, \
      0.025, 0.00, 1.00, 0.00, 0.00, \
      0.03,  0.00, 1.00, 1.00, 0.20, \
      0.035, 0.00, 0.00, 1.00, 0.00, \
    ])\n')


for i in range(0,5):
    
            name="H-"+'%s' % (i)
            full_name = name+'.cube'
            vol_name = name+'_v'
            f2.write('load /Users/Nataliyadi/Downloads/Golovin/4hw/%s, %s\n' % (full_name, name)) 
            f2.write('volume %s, %s\n' % (vol_name, name))
            f2.write('volume_color %s, ramp008\n' % (vol_name)) 
f2.close()

image_set2=[]
for i in range(0,5):
            name="H-"+'%s' % (i)
            full_name = '/Users/Nataliyadi/Downloads/Golovin/4hw/cube/'+name+'.png'
            image_set2.append(full_name)

image_set2

['/Users/Nataliyadi/Downloads/Golovin/4hw/cube/H-0.png',
 '/Users/Nataliyadi/Downloads/Golovin/4hw/cube/H-1.png',
 '/Users/Nataliyadi/Downloads/Golovin/4hw/cube/H-2.png',
 '/Users/Nataliyadi/Downloads/Golovin/4hw/cube/H-3.png',
 '/Users/Nataliyadi/Downloads/Golovin/4hw/cube/H-4.png']

for i in image_set2:
    display(Image(i))

В результате были получены очень схожие изображения. Немного смущает, что три последних на глаз совершенно одинаковые.. Что касается скорости точных подсчетов и с помощью программы orca, для меня разница не была значительна и ощутима.