Обитали в Pymol

import numpy as np
import pandas as pd
import scipy.special
import scipy.misc
import math
import matplotlib.pyplot as plt
import time

def npy2cube(grid, start, step, cube_f):

    bohr_to_angs = 0.529177210859  # const
    start = [x / bohr_to_angs for x in start]
    step = [x / bohr_to_angs for x in step]

    with open(cube_f, 'w') as cube_file:
        cube_file.write('CPMD CUBE FILE.\n'
                        'OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z\n')
        cube_file.write("    1 %f %f %f\n" % (start[0], start[1], start[2]))
        cube_file.write("  %i    %f   0.000000     0.000000\n" %
                        (grid.shape[0], step[0]))
        cube_file.write("  %i    0.000000    %f    0.000000\n" %
                        (grid.shape[1], step[1]))
        cube_file.write("  %i    0.000000    0.000000    %f\n" %
                        (grid.shape[2], step[2]))
        cube_file.write("    1    0.000000    %f %f %f\n" %
                        (start[0], start[1], start[2]))


        i = 0
        for grid_item in np.nditer(grid, order='C'):
            if i < 5:
                cube_file.write('%f ' % (float(grid_item)))
                i += 1
            elif i == 5:
                cube_file.write('%f\n' % (float(grid_item)))
                i = 0

#n,l,m – квантовые числа, d – шаг
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.arctan2(y, x) + math.pi
    
#Боровский радиус
    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)
    
#Нормированный коэффициент
    k = lambda n,l: np.sqrt((2/n/a0)**3 * math.factorial(n-l-1)/(2*n*math.factorial(n+l)))
    
#Волновая фунция
    WF = lambda r,theta,phi,n,l,m: k(n, l) * R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
    
#Плотность вероятности WF   
    absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2
    
    result = WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)    
    return np.real(result) + np.imag(result)

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

from xmlrpc.client import ServerProxy
from pymol import cmd,stored
cmd = ServerProxy(uri='http://localhost:9123/RPC2')

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

out.write('cmd.volume_ramp_new(\'ramp007\', [\
     -0.015, 1.00, 0.00, 0.00, 0.00, \
      -0.01,  1.00, 1.00, 0.00, 0.20, \
      -0.005, 0.00, 0.00, 1.00, 0.00, \
      0.005, 0.00, 0.00, 1.00, 0.00, \
      0.01,  0.00, 1.00, 1.00, 0.20, \
      0.015, 0.00, 0.00, 1.00, 0.00, \
    ])\n')

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)
            lname = name+'.cube'
            vname = name+'_vol'
            out.write('load %s, %s\n' % (lname, name))
            out.write('volume %s, %s\n' % (vname, name))
            out.write('volume_color %s, ramp007\n' % (vname))
out.close()

cmd.load("volume_color.pml")

display(Image(filename='310.png'))
display(Image(filename='210.png'))
display(Image(filename='322.png'))

Орбитали в ORCA

with open('h.inp', 'w') as f:
    f.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);
end
* xyz 0 4
 H 0 0 0
*''')

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

out = open ('volume_color_orca.pml', 'w')

out.write('cmd.volume_ramp_new(\'ramp007\', [\
	  -0.015, 1.00, 0.00, 0.00, 0.00, \
      -0.01,  1.00, 1.00, 0.00, 0.20, \
      -0.005, 0.00, 0.00, 1.00, 0.00, \
      0.005, 0.00, 0.00, 1.00, 0.00, \
      0.01,  0.00, 1.00, 1.00, 0.20, \
      0.015, 0.00, 0.00, 1.00, 0.00, \
    ])\n')

for file in range(1,5):
    name='%s-%s' % ('H', file)
    lname = name+'.cube'
    vname = name+'_vol'
    out.write('load %s, %s\n' % (lname, name))
    out.write('volume %s, %s\n' % (vname, name))
    out.write('volume_color %s, ramp007\n' % (vname))
out.close()

display(Image(filename='H1.png'))
display(Image(filename='H2.png'))