Практикум 4

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
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
    
#Задаем боровский радиус (расстояние от ядра до ближайшей орбиты электрона) = 1
    a0 = 1
    
#Прописываем радиальную функцию R (зависящую от n, l и r) (genlaguerre - подсчет полинома Лагерра)
    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)
    
#И не забываем про нормировочный коэффициент
    norm = lambda n,l: np.sqrt((2/n/a0)**3 * math.factorial(n-l-1)/(2*n*math.factorial(n+l)))
    
#Итоговая волновая функция (sph_harm – подсчет сферической гармоники)
    WF = lambda r,theta,phi,n,l,m: norm(n, l) * R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
    
#Итоговая плотность вероятности WF будет высчитываться как квадрат модуля
#Однако если выводить именно квадрат модуля, то итоговые электронные орбитали получаются неверными 
#Опираясь на опыт людей, уже выполнивших данный практикум, было решено выводить сумму действительной и мнимой частей WF
    result = WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)    
    return np.real(result) + np.imag(result)

#Для каждого главного числа n от 1 до 3 перебираем орбитальные числа l от 0 до n-1, и для них перебираем магнитные m от -l до l
#Для каждой тройки таких чисел применяем скрипт
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 = f'{n}_{l}_{m}'
            npy2cube(grid, (-d,-d,-d), (step,step,step), f'{name}.cube')

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

import pymol
pymol.finish_launching()
from pymol import cmd,stored

cmd.reinitialize()
#Зададим цвет орбиталей (volume), например красный и зеленый
cmd.volume_ramp_new('ramp007', [-0.015, 0.00, 0.00, 0.00, 0.00, 
      -0.01,  1.00, 0.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, 0.00, 0.20, 
      0.015, 0.00, 0.00, 0.00, 0.00])

#Для каждого главного числа n от 1 до 3 перебираем орбитальные числа l от 0 до n-1, и для них перебираем магнитные m от -l до l
#Для каждой тройки таких чисел загружаем соответствующий файл cube, визуализируем его (либо срез, либо поверхность) и сохраняем
for n in range(1, 4):
    for l in range(0, n):
        for m in range(-l, l+1):
            name = f'{n}_{l}_{m}'
            cmd.load(f'{name}.cube')
            cmd.volume(f'{name}-vol', name, ramp='ramp007')
            cmd.hide()
            cmd.show('volume', f'{name}-vol')
            cmd.reset()
            cmd.turn('x', 90)
            cmd.clip('slab', 0)
            time.sleep(0.3)
            cmd.do(f'png clip_{name}, width={480}, height={270}')
for n in range(1, 4):
    for l in range(0, n):
        for m in range(-l, l+1):
            name = f'{n}_{l}_{m}'
            cmd.load(f'{name}.cube')
            cmd.volume(f'{name}-vol', name, ramp='ramp007')
            cmd.hide()
            cmd.show('volume', f'{name}-vol')
            cmd.reset()
            cmd.turn('x', 45)
            cmd.turn('y', 45)
            cmd.turn('z', 45)
            time.sleep(0.3)
            cmd.do(f'png {name}, width={480}, height={270}')

Покажем срезы орбителей

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'clip_{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.save("/home/masha/Документы/Головин/прак4/im1.png")

display(Image("/home/masha/Документы/Головин/прак4/im1.png", format='png'))

Покажем поверхности орбителей

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

img_paths = [f'{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.save("/home/masha/Документы/Головин/прак4/im2.png")

display(Image("/home/masha/Документы/Головин/прак4/im2.png", format='png'))

Чтобы проверить достоверность рассчитанных орбиталей, сравним их с выдачей программы Orca. Для этого запишем команды для подсчета первых 14 орбиталей в файл hydrogen.inp и выполним команду orca hydrogen.inp на сервере kodomo. С выдачей orca проводим все те же операции, что и с выдачей npy2cube

with open('hydrogen.inp', 'w') as f:
    f.write('''! UHF QZVPP XYZFile
%plots Format Cube
MO("H-0.cube",0,0);
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",5,0);
MO("H-6.cube",6,0);
MO("H-7.cube",7,0);
MO("H-8.cube",8,0);
MO("H-9.cube",9,0);
MO("H-10.cube",10,0);
MO("H-11.cube",11,0);
MO("H-12.cube",12,0);
MO("H-13.cube",13,0);
end
* xyz 0 4
H 0 0 0
* ''')

cmd.reinitialize()
cmd.volume_ramp_new('ramp007', [-0.015, 0.00, 0.00, 0.00, 0.00, 
      -0.01,  1.00, 0.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, 0.00, 0.20, 
      0.015, 0.00, 0.00, 0.00, 0.00])

for n in range(14):
    cmd.load(f'H-{n}.cube')
    cmd.volume(f'H-{n}-vol', f'H-{n}', ramp='ramp007')
    cmd.hide()
    cmd.show('volume', f'H-{n}-vol')
    cmd.reset()
    cmd.turn('x', 90)
    cmd.clip('slab', 0)
    time.sleep(0.5)
    cmd.do(f'png clip_orca_{n}, width={480}, height={270}')
    
for n in range(14):
    cmd.load(f'H-{n}.cube')
    cmd.volume(f'H-{n}-vol', f'H-{n}', ramp='ramp007')
    cmd.hide()
    cmd.show('volume', f'H-{n}-vol')
    cmd.reset()
    cmd.turn('x', 45)
    cmd.turn('y', 45)
    cmd.turn('z', 45)
    time.sleep(0.4)
    cmd.do(f'png orca_{n}, width={480}, height={270}')

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

img_paths = [f'clip_orca_{i}.png' for i in range(14)]
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.save("/home/masha/Документы/Головин/прак4/im3.png")

display(Image("/home/masha/Документы/Головин/прак4/im3.png", format='png'))

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

img_paths = [f'orca_{i}.png' for i in range(14)]
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.save("/home/masha/Документы/Головин/прак4/im4.png")

display(Image("/home/masha/Документы/Головин/прак4/im4.png", format='png'))

В итоге мы видим, что орбитали Orca если и отличаются от рассчитанных нами, то совсем немного. Это может быть связано с тем, что в нашем случае мы взяли сумму действительной и мнимой частей волновой функции (но квадрат модуля почему-то работать не захотел, точнее выдавал странные орбитали). Также мы приняли a0 за 1, хотя конечно эту величиную можно уточнить, и т.д. С другой стороны, наш вариант орбиталей выглядит симпатичнее)