Импортируем все необходимое

import npy2cube

import numpy as np
import scipy.special
import scipy.misc
import math

conda install -c schrodinger pymol

Collecting package metadata (current_repodata.json): done
Solving environment: done

# All requested packages already installed.


Note: you may need to restart the kernel to use updated packages.

import __main__
__main__.pymol_argv = [ 'pymol', '-c' ]
import pymol
pymol.finish_launching()
from pymol import cmd,stored

 Detected 8 CPU cores.  Enabled multithreaded rendering.

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
    
    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)  
    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: np.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)

# Для каждого набора квантовых чисел подсчитаем орбитали

for n in range(1,4):
    for l in range(0,n):
        for m in range(-l,l+1,1):
            d = 10 * n
            step=2*d/(30-1)
            grid=w(n, l, m, d)
            name='%s_%s_%s' % (n,l,m)
            npy2cube.npy2cube(grid, (-d,)*3, (step,)*3, 'cube/'+name+'.cube')

Далее рассчитагнные .cube файлы вручную визуализировала в PyMol

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'cube/{i}.png'for i in numbers]
images = []
for img_path in img_paths:
    images.append(mpimg.imread(img_path))

orbitals = ['s1','s2','p2','p2','p2','s3','p3','p3','p3','d3','d3','d3','d3','d3']
f = plt.figure(figsize=(20,10))
columns = 4
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'{orbitals[i]}')

plt.show()

/opt/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:15: MatplotlibDeprecationWarning: Passing non-integers as three-element position specification is deprecated since 3.3 and will be removed two minor releases later.
  from ipykernel import kernelapp as app

Файл для orca

print('''
h.inp:
! 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);
end
* xyz 0 4
 H 0 0 0
*
''')

h.inp:
! 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);
end
* xyz 0 4
 H 0 0 0
*

Полученные рассчитанные орбитали визуализируем в PyMol

import glob
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

numbers = ['H-1', 'H-2', 'H-3', 'H-4']
img_paths = [f'cube/{i}.png' for i in numbers]
images = []
for img_path in img_paths:
    images.append(mpimg.imread(img_path))

orbitals = ['s1','s2','p2','p2','p2','s3','p3','p3','p3','d3','d3','d3','d3','d3']
f = plt.figure(figsize=(20,10))
columns = 4
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'{orbitals[i]}')

plt.show()

/opt/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:15: MatplotlibDeprecationWarning: Passing non-integers as three-element position specification is deprecated since 3.3 and will be removed two minor releases later.
  from ipykernel import kernelapp as app

Надо сказать, что орбитали не похожи на реальные ни в случае Orca, ни в случае локального рассчета. Не понимаю почему такая картина - возможно проблема в визуализации PyMol, которую нужнол подкрутить.