{ "cells": [ { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "#### Hydrogen atomic orbitals: Part 2\n", "In this notebook we will introduce 2D and 3D plotting. Let's start by copying the function for evaluating the hydrogen wavefunction from the previous notebook." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import scipy as sp\n", "import scipy.special as sps\n", "import matplotlib.pyplot as plt\n", "\n", "def h_orbital(r,theta,phi,n,l,m):\n", " pf = ( (2./n)**3. * sps.factorial(n-l-1) / (2. * n * sps.factorial(n+l)) )**0.5\n", " return pf * np.exp(-r/2.) *r**l * sps.sph_harm(m,l,phi,theta) * sps.eval_genlaguerre(n-l-1,2*l+1,r)" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "#### Conversion of Cartesian to spherical coordinates\n", "\n", "Before we get started on making 2D and 3D plots, we need a function that can convert Cartesian coordinates into spherical coordinates. Using trigonometric functions, it is easy to derive the expressions for the spherical coordinates from $x, y, z$ as follows:\n", "\n", "$$ r = \\sqrt{x^2 + y^2 + z^2} ; \\quad \\theta = \\cos^{-1} \\frac{z}{r} ; \\quad\\phi = \\tan^{-1} \\frac{y}{x} $$\n", "\n", "However, when applying these formulas there is quickly a problem. Notice that for $\\theta$, there will be a divide-by-zero error at the origin. Even worse, the formula for $\\phi$ returns the same value for $(-x, -y)$ as for $(x, y)$, even though we know the actual $\\phi$ angle differs by $\\pi$ for these two points. \n", "\n", "Both problems can be solved by using the `np.arctan2(y, x)` function. This function takes in two arguments that would normally be divided in the conventional arctangent, thus giving it the ability to return values in the full range of $\\phi$. Also, this function returns standard values when one or both inputs are equal to zero; if we had used `np.arctan`, this would have resulted in errors. \n", "\n", "The boardwork below shows how `np.arctan2` is used in place of `np.arccos` and `np.arcsin` to give a robust conversion from Cartesian to spherical coordinates.\n", "\n", "
\n", "\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def to_spherical(x, y, z):\n", " \"\"\"Converts the Cartesian coordinates x,y,z to r,theta,phi\n", "\n", " This function expects x, y, and z to be 3D arrays,\n", " such as those generated by np.meshgrid. The output arrays are\n", " always dense 3D grids, and correspond to r, theta, and phi\n", "\n", " Parameters\n", " ----------\n", " x, y, z : np.ndarray\n", " Cartesian x / y / z values, in a.u.\n", "\n", " Returns\n", " -------\n", " tuple of 3x np.ndarrays\n", " 3D grid of r, theta, phi values\n", " \"\"\"\n", "\n", " r = ( x**2 + y**2 + z**2 )**0.5\n", " theta = np.arctan2((x**2 + y**2)**0.5, z)\n", " phi = np.arctan2(y, x)\n", "\n", " return r, theta, phi" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "#### Plotting heat maps of orbital slices in 2D\n", "\n", "In what follows, we will create a 2-D grid where two out of the three $(x, y, z)$ coordinates vary and the other coordinate is set to a constant. Note that we still have three input arrays for `to_spherical()` because all three Cartesian coordinates are needed to calculate the corresponding spherical coordinates. The shape of the 3-D arrays in the following cells is (101, 1, 101). The spherical coordinates are then used as input to the orbital function." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create range of values along each axis.\n", "x = np.linspace(-10, 10, 101)\n", "y = 0.0\n", "z = np.linspace(-10, 10, 101)\n", "\n", "# Create 3D arrays containing x, y, z-coords of each point.\n", "# The 'indexing' keyword argument is needed, otherwise the \n", "# ordering of our dimensions will be messed up.\n", "xg, yg, zg = np.meshgrid(x, y, z, indexing='ij')\n", "\n", "print(\"The shape of xg is: \", xg.shape)\n", "\n", "# Calculate r, theta, and phi values at each grid point.\n", "rg, tg, pg = to_spherical(xg, yg, zg)\n", "\n", "# Calculate the hydrogen orbital wavefunction.\n", "psi_3p = h_orbital(rg,tg,pg,3,1,0)\n", "\n", "# The output of h_orbital is a complex type, even though the chosen \n", "# 3pz orbital is real (imaginary part is zero). Matplotlib can only plot\n", "# real numbers, so we take the real part here.\n", "psi_3p = psi_3p.real" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "##### Choosing a colormap for plotting\n", "\n", "The following cell is a customization of the color map before starting to draw plots in 2D. Blue and red are used for oppositely signed values. The full list of color maps is available here:\n", "\n", "[Colormaps in Matplotlib by name](https://matplotlib.org/stable/tutorials/colors/colormaps.html)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib \n", "matplotlib.rc('image', cmap='seismic')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Plotting the heat map\n", "\n", "To make the heat map, we are using [`Axes.pcolormesh`](https://matplotlib.org/api/_as_gen/matplotlib.axes.Axes.pcolormesh.html#matplotlib.axes.Axes.pcolormesh). Its main function signature is `pcolormesh(X,Y,C)` where the arguments are the x values, the y values, and the data to be plotted. By passing the same arguments into `Axes.contour()` we can also create level curves.\n", "\n", "Our coordinate and function data are stored in 3D arrays, but `pcolormesh()` can only plot 2D data. Therefore we need to use indexing to extract a 2D slice out of each 3D array. This is done using the syntax `xg[:,0,:], zg[:,0,:], psi_3p[:,0,:]` to extract all values in the $x$ and $z$ dimensions but keep only the zeroth element for $y$. (In this example we only had one value in the $y$ dimension, but reducing the array dimension is still essential.)\n", "\n", "We also called `ax.set_aspect('equal', 'box')`. With this, the size and extent of the horizontal and vertical axes will be set equal and the plot area will be made square. Otherwise the plot will appear squashed in one direction or another. This comes in handy when both axes correspond to Cartesian coordinates, and various other situations.\n", "\n", "Finally, a color scale is useful for knowing the correspondence between color and function value. The color scale is the return value of `pcolormesh()`, and it is passed to `Figure.colorbar()` to add a special Axes object containing the color bar." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create a Figure and Axes object\n", "fig, ax = plt.subplots()\n", "# This sets the sets the x- and y- axis limits to be equal\n", "# and makes the plotting area square. \n", "ax.set_aspect('equal', 'box')\n", "cbar = ax.pcolormesh(xg[:,0,:], zg[:,0,:], psi_3p[:,0,:])\n", "ax.contour(xg[:,0,:], zg[:,0,:], psi_3p[:,0,:], colors='k', \n", " linewidths=1, levels=[-0.01, 0.0, 0.01])\n", "fig.colorbar(cbar)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Plotting several heat maps in a single figure\n", "\n", "Here the calculations and plotting codes are organized into a function, making it easy to create different plots by changing the quantum numbers. Four `Axes` objects are created in a 2x2 grid, making it easy to compare the $xz-$ cross sections of orbitals that differ by only one quantum number; here the quantum numbers are $(4, 0, 0); (4, 1, 0); (4, 2, 0); (4, 3, 0)$. These correspond to the $4s, 4p_z, 4d_{z^2}$ and $4f_{z^3}$ orbitals.\n", "\n", "In these plots we also used the keyword arguments `vmin, vmax` to manually set the color scale (otherwise it would not always be centered around zero). Normally you need to know the range of the values before setting a reasonable color scale. The contour lines are drawn only at $\\psi = 0.0$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_orbital_xz(x, y, z, n, l, m, fig, ax):\n", " xg, yg, zg = np.meshgrid(x, y, z, indexing='ij')\n", " rg, tg, pg = to_spherical(xg, yg, zg)\n", " psi = h_orbital(rg,tg,pg,n,l,m)\n", " psi = psi.real\n", " ax.set_aspect('equal', 'box')\n", " ax.pcolormesh(xg[:,0,:], zg[:,0,:], psi[:,0,:], vmin=-0.02, vmax=0.02)\n", " ax.contour(xg[:,0,:], zg[:,0,:], psi[:,0,:], colors='k', \n", " linewidths=1, levels=[0.0])\n", " #fig.colorbar(cbar)\n", " \n", "fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2)\n", "x = np.linspace(-10, 10, 101)\n", "y = 0.0\n", "z = np.linspace(-10, 10, 101)\n", "plot_orbital_xz(x, y, z, 4, 0, 0, fig, ax1)\n", "plot_orbital_xz(x, y, z, 4, 1, 0, fig, ax2)\n", "plot_orbital_xz(x, y, z, 4, 2, 0, fig, ax3)\n", "plot_orbital_xz(x, y, z, 4, 3, 0, fig, ax4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Plotting the real and imaginary part\n", "\n", "We can also plot orbital cross sections in different planes. Here I've chosen the $xy-$ plane, and chosen the quantum numbers $(3, 2, 2)$. This is the first time we are using nonzero $m$, and these orbitals are complex. The code below plots the real and imaginary part separately. The real part is proportional to $3d_{x^2-y^2}$ and the imaginary part is proportional to $3d_{xy}$. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, (ax1, ax2) = plt.subplots(1,2)\n", "x = np.linspace(-10, 10, 101)\n", "y = np.linspace(-10, 10, 101)\n", "z = 0.0\n", "n = 3\n", "l = 2\n", "m = 2\n", "\n", "xg, yg, zg = np.meshgrid(x, y, z, indexing='ij')\n", "rg, tg, pg = to_spherical(xg, yg, zg)\n", "psi = h_orbital(rg,tg,pg,n,l,m)\n", "\n", "ax1.set_aspect('equal', 'box')\n", "ax1.pcolormesh(xg[:,:,0], yg[:,:,0], psi[:,:,0].real, vmin=-0.02, vmax=0.02)\n", "ax1.contour(xg[:,:,0], yg[:,:,0], psi[:,:,0].real, colors='k', \n", " linewidths=1, levels=[0.0])\n", "\n", "ax2.set_aspect('equal', 'box')\n", "ax2.pcolormesh(xg[:,:,0], yg[:,:,0], psi[:,:,0].imag, vmin=-0.02, vmax=0.02)\n", "ax2.contour(xg[:,:,0], yg[:,:,0], psi[:,:,0].imag, colors='k', \n", " linewidths=1, levels=[0.0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Interactive 3D plots of isosurfaces\n", "\n", "Isosurfaces are not a standard feature in current versions of Matplotlib. Although it is possible to do, one needs to separately calculate the isosurface, and then pass it to Matplotlib as a collection of triangular faces. On the other hand, the plotly package does support isosurfaces, and the resulting plots are interactive to boot. Here is an example code for creating a 3D isosurface in Plotly.\n", "\n", "Check out the [plotly `Isosurface` documentation](https://plotly.com/python-api-reference/generated/plotly.graph_objects.Isosurface.html). The most important keyword arguments are `x`, `y`, `z`, and `value`, representing the coordinates and function value on a grid of points. Unlike Matplotlib's `pcolormesh`, the coordinate arrays and values are all expected to be 1-D arrays, which can be easily done using the `array.flatten()` method. This is because the `Isosurface` class does not require points to be on a regular 3-D grid.\n", "\n", "Other arguments to `Isosurface` include specifying the color scale (`colorscale`), and plotting multiple surfaces by specifying a minimum and maximum value (`isomin` and `isomax`), as well as how many isosurfaces we want (`surface_count`). The `caps` allows the edges of the 3D plot to show the interior of the isosurface (see [Matlab's documentation](https://www.mathworks.com/help/matlab/visualize/isocaps-add-context-to-visualizations.html) for a conceptual explanation, because plotly's documentation is unclear). Here we simplify things by turning them off." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import plotly.graph_objects as go\n", "\n", "x = np.linspace(-15, 15, 51)\n", "y = np.linspace(-15, 15, 51)\n", "z = np.linspace(-15, 15, 51)\n", "n = 3\n", "l = 2\n", "m = 0\n", "\n", "xg, yg, zg = np.meshgrid(x, y, z, indexing='ij')\n", "rg, tg, pg = to_spherical(xg, yg, zg)\n", "psi = h_orbital(rg,tg,pg,n,l,m)\n", "psi = psi.real\n", "\n", "psimax = psi.max()\n", "\n", "iso = go.Isosurface(x=xg.flatten(), y=yg.flatten(), z=zg.flatten(),\n", " value=psi.flatten(), colorscale='RdBu',\n", " isomin=-.75*psimax, isomax=.75*psimax,\n", " surface_count=7,opacity=0.5, caps=dict(x_show=False,y_show=False,z_show=False))\n", "\n", "fig= go.Figure(data=iso)\n", "\n", "# I thought the default behavior of using scroll to zoom was distracting so I turned it off.\n", "fig.show(config={'scrollZoom': False})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Baby steps toward molecular orbitals\n", "\n", "The creation of molecular orbitals from linear combinations of atomic orbitals (LCAO-MO) is a foundation of modern quantum chemistry. At its core is a very simple concept: MOs are formed by taking linear combinations of AOs - essentially hydrogen orbitals - centered on different atoms in the molecule. \n", "\n", "The math required to describe a MO is a simple extension of the hydrogen orbitals. We need a way to move the AO center to different locations, and this could simply be done by adding an `origin` keyword argument to the `to_spherical()` function that computes the spherical coordinates. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def to_spherical(x, y, z, origin=np.array([0, 0, 0])):\n", " \"\"\"Converts the Cartesian coordinates x,y,z to r,theta,phi\n", "\n", " This function expects x, y, and z to be 3D arrays,\n", " such as those generated by np.meshgrid. The output arrays are\n", " always dense 3D grids, and correspond to r, theta, and phi\n", "\n", " Parameters\n", " ----------\n", " x, y, z : np.ndarray\n", " Cartesian x / y / z values, in a.u.\n", "\n", " Returns\n", " -------\n", " tuple of 3x np.ndarrays\n", " 3D grid of r, theta, phi values\n", " \"\"\"\n", " dx = x - origin[0]\n", " dy = y - origin[1]\n", " dz = z - origin[2]\n", " \n", " r = ( dx**2 + dy**2 + dz**2 )**0.5\n", " theta = np.arctan2((dx**2 + dy**2)**0.5, dz)\n", " phi = np.arctan2(dy, dx)\n", "\n", " return r, theta, phi\n", "\n", "# Convenience function so we don't have to enter \n", "# the keyword arguments repeatedly. Setting an odd number \n", "# of surfaces will result in showing the nodal surfaces (isovalue 0).\n", "def plotly_isosurface(xg, yg, zg, psi):\n", " psimax = psi.max()\n", " return go.Isosurface(x=xg.flatten(), y=yg.flatten(), z=zg.flatten(),\n", " value=psi.flatten(), colorscale='RdBu',\n", " isomin=-.75*psimax, isomax=.75*psimax,\n", " surface_count=7,opacity=0.5, \n", " caps=dict(x_show=False,y_show=False,z_show=False))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's an orbital plot with the origin shifted to $(0, 0, -5)$ on the $z-$axis. Again we are taking the real part of the orbital, which is proportional to the $3d_{yz}$ orbital in this case. More rigorously, the $3d_{xz}$ and $3d_{yz}$ orbitals can be obtained by taking linear combinations of the complex $3d$ orbitals with $m=-1$ and $m=+1$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "rg, tg, pg = to_spherical(xg, yg, zg, origin=np.array([0, 0, -5]))\n", "psi = h_orbital(rg,tg,pg,3,2,1)\n", "psi = psi.real\n", "\n", "iso = plotly_isosurface(xg, yg, zg, psi)\n", "\n", "fig= go.Figure(data=iso)\n", "fig.show(config={'scrollZoom': False})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By adding two orbitals with different centers in an in-phase way, we can make a visualization of a bonding molecular orbital. In actual QM calculations, the atomic orbital sizes would be affected by the nuclear charge and core electron shielding, and sums of Gaussian functions are typically used to approximate the exponential decay of the orbital with increasing $r$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "rg1, tg1, pg1 = to_spherical(xg, yg, zg, origin=np.array([0, 0, -5]))\n", "psi1 = h_orbital(rg1,tg1,pg1,3,2,1)\n", "psi1 = psi1.real\n", "\n", "rg2, tg2, pg2 = to_spherical(xg, yg, zg, origin=np.array([0, 0, 6]))\n", "psi2 = h_orbital(rg2,tg2,pg2,2,1,1)\n", "psi2 = psi2.real * 0.5 # You can tune this constant to see how it affects the orbital shapes.\n", "\n", "psi = psi1 + psi2\n", "\n", "iso = plotly_isosurface(xg, yg, zg, psi)\n", "\n", "fig= go.Figure(data=iso)\n", "fig.show(config={'scrollZoom': False})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Does plotly support having multiple subplots in the same figure? Yes, it does. Here we draw the two AOs that contribute to form the bonding and antibonding MO." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from plotly.subplots import make_subplots\n", "# Initialize figure with 4 3D subplots\n", "fig = make_subplots(\n", " rows=1, cols=4,\n", " specs=[[{'type': 'surface'}, {'type': 'surface'},\n", " {'type': 'surface'}, {'type': 'surface'}]])\n", "\n", "rg1, tg1, pg1 = to_spherical(xg, yg, zg, origin=np.array([0, 0, -4]))\n", "psi1 = h_orbital(rg1,tg1,pg1,3,2,1)\n", "psi1 = psi1.real\n", "\n", "rg2, tg2, pg2 = to_spherical(xg, yg, zg, origin=np.array([0, 0, 4]))\n", "psi2 = h_orbital(rg2,tg2,pg2,2,1,1)\n", "psi2 = psi2.real * 0.5\n", "\n", "psi_bond = psi1 + psi2\n", "psi_anti = psi1 - psi2\n", "\n", "iso1 = plotly_isosurface(xg, yg, zg, psi1)\n", "iso2 = plotly_isosurface(xg, yg, zg, psi2)\n", "iso_bond = plotly_isosurface(xg, yg, zg, psi_bond)\n", "iso_anti = plotly_isosurface(xg, yg, zg, psi_anti)\n", "\n", "# adding surfaces to subplots.\n", "fig.add_trace(iso1, row=1, col=1)\n", "fig.add_trace(iso2, row=1, col=2)\n", "fig.add_trace(iso_bond, row=1, col=3)\n", "fig.add_trace(iso_anti, row=1, col=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are many possible ways to extend this code. You can write a function that makes it convenient to draw multiple orbitals, each having a different set of three quantum numbers and a different origin. Ultimately you could make visualizations of the MO diagrams that you see in organic and inorganic chemistry textbooks. The possibilities are endless!" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.18" } }, "nbformat": 4, "nbformat_minor": 4 }