{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Домашнее задание 4\n", "**Выполнил: Ким Адамейко, группа мАДБМ16**" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "! wget -nv -nc http://kodomo.fbb.msu.ru/~golovin/ipynb/npy2cube.py\n", "from npy2cube import *" ] }, { "cell_type": "code", "execution_count": 186, "metadata": {}, "outputs": [], "source": [ "import time\n", "import __main__\n", "import os\n", "\n", "import numpy as np\n", "import scipy.special\n", "import scipy.misc\n", "\n", "from IPython.display import Image, display, HTML\n", "from __future__ import print_function, division\n", "\n", "from xmlrpclib import ServerProxy\n", "cmd = ServerProxy(uri=\"http://localhost:9123/RPC2\")\n", "cmd.do(\"cd \" + os.getcwd())\n", "\n", "def img_row(img_data): # tuples of (filename, title)\n", " img_files, titles = zip(*img_data)\n", " display(HTML(\"\" + \n", " ('' + ''.join(['' \n", " for title in titles]) + '' if titles else '') +\n", " '' + ''.join(['' \n", " for img in img_files])+'
' + title + '
'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Часть 1. Волновая функция для атома водорода\n", "\n", "В полярных координатах $(r,\\vartheta, \\varphi)$:\n", "\n", "$${\\large \\psi_{n\\,l\\,m}(r,\\vartheta, \\varphi)} = \\sqrt{\\left(\\dfrac{2}{n a_{0}}\\right)^{\\!\\!3} \\dfrac{(n-l-1)!}{2n\\cdot(n+l)!}}\\ {\\large e^{-\\frac{p}{2}}\\ p^{l}\\ L^{2l+1}_{n-l-1}(p)\\ Y^{m}_{l}(\\vartheta, \\varphi),\\quad p = \\dfrac{2r}{na_{0}}}$$ \n", "где:\n", "\n", "* $L^{2l+1}_{n-l-1}(p) -$ [обобщенный полином Лагерра](https://en.wikipedia.org/wiki/Laguerre_polynomials#Generalized_Laguerre_polynomials) степени $n-l-1;$\n", "* $ Y^{m}_{l}(\\vartheta, \\varphi) -$ [сферическая гармоника](https://en.wikipedia.org/wiki/Spherical_harmonics) степени $l$ и порядка $m$;\n", "* $a_0 = \\dfrac{4 \\pi \\varepsilon_0 \\hbar^2}{\\mu e^2} -$ приведенный [радиус Бора](https://en.wikipedia.org/wiki/Bohr_radius);\n", " * $\\varepsilon_0 - $ электрическая постоянная;\n", " * $e - $ элементарный заряд;\n", "\n", "\n", "* $n, l, m\\ - $ основные квантовые числа\n", " * $n$ – основное число $(1,2,3, \\dots)$\n", " * $l$ – орбитальное число $(0,1,2,\\dots n-1)$\n", " * $m$ – магнитное число $(-l, \\dots, +l)$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Часть 2. Ручной расчёт орбиталей \n", "Определим необходимые функции для расчётов и воспользуемся функцией() `npy2cube` из соответствующего файла, предоставленного нам заранее:" ] }, { "cell_type": "code", "execution_count": 229, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def w(n, l, m, d):\n", " \n", " # Неразреженный 3D массив точек, величина компл.ч. (30j) задаёт число точек между границами (такой синтаксис)\n", " x,y,z = np.mgrid[-d:d:30j, -d:d:30j, -d:d:30j]\n", " \n", " # Переход к сферическим координатам\n", " r = lambda x,y,z: np.sqrt(x**2+y**2+z**2)\n", " theta = lambda x,y,z: np.arccos(z/r(x,y,z))\n", " phi = lambda x,y,z: np.arctan(y/x)\n", " \n", " a0 = 1.0\n", " \n", " # Радиальная часть волновой функции\n", " R = lambda r,n,l: (2.0*r/n/a0)**l * np.exp(-r/n/a0) * scipy.special.genlaguerre(n-l-1,2*l+1)(2.0*r/n/a0)\n", " # Умножение радиальной части на угловую\n", " WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)\n", " # Абсолютное значение волновой функции (плотность вероятности нахождения электрона в конкретном состоянии)\n", " absWF = lambda r,theta,phi,n,l,m: np.absolute(WF(r,theta,phi,n,l,m))**2\n", " \n", " return WF(r(x,y,z), theta(x,y,z), phi(x,y,z), n, l, m)" ] }, { "cell_type": "code", "execution_count": 236, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1-0-0.cube 2-1-0.cube\t3-0-0.cube 3-1-1.cube\t3-2-1.cube\r\n", "2-0-0.cube 2-1-1.cube\t3-1-0.cube 3-2-0.cube\t3-2-2.cube\r\n" ] } ], "source": [ "d = 30\n", "step= float(2.*d/29)\n", "\n", "cdir = './cube/'\n", "! mkdir -p {cdir}\n", "qns = []\n", "for n in range(1, 4):\n", " for l in range(0, n):\n", " for m in range(0, l+1):\n", " data = w(n, l, m, d)\n", " qns.append('%d-%d-%d' % (n, l, m))\n", " filename = '%s%d-%d-%d.cube' % (cdir, n, l, m)\n", " npy2cube(data, (-d,-d, -d), (step,step,step), filename) # парам: коорд. \"самого отриц.точки\" куба и шаг\n", "\n", "!ls {cdir}" ] }, { "cell_type": "code", "execution_count": 237, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " CPMD CUBE FILE.\r\n", "OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z\r\n", " 1 -56.691784 -56.691784 -56.691784\r\n", " 30 3.909778 0.000000 0.000000\r\n", " 30 0.000000 3.909778 0.000000\r\n", " 30 0.000000 0.000000 3.909778\r\n", " 1 0.000000 -56.691784 -56.691784 -56.691784\r\n", "-0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000\r\n", "-0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000\r\n", "-0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000\r\n" ] } ], "source": [ "!head {cdir}2-1-1.cube" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Сделаем необходимые настройки и осуществим рендеринг объемов с помощью Pymol (на локальном компьютере)" ] }, { "cell_type": "code", "execution_count": 269, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# настройка цветового градиента, для правильного отображения он должен быть симметричен относительно нуля\n", "# числа: x-координата, RGB-код цвета, прозрачность\n", "cmd.do( \n", "\"\"\"\n", "cmd.volume_ramp_new('ramp', [\\\n", " -0.05, 0.36, 0.70, 1.00, 0.0, \\\n", " -0.03, 0.50, 0.00, 1.00, 0.05, \\\n", " -0.02, 0.92, 0.00, 0.80, 0.00, \\\n", " 0.02, 0.92, 0.00, 0.80, 0.00, \\\n", " 0.03, 0.50, 0.00, 1.00, 0.05, \\\n", " 0.05, 0.36, 0.70, 1.00, 0.0, \\\n", " ])\n", "\"\"\")" ] }, { "cell_type": "code", "execution_count": 270, "metadata": {}, "outputs": [], "source": [ "imgdir = './images/manualcomp/'\n", "!mkdir -p {imgdir}\n", "\n", "for qn in qns:\n", " cmd.reinitialize()\n", " cmd.do('load ./cube/{0}.cube, {0}'.format(qn))\n", " cmd.do('volume {0}_vol, {0}'.format(qn))\n", " cmd.do('volume_color {0}_vol, ramp'.format(qn))\n", " time.sleep(0.5)\n", " if qn in ['2-1-0', '3-1-0', '3-2-0']: # некоторые картинки нужно повернуть для лучшего отражения ситуации\n", " cmd.turn('y', 90)\n", " cmd.turn('x', -45)\n", " cmd.do('png {0}{1}.png, width={2}, height={3}'.format(imgdir, qn, 400, 300))\n", " time.sleep(0.5)\n", " cmd.delete('all')" ] }, { "cell_type": "code", "execution_count": 271, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
1-0-02-0-0
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
2-1-02-1-1
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
3-0-03-1-0
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
3-1-13-2-0
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
3-2-13-2-2
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "for i in range(0, len(qns), 2):\n", " img_row(zip(map(lambda x: imgdir + x, qns[i:i+2]), qns[i:i+2]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Часть 3. Расчёт орбиталей в ORCA\n", "Следующие две ячейки запускались на сервере. Затем результаты расчётов в виде cube-файлов архивировались и загружались на локальный компьютер, где снова отрисовывались в Pymol. К сожалению, удалось получить корректные изображения только для первых четырёх конфигураций квантовых чисел. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%%writefile h.in\n", "! UHF SVP XYZFile\n", "%plots Format Cube\n", " MO(\"./cube-orca/1-0-0.cube\",0,0);\n", " MO(\"./cube-orca/2-0-0.cube\",1,0);\n", " MO(\"./cube-orca/2-1-0.cube\",2,0);\n", " MO(\"./cube-orca/2-1-1.cube\",3,0);\n", "end\n", "* xyz 0 4\n", " H 0 0 0\n", "*" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import os\n", "curpath = os.getcwd()\n", "! mkdir -p './cube-orca' \n", "! export PATH=${PATH}:/{curpath}\n", "! orca h.in > h.out\n", "! tail h.out -n 2\n", "! tar -zcf cube-orca.tar.gz cube-orca" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Следующие ячейки запускались на локальном компьютере" ] }, { "cell_type": "code", "execution_count": 212, "metadata": {}, "outputs": [], "source": [ "cdir2 = './cube-orca/'\n", "! mkdir -p {cdir2}\n", "! tar -zxf cube-orca.tar.gz -C {cdir2} --strip-components=1" ] }, { "cell_type": "code", "execution_count": 282, "metadata": {}, "outputs": [], "source": [ "imgdir2 = './images/orca/'\n", "!mkdir -p {imgdir2}\n", "\n", "for qn in qns[:4]:\n", " cmd.reinitialize()\n", " cmd.do('load ./cube-orca/{0}.cube, {0}'.format(qn))\n", " cmd.do('volume {0}_vol, {0}'.format(qn))\n", " cmd.do('volume_color {0}_vol, ramp'.format(qn))\n", " time.sleep(0.5)\n", " if qn in ['2-1-0', '2-1-1']:\n", " cmd.turn('y', 90)\n", " cmd.turn('x', 45)\n", " cmd.do('png {0}{1}.png, width={2}, height={3}'.format(imgdir2, qn, 400, 300))\n", " time.sleep(0.5)\n", " cmd.delete('all')" ] }, { "cell_type": "code", "execution_count": 283, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
1-0-02-0-0
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
2-1-02-1-1
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "for i in range(0, 4, 2):\n", " img_row(zip(map(lambda x: imgdir2 + x, qns[i:i+2]), qns[i:i+2]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Наблюдается некоторое сходство с картинками, полученными с помощью \"ручных\" расчётов. Также заметна разница в масштабах. Возможно, это можно настроить с помощью каких-либо параметров ORCA, но я не успел разобраться в этом." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Источники \n", "\n", "1. [Wave function: Hydrogen atom -- Wikipedia](https://en.wikipedia.org/wiki/Wave_function#Hydrogen_atom)\n", "1. [Volume color -- PymolWiki](https://pymolwiki.org/index.php/Volume_color)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.14" } }, "nbformat": 4, "nbformat_minor": 2 }