diff --git a/Tutorials/14_Visualization/vizualize.ipynb b/Tutorials/14_Visualization/vizualize.ipynb new file mode 100644 index 00000000..268faae5 --- /dev/null +++ b/Tutorials/14_Visualization/vizualize.ipynb @@ -0,0 +1,1154 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Visualization\n", + "\n", + "Psi4 can produce cube files containing various properties. Here's how to visualize the information in Jupyter." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "%matplotlib notebook\n", + "\n", + "import numpy as np\n", + "import psi4\n", + "from skimage import measure\n", + "import matplotlib.pyplot as plt\n", + "from mpl_toolkits.mplot3d import Axes3D\n", + "from ipywidgets import interact, interactive, widgets" + ] + }, + { + "cell_type": "code", + "execution_count": 65, + "metadata": {}, + "outputs": [], + "source": [ + "psi4.set_memory(int(5e8))\n", + "numpy_memory = 2\n", + "\n", + "psi4.core.set_output_file('output.dat', False)\n", + "\n", + "mol = psi4.geometry(\"\"\"\n", + "1 2\n", + "O\n", + "H 1 1.1\n", + "H 1 1.1 2 104\n", + "\"\"\")\n", + "\n", + "# Set computation options\n", + "psi4.set_options({'basis': 'cc-pvdz',\n", + " 'reference' : 'uhf',\n", + " 'e_convergence': 1e-8})\n", + "\n", + "E,wfn = psi4.energy('scf', return_wfn=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 85, + "metadata": {}, + "outputs": [], + "source": [ + "class Visualizer:\n", + " def parse_cube(self, filename):\n", + " \"\"\" Parses a cube file, returning a dict of the information contained.\n", + " The cubefile itself is stored in a numpy array. \"\"\"\n", + " with open(filename) as fp:\n", + " results = {}\n", + "\n", + " # skip over the title\n", + " fp.readline()\n", + " fp.readline()\n", + "\n", + " origin = fp.readline().split()\n", + " natoms = int(origin[0])\n", + " results['minx'] = minx = float(origin[1])\n", + " results['miny'] = miny = float(origin[2])\n", + " results['minz'] = minz = float(origin[3])\n", + "\n", + " infox = fp.readline().split()\n", + " numx = int(infox[0])\n", + " incx = float(infox[1])\n", + " results['incx'] = incx\n", + " results['numx'] = numx\n", + " results['maxx'] = minx + incx * numx\n", + "\n", + " infoy = fp.readline().split()\n", + " numy = int(infoy[0])\n", + " incy = float(infoy[2])\n", + " results['incy'] = incy\n", + " results['numy'] = numy\n", + " results['maxy'] = miny + incy * numy\n", + "\n", + " infoz = fp.readline().split()\n", + " numz = int(infoz[0])\n", + " incz = float(infoz[3])\n", + " results['incz'] = incz\n", + " results['numz'] = numz\n", + " results['maxz'] = minz + incz * numz\n", + "\n", + " atnums = []\n", + " coords = []\n", + " for atom in range(natoms):\n", + " coordinfo = fp.readline().split()\n", + " atnums.append(int(coordinfo[0]))\n", + " coords.append(list(map(float, coordinfo[2:])))\n", + " results['atom_numbers'] = np.array(atnums)\n", + " results['atom_coords'] = np.array(coords)\n", + "\n", + " data = np.array([ float(entry) for line in fp for entry in line.split() ])\n", + " if len(data) != numx*numy*numz:\n", + " raise Exception(\"Amount of parsed data is inconsistent with header in Cube file!\")\n", + " results['data'] = data.reshape((numx,numy,numz))\n", + "\n", + " return results\n", + "\n", + " \n", + " def run_cubeprop(self):\n", + " import tempfile, shutil\n", + " tmpdir = tempfile.mkdtemp()\n", + " psi4.set_options({'cubeprop_filepath' : tmpdir,\n", + " 'cubic_grid_spacing' : [self.spacing]*3})\n", + " psi4.cubeprop(self.wfn)\n", + " if self.type_cubefile[0] == \"potential\":\n", + " data = self.parse_cube(tmpdir + \"/Dt.cube\")\n", + " espdata = self.parse_cube(tmpdir + \"/\" + self.type_cubefile[1])\n", + " else:\n", + " data = self.parse_cube(tmpdir + \"/\" + self.type_cubefile[1])\n", + " espdata = None\n", + " shutil.rmtree(tmpdir)\n", + " return data, espdata\n", + "\n", + " \n", + " def process_cubefile(self):\n", + " if self.type_cubefile[0] == \"orbital\":\n", + " import re\n", + " orbitalre = re.compile(r'Psi_([ab])_(\\d+)_.*.cube')\n", + " matchobj = orbitalre.match(self.type_cubefile[1])\n", + " if not matchobj:\n", + " raise Exception(\"Unable to parse orbital cube file name.\")\n", + " orbnum = int(matchobj.group(2))\n", + " if matchobj.group(1) == 'b': orbnum = -orbnum\n", + " psi4.set_options({'cubeprop_tasks' : ['orbitals']})\n", + " psi4.set_options({'cubeprop_orbitals' : [orbnum]})\n", + " elif self.type_cubefile[0] == \"density\":\n", + " psi4.set_options({'cubeprop_tasks' : ['density']})\n", + " elif self.type_cubefile[0] == \"potential\":\n", + " psi4.set_options({'cubeprop_tasks' : ['ESP']})\n", + " else:\n", + " raise Exception(\"Unrecognized cube file type.\")\n", + " return self.run_cubeprop()\n", + " \n", + " \n", + " def plot(self, isovalue, grid_spacing, type_cubefile, surface_opacity, mesh_opacity, show_atoms):\n", + " \n", + " # Recompute cube data if needed. This is quite expensive, so only update when necessary.\n", + " if not self.data or grid_spacing != self.spacing or type_cubefile != self.type_cubefile:\n", + " self.spacing = grid_spacing\n", + " self.type_cubefile = type_cubefile\n", + " self.data, self.espdata = self.process_cubefile()\n", + " \n", + " spacetuple = (self.data['incx'], self.data['incy'], self.data['incz'])\n", + " mclkwargs = {\"level\" : isovalue, \"spacing\" : spacetuple}\n", + " pospltkwargs = {\"color\" : (0, 0, 1, surface_opacity/100),\n", + " \"edgecolors\" : (0,0,0,mesh_opacity/100),\n", + " \"lw\" : 0.02}\n", + " negpltkwargs = {\"color\" : (1, 0, 0, surface_opacity/100),\n", + " \"edgecolors\" : (0,0,0,mesh_opacity/100),\n", + " \"lw\" : 0.02}\n", + " \n", + " self.ax.cla()\n", + " \n", + " mcl = measure.marching_cubes_lewiner\n", + " # the positive plot info; needed for all plot types\n", + " verts, faces, normals, values = mcl(self.data['data'], **mclkwargs)\n", + " \n", + " if type_cubefile[0] == \"potential\":\n", + " raise Exception(\"ESP plot doesn't work yet!\")\n", + " # Matplotlib unfortunately colors only using the z coordinate for a colormap\n", + " # so we can't pass in an array of colors (like espcolors, below) to project\n", + " # a function onto a surface. The other 3D plotting routines do allow this, so\n", + " # I'll leave this code in place in the hopes that the MPL folks add the trisurf\n", + " # custom face colors. If not, I'll add it myself if there's sufficient demand.\n", + " \n", + " # plot the positive lobe, computed above\n", + " ndata = verts[:].shape[0]\n", + " espcolors = np.zeros(ndata)\n", + " espdata = self.espdata['data']\n", + " for point in range(ndata):\n", + " xyz = verts[point]\n", + " xval = int(xyz[0] / self.data['incx'])\n", + " yval = int(xyz[1] / self.data['incy'])\n", + " zval = int(xyz[2] / self.data['incz'])\n", + " espcolors[point] = espdata[xval,yval,zval]\n", + " \n", + " self.ax.plot_trisurf(verts[:, 0], verts[:,1], faces, verts[:, 2], \n", + " cmap='bwr',\n", + " edgecolors=(0,0,0,mesh_opacity/100), lw=0.02)\n", + " elif type_cubefile[0] in [\"density\", \"orbital\"]:\n", + " # plot the positive lobe, computed above\n", + " self.ax.plot_trisurf(verts[:, 0], verts[:,1], faces, verts[:, 2], **pospltkwargs)\n", + " try:\n", + " # compute and plot the negative lobe\n", + " verts, faces, normals, values = mcl(-self.data['data'], **mclkwargs)\n", + " self.ax.plot_trisurf(verts[:, 0], verts[:,1], faces, verts[:, 2], **negpltkwargs)\n", + " except:\n", + " # Some orbitals/densities are positive all over, so skip\n", + " pass\n", + " else:\n", + " raise Exception(\"Unknow cube type in plot() function\")\n", + "\n", + " self.ax.set_xticks([])\n", + " self.ax.set_yticks([])\n", + " self.ax.set_zticks([])\n", + "\n", + " if show_atoms:\n", + " xs = self.data['atom_coords'][:,0]-self.data['minx']\n", + " ys = self.data['atom_coords'][:,1]-self.data['miny']\n", + " zs = self.data['atom_coords'][:,2]-self.data['minz']\n", + " self.ax.scatter(xs, ys, zs, c='g', s=15.0*self.data['atom_numbers'])\n", + "\n", + " plt.axis('off')\n", + " plt.show()\n", + " \n", + " def __init__(self, wfn):\n", + " self.data = None\n", + " self.espdata = None\n", + " self.wfn = wfn\n", + " self.spacing = 0.3\n", + " self.fig = plt.figure()\n", + " self.ax = self.fig.add_subplot(111, projection='3d')\n", + " self.type_cubefile = (\"density\", \"Dt.cube\")\n", + " \n", + " orbsyms = [ irrep for irrep in range(wfn.nirrep()) for sym in range(wfn.nmopi()[irrep]) ]\n", + " aevals = wfn.epsilon_a().nph\n", + " bevals = wfn.epsilon_b().nph\n", + " aorbenergies = [ aevals[irrep][orb] for irrep in range(wfn.nirrep()) for orb in range(wfn.nmopi()[irrep]) ]\n", + " borbenergies = [ bevals[irrep][orb] for irrep in range(wfn.nirrep()) for orb in range(wfn.nmopi()[irrep]) ]\n", + " aevals = sorted(zip(aorbenergies, orbsyms))\n", + " bevals = sorted(zip(borbenergies, orbsyms))\n", + " def process_evals(label, nocc, evals):\n", + " irrep_labels = wfn.molecule().irrep_labels()\n", + " symcount = [ 1 for irrep in range(wfn.nirrep()) ]\n", + " orbitaldict = {}\n", + " for orbital,evaltuple in enumerate(evals):\n", + " evalue = evaltuple[0]\n", + " irrep = evaltuple[1]\n", + " symlabel = irrep_labels[irrep]\n", + " filename = 'Psi_{}_{}_{}-{}.cube'.format(label[0], orbital+1, symcount[irrep], symlabel)\n", + " labelstr = '{} orbital {} ({}{}) {}{}, eigenvalue = {:.6}'\n", + " if orbital < nocc: \n", + " orblabel = labelstr.format(\n", + " label, orbital+1, symcount[irrep], symlabel, \"HOMO\", str(orbital-nocc), evalue)\n", + " elif orbital == nocc:\n", + " orblabel = labelstr.format(\n", + " label, orbital+1, symcount[irrep], symlabel, \"HOMO\", \"\", evalue)\n", + " elif orbital == nocc+1:\n", + " orblabel = labelstr.format(\n", + " label, orbital+1, symcount[irrep], symlabel, \"LUMO\", \"\", evalue)\n", + " else:\n", + " orblabel = labelstr.format(\n", + " label, orbital+1, symcount[irrep], symlabel, \"LUMO\", \"+\"+str(orbital-nocc-1), evalue)\n", + " symcount[irrep] += 1\n", + " orbitaldict[orblabel] = (\"orbital\", filename)\n", + " return orbitaldict\n", + " plotdict = {\"total density\" : (\"density\", \"Dt.cube\"),\n", + " \"alpha density\" : (\"density\", \"Da.cube\"),\n", + " \"beta density\" : (\"density\", \"Db.cube\"),\n", + " \"spin density\" : (\"density\", \"Ds.cube\"),\n", + " # disabled: see note in plot function \"potential\" : (\"potential\", \"ESP.cube\")\n", + " }\n", + " plotdict.update(process_evals('alpha', wfn.nalpha(), aevals))\n", + " plotdict.update(process_evals('beta', wfn.nbeta(), bevals))\n", + " interact(self.plot,\n", + " show_atoms=widgets.Checkbox(value=True),\n", + " isovalue=widgets.FloatSlider(min=0,max=0.1,step=.01,value=0.01,\n", + " description=\"Isovalue: \"),\n", + " mesh_opacity=widgets.FloatSlider(min=0.0,max=100.0,step=10.0,value=40.0,\n", + " description=\"Mesh Opacity: \"),\n", + " surface_opacity=widgets.FloatSlider(min=0.0,max=100.0,step=10.0,value=40.0,\n", + " description=\"Surface Opacity: \"),\n", + " grid_spacing=widgets.FloatSlider(min=0.05,max=0.5,step=.05,value=self.spacing,\n", + " description=\"Grid spacing: \"),\n", + " type_cubefile=widgets.Dropdown(options=plotdict, value=self.type_cubefile,\n", + " description=\"Plot type: \"))" + ] + }, + { + "cell_type": "code", + "execution_count": 86, + "metadata": {}, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " if (mpl.ratio != 1) {\n", + " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", + " }\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " this.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var backingStore = this.context.backingStorePixelRatio ||\n", + "\tthis.context.webkitBackingStorePixelRatio ||\n", + "\tthis.context.mozBackingStorePixelRatio ||\n", + "\tthis.context.msBackingStorePixelRatio ||\n", + "\tthis.context.oBackingStorePixelRatio ||\n", + "\tthis.context.backingStorePixelRatio || 1;\n", + "\n", + " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width * mpl.ratio);\n", + " canvas.attr('height', height * mpl.ratio);\n", + " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('