{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 01 Molecular Geometry Analysis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The purpose of this project is to introduce you to fundamental Python programming techniques in the context of a scientific problem, viz. the calculation of the internal coordinates (bond lengths, bond angles, dihedral angles), moments of inertia, and rotational constants of a polyatomic molecule. A concise set of instructions for this project may be found [here](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2301/project1-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": "markdown", "metadata": { "tags": [ "hide-cell" ] }, "source": [ "Solution of this project could be {download}`solution_01.py`." ] }, { "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_01\":\n", " os.chdir(\"source/Project_01\")\n", "from solution_01 import Molecule as SolMol" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 1: Read the Coordinate Data from Input" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The input to the program is the set of Cartesian coordinates of the atoms (**in bohr**) and their associated atomic numbers. A sample molecule (acetaldehyde) to use as input to the program is:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " 7\n", " 6 0.000000000000 0.000000000000 0.000000000000\n", " 6 0.000000000000 0.000000000000 2.845112131228\n", " 8 1.899115961744 0.000000000000 4.139062527233\n", " 1 -1.894048308506 0.000000000000 3.747688672216\n", " 1 1.942500819960 0.000000000000 -0.701145981971\n", " 1 -1.007295466862 -1.669971842687 -0.705916966833\n", " 1 -1.007295466862 1.669971842687 -0.705916966833" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first line above is the number of atoms (an integer), while the remaining lines contain the z-values (atom charges) and x-, y-, and z-coordinates of each atom (one integer followed by three double-precision floating-point numbers). This input file({download}`input/acetaldehyde.dat`) along with a few other test cases can be found in this repository in the input directory." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After downloading the file to your computer (to a file called `geom.dat`, for example), you must open the file, read the data from each line into appropriate variables, and finally close the file." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": { "tags": [ "hide_cell", "thebe-init" ] }, "source": [ "````{admonition} Hint 1: Opening and closing the file\n", ":class: dropdown\n", "\n", "In Python, using `open` and `close` can easily handle file I/Os:\n", "\n", "```python\n", ">>> f = open(\"input/acetaldehyde.dat\", \"r\")\n", ">>> # handle file `f`\n", ">>> f.close()\n", "```\n", "\n", "However, there's more elegent way to do this by `with`:\n", "\n", "```python\n", ">>> with open(\"input/acetaldehyde.dat\", \"r\") as f:\n", ">>> pass # handle file `f`\n", "```\n", "\n", "The latter code is prefered. Using `with` would automatically close file, avoiding possibilities that programmer forget inserting a line to close file.\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 2: Reading the numbers\n", ":class: dropdown\n", "\n", "Handling multiple words (decimals, floats as strings) is relatively easy in Python. For example,\n", "\n", "```python\n", ">>> \"6 0.000000000000 0.000000000000 0.000000000000\".split()\n", "['6', '0.000000000000', '0.000000000000', '0.000000000000']\n", "```\n", "\n", "However, these are string types. One need to transform to proper types to enable calculation afterwards.\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 3: Storing the z-values (atom charges) and the coordinates by NumPy\n", ":class: dropdown\n", "\n", "There are two ways in python.\n", "\n", "An intutive way to do this is generating a list to store z-values (as array), and list of list to store coordinates (as matrix). Notice do not use string type, for the sake of convenience of calculation tasks afterwards. It is also the most used way to store array/matrix in C++ and fortran90+. In FORTRAN77, one must use array to store matrix, along with it's leading dimension.\n", "\n", "Another way is to use numpy. For example,\n", "\n", "```python\n", ">>> arr = np.array([\n", "... [\"0.0\", \"0.1\"],\n", "... [\"2.5\", \"4.6\"],\n", "... ], dtype=float)\n", ">>> arr\n", "array([[0. , 0.1],\n", " [2.5, 4.6]])\n", "```\n", "\n", "Numpy can be extremely powerful, not only it wraps many array/matrix properties (raw data, dimension, type, etc.) into one class (`numpy.ndarray`), but also offers powerful operations in very concise way (matrix multiplication, SVD, tensor contraction, etc.).\n", "\n", "Learning some essentials of NumPy can prove to be useful. In this document, NumPy will be the default array/matrix/tensor storage structure. Other robust array/matrix/tensor engines could be PyTorch/TensorFlow, etc.\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 4: Wrapping molecule information into class\n", ":class: dropdown\n", "\n", "What defines a molecule is somehow complicated. Not only it's atom charges, but also its coordinates. We can derive bond lengths, angles, etc. only from molecule information. We may term atom charges and coordinates as *Attribiutes*; they are variables and information describing molecule, like adjectives in gramma. We may term how we generate bond lengths, angles, etc. as *Methods*; they are functions, mapping from molecule information to some properties we desire, like verbs in gramma.\n", "\n", "An example to construct molecule class may be\n", "\n", "```python\n", "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", "\n", " def construct_from_dat_file(self, file_path: str):\n", " # Input: Read file from `file_path`\n", " # Attribute modification: Obtain atom charges to `self.atom_charges`\n", " # Attribute modification: Obtain coordinates to `self.atom_coords`\n", " # Attribute modification: Obtain atom numbers in molecule to `self.natm`\n", " raise NotImplementedError(\"Reader need to fill this method function!\")\n", "```\n", "\n", "For convenience, redundant variable `self.natm` (number of atoms in molecule) is also included as *Attribute*. It can be easily obtained by counting length of atom charges array `self.atom_charges`.\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implementation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- If `mole` is an instance of `Molecule` class, and it's acetaldehyde molecule, then\n", "\n", " - `mole.atom_charges` should be a NumPy array with shape (7, )\n", " \n", " - `mole.atom_coords` should be a NumPy array with shape (7, 3) (we use *array* instead of matrix, just convention; like any array/matrix in TensorFlow/PyTorch is termed *tensor*)\n", " \n", " - `mole.natm` should be an integer" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Reader should fill all `NotImplementedError` in the following code (notice that `NotImplemented` is coded intentionally, not the same to `NotImplementedError`):" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "import numpy as np\n", "\n", "\n", "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", "\n", " def construct_from_dat_file(self, file_path: str):\n", " # Input: Read file from `file_path`\n", " # Attribute modification: Obtain atom charges to `self.atom_charges`\n", " # Attribute modification: Obtain coordinates to `self.atom_coords`\n", " # Attribute modification: Obtain atom numbers in molecule to `self.natm`\n", " raise NotImplementedError(\"About 5~15 lines of code\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Following is reference output. Reader may repeat the output results with his/hers own implementation." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[6 6 8 1 1 1 1]\n", "[[ 0. 0. 0. ]\n", " [ 0. 0. 2.84511213]\n", " [ 1.89911596 0. 4.13906253]\n", " [-1.89404831 0. 3.74768867]\n", " [ 1.94250082 0. -0.70114598]\n", " [-1.00729547 -1.66997184 -0.70591697]\n", " [-1.00729547 1.66997184 -0.70591697]]\n", "7\n" ] } ], "source": [ "sol_mole = SolMol()\n", "sol_mole.construct_from_dat_file(\"input/acetaldehyde.dat\")\n", "print(sol_mole.atom_charges)\n", "print(sol_mole.atom_coords)\n", "print(sol_mole.natm)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 2: Bond Lengths" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate the interatomic distances using the expression:\n", "\n", "$$\n", "R_{ij} = \\sqrt{(x_i - x_j)^2 + (y_i - y_j)^2 + (z_i - z_j)^2}\n", "$$\n", "\n", "where $x$, $y$, and $z$ are Cartesian coordinates and $i$ and $j$ denote atomic indices.\n", "\n", "Print all bond lengths $R_{ij}$ in Bohr." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 1: Bond length method function signiture\n", ":class: dropdown\n", "\n", "Before writing program, we need to ask ourselfs, that what information do we need to definitely define bond length. Actually, we need the two atoms ($i, j$ as subscripts of $R_{ij}$), and molecule coordinates ($x_i, y_i, z_i$). So we define the program as\n", "\n", "```python\n", "def bond_length(mole: Molecule, i: int, j: int) -> float:\n", " # Input: `i`, `j` index of molecule's atom\n", " # Output: Bond length from atom `i` to atom `j`\n", " raise NotImplementedError(\"Reader need to fill this method function!\")\n", "```\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 2: $L_2$ norm\n", ":class: dropdown\n", "\n", "Directly programming by the expression of $R_{ij}$ is very intutive. However, this expression can also be written as\n", "\n", "$$\n", "R_{ij} = |\\boldsymbol{r}_i - \\boldsymbol{r}_j | := \\Vert \\boldsymbol{r}_i - \\boldsymbol{r}_j \\Vert_2\n", "$$\n", "\n", "where $\\boldsymbol{r}_i = (x_i, y_i, z_i)$. $R_{ij}$ could be also termed as $L_2$ (Frobenius) norm of vector $\\boldsymbol{r}_i - \\boldsymbol{r}_j$, and $L_2$ norm of vector is implemented by NumPy (`numpy.linalg.norm`). See [NumPy API](https://numpy.org/doc/stable/reference/generated/numpy.linalg.norm.html) for its usage.\n", "\n", "Using `numpy.linalg.norm` is preferred, but not compulsory at this stage.\n", "\n", "If reader is not ready to use NumPy, python's standard library `math` is prefered to make sqrt calculation. For example, norm of vector $(0, 3, 4)$ could be calculated as\n", "\n", "```python\n", ">>> import math\n", ">>> math.sqrt(0**2 + 3**2 + 4**2)\n", "5.0\n", "```\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 3: Pretty printing\n", ":class: dropdown\n", "\n", "We may use python's default printing utilities for debugging:\n", "\n", "```python\n", ">>> print(0, 1, 2.845112131228)\n", "0 1 2.845112131228\n", "```\n", "\n", "However, if we want to format printing string, we may use `str.format` method function. Reader may search `python string format` in search engines for more details. For example, if we want to print atom indexes in 3-char-long decimals and bond length in 10-char-long-5-after-decimal-point, we may write program as \n", "\n", "```python\n", ">>> print(\"{:3d} - {:3d}: {:10.5f} Bohr\".format(0, 1, 2.845112131228))\n", " 0 - 1: 2.84511 Bohr\n", "```\n", "\n", "You may wish to implement another function `print_bond_length` to realize this feature.\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 4: Loop structure\n", ":class: dropdown\n", "\n", "We need to print all bond length information in function `print_bond_length`, as Step 2 Hint 3 stated. Note that bond between atom 1 - atom 2 is the same to atom 2 - atom 1. So, reader should notice the range for each loop. In the following code, atom index `j` is always larger than `i`:\n", "\n", "```python\n", "for i in range(mole.natm):\n", " for j in range(i + 1, mole.natm):\n", " raise NotImplementedError # Reader should complete printing code\n", "```\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 5: Extending molecule class\n", ":class: dropdown\n", "\n", "We have discussed how to construct `Molecule` class in Step 1 Hint 4. However, we could regard `bond_length` function defined in Step 2 Hint 1 as method function of `Molecule` class. Though this is somehow not \"pythonic\", one can add method function (or attribute) after class construction, in fashion like C++. But we note that python is more flexable than C++ when handling methods in class, since in C++, although method implementation is not required when constructing class, one still need to declare signiture of function. In python, there's no need declaring function, so we just do implementations.\n", "\n", "In this case, if we have defined function `bond_length` in Step 2 Hint 1, we just make the following statement to add this function to `Molecule` class:\n", "\n", "```python\n", "Molecule.bond_length = bond_length\n", "```\n", "\n", "We note that when calling `bond_length`, we need use\n", "\n", "```python\n", "bond_length(mole, i, j)\n", "```\n", "\n", "but when calling `Molecule.bond_length`, we need use\n", "\n", "```python\n", "mole.bond_length(i, j)\n", "```\n", "\n", "The first parameter `mole: Molecule` in `Molecule.bond_length`, like other (non-static) method functions, is equivalent to `self`.\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": 4, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def bond_length(mole: Molecule, i: int, j: int) -> float:\n", " # Input: `i`, `j` index of molecule's atom\n", " # Output: Bond length from atom `i` to atom `j`\n", " raise NotImplementedError(\"About 1~3 lines of code\")\n", "\n", "Molecule.bond_length = bond_length" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def print_bond_length(mole: Molecule):\n", " # Usage: Print all bond length\n", " print(\"=== Bond Length ===\")\n", " for i in range(mole.natm):\n", " for j in range(i + 1, mole.natm):\n", " raise NotImplementedError(\"Exactly 1 line of code\")\n", "\n", "Molecule.print_bond_length = print_bond_length" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Bond Length ===\n", " 0 - 1: 2.84511 Bohr\n", " 0 - 2: 4.55395 Bohr\n", " 0 - 3: 4.19912 Bohr\n", " 0 - 4: 2.06517 Bohr\n", " 0 - 5: 2.07407 Bohr\n", " 0 - 6: 2.07407 Bohr\n", " 1 - 2: 2.29803 Bohr\n", " 1 - 3: 2.09811 Bohr\n", " 1 - 4: 4.04342 Bohr\n", " 1 - 5: 4.05133 Bohr\n", " 1 - 6: 4.05133 Bohr\n", " 2 - 3: 3.81330 Bohr\n", " 2 - 4: 4.84040 Bohr\n", " 2 - 5: 5.89151 Bohr\n", " 2 - 6: 5.89151 Bohr\n", " 3 - 4: 5.87463 Bohr\n", " 3 - 5: 4.83836 Bohr\n", " 3 - 6: 4.83836 Bohr\n", " 4 - 5: 3.38971 Bohr\n", " 4 - 6: 3.38971 Bohr\n", " 5 - 6: 3.33994 Bohr\n" ] } ], "source": [ "sol_mole.print_bond_length()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 3: Bond Angles" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate bond angles. For example, the angle, $\\phi_{ijk}$, between atoms $i$-$j$-$k$, where $j$ is the central atom is given by:\n", "\n", "$$\n", "\\cos \\phi_{ijk} = \\boldsymbol{e}_{ji} \\cdot \\boldsymbol{e}_{jk}\n", "$$\n", "\n", "where the $\\boldsymbol{e}_{ij}$ are unit vectors between the atoms, e.g.,\n", "\n", "$$\n", "e_{ij}^x = - \\frac{x_i - x_j}{R_{ij}}, \\quad e_{ij}^y = - \\frac{y_i - y_j}{R_{ij}}, \\quad e_{ij}^z = - \\frac{z_i - z_j}{R_{ij}}\n", "$$\n", "\n", "Print all valid bond angles $\\phi_{ijk}$ (with bond length $R_{ij} < 3 \\, \\mathsf{Bohr}$ and $R_{jk} < 3 \\, \\mathsf{Bohr}$, and $\\phi_{ijk} > 90^\\circ$) in degree.\n", "\n", "Original C++ project suggests storing unit vectors in memory. However, this python project will generate unit vectors on-the-fly. Also, the final result of this step is a little different to the original project." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 1: Unit vector generation\n", ":class: dropdown\n", "\n", "From the definition above, we may state that\n", "\n", "$$\n", "\\boldsymbol{e}_{ij} = \\frac{\\boldsymbol{r}_j - \\boldsymbol{r}_i}{| \\boldsymbol{r}_j - \\boldsymbol{r}_i |}\n", "$$\n", "\n", "It is a combined opeation of vector substraction, $L_2$ norm, and vector divided by numerical value. In NumPy, they are quite easy to achieve.\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 2: Vector inner product\n", ":class: dropdown\n", "\n", "We used inner product between $\\boldsymbol{e}_{ji}$ and $\\boldsymbol{e}_{jk}$ in the fomular above:\n", "\n", "$$\n", "\\cos \\phi_{ijk} = \\boldsymbol{e}_{ji} \\cdot \\boldsymbol{e}_{jk} = e_{ji}^x e_{jk}^x + e_{ji}^y e_{jk}^y + e_{ji}^z e_{jk}^z\n", "$$\n", "\n", "In NumPy, there are multiple ways to achieve vector inner product. An intutive way is using `numpy.inner` (see [NumPy API](https://numpy.org/doc/stable/reference/generated/numpy.inner.html)). For example,\n", "\n", "```python\n", ">>> np.inner([-1, 1, 2], [3, 9, 5])\n", "16\n", "```\n", "\n", "Other ways could be `numpy.dot`, `numpy.tensordot`, `numpy.einsum`, `np.multiply` with `np.sum`.\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 3: Inverse trigonometric functions\n", ":class: dropdown\n", "\n", "We need calculate arccos value of vector inner product. In NumPy, `numpy.arccos` can do this job. Without NumPy, there's still `math.acos`.\n", "\n", "Note that we need to print angle in degree, so we need to multiply result by $180 / \\pi$.\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 4: Loop structure\n", ":class: dropdown\n", "\n", "From the description, we should know that bond angle 1-2-3 is the same to 3-2-1. However, we regard 1-2-3 as different bond angle to 2-1-3. If bond length of 1-2 is larger than 4 Angstrom, then we silent all bond angles 1-2-\\*; but if bond length of 1-3 is larger than 3 Angstrom, as well as 1-2, 2-3 shorter than 3 Angstrom, we still print bond angle 1-2-3.\n", "\n", "Utilizing that the sum of interior angles of a triangle is 180°, we know that for any three atoms, they have utmost one angle larger than 90°.\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": 7, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def bond_unit_vector(mole: Molecule, i: int, j: int) -> np.ndarray:\n", " # Input: `i`, `j` index of molecule's atom\n", " # Output: Unit vector of bond from atom `i` to atom `j`\n", " raise NotImplementedError(\"About 1~5 lines of code\")\n", "\n", "Molecule.bond_unit_vector = bond_unit_vector" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def bond_angle(mole: Molecule, i: int, j: int, k: int) -> float:\n", " # Input: `i`, `j`, `k` index of molecule's atom; where `j` is the central atom\n", " # Output: Bond angle for atoms `i`-`j`-`k`\n", " raise NotImplementedError(\"About 3~5 lines of code\")\n", "\n", "Molecule.bond_angle = bond_angle" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def is_valid_angle(mole: Molecule, i: int, j: int, k: int) -> bool:\n", " # Input: `i`, `j`, `k` index of molecule's atom; where `j` is the central atom\n", " # Output: Test if `i`-`j`-`k` is a valid bond angle\n", " # if i != j != k\n", " # and if i-j and j-k bond length smaller than 3 Angstrom,\n", " # and if angle i-j-k > 90 degree\n", " raise NotImplementedError(\"About 1~5 lines of code\")\n", "\n", "Molecule.is_valid_angle = is_valid_angle" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def print_bond_angle(mole: Molecule):\n", " # Usage: Print all bond angle i-j-k which is considered as valid angle\n", " raise NotImplementedError(\"About 5~20 lines of code\")\n", "\n", "Molecule.print_bond_angle = print_bond_angle" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For this solution, there could be different outputs, such as angle `0-1-2` should be equal to `2-1-0`. Such kind of minor differences could be tolerated. However, angle values, number of valid angles, etc. should be exactly the same." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Bond Angle ===\n", " 0 - 1 - 2: 124.26831 Degree\n", " 0 - 1 - 3: 115.47934 Degree\n", " 1 - 0 - 4: 109.84706 Degree\n", " 1 - 0 - 5: 109.89841 Degree\n", " 1 - 0 - 6: 109.89841 Degree\n", " 4 - 0 - 5: 109.95368 Degree\n", " 4 - 0 - 6: 109.95368 Degree\n", " 5 - 0 - 6: 107.25265 Degree\n", " 2 - 1 - 3: 120.25235 Degree\n" ] } ], "source": [ "sol_mole.print_bond_angle()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 4: Out-of-Plane Angles" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate out-of-plane angles. For example, the angle $\\theta_{ijkl}$ for atom $i$ out of the plane containing atoms $j$-$k$-$l$ (with $k$ as the central atom, connected to $i$) is given by:\n", "\n", "$$\n", "\\sin \\theta_{ijkl} = \\frac{\\boldsymbol{e}_{kj} \\times \\boldsymbol{e}_{kl}}{\\sin \\phi_{jkl}} \\cdot \\boldsymbol{e}_{ki}\n", "$$\n", "\n", "Print all valid out-of-plane angles (bond angle $\\phi_{jkl}$ should be valid bond angle, and bond length $R_{ki} < 3 \\, \\mathsf{Bohr}$) in degree. Note that situation of $\\sin \\phi_{jkl} = 0$ could occasionally (especially Molecule HCN or H2C=C=CH2), so excluding this kind of situation is encouraged." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 1: Cross products\n", ":class: dropdown\n", "\n", "Cross product of two arrays (length of array is 3) $\\boldsymbol{c} = \\boldsymbol{a} \\times \\boldsymbol{b}$ is defined as\n", "\n", "$$\n", "\\begin{align}\n", "c_x &= a_y b_z - a_z b_y \\\\\n", "c_y &= a_z b_x - a_x b_z \\\\\n", "c_z &= a_x b_y - a_y b_x\n", "\\end{align}\n", "$$\n", "\n", "In NumPy, `numpy.cross` ([NumPy API](https://numpy.org/doc/stable/reference/generated/numpy.cross.html)) could fulfill this feature.\n", "\n", "```python\n", ">>> np.cross([5, 2, 3], [7, 8, 9])\n", "array([ -6, -24, 26])\n", "```\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 2: Numerical precision\n", ":class: dropdown\n", "\n", "Actually, when implementing arcsin operation, one may encounter `RuntimeWarning` from NumPy:\n", "\n", "```python\n", ">>> np.arcsin(-1.000000000000001)\n", ":1: RuntimeWarning: invalid value encountered in arcsin\n", "nan\n", "```\n", "\n", "To avoid such kind of error, one need to make sure the value passed into `numpy.arcsin` is not larger than 1 nor smaller than -1, strictly. One may also assert that the passed value is not too large or too small, in case of internal programming mistakes.\n", "\n", "```python\n", ">>> theta = -1.000000000000001\n", ">>> assert(np.abs(theta) < 1 + 1e-7) # make sure that theta is not too wrong, in case of programming mistakes\n", ">>> theta = np.sign(theta) if np.abs(theta) > 1 else theta # make theta bounded between -1 to 1\n", ">>> np.arcsin(theta)\n", "-1.5707963267948966\n", "```\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": 12, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def out_of_plane_angle(mole: Molecule, i: int, j: int, k: int, l: int) -> float:\n", " # Input: `i`, `j`, `k`, `l` index of molecule's atom; where `k` is the central atom, and angle is i - j-k-l\n", " # Output: Out-of-plane bond angle for atoms `i`-`j`-`k`-`l`\n", " raise NotImplementedError(\"About 3~8 lines of code\")\n", "\n", "Molecule.out_of_plane_angle = out_of_plane_angle" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def is_valid_out_of_plane_angle(mole: Molecule, i: int, j: int, k: int, l: int) -> bool:\n", " # Input: `i`, `j`, `k`, `l` index of molecule's atom; where `k` is the central atom, and angle is i - j-k-l\n", " # Output: Test if `i`-`j`-`k` is a valid out-of-plane bond angle\n", " # if i != j != k != l\n", " # and if angle j-k-l is valid bond angle\n", " # and if i-k bond length smaller than 3 Angstrom\n", " # and bond angle of j-k-l is not linear\n", " raise NotImplementedError(\"About 1~5 lines of code\")\n", "\n", "Molecule.is_valid_out_of_plane_angle = is_valid_out_of_plane_angle" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def print_out_of_plane_angle(mole: Molecule):\n", " # Usage: Print all out-of-plane bond angle i-j-k-l which is considered as valid\n", " raise NotImplementedError(\"About 6~15 lines of code\")\n", "\n", "Molecule.print_out_of_plane_angle = print_out_of_plane_angle" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Difference of output may also be tolerated in this solution." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Out-of-Plane Angle ===\n", " 3 - 0 - 1 - 2: 0.00000 Degree\n", " 2 - 0 - 1 - 3: 0.00000 Degree\n", " 5 - 1 - 0 - 4: -53.62632 Degree\n", " 6 - 1 - 0 - 4: 53.62632 Degree\n", " 4 - 1 - 0 - 5: 53.65153 Degree\n", " 6 - 1 - 0 - 5: -56.27711 Degree\n", " 4 - 1 - 0 - 6: -53.65153 Degree\n", " 5 - 1 - 0 - 6: 56.27711 Degree\n", " 1 - 4 - 0 - 5: -53.67878 Degree\n", " 6 - 4 - 0 - 5: 56.19462 Degree\n", " 1 - 4 - 0 - 6: 53.67878 Degree\n", " 5 - 4 - 0 - 6: -56.19462 Degree\n", " 1 - 5 - 0 - 6: -54.97706 Degree\n", " 4 - 5 - 0 - 6: 54.86999 Degree\n", " 0 - 2 - 1 - 3: 0.00000 Degree\n" ] } ], "source": [ "sol_mole.print_out_of_plane_angle()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 5: Torsion/Dihedral Angles" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate torsional angles. For example, the torsional angle $\\tau_{ijkl}$ for the atom connectivity $i$-$j$-$k$-$l$ is given by:\n", "\n", "$$\n", "\\cos \\tau_{ijkl} = \\frac{(\\boldsymbol{e}_{ji} \\times \\boldsymbol{e}_{jk}) \\cdot (\\boldsymbol{e}_{kj} \\times \\boldsymbol{e}_{kl})}{\\sin \\phi_{ijk} \\cdot \\sin \\phi_{jkl}}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Print all valid dihedral angles ($i, j, k$ and $j, k, l$ constructs a valid angle, where $j, k$ should be bonded or $R_{jk} < 3 \\, \\mathsf{Bohr}$) in degree." ] }, { "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": 16, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def dihedral_angle(mole: Molecule, i: int, j: int, k: int, l: int) -> float:\n", " # Input: `i`, `j`, `k`, `l` index of molecule's atom; where `k` is the central atom, and angle is i - j-k-l\n", " # Output: Dihedral angle for atoms `i`-`j`-`k`-`l`\n", " raise NotImplementedError(\"About 3~8 lines of code\")\n", "\n", "Molecule.dihedral_angle = dihedral_angle" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def is_valid_dihedral_angle(mole: Molecule, i: int, j: int, k: int, l: int) -> bool:\n", " # Input: `i`, `j`, `k`, `l` index of molecule's atom; where `k` is the central atom, and angle is i - j-k-l\n", " # Output: Test if `i`-`j`-`k` is a valid dihedral bond angle\n", " # if i != j != k != l\n", " # and if i, j, k construct a valid bond angle (with j-k bonded)\n", " # and if j, k, l construct a valid bond angle (with j-k bonded)\n", " raise NotImplementedError(\"About 1~5 lines of code\")\n", "\n", "Molecule.is_valid_dihedral_angle = is_valid_dihedral_angle" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def print_dihedral_angle(mole: Molecule):\n", " # Usage: Print all dihedral bond angle i-j-k-l which is considered as valid\n", " raise NotImplementedError(\"About 6~15 lines of code\")\n", "\n", "Molecule.print_dihedral_angle = print_dihedral_angle" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Difference of output may also be tolerated in this solution." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Dihedral Angle ===\n", " 2 - 0 - 1 - 3: 180.00000 Degree\n", " 2 - 0 - 1 - 4: 0.00000 Degree\n", " 2 - 0 - 1 - 5: 121.09759 Degree\n", " 2 - 0 - 1 - 6: 121.09759 Degree\n", " 3 - 0 - 1 - 4: 180.00000 Degree\n", " 3 - 0 - 1 - 5: 58.90241 Degree\n", " 3 - 0 - 1 - 6: 58.90241 Degree\n", " 4 - 0 - 1 - 5: 121.09759 Degree\n", " 4 - 0 - 1 - 6: 121.09759 Degree\n", " 5 - 0 - 1 - 6: 117.80483 Degree\n", " 1 - 0 - 4 - 5: 121.06434 Degree\n", " 1 - 0 - 4 - 6: 121.06434 Degree\n", " 5 - 0 - 4 - 6: 117.87131 Degree\n", " 1 - 0 - 5 - 4: 121.03351 Degree\n", " 1 - 0 - 5 - 6: 119.43447 Degree\n", " 4 - 0 - 5 - 6: 119.53201 Degree\n", " 1 - 0 - 6 - 4: 121.03351 Degree\n", " 1 - 0 - 6 - 5: 119.43447 Degree\n", " 4 - 0 - 6 - 5: 119.53201 Degree\n", " 0 - 1 - 2 - 3: 180.00000 Degree\n", " 0 - 1 - 3 - 2: 180.00000 Degree\n" ] } ], "source": [ "sol_mole.print_dihedral_angle()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 6: Center-of-Mass Translation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Find the center of mass of the molecule:\n", "\n", "$$\n", "\\boldsymbol{r}_\\mathrm{c.m.} = \\left. \\sum_{i} m_i \\boldsymbol{r}_i \\right/ \\sum_{i} m_i\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "or in its components\n", "\n", "$$\n", "x_\\mathrm{c.m.} = \\left. \\sum_i m_i x_i \\right/ \\sum_i m_i , \\quad\n", "y_\\mathrm{c.m.} = \\left. \\sum_i m_i y_i \\right/ \\sum_i m_i , \\quad\n", "z_\\mathrm{c.m.} = \\left. \\sum_i m_i z_i \\right/ \\sum_i m_i\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $m_i$ is the mass of atom $i$ and the summation runs over all atoms in the molecule.\n", "\n", "Translate the input coordinates of the molecule to the center-of-mass.\n", "\n", "We use isotopic-averaged atom mass, so the final results may be different to the original Crawford's project at $10^{-3} \\, \\mathsf{a.u.}$ level." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 1: Atomic masses\n", ":class: dropdown\n", "\n", "An excellent source for atomic masses and other physical constants is the [National Institute of Standard and Technology (NIST) website](https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl).\n", "\n", "For python users, one may find package [molmass](https://github.com/cgohlke/molmass) useful.\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implementation" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def center_of_mass(mole: Molecule) -> np.ndarray or list:\n", " # Output: Center of mass for this molecule\n", " raise NotImplementedError(\"About 5~10 lines of code\")\n", "\n", "Molecule.center_of_mass = center_of_mass" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "array([0.64475065, 0. , 2.31636762])" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.center_of_mass()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 7: Principal Moments of Inertia" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate elements of the [moment of inertia tensor](http://en.wikipedia.org/wiki/Moment_of_inertia_tensor).\n", "\n", "Diagonal:\n", "\n", "$$\n", "I_{xx} = \\sum_i m_i (\\tilde y_i^2 + \\tilde z_i^2), \\quad I_{yy} = \\sum_i m_i (\\tilde z_i^2 + \\tilde x_i^2) \\, \\quad I_{zz} = \\sum_i m_i (\\tilde x_i^2 + \\tilde y_i^2)\n", "$$\n", "\n", "Off-diagonal (add a negative sign):\n", "\n", "$$\n", "I_{xy} = I_{yx} = - \\sum_i m_i \\tilde x_i \\tilde y_i, \\quad I_{yz} = I_{zy} = - \\sum_i m_i \\tilde y_i \\tilde z_i, \\quad I_{zx} = I_{xz} = - \\sum_i m_i \\tilde z_i \\tilde x_i\n", "$$\n", "\n", "Calculate eigenvalue of inertia tensor to obtain the principal moments of inertia:\n", "\n", "$$\n", "I_a \\leqslant I_b \\leqslant I_c\n", "$$\n", "\n", "Report the moments of inertia in $\\mathsf{amu \\cdot Bohr^2}$, $\\mathsf{amu \\cdot {\\buildrel_{\\circ} \\over{\\mathsf{A}}}{}^2}$, $\\mathsf{g \\cdot cm^2}$.\n", "\n", "Note that $\\tilde x_i, \\tilde y_i, \\tilde z_i$ in the equations above are translated by center of mass:\n", "\n", "$$\n", "\\tilde x_i = x_i - x_\\mathrm{c.m.}, \\quad \\tilde y_i = y_i - y_\\mathrm{c.m.}, \\quad \\tilde z_i = z_i - z_\\mathrm{c.m.}\n", "$$\n", "\n", "Based on the relative values of the principal moments, determine the [molecular rotor type](http://en.wikipedia.org/wiki/Rotational_spectroscopy). Criterion of equality and far larger/smaller could be set to $10^{-4} \\mathsf{amu \\cdot Bohr^2}$.\n", "\n", "- Spherical: $I_a = I_b = I_c$\n", "- Linear: $I_a \\ll I_b = I_c$\n", "- Prolate: $I_a < I_b = I_c$\n", "- Oblate: $I_a = I_b < I_c$\n", "\n", "- Asymmetric: $I_a \\neq I_b \\neq I_c$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 1: Eigenvalue of symmetric matrix\n", ":class: dropdown\n", "\n", "Eigenvalue $\\lambda$ of moments of inertia matrix can be obtained as\n", "\n", "$$\n", "\\left\\vert\n", "\\begin{matrix}\n", "I_{xx} - \\lambda & I_{xy} & I_{xz} \\\\\n", "I_{yx} & I_{yy} - \\lambda & I_{yz} \\\\\n", "I_{zx} & I_{zy} & I_{zz} - \\lambda\n", "\\end{matrix}\n", "\\right\\vert = 0\n", "$$\n", "\n", "This leads to a cubic equation in $\\lambda$, which one can solve directly.\n", "\n", "However, solving symmetric matrix by program is usually called *diagonalization*. This can be done by `numpy.linalg.eigh` which returns eigenvalues and eigenvectors ([NumPy API](https://numpy.org/doc/stable/reference/generated/numpy.linalg.eigh.html)) or `numpy.linalg.eigvalsh` which only returns eigenvalues ([NumPy API](https://numpy.org/doc/stable/reference/generated/numpy.linalg.eigvalsh.html)).\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 2: Physical constants\n", ":class: dropdown\n", "\n", "Lots of useful and precise physical constants are available at the [National Institute of Standards and Technology website](https://physics.nist.gov/cuu/Constants/index.html?/codata86.html).\n", "\n", "For python, `scipy.constants` ([SciPy API](https://docs.scipy.org/doc/scipy/reference/constants.html)) offers lots of constants and unit conversion utilities.\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": 22, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def moment_of_inertia(mole: Molecule) -> np.ndarray or list:\n", " # Output: Principle of moments of intertia\n", " raise NotImplementedError(\"About 4~25 lines of code\")\n", "\n", "Molecule.moment_of_inertia = moment_of_inertia" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def print_moment_of_interia(mole: Molecule):\n", " # Output: Print moments of inertia in amu bohr2, amu Å2, and g cm2\n", " raise NotImplementedError(\"About 3~15 lines of code\")\n", "\n", "Molecule.print_moment_of_interia = print_moment_of_interia" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def type_of_moment_of_interia(mole: Molecule) -> str:\n", " # Output: Judge which type of moment of interia is\n", " raise NotImplementedError(\"About 7~11 lines of code\")\n", "\n", "Molecule.type_of_moment_of_interia = type_of_moment_of_interia" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Value of principles of moment of interia:" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "In amu bohr^2: 3.19751875e+01 1.78736281e+02 1.99467661e+02\n", "In amu angstrom^2: 8.95396445e+00 5.00512562e+01 5.58566339e+01\n", "In g cm^2: 1.48684078e-39 8.31120663e-39 9.27521227e-39\n" ] } ], "source": [ "sol_mole.print_moment_of_interia()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Type of interia:" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "'Asymmetric'" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.type_of_moment_of_interia()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 8: Rotational Constants" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute the rotational constants in $\\mathsf{cm}^{-1}$ and $\\mathsf{MHz}$:\n", "\n", "$$\n", "A = \\frac{h}{8 \\pi^2 c I_a} , \\quad B = \\frac{h}{8 \\pi^2 c I_b} , \\quad C = \\frac{h}{8 \\pi^2 c I_c}\n", "$$\n", "\n", "where $A \\geqslant B \\geqslant C$." ] }, { "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": 27, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def rotational_constants(mole: Molecule) -> np.ndarray or list:\n", " # Output: Rotational constant in cm^-1\n", " raise NotImplementedError(\"About 1~10 lines of code\")\n", "\n", "Molecule.rotational_constants = rotational_constants" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute the rotational constants in $\\mathsf{cm}^{-1}$" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "array([1.88270004, 0.33680731, 0.30180174])" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.rotational_constants()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Test Cases" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Acetaldehyde**\n", "- Input file: {download}`input/acetaldehyde.dat`" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Atom Charges ===\n", "[6 6 8 1 1 1 1]\n", "=== Coordinates ===\n", "[[ 0. 0. 0. ]\n", " [ 0. 0. 2.84511213]\n", " [ 1.89911596 0. 4.13906253]\n", " [-1.89404831 0. 3.74768867]\n", " [ 1.94250082 0. -0.70114598]\n", " [-1.00729547 -1.66997184 -0.70591697]\n", " [-1.00729547 1.66997184 -0.70591697]]\n", "=== Bond Length ===\n", " 0 - 1: 2.84511 Bohr\n", " 0 - 2: 4.55395 Bohr\n", " 0 - 3: 4.19912 Bohr\n", " 0 - 4: 2.06517 Bohr\n", " 0 - 5: 2.07407 Bohr\n", " 0 - 6: 2.07407 Bohr\n", " 1 - 2: 2.29803 Bohr\n", " 1 - 3: 2.09811 Bohr\n", " 1 - 4: 4.04342 Bohr\n", " 1 - 5: 4.05133 Bohr\n", " 1 - 6: 4.05133 Bohr\n", " 2 - 3: 3.81330 Bohr\n", " 2 - 4: 4.84040 Bohr\n", " 2 - 5: 5.89151 Bohr\n", " 2 - 6: 5.89151 Bohr\n", " 3 - 4: 5.87463 Bohr\n", " 3 - 5: 4.83836 Bohr\n", " 3 - 6: 4.83836 Bohr\n", " 4 - 5: 3.38971 Bohr\n", " 4 - 6: 3.38971 Bohr\n", " 5 - 6: 3.33994 Bohr\n", "=== Bond Angle ===\n", " 0 - 1 - 2: 124.26831 Degree\n", " 0 - 1 - 3: 115.47934 Degree\n", " 1 - 0 - 4: 109.84706 Degree\n", " 1 - 0 - 5: 109.89841 Degree\n", " 1 - 0 - 6: 109.89841 Degree\n", " 4 - 0 - 5: 109.95368 Degree\n", " 4 - 0 - 6: 109.95368 Degree\n", " 5 - 0 - 6: 107.25265 Degree\n", " 2 - 1 - 3: 120.25235 Degree\n", "=== Out-of-Plane Angle ===\n", " 3 - 0 - 1 - 2: 0.00000 Degree\n", " 2 - 0 - 1 - 3: 0.00000 Degree\n", " 5 - 1 - 0 - 4: -53.62632 Degree\n", " 6 - 1 - 0 - 4: 53.62632 Degree\n", " 4 - 1 - 0 - 5: 53.65153 Degree\n", " 6 - 1 - 0 - 5: -56.27711 Degree\n", " 4 - 1 - 0 - 6: -53.65153 Degree\n", " 5 - 1 - 0 - 6: 56.27711 Degree\n", " 1 - 4 - 0 - 5: -53.67878 Degree\n", " 6 - 4 - 0 - 5: 56.19462 Degree\n", " 1 - 4 - 0 - 6: 53.67878 Degree\n", " 5 - 4 - 0 - 6: -56.19462 Degree\n", " 1 - 5 - 0 - 6: -54.97706 Degree\n", " 4 - 5 - 0 - 6: 54.86999 Degree\n", " 0 - 2 - 1 - 3: 0.00000 Degree\n", "=== Dihedral Angle ===\n", " 2 - 0 - 1 - 3: 180.00000 Degree\n", " 2 - 0 - 1 - 4: 0.00000 Degree\n", " 2 - 0 - 1 - 5: 121.09759 Degree\n", " 2 - 0 - 1 - 6: 121.09759 Degree\n", " 3 - 0 - 1 - 4: 180.00000 Degree\n", " 3 - 0 - 1 - 5: 58.90241 Degree\n", " 3 - 0 - 1 - 6: 58.90241 Degree\n", " 4 - 0 - 1 - 5: 121.09759 Degree\n", " 4 - 0 - 1 - 6: 121.09759 Degree\n", " 5 - 0 - 1 - 6: 117.80483 Degree\n", " 1 - 0 - 4 - 5: 121.06434 Degree\n", " 1 - 0 - 4 - 6: 121.06434 Degree\n", " 5 - 0 - 4 - 6: 117.87131 Degree\n", " 1 - 0 - 5 - 4: 121.03351 Degree\n", " 1 - 0 - 5 - 6: 119.43447 Degree\n", " 4 - 0 - 5 - 6: 119.53201 Degree\n", " 1 - 0 - 6 - 4: 121.03351 Degree\n", " 1 - 0 - 6 - 5: 119.43447 Degree\n", " 4 - 0 - 6 - 5: 119.53201 Degree\n", " 0 - 1 - 2 - 3: 180.00000 Degree\n", " 0 - 1 - 3 - 2: 180.00000 Degree\n", "=== Center of Mass ===\n", " 0.64475 0.00000 2.31637\n", "=== Moments of Inertia ===\n", "In amu bohr^2: 3.19751875e+01 1.78736281e+02 1.99467661e+02\n", "In amu angstrom^2: 8.95396445e+00 5.00512562e+01 5.58566339e+01\n", "In g cm^2: 1.48684078e-39 8.31120663e-39 9.27521227e-39\n", "Type: Asymmetric\n", "=== Rotational Constants ===\n", "In cm^-1: 1.88270 0.33681 0.30180\n", "In MHz: 56441.92712 10097.22927 9047.78849\n" ] } ], "source": [ "sol_mole = SolMol()\n", "sol_mole.construct_from_dat_file(\"input/acetaldehyde.dat\")\n", "sol_mole.print_solution_01()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Allene**\n", "- Input file: {download}`input/allene.dat`" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Atom Charges ===\n", "[6 6 6 1 1 1 1]\n", "=== Coordinates ===\n", "[[ 0. 0. 1.88972599]\n", " [ 2.55113008 0. 1.88972599]\n", " [-2.55113008 0. 1.88972599]\n", " [ 3.49599308 1.15721611 0.73250988]\n", " [ 3.49599308 -1.15721611 3.0469421 ]\n", " [-3.49599308 -1.15721611 0.73250988]\n", " [-3.49599308 1.15721611 3.0469421 ]]\n", "=== Bond Length ===\n", " 0 - 1: 2.55113 Bohr\n", " 0 - 2: 2.55113 Bohr\n", " 0 - 3: 3.86009 Bohr\n", " 0 - 4: 3.86009 Bohr\n", " 0 - 5: 3.86009 Bohr\n", " 0 - 6: 3.86009 Bohr\n", " 1 - 2: 5.10226 Bohr\n", " 1 - 3: 1.88973 Bohr\n", " 1 - 4: 1.88973 Bohr\n", " 1 - 5: 6.26466 Bohr\n", " 1 - 6: 6.26466 Bohr\n", " 2 - 3: 6.26466 Bohr\n", " 2 - 4: 6.26466 Bohr\n", " 2 - 5: 1.88973 Bohr\n", " 2 - 6: 1.88973 Bohr\n", " 3 - 4: 3.27310 Bohr\n", " 3 - 5: 7.36508 Bohr\n", " 3 - 6: 7.36508 Bohr\n", " 4 - 5: 7.36508 Bohr\n", " 4 - 6: 7.36508 Bohr\n", " 5 - 6: 3.27310 Bohr\n", "=== Bond Angle ===\n", " 1 - 0 - 2: 180.00000 Degree\n", " 0 - 1 - 3: 120.00000 Degree\n", " 0 - 1 - 4: 120.00000 Degree\n", " 0 - 2 - 5: 120.00000 Degree\n", " 0 - 2 - 6: 120.00000 Degree\n", " 3 - 1 - 4: 120.00000 Degree\n", " 5 - 2 - 6: 120.00000 Degree\n", "=== Out-of-Plane Angle ===\n", " 4 - 0 - 1 - 3: 0.00000 Degree\n", " 3 - 0 - 1 - 4: 0.00000 Degree\n", " 6 - 0 - 2 - 5: 0.00000 Degree\n", " 5 - 0 - 2 - 6: 0.00000 Degree\n", " 0 - 3 - 1 - 4: 0.00000 Degree\n", " 0 - 5 - 2 - 6: 0.00000 Degree\n", "=== Dihedral Angle ===\n", " 2 - 0 - 1 - 3: 90.00000 Degree\n", " 2 - 0 - 1 - 4: 90.00000 Degree\n", " 3 - 0 - 1 - 4: 180.00000 Degree\n", " 1 - 0 - 2 - 5: 90.00000 Degree\n", " 1 - 0 - 2 - 6: 90.00000 Degree\n", " 5 - 0 - 2 - 6: 180.00000 Degree\n", " 0 - 1 - 3 - 4: 180.00000 Degree\n", " 0 - 1 - 4 - 3: 180.00000 Degree\n", " 0 - 2 - 5 - 6: 180.00000 Degree\n", " 0 - 2 - 6 - 5: 180.00000 Degree\n", "=== Center of Mass ===\n", " 0.00000 0.00000 1.88973\n", "=== Moments of Inertia ===\n", "In amu bohr^2: 1.07982664e+01 2.11013373e+02 2.11013373e+02\n", "In amu angstrom^2: 3.02382256e+00 5.90897626e+01 5.90897626e+01\n", "In g cm^2: 5.02117550e-40 9.81208592e-39 9.81208592e-39\n", "Type: Prolate\n", "=== Rotational Constants ===\n", "In cm^-1: 5.57494 0.28529 0.28529\n", "In MHz: 167132.49483 8552.73379 8552.73379\n" ] } ], "source": [ "sol_mole = SolMol()\n", "sol_mole.construct_from_dat_file(\"input/allene.dat\")\n", "sol_mole.print_solution_01()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Benzene**\n", "- Input file: {download}`input/benzene.dat`" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "scrolled": false, "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Atom Charges ===\n", "[6 6 6 6 6 6 1 1 1 1 1 1]\n", "=== Coordinates ===\n", "[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [ 0.00000000e+00 0.00000000e+00 2.61644846e+00]\n", " [ 2.26591048e+00 0.00000000e+00 3.92467249e+00]\n", " [ 4.53182108e+00 0.00000000e+00 2.61644839e+00]\n", " [ 4.53182108e+00 0.00000000e+00 -1.32281000e-07]\n", " [ 2.26591069e+00 0.00000000e+00 -1.30822415e+00]\n", " [ 2.26591069e+00 -0.00000000e+00 -3.33421820e+00]\n", " [ 6.28638327e+00 0.00000000e+00 -1.01299708e+00]\n", " [ 6.28638324e+00 0.00000000e+00 3.62944532e+00]\n", " [ 2.26591048e+00 0.00000000e+00 5.95066654e+00]\n", " [-1.75456219e+00 0.00000000e+00 3.62944541e+00]\n", " [-1.75456216e+00 -0.00000000e+00 -1.01299693e+00]]\n", "=== Bond Length ===\n", " 0 - 1: 2.61645 Bohr\n", " 0 - 2: 4.53182 Bohr\n", " 0 - 3: 5.23290 Bohr\n", " 0 - 4: 4.53182 Bohr\n", " 0 - 5: 2.61645 Bohr\n", " 0 - 6: 4.03130 Bohr\n", " 0 - 7: 6.36748 Bohr\n", " 0 - 8: 7.25889 Bohr\n", " 0 - 9: 6.36748 Bohr\n", " 0 - 10: 4.03130 Bohr\n", " 0 - 11: 2.02599 Bohr\n", " 1 - 2: 2.61645 Bohr\n", " 1 - 3: 4.53182 Bohr\n", " 1 - 4: 5.23290 Bohr\n", " 1 - 5: 4.53182 Bohr\n", " 1 - 6: 6.36748 Bohr\n", " 1 - 7: 7.25889 Bohr\n", " 1 - 8: 6.36748 Bohr\n", " 1 - 9: 4.03130 Bohr\n", " 1 - 10: 2.02599 Bohr\n", " 1 - 11: 4.03130 Bohr\n", " 2 - 3: 2.61645 Bohr\n", " 2 - 4: 4.53182 Bohr\n", " 2 - 5: 5.23290 Bohr\n", " 2 - 6: 7.25889 Bohr\n", " 2 - 7: 6.36748 Bohr\n", " 2 - 8: 4.03130 Bohr\n", " 2 - 9: 2.02599 Bohr\n", " 2 - 10: 4.03130 Bohr\n", " 2 - 11: 6.36748 Bohr\n", " 3 - 4: 2.61645 Bohr\n", " 3 - 5: 4.53182 Bohr\n", " 3 - 6: 6.36748 Bohr\n", " 3 - 7: 4.03130 Bohr\n", " 3 - 8: 2.02599 Bohr\n", " 3 - 9: 4.03130 Bohr\n", " 3 - 10: 6.36748 Bohr\n", " 3 - 11: 7.25889 Bohr\n", " 4 - 5: 2.61645 Bohr\n", " 4 - 6: 4.03130 Bohr\n", " 4 - 7: 2.02599 Bohr\n", " 4 - 8: 4.03130 Bohr\n", " 4 - 9: 6.36748 Bohr\n", " 4 - 10: 7.25889 Bohr\n", " 4 - 11: 6.36748 Bohr\n", " 5 - 6: 2.02599 Bohr\n", " 5 - 7: 4.03130 Bohr\n", " 5 - 8: 6.36748 Bohr\n", " 5 - 9: 7.25889 Bohr\n", " 5 - 10: 6.36748 Bohr\n", " 5 - 11: 4.03130 Bohr\n", " 6 - 7: 4.64244 Bohr\n", " 6 - 8: 8.04095 Bohr\n", " 6 - 9: 9.28488 Bohr\n", " 6 - 10: 8.04095 Bohr\n", " 6 - 11: 4.64244 Bohr\n", " 7 - 8: 4.64244 Bohr\n", " 7 - 9: 8.04095 Bohr\n", " 7 - 10: 9.28488 Bohr\n", " 7 - 11: 8.04095 Bohr\n", " 8 - 9: 4.64244 Bohr\n", " 8 - 10: 8.04095 Bohr\n", " 8 - 11: 9.28488 Bohr\n", " 9 - 10: 4.64244 Bohr\n", " 9 - 11: 8.04095 Bohr\n", " 10 - 11: 4.64244 Bohr\n", "=== Bond Angle ===\n", " 0 - 1 - 2: 120.00000 Degree\n", " 1 - 0 - 5: 120.00000 Degree\n", " 0 - 1 - 10: 120.00000 Degree\n", " 1 - 0 - 11: 120.00000 Degree\n", " 0 - 5 - 4: 120.00000 Degree\n", " 0 - 5 - 6: 120.00000 Degree\n", " 5 - 0 - 11: 120.00000 Degree\n", " 1 - 2 - 3: 120.00000 Degree\n", " 1 - 2 - 9: 120.00000 Degree\n", " 2 - 1 - 10: 120.00000 Degree\n", " 2 - 3 - 4: 120.00000 Degree\n", " 2 - 3 - 8: 120.00000 Degree\n", " 3 - 2 - 9: 120.00000 Degree\n", " 3 - 4 - 5: 120.00000 Degree\n", " 3 - 4 - 7: 120.00000 Degree\n", " 4 - 3 - 8: 120.00000 Degree\n", " 4 - 5 - 6: 120.00000 Degree\n", " 5 - 4 - 7: 120.00000 Degree\n", "=== Out-of-Plane Angle ===\n", " 10 - 0 - 1 - 2: 0.00000 Degree\n", " 11 - 1 - 0 - 5: 0.00000 Degree\n", " 2 - 0 - 1 - 10: 0.00000 Degree\n", " 5 - 1 - 0 - 11: 0.00000 Degree\n", " 6 - 0 - 5 - 4: 0.00000 Degree\n", " 4 - 0 - 5 - 6: 0.00000 Degree\n", " 1 - 5 - 0 - 11: 0.00000 Degree\n", " 9 - 1 - 2 - 3: 0.00000 Degree\n", " 3 - 1 - 2 - 9: 0.00000 Degree\n", " 0 - 2 - 1 - 10: 0.00000 Degree\n", " 8 - 2 - 3 - 4: 0.00000 Degree\n", " 4 - 2 - 3 - 8: 0.00000 Degree\n", " 1 - 3 - 2 - 9: 0.00000 Degree\n", " 7 - 3 - 4 - 5: 0.00000 Degree\n", " 5 - 3 - 4 - 7: 0.00000 Degree\n", " 2 - 4 - 3 - 8: 0.00000 Degree\n", " 0 - 4 - 5 - 6: 0.00000 Degree\n", " 3 - 5 - 4 - 7: 0.00000 Degree\n", "=== Dihedral Angle ===\n", " 2 - 0 - 1 - 5: 0.00000 Degree\n", " 2 - 0 - 1 - 10: 180.00000 Degree\n", " 2 - 0 - 1 - 11: 180.00000 Degree\n", " 5 - 0 - 1 - 10: 180.00000 Degree\n", " 5 - 0 - 1 - 11: 180.00000 Degree\n", " 10 - 0 - 1 - 11: 0.00000 Degree\n", " 1 - 0 - 5 - 4: 0.00000 Degree\n", " 1 - 0 - 5 - 6: 180.00000 Degree\n", " 1 - 0 - 5 - 11: 180.00000 Degree\n", " 4 - 0 - 5 - 6: 180.00000 Degree\n", " 4 - 0 - 5 - 11: 180.00000 Degree\n", " 6 - 0 - 5 - 11: 0.00000 Degree\n", " 1 - 0 - 11 - 5: 180.00000 Degree\n", " 0 - 1 - 2 - 3: 0.00000 Degree\n", " 0 - 1 - 2 - 9: 180.00000 Degree\n", " 0 - 1 - 2 - 10: 180.00000 Degree\n", " 3 - 1 - 2 - 9: 180.00000 Degree\n", " 3 - 1 - 2 - 10: 180.00000 Degree\n", " 9 - 1 - 2 - 10: 0.00000 Degree\n", " 0 - 1 - 10 - 2: 180.00000 Degree\n", " 1 - 2 - 3 - 4: 0.00000 Degree\n", " 1 - 2 - 3 - 8: 180.00000 Degree\n", " 1 - 2 - 3 - 9: 180.00000 Degree\n", " 4 - 2 - 3 - 8: 180.00000 Degree\n", " 4 - 2 - 3 - 9: 180.00000 Degree\n", " 8 - 2 - 3 - 9: 0.00000 Degree\n", " 1 - 2 - 9 - 3: 180.00000 Degree\n", " 2 - 3 - 4 - 5: 0.00000 Degree\n", " 2 - 3 - 4 - 7: 180.00000 Degree\n", " 2 - 3 - 4 - 8: 180.00000 Degree\n", " 5 - 3 - 4 - 7: 180.00000 Degree\n", " 5 - 3 - 4 - 8: 180.00000 Degree\n", " 7 - 3 - 4 - 8: 0.00000 Degree\n", " 2 - 3 - 8 - 4: 180.00000 Degree\n", " 0 - 4 - 5 - 3: 0.00000 Degree\n", " 0 - 4 - 5 - 6: 180.00000 Degree\n", " 0 - 4 - 5 - 7: 180.00000 Degree\n", " 3 - 4 - 5 - 6: 180.00000 Degree\n", " 3 - 4 - 5 - 7: 180.00000 Degree\n", " 6 - 4 - 5 - 7: 0.00000 Degree\n", " 3 - 4 - 7 - 5: 180.00000 Degree\n", " 0 - 5 - 6 - 4: 180.00000 Degree\n", "=== Center of Mass ===\n", " 2.26591 0.00000 1.30822\n", "=== Moments of Inertia ===\n", "In amu bohr^2: 3.11839639e+02 3.11839704e+02 6.23679344e+02\n", "In amu angstrom^2: 8.73239929e+01 8.73240110e+01 1.74648004e+02\n", "In g cm^2: 1.45004902e-38 1.45004932e-38 2.90009833e-38\n", "Type: Oblate\n", "=== Rotational Constants ===\n", "In cm^-1: 0.19305 0.19305 0.09652\n", "In MHz: 5787.40152 5787.40032 2893.70046\n" ] } ], "source": [ "sol_mole = SolMol()\n", "sol_mole.construct_from_dat_file(\"input/benzene.dat\")\n", "sol_mole.print_solution_01()" ] }, { "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 }