docs/doc-source/ePS_ePSproc_tutorial/ePS_tutorial_080520.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ePolyScat basics tutorial\n",
    "08/05/20\n",
    "\n",
    "Paul Hockett\n",
    "\n",
    "Disclaimer: I am an enthusiastic ePolyScat user for photoionization calculations, but not an expert on the code. Nonetheless, this tutorial aims to go over some of the key features/uses of ePS for such problems - as far as my own usage goes - and provide an introduction and resource to new users."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview & resources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* [ePolyScat](https://epolyscat.droppages.com) (ePS) is an open-source tool for numerical computation of electron-molecule scattering & photoionization by Lucchese & coworkers. For more details: \n",
    "    \n",
    "    * The [ePolyScat website and manual](https://epolyscat.droppages.com). *Note* that the manual uses a frames-based layout, use the menu on the page to navigate to sub-sections (direct links break the menu).\n",
    "    * *Calculation of low-energy elastic cross sections for electron-CF4 scattering*, F. A. Gianturco, R. R. Lucchese, and N. Sanna, J. Chem. Phys. 100, 6464 (1994), http://dx.doi.org/10.1063/1.467237\n",
    "    * *Cross section and asymmetry parameter calculation for sulfur 1s photoionization of SF6*, A. P. P. Natalense and R. R. Lucchese, J. Chem. Phys. 111, 5344 (1999), http://dx.doi.org/10.1063/1.479794\n",
    "    * *Applications of the Schwinger variational principle to electron-molecule collisions and molecular photoionization.* Lucchese, R. R., Takatsuka, K., & McKoy, V. (1986).  Physics Reports, 131(3), 147–221. https://doi.org/10.1016/0370-1573(86)90147-X (comprehensive discussion of the theory and methods underlying the code).\n",
    "    \n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* [ePSproc](https://epsproc.readthedocs.io) is an open-source tool for post-processing & visualisation of ePS results, aimed primarily at photoionization studies.\n",
    "    \n",
    "    * Ongoing documentation is on [Read the Docs](https://epsproc.readthedocs.io).\n",
    "    * Source code is [available on Github](https://github.com/phockett/ePSproc).\n",
    "    * For more background, see the software metapaper for the original release of ePSproc (Aug. 2016): *ePSproc: Post-processing suite for ePolyScat electron-molecule scattering calculations*, on [Authorea](https://www.authorea.com/users/71114/articles/122402/_show_article) or [arXiv 1611.04043](https://arxiv.org/abs/1611.04043).\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* [ePSdata](https://phockett.github.io/ePSdata/index.html) is an open-data/open-science collection of ePS + ePSproc results.\n",
    "\n",
    "    * ePSdata collects ePS datasets, post-processed via ePSproc (Python) in [Jupyter notebooks](https://jupyter.org), for a full open-data/open-science transparent pipeline.\n",
    "    * ePSdata is currently (Jan 2020) collecting existing calculations from 2010 - 2019, from the [femtolabs at NRC](http://femtolab.ca), with one notebook per ePS job.\n",
    "    * In future, ePSdata pages will be automatically generated from ePS jobs (via the ePSman toolset, currently in development), for immediate dissemination to the research community.\n",
    "    * Source notebooks are available on the [Github project pages](https://github.com/phockett/ePSdata/), and notebooks + datasets via [Zenodo repositories](https://about.zenodo.org) (one per dataset). Each notebook + dataset is given a Zenodo DOI for full traceability, and notebooks are versioned on Github.\n",
    "    * Note: ePSdata may also be linked or mirrored on the existing [ePolyScat Collected Results OSF project](https://osf.io/psjxt/), but will effectively supercede those pages.\n",
    "    * All results are released under <a href=\"https://creativecommons.org/licenses/by-nc-sa/4.0/\">Creative Commons Attribution-NonCommercial-ShareAlike 4.0 (CC BY-NC-SA 4.0)</a> license, and are part of our ongoing [Open Science initiative](http://femtolab.ca/?p=877)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Workflow\n",
    "\n",
    "The general workflow for photoionization calculations plus post-processing is shown below. This pipeline involves a range of code suites, as shown in the main workflow; some additional details are also illustrated. This tutorial will only discuss ePS."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![ePS workflow](https://phockett.github.io/ePSdata/_images/ePSworkflow_170220.gv.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Theoretical background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For general scattering theory notes, try these textbooks:\n",
    "\n",
    "* *Quantum Mechanics Volume I.* Messiah, A. (1970).  North-Holland Publishing Company.\n",
    "* *Introduction to the Quantum Theory of Scattering.* Rodberg, L. S., & Thaler, R. M. (1967).  Academic Press.\n",
    "\n",
    "For notes focussed on photoionization problems, try:\n",
    "\n",
    "* [Photoelectron angular distributions seminar series](https://doi.org/10.6084/m9.figshare.c.3511731)\n",
    "* [Brief notes on scattering theory for photoionization](https://github.com/phockett/ePSproc/blob/master/docs/additional/ePSproc_scattering_theory_ePS_notes_011015.pdf)\n",
    "* *Quantum Metrology with Photoelectrons, Volume 1 Foundations.* Hockett, P. (2018). IOP Publishing. https://doi.org/10.1088/978-1-6817-4684-5\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For ePS, the details can be found in:\n",
    "\n",
    "* *Applications of the Schwinger variational principle to electron-molecule collisions and molecular photoionization.* Lucchese, R. R., Takatsuka, K., & McKoy, V. (1986).  Physics Reports, 131(3), 147–221. https://doi.org/10.1016/0370-1573(86)90147-X (comprehensive discussion of the theory and methods underlying the code).\n",
    "* *Calculation of low-energy elastic cross sections for electron-CF4 scattering*, F. A. Gianturco, R. R. Lucchese, and N. Sanna, J. Chem. Phys. 100, 6464 (1994), http://dx.doi.org/10.1063/1.467237\n",
    "* *Cross section and asymmetry parameter calculation for sulfur 1s photoionization of SF6*, A. P. P. Natalense and R. R. Lucchese, J. Chem. Phys. 111, 5344 (1999), http://dx.doi.org/10.1063/1.479794\n",
    "\n",
    "\n",
    "**The Method section from the latter paper is reproduced here for reference:**\n",
    "\n",
    "\n",
    "---\n",
    "\n",
    "Both the initial neutral molecule electronic wave function and the\n",
    "final ionized molecule electronic wave function are represented by\n",
    "single determinants constructed using the Hartree--Fock orbitals\n",
    "of the initial neutral state. The final $N$--electron state (continuum\n",
    "photoelectron + molecular ion) can then be variationally obtained\n",
    "from:\n",
    "\n",
    "\\begin{equation}\n",
    "\\langle\\delta\\Psi_{\\mathbf{k}}|H-E|\\Psi_{\\mathbf{k}}\\rangle=0,\n",
    "\\end{equation}\n",
    "\n",
    "where $\\delta\\Psi_{\\mathbf{k}}$ represent variations in the final\n",
    "state due to variations on the continuum wave function $\\psi_{\\mathbf{k}}(\\mathbf{r})$,\n",
    "and $\\mathbf{k}$ is the momentum of the ejected electron. The final\n",
    "photoionization problem is then reduced to solving the problem of\n",
    "an electron under the action of the potential of the ion. Thus, we\n",
    "do not consider many-electron effects. The following scattering equation\n",
    "can be obtained from Eq. (1) (Ref. 19) (in atomic units),\n",
    "\n",
    "\\begin{equation}\n",
    "\\left[-\\frac{1}{2}\\nabla^{2}+V(\\mathbf{r})-\\frac{k^{2}}{2}\\right]\\psi_{\\mathbf{k}}^{(\\pm)}(\\mathbf{r})=0,\n",
    "\\end{equation}\n",
    "\n",
    "where $V(\\mathbf{r})$ is the static exchange potential given by:\n",
    "\n",
    "\\begin{equation}\n",
    "V(\\mathbf{r})=-\\sum_{\\gamma=1}^{M}Z_{\\gamma}|\\mathbf{r}-\\mathbf{R}_{\\gamma}|^{-1}+\\sum_{i=1}^{n_{occ}}(2\\hat{J}_{i}-\\hat{K}_{i})+\\hat{J}_{\\perp5}+\\hat{K}_{\\perp5}\n",
    "\\end{equation}\n",
    "\n",
    "for $M$ nuclei of charge $Z_{\\gamma}$ located at $\\mathbf{R}\\gamma$\n",
    ", where $n_{occ}$ is the number of doubly occupied orbitals. $\\hat{J}_{i}$\n",
    "is the Coulomb operator,\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{J}_{i}(\\mathbf{r}_{1})=\\int\\frac{\\phi_{i}^{*}(\\mathbf{r}_{2})\\phi_{i}(\\mathbf{r}_{2})}{r_{12}}d^{3}r_{2},\n",
    "\\end{equation}\n",
    "\n",
    "and $\\hat{K}_{i}$ is the nonlocal exchange potential operator\n",
    "\n",
    "\\begin{equation}\n",
    "(\\hat{K}_{i}\\psi)(\\mathbf{r}_{1})=\\phi_{i}(\\mathbf{r}_{2})\\int\\frac{\\phi_{i}^{*}(\\mathbf{r}_{2})\\psi(\\mathbf{r}_{2})}{r_{12}}d^{3}r_{2}\n",
    "\\end{equation}\n",
    "\n",
    "Correlation and polarization effects can be included in the calculations\n",
    "through the addition of a local, energy independent model correlation\n",
    "polarization potential, such as described in Ref. 10 {[}...{]}\n",
    "\n",
    "To solve the scattering problem, we use the single center expansion\n",
    "(SCE) method, where all three dimensional functions are expanded on\n",
    "a set of angular functions $X_{lh}^{p\\mu}(\\theta,\\phi)$, which are\n",
    "symmetry adapted, according to the irreducible representations of\n",
    "the molecular point group. An arbitrary three dimensional function\n",
    "$F^{p\\mu}(r,\\theta,\\phi)$ is then expanded as\n",
    "\n",
    "\\begin{equation}\n",
    "F^{p\\mu}(r,\\theta,\\phi)=\\sum_{lh}r^{-1}f_{lh}^{p\\mu}(r)X_{lh}^{p\\mu}(\\theta,\\phi),\n",
    "\\end{equation}\n",
    "\n",
    "where\n",
    "\n",
    "\\begin{equation}\n",
    "X_{lh}^{p\\mu}(\\theta,\\phi)=\\sum_{m}b_{lhm}^{p\\mu}Y_{lm}(\\theta,\\phi),\n",
    "\\end{equation}\n",
    "\n",
    "and $p$ is one of the irreducible representations of the molecular\n",
    "point group, $m$ is a component of the representation $p$, and $h$\n",
    "indexes all possible $X_{lh}^{p\\mu}$ belonging to the same irreducible\n",
    "representation ($p\\mu$) with the same value of $l$. \\emph{The radial\n",
    "functions $f_{lh}^{p\\mu}(r)$ are represented on a numerical grid}.\n",
    "When solving the scattering equations we enforce orthogonality between\n",
    "the continuum solutions and the occupied orbitals.$^{19}$\n",
    "\n",
    "The matrix elements of the dipole operator are\n",
    "\n",
    "\\begin{equation}\n",
    "I_{\\mathbf{k},\\hat{n}}^{L}=(k)^{1/2}\\langle\\Psi_{i}|\\mathbf{r}.\\hat{n}|\\Psi_{f,\\mathbf{k}}^{(-)}\\rangle\n",
    "\\end{equation}\n",
    "\n",
    "for the dipole length form, and\n",
    "\n",
    "\\begin{equation}\n",
    "I_{\\mathbf{k},\\hat{n}}^{V}=\\frac{(k)^{1/2}}{E}\\langle\\Psi_{i}|\\nabla.\\hat{n}|\\Psi_{f,\\mathbf{k}}^{(-)}\\rangle\n",
    "\\end{equation}\n",
    "\n",
    "for the dipole velocity form, where $|\\Psi_{i}\\rangle$ is the initial\n",
    "bound state, $|\\Psi_{f,\\mathbf{k}}^{(-)}\\rangle$ is the final continuum\n",
    "state, $E$ is the photon energy, $\\mathbf{k}$ is the momentum of\n",
    "the photoelectron, and $\\hat{n}$ is the direction of polarization\n",
    "of the light, which is assumed to be linearly polarized.\n",
    "\n",
    "The matrix elements $I_{\\mathbf{k},\\hat{n}}^{(L,V)}$ of Eqs. (8)\n",
    "and (9) can be expanded in terms of the $X_{lh}^{p\\mu}$ functions\n",
    "of Eq. (7) as$^{14}$\n",
    "\n",
    "\\begin{equation}\n",
    "I_{\\mathbf{k},\\hat{n}}^{(L,V)}=\\left[\\frac{4\\pi}{3}\\right]^{1/2}\\sum_{p\\mu lhv}I_{lhv}^{p\\mu(L,V)}X_{lh}^{p\\mu}(\\hat{k})X_{1v}^{p_{v}\\mu_{v}}(\\hat{n}).\n",
    "\\end{equation}\n",
    "\n",
    "{[}Note here the final term gives polarization (dipole) terms, with\n",
    "$l=1$, $h=v$, corresponding to a photon with one unit of angular\n",
    "momentum and projections $v=-1,0,1$, correlated with irreducible\n",
    "representations $p_{v}\\mu_{v}$.{]}\n",
    "\n",
    "The differential cross section is given by\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{d\\sigma^{L,V}}{d\\Omega_{\\mathbf{k}}}=\\frac{\\sigma^{L,V}}{4\\pi}[1+\\beta_{\\mathbf{k}}^{L,V}P_{2}(\\cos\\theta)],\n",
    "\\end{equation}\n",
    "\n",
    "where the asymmetry parameter can be written as$^{14}$\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\beta_{\\mathbf{k}}^{L,V} & = & \\frac{3}{5}\\frac{1}{\\sum_{p\\mu lhv}|I_{\\mathbf{k},\\hat{n}}^{(L,V)}|^{2}}\\sum_{\\stackrel{p\\mu lhvmm_{v}}{p'\\mu'l'h'v'm'm'_{v}}}(-1)^{m'-m_{v}}I_{\\mathbf{k},\\hat{n}}^{(L,V)}\\nonumber \\\\\n",
    " & \\times & \\left(I_{\\mathbf{k},\\hat{n}}^{(L,V)}\\right)^{*}b_{lhm}^{p\\mu}b_{l'h'm'}^{p'\\mu'*}b_{1vm_{v}}^{p_{v}\\mu_{v}}b_{1v'm'_{v}}^{p'_{v}\\mu'_{v}*}\\nonumber \\\\\n",
    " & \\times & [(2l+1)(2l'+1)]^{1/2}(1100|20)(l'l00|20)\\nonumber \\\\\n",
    " & \\times & (11-m'_{v}m_{v}|2M')(l'l-m'm|2-M'),\n",
    "\\end{eqnarray}\n",
    "\n",
    "and the $(l'lm'm|L'M')$ are the usual Clebsch--Gordan coefficients.\n",
    "The total cross section is\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma^{L,V}=\\frac{4\\pi^{2}}{3c}E\\sum_{p\\mu lhv}|I_{\\mathbf{k},\\hat{n}}^{(L,V)}|^{2},\n",
    "\\end{equation}\n",
    "\n",
    "where c is the speed of light.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical background\n",
    "\n",
    "Numerically, ePS proceeds approximately as in the flow-chart above: \n",
    "\n",
    "- Molecular structure is parsed from the input quantum chemistry results;\n",
    "- Various functions process this input, plus additional user commands, to set up the problem (e.g. ionizing orbital selection, symmetry-adapted single point expansion, and subsequent determination of the scattering potential $V(\\mathbf{r})$);\n",
    "- The scattering solution is determined variationally (i.e. solving $\\langle\\delta\\Psi_{\\mathbf{k}}|H-E|\\Psi_{\\mathbf{k}}\\rangle=0$) - this is essentially the [Schwinger variational method](https://en.wikipedia.org/wiki/Schwinger_variational_principle), see [ePS references at head of section for further details](#Theoretical-background);\n",
    "- Various outputs are written to file - some examples are discussed below.\n",
    "\n",
    "It's the variational solution for the continuum wavefunction  that (usually) takes most of the computational effort, and is the core of an ePS run. This procedure can take anywhere from seconds to minutes to hours for a single point (== single input geometry, energy and symmetry), depending on the problem at hand, and the hardware. \n",
    "\n",
    "In terms of the code, ePS is written in fortran 90, with MPI for parallelism, and LAPACK and BLAS libraries for numerical routines - see the [intro section of the manual](https://epolyscat.droppages.com/IntroREADME2) for more details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ePolyScat basic example calculation: N2 $3\\sigma_g^{-1}$ photoionization\n",
    "\n",
    "ePS uses a similar input format to standard quantum chemistry codes, with a series of data records and commands set by the user for a specific job.\n",
    "\n",
    "As a basic intro, here is `test12` from the [ePS sample jobs](https://epolyscat.droppages.com/TestsREADME):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "#\n",
    "# input file for test12\n",
    "#\n",
    "# N2 molden SCF, (3-sigma-g)^-1 photoionization\n",
    "#\n",
    "  LMax   22     # maximum l to be used for wave functions\n",
    "  EMax  50.0    # EMax, maximum asymptotic energy in eV\n",
    "  FegeEng 13.0   # Energy correction (in eV) used in the fege potential\n",
    "\n",
    "  ScatEng  10.0   # list of scattering energies\n",
    "\n",
    " InitSym 'SG'      # Initial state symmetry\n",
    " InitSpinDeg 1     # Initial state spin degeneracy\n",
    " OrbOccInit 2 2 2 2 2 4  # Orbital occupation of initial state\n",
    " OrbOcc     2 2 2 2 1 4  # occupation of the orbital groups of target\n",
    " SpinDeg 1         # Spin degeneracy of the total scattering state (=1 singlet)\n",
    " TargSym 'SG'      # Symmetry of the target state\n",
    " TargSpinDeg 2     # Target spin degeneracy\n",
    " IPot 15.581    # ionization potentail\n",
    "\n",
    "Convert '$pe/tests/test12.molden' 'molden'\n",
    "GetBlms\n",
    "ExpOrb\n",
    "\n",
    " ScatSym     'SU'  # Scattering symmetry of total final state\n",
    " ScatContSym 'SU'  # Scattering symmetry of continuum electron\n",
    "\n",
    "FileName 'MatrixElements' 'test12SU.idy' 'REWIND'\n",
    "GenFormPhIon\n",
    "DipoleOp\n",
    "GetPot\n",
    "PhIon\n",
    "GetCro\n",
    "#\n",
    " ScatSym     'PU'  # Scattering symmetry of total final state\n",
    " ScatContSym 'PU'  # Scattering symmetry of continuum electron\n",
    "\n",
    "FileName 'MatrixElements' 'test12PU.idy' 'REWIND'\n",
    "GenFormPhIon\n",
    "DipoleOp\n",
    "GetPot\n",
    "PhIon\n",
    "GetCro\n",
    "#\n",
    "GetCro 'test12PU.idy' 'test12SU.idy'\n",
    "#\n",
    "#\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's breakdown the various segments to this input file..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### System & job definition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Basic file header info, comments with #.\n",
    "\n",
    "```\n",
    "#\n",
    "# input file for test12\n",
    "#\n",
    "# N2 molden SCF, (3-sigma-g)^-1 photoionization\n",
    "#\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Defined some numerical values for the calculation - the values here are defined in the `Data Records` [section of the manual](https://epolyscat.droppages.com/SourceGetDataRecordDef_f90), and may have defaults if not set.\n",
    "```\n",
    "  LMax   22     # maximum l to be used for wave functions\n",
    "  EMax  50.0    # EMax, maximum asymptotic energy in eV\n",
    "  FegeEng 13.0   # Energy correction (in eV) used in the fege potential\n",
    "\n",
    "  ScatEng  10.0   # list of scattering energies\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define symmetries and orbital occupations. Note that the orbitals are grouped by degeneracy in ePS, so the numbering here may be different from that in the raw computational chemistry file output (which is typically *not* grouped).\n",
    "\n",
    "```\n",
    " InitSym 'SG'      # Initial state symmetry\n",
    " InitSpinDeg 1     # Initial state spin degeneracy\n",
    " OrbOccInit 2 2 2 2 2 4  # Orbital occupation of initial state\n",
    " OrbOcc     2 2 2 2 1 4  # occupation of the orbital groups of target\n",
    " SpinDeg 1         # Spin degeneracy of the total scattering state (=1 singlet)\n",
    " TargSym 'SG'      # Symmetry of the target state\n",
    " TargSpinDeg 2     # Target spin degeneracy\n",
    " IPot 15.581    # ionization potentail\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the symmetries set here correspond to specific orbital and continuum sets, so there may be multiple symmetries for a given problem - more on this later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Init job & run calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define electronic structure (in this case, from a Molden file), and run a single-point expansion in symmetry adapted harmonics (see `Commands` [manual pages for more details](https://epolyscat.droppages.com/SourceePolyScat_f90)).\n",
    "\n",
    "```\n",
    "Convert '$pe/tests/test12.molden' 'molden'   # Read electronic structure file\n",
    "GetBlms\n",
    "ExpOrb\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The meat of ePS is running scattering computations, and calculating associated matrix elements/parameters/observables. In this case, there are two (continuum) symmetries, (SU, PU). In each case, the data records are set, and a sequence of commands run the computations to numerically determine the continuum (scattering) wavefunction (again, more details can be found via the `Commands` [manual pages](https://epolyscat.droppages.com/SourceePolyScat_f90), comments added below are the basic command descriptions).\n",
    "\n",
    "```\n",
    " ScatSym     'SU'  # Scattering symmetry of total final state\n",
    " ScatContSym 'SU'  # Scattering symmetry of continuum electron\n",
    "\n",
    "FileName 'MatrixElements' 'test12SU.idy' 'REWIND'   # Set file for matrix elements\n",
    "GenFormPhIon   # Generate potential formulas for photoioniztion.\n",
    "DipoleOp       # Compute the dipole operator onto an orbital. \n",
    "GetPot         # Calculate electron density, static potential, and V(CP) potential. \n",
    "PhIon          # Calculate photionization dipole matrix elements.\n",
    "GetCro         # Compute photoionization cross sections from output of scatstab or from a file of dynamical coefficients.\n",
    "\n",
    "# Repeat for 2nd symmetry\n",
    "#\n",
    " ScatSym     'PU'  # Scattering symmetry of total final state\n",
    " ScatContSym 'PU'  # Scattering symmetry of continuum electron\n",
    "\n",
    "FileName 'MatrixElements' 'test12PU.idy' 'REWIND'\n",
    "GenFormPhIon\n",
    "DipoleOp\n",
    "GetPot\n",
    "PhIon\n",
    "GetCro\n",
    "\n",
    "# Run a final GetCro command including both symmetries, set by specifying the matrix element files to use.\n",
    "#\n",
    "GetCro 'test12PU.idy' 'test12SU.idy'\n",
    "#\n",
    "#\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running ePolyScat\n",
    "\n",
    "Assuming you have ePS compiled/installed, then running the above file is a simple case of passing the file to the ePolyScat executable, e.g. `/opt/ePolyScat.E3/bin/ePolyScat inputFile`.\n",
    "\n",
    "For a general [code overview see the manual](https://epolyscat.droppages.com), contact R. R. Lucchese for source code.\n",
    "\n",
    "Performance will depend on the machine, but for a decent multi-core workstation expect minutes to hours per energy point depending on the size and symmetry of the problem (again, see the [test/example jobs](https://epolyscat.droppages.com/TestsREADME) for some typical timings.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results\n",
    "Again following quantum chemistry norms, the main output file is an ASCII file with various sections, separated by keyword headers, corresponding different steps in the computations. Additional files - e.g. the matrix elements listed above - may also be of use, although are typically not user-readable/interpretable, but rather provided for internal use (e.g. running further calculations on matrix elements, [used with utility programmes](https://epolyscat.droppages.com/UtilityProgs), etc.).\n",
    "\n",
    "The [full test12 output file can be found in the manual](https://epolyscat.droppages.com/Tests_test12_ostnd).  We'll just look at the `CrossSection` segments, correlated with the `GetCro` commands in the input, which give the calculated $\\sigma^{L,V}$ and $\\beta_{\\mathbf{k}}^{L,V}$ values (as [per definitions above](#Theoretical-background))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the first exampe, for SU symmetry\n",
    "\n",
    "```\n",
    "----------------------------------------------------------------------\n",
    "CrossSection - compute photoionization cross section\n",
    "----------------------------------------------------------------------\n",
    "\n",
    "Ionization potential (IPot) =     15.5810 eV\n",
    "Label -\n",
    "Cross section by partial wave      F\n",
    "Cross Sections for\n",
    "\n",
    "     Sigma LENGTH   at all energies\n",
    "      Eng\n",
    "    25.5810  0.57407014E+01\n",
    "\n",
    "     Sigma MIXED    at all energies\n",
    "      Eng\n",
    "    25.5810  0.53222735E+01\n",
    "\n",
    "     Sigma VELOCITY at all energies\n",
    "      Eng\n",
    "    25.5810  0.49344024E+01\n",
    "\n",
    "     Beta LENGTH   at all energies\n",
    "      Eng\n",
    "    25.5810  0.47433563E+00\n",
    "\n",
    "     Beta MIXED    at all energies\n",
    "      Eng\n",
    "    25.5810  0.47537873E+00\n",
    "\n",
    "     Beta VELOCITY at all energies\n",
    "      Eng\n",
    "    25.5810  0.47642429E+00\n",
    "\n",
    "          COMPOSITE CROSS SECTIONS AT ALL ENERGIES\n",
    "         Energy  SIGMA LEN  SIGMA MIX  SIGMA VEL   BETA LEN   BETA MIX   BETA VEL\n",
    "EPhi     25.5810     5.7407     5.3223     4.9344     0.4743     0.4754     0.4764\n",
    "Time Now =        19.9228  Delta time =         0.0067 End CrossSection\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This provides a listing of the photoionization cross-section(s) $\\sigma$ (in [mega barns]((https://en.wikipedia.org/wiki/Barn_(unit))), 1 Mb = $10^{-18} cm^2$), and associated anisotropy parameter(s) $\\beta$ (dimensionless) - in this case, just for a single energy point. Note that there are multiple values, these correspond to calculations in $L$ or $V$ gauge. The values are also tabulated at the end of the output, along with a note on the computational time. \n",
    "\n",
    "Note that these values correspond to calculations for an isotropic ensemble, i.e. what would usually be measured in the lab frame for a 1-photon ionization from a gas sample, as given by the $\\beta_{\\mathbf{k}}^{L,V}$ values [defined earlier](#Theoretical-background)). For these normalised values, the differential cross section (photoelectron flux vs. angle) $\\frac{d\\sigma^{L,V}}{d\\Omega_{\\mathbf{k}}}$ [defined earlier](#Theoretical-background) can be written in a slightly simplified form, $I(\\theta) = \\sigma(1 + \\beta \\cos^2(\\theta))$ (note this is, by definition, cylindrically symmetric)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the values corresponding to PU symmetry:\n",
    "\n",
    "```\n",
    "----------------------------------------------------------------------\n",
    "CrossSection - compute photoionization cross section\n",
    "----------------------------------------------------------------------\n",
    "\n",
    "Ionization potential (IPot) =     15.5810 eV\n",
    "Label -\n",
    "Cross section by partial wave      F\n",
    "Cross Sections for\n",
    "\n",
    "     Sigma LENGTH   at all energies\n",
    "      Eng\n",
    "    25.5810  0.30571336E+01\n",
    "\n",
    "     Sigma MIXED    at all energies\n",
    "      Eng\n",
    "    25.5810  0.27828683E+01\n",
    "\n",
    "     Sigma VELOCITY at all energies\n",
    "      Eng\n",
    "    25.5810  0.25404131E+01\n",
    "\n",
    "     Beta LENGTH   at all energies\n",
    "      Eng\n",
    "    25.5810  0.12686831E+01\n",
    "\n",
    "     Beta MIXED    at all energies\n",
    "      Eng\n",
    "    25.5810  0.13002321E+01\n",
    "\n",
    "     Beta VELOCITY at all energies\n",
    "      Eng\n",
    "    25.5810  0.13302595E+01\n",
    "\n",
    "          COMPOSITE CROSS SECTIONS AT ALL ENERGIES\n",
    "         Energy  SIGMA LEN  SIGMA MIX  SIGMA VEL   BETA LEN   BETA MIX   BETA VEL\n",
    "EPhi     25.5810     3.0571     2.7829     2.5404     1.2687     1.3002     1.3303\n",
    "Time Now =        35.0804  Delta time =         0.0067 End CrossSection\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and the final set of values for both continua:\n",
    "    \n",
    "```\n",
    "----------------------------------------------------------------------\n",
    "CrossSection - compute photoionization cross section\n",
    "----------------------------------------------------------------------\n",
    "\n",
    "Ionization potential (IPot) =     15.5810 eV\n",
    "Label -\n",
    "Cross section by partial wave      F\n",
    "Cross Sections for\n",
    "\n",
    "     Sigma LENGTH   at all energies\n",
    "      Eng\n",
    "    25.5810  0.87978350E+01\n",
    "\n",
    "     Sigma MIXED    at all energies\n",
    "      Eng\n",
    "    25.5810  0.81051418E+01\n",
    "\n",
    "     Sigma VELOCITY at all energies\n",
    "      Eng\n",
    "    25.5810  0.74748155E+01\n",
    "\n",
    "     Beta LENGTH   at all energies\n",
    "      Eng\n",
    "    25.5810  0.10200847E+01\n",
    "\n",
    "     Beta MIXED    at all energies\n",
    "      Eng\n",
    "    25.5810  0.10432056E+01\n",
    "\n",
    "     Beta VELOCITY at all energies\n",
    "      Eng\n",
    "    25.5810  0.10660644E+01\n",
    "\n",
    "          COMPOSITE CROSS SECTIONS AT ALL ENERGIES\n",
    "         Energy  SIGMA LEN  SIGMA MIX  SIGMA VEL   BETA LEN   BETA MIX   BETA VEL\n",
    "EPhi     25.5810     8.7978     8.1051     7.4748     1.0201     1.0432     1.0661\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the SU continua has a larger XS (~5.7 Mb) and a larger $\\beta$ value, hence will have a more anisotropic angular scattering distribution. The values for the different gauges are similar, indicating that these results should be (numerically) accurate - generally speaking, the differences in results between different gauges can be regarded as indicative of numerical accuracy and stability. (See, for example [*Atoms and Molecules in Intense Laser Fields:  Gauge Invariance ofTheory and Models*, A. D. Bandrauk, F. Fillion-Gourdeau and E. Lorin, arXiv 1302.2932](https://arxiv.org/pdf/1302.2932.pdf) (best ref...?).)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ePolyScat: preparing an input file\n",
    "\n",
    "In most cases, this is relatively simple, assuming that there is a suitable test example to use as a template. A few caveats...\n",
    "\n",
    "- There are quite a few computational parameters which can, and should, be played with. (I have *not* done enough of this, to be honest.)\n",
    "- There are lots of commands which are not really covered in the test examples, so exploring the manual is worthwhile.\n",
    "- One aspect which is non-trivial is the assignment of correct symmetries. To the best of my knowledge (which doesn't go very far here), there is no easy way to auto-generate these (although ePS will tell you if you got things wrong...), so maybe this can be considered as the necessary barrier to entry...! Some examples are given below.\n",
    "\n",
    "(Some autogeneration for input files & general job management is currently in development.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Symmetry selection rules in photoionization\n",
    "\n",
    "Essentially, the relevant direct products must contain the totally symmetric representation of the point group to constitute an allowed combination:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\Gamma_{\\mathrm{ion}}\\otimes\\Gamma_{\\mathrm{electron}}\\otimes\\Gamma_{\\mathrm{dipole}}\\otimes\\Gamma_{\\mathrm{neutral}} & \\supseteq & \\mathrm{A_{1}}\\\\\n",
    "\\end{eqnarray}\n",
    "\n",
    "Where $\\Gamma$ is the character of the quantity of interest. For the neutral and ion this will be the direct product of unpaired electrons - hence, for a closed-shell system, this is usually totally symmetric for the neutral, and identical to the character of the ionizing orbital for the ion. The dipole character corresponds to the $(x,y,z)$ operators in the point group, and the electron (continuum) character is what one needs to work out."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more discussion, see, for example, [Signorell, R., & Merkt, F. (1997). General symmetry selection rules for the photoionization of polyatomic molecules. Molecular Physics, 92(5), 793–804. DOI: 10.1080/002689797169745](https://doi.org/10.1080/002689797169745)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Worked example: N2 $3\\sigma_g^{-1}$\n",
    "\n",
    "This corresponds to the test case above. We'll need the [character table](http://symmetry.jacobs-university.de/cgi-bin/group.cgi?group=1001&option=4) for $D_{\\infty h}$ (see [also Oxford materials page from Atkins, Child & Philips](https://global.oup.com/uk/orc/chemistry/qchem2e/student/tables/), which includes the direct product tables ([PDF version](https://fdslive.oup.com/www.oup.com/orc/resources/chemistry/qchem2e/student/tables/tables_for_group_theory.pdf)). Then just plug in and work through..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{eqnarray}\n",
    "\\Gamma_{\\mathrm{ion}}\\otimes\\Gamma_{\\mathrm{electron}}\\otimes\\Gamma_{\\mathrm{dipole}}\\otimes\\Gamma_{\\mathrm{neutral}} & \\supseteq & \\Sigma_{g}^{+}\\\\\n",
    "\\Sigma_{g}^{+}\\otimes\\Gamma_{\\mathrm{electron}}\\otimes\\begin{array}{c}\n",
    "\\Pi_{u}(x,y)\\\\\n",
    "\\Sigma_{u}^{+}(z)\n",
    "\\end{array}\\otimes\\Sigma_{g}^{+} & \\supseteq & \\Sigma_{g}^{+}\\\\\n",
    "\\Sigma_{g}^{+}\\otimes\\Gamma_{\\mathrm{electron}}\\otimes\\begin{array}{c}\n",
    "\\Pi_{u}(x,y)\\\\\n",
    "\\Sigma_{u}^{+}(z)\n",
    "\\end{array} & \\supseteq & \\Sigma_{g}^{+}\\\\\n",
    "\\Gamma_{\\mathrm{electron}}\\otimes\\begin{array}{c}\n",
    "\\Pi_{u}(x,y)\\\\\n",
    "\\Sigma_{u}^{+}(z)\n",
    "\\end{array} & \\supseteq & \\Sigma_{g}^{+}\n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence:\n",
    "\n",
    "\\begin{equation}\n",
    "\\Gamma_{\\mathrm{electron}}=\\begin{array}{c}\n",
    "\\Pi_{u}(x,y)\\\\\n",
    "\\Sigma_{u}^{+}(z)\n",
    "\\end{array}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we also need to specify the total scattering symmetry $\\Gamma_{\\mathrm{scat}}=\\Gamma_{\\mathrm{ion}}\\otimes\\Gamma_{\\mathrm{electron}}$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\Gamma_{\\mathrm{scat}}=\\Sigma_{g}^{+}\\otimes\\begin{array}{c}\n",
    "\\Pi_{u}(x,y)\\\\\n",
    "\\Sigma_{u}^{+}(z)\n",
    "\\end{array}=\\begin{array}{c}\n",
    "\\Pi_{u}(x,y)\\\\\n",
    "\\Sigma_{u}^{+}(z)\n",
    "\\end{array}\n",
    "\\end{equation}\n",
    "\n",
    "... which is identical to $\\Gamma_{\\mathrm{electron}}$ in this simple case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These symmetries correspond to the SU and PU cases set in the ePS input earlier. (A [full list of ePS supported symmetries is given in the manual](https://epolyscat.droppages.com/SymmetryLabels).) It is worth noting that the continua correspond to the polarisation of the electric field in the molecular frame, hence which (Cartesian) dipole component is selected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: NO2 with autogeneration\n",
    "\n",
    "(To follow - part of epsman development...!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multi-energy calculations\n",
    "\n",
    "The most common case (for photoionization problems) is likely to be setting multiple energy points for a calculation. An example is shown in `test04` ([full output](https://epolyscat.droppages.com/Tests_test04_ostnd)), which basically uses the `ScatEng` record:\n",
    "\n",
    "```\n",
    "ScatEng 0.5 10.0 15.0      # list of scattering energies\n",
    "```\n",
    "\n",
    "In this case, later use of `GetCro` will output properties at the energies set here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ePolyScat: advanced usage\n",
    "\n",
    "See pt. 2 notebook!\n",
    "\n",
    "- Working directly with matrix elements.\n",
    "- Post-processing with ePSproc."
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}