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

Цель занятия: построить электронную плотность на основе уравнения в одно-электронном атоме и сравнить с известными программами

In [1]:
import numpy as np
import math
import scipy.special
import scipy.misc
from npy2cube import npy2cube
import time
In [ ]:
Волновая функция:
In [2]:
def w(n, l, m, d):
    """
    Вычисляет волновую функцию атома водорода в трехмерном пространстве.
    
    Параметры:
    n (int): Главное квантовое число (n >= 1)
    l (int): Орбитальное квантовое число (0 <= l < n)
    m (int): Магнитное квантовое число (-l <= m <= l)
    d (float): Размер области вычисления в пространстве [-d, d] по всем осям
    
    Возвращает:
    numpy.ndarray: Значения волновой функции в точках 3D сетки
    """
    
    # Создание 3D сетки в декартовых координатах
    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))  # Полярный угол (0 до π)
    phi = lambda x, y, z: np.arctan2(y, x)           # Азимутальный угол (0 до 2π)
    
    # Боровский радиус (в атомных единицах равен 1)
    a0 = 1.0
    
    # Радиальная часть волновой функции
    # Использует обобщенные полиномы Лагерра
    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))
    
    # Нормировочный коэффициент для волновой функции
    base = 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: base(n, l) * R(r, n, l) * scipy.special.sph_harm(m, l, phi, theta)
    
    # Для визуализации используем сумму действительной и мнимой частей,
    # что дает более точное представление орбиталей
    res = WF(r(x, y, z), theta(x, y, z), phi(x, y, z), n, l, m)
    return np.real(res) + np.imag(res)

Рассчёт значения для первых трех уровней. Функция w выдает трехмерный массив из 303030 элементов с неким шагом (или grid).

In [3]:
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,1):
            grid= w(n, l, m, d)
            name='%s_%s_%s' % (n,l,m)
            # для сохранения нужно задать координаты старта grid и шаг по каждому направлению
            npy2cube(grid, (-d,-d,-d), (step,step,step), f'{name}.cube')

В результате получились .cube файлы. Визуализируем один из них в py3Dmol:

In [4]:
import py3Dmol
view = py3Dmol.view()
alpha = open('1_0_0.cube','r').read()
view.addVolumetricData(alpha, "cube", {'isoval': 0.04, 'color': "red", 'opacity': 0.5})
view.zoomTo()
view.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

Визуализируем всё через pyMol

In [5]:
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]

import pymol
pymol.finish_launching()
from pymol import cmd,stored
In [6]:
cmd.delete('all')
cmd.volume_ramp_new('colors', [-0.015, 0.00, 0.00, 0.00, 0.00, 
      -0.01,  1.00, 0.00, 0.00, 0.20, 
      -0.005, 0.00, 1.00, 0.00, 0.00, 
      0.005, 0.00, 0.00, 1.00, 0.00, 
      0.01,  0.50, 1.00, 0.00, 0.20, 
      0.015, 0.50, 0.00, 0.50, 0.00])

names = []
for n in range(1, 4):
    for l in range(0, n):
        for m in range(-l, l+1):
            name = f'{n}_{l}_{m}'
            names.append(name)
            cmd.load(f'{name}.cube')
            cmd.volume(f'{name}-V', name, ramp='colors')
            cmd.hide()
            cmd.show('volume', f'{name}-V')
            cmd.reset()
            cmd.turn('x', 55)
            cmd.turn('y', 35)
            cmd.turn('z', 45)
            cmd.do(f'''
            png {name}.png, dpi=300
            ''')

Орбитали:

In [7]:
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

items = ['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']
orb = ['s1', 's2', 'p2', 'p2', 'p2',
            's3', 'p3', 'p3', 'p3',
            'd3', 'd3', 'd3', 'd3', 'd3']

images = [mpimg.imread(f'{num}.png') for num in items]
cols = 5  
rows = (len(images) + cols - 1) // cols  
plt.figure(figsize=(20, 4 * rows)) 
for i, (img, o) in enumerate(zip(images, orb), 1):
    plt.subplot(rows, cols, i)
    plt.imshow(img)
    plt.title(o, fontsize=12, pad=10)
    plt.axis('off')

plt.tight_layout(pad=2.0, h_pad=3.0, w_pad=3.0)  
plt.show()
No description has been provided for this image

Работа с программой psi4 осуществлялась в другом ноутбуке на сервере со следующим кодом

import psi4

import numpy as np

psi4.core.set_output_file('output.dat')

Зададим геометрию

psi4.core.set_output_file('output.dat')

g = '''

H

'''

Расчитаем энергию b орбитали

m = psi4.geometry(g)

psi4.set_options({"maxiter": 200, "fail_on_maxiter" : True, 'reference':'uhf'})

ener, wave=psi4.energy('scf/cc-pvtz', return_wfn = True, molecule = m )

psi4.cubeprop(wave)

Примеры полученных файлов ниже

Для дальнейшей работы файлы были скачаны в рабочую дирректорию этого ноутбука.

In [9]:
import os
all=[i for i in os.listdir() if 'Psi_a' in i]
In [10]:
all
Out[10]:
['Psi_a_10_4-Ag.cube',
 'Psi_a_11_5-Ag.cube',
 'Psi_a_12_2-B1u.cube',
 'Psi_a_13_2-B3u.cube',
 'Psi_a_14_2-B2u.cube',
 'Psi_a_1_1-Ag.cube',
 'Psi_a_2_2-Ag.cube',
 'Psi_a_3_1-B1u.cube',
 'Psi_a_4_1-B2u.cube',
 'Psi_a_5_1-B3u.cube',
 'Psi_a_6_3-Ag.cube',
 'Psi_a_7_1-B2g.cube',
 'Psi_a_8_1-B1g.cube',
 'Psi_a_9_1-B3g.cube']

Ниже визуализация и сравнение

In [11]:
cmd.delete('all')
cmd.volume_ramp_new('colors', [-0.015, 0.00, 0.00, 0.00, 0.00, 
      -0.01,  1.00, 0.00, 0.00, 0.20, 
      -0.005, 0.00, 1.00, 0.00, 0.00, 
      0.005, 0.00, 0.00, 1.00, 0.00, 
      0.01,  0.50, 1.00, 0.00, 0.20, 
      0.015, 0.50, 0.00, 0.50, 0.00])

names = []
for n in all:
    name=n.split('.')[0]
    cmd.load(n)
    cmd.volume(f'{name}_vol', name, ramp='colors')
    cmd.hide()
    cmd.show('volume', f'{name}_vol')
    cmd.reset()
    cmd.turn('x', 55)
    cmd.turn('y', 35)
    cmd.turn('z', 45)
    cmd.do(f'''
    png {name}.png, dpi=300
    ''')
In [12]:
import glob
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

pngs=[i.split('.')[0]+'.png' for i in all]
gets = [mpimg.imread(img_path) for img_path in pngs]

columns = 5
rows = (len(gets) + columns - 1) // columns 
f = plt.figure(figsize=(30, 5*rows)) 
for i, image in enumerate(gets):
    f.add_subplot(rows, columns, i + 1)
    plt.imshow(image)
    plt.axis('off') 
    plt.tight_layout() 

plt.show()
No description has been provided for this image

В целом похоже, но есть странности.

Некоторые расчётные орбитали отличаются:

2p-орбитали имеют провалы электронной плотности,

3s-орбитали демонстрируют избыточную плотность в неожиданных областях