♂♂♂Вычисление орбиталей водорода♂♂♂

import numpy
import scipy.special
import scipy.misc
import npy2cube
import os
import IPython
import math
import matplotlib.pyplot as van
#import ♂finger♂ into ♂ass♂

После импорта необходимых библиотек (файл npy2cube.py, полученный с ♂website♂ http://kodomo.fbb.msu.ru/~golovin/ipynb/npy2cube.py, был адаптирован для Python 3) я написал комметарии к функциям расчета орбиталей.

#from ♂dungeon♂ import ♂slave♂
def weewee(aniki, leatherman, m, d):
    
    # Задаем декартову кубическую ортогональную систему координат с шагом сетки 30
    x, y, z = numpy.mgrid[-d:d:30j, -d:d:30j, -d:d:30j]


    # Переходим в сферические полярные координаты.
    # r -- радиус-вектор, расстояние от центра ядра
    # theta -- зенит
    # phi -- азимут
    r = numpy.sqrt(x**2 + y**2 + z**2)
    theta = numpy.arccos(z / r)
    phi = numpy.arctan2(y, x)

    # Волновую функцию можно разделить на произведение компонент, зависящих только от радиус-вектора и только от углов
    
    # a0 -- Боров радиус - наиболее вероятное расстояние от ядра до элеткрона
    a0 = 1.0
    # численный множитель для волновой функции
    norm_coef = math.sqrt((2/aniki/a0)**3 * math.factorial(aniki-leatherman-1)/(2*aniki*math.factorial(aniki+leatherman)))

    # добавляем радиальные функции 
    R = (2*r/aniki/a0)**l * numpy.exp(-r/aniki/a0) * scipy.special.genlaguerre(aniki-leatherman-1, 2*leatherman+1)(2*r/aniki/a0)

    # Сферические гармоники
    WF = norm_coef * R * scipy.special.sph_harm(m, leatherman, phi, theta)

    # Квадрат нормы волновой функции - плотность вероятности нахождения электрона в точке
    absWF = numpy.absolute(WF)**2

    return numpy.real(WF) + numpy.imag(WF)

# Переберем квантовые числа в цикле, задав шаг и размер ячейки

            
leatherclub = os.getcwd()

for harrington in range(1,4):
    daddy = 6.5*harrington
    fist = 2*daddy/30
    for milos in range(0,harrington):
        for darkholme in range(-milos,milos+1):
            latexglove = weewee(harrington,milos,darkholme,daddy)
            fantasy = '{}_{}_{}'.format(harrington, milos, darkholme)
            gym_path = os.path.join(leatherclub, 'wf_cube', fantasy + '.cube')
            npy2cube.npy2cube(latexglove, (-daddy,-daddy,-daddy), (fist,fist,fist), gym_path)

Запустим PyMol локально для визуализации орбиталей и who's ♂BOSS♂ of this ♂gym♂

import __main__

__main__.pymol_argv = [ 'pymol', '-x' ]

import pymol 
pymol.finish_launching()
from pymol import cmd

import time
get_your_whole_arm_up_there = ['3_1_0', '3_1_1', '3_1_-1', '3_2_2', '3_2_-2']

cmd.reinitialize()
cmd.bg_color('white')

cmd.volume_ramp_new('rip_the_skin', [-0.015, 1.00, 1.00, 0.00, 0.50, 
                                -0.01,  0.20, 0.80, 0.00, 0.20, 
                                -0.005, 0.60, 0.60, 0.00, 0.00, 
                                0.005, 0.80, 0.40, 0.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(1,4):
    for l in range(0,n):
        for master in range(-l,l+1):
            name='{}_{}_{}'.format(n, l, master)
            fetishname = 'wf_cube/' + name + '.cube'
            vname = name + '_vol'
            cmd.load(fetishname, name)
            cmd.volume(vname, name, ramp='rip_the_skin')
            if get_your_whole_arm_up_there[master] in vname:
                cmd.rotate('y', 90)
            cmd.hide('all')
            cmd.show('volume', vname)
            time.sleep(0.1) #иначе изображения кривые
            cmd.draw(720, 720, antialias=2)
            cmd.png('wf_cube/nextdoor/' + name)

f, ax = van.subplots(5, 3, figsize=(10, 10))

deep, dark = 0, 0
for n in range(1, 4):
    for l in range(0, n):
        for m in range(-l, l+1):
            name = '{}_{}_{}'.format(n, l, m)
            pic_path = os.path.join('wf_cube', 'nextdoor', name)
            f = van.imread('{}.png'.format(pic_path))
            ax[deep, dark].imshow(f)
            ax[deep, dark].axis('off')
            ax[deep, dark].set_title('{}'.format(name))
            dark += 1
            if dark == 3: deep += 1; dark = 0

ax[4,2].set_visible(False)    
van.tight_layout()

Создадим входные файлы для ORCA и потом их визуализируем - прямо как ♂Performance Artist♂

with open('orca_file.inp', 'w') as take_it_boy:
    take_it_boy.write('''! UHF SVP XYZFile
%plots Format Cube
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",1,1);
MO("H_6.cube",2,1);
MO("H_7.cube",3,1);
MO("H_8.cube",4,1);
end
* xyz 0 2
H 0 0 0
*''')

%%bash
export PATH=${PATH}:/srv/databases/orca
orca orca_file.inp

import time
cmd.reinitialize()
cmd.bg_color('white')

cmd.volume_ramp_new('bondage', [-0.015, 1.00, 1.00, 0.00, 0.50, 
                                -0.01,  0.20, 0.80, 0.00, 0.20, 
                                -0.005, 0.60, 0.60, 0.00, 0.00, 
                                0.005, 0.80, 0.40, 0.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(1,9):
    name='H_{}'.format(n)
    fetishname = name + '.cube'
    vname = name + '_vol'
    cmd.load(fetishname, name)
    cmd.volume(vname, name, ramp='bondage')
    cmd.hide('all')
    cmd.show('volume', vname)
    time.sleep(0.1)
    cmd.draw(720, 720, antialias=2)
    cmd.png('suction/' + name)

lube, spank = van.subplots(2, 4, figsize=(10, 3))

row, col = 0, 0
for i in range(1,9):
    name = 'H_{}'.format(i)
    pic_path = os.path.join('suction', name)
    lube = van.imread('{}.png'.format(pic_path))
    spank[row, col].imshow(lube)
    spank[row, col].axis('off')
    spank[row, col].set_title('{}'.format(name))
    col += 1
    if col == 4: row += 1; col = 0


van.tight_layout()

На одноэлектронной модели существенной разницы между приближенным вычислением и специализированной программой не видно. Заранее ♂Thank♂you♂Sir♂ за проверку