ePSproc demo¶
Basic IO¶
[1]:
import sys
import os
import numpy as np
# For module testing, include path to module here
modPath = r'D:\code\github\ePSproc'
sys.path.append(modPath)
import epsproc as ep
* pyevtk not found, VTK export not available.
[2]:
# Load data from modPath\data
dataPath = os.path.join(modPath, 'data')
# Scan data dir
dataSet = ep.readMatEle(fileBase = dataPath)
*** ePSproc readMatEle(): scanning files for DumpIdy segments (matrix elements)
*** Scanning dir
D:\code\github\ePSproc\data
Found 2 .out file(s)
*** Reading ePS output file: D:\code\github\ePSproc\data\n2_3sg_0.1-50.1eV_A2.inp.out
Expecting 51 energy points.
Expecting 2 symmetries.
Expecting 102 dumpIdy segments.
Found 102 dumpIdy segments (sets of matrix elements).
Processing segments to Xarrays...
Processed 102 sets of matrix elements (0 blank)
*** Reading ePS output file: D:\code\github\ePSproc\data\no2_demo_ePS.out
Expecting 1 energy points.
Expecting 3 symmetries.
Expecting 3 dumpIdy segments.
Found 3 dumpIdy segments (sets of matrix elements).
Processing segments to Xarrays...
Processed 3 sets of matrix elements (0 blank)
Structure¶
Data is read and sorted into Xarrays, currently one Xarray per input file. The full dimensionality is maintained here.
Calling the array will provide some output…
[3]:
dataSet[1]
[3]:
<xarray.DataArray 'no2_demo_ePS.out' (LM: 110, Eke: 1, Sym: 3, mu: 3, it: 1, Type: 2)>
array([[[[[[ nan +nanj, nan +nanj]],
...,
[[ nan +nanj, nan +nanj]]],
...,
[[[-0.006893+0.212752j, -0.002009+0.126392j]],
...,
[[-0.006893+0.212752j, -0.002009+0.126392j]]]]],
...,
[[[[[ nan +nanj, nan +nanj]],
...,
[[ nan +nanj, nan +nanj]]],
...,
[[[ nan +nanj, nan +nanj]],
...,
[[ nan +nanj, nan +nanj]]]]]])
Coordinates:
* LM (LM) MultiIndex
- l (LM) int64 1 1 2 2 2 2 3 3 3 3 3 ... 10 10 10 10 10 10 10 10 10 10
- m (LM) int64 -1 1 -2 -1 1 2 -3 -2 -1 1 2 ... -1 1 2 3 4 5 6 7 8 9 10
* mu (mu) int64 -1 0 1
* it (it) int64 1
* Type (Type) object 'L' 'V'
* Eke (Eke) float64 0.81
* Sym (Sym) MultiIndex
- Cont (Sym) object 'A2' 'B2' 'B1'
- Targ (Sym) object 'A2' 'A2' 'A2'
- Total (Sym) object 'A1' 'B1' 'B2'
Attributes:
E: 0.81
Ehv: 14.402
SF: (2.0663304+3.9041597j)
Lmax: 12
Targ: A2
QNs: ['m', 'l', 'mu', 'ip', 'it', 'Value']
file: no2_demo_ePS.out
fileBase: D:\code\github\ePSproc\data
… and sub-selection can provide sets of matrix elements as a function of energy, symmetry and type.
[4]:
inds = {'Type':'L','Cont':'A2','mu':0}
dataSet[1].sel(inds).squeeze()
[4]:
<xarray.DataArray 'no2_demo_ePS.out' (LM: 110)>
array([ nan +nanj, nan +nanj,
-1.197776e-01-1.466388e-02j, nan +nanj,
nan +nanj, 1.197776e-01+1.466388e-02j,
nan +nanj, -7.045917e-01+5.263324e-01j,
nan +nanj, nan +nanj,
7.045917e-01-5.263324e-01j, nan +nanj,
1.411549e-03+5.641442e-03j, nan +nanj,
1.445839e-02-3.747604e-02j, nan +nanj,
nan +nanj, -1.445839e-02+3.747604e-02j,
nan +nanj, -1.411549e-03-5.641442e-03j,
nan +nanj, -1.026345e-02-9.266936e-03j,
nan +nanj, -4.208815e-03-3.131550e-03j,
nan +nanj, nan +nanj,
4.208815e-03+3.131550e-03j, nan +nanj,
1.026345e-02+9.266936e-03j, nan +nanj,
-5.931944e-05-1.768185e-06j, nan +nanj,
4.142273e-04+8.937456e-05j, nan +nanj,
2.414040e-04+8.185349e-05j, nan +nanj,
nan +nanj, -2.414040e-04-8.185349e-05j,
nan +nanj, -4.142273e-04-8.937456e-05j,
nan +nanj, 5.931944e-05+1.768185e-06j,
nan +nanj, 5.544370e-05-6.463095e-05j,
nan +nanj, 2.389949e-05-2.418784e-05j,
nan +nanj, 3.925946e-06-8.804978e-07j,
nan +nanj, nan +nanj,
-3.925946e-06+8.804978e-07j, nan +nanj,
-2.389949e-05+2.418784e-05j, nan +nanj,
-5.544370e-05+6.463095e-05j, nan +nanj,
-6.068016e-08-1.736560e-07j, nan +nanj,
2.439349e-06+5.382094e-06j, nan +nanj,
3.867429e-06+7.654976e-06j, nan +nanj,
2.323386e-06+5.515427e-06j, nan +nanj,
nan +nanj, -2.323386e-06-5.515427e-06j,
nan +nanj, -3.867429e-06-7.654976e-06j,
nan +nanj, -2.439349e-06-5.382094e-06j,
nan +nanj, 6.068016e-08+1.736560e-07j,
nan +nanj, 3.623968e-08-1.024486e-07j,
nan +nanj, 1.568559e-07-1.733682e-07j,
nan +nanj, 9.961257e-08-4.818728e-08j,
nan +nanj, -3.466612e-08+1.962320e-08j,
nan +nanj, nan +nanj,
3.466612e-08-1.962320e-08j, nan +nanj,
-9.961257e-08+4.818728e-08j, nan +nanj,
-1.568559e-07+1.733682e-07j, nan +nanj,
-3.623968e-08+1.024486e-07j, nan +nanj,
4.439372e-09-1.080641e-08j, nan +nanj,
7.202754e-09-1.264103e-08j, nan +nanj,
1.078749e-08-1.425429e-08j, nan +nanj,
1.810095e-08+7.838436e-09j, nan +nanj,
1.051935e-08+1.015762e-08j, nan +nanj,
nan +nanj, -1.051935e-08-1.015762e-08j,
nan +nanj, -1.810095e-08-7.838436e-09j,
nan +nanj, -1.078749e-08+1.425429e-08j,
nan +nanj, -7.202754e-09+1.264103e-08j,
nan +nanj, -4.439372e-09+1.080641e-08j])
Coordinates:
* LM (LM) MultiIndex
- l (LM) int64 1 1 2 2 2 2 3 3 3 3 3 ... 10 10 10 10 10 10 10 10 10 10
- m (LM) int64 -1 1 -2 -1 1 2 -3 -2 -1 1 2 ... -1 1 2 3 4 5 6 7 8 9 10
mu int64 0
it int64 1
Type <U1 'L'
Eke float64 0.81
Sym object ('A2', 'A1')
Attributes:
E: 0.81
Ehv: 14.402
SF: (2.0663304+3.9041597j)
Lmax: 12
Targ: A2
QNs: ['m', 'l', 'mu', 'ip', 'it', 'Value']
file: no2_demo_ePS.out
fileBase: D:\code\github\ePSproc\data
The matEleSelector function does the same thing, and also includes thresholding on abs values:
[5]:
# Set sq = True to squeeze on singleton dimensions
ep.matEleSelector(dataSet[1], thres=1e-2, inds = inds, sq = True)
[5]:
<xarray.DataArray 'no2_demo_ePS.out' (LM: 8)>
array([-0.119778-0.014664j, 0.119778+0.014664j, -0.704592+0.526332j,
0.704592-0.526332j, 0.014458-0.037476j, -0.014458+0.037476j,
-0.010263-0.009267j, 0.010263+0.009267j])
Coordinates:
* LM (LM) MultiIndex
- l (LM) int64 2 2 3 3 4 4 5 5
- m (LM) int64 -2 2 -2 2 -2 2 -4 4
mu int64 0
it int64 1
Type <U1 'L'
Eke float64 0.81
Sym object ('A2', 'A1')
Attributes:
E: 0.81
Ehv: 14.402
SF: (2.0663304+3.9041597j)
Lmax: 12
Targ: A2
QNs: ['m', 'l', 'mu', 'ip', 'it', 'Value']
file: no2_demo_ePS.out
fileBase: D:\code\github\ePSproc\data
Basic plotting from Xarray¶
[6]:
# Plot matrix elements using Xarray functionality
daPlot = dataSet[0].sum('mu').sum('Sym').sel({'Type':'L'}).squeeze()
daPlot.pipe(np.abs).plot.line(x='Eke')
[6]:
[<matplotlib.lines.Line2D at 0x13658d327f0>,
<matplotlib.lines.Line2D at 0x13658d32588>,
<matplotlib.lines.Line2D at 0x13658d32780>,
<matplotlib.lines.Line2D at 0x13658d32f28>,
<matplotlib.lines.Line2D at 0x13658d322b0>,
<matplotlib.lines.Line2D at 0x13658d32320>,
<matplotlib.lines.Line2D at 0x13658d32908>,
<matplotlib.lines.Line2D at 0x13658d32518>,
<matplotlib.lines.Line2D at 0x136597cfda0>,
<matplotlib.lines.Line2D at 0x136597cf588>,
<matplotlib.lines.Line2D at 0x13659b91978>,
<matplotlib.lines.Line2D at 0x136597cfc18>,
<matplotlib.lines.Line2D at 0x136597cfef0>,
<matplotlib.lines.Line2D at 0x136597cffd0>,
<matplotlib.lines.Line2D at 0x136597cf400>,
<matplotlib.lines.Line2D at 0x136597cfc50>,
<matplotlib.lines.Line2D at 0x136597cfbe0>,
<matplotlib.lines.Line2D at 0x136597cfa90>]
[7]:
# Plot with faceting on type
daPlot = dataSet[0].sum('mu').sum('Sym').squeeze()
daPlot.pipe(np.abs).plot.line(x='Eke', col='Type')
# Plot with faceting on symmetry
daPlot = dataSet[0].sum('mu').squeeze()
daPlot.pipe(np.abs).plot.line(x='Eke', col='Sym', row='Type')
[7]:
<xarray.plot.facetgrid.FacetGrid at 0x136588d90f0>
Calculate MFPADs¶
Calculate MFPADs, as given by:
\begin{equation}
I_{\mu_{0}}(\theta_{\hat{k}},\phi_{\hat{k}},\theta_{\hat{n}},\phi_{\hat{n}})=\frac{4\pi^{2}E}{cg_{p_{i}}}\sum_{\mu_{i},\mu_{f}}|T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\theta_{\hat{k}},\phi_{\hat{k}},\theta_{\hat{n}},\phi_{\hat{n}})|^{2}\label{eq:MFPAD}
\end{equation}
\begin{equation}
T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\theta_{\hat{k}},\phi_{\hat{k}},\theta_{\hat{n}},\phi_{\hat{n}})=\sum_{l,m,\mu}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)Y_{lm}^{*}(\theta_{\hat{k}},\phi_{\hat{k}})D_{\mu,\mu_{0}}^{1}(R_{\hat{n}})\label{eq:TMF}
\end{equation}
\begin{equation}
I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)=\langle\Psi_{i}^{p_{i},\mu_{i}}|\hat{d_{\mu}}|\Psi_{f}^{p_{f},\mu_{f}}\varphi_{klm}^{(-)}\rangle\label{eq:I}
\end{equation}
In this formalism:
- \(I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)\) is the radial part of the dipole matrix element, determined from the initial and final state electronic wavefunctions \(\Psi_{i}^{p_{i},\mu_{i}}\)and \(\Psi_{f}^{p_{f},\mu_{f}}\), photoelectron wavefunction \(\varphi_{klm}^{(-)}\) and dipole operator \(\hat{d_{\mu}}\). Here the wavefunctions are indexed by irreducible representation (i.e. symmetry) by the labels \(p_{i}\) and \(p_{f}\), with components \(\mu_{i}\) and \(\mu_{f}\) respectively; \(l,m\) are angular momentum components, \(\mu\) is the projection of the polarization into the MF (from a value \(\mu_{0}\) in the LF). Each energy and irreducible representation corresponds to a calculation in ePolyScat.
- \(T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\theta_{\hat{k}},\phi_{\hat{k}},\theta_{\hat{n}},\phi_{\hat{n}})\) is the full matrix element (expanded in polar coordinates) in the MF, where \(\hat{k}\) denotes the direction of the photoelectron \(\mathbf{k}\)-vector, and \(\hat{n}\) the direction of the polarization vector \(\mathbf{n}\) of the ionizing light. Note that the summation over components \(\{l,m,\mu\}\) is coherent, and hence phase sensitive.
- \(Y_{lm}^{*}(\theta_{\hat{k}},\phi_{\hat{k}})\) is a spherical harmonic.
- \(D_{\mu,\mu_{0}}^{1}(R_{\hat{n}})\) is a Wigner rotation matrix element, with a set of Euler angles \(R_{\hat{n}}=(\phi_{\hat{n}},\theta_{\hat{n}},\chi_{\hat{n}})\), which rotates/projects the polarization into the MF .
- \(I_{\mu_{0}}(\theta_{\hat{k}},\phi_{\hat{k}},\theta_{\hat{n}},\phi_{\hat{n}})\) is the final (observable) MFPAD, for a polarization \(\mu_{0}\) and summed over all symmetry components of the initial and final states, \(\mu_{i}\) and \(\mu_{f}\). Note that this sum can be expressed as an incoherent summation, since these components are (by definition) orthogonal.
- \(g_{p_{i}}\) is the degeneracy of the state \(p_{i}\).
See: Toffoli, D., Lucchese, R. R., Lebech, M., Houver, J. C., & Dowek, D. (2007). Molecular frame and recoil frame photoelectron angular distributions from dissociative photoionization of NO2. The Journal of Chemical Physics, 126(5), 054307. https://doi.org/10.1063/1.2432124
[8]:
print('MFPADs for test NO2 dataset (single energy, (z,x,y) pol states)')
TX, TlmX = ep.mfpad(dataSet[1])
# Plot for each pol geom (symmetry)
for n in range(0,3):
ep.sphSumPlotX(TX[n].sum('Sym').squeeze(), pType = 'a')
MFPADs for test NO2 dataset (single energy, (z,x,y) pol states)
[9]:
# Plot abs(TX) images using Xarray functionality
TX.squeeze().pipe(np.abs).plot(x='Theta',y='Phi', col='Euler', row='Sym')
[9]:
<xarray.plot.facetgrid.FacetGrid at 0x1365e8b7710>
[10]:
# Plot MFPAD surfaces vs E
print('N2 test data, abs(TX) vs E and (z,x,y) pol geom')
TX, TlmX = ep.mfpad(dataSet[0])
TXplot = TX.sum('Sym').squeeze().isel(Eke=slice(0,50,10))
TXplot.pipe(np.abs).plot(x='Theta',y='Phi', row='Euler', col='Eke')
N2 test data, abs(TX) vs E and (z,x,y) pol geom
[10]:
<xarray.plot.facetgrid.FacetGrid at 0x1365d5980f0>
[11]:
# Try Plotly with looping functionality... this gives 3D interactive surf plots.
# Note this is currently set to expect 3D data only, and loop over 3rd dim.
# This is a work in progress...!
TX, TlmX = ep.mfpad(dataSet[0])
TXplot = TX.sum('Sym').squeeze().isel(Eke=slice(0,50,10))
ep.mfpadPlotPL(np.abs(TXplot[0]), rc = [1,5])