Вычисление атомных орбиталей водорода: построить электронную плотность в одно-электронном атоме и сравнить с известными программами.

import glob
import numpy
import scipy.special
import scipy.misc

%%bash

wget http://kodomo.fbb.msu.ru/~golovin/ipynb/npy2cube.py

import npy2cube

Задаем волновую функцию Ψ(r,θ,φ) = R(r) * Y(θ,φ)

r,θ,φ - сферические координаты n (главное), l (орбитальное), m (магнитное) - квантовые числа

def w(n,l,m,d):
    
    # пространственные координаты системы (задаются параметром 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)
 
    # радиус Бора ближайшей к ядру орбиты электрона атома водорода
    # значение исправлено на 0.53 А
    a0 = 0.53

    # определение нормированной радиальной функции через присоединенные полиномы Лагерра
    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)
    
    # разделение волновой функции на радиальную и угловую части
    # задаем угловую часть (сферические гармоники)
    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=15
step=2*d

# цикл по перебору квантовых чисел
# для первых трех энергетических уровней
# для основного состояния атома водорода n=1, l=0, m=0

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)
            # для сохранения .cube задаем координаты старта grid и шаг по каждому направлению
            npy2cube.npy2cube(grid,(-d,-d,-d),(step, step, step),(name+'.cube'))

# launch PyMOL 
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]

import pymol
from pymol import cmd
pymol.finish_launching()

# set image
import time
from IPython.display import Image

# set custom color ramp 
cmd.set("ray_volume")
cmd.volume_ramp_new('ramp007', [\
      -0.009, 1.00, 1.00, 0.00, 0.005, \
      -0.005, 1.00, 0.00, 0.00, 0.007, \
      -0.001, 0.00, 0.00, 1.00, 0.0001, \
      0.004, 1.00, 1.00, 0.00, 0.005, \
      0.005, 0.00, 0.00, 0.00, 0.007, \
      0.006, 1.00, 1.00, 0.00, 0.005, \
    ])

# строим электронную плотность 
for file in glob.glob('*.cube'):
    cmd.load(file)
    cmd.volume("vol", file.rsplit('.')[0], ramp="ramp007")
    time.sleep(3)
    cmd.ray()
    cmd.png('%s.png' % file.rsplit('.')[0])

for image in sorted(glob.glob('*.cube.png'), key=str):
    display(Image(image))

%%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
*

%%bash 

# cp -r ../../progs/bin/orca .
orca h.in

# строим электронную плотность для .cube-файлов, полученных с orca
for file in glob.glob('H*.cube'):
    cmd.load(file)
    cmd.volume("vol", file.rsplit('.')[0], ramp="ramp007")
    time.sleep(3)
    cmd.ray()
    cmd.png('%s-orca.png' % file.rsplit('.')[0])

for image in sorted(glob.glob('*-orca.png'), key=str):
    display(Image(image))

%%bash
jupyter nbconvert --to html hw4_filippova.ipynb