{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 02 Harmonic Vibrational Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The purpose of this project is to extend your fundamental Python language programming techniques through a normal coordinate/harmonic vibrational frequency calculation. The theoretical background and a concise set of instructions for this project may be found [here](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2302/project2-instructions.pdf).\n",
    "\n",
    "Original authors (Crawford, et al.) thank Dr. Yukio Yamaguchi of the University of Georgia for the original version of this project."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": [
     "thebe-init",
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "# Following os.chdir code is only for thebe (live code), since only in thebe default directory is /home/jovyan\n",
    "import os\n",
    "if os.getcwd().split(\"/\")[-1] != \"Project_02\":\n",
    "    os.chdir(\"source/Project_02\")\n",
    "from solution_02 import Molecule as SolMol\n",
    "# Following code is called for numpy array pretty printing\n",
    "import numpy as np\n",
    "np.set_printoptions(precision=7, linewidth=120, suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: Read the Coordinate Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coordinate data are given in a format identical to that for [Project #1](../Project_01/Project_01.ipynb). The test case for the remainder of this project is the water molecule ({download}`input/h2o_geom.txt`), optimized at the SCF/DZP level of theory. You can find the coordinates (in bohr) in the input directory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this project, `Molecule` object can be initialized as the following toggled code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": [
     "thebe-init",
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "class Molecule:\n",
    "\n",
    "    def __init__(self):\n",
    "        self.atom_charges = NotImplemented  # type: np.ndarray\n",
    "        self.atom_coords = NotImplemented  # type: np.ndarray\n",
    "        self.natm = NotImplemented  # type: int\n",
    "        self.hess = NotImplemented  # type: np.ndarray\n",
    "\n",
    "    def construct_from_dat_file(self, file_path: str):\n",
    "        # Same to Project 01\n",
    "        with open(file_path, \"r\") as f:\n",
    "            dat = np.array([line.split() for line in f.readlines()][1:])\n",
    "            self.atom_charges = np.array(dat[:, 0], dtype=float).astype(int)\n",
    "            self.atom_coords = np.array(dat[:, 1:4], dtype=float)\n",
    "            self.natm = self.atom_charges.shape[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Read the Cartesian Hessian Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The primary input data for the harmonic vibrational calculation is the Hessian matrix, which consists of second derivatives of the energy with respect to atomic positions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F_{r_i s_j} = \\frac{\\partial^2 E}{\\partial q_{r_i} \\partial q_{s_j}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notations:\n",
    "\n",
    "- $E$: Total energy of molecule\n",
    "- $i, j$: refer to index of atom (0, 1, 2, ...)\n",
    "- $r, s$: refer to coordinate components ($x$, $y$, $z$)\n",
    "- $q_{x_0}$: $x$ coordinate component of atom 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Hessian matrix (in units of $E_\\mathrm{h} / a_0^2$, where $E_\\mathrm{h}$ stands for Hartree energy, and $a_0$ Bohr radius) can be downloaded here ({download}`input/h2o_hessian.txt`) for the H2O test case. The first integer in the file is the number of atoms (which you may compare to the corresponding value from the geometry file as a test of consistency), while the remaining values have the following format:\n",
    "\n",
    "$$\n",
    "\\begin{matrix}\n",
    "F_{x_1, x_1} & F_{x_1, y_1} & F_{x_1, z_1} \\\\\n",
    "F_{x_1, x_2} & F_{x_1, y_2} & F_{x_1, z_2} \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "F_{x_2, x_1} & F_{x_2, y_1} & F_{x_2, z_1} \\\\\n",
    "\\vdots & \\vdots & \\vdots\n",
    "\\end{matrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "````{admonition} Hint 1: Reading and storing Hessian\n",
    ":class: dropdown\n",
    "\n",
    "The Hessian stored in memory should be a *symmetric* matrix, while the format of the input file is rectangular. Understanding the translation between the two takes a bit of thinking.\n",
    "\n",
    "One may take advantage of `numpy.reshape` ([NumPy API](https://numpy.org/doc/stable/reference/generated/numpy.reshape.html)) to compactly convert rectangular matrix to symmetric matrix.\n",
    "\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reader should fill all `NotImplementedError` in the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "def obtain_hessian(mole: Molecule, file_path: str):\n",
    "    # Input: Read Hessian file from `file_path`\n",
    "    # Attribute modification: Obtain raw Hessian to `mole.hess`\n",
    "    raise NotImplementedError(\"About 2~15 lines of code\")\n",
    "\n",
    "Molecule.obtain_hessian = obtain_hessian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.0927643,  0.       ,  0.       , -0.0463822,  0.       ,  0.       , -0.0463822,  0.       ,  0.       ],\n",
       "       [ 0.       ,  0.3171327,  0.       ,  0.       , -0.1585664,  0.0800202,  0.       , -0.1585664, -0.0800202],\n",
       "       [ 0.       ,  0.       ,  0.2800907,  0.       ,  0.0347766, -0.1400454,  0.       , -0.0347766, -0.1400454],\n",
       "       [-0.0463822,  0.       ,  0.       ,  0.0514668,  0.       ,  0.       , -0.0050847,  0.       ,  0.       ],\n",
       "       [ 0.       , -0.1585664,  0.0347766,  0.       ,  0.1730076, -0.0573984,  0.       , -0.0144412,  0.0226218],\n",
       "       [ 0.       ,  0.0800202, -0.1400454,  0.       , -0.0573984,  0.1268373,  0.       , -0.0226218,  0.013208 ],\n",
       "       [-0.0463822,  0.       ,  0.       , -0.0050847,  0.       ,  0.       ,  0.0514668,  0.       ,  0.       ],\n",
       "       [ 0.       , -0.1585664, -0.0347766,  0.       , -0.0144412, -0.0226218,  0.       ,  0.1730076,  0.0573984],\n",
       "       [ 0.       , -0.0800202, -0.1400454,  0.       ,  0.0226218,  0.013208 ,  0.       ,  0.0573984,  0.1268373]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol_mole = SolMol()\n",
    "sol_mole.construct_from_dat_file(\"input/h2o_geom.txt\")\n",
    "sol_mole.obtain_hessian(\"input/h2o_hessian.txt\")\n",
    "sol_mole.hess"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Mass-Weight the Hessian Matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Divide each element of the Hessian matrix by the product of square-roots of the masses of the atoms associated with the given coordinates:\n",
    "\n",
    "$$\n",
    "F_{r_i s_j}^\\mathrm{M} = \\frac{F_{r_i s_i}}{\\sqrt{m_i m_j}}\n",
    "$$\n",
    "\n",
    "where $m_i$ represents the mass of the atom corresponding to atom $i$. Use atomic mass units ($\\mathsf{amu}$) for the masses, just as for [Project 01](../Project_01/Project_01.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reader should fill all `NotImplementedError` in the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "def mass_weighted_hess(mole: Molecule) -> np.ndarray or list:\n",
    "    # Output: Mass-weighted Hessian matrix (in unit Eh/(amu*a0^2))\n",
    "    raise NotImplementedError(\"About 2~10 lines of code\")\n",
    "\n",
    "Molecule.mass_weighted_hess = mass_weighted_hess"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.005798 ,  0.       ,  0.       , -0.01155  ,  0.       ,  0.       , -0.01155  ,  0.       ,  0.       ],\n",
       "       [ 0.       ,  0.0198215,  0.       ,  0.       , -0.0394859,  0.0199265,  0.       , -0.0394859, -0.0199265],\n",
       "       [ 0.       ,  0.       ,  0.0175063,  0.       ,  0.00866  , -0.0348738,  0.       , -0.00866  , -0.0348738],\n",
       "       [-0.01155  ,  0.       ,  0.       ,  0.0510613,  0.       ,  0.       , -0.0050446,  0.       ,  0.       ],\n",
       "       [ 0.       , -0.0394859,  0.00866  ,  0.       ,  0.1716445, -0.0569462,  0.       , -0.0143274,  0.0224436],\n",
       "       [ 0.       ,  0.0199265, -0.0348738,  0.       , -0.0569462,  0.1258381,  0.       , -0.0224436,  0.013104 ],\n",
       "       [-0.01155  ,  0.       ,  0.       , -0.0050446,  0.       ,  0.       ,  0.0510613,  0.       ,  0.       ],\n",
       "       [ 0.       , -0.0394859, -0.00866  ,  0.       , -0.0143274, -0.0224436,  0.       ,  0.1716445,  0.0569462],\n",
       "       [ 0.       , -0.0199265, -0.0348738,  0.       ,  0.0224436,  0.013104 ,  0.       ,  0.0569462,  0.1258381]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol_mole.mass_weighted_hess()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: Diagonalize the Mass-Weighted Hessian Matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the eigenvalues of the mass-weighted Hessian:\n",
    "\n",
    "$$\n",
    "\\mathbf{F}^\\mathrm{M} \\mathbf{L} = \\mathbf{L} \\mathbf{\\Lambda}\n",
    "$$\n",
    "\n",
    "where $\\mathbf{\\Lambda}$ is a diagonal matrix containing eigenvalues, and $\\mathbf{L}$ contains eigenvectors. You should consider using the same canned diagonalization function you used in [Project 01](../Project_01/Project_01.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reader should fill all `NotImplementedError` in the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "def eig_mass_weight_hess(mole: Molecule) -> np.ndarray or list:\n",
    "    # Output: Eigenvalue of mass-weighted Hessian matrix (in unit Eh/(amu*a0^2))\n",
    "    raise NotImplementedError(\"Exactly 1 line of code using numpy\")\n",
    "\n",
    "Molecule.eig_mass_weight_hess = eig_mass_weight_hess"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.       ,  0.       ,  0.       ,  0.0518147,  0.0547485,  0.0561059,  0.1317344,  0.2106858,  0.2351242])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol_mole.eig_mass_weight_hess()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5: Compute the Harmonic Vibrational Frequencies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The vibrational frequencies are proportional to the squareroot of the eigenvalues of the mass-weighted Hessian:\n",
    "\n",
    "$$\n",
    "\\tilde \\omega_i = \\mathrm{constant} \\times \\sqrt{\\lambda_i}\n",
    "$$\n",
    "\n",
    "The most common units to use for vibrational frequencies is $\\mathsf{cm}^{-1}$ ([spectroscopy wavenumber](https://en.wikipedia.org/wiki/Wavenumber))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "````{admonition} Hint 1: Unit dimension analysis\n",
    ":class: dropdown\n",
    "\n",
    "It should be relatively appearant that unit of eigenvalues $\\sqrt{\\lambda_i}$ is the same to $F_{r_i s_j}^\\mathrm{M}$, i.e. $E_\\mathrm{h} / (\\mathsf{amu} \\cdot a_0^2)$.\n",
    "\n",
    "If we regard $E_\\mathrm{h}$, $\\mathrm{amu}$, $a_0$ as unit conversion constant, then we may have\n",
    "\n",
    "$$\n",
    "\\frac{E_\\mathrm{h}}{\\mathrm{amu} \\, a_0^2} \\frac{\\mathsf{J}}{\\mathsf{kg} \\cdot \\mathsf{m}^2} = \\frac{E_\\mathrm{h}}{\\mathrm{amu} \\, a_0^2} \\, \\mathsf{s}^{-2}\n",
    "$$\n",
    "\n",
    "So unit of $\\sqrt{\\lambda_i}$ is exactly $\\mathsf{s}^{-1} = \\mathsf{Hz}$. In spectroscopy, wavenumber refers to a frequency which is divided by speed of light in vacuum:\n",
    "\n",
    "$$\n",
    "\\tilde \\omega_i = \\frac{\\sqrt{\\lambda_i}}{2 \\pi c}\n",
    "$$\n",
    "\n",
    "So finally, the unit conversion should be\n",
    "\n",
    "$$\n",
    "\\frac{\\mathrm{centi}}{2 \\pi c} \\sqrt{\\frac{E_\\mathrm{h}}{\\mathrm{amu} \\, a_0^2}} \\, \\mathsf{cm}^{-1}\n",
    "$$\n",
    "\n",
    "All the unit conversion constants in the fomular above could be simply obtained from `scipy.constants` ([SciPy API](https://docs.scipy.org/doc/scipy/reference/constants.html)).\n",
    "\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "````{admonition} Hint 2: Imaginary frequency\n",
    ":class: dropdown\n",
    "\n",
    "In some cases, you may find some $\\lambda_i$ could be smaller than zero, making $\\sqrt{\\lambda_i}$ becoming an imaginary value. This \"frequency\" is called imaginary frequency. If these values are far from zero, they can be some indication that the molecule is far from optimized geometry or near transition state.\n",
    "\n",
    "We do not discuss imaginary frequency in detail here. Usually, these imaginary values is expressed in minus values in common quantum chemistry softwares (so the reader should bear in mind that minus frequencies are actually imaginary). So, the code is something like\n",
    "\n",
    "```python\n",
    "return np.sign(eigs) * np.sqrt(np.abs(eigs)) * unit_conversion\n",
    "```\n",
    "\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "````{admonition} Hint 3: \"Rotational\" Frequencies\n",
    ":class: dropdown\n",
    "\n",
    "For a molecule which is fully geometry optimized, there should be three lowest frequencies *exactly* to be zero. However, you may find there exist three frequencies *near to* zero. So the structure used in the computation is not exactly the stationary point on the potential energy surface (PES), although really close to that.\n",
    "\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reader should fill all `NotImplementedError` in the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "def harmonic_vib_freq(mole: Molecule) -> np.ndarray or list:\n",
    "    # Output: Harmonic vibrational frequencies (in unit cm^-1)\n",
    "    raise NotImplementedError(\"About 2~15 lines of code\")\n",
    "\n",
    "Molecule.harmonic_vib_freq = harmonic_vib_freq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  -0.0000358,    0.       ,    0.020977 , 1170.1214147, 1202.791828 , 1217.611377 , 1865.752138 , 2359.5107516,\n",
       "       2492.6021517])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol_mole.harmonic_vib_freq()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test Cases"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Water**\n",
    "- Geometry file: {download}`input/h2o_geom.txt`\n",
    "- Hessian file: {download}`input/h2o_hessian.txt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== Mass-Weighted Hessian Matrix (in unit Eh/(amu*a0^2)) ===\n",
      "[[ 0.005798   0.         0.        -0.01155    0.         0.        -0.01155    0.         0.       ]\n",
      " [ 0.         0.0198215  0.         0.        -0.0394859  0.0199265  0.        -0.0394859 -0.0199265]\n",
      " [ 0.         0.         0.0175063  0.         0.00866   -0.0348738  0.        -0.00866   -0.0348738]\n",
      " [-0.01155    0.         0.         0.0510613  0.         0.        -0.0050446  0.         0.       ]\n",
      " [ 0.        -0.0394859  0.00866    0.         0.1716445 -0.0569462  0.        -0.0143274  0.0224436]\n",
      " [ 0.         0.0199265 -0.0348738  0.        -0.0569462  0.1258381  0.        -0.0224436  0.013104 ]\n",
      " [-0.01155    0.         0.        -0.0050446  0.         0.         0.0510613  0.         0.       ]\n",
      " [ 0.        -0.0394859 -0.00866    0.        -0.0143274 -0.0224436  0.         0.1716445  0.0569462]\n",
      " [ 0.        -0.0199265 -0.0348738  0.         0.0224436  0.013104   0.         0.0569462  0.1258381]]\n",
      "=== Eigenvalue of Mass-Weighted Hessian Matrix (in unit Eh/(amu*a0^2)) ===\n",
      "[-0.         0.         0.         0.0518147  0.0547485  0.0561059  0.1317344  0.2106858  0.2351242]\n",
      "=== Harmonic Vibrational Frequencies (in unit cm^-1) ===\n",
      "[  -0.0000358    0.           0.020977  1170.1214147 1202.791828  1217.611377  1865.752138  2359.5107516 2492.6021517]\n"
     ]
    }
   ],
   "source": [
    "sol_mole = SolMol()\n",
    "sol_mole.construct_from_dat_file(\"input/h2o_geom.txt\")\n",
    "sol_mole.obtain_hessian(\"input/h2o_hessian.txt\")\n",
    "sol_mole.print_solution_02()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Benzene**\n",
    "- Geometry file: {download}`input/benzene_geom.txt`\n",
    "- Hessian file: {download}`input/benzene_hessian.txt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== Mass-Weighted Hessian Matrix (in unit Eh/(amu*a0^2)) ===\n",
      "[[ 0.0696916  0.         0.        ... -0.0234218  0.         0.       ]\n",
      " [ 0.         0.0855755  0.        ...  0.        -0.136488   0.       ]\n",
      " [ 0.         0.         0.0156863 ...  0.         0.        -0.0153281]\n",
      " ...\n",
      " [-0.0234218  0.         0.        ...  0.0763606  0.         0.       ]\n",
      " [ 0.        -0.136488   0.        ...  0.         0.4833777  0.       ]\n",
      " [ 0.         0.        -0.0153281 ...  0.         0.         0.0365666]]\n",
      "=== Eigenvalue of Mass-Weighted Hessian Matrix (in unit Eh/(amu*a0^2)) ===\n",
      "[-0.        -0.         0.         0.         0.         0.         0.0086435  0.0086435  0.0185369  0.0185369\n",
      "  0.0248907  0.0267604  0.0407475  0.0407475  0.0503805  0.0505931  0.051977   0.0536311  0.0536311  0.0558216\n",
      "  0.0568439  0.0568439  0.0711535  0.0711535  0.0717443  0.0962824  0.1188397  0.1188397  0.1412093  0.1412093\n",
      "  0.5192258  0.5244019  0.5244019  0.5282622  0.5282622  0.5313338]\n",
      "=== Harmonic Vibrational Frequencies (in unit cm^-1) ===\n",
      "[  -0.0164509   -0.0116326    0.0246764    0.6680535    0.6680535    0.6691684  477.9142727  477.9142727  699.8795264\n",
      "  699.8795269  811.0042569  840.9121042 1037.65959   1037.65959   1153.81349   1156.2446375 1171.9516278 1190.4538092\n",
      " 1190.4538092 1214.5225828 1225.5931503 1225.5931506 1371.2049457 1371.2049459 1376.8857923 1595.0623276 1772.087214\n",
      " 1772.0872146 1931.6832235 1931.6832236 3704.0974408 3722.514595  3722.5145955 3736.1908987 3736.1908987 3747.0369861]\n"
     ]
    }
   ],
   "source": [
    "sol_mole = SolMol()\n",
    "sol_mole.construct_from_dat_file(\"input/benzene_geom.txt\")\n",
    "sol_mole.obtain_hessian(\"input/benzene_hessian.txt\")\n",
    "sol_mole.print_solution_02()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3-chloro-1-butene**\n",
    "- Geometry file: {download}`input/3c1b_geom.txt`\n",
    "- Hessian file: {download}`input/3c1b_hessian.txt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== Mass-Weighted Hessian Matrix (in unit Eh/(amu*a0^2)) ===\n",
      "[[ 0.0720844 -0.0016129  0.0003104 ... -0.0009051 -0.0017714  0.0003552]\n",
      " [-0.0016129  0.0597496  0.0045083 ... -0.0003417 -0.0006569  0.0001635]\n",
      " [ 0.0003104  0.0045083  0.0707321 ... -0.0003336 -0.0005921  0.0002486]\n",
      " ...\n",
      " [-0.0009051 -0.0003417 -0.0003336 ...  0.0073687  0.0026435  0.0017749]\n",
      " [-0.0017714 -0.0006569 -0.0005921 ...  0.0026435  0.0022558  0.0007853]\n",
      " [ 0.0003552  0.0001635  0.0002486 ...  0.0017749  0.0007853  0.0013844]]\n",
      "=== Eigenvalue of Mass-Weighted Hessian Matrix (in unit Eh/(amu*a0^2)) ===\n",
      "[-0.000414   0.         0.         0.         0.         0.0000001  0.0000001  0.0022826  0.0027992  0.0046311\n",
      "  0.0065116  0.0088534  0.0213522  0.0283228  0.0409907  0.0498683  0.0511909  0.0569821  0.0595559  0.0640405\n",
      "  0.0690233  0.0835733  0.089757   0.0940067  0.110188   0.1124069  0.1240443  0.1246022  0.1593426  0.4844511\n",
      "  0.4991677  0.5062146  0.5254742  0.5350153  0.537149   0.5503098]\n",
      "=== Harmonic Vibrational Frequencies (in unit cm^-1) ===\n",
      "[-104.5906768    0.0099111    0.0289641    0.0304254    0.3665457    1.5129799    1.7406089  245.5941566  271.9691098\n",
      "  349.8231157  414.8104445  483.6819181  751.1480607  865.1127581 1040.7507835 1147.9332085 1163.0562807 1227.0817046\n",
      " 1254.488878  1300.863303  1350.5239405 1486.0657982 1540.0631407 1576.0995807 1706.3629779 1723.4581894 1810.4757465\n",
      " 1814.5421487 2051.9665893 3577.9088302 3631.8466432 3657.3927489 3726.3186775 3759.995898  3767.4862956 3813.3608088]\n"
     ]
    }
   ],
   "source": [
    "sol_mole = SolMol()\n",
    "sol_mole.construct_from_dat_file(\"input/3c1b_geom.txt\")\n",
    "sol_mole.obtain_hessian(\"input/3c1b_hessian.txt\")\n",
    "sol_mole.print_solution_02()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Wilson, E. B.; Decius, J. C.; Cross, P. C. *Molecular Vibrations* Dover Publication Inc., 1980.\n",
    "\n",
    "    ISBN-13: 978-0486639413"
   ]
  }
 ],
 "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.8.1"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}