{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running vibronic calculations from a Molcas RASSI set of calculations\n",
    "\n",
    "Requirements:\n",
    "\n",
    " * Spin-free hamiltonian for every structure displaced along the normal modes for the numerical first derivative.\n",
    " * Zero-order file with the following spin-free properties for each of the vibronic coupling calculations:\n",
    "  - Electric dipole moments (`'electric_dipole'`)\n",
    "  - Magnetic dipole moments (`'magnetic_dipole'`)\n",
    "  - Electric quadrupole moments (`'electric_quadrupole'`)\n",
    " * Zero-order file with the spin-free and spin-orbit energies.\n",
    " \n",
    "The Zero-order file referenced above is the calculation at the equilibrium geometry with the `EJOB` option in the molcas input file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some imports for the resource to deal with the data files for the purpose of this example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vibrav.base import resource\n",
    "import tarfile\n",
    "import os\n",
    "import glob\n",
    "import shutil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "tarball = resource('molcas-ucl6-2minus-vibronic-coupling.tar.xz')\n",
    "with tarfile.open(tarball, 'r:xz') as tar:\n",
    "    tar.extractall()\n",
    "parent = os.getcwd()\n",
    "os.chdir('molcas-ucl6-2minus-vibronic-coupling')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Begin the vibronic coupling calculations\n",
    "\n",
    "The configuration file must be give as it will tell the program where to look for all of the things it needs. Refer to the `vibrav.vibronic.vibronic.Vibronic` documentation for more information about the required and default arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vibrav.vibronic import Vibronic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "vib = Vibronic(config_file='va.conf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Brief rundown of the most important input arguments\n",
    "\n",
    " - `property`: Tell the program which vibronic property to calculate.\n",
    " - `write_oscil`: If set to `True` the program will write the ocillator strengths for the x, y and z directions along with the averaged values in the `vibronic-outputs` directory under the filenames `oscillators-1.txt`, `oscillators-2.txt`, `oscillators-3.txt` and `oscillators-0.txt`, respectively.\n",
    " - `write_property`: If set to `True` the program will write the spin orbit vibronic properties that are calculated into the directories under the names `vib???` for the respective vibrational mode. It will then create two directories `plus` and `minus` referring to the transitions of `dE + E_vib` and `dE - E_vib`, respectively. Inside those directories it will write the spin-orbit property values for the three catrtesian directions.\n",
    "\n",
    "Refer to the `vibrav.vibronic.vibronic.Vibronic.vibronic_coupling` documentation for more information on the input parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/herbertl/github/vibrav/vibrav/numerical/boltzmann.py:114: Warning: Calculating only the first 2 states for the Boltzmann distribution.\n",
      "  warnings.warn(\"Calculating only the first {} \".format(states) \\\n"
     ]
    }
   ],
   "source": [
    "vib.vibronic_coupling(property='electric-dipole', print_stdout=False, temp=298,\n",
    "                      write_property=True, write_oscil=True, boltz_states=2,\n",
    "                      verbose=False, eq_cont=False, select_fdx=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Directory listing of the `vibronic-outputs` directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['oscillators-0.txt',\n",
       " 'oscillators-1.txt',\n",
       " 'oscillators-3.txt',\n",
       " 'boltzmann-populations.csv',\n",
       " 'oscillators-2.txt',\n",
       " 'alpha.txt']"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.listdir('vibronic-outputs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The different directories created for each vibrational mode. Note, this uses a one based index as we typically reserve 0 for the normal mode coordinate at the equilibrium structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['vib001',\n",
       " 'vib002',\n",
       " 'vib003',\n",
       " 'vib004',\n",
       " 'vib005',\n",
       " 'vib006',\n",
       " 'vib007',\n",
       " 'vib008',\n",
       " 'vib009',\n",
       " 'vib010',\n",
       " 'vib011',\n",
       " 'vib012',\n",
       " 'vib013',\n",
       " 'vib014',\n",
       " 'vib015']"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sorted(glob.glob('vib???'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Directory listing of one of the vibronic coupling property output directories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['dipole-1.txt', 'dipole-2.txt', 'dipole-3.txt', 'energies.txt']"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.listdir(os.path.join('vib001', 'plus'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading the `.txt` files that are generated\n",
    "\n",
    "There is a utility in the vibrav code that allows the user to easily open and read the `.txt` files that are generated by both molcas and vibrav. You cna find more information in the documentation for `vibrav.util.io.open_txt`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vibrav.util.io import open_txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reading all of the oscillator files that are generated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>nrow</th>\n",
       "      <th>ncol</th>\n",
       "      <th>oscil</th>\n",
       "      <th>energy</th>\n",
       "      <th>freqdx</th>\n",
       "      <th>sign</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>2.931277e-11</td>\n",
       "      <td>0.002609</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1.970756e-11</td>\n",
       "      <td>0.002609</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>3.586465e-11</td>\n",
       "      <td>0.002609</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>1.083229e-10</td>\n",
       "      <td>0.004795</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>2.710447e-11</td>\n",
       "      <td>0.001807</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   nrow  ncol         oscil    energy  freqdx   sign\n",
       "0     0     1  2.931277e-11  0.002609       0  minus\n",
       "1     0     2  1.970756e-11  0.002609       0  minus\n",
       "2     0     3  3.586465e-11  0.002609       0  minus\n",
       "3     0     4  1.083229e-10  0.004795       0  minus\n",
       "4     1     4  2.710447e-11  0.001807       0  minus"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "iso_oscil = open_txt(os.path.join('vibronic-outputs', 'oscillators-0.txt'),\n",
    "                     rearrange=False)\n",
    "iso_oscil.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>nrow</th>\n",
       "      <th>ncol</th>\n",
       "      <th>oscil</th>\n",
       "      <th>energy</th>\n",
       "      <th>freqdx</th>\n",
       "      <th>sign</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>5.378963e-12</td>\n",
       "      <td>0.002609</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1.240604e-11</td>\n",
       "      <td>0.002609</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>4.269799e-13</td>\n",
       "      <td>0.002609</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>9.580531e-11</td>\n",
       "      <td>0.004795</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>3.322260e-11</td>\n",
       "      <td>0.001807</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   nrow  ncol         oscil    energy  freqdx   sign\n",
       "0     0     1  5.378963e-12  0.002609       0  minus\n",
       "1     0     2  1.240604e-11  0.002609       0  minus\n",
       "2     0     3  4.269799e-13  0.002609       0  minus\n",
       "3     0     4  9.580531e-11  0.004795       0  minus\n",
       "4     1     4  3.322260e-11  0.001807       0  minus"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_oscil = open_txt(os.path.join('vibronic-outputs', 'oscillators-1.txt'),\n",
    "                   rearrange=False)\n",
    "x_oscil.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>nrow</th>\n",
       "      <th>ncol</th>\n",
       "      <th>oscil</th>\n",
       "      <th>energy</th>\n",
       "      <th>freqdx</th>\n",
       "      <th>sign</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>6.395906e-12</td>\n",
       "      <td>0.002609</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>9.883249e-15</td>\n",
       "      <td>0.002609</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>6.191645e-11</td>\n",
       "      <td>0.002609</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>7.761667e-11</td>\n",
       "      <td>0.004795</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>1.647976e-12</td>\n",
       "      <td>0.001807</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   nrow  ncol         oscil    energy  freqdx   sign\n",
       "0     0     1  6.395906e-12  0.002609       0  minus\n",
       "1     0     2  9.883249e-15  0.002609       0  minus\n",
       "2     0     3  6.191645e-11  0.002609       0  minus\n",
       "3     0     4  7.761667e-11  0.004795       0  minus\n",
       "4     1     4  1.647976e-12  0.001807       0  minus"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_oscil = open_txt(os.path.join('vibronic-outputs', 'oscillators-2.txt'),\n",
    "                   rearrange=False)\n",
    "y_oscil.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>nrow</th>\n",
       "      <th>ncol</th>\n",
       "      <th>oscil</th>\n",
       "      <th>energy</th>\n",
       "      <th>freqdx</th>\n",
       "      <th>sign</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>7.616344e-11</td>\n",
       "      <td>0.002609</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>4.670676e-11</td>\n",
       "      <td>0.002609</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>4.525051e-11</td>\n",
       "      <td>0.002609</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>1.515468e-10</td>\n",
       "      <td>0.004795</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>4.644282e-11</td>\n",
       "      <td>0.001807</td>\n",
       "      <td>0</td>\n",
       "      <td>minus</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   nrow  ncol         oscil    energy  freqdx   sign\n",
       "0     0     1  7.616344e-11  0.002609       0  minus\n",
       "1     0     2  4.670676e-11  0.002609       0  minus\n",
       "2     0     3  4.525051e-11  0.002609       0  minus\n",
       "3     0     4  1.515468e-10  0.004795       0  minus\n",
       "4     1     4  4.644282e-11  0.001807       0  minus"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z_oscil = open_txt(os.path.join('vibronic-outputs', 'oscillators-3.txt'),\n",
    "                   rearrange=False)\n",
    "z_oscil.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cleaning up for the purpose of the example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(parent)\n",
    "shutil.rmtree('molcas-ucl6-2minus-vibronic-coupling')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}