{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Vibrational Modes of a Molecule\n\nLet's calculate the vibrational modes of :mol:`H_2O`.\n\nConsider the molecule at its equilibrium positions. If we displace the\natoms slightly, the energy $E$ will increase, and restoring\nforces will make the atoms oscillate in some pattern around the equilibrium.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can Taylor expand the energy with respect to the 9 coordinates\n(generally $3N$ coordinates for a molecule with\n$N$ atoms), $u_i$:\n\n\\begin{align}E = E_0 + \\frac{1}{2}\\sum_{i}^{3N} \\sum_{j}^{3N}\n \\frac{\\partial^2 E}{\\partial u_{i}\\partial u_{j}}\\bigg\\rvert_0\n (u_i - u_{i0}) (u_j - u_{j0}) + \\cdots\\end{align}\n\nSince we are expanding around the equilibrium positions, the energy\nshould be stationary and we can omit linear contributions.\n\nThe matrix of all the second derivatives is called the Hessian,\n$\\mathbf H$, and it expresses a linear system of differential equations\n\n\\begin{align}\\end{align}\n\n \\mathbf{Hu}_k = \\omega_k^2\\mathbf{Mu}_k\n\nfor the vibrational eigenmodes $u_k$ and their frequencies\n$\\omega_k$ that will characterise the collective movement of the atoms.\nIn short, we need the eigenvalues and eigenvectors of the Hessian.\n\nThe elements of the Hessian can be approximated as\n\n\\begin{align}\\end{align}\n\n H_{ij} = \\frac{\\partial^2 E}{\\partial u_{i}\\partial u_{j}}\\bigg\\rvert_0\n = -\\frac{\\partial F_{j}}{\\partial u_{i}},\n\nwhere $F_j$ are the forces. Hence we calculate the derivative\nof the forces using finite differences.\nWe need to displace each atom back and forth along each Cartesian direction,\ncalculating forces at each configuration to establish\n$H_{ij} \\approx \\Delta F_{j} / \\Delta u_{i}$,\nthen get eigenvalues and vectors of that.\n\n\nASE provides the :class:`~ase.vibrations.Vibrations` class for this\npurpose. Note how the linked documentation contains an example for\nthe :mol:`N_2` molecule, which means we almost don't have to do any\nwork ourselves. We just scavenge the online ASE\ndocumentation like we always do, then hack as necessary until the thing runs.\n\n\n.. admonition:: Exercise\n\n Calculate the vibrational frequencies of :mol:`H_2O` using GPAW in\n LCAO mode, saving the modes to trajectory files. What are the\n frequencies, and what do the eigenmodes look like?\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution\nFirst, we import the necessary examples and build a water structure.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from math import cos, pi, sin\n\nfrom gpaw import GPAW\n\nfrom ase import Atoms\nfrom ase.build import molecule\nfrom ase.optimize import QuasiNewton\nfrom ase.vibrations import Vibrations\n\n# Water molecule:\nh2o = molecule('H2O', vacuum=3.5)\nd = 0.9575\nt = pi / 180 * 104.51\n\nh2o = Atoms(\n 'H2O', positions=[(0, 0, 0), (d, 0, 0), (d * cos(t), d * sin(t), 0)]\n)\n\nh2o.center(vacuum=3.5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we attach a calculator (here: GPAW) and compute the vibrational\nmodes of the water molecule.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "h2o.calc = GPAW(txt='h2o.txt', mode='lcao', basis='dzp', symmetry='off')\n\nQuasiNewton(h2o).run(fmax=0.05)\n\n\n\"\"\"Calculate the vibrational modes of a H2O molecule.\"\"\"\n\n# Create vibration calculator\nvib = Vibrations(h2o)\nvib.run()\nvib.summary(method='frederiksen')\n\n# Make trajectory files to visualize normal modes:\nfor mode in range(9):\n vib.write_mode(mode)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since there are nine coordinates, we get nine eigenvalues and\ncorresponding modes. However the three translational and three\nrotational degrees of freedom will contribute six \"modes\" that do not\ncorrespond to true vibrations. In principle there are no restoring\nforces if we translate or rotate the molecule, but these will\nnevertheless have different energies (often imaginary) because of\nvarious artifacts of the simulation such as the grid used to represent\nthe density, or effects of the simulation box size.\n\nLet's visualize the last vibrational mode of water, which corresponds\nto assymetric stretching of the H-O bonds.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.animation as animation\nimport matplotlib.pyplot as plt\n\nfrom ase.io import read\nfrom ase.visualize.plot import plot_atoms\n\n\nconfigs = read('vib.8.traj', ':')\nfig, ax = plt.subplots()\nplot_atoms(h2o, ax, rotation=('0x,0y,0z'))\nax.set_axis_off()\n\n\ndef animate(i):\n # Remove the previous atomic plot\n [p.remove() for p in ax.patches]\n plot_atoms(configs[i], ax, rotation=('0x,0y,0z'), show_unit_cell=1)\n ax.set_xlim(-4, 7)\n ax.set_ylim(-4, 7)\n ax.set_axis_off()\n return (ax,)\n\n\nani = animation.FuncAnimation(\n fig, animate, repeat=True, frames=len(configs) - 1, interval=200\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This solution is based on an example from the GPAW web page,\nwhere also other comments to this exercise can be found:\n\n[GPAW website](https://gpaw.readthedocs.io/tutorialsexercises/vibrational/vibrations/vibrations.html)_\n\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.13.7" } }, "nbformat": 4, "nbformat_minor": 0 }