Tutorial_CH4_molecule_IR_Raman_symmetry.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "880e2b4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "import os\n",
    "import pickle\n",
    "import random\n",
    "import sys\n",
    "import pickle\n",
    "import jraph\n",
    "import tqdm\n",
    "\n",
    "from typing import Callable, Dict, List, Optional\n",
    "import haiku as hk\n",
    "import ase\n",
    "import ase.io\n",
    "import jax\n",
    "import jax.numpy as jnp\n",
    "import numpy as np\n",
    "import optax\n",
    "import yaml\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from functools import partial\n",
    "\n",
    "import sys\n",
    "import os\n",
    "\n",
    "\n",
    "from phonax.datasets import datasets\n",
    "from phonax.utils import (\n",
    "    create_directory_with_random_name,\n",
    "    compute_avg_num_neighbors,\n",
    ")\n",
    "from phonax.data_utils import (\n",
    "    get_atomic_number_table_from_zs,\n",
    "    compute_average_E0s,\n",
    ")\n",
    "from phonax.predictors import predict_energy_forces_stress\n",
    "from phonax.optimizer import optimizer\n",
    "from phonax.energy_force_train import energy_force_train\n",
    "from phonax.loss import WeightedEnergyFrocesStressLoss\n",
    "\n",
    "\n",
    "from phonax.trained_models import NequIP_JAXMD_molecule_model\n",
    "\n",
    "from phonax.phonons import (\n",
    "    atoms_to_ext_graph,\n",
    "    predict_molecular_spectra_unfilter,\n",
    ")\n",
    "\n",
    "from posym import SymmetryNormalModes\n",
    "\n",
    "jax.config.update(\"jax_debug_nans\", True)\n",
    "jax.config.update(\"jax_debug_infs\", True)\n",
    "np.set_printoptions(precision=3, suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97cc7113",
   "metadata": {},
   "source": [
    "# Molecular symmetry analysis for IR / Raman active modes\n",
    "\n",
    "In this tutorial, we demonstrate how to perform symmetry analysis for vibrational modes of the molecules and determine the IR / Raman active modes.\n",
    "We use CH$_4$ molecule as the example (as used in the book on group theory, application to the physics of condensed matter, by M. Dresselhaus, G. Dresselhaus and A. Jorio.)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7afa2c2",
   "metadata": {},
   "source": [
    "## Pre-trained NequIP energy model for molecules\n",
    "\n",
    "We we load the pre-trained NequIP model for molecules.\n",
    "This energy model is trained with energy and force data selected from universal IAP dataset which contain only H, C, N, O and F elements that are common for organic molecules. The DFT calculations are done with the PBE exchange-correlation functionals.\n",
    "This simplified model is used to discuss the symmetry analysis for IR / Raman active modes, and the model can be improved by more advanced DFT (hybrid functional) or quantum chemistry methods (such as CCSD).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "97253f76",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Create NequIP (JAX-MD version) with parameters {'use_sc': True, 'graph_net_steps': 2, 'hidden_irreps': '96x0e+ 96x0o  + 64x1e +64x1o +64x2e +64x2o', 'nonlinearities': {'e': 'swish', 'o': 'tanh'}, 'r_max': 5.0, 'avg_num_neighbors': 53.2808396682792, 'avg_r_min': None, 'num_species': 100, 'radial_basis': <function bessel_basis at 0x7fd0775145e0>, 'radial_envelope': <function soft_envelope at 0x7fd077514540>}\n"
     ]
    }
   ],
   "source": [
    "model_fn, params, num_message_passing, r_max= NequIP_JAXMD_molecule_model(os.path.join(os.getcwd(), 'trained-models'))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a647c64",
   "metadata": {},
   "source": [
    "Next, we load the molecular graph data and compute the molecular vibrational modes and spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e887a141",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Molecular vibration mode energies (in cm^-1):\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Array([  -0.17 ,    0.061,    0.279,    8.22 ,    8.23 ,    8.232,\n",
       "       1246.004, 1246.01 , 1342.545, 1342.547, 1342.55 , 3277.403,\n",
       "       3277.416, 3277.425, 3305.425], dtype=float32)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mol_atoms = ase.io.read('data/molecules/CH4-symmetry.vasp')\n",
    "\n",
    "mol_graph = atoms_to_ext_graph(mol_atoms, r_max, num_message_passing=1)\n",
    "mol_masses = ase.data.atomic_masses[mol_graph.nodes.species][mol_graph.nodes.mask_primitive]\n",
    "gamma_eigs, gamma_vecs, D_gamma = predict_molecular_spectra_unfilter(params, model_fn,mol_graph,mol_masses)\n",
    "\n",
    "print('Molecular vibration mode energies (in cm^-1):')\n",
    "# in cm^-1\n",
    "gamma_eigs\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ca34ac4",
   "metadata": {},
   "source": [
    "## Molecular point group symmetry\n",
    "\n",
    "Now we can use posym tool to figure out the molecular point group symmetry.\n",
    "\n",
    "https://github.com/abelcarreras/posym\n",
    "\n",
    "\n",
    "To be confirmed: the output Continuous symmetry measures (CSM) probes the amount of symmetry break with respect to a point group symmetry? We should select the point group (highest symmetry) with small measure?\n",
    "Any other more direct way to get the point group symmetry? (pymatgen might have a relevant tool for this, but needs to figure out how to convert the molecular structure object.)\n",
    "\n",
    "Response from Abel Carreras confirms this:\n",
    "\n",
    "Using Posym, indeed you can also determine the symmetry point group a molecule belongs by computing the CSM\n",
    "and looking for the group that gives a value close to zero. However, in general this requires to analyze\n",
    "all the possible groups one by one which can be quite inefficient compared to other methods.\n",
    "\n",
    "Also yes, you will find zero CMS value for all subgroups (since the molecule belongs to all of them),\n",
    "then is your criteria to choose the highest symmetry group among them. The advantage of using posym\n",
    "in this case is that you have a clear measure (the CSM) that you can use as a criteria to determine if the\n",
    "molecule belongs or not to the group given a certain tolerance (concept of quasi-symmetry).\n",
    "\n",
    "There are many ways (and software) to determine the point group.\n",
    "Usually they are based in the determination of the principal axis of inertia to reduce rule out\n",
    "different symmetry groups and then looking for specific symmetry operations to restrict the\n",
    "search until the group is found. Very similar to the usual scheme that you can find in chemistry books.\n",
    "\n",
    "Below, we can also use another method from his repository:\n",
    "\n",
    "https://github.com/abelcarreras/pointgroup\n",
    "\n",
    "\n",
    "\n",
    "Regarding the proper citation(s) for his work:\n",
    "\n",
    "At this moment I don’t have any paper specific to this library. That’s why I created a citable zenodo\n",
    "repository:\n",
    "\n",
    "https://zenodo.org/record/7261326\n",
    "\n",
    "Here I list some papers related to the methodology employed in this library:\n",
    "\n",
    "Pinsky M, Dryzun C, Casanova D, Alemany P, Avnir D, J Comput Chem. 29:2712-21 (2008) \n",
    "Pinsky M, Casanova D, Alemany P, Alvarez S, Avnir D, Dryzun C, Kizner Z, Sterkin A. J Comput Chem. 29:190-7 (2008)\n",
    "Casanova D, Alemany P. Phys Chem Chem Phys. 12(47):15523–9 (2010)\n",
    "Casanova D, Alemany P, Falceto A, Carreras A, Alvarez S. J Comput Chem 34(15):1321–31 (2013)\n",
    "A. Carreras, E. Bernuz, X. Marugan, M. Llunell, P. Alemany, Chem. Eur. J. 25, 673 – 691 (2019) \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "be3be14c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Symmetry measure Td :  4.3256913606626313e-07\n",
      "Symmetry measure C3v :  4.022350230670213e-08\n",
      "Symmetry measure C4v :  10.86539664434758\n",
      "Symmetry measure D6h :  24.602044733638294\n"
     ]
    }
   ],
   "source": [
    "from posym import SymmetryMolecule\n",
    "pos_CM = mol_atoms.get_center_of_mass()\n",
    "coordinates = mol_atoms.get_positions() - pos_CM\n",
    "\n",
    "symbols = mol_atoms.get_chemical_symbols()\n",
    "\n",
    "sym_geom = SymmetryMolecule(group='Td', coordinates=coordinates, symbols=symbols)\n",
    "print('Symmetry measure Td : ', sym_geom.measure)\n",
    "\n",
    "# C3v as a subgroup of D6h, the measure will be small as well.\n",
    "sym_geom = SymmetryMolecule(group='C3v', coordinates=coordinates, symbols=symbols)\n",
    "print('Symmetry measure C3v : ', sym_geom.measure)\n",
    "\n",
    "sym_geom = SymmetryMolecule(group='C4v', coordinates=coordinates, symbols=symbols)\n",
    "print('Symmetry measure C4v : ', sym_geom.measure)\n",
    "\n",
    "# this is our symmetry for the benzene\n",
    "sym_geom = SymmetryMolecule(group='D6h', coordinates=coordinates, symbols=symbols)\n",
    "print('Symmetry measure D6h : ', sym_geom.measure)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80428d0f",
   "metadata": {},
   "source": [
    "Alternatively, we can use the pointgroup library by the same author of the library, and we confirm a CH$_4$ molecule has a Td point group symmetry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f5149a83",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Point group:  Td\n"
     ]
    }
   ],
   "source": [
    "from pointgroup import PointGroup\n",
    "\n",
    "pos_CM = mol_atoms.get_center_of_mass()\n",
    "\n",
    "coordinates = mol_atoms.get_positions() - pos_CM\n",
    "\n",
    "symbols = mol_atoms.get_chemical_symbols()\n",
    "\n",
    "pg = PointGroup(positions=coordinates, \n",
    "                symbols=symbols)\n",
    "\n",
    "print('Point group: ', pg.get_point_group())\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9a5e54d",
   "metadata": {},
   "source": [
    "## Translational / Rotational mode projections\n",
    "\n",
    "In the following, I would like to show how one can perform the mode projections to filter out the pure traslational / rotational modes (with respect to the center of mass)\n",
    "\n",
    "Note with the point group symmetry analysis above, the translational modes would carry the vector like symmetry irreducible symmetry representations (irreps), while the rotational modes would be of the anti-symmetric second rank tensor type irreps. These will be relevant later for the IR / Raman active mode analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b801992e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# here we would like to project out the translation + rotation modes by defining the projector\n",
    "\n",
    "# function to orthogonalize the modes we select\n",
    "def vector_gram_schmidt(vecs_in):\n",
    "    vecs_out = []\n",
    "    for vec in vecs_in:\n",
    "        if len(vecs_out) == 0:\n",
    "            vecs_out.append(vec/jnp.linalg.norm(vec))\n",
    "        else:\n",
    "            for vec_proj in vecs_out:\n",
    "                vec -= jnp.sum(vec*vec_proj) *vec_proj\n",
    "            vecs_out.append(vec/jnp.linalg.norm(vec))\n",
    "    return vecs_out\n",
    "\n",
    "num_atoms = len(mol_atoms)\n",
    "\n",
    "mol_cm = mol_atoms.get_center_of_mass()\n",
    "shift_pos = mol_atoms.get_positions() - mol_cm\n",
    "\n",
    "# write down the mode vector, note the subtle sqrt(mass) factor here!\n",
    "tr_mode1 = jnp.zeros((num_atoms,3))\n",
    "tr_mode1 = tr_mode1.at[:,0].add(1.0) * jnp.sqrt(mol_masses)[:,None]\n",
    "tr_mode2 = jnp.zeros((num_atoms,3))\n",
    "tr_mode2 = tr_mode2.at[:,1].add(1.0)* jnp.sqrt(mol_masses)[:,None]\n",
    "tr_mode3 = jnp.zeros((num_atoms,3))\n",
    "tr_mode3 = tr_mode3.at[:,2].add(1.0)* jnp.sqrt(mol_masses)[:,None]\n",
    "\n",
    "\n",
    "tr_mode1 /= jnp.linalg.norm(tr_mode1)\n",
    "tr_mode2 /= jnp.linalg.norm(tr_mode2)\n",
    "tr_mode3 /= jnp.linalg.norm(tr_mode3)\n",
    "\n",
    "[tr_mode1,tr_mode2,tr_mode3] = vector_gram_schmidt([tr_mode1,tr_mode2,tr_mode3])\n",
    "\n",
    "x_axis = jnp.array([1.0,0.0,0.0])\n",
    "y_axis = jnp.array([0.0,1.0,0.0])\n",
    "z_axis = jnp.array([0.0,0.0,1.0])\n",
    "\n",
    "rot_mode1 = jnp.cross(x_axis, shift_pos) * jnp.sqrt(mol_masses)[:,None]\n",
    "rot_mode1 /= jnp.linalg.norm(rot_mode1)\n",
    "rot_mode2 = jnp.cross(y_axis, shift_pos) * jnp.sqrt(mol_masses)[:,None]\n",
    "rot_mode2 /= jnp.linalg.norm(rot_mode2)\n",
    "rot_mode3 = jnp.cross(z_axis, shift_pos) * jnp.sqrt(mol_masses)[:,None]\n",
    "rot_mode3 /= jnp.linalg.norm(rot_mode3)\n",
    "\n",
    "[rot_mode1,rot_mode2,rot_mode3] = vector_gram_schmidt([rot_mode1,rot_mode2,rot_mode3])\n",
    "\n",
    "\n",
    "proj_modes = [tr_mode1,tr_mode2,tr_mode3,rot_mode1,rot_mode2,rot_mode3]\n",
    "\n",
    "\n",
    "proj_mat = jnp.zeros((3*num_atoms, 3*num_atoms))\n",
    "for pmode in proj_modes:\n",
    "    pmode = pmode.reshape(-1)\n",
    "    proj_mat += jnp.einsum('i,j->ij',pmode,pmode)\n",
    "    \n",
    "gamma_proj_norm = []\n",
    "\n",
    "for gamma_vec in gamma_vecs.T:\n",
    "    proj_vec = proj_mat @ gamma_vec.real\n",
    "    gamma_proj_norm.append(jnp.linalg.norm(proj_vec).item())\n",
    "    \n",
    "gamma_proj_norm = np.array(gamma_proj_norm) **2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a91d151",
   "metadata": {},
   "source": [
    "With the total 6 modes for translations and rotations, we can compute the projection weight of these states in the molecular vibrational spectrum. With the projection weight distribution, we find these modes are concentrated at the lowest $E \\sim 0$ zero energy modes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "aace4599",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'projected norm')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9R0lEQVR4nO3deXhU5fnG8Xuy75ElKwZBZN8FwYDWLT9pbUFqtQgIiLhRkCXVItUQqdYgakpVCoW6tlqxdauCIKaggqwJUFB2QayShIAwJEC2Ob8/cAanBMiEyZyZM9/PdeW6zJkz8z4T4nDznPc9r80wDEMAAAAWEWJ2AQAAAN5EuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJYSZnYBvuZwOPTtt98qPj5eNpvN7HIAAEA9GIaho0ePKj09XSEhZ+/NBF24+fbbb5WRkWF2GQAAoAG+/vprXXjhhWc9J+jCTXx8vKSTP5yEhASTqwEAAPVht9uVkZHh+nv8bIIu3DgvRSUkJBBuAAAIMPWZUsKEYgAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYClBd4fixlLrMLR2zyGVHj2h5Pgo9WndVKEhvtuYk/HNHd8fajB7fADwF6aGm08++URPPvmkCgsLtX//fr399tsaPHjwWZ+zfPlyZWdn6/PPP1dGRoYefvhh3X777T6p90wWb9mv6e99of1HTriOpSVGKXdgJ/24SxrjW3x8f6jB7PEBwJ+YelmqoqJC3bt31+zZs+t1/p49e/TTn/5U11xzjTZu3KhJkybpzjvv1JIlSxq50jNbvGW/xv6tyO0vFUkqPnJCY/9WpMVb9jO+hcf3hxrMHh8A/I3NMAzD7CKkkxthnatzM2XKFC1cuFBbtmxxHbv11lt1+PBhLV68uF7j2O12JSYm6siRI+e9cWatw9AVT/z7tL9UfigpPlJv3ttPIY0QIx0O6RdzP9OBo5WMb8L4/lDDuca3SUpNjNKKKddyiQpAQPPk7++AmnOzatUqZWVluR0bMGCAJk2adMbnVFZWqrLy1Ae/3W73Wj1r9xw6a7CRpANHK/WjJ5d5bUxPMb6545tdgyFp/5ETWrvnkDLbNDOlBgDwtYAKN8XFxUpJSXE7lpKSIrvdruPHjys6Ovq05+Tl5Wn69OmNUk/p0bMHG6fwUJtC6rFFu6cchqHq2nM33hi/ccb3hxrqO359f1cBwAoCKtw0xNSpU5Wdne363m63KyMjwyuvnRwfVa/zXrmjb6P8q3nV7oMaOn8145s0vj/UUN/x6/u7CgBWEFD3uUlNTVVJSYnbsZKSEiUkJNTZtZGkyMhIJSQkuH15S5/WTZWWGKUz/XvcppMrVvq0buq1MRnff8b3hxrMHh8A/FFAhZvMzEwVFBS4HVu6dKkyMzNNqSc0xKbcgZ0k6bS/XJzf5w7s1GgTORnf3PH9oQazxwcAf2RquCkvL9fGjRu1ceNGSSeXem/cuFH79u2TdPKS0siRI13n33vvvfryyy/1m9/8Rtu2bdOf/vQnvfHGG5o8ebIZ5UuSftwlTXNuu1Spie5t/9TEKM257dJGv8cI45s7vj/UYPb4AOBvTF0Kvnz5cl1zzTWnHR81apReeukl3X777dq7d6+WL1/u9pzJkyfriy++0IUXXqicnByPbuLnzaXgP2T23WEZ3/y785pdQ63D0G3Pr9Gq3Qc14vKL9MigznRsAFiGJ39/+819bnylscIN4A9+994XemHlHt3zo4s19YaOZpcDAF7jyd/fATXnBsDZJSdEStJZbyoIAFZHuAEsJCnuZLgpJdwACGKEG8BC6NwAAOEGsJSkeGfnhjsSAwhehBvAQpx3Iv7uWLWqahwmVwMA5iDcABZyQXS4wr5f/l1WzqUpAMGJcANYSEiIzXVpink3AIIV4QawmFPzbgg3AIIT4QawmGQ6NwCCHOEGsBhWTAEIdoQbwGKSvl8xxWUpAMGKcANYDJelAAQ7wg1gMUwoBhDsCDeAxTg7N2WEGwBBinADWMwP73NjGIbJ1QCA7xFuAItxhpuqWoeOHK82uRoA8D3CDWAxkWGhSowOl8S8GwDBiXADWBArpgAEM8INYEHcyA9AMCPcABZE5wZAMCPcABbk6tzYCTcAgg/hBrCg5O+3YDhQTrgBEHwIN4AF0bkBEMwIN4AFJTOhGEAQI9wAFpScwIRiAMGLcANYUFLcyTk39hM1OlFda3I1AOBbhBvAghKiwxQRdvJ/b7o3AIIN4QawIJvNpqQ457wbwg2A4EK4ASyKeTcAghXhBrAoZ+fmACumAAQZwg1gUXRuAAQrwg1gUc4VU8y5ARBsCDeARdG5ARCsCDeARbFaCkCwItwAFkXnBkCwItwAFuXcPLOsvFIOh2FyNQDgO4QbwKKax0XKZpNqHIYOHasyuxwA8BnCDWBR4aEhahoTIYlLUwCCC+EGsDDnpSkmFQMIJoQbwMKc4YbODYBgQrgBLOxU54YtGAAED8INYGHJ8SfvUkznBkAwIdwAFsacGwDBiHADWFgyc24ABCHCDWBhTCgGEIwIN4CF0bkBEIwIN4CFOTs35ZU1OlZVY3I1AOAbhBvAwuIiwxQdHiqJ7g2A4EG4ASzMZrOxYgpA0CHcABbnnHdTaifcAAgOhBvA4pITnJOKuUsxgOBAuAEsLimOy1IAggvhBrC45AS2YAAQXAg3gMXRuQEQbAg3gMUlJXAjPwDBhXADWBydGwDBhnADWJxztdShikrVOgyTqwGAxmd6uJk9e7ZatWqlqKgo9e3bV2vXrj3r+bNmzVL79u0VHR2tjIwMTZ48WSdOsMQVOJNmsZEKsUkOQzpYTvcGgPWZGm4WLFig7Oxs5ebmqqioSN27d9eAAQNUWlpa5/mvvfaaHnzwQeXm5mrr1q16/vnntWDBAv32t7/1ceVA4AgNsakZl6YABBFTw01+fr7uuusujR49Wp06ddLcuXMVExOjF154oc7zP/vsM/Xv31/Dhg1Tq1atdP3112vo0KFn7fZUVlbKbre7fQHBxjnvhknFAIKBaeGmqqpKhYWFysrKOlVMSIiysrK0atWqOp/Tr18/FRYWusLMl19+qUWLFumGG2444zh5eXlKTEx0fWVkZHj3jQABIJkVUwCCSJhZA5eVlam2tlYpKSlux1NSUrRt27Y6nzNs2DCVlZXpiiuukGEYqqmp0b333nvWy1JTp05Vdna263u73U7AQdA5tWKK+WkArM/0CcWeWL58uR5//HH96U9/UlFRkd566y0tXLhQjz766BmfExkZqYSEBLcvINg4OzfMuQEQDEzr3DRv3lyhoaEqKSlxO15SUqLU1NQ6n5OTk6MRI0bozjvvlCR17dpVFRUVuvvuu/XQQw8pJCSgshrgM8nxbMEAIHiYlgYiIiLUq1cvFRQUuI45HA4VFBQoMzOzzuccO3bstAATGhoqSTIM7t8BnElSPJ0bAMHDtM6NJGVnZ2vUqFHq3bu3+vTpo1mzZqmiokKjR4+WJI0cOVItWrRQXl6eJGngwIHKz89Xz5491bdvX+3atUs5OTkaOHCgK+QAOF1yPBOKAQQPU8PNkCFDdODAAU2bNk3FxcXq0aOHFi9e7JpkvG/fPrdOzcMPPyybzaaHH35Y33zzjZKSkjRw4ED9/ve/N+stAAHhVOfmhAzDkM1mM7kiAGg8NiPIrufY7XYlJibqyJEjTC5G0DhWVaNO05ZIkjY/cr3io8JNrggAPOPJ39/MwAWCQExEmOIiTzZqmXcDwOoIN0CQYN4NgGBBuAGCRHNWTAEIEoQbIEjQuQEQLAg3QJD44YopALAywg0QJLhLMYBgQbgBgkQSl6UABAnCDRAknHNuSu2EGwDWRrgBgoRzZ/AD5YQbANZGuAGCRFLcyXBzqKJKVTUOk6sBgMZDuAGCRJOYCIWFnNxT6mAF3RsA1kW4AYJESIhNzeOYdwPA+gg3QBBxzbthxRQACyPcAEHEOe+GLRgAWBnhBggidG4ABAPCDRBETnVu2IIBgHURboAgkpTAFgwArI9wAwQR5twACAaEGyCIMOcGQDAg3ABBxNm5OXC0UoZhmFwNADQOwg0QRJw7g1fVOnTkeLXJ1QBA4yDcAEEkKjxUidHhkrg0BcC6CDdAkHF2b5hUDMCqCDdAkEmOZ1IxAGsj3ABB5lTnhhv5AbAmwg0QZOjcALA6wg0QZJhzA8DqCDdAkEmOZwsGANZGuAGCDJ0bAFZHuAGCDHNuAFgd4QYIMs7OzZHj1TpRXWtyNQDgfYQbIMgkRocrIvTk//p0bwBYEeEGCDI2m83VvTlQTrgBYD2EGyAIuSYV2wk3AKyHcAMEoWQ6NwAsjHADBCHXZSk7WzAAsB7CDRCEXDfyo3MDwIIIN0AQYs4NACsj3ABBiDk3AKyMcAMEITo3AKyMcAMEoeSEk+GmrLxSDodhcjUA4F2EGyAINYs9GW5qHIa+O1ZlcjUA4F1hDXnSiRMn9J///EelpaVyOBxujw0aNMgrhQFoPBFhIWoaG6FDFVU6UF6pZnGRZpcEAF7jcbhZvHixRo4cqbKystMes9lsqq1lIz4gECTFRepQRZVK7ZXqkGp2NQDgPR5flrrvvvt0yy23aP/+/XI4HG5fBBsgcDjn3ZSyeSYAi/E43JSUlCg7O1spKSmNUQ8AH3HdpZhwA8BiPA43N998s5YvX94IpQDwJddy8KNswQDAWjyec/Pcc8/plltu0aeffqquXbsqPDzc7fEJEyZ4rTgAjce1BQOdGwAW43G4+fvf/64PP/xQUVFRWr58uWw2m+sxm81GuAECxKnODeEGgLV4HG4eeughTZ8+XQ8++KBCQrhNDhConFswlBFuAFiMx+mkqqpKQ4YMIdgAAY7ODQCr8jihjBo1SgsWLGiMWgD4kLNzU15Zo2NVNSZXAwDe4/FlqdraWs2cOVNLlixRt27dTptQnJ+f77XiADSeuMgwRYWH6ES1QweOVuqiZg26YTkA+B2PP802b96snj17SpK2bNni9tgPJxcD8G82m03J8VHad+jY9+Em1uySAMArPAo3tbW1mj59urp27aomTZo0Vk0AfCQpPlL7Dh1j3g0AS/Fozk1oaKiuv/56HT58uJHKAeBLydylGIAFeTyhuEuXLvryyy+9VsDs2bPVqlUrRUVFqW/fvlq7du1Zzz98+LDGjRuntLQ0RUZGql27dlq0aJHX6gGCCXcpBmBFHoebxx57TPfff7/ef/997d+/X3a73e3LEwsWLFB2drZyc3NVVFSk7t27a8CAASotLa3z/KqqKv3f//2f9u7dq3/+85/avn275s+frxYtWnj6NgDoVOem1E7nBoB1eDyh+IYbbpAkDRo0yG0CsWEYstlsHu0Mnp+fr7vuukujR4+WJM2dO1cLFy7UCy+8oAcffPC081944QUdOnRIn332mWuVVqtWrc46RmVlpSorT31wexrAACtzbcFQTrgBYB0eh5tly5Z5ZeCqqioVFhZq6tSprmMhISHKysrSqlWr6nzOv/71L2VmZmrcuHF69913lZSUpGHDhmnKlCkKDQ2t8zl5eXmaPn26V2oGrCaJzg0AC/I43Fx11VVeGbisrEy1tbVKSUlxO56SkqJt27bV+Zwvv/xS//73vzV8+HAtWrRIu3bt0q9+9StVV1crNze3zudMnTpV2dnZru/tdrsyMjK88h6AQOcMN3RuAFhJg+7adfjwYT3//PPaunWrJKlz58664447lJiY6NXi/pfD4VBycrLmzZun0NBQ9erVS998842efPLJM4abyMhIRUZGNmpdQKByzrk5WF6pWoeh0BDuVQUg8Hk8oXj9+vVq06aN/vCHP+jQoUM6dOiQ8vPz1aZNGxUVFdX7dZo3b67Q0FCVlJS4HS8pKVFqamqdz0lLS1O7du3cLkF17NhRxcXFqqqq8vStAEGvWVykQmySw5AOVtC9AWANHoebyZMna9CgQdq7d6/eeustvfXWW9qzZ49+9rOfadKkSfV+nYiICPXq1UsFBQWuYw6HQwUFBcrMzKzzOf3799euXbvkcDhcx3bs2KG0tDRFRER4+laAoBcaYlPTWObdALCWBnVupkyZorCwU1e0wsLC9Jvf/Ebr16/36LWys7M1f/58vfzyy9q6davGjh2riooK1+qpkSNHuk04Hjt2rA4dOqSJEydqx44dWrhwoR5//HGNGzfO07cB4HvJzLsBYDEez7lJSEjQvn371KFDB7fjX3/9teLj4z16rSFDhujAgQOaNm2aiouL1aNHDy1evNg1yXjfvn0KCTmVvzIyMrRkyRJNnjxZ3bp1U4sWLTRx4kRNmTLF07cB4HtJ8ZHSfukAnRsAFuFxuBkyZIjGjBmjp556Sv369ZMkrVy5Ug888ICGDh3qcQHjx4/X+PHj63xs+fLlpx3LzMzU6tWrPR4HQN3o3ACwGo/DzVNPPSWbzaaRI0eqpqZGkhQeHq6xY8dqxowZXi8QQOM6da8btmAAYA0eh5uIiAj98Y9/VF5ennbv3i1JatOmjWJiYrxeHIDGR+cGgNU06D43khQTE6OuXbt6sxYAJkj6fgsGVksBsAqPw01FRYVmzJihgoIClZaWui3LluTVHcMBNL7kBOfO4IQbANbgcbi588479fHHH2vEiBFKS0tz2zwTQOBxXZY6WunaABcAApnH4eaDDz7QwoUL1b9//8aoB4CPOScUH6+uVXlljeKjwk2uCADOj8c38WvSpImaNm3aGLUAMEFMRJjiIk/+O+cAl6YAWIDH4ebRRx/VtGnTdOzYscaoB4AJXMvBCTcALMDjy1JPP/20du/erZSUFLVq1Urh4e4tbE82zwTgH5LiI7WnrILODQBL8DjcDB48uBHKAGAmOjcArMTjcJObm9sYdQAw0Q9XTAFAoPN4zg0A6znVuWELBgCBj3ADQMnf36WYzg0AKyDcAHB1bgg3AKyAcAOAOTcALIVwA8DVuTlYUaXqWsc5zgYA/1av1VLZ2dn1fsH8/PwGFwPAHE1jIhQWYlONw1BZeaXSEqPNLgkAGqxe4WbDhg1u3xcVFammpkbt27eXJO3YsUOhoaHq1auX9ysE0OhCQmxqHhepYvsJHThKuAEQ2OoVbpYtW+b67/z8fMXHx+vll19WkyZNJEnfffedRo8erSuvvLJxqgTQ6JLiT4abUjvzbgAENo/n3Dz99NPKy8tzBRvp5Gaajz32mJ5++mmvFgfAd1yTissJNwACm8fhxm6368CBA6cdP3DggI4ePeqVogD4nutGfnRuAAQ4j8PNz3/+c40ePVpvvfWW/vvf/+q///2v3nzzTY0ZM0Y33XRTY9QIwAdOdW64SzGAwObx3lJz587V/fffr2HDhqm6uvrki4SFacyYMXryySe9XiAA36BzA8AqPA43MTEx+tOf/qQnn3xSu3fvliS1adNGsbGxXi8OgO8kObdgYM4NgADX4Jv47d+/X/v371fbtm0VGxsrwzC8WRcAH6NzA8AqPA43Bw8e1HXXXad27drphhtu0P79+yVJY8aM0a9//WuvFwjAN364Wop/rAAIZB6Hm8mTJys8PFz79u1TTEyM6/iQIUO0ePFirxYHwHecnZuqGofsx2tMrgYAGs7jOTcffvihlixZogsvvNDteNu2bfXVV195rTAAvhUVHqqEqDDZT9ToQPkJJcaEm10SADSIx52biooKt46N06FDhxQZGemVogCYg3k3AKzA43Bz5ZVX6pVXXnF9b7PZ5HA4NHPmTF1zzTVeLQ6AbyV/v2Kq9CjhBkDg8viy1MyZM3Xddddp/fr1qqqq0m9+8xt9/vnnOnTokFauXNkYNQLwkeSE7ycVE24ABDCPOzddunTRjh07dMUVV+jGG29URUWFbrrpJm3YsEFt2rRpjBoB+EhS3PeXpY5yl2IAgcvjzs2+ffuUkZGhhx56qM7HWrZs6ZXCAPgenRsAVuBx56Z169Z1bpx58OBBtW7d2itFATCHa0Ix4QZAAPM43BiGIZvNdtrx8vJyRUVFeaUoAOZwTiimcwMgkNX7slR2drakk6ujcnJy3JaD19bWas2aNerRo4fXCwTgO3RuAFhBvcPNhg0bJJ3s3GzevFkRERGuxyIiItS9e3fdf//93q8QgM84t2A4crxalTW1igwLNbkiAPBcvcPNsmXLJEmjR4/WH//4RyUkJDRaUQDMkRgdrojQEFXVOnTgaKUubHL6DTsBwN95POdm1qxZqqk5fd+ZQ4cOyW63e6UoAOaw2WyuS1PMuwEQqDwON7feeqtef/31046/8cYbuvXWW71SFADzNGfeDYAA53G4WbNmTZ3bLFx99dVas2aNV4oCYJ5kwg2AAOdxuKmsrKzzslR1dbWOHz/ulaIAmIfLUgACncfhpk+fPpo3b95px+fOnatevXp5pSgA5kl2hRu2YAAQmDzefuGxxx5TVlaWNm3apOuuu06SVFBQoHXr1unDDz/0eoEAfIsb+QEIdB53bvr3769Vq1bpwgsv1BtvvKH33ntPl1xyif7zn//oyiuvbIwaAfgQN/IDEOg87txIUo8ePfTaa695uxYAfiCZOTcAApzHnRtJ2r17tx5++GENGzZMpaWlkqQPPvhAn3/+uVeLA+B7P5xQ7HAYJlcDAJ7zONx8/PHH6tq1q9asWaM333xT5eXlkqRNmzYpNzfX6wUC8K3mcSfDTY3D0OHj1SZXAwCe8zjcPPjgg3rssce0dOlSt/2lrr32Wq1evdqrxQHwvYiwEDWJCZcklbJiCkAA8jjcbN68WT//+c9PO56cnKyysjKvFAXAXKyYAhDIPA43F1xwgfbv33/a8Q0bNqhFixZeKQqAuVwrpuyEGwCBp0F7S02ZMkXFxcWy2WxyOBxauXKl7r//fo0cObIxagTgY64VU+WEGwCBx+Nw8/jjj6tDhw7KyMhQeXm5OnXqpB/96Efq16+fHn744caoEYCP0bkBEMg8vs9NRESE5s+fr5ycHG3ZskXl5eXq2bOn2rZt2xj1ATDBqRv5MaEYQOBp0E38JKlly5Zq2bKlN2sB4CfYPBNAIKtXuMnOztajjz6q2NhYZWdnn/XcuLg4de7cWTfffLNCQ0O9UiQA32K1FIBAVq85Nxs2bFB1dbXrv8/29dFHH+nuu+/WHXfcUe8iZs+erVatWikqKkp9+/bV2rVr6/W8119/XTabTYMHD673WADOLTmBzg2AwFWvzs2yZcvq/O8zWb9+vWvH8HNZsGCBsrOzNXfuXPXt21ezZs3SgAEDtH37diUnJ5/xeXv37tX999/PZp1AI3BeljpaWaPjVbWKjqALCyBwNGhvKSfDMGQYp+89061bN73yyiv1eo38/HzdddddGj16tDp16qS5c+cqJiZGL7zwwhmfU1tbq+HDh2v69Om6+OKLG1w/gLrFR4YpKvzkxwPdGwCBpkHh5pVXXlHXrl0VHR2t6OhodevWTX/9619dj0dEROjGG2885+tUVVWpsLBQWVlZpwoKCVFWVpZWrVp1xuf97ne/U3JyssaMGXPOMSorK2W3292+AJydzWZjxRSAgOXxaqn8/Hzl5ORo/Pjx6t+/vyRpxYoVuvfee1VWVqbJkyfX+7XKyspUW1urlJQUt+MpKSnatm1bnc9ZsWKFnn/+eW3cuLFeY+Tl5Wn69On1rgnAScnxUfr60HE6NwACjsfh5tlnn9WcOXPc7kY8aNAgde7cWY888ohH4cZTR48e1YgRIzR//nw1b968Xs+ZOnWq2wovu92ujIyMxioRsIykOGfnhnADILB4HG7279+vfv36nXa8X79+de45dTbNmzdXaGioSkpK3I6XlJQoNTX1tPN3796tvXv3auDAga5jDodDkhQWFqbt27erTZs2bs+JjIxUZGSkR3UBYMUUgMDl8ZybSy65RG+88cZpxxcsWODxXYojIiLUq1cvFRQUuI45HA4VFBQoMzPztPM7dOigzZs3a+PGja6vQYMG6ZprrtHGjRvpyABedKpzw5wbAIHF487N9OnTNWTIEH3yySeuOTcrV65UQUFBnaHnXLKzszVq1Cj17t1bffr00axZs1RRUaHRo0dLkkaOHKkWLVooLy9PUVFR6tKli9vzL7jgAkk67TiA80PnBkCg8jjc/OIXv9DatWuVn5+vd955R5LUsWNHrV27Vj179vS4gCFDhujAgQOaNm2aiouL1aNHDy1evNg1yXjfvn0KCTmvFesAGuDUainCDYDAYjPqulHNGVRXV+uee+5RTk6OWrdu3Zh1NRq73a7ExEQdOXJECQkJZpcD+K0t3xzRz55doaT4SK17KOvcTwCARuTJ398etUTCw8P15ptvnldxAAKDs3NzsLxStY56/xsIAEzn8fWewYMHuy5HAbCuZrERstkkhyEdrODSFIDA4fGcm7Zt2+p3v/udVq5cqV69eik2Ntbt8QkTJnitOADmCQsNUbPYSJWVV+rA0UrXTuEA4O88DjfPP/+8LrjgAhUWFqqwsNDtMZvNRrgBLCQp/mS4KT1aqc5mFwMA9eRxuNmzZ09j1AHADyXHR2rrfpaDAwgsjbIrOABrcE4qJtwACCQNCjfPP/+8unTpoqioKNeN9f7yl794uzYAJksm3AAIQB5flpo2bZry8/N13333ubZIWLVqlSZPnqx9+/bpd7/7ndeLBGCOUzfyYwsGAIHD43AzZ84czZ8/X0OHDnUdGzRokLp166b77ruPcANYiHOFFJ0bAIHE48tS1dXV6t2792nHe/XqpZqaGq8UBcA/sAUDgEDkcbgZMWKE5syZc9rxefPmafjw4V4pCoB/YM4NgEDk8WUp6eSE4g8//FCXX365JGnNmjXat2+fRo4cqezsbNd5+fn53qkSgCmcnZtjVbUqr6xRXGSDPjIAwKc8/qTasmWLLr30UknS7t27JUnNmzdX8+bNtWXLFtd5NpvNSyUCMEtsZJhiI0JVUVWrUvsJxSXFmV0SAJyTx+Fm2bJljVEHAD+VFB+pioPHdOBopS4m3AAIAOd1Ez8A1udcMcWkYgCBgnAD4KySEphUDCCwEG4AnFVSHMvBAQQWwg2As0qmcwMgwBBuAJzVqc4NWzAACAyEGwBnlZzAFgwAAgvhBsBZOTs3hBsAgYJwA+CsnHNuDh2rUnWtw+RqAODcCDcAzqppTIRCQ2wyDOlgeZXZ5QDAORFuAJxVSIhNzeMiJHFpCkBgINwAOCfnBpqsmAIQCAg3AM6JLRgABBLCDYBzYsUUgEBCuAFwTs4VU1yWAhAICDcAzik5ns4NgMBBuAFwTqcmFBNuAPg/wg2Ac0qKZwsGAIGDcAPgnJJ/0LkxDMPkagDg7Ag3AM7JeVmqqsYh+4kak6sBgLMj3AA4p6jwUMVHhUmSDrBiCoCfI9wAqJdkJhUDCBCEGwD1ksRycAABgnADoF6SWTEFIEAQbgDUC/e6ARAoCDcA6sU158bOhGIA/o1wA6BeXHNuyuncAPBvhBsA9eKcc1NqJ9wA8G+EGwD14twZnM4NAH9HuAFQL0lxJ8PN4WPVqqypNbkaADgzwg2AerkgJlzhoTZJUll5lcnVAMCZEW4A1IvNZnN1b1gxBcCfEW4A1FtSAjfyA+D/CDcA6s3VuSHcAPBjhBsA9eZaMUW4AeDHCDcA6o3ODYBAQLgBUG90bgAEAsINgHpzdm4OHGW1FAD/RbgBUG/J36+W4rIUAH9GuAFQb86dwcvKK+VwGCZXAwB1I9wAqLfm31+Wqq41dPh4tcnVAEDdCDcA6i0iLERNYsIlMakYgP8i3ADwSFK8czk4k4oB+Ce/CDezZ89Wq1atFBUVpb59+2rt2rVnPHf+/Pm68sor1aRJEzVp0kRZWVlnPR+AdyXHswUDAP9merhZsGCBsrOzlZubq6KiInXv3l0DBgxQaWlpnecvX75cQ4cO1bJly7Rq1SplZGTo+uuv1zfffOPjyoHgdKpzQ7gB4J9MDzf5+fm66667NHr0aHXq1Elz585VTEyMXnjhhTrPf/XVV/WrX/1KPXr0UIcOHfSXv/xFDodDBQUFPq4cCE7OFVN0bgD4K1PDTVVVlQoLC5WVleU6FhISoqysLK1ataper3Hs2DFVV1eradOmdT5eWVkpu93u9gWg4ejcAPB3poabsrIy1dbWKiUlxe14SkqKiouL6/UaU6ZMUXp6ultA+qG8vDwlJia6vjIyMs67biCYJcVzl2IA/s30y1LnY8aMGXr99df19ttvKyoqqs5zpk6dqiNHjri+vv76ax9XCVgLnRsA/i7MzMGbN2+u0NBQlZSUuB0vKSlRamrqWZ/71FNPacaMGfroo4/UrVu3M54XGRmpyMhIr9QLgNVSAPyfqZ2biIgI9erVy20ysHNycGZm5hmfN3PmTD366KNavHixevfu7YtSAXzP2bk5eqJGJ6prTa4GAE5n+mWp7OxszZ8/Xy+//LK2bt2qsWPHqqKiQqNHj5YkjRw5UlOnTnWd/8QTTygnJ0cvvPCCWrVqpeLiYhUXF6u8vNystwAElYSoMEWGnfzoKLXTvQHgf0y9LCVJQ4YM0YEDBzRt2jQVFxerR48eWrx4sWuS8b59+xQSciqDzZkzR1VVVbr55pvdXic3N1ePPPKIL0sHgpLNZlNyQqS+PnRcB8pPqGWzGLNLAgA3pocbSRo/frzGjx9f52PLly93+37v3r2NXxCAs0qKOxlu6NwA8EemX5YCEHhck4rLCTcA/A/hBoDHXMvB6dwA8EOEGwAeYwsGAP6McAPAY6du5MddigH4H8INAI8lJ3zfuWHODQA/RLgB4LGkuJMTiplzA8AfEW4AeMzZuTlYUaVah2FyNQDgjnADwGPNYiNks0m1DkOHKqrMLgcA3BBuAHgsLDREzWIjJDGpGID/IdwAaJDmcSwHB+CfCDcAGiQ54ftJxYQbAH6GcAOgQbiRHwB/RbgB0CBJhBsAfopwA6BB6NwA8FeEGwANwhYMAPwV4QZAgyTHn5xQTOcGgL8h3ABokFOdG8INAP9CuAHQIM45N8eqalVRWWNyNQBwCuEGQIPERoYpJiJUEt0bAP6FcAOgwVgxBcAfEW4ANBgrpgD4I8INgAZzrpgqtdO5AeA/CDcAGsx1l+Jywg0A/0G4AdBgrstSdG4A+BHCDYAGS6ZzA8APEW4ANNipzg0TigH4D8INgAZzTiguo3MDwI8QbgA0mLNzc7CiSjW1DpOrAYCTCDcAGqxpbIRCQ2wyjJMBBwD8AeEGQIOFhtjULDZCEiumAPgPwg2A85Kc4FwxxaRiAP6BcAPgvCTFca8bAP6FcAPgvDhXTLF5JgB/QbgBcF5ObZ5JuAHgHwg3AM6Lc84NO4MD8BeEGwDnxTnnhstSAPwF4QbAeTnVuSHcAPAPhBsA5+WHE4oNwzC5GgAg3AA4T84JxZU1DtlP1JhcDQAQbgCcp6jwUMVHhUli3g0A/0C4AXDeTi0HZ8UUAPMRbgCct+R4VkwB8B+EGwDnLYm7FAPwI4QbAOeNzg0Af0K4AXDe2IIBgD8h3AA4b3RuAPgTwg2A88ZqKQD+hHAD4Lw571LMZSkA/oBwA+C8OTs3h49Vq7Km1uRqAAQ7wg2A89YkJlzhoTZJUll5lcnVAAh2hBsA581msykpjknFAPwD4QaAV7gmFduZVAzAXIQbAF7huktxOZ0bAOYi3ADwilOdG8INAHMRbgB4hetGfnRuAJgszOwCAFhD87gISdJ//ntYq3YfVJ/WTRUaYvPZ+LUOQ2v3HFLp0RNKjo8KuvH9oQbG53fA7PGd/CLczJ49W08++aSKi4vVvXt3Pfvss+rTp88Zz//HP/6hnJwc7d27V23bttUTTzyhG264wYcVA/ihxVv2K3/pDknSlm/sGjp/tdISo5Q7sJN+3CXNJ+NPf+8L7T9yajJzMI3vDzUwPr8DZo//Q6ZfllqwYIGys7OVm5uroqIide/eXQMGDFBpaWmd53/22WcaOnSoxowZow0bNmjw4MEaPHiwtmzZ4uPKAUgnP9DG/q1I3x2rdjtefOSExv6tSIu37PfJ+D/8QA2m8f2hBsbnd8Ds8f+XzTAMw6cj/o++ffvqsssu03PPPSdJcjgcysjI0H333acHH3zwtPOHDBmiiooKvf/++65jl19+uXr06KG5c+eeczy73a7ExEQdOXJECQkJ3nsjQBCqdRi64ol/n/aB9kNJ8ZF6895+CmmEf0o5HNIv5n521nvrWHl8f6iB8fkdONf4NkmpiVFaMeXa87pE5cnf36aGm6qqKsXExOif//ynBg8e7Do+atQoHT58WO++++5pz2nZsqWys7M1adIk17Hc3Fy988472rRp02nnV1ZWqrLy1A/cbrcrIyODcAN4wardBzV0/mqzywAQAP5+1+XKbNOswc/3JNyYOuemrKxMtbW1SklJcTuekpKibdu21fmc4uLiOs8vLi6u8/y8vDxNnz7dOwUDcFPfXcDDQ20KsXl/UqHDMFRde+5/n1l1fH+ogfH5Hajv+PX9vPAGv5hQ3JimTp2q7Oxs1/fOzg2A8+fcDfxcXrmj73n9i+1M6ts5sur4/lAD4/M7UN/x6/t54Q2mTihu3ry5QkNDVVJS4na8pKREqampdT4nNTXVo/MjIyOVkJDg9gXAO/q0bqq0xCid6d+CNp1cLdGndVPGbyRm18D4/A6YPX5dTA03ERER6tWrlwoKClzHHA6HCgoKlJmZWedzMjMz3c6XpKVLl57xfACNJzTEptyBnSTptA825/e5Azs12n0ugn18f6iB8fkdMHv8upi+FDw7O1vz58/Xyy+/rK1bt2rs2LGqqKjQ6NGjJUkjR47U1KlTXedPnDhRixcv1tNPP61t27bpkUce0fr16zV+/Hiz3gIQ1H7cJU1zbrtUqYnuLefUxCjNue3SRr+/RbCP7w81MD6/A2aP/79MXwouSc8995zrJn49evTQM888o759+0qSrr76arVq1UovvfSS6/x//OMfevjhh1038Zs5c2a9b+LHUnCgcZh9Z9JgH98famB8fgcac/yAWQpuBsINAACBx5O/v02/LAUAAOBNhBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGApYWYX4GvOGzLb7XaTKwEAAPXl/Hu7PhsrBF24OXr0qCQpIyPD5EoAAICnjh49qsTExLOeE3R7SzkcDn377beKj4+XzebdzcTsdrsyMjL09ddfB+W+VcH+/iV+Brz/4H7/Ej+DYH//UuP9DAzD0NGjR5Wenq6QkLPPqgm6zk1ISIguvPDCRh0jISEhaH+pJd6/xM+A9x/c71/iZxDs719qnJ/BuTo2TkwoBgAAlkK4AQAAlkK48aLIyEjl5uYqMjLS7FJMEezvX+JnwPsP7vcv8TMI9vcv+cfPIOgmFAMAAGujcwMAACyFcAMAACyFcAMAACyFcAMAACyFcOMls2fPVqtWrRQVFaW+fftq7dq1ZpfkM3l5ebrssssUHx+v5ORkDR48WNu3bze7LNPMmDFDNptNkyZNMrsUn/rmm2902223qVmzZoqOjlbXrl21fv16s8vyidraWuXk5Kh169aKjo5WmzZt9Oijj9ZrD5xA9cknn2jgwIFKT0+XzWbTO++84/a4YRiaNm2a0tLSFB0draysLO3cudOcYhvB2d5/dXW1pkyZoq5duyo2Nlbp6ekaOXKkvv32W/MK9rJz/fn/0L333iubzaZZs2b5rD7CjRcsWLBA2dnZys3NVVFRkbp3764BAwaotLTU7NJ84uOPP9a4ceO0evVqLV26VNXV1br++utVUVFhdmk+t27dOv35z39Wt27dzC7Fp7777jv1799f4eHh+uCDD/TFF1/o6aefVpMmTcwuzSeeeOIJzZkzR88995y2bt2qJ554QjNnztSzzz5rdmmNpqKiQt27d9fs2bPrfHzmzJl65plnNHfuXK1Zs0axsbEaMGCATpw44eNKG8fZ3v+xY8dUVFSknJwcFRUV6a233tL27ds1aNAgEyptHOf683d6++23tXr1aqWnp/uosu8ZOG99+vQxxo0b5/q+trbWSE9PN/Ly8kysyjylpaWGJOPjjz82uxSfOnr0qNG2bVtj6dKlxlVXXWVMnDjR7JJ8ZsqUKcYVV1xhdhmm+elPf2rccccdbsduuukmY/jw4SZV5FuSjLffftv1vcPhMFJTU40nn3zSdezw4cNGZGSk8fe//92EChvX/77/uqxdu9aQZHz11Ve+KcqHzvT+//vf/xotWrQwtmzZYlx00UXGH/7wB5/VROfmPFVVVamwsFBZWVmuYyEhIcrKytKqVatMrMw8R44ckSQ1bdrU5Ep8a9y4cfrpT3/q9rsQLP71r3+pd+/euuWWW5ScnKyePXtq/vz5ZpflM/369VNBQYF27NghSdq0aZNWrFihn/zkJyZXZo49e/aouLjY7f+FxMRE9e3bN6g/F202my644AKzS/EJh8OhESNG6IEHHlDnzp19Pn7QbZzpbWVlZaqtrVVKSorb8ZSUFG3bts2kqszjcDg0adIk9e/fX126dDG7HJ95/fXXVVRUpHXr1pldiim+/PJLzZkzR9nZ2frtb3+rdevWacKECYqIiNCoUaPMLq/RPfjgg7Lb7erQoYNCQ0NVW1ur3//+9xo+fLjZpZmiuLhYkur8XHQ+FkxOnDihKVOmaOjQoUGzmeYTTzyhsLAwTZgwwZTxCTfwqnHjxmnLli1asWKF2aX4zNdff62JEydq6dKlioqKMrscUzgcDvXu3VuPP/64JKlnz57asmWL5s6dGxTh5o033tCrr76q1157TZ07d9bGjRs1adIkpaenB8X7x5lVV1frl7/8pQzD0Jw5c8wuxycKCwv1xz/+UUVFRbLZbKbUwGWp89S8eXOFhoaqpKTE7XhJSYlSU1NNqsoc48eP1/vvv69ly5bpwgsvNLscnyksLFRpaakuvfRShYWFKSwsTB9//LGeeeYZhYWFqba21uwSG11aWpo6derkdqxjx47at2+fSRX51gMPPKAHH3xQt956q7p27aoRI0Zo8uTJysvLM7s0Uzg/+4L9c9EZbL766istXbo0aLo2n376qUpLS9WyZUvXZ+JXX32lX//612rVqpVPaiDcnKeIiAj16tVLBQUFrmMOh0MFBQXKzMw0sTLfMQxD48eP19tvv61///vfat26tdkl+dR1112nzZs3a+PGja6v3r17a/jw4dq4caNCQ0PNLrHR9e/f/7Tl/zt27NBFF11kUkW+dezYMYWEuH+choaGyuFwmFSRuVq3bq3U1FS3z0W73a41a9YEzeeiM9js3LlTH330kZo1a2Z2ST4zYsQI/ec//3H7TExPT9cDDzygJUuW+KQGLkt5QXZ2tkaNGqXevXurT58+mjVrlioqKjR69GizS/OJcePG6bXXXtO7776r+Ph41zX1xMRERUdHm1xd44uPjz9tflFsbKyaNWsWNPOOJk+erH79+unxxx/XL3/5S61du1bz5s3TvHnzzC7NJwYOHKjf//73atmypTp37qwNGzYoPz9fd9xxh9mlNZry8nLt2rXL9f2ePXu0ceNGNW3aVC1bttSkSZP02GOPqW3btmrdurVycnKUnp6uwYMHm1e0F53t/aelpenmm29WUVGR3n//fdXW1ro+F5s2baqIiAizyvaac/35/2+YCw8PV2pqqtq3b++bAn22Lsvinn32WaNly5ZGRESE0adPH2P16tVml+Qzkur8evHFF80uzTTBthTcMAzjvffeM7p06WJERkYaHTp0MObNm2d2ST5jt9uNiRMnGi1btjSioqKMiy++2HjooYeMyspKs0trNMuWLavz//tRo0YZhnFyOXhOTo6RkpJiREZGGtddd52xfft2c4v2orO9/z179pzxc3HZsmVml+4V5/rz/1++XgpuMwwL30ITAAAEHebcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAEGgVatWmjVrltll+L3ly5fLZrPp8OHDDX6NRx55RD169DjvWvgzAxqOcAMAXnT//fe7bRgJwPfYOBMAvCguLk5xcXFmlwEENTo3QIC7+uqrNX78eI0fP16JiYlq3ry5cnJy9L/bxh07dkx33HGH4uPj1bJly9N27N68ebOuvfZaRUdHq1mzZrr77rtVXl7uevz222/X4MGD9dRTTyktLU3NmjXTuHHjVF1d7Trnu+++08iRI9WkSRPFxMToJz/5iXbu3Ol6/KWXXtIFF1ygJUuWqGPHjoqLi9OPf/xj7d+//4zvz3mpaMmSJerZs6eio6N17bXXqrS0VB988IE6duyohIQEDRs2TMeOHXM9r7KyUhMmTFBycrKioqJ0xRVXaN26dW6vvWjRIrVr107R0dG65pprtHfv3tPGX7Fiha688kpFR0crIyNDEyZMUEVFxRnr/d/LUvX5uZWWlmrgwIGKjo5W69at9eqrr572uocPH9add96ppKQkJSQk6Nprr9WmTZskSQcOHFBqaqoef/xx1/mfffaZIiIi6CIhOPlsi04AjeKqq64y4uLijIkTJxrbtm0z/va3vxkxMTFuu3JfdNFFRtOmTY3Zs2cbO3fuNPLy8oyQkBBj27ZthmEYRnl5uZGWlmbcdNNNxubNm42CggKjdevWbjv8jho1ykhISDDuvfdeY+vWrcZ777132jiDBg0yOnbsaHzyySfGxo0bjQEDBhiXXHKJUVVVZRiGYbz44otGeHi4kZWVZaxbt84oLCw0OnbsaAwbNuyM78+5+/Dll19urFixwigqKjIuueQS46qrrjKuv/56o6ioyPjkk0+MZs2aGTNmzHA9b8KECUZ6erqxaNEi4/PPPzdGjRplNGnSxDh48KBhGIaxb98+IzIy0sjOznb93FJSUgxJxnfffWcYhmHs2rXLiI2NNf7whz8YO3bsMFauXGn07NnTuP32289Yb25urtG9e3ePfm4/+clPjO7duxurVq0y1q9fb/Tr18+Ijo5220U5KyvLGDhwoLFu3Tpjx44dxq9//WujWbNmrvezcOFCIzw83Fi3bp1ht9uNiy++2Jg8efIZ6wSsjHADBLirrrrK6Nixo+FwOFzHpkyZYnTs2NH1/UUXXWTcdtttru8dDoeRnJxszJkzxzAMw5g3b57RpEkTo7y83HXOwoULjZCQEKO4uNgwjJN/SV900UVGTU2N65xbbrnFGDJkiGEYhrFjxw5DkrFy5UrX42VlZUZ0dLTxxhtvGIZxMtxIMnbt2uU6Z/bs2UZKSsoZ358z3Hz00UeuY3l5eYYkY/fu3a5j99xzjzFgwADDME6GtfDwcOPVV191PV5VVWWkp6cbM2fONAzDMKZOnWp06tTJbawpU6a4hZsxY8YYd999t9s5n376qRESEmIcP368znrrCjdn+7lt377dkGSsXbvW9fjWrVsNSa5w8+mnnxoJCQnGiRMn3MZq06aN8ec//9n1/a9+9SujXbt2xrBhw4yuXbuedj4QLLgsBVjA5ZdfLpvN5vo+MzNTO3fuVG1tretYt27dXP9ts9mUmpqq0tJSSdLWrVvVvXt3xcbGus7p37+/HA6Htm/f7jrWuXNnhYaGur5PS0tze42wsDD17dvX9XizZs3Uvn17bd261XUsJiZGbdq0qfM1zuaH9aekpCgmJkYXX3yx2zHn6+zevVvV1dXq37+/6/Hw8HD16dPHVcvWrVvdapVO/tx+aNOmTXrppZdc82ji4uI0YMAAORwO7dmz55w1O9Xn59arVy/X4x06dNAFF1zgVkd5ebmaNWvmVsuePXu0e/du13lPPfWUampq9I9//EOvvvqqIiMj610jYCVMKAaCRHh4uNv3NptNDofDL17D+J/5Qed6ns1m80ot51JeXq577rlHEyZMOO2xli1b1vt1zrfW8vJypaWlafny5ac99sMQtHv3bn377bdyOBzau3evunbtWu8xACsh3AAWsGbNGrfvV69erbZt27p1C86mY8eOeumll1RRUeHq3qxcuVIhISFq3759vV+jpqZGa9asUb9+/SRJBw8e1Pbt29WpUycP3s35a9OmjSIiIrRy5UpddNFFkqTq6mqtW7dOkyZNctX7r3/9y+15q1evdvv+0ksv1RdffKFLLrmk0Wrt0KGDampqVFhYqMsuu0yStH37drd77Vx66aUqLi5WWFiYWrVqVefrVFVV6bbbbtOQIUPUvn173Xnnndq8ebOSk5MbrXbAX3FZCrCAffv2KTs7W9u3b9ff//53Pfvss5o4cWK9nz98+HBFRUVp1KhR2rJli5YtW6b77rtPI0aMUEpKSr1eo23btrrxxht11113acWKFdq0aZNuu+02tWjRQjfeeGND31qDxMbGauzYsXrggQe0ePFiffHFF7rrrrt07NgxjRkzRpJ07733aufOnXrggQe0fft2vfbaa3rppZfcXmfKlCn67LPPNH78eG3cuFE7d+7Uu+++q/Hjx3ut1vbt2+vHP/6x7rnnHq1Zs0aFhYW68847FR0d7TonKytLmZmZGjx4sD788EPt3btXn332mR566CGtX79ekvTQQw/pyJEjeuaZZzRlyhS1a9dOd9xxh9fqBAIJ4QawgJEjR+r48ePq06ePxo0bp4kTJ+ruu++u9/NjYmK0ZMkSHTp0SJdddpluvvlmXXfddXruuec8quPFF19Ur1699LOf/UyZmZkyDEOLFi067bKML8yYMUO/+MUvNGLECF166aXatWuXlixZoiZNmkg6eVnpzTff1DvvvKPu3btr7ty5bkuppZPzfD7++GPt2LFDV155pXr27Klp06YpPT3dq7W++OKLSk9P11VXXaWbbrpJd999t1vHxWazadGiRfrRj36k0aNHq127drr11lv11VdfKSUlRcuXL9esWbP017/+VQkJCQoJCdFf//pXffrpp5ozZ45XawUCgc2oz8VuAH7r6quvVo8ePbhVPwB8j84NAACwFMINAACwFC5LAQAAS6FzAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALOX/AcWx2gLZHLs7AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(gamma_proj_norm)\n",
    "plt.scatter(x=np.arange(len(gamma_proj_norm)),y=gamma_proj_norm)\n",
    "plt.xlabel('phonon mode index')\n",
    "plt.ylabel('projected norm')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51ca6bab",
   "metadata": {},
   "source": [
    "With the external library to perform symmetry analysis, we can determine the vibrational mode symmetry irrep as below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "96252324",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Symmetry irrep:  T2\n"
     ]
    }
   ],
   "source": [
    "# In this example, they are translational modes, which transform as (x,y,z)\n",
    "test_modes = np.array(gamma_vecs.T[:3].real)\n",
    "test_modes = test_modes.reshape(3,5,3)\n",
    "\n",
    "sym_modes_gs = SymmetryNormalModes(group='Td', coordinates=coordinates, modes=test_modes, symbols=symbols)\n",
    "print('Symmetry irrep: ', sym_modes_gs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "f9ca68c5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Symmetry irrep:  T1\n"
     ]
    }
   ],
   "source": [
    "# In this example, they are rotational modes, which transform as (Rx,Ry,Rz)\n",
    "test_modes = np.array(gamma_vecs.T[3:6].real)\n",
    "test_modes = test_modes.reshape(3,5,3)\n",
    "\n",
    "sym_modes_gs = SymmetryNormalModes(group='Td', coordinates=coordinates, modes=test_modes, symbols=symbols)\n",
    "print('Symmetry irrep: ', sym_modes_gs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "608531ba",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Symmetry irrep:  E\n"
     ]
    }
   ],
   "source": [
    "# one can check any (degenerate) mode(s) to find their symmetry irrep\n",
    "test_modes = np.array(gamma_vecs.T[6:8].real)\n",
    "test_modes = test_modes.reshape(2,5,3)\n",
    "\n",
    "sym_modes_gs = SymmetryNormalModes(group='Td', coordinates=coordinates, modes=test_modes, symbols=symbols)\n",
    "print('Symmetry irrep: ', sym_modes_gs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0de9034d",
   "metadata": {},
   "source": [
    "## IR / Raman activity\n",
    "\n",
    "The point group symmetry imposes constraints of the IR / Raman active modes.\n",
    "We fwllow the symmetry arguments in the group theory book relevant for the IR / Raman transitions with phonons.\n",
    "In short summary, the IR modes would transform as vectors, while the Raman modes would transform as symmetric second rank tensors.\n",
    "In the following, we will implement these symmetry selection rules:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6b302186",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Molecular point group:  Td\n",
      "Translation mode irreps:  T2\n",
      "Rotational mode irreps:  T1\n"
     ]
    }
   ],
   "source": [
    "from posym import SymmetryNormalModes\n",
    "\n",
    "from pointgroup import PointGroup\n",
    "\n",
    "pos_CM = mol_atoms.get_center_of_mass()\n",
    "\n",
    "coordinates = mol_atoms.get_positions() - pos_CM\n",
    "\n",
    "symbols = mol_atoms.get_chemical_symbols()\n",
    "\n",
    "pg = PointGroup(positions=coordinates, \n",
    "                symbols=symbols)\n",
    "\n",
    "\n",
    "print('Molecular point group: ',pg.get_point_group())\n",
    "\n",
    "tran_modes_all = [tr_mode1,tr_mode2,tr_mode3]\n",
    "sym_modes_trans = SymmetryNormalModes(group=pg.get_point_group(), coordinates=coordinates, modes=tran_modes_all, symbols=symbols)\n",
    "\n",
    "\n",
    "print('Translation mode irreps: ',sym_modes_trans)\n",
    "\n",
    "\n",
    "rot_modes_all = np.array([rot_mode1,rot_mode2,rot_mode3])\n",
    "\n",
    "sym_modes_gs_rot = SymmetryNormalModes(group=pg.get_point_group(), coordinates=coordinates, modes=rot_modes_all, symbols=symbols)\n",
    "print('Rotational mode irreps: ',sym_modes_gs_rot)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0fc5a4d",
   "metadata": {},
   "source": [
    "Because of the symmetry selection rules, the active modes:\n",
    "\n",
    "IR active modes: They will transform as vector irreps found in the translational modes at zero energy.\n",
    "\n",
    "Raman active modes: Since they transform as symmetric second rank tensors, we derive the symmetry irreps by finding the irreps for the second rank tensor first (square of the vector irreps), then removing the anti-symmetric component (derived in the rotational modes)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "9f31178d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "point group:  Td\n",
      "IR active modes:  T2\n",
      "Raman active modes:  A1 + E + T2\n"
     ]
    }
   ],
   "source": [
    "set_pg = pg.get_point_group()\n",
    "\n",
    "\n",
    "print('point group: ',set_pg)\n",
    "\n",
    "def decomp_irreps(sym_irreps):\n",
    "    ir_rep = np.dot(sym_irreps._pg.trans_matrix_inv, sym_irreps.get_reduced_op_representation().values)\n",
    "    index=sym_irreps._pg.ir_labels\n",
    "    in_irreps = []\n",
    "    for ir_, idx_label in zip(ir_rep,index):\n",
    "        if ir_ > 0.5:\n",
    "            in_irreps.append(idx_label)\n",
    "    return in_irreps\n",
    "            \n",
    "    \n",
    "trans_irs = decomp_irreps(sym_modes_trans)\n",
    "\n",
    "total_r2_irreps = None\n",
    "\n",
    "from posym import PointGroup, SymmetryObject\n",
    "\n",
    "for ir1_ in trans_irs:\n",
    "    ir1 = SymmetryObject(group=set_pg, rep=ir1_)\n",
    "    for ir2_ in trans_irs:\n",
    "        ir2 = SymmetryObject(group=set_pg, rep=ir2_)\n",
    "\n",
    "        if total_r2_irreps is None:\n",
    "            total_r2_irreps = ir1 * ir2\n",
    "        else:\n",
    "            total_r2_irreps = total_r2_irreps + ir1*ir2\n",
    "            \n",
    "# now remove the rotational ones:\n",
    "rot_irs = decomp_irreps(sym_modes_gs_rot)\n",
    "\n",
    "for ir_ in rot_irs:\n",
    "    ir = SymmetryObject(group=set_pg, rep=ir_)\n",
    "    total_r2_irreps = total_r2_irreps + (-1*ir)\n",
    "    \n",
    "print('IR active modes: ',sym_modes_trans)\n",
    "print('Raman active modes: ',total_r2_irreps)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9c081aa",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76f2e08d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}