orca-ir.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf560467-5f07-4318-ab26-ed44e44e1c79",
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "'''\n",
    "orca-ir Jupyter edition\n",
    "\n",
    "1. Upload 'orca.out' to the same folder\n",
    "2. Insert the filename here: imp_data('orca.out')\n",
    "3. Adjust the plot; use the interactive Matplotlib to select a specific region\n",
    "4. Save the spectrum by clicking on the 'disk' symbol of the interactive Matplotlib window\n",
    "\n",
    "'''\n",
    "\n",
    "# the matplotlip widget in IPython \n",
    "%matplotlib widget\n",
    "\n",
    "import re                                      # regular expressions\n",
    "import numpy as np                             # summation\n",
    "import matplotlib.pyplot as plt                # plots\n",
    "from scipy.signal import find_peaks            # peak detection\n",
    "from ipywidgets import interactive, HBox, VBox # interactive ui elements\n",
    "\n",
    "def imp_data(file):\n",
    "    specstring_start = 'IR SPECTRUM'          #check orca.out from here\n",
    "    specstring_end = 'The first'              #stop reading orca.out from here\n",
    "    # 3 empty lists (first one is not really necessary)\n",
    "    modelist = []\n",
    "    freqlist = []\n",
    "    intenslist = []\n",
    "    try:\n",
    "        with open(file, \"r\") as input_file:\n",
    "            for line in input_file:\n",
    "                #start exctract text \n",
    "                if \"Program Version 6\" in line:\n",
    "                    intens_column=3\n",
    "                if \"Program Version 5\" in line:\n",
    "                    intens_column=3\n",
    "                if \"Program Version 4\" in line:\n",
    "                    intens_column=2\n",
    "                if \"Program Version 3\" in line:\n",
    "                    intens_column=2\n",
    "                if line.startswith(specstring_start):\n",
    "                    #found IR data in orca.out\n",
    "                    for line in input_file:\n",
    "                        #stop exctract text \n",
    "                        if line.startswith(specstring_end):\n",
    "                            break\n",
    "                        #only recognize lines that start with number\n",
    "                        #split line into 3 lists mode, frequencies, intensities\n",
    "                        #line should start with a number\n",
    "                        if re.search(\"\\d:\",line): \n",
    "                            modelist.append(int(line.strip().split(\":\")[0])) \n",
    "                            freqlist.append(float(line.strip().split()[1]))\n",
    "                            intenslist.append(float(line.strip().split()[intens_column]))\n",
    "    # file not found -> message\n",
    "    except FileNotFoundError:\n",
    "        print(f\"'{file}'\" + \" not found\")\n",
    "    if modelist:\n",
    "        return modelist, freqlist, intenslist\n",
    "    # IR part not found -> message\n",
    "    else:\n",
    "        print(f\"'{specstring_start}'\" + \" not found in\" + f\"'{file}'\")\n",
    "    \n",
    "def roundup(x):\n",
    "    #round to next 100\n",
    "    return x if x % 100 == 0 else x + 100 - x % 100\n",
    "\n",
    "def gauss(a,m,x,w):\n",
    "    # calculation of the Gaussian line shape\n",
    "    # a = amplitude (max y, intensity)\n",
    "    # x = position\n",
    "    # m = maximum/meadian (stick position in x, wave number)\n",
    "    # w = line width, FWHM\n",
    "    return a*np.exp(-(np.log(2)*((m-x)/w)**2))\n",
    "\n",
    "# plot the IR spectrum\n",
    "def plot_ir(w = 15,                   \n",
    "            transm_style = True, \n",
    "            high_to_low_wn = True, \n",
    "            show_stick = True, \n",
    "            label_peaks = True,\n",
    "            show_grid = False,\n",
    "            show_single_gauss = False):\n",
    "    # w = line width\n",
    "    # high_to_low_wn = go from high to low wn \n",
    "    # show_stick = show the stick spectrum\n",
    "    # labe_peaks = label the peaks\n",
    "    # show_grid = show a grid\n",
    "    # show_single_gauss = show the single gauss functions\n",
    "\n",
    "    # always close before draw a new one\n",
    "    plt.close()\n",
    "\n",
    "    # list for sum of single gauss functions\n",
    "    gauss_sum = []\n",
    "    # position\n",
    "    freqlist = ir_data[1]\n",
    "    # intensity\n",
    "    intenslist = ir_data[2]\n",
    "    # add some extra space in x \n",
    "    wn_add = 120 + w * 2\n",
    "    # prepare plot\n",
    "    fig, ax = plt.subplots()\n",
    "    # plot range in x\n",
    "    plt_range_x=np.arange(0, max(freqlist) + wn_add, 1)\n",
    "\n",
    "    # plot single gauss function for every frequency freq if TRUE\n",
    "    # generate summation of single gauss functions\n",
    "    for index, freq in enumerate(freqlist):\n",
    "        if show_single_gauss:\n",
    "            ax.plot(plt_range_x,gauss(intenslist[index], plt_range_x, freq, w), alpha=0.5) \n",
    "            ax.fill_between(plt_range_x,gauss(intenslist[index], plt_range_x, freq, w), alpha=0.5)\n",
    "        gauss_sum.append(gauss(intenslist[index], plt_range_x, freq, w))\n",
    "    \n",
    "    # y values of the gauss summation\n",
    "    plt_range_gauss_sum_y = np.sum(gauss_sum, axis=0)\n",
    "    \n",
    "    # find peaks scipy function, change height for level of detection\n",
    "    peaks , _ = find_peaks(plt_range_gauss_sum_y, height = 0.1)\n",
    "\n",
    "    # plot the spectrum\n",
    "    ax.plot(plt_range_x, plt_range_gauss_sum_y, color=\"black\", linewidth=0.8)\n",
    "\n",
    "    # show stick spectrum if TRUE\n",
    "    if show_stick:\n",
    "        ax.stem(freqlist, intenslist, linefmt=\"dimgrey\", markerfmt=\" \", basefmt=\" \")\n",
    "\n",
    "    # title, label, ticks, linear locator\n",
    "    ax.set_title('IR spectrum', fontweight = 'bold')\n",
    "    ax.set_xlabel(r'$\\tilde{\\nu}$ /cm$^{-1}$')  \n",
    "    ax.get_yaxis().set_ticks([])   \n",
    "    ax.minorticks_on()\n",
    "    ax.xaxis.set_major_locator(plt.LinearLocator())\n",
    "\n",
    "    # transmission or absorption style\n",
    "    if transm_style:\n",
    "        plt.ylim(max(plt_range_gauss_sum_y)+max(plt_range_gauss_sum_y)*0.1,0)\n",
    "        ax.set_ylabel('transmittance')\n",
    "    else:\n",
    "        plt.ylim(0,max(plt_range_gauss_sum_y)+max(plt_range_gauss_sum_y)*0.1)\n",
    "        ax.set_ylabel('intensity')\n",
    "        label_rel_pos_y=5\n",
    "\n",
    "    # high to low wn or low to high wn\n",
    "    if high_to_low_wn:\n",
    "        plt.xlim(roundup(max(plt_range_x)),0) # round to next 100\n",
    "    else:\n",
    "        plt.xlim(0,roundup(max(plt_range_x)))\n",
    "\n",
    "    # show grid if TRUE\n",
    "    if show_grid:\n",
    "        ax.grid(True,which='major', axis='x',color='black',linestyle='dotted', linewidth=0.5)\n",
    "\n",
    "    # show labels on peaks if TRUE\n",
    "    if label_peaks:\n",
    "        for index, txt in enumerate(peaks):\n",
    "            \n",
    "            if transm_style:\n",
    "                #corr_factor - maintain distance from label to peak in transmittance style\n",
    "                #sensitive to peak label font size\n",
    "                corr_factor = (4 - len(str(peaks[index])))*3.75\n",
    "                label_rel_pos_y = -15\n",
    "                #if one does not care:\n",
    "                #corr_factor = 0 \n",
    "            else:\n",
    "                # not necessary for absorption style\n",
    "                label_rel_pos_y=5  \n",
    "                corr_factor = 0\n",
    "                \n",
    "            ax.annotate(peaks[index], xy=(peaks[index], plt_range_gauss_sum_y[peaks[index]]),\n",
    "                        ha = \"center\", rotation = 90, size = 6,\n",
    "                        xytext=(0,label_rel_pos_y+corr_factor), textcoords = 'offset points')\n",
    "\n",
    "    # resize the Matplotlib window\n",
    "    N = 1.5\n",
    "    params = plt.gcf()\n",
    "    plSize = params.get_size_inches()\n",
    "    params.set_size_inches((plSize[0]*N, plSize[1]*N))\n",
    "\n",
    "    # tight layout and show\n",
    "    plt.tight_layout() \n",
    "    plt.show()\n",
    "\n",
    "\n",
    "# upload orca.out in the same folder (location)\n",
    "# insert the filename; imp_data('orca.out)\n",
    "ir_data = imp_data('bf43.out')\n",
    "\n",
    "# make it interactive\n",
    "interactive_plot = interactive(plot_ir, w = (5, 100, 5), \n",
    "                               show_stick = True, \n",
    "                               transm_style = True, \n",
    "                               high_to_low_wn = True,\n",
    "                               label_peaks = True,\n",
    "                               show_grid = False,\n",
    "                               show_single_gauss = False)\n",
    "\n",
    "# ui elements (all to last -1)\n",
    "ui_el = interactive_plot.children[:-1]\n",
    "# plot window (last)\n",
    "plot_win = interactive_plot.children[-1]\n",
    "\n",
    "# ui elements in rows\n",
    "row1 = HBox(ui_el[:1])\n",
    "row2 = HBox(ui_el[1:3])\n",
    "row3 = HBox(ui_el[3:5])\n",
    "row4 = HBox(ui_el[5:])\n",
    "\n",
    "# display\n",
    "interactive_plot = VBox([row1, row2, row3, row4, plot_win])\n",
    "interactive_plot"
   ]
  }
 ],
 "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}