import numpy
import math
import scipy.special
import scipy.misc
import os
import py3Dmol
import psi4

import npy2cube

import matplotlib.pyplot as plt

def w(n,l,m,d):
    
    x,y,z = numpy.mgrid[-d:d:50j,-d:d:50j,-d:d:50j]
    
    r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2)
#Distance from cener
    theta = lambda x,y,z: numpy.arccos(z/r(x,y,z))
#Angle theta
    phi = lambda x,y,z: numpy.arctan(y/x)
#Angle phi
    
    a0 = 0.529177
#Bohr radius 
    
    R = lambda r,n,l: ((((2/n/a0)**3)*math.factorial(n-l-1)/2/n/math.factorial(n+l))**0.5)*(2*r/n/a0)**l*numpy.exp(-r/n/a0)*scipy.special.genlaguerre(n-l-1,2*l+1)(2*r/n/a0)
#Radial part

##########
    WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
#Wave function
    absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2
#Probability in voxel?????
#########

    prob = R(r(x,y,z),n,l) * (numpy.real(scipy.special.sph_harm(m,l,phi(x,y,z),theta(x,y,z))) + numpy.imag(scipy.special.sph_harm(m,l,phi(x,y,z),theta(x,y,z))))
#Get probability sqr in voxel
    return prob**2


def w2(n,l,m,d):
    
    x,y,z = numpy.mgrid[-d:d:50j,-d:d:50j,-d:d:50j]
    
    r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2)
#Distance from cener
    theta = lambda x,y,z: numpy.arccos(z/r(x,y,z))
#Angle theta
    phi = lambda x,y,z: numpy.arctan(y/x)
#Angle phi
    
    a0 = 0.529177
#Bohr radius 
    
    R = lambda r,n,l: ((((2/n/a0)**3)*math.factorial(n-l-1)/2/n/math.factorial(n+l))**0.5)*(2*r/n/a0)**l*numpy.exp(-r/n/a0)*scipy.special.genlaguerre(n-l-1,2*l+1)(2*r/n/a0)
#Radial part

    WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
#Wave function
    absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2
#Probability in voxel?????

    return absWF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)

try:
    os.remove(name+'.cube')
except:
    pass
# определите шаг grid при заданном диапозоне от -d до d
d=15
step=150

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

for n in range(0,4):
    for l in range(0,n):
        for m in range(0,l+1,1):
            grid=w(n,l,m,d) 
            name='%s-%s-%s' % (n,l,m)
            # для сохранения нужно задать координаты старта grid и шаг по каждому направлению
            npy2cube.npy2cube(grid,(-d,-d,-d),(step,step,step),name+'.cube')

view = py3Dmol.view()
alpha = open('1-0-0.cube','r').read()
view.addVolumetricData(alpha, "cube", {'isoval': 0.001, 'color': "red", 'opacity': 0.5})
view.zoomTo()
view.show()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

import __main__

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

import pymol
 
pymol.finish_launching()
 
from pymol import cmd

from IPython.display import Image

#color scheme
cmd.volume_ramp_new('ramp', [\
      -0.005, 0.00, 0.00, 1.00, 0.00, \
      -0.001,  0.00, 0.00, 1.00, 0.20, \
      -0.0015, 0.00, 0.00, 1.00, 0.00, \
      0.005, 0.00, 1.00, 1.00, 0.10, \
      0.001,  3.00, 0.00, 0.00, 0.20, \
      0.0015, 0.00, 0.00, 1.00, 0.00, \
    ])

#make pics
for n in range(0,4):
    for l in range(0,n):
        for m in range(0,l+1,1):
            name = '%s-%s-%s' % (n,l,m)
            cmd.do('''
            reinitialize
            load ./{0}.cube, {0}
            volume {0}_vol, {0}
            volume_color {0}_vol, ramp
            turn x, 45
            turn y, 45
            png {0}.png, width=400, height=300'''.format(name))

for n in range(1,4):
    fig = plt.figure(figsize = (2*n, 2*n))
    for l in range(0,n):
        for m in range(0, l+1): 
            ax = fig.add_subplot(n, n, 1+l*n+m)
            ax.imshow(plt.imread(f'{n}-{l}-{m}.png'))
            ax.set_title(f'{n}-{l}-{m}')
            ax.set_xticks([])
            ax.set_yticks([])
    fig.tight_layout()
plt.show()

!pip install -q condacolab
import condacolab
condacolab.install()

%%capture
! mamba install -c anaconda intel-openmp
! mamba install -c psi4 psi4
! mamba install -c conda-forge psiresp
!pip install py3Dmol

import psi4
import numpy as np
psi4.core.set_output_file('output.dat')

# as https://sharifsuliman.medium.com/visualizing-frontier-orbitals-with-psi4-and-python-c01f335afaad

m = psi4.geometry('''
         0 1
         C     0.000000     0.000000     0.000000
''')
psi4.set_options({'reference': 'uhf'})
eng, wave_function = psi4.energy('b3lyp/6-311G(2DF,P)', molecule=m, return_wfn=True)

cubeprop_tasks = []
cubeprop_tasks.append('ORBITALS')

#возьмем орбиталей с запасом
psi4.set_options(
{
  'scf_type': 'df',
  'g_convergence': 'gau_tight',
  'freeze_core': 'true',
  'CUBEPROP_TASKS': cubeprop_tasks,
  'CUBEPROP_FILEPATH': './',
  'CUBEPROP_ORBITALS': list(range(1, 17)),
}
)

psi4.cubeprop(wave_function)

import __main__

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

import pymol
 
pymol.finish_launching()
 
from pymol import cmd

from IPython.display import Image

#color scheme
cmd.volume_ramp_new('ramp', [\
      -0.005, 0.00, 0.00, 1.00, 0.00, \
      -0.001,  0.00, 0.00, 1.00, 0.20, \
      -0.0015, 0.00, 0.00, 1.00, 0.10, \
      0.005, 0.00, 1.00, 1.00, 0.10, \
      0.001,  3.00, 0.00, 0.00, 0.20, \
      0.0015, 0.00, 0.00, 1.00, 0.00, \
    ])

#make pics
for n in range(1,17):
    cmd.do('''
    reinitialize
    load ./Psi_a_{0}.cube, {0}
    volume {0}_vol, {0}
    volume_color {0}_vol, ramp
    turn x, 45
    turn y, 45
    png {0}.png, width=400, height=300'''.format(n))

fig = plt.figure(figsize = (2*3, 2*5))
for n in range(1,17):
        ax = fig.add_subplot(4, 4, n)
        ax.imshow(plt.imread(f'{n}.png'))
        ax.set_title(f'Psi-{n}')
        ax.set_xticks([])
        ax.set_yticks([])
fig.tight_layout()
plt.show()

Вышли какие-то половинки орбиталей. Может быть, это части решения уравнения Шредингера с разными знаками частей?..