Symmetrized harmonics demo pt III: Working with symmetry, direct products and photoionization selection rules

21/04/23

This notebook demonstrates functionality for working with symmetry, making use of libmsym for the core symmetry functionality.

For more background, see the Symmetriezed Haromnics demo.

For more on point groups and symmetry, see http://www.gernot-katzers-spice-pages.com/character_tables/index.html, which includes full character tables and associated data.

Computational refs:

  • Boiko, D.L., Féron, P. and Besnard, P. (2006) ‘Simple method to obtain symmetry harmonics of point groups’, The Journal of Chemical Physics, 125(9), p. 094104. doi:10.1063/1.2338040.

  • Johansson, M. and Veryazov, V. (2017) ‘Automatic procedure for generating symmetry adapted wavefunctions’, Journal of Cheminformatics, 9(1), p. 8. doi:10.1186/s13321-017-0193-3. (Manuscript for libmsym library.)

Generating basis sets and direct products

With the symHarm class, a basis set of symmetrized harmonics can be created.

As of commit d470855dae78245257a960a0306bbd85ddcd0160, the class also computes direct products, including full direct product tables and dipole-allowed terms.

Direct products are computed as array multiplications and vector projections from the libmsym objects, see https://github.com/mcodev31/libmsym/issues/29 for the basic method, and special thanks to mcodev31 for the library and original methods.

[1]:

# Import class and generate basis
from pemtk.sym.symHarm import symHarm

sym = 'D10h'
symBasis = symHarm(PG = sym, lmax=4)

OMP: Info #273: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

* sparse not found, sparse matrix forms not available.
* natsort not found, some sorting functions not available.

* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.
* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available.
*** Mapping coeffs to ePSproc dataType = matE
Remapped dims: {'C': 'Cont', 'mu': 'it'}
Added dim Eke
Added dim Targ
Added dim Total
Added dim mu
Added dim Type
Found dipole symmetries:
{'A2u': {'m': [0], 'pol': ['z']}, 'E1u': {'m': [-1, 1, -1, 1], 'pol': ['x', 'y']}}

[2]:

# The full direct product table is stored in the class, as a Pandas DataFrame
symBasis.directProductTable

[2]:
A1g A2g A1u A2u B1g B1u B2g B2u E1u E1g E2g E2u E3u E3g E4g E4u
A1g [A1g] [A2g] [A1u] [A2u] [B1g] [B1u] [B2g] [B2u] [E1u] [E1g] [E2g] [E2u] [E3u] [E3g] [E4g] [E4u]
A2g [A2g] [A1g] [A2u] [A1u] [B2g] [B2u] [B1g] [B1u] [E1u] [E1g] [E2g] [E2u] [E3u] [E3g] [E4g] [E4u]
A1u [A1u] [A2u] [A1g] [A2g] [B1u] [B1g] [B2u] [B2g] [E1g] [E1u] [E2u] [E2g] [E3g] [E3u] [E4u] [E4g]
A2u [A2u] [A1u] [A2g] [A1g] [B2u] [B2g] [B1u] [B1g] [E1g] [E1u] [E2u] [E2g] [E3g] [E3u] [E4u] [E4g]
B1g [B1g] [B2g] [B1u] [B2u] [A1g] [A1u] [A2g] [A2u] [E4u] [E4g] [E3g] [E3u] [E2u] [E2g] [E1g] [E1u]
B1u [B1u] [B2u] [B1g] [B2g] [A1u] [A1g] [A2u] [A2g] [E4g] [E4u] [E3u] [E3g] [E2g] [E2u] [E1u] [E1g]
B2g [B2g] [B1g] [B2u] [B1u] [A2g] [A2u] [A1g] [A1u] [E4u] [E4g] [E3g] [E3u] [E2u] [E2g] [E1g] [E1u]
B2u [B2u] [B1u] [B2g] [B1g] [A2u] [A2g] [A1u] [A1g] [E4g] [E4u] [E3u] [E3g] [E2g] [E2u] [E1u] [E1g]
E1u [E1u] [E1u] [E1g] [E1g] [E4u] [E4g] [E4u] [E4g] [A1g, A2g, E2g] [A1u, A2u, E2u] [E1u, E3u] [E1g, E3g] [E2g, E4g] [E2u, E4u] [B1u, B2u, E3u] [B1g, B2g, E3g]
E1g [E1g] [E1g] [E1u] [E1u] [E4g] [E4u] [E4g] [E4u] [A1u, A2u, E2u] [A1g, A2g, E2g] [E1g, E3g] [E1u, E3u] [E2u, E4u] [E2g, E4g] [B1g, B2g, E3g] [B1u, B2u, E3u]
E2g [E2g] [E2g] [E2u] [E2u] [E3g] [E3u] [E3g] [E3u] [E1u, E3u] [E1g, E3g] [A1g, A2g, E4g] [A1u, A2u, E4u] [B1u, B2u, E1u] [B1g, B2g, E1g] [E2g, E4g] [E2u, E4u]
E2u [E2u] [E2u] [E2g] [E2g] [E3u] [E3g] [E3u] [E3g] [E1g, E3g] [E1u, E3u] [A1u, A2u, E4u] [A1g, A2g, E4g] [B1g, B2g, E1g] [B1u, B2u, E1u] [E2u, E4u] [E2g, E4g]
E3u [E3u] [E3u] [E3g] [E3g] [E2u] [E2g] [E2u] [E2g] [E2g, E4g] [E2u, E4u] [B1u, B2u, E1u] [B1g, B2g, E1g] [A1g, A2g, E4g] [A1u, A2u, E4u] [E1u, E3u] [E1g, E3g]
E3g [E3g] [E3g] [E3u] [E3u] [E2g] [E2u] [E2g] [E2u] [E2u, E4u] [E2g, E4g] [B1g, B2g, E1g] [B1u, B2u, E1u] [A1u, A2u, E4u] [A1g, A2g, E4g] [E1g, E3g] [E1u, E3u]
E4g [E4g] [E4g] [E4u] [E4u] [E1g] [E1u] [E1g] [E1u] [B1u, B2u, E3u] [B1g, B2g, E3g] [E2g, E4g] [E2u, E4u] [E1u, E3u] [E1g, E3g] [A1g, A2g, E2g] [A1u, A2u, E2u]
E4u [E4u] [E4u] [E4g] [E4g] [E1u] [E1g] [E1u] [E1g] [B1g, B2g, E3g] [B1u, B2u, E3u] [E2u, E4u] [E2g, E4g] [E1g, E3g] [E1u, E3u] [A1u, A2u, E2u] [A1g, A2g, E2g] 
[3]:

# Direct products can be determined for pairs of irreps from the table elements via a Pandas lookup
print(symBasis.directProductTable.loc['A1g','B1u'])
print(symBasis.directProductTable.loc['E1u','E4u'])

['B1u']
['B1g', 'B2g', 'E3g']

[4]:

# Arb direct products are also computable with the directProduct method - this computes the direct product of all the terms passed
print(symBasis.directProduct(['A1g','B1u']))
print(symBasis.directProduct(['E1u','E4u']))
print(symBasis.directProduct(['A1g','A1u','B1u', 'B1u']))
print(symBasis.directProduct(['A1g','A1u','B1u', 'B1u','E1u','E4u']))

['B1u']
['B1g', 'B2g', 'E3g']
['A1u']
['B1u', 'B2u', 'E3u']

Dipole terms

For a given point group, the allow dipole terms correspond to \(l=1\). Futher, \(l=1, m=0\) corresponds to the Cartesian \(z\) component, and linear combinations of \(m=\pm1\) the \(x,y\) components. Depending on the point group the \(x,y\) terms may have the same symmetry (doubly-degenerate), or may not.

[5]:

# Some corresponding outputs are in the dipole dictionary
symBasis.dipole.keys()

[5]:

dict_keys(['dipoleSyms', 'dipolePD', 'dipoleXR', 'coeffsXR', 'dipContList', 'prodTerms', 'prodValid'])

[6]:

# Allowed terms and mappings are given in 'dipoleSyms'
symBasis.dipole['dipoleSyms']

[6]:

{'A2u': {'m': [0], 'pol': ['z']},
 'E1u': {'m': [-1, 1, -1, 1], 'pol': ['x', 'y']}}

[7]:

# A subset of the direct product table for only the dipole terms is given in 'prodTerms'
symBasis.dipole['prodTerms']

[7]:
A2u E1u
A1g [A2u] [E1u]
A2g [A1u] [E1u]
A1u [A2g] [E1g]
A2u [A1g] [E1g]
B1g [B2u] [E4u]
B1u [B2g] [E4g]
B2g [B1u] [E4u]
B2u [B1g] [E4g]
E1u [E1g] [A1g, A2g, E2g]
E1g [E1u] [A1u, A2u, E2u]
E2g [E2u] [E1u, E3u]
E2u [E2g] [E1g, E3g]
E3u [E3g] [E2g, E4g]
E3g [E3u] [E2u, E4u]
E4g [E4u] [B1u, B2u, E3u]
E4u [E4g] [B1g, B2g, E3g]

Mapping allowed terms to photoionization matrix elements

With knowledge of the dipole terms and direct products, the allowed set of (symmetrized) matrix elements for a given case can be determined…

Symmetry selection rules in photoionization

Essentially, the relevant direct products must contain the totally symmetric representation of the point group to constitute an allowed combination:

\begin{eqnarray} \Gamma_{\mathrm{ion}}\otimes\Gamma_{\mathrm{electron}}\otimes\Gamma_{\mathrm{dipole}}\otimes\Gamma_{\mathrm{neutral}} & \supseteq & \Gamma_{\mathrm{symmetric}}\\ \end{eqnarray}

Where \(\Gamma\) is the character of the quantity of interest, and \(\Gamma_{\mathrm{symmetric}}\) denotes the totally symmetric irrep in the group, e.g. \(\mathrm{A_{1}}\). 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.

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

Worked example: \(N_2~3\sigma_g^{-1}\) ionization

By hand

To do this by hand, we’ll need the character table for \(D_{\infty h}\) (see also Oxford materials page from Atkins, Child & Philips, which includes the direct product tables (PDF version).

For this example, \(\Gamma_{\mathrm{ion}} = \Gamma_{\mathrm{neutral}} = \Gamma_{\mathrm{symmetric}} = \Sigma_{g}^{+}\). The dipole symmetries correspond to the Cartesian axes/translations given in the character tables, hence we find \(\Sigma_{u}^{+} = (z)\) and \(\Pi_{u} = (x,y)\)

Then just plug in and work through to determine the allowed continuum symmetries…

\begin{eqnarray} \Gamma_{\mathrm{ion}}\otimes\Gamma_{\mathrm{electron}}\otimes\Gamma_{\mathrm{dipole}}\otimes\Gamma_{\mathrm{neutral}} & \supseteq & \Sigma_{g}^{+}\\ \Sigma_{g}^{+}\otimes\Gamma_{\mathrm{electron}}\otimes\begin{array}{c} \Pi_{u}(x,y)\\ \Sigma_{u}^{+}(z) \end{array}\otimes\Sigma_{g}^{+} & \supseteq & \Sigma_{g}^{+}\\ \Sigma_{g}^{+}\otimes\Gamma_{\mathrm{electron}}\otimes\begin{array}{c} \Pi_{u}(x,y)\\ \Sigma_{u}^{+}(z) \end{array} & \supseteq & \Sigma_{g}^{+}\\ \Gamma_{\mathrm{electron}}\otimes\begin{array}{c} \Pi_{u}(x,y)\\ \Sigma_{u}^{+}(z) \end{array} & \supseteq & \Sigma_{g}^{+} \end{eqnarray}

Hence:

\begin{equation} \Gamma_{\mathrm{electron}}=\begin{array}{c} \Pi_{u}(x,y)\\ \Sigma_{u}^{+}(z) \end{array} \end{equation}

The total scattering symmetry is given by \(\Gamma_{\mathrm{scat}}=\Gamma_{\mathrm{ion}}\otimes\Gamma_{\mathrm{electron}}\):

\begin{equation} \Gamma_{\mathrm{scat}}=\Sigma_{g}^{+}\otimes\begin{array}{c} \Pi_{u}(x,y)\\ \Sigma_{u}^{+}(z) \end{array}=\begin{array}{c} \Pi_{u}(x,y)\\ \Sigma_{u}^{+}(z) \end{array} \end{equation}

… which is identical to \(\Gamma_{\mathrm{electron}}\) in this simple case.

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.) 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.

Computationally

The cylindrically-symmetric \(\infty\) groups can be approximated by a high-order group, e.g. \(D_{\infty h} \approx D_{10h}\), as set above. (For a cross-check, see the full tables and direct products online.) Here the notational convention is \(A = \Sigma\), \(E1 = \Pi\), \(E2 = \Delta\) and so forth (see the :math:`D_{infty h} character table <http://symmetry.jacobs-university.de/cgi-bin/group.cgi?group=1001&option=4>`__ for more details).

A table of values \(\Gamma_{x,y,z}\otimes\Gamma\) are pre-computed for the point group and set to self.dipole['prodTerms'], as shown above. The variable self.dipole['prodValid'] additionally gives the cases only fulfilling the condition \(\Gamma_{product} \supseteq \Gamma_{sym}\), this may be useful in some simple cases, such as the example here.

[8]:

# A subset of the direct product table for only the dipole terms, and if the results contain the totally symmetric representation, is given in 'prodValid'
symBasis.dipole['prodValid']

[8]:
A2u E1u
A2u True
E1u True

For the basic case here, this provides a quick look-up, with valid cases \(A2u\) and \(E1u\) corresponding to the allowed electron terms. Note, however, that this is only the result for a direct product from the pair shown, so for other choices of \(\Gamma_{\mathrm{ion}}\) and \(\Gamma_{\mathrm{neutral}}\) these values may not correspond to the desired results.

Worked example: \(N_2~3\pi_u^{-1}\) ionization

For more complicated cases,

  1. the full set of terms can be specified, and the dipole-product terms computed as defined by the point group with directProductDipole.

  2. Alternatively, the direct products for all possible continua can be determined using directProductContinuum.

Both cases are illustrated below for the choice of \(\Gamma_{\mathrm{ion}} = \Pi_u\), \(\Gamma_{\mathrm{neutral}} = \Sigma_{g}^{+}\), using \(D_{10 h}\). As previously, allowed terms are defined to be those fulfilling \(\Gamma_{product} \supseteq \Gamma_{sym}\).

[9]:

# Specify full terms - this will indicate if a given selection fulfils the dipole selection rules for the dipole terms in the given point group.
sNeutral = 'A1g'
sIon = 'E1u'
sContinuum = 'A1g'

# The `directProductDipole` method outputs results with the dipole symmetries included
symBasis.directProductDipole([sNeutral, sIon, sContinuum])
symBasis.dipole['dipoleProducts']

[9]:

{'A2u': {'m': [0],
  'pol': ['z'],
  'dipSym': 'A2u',
  'result': ['E1g'],
  'allowed': False,
  'terms': ['A1g', 'E1u', 'A1g']},
 'E1u': {'m': [-1, 1, -1, 1],
  'pol': ['x', 'y'],
  'dipSym': 'E1u',
  'result': ['A1g', 'A2g', 'E2g'],
  'allowed': True,
  'terms': ['A1g', 'E1u', 'A1g']}}

For this case we see that the set of states [sNeutral, sIon, sContinuum] is allowed for the \(E_{1g}(x,y)\) dipole component.

[10]:

# For cases where the continuum symmetry is not known/specified, the `directProductContinuum` method can be used, and will output results looped over all irreps
symBasis.directProductContinuum([sNeutral, sIon])
    allowed m pol result terms
Dipole Target          
A2u A1g 0 [0] ['z'] ['E1g'] ['A1g', 'E1u']
A2g 0 [0] ['z'] ['E1g'] ['A1g', 'E1u']
A1u 0 [0] ['z'] ['E1u'] ['A1g', 'E1u']
A2u 0 [0] ['z'] ['E1u'] ['A1g', 'E1u']
B1g 0 [0] ['z'] ['E4g'] ['A1g', 'E1u']
B1u 0 [0] ['z'] ['E4u'] ['A1g', 'E1u']
B2g 0 [0] ['z'] ['E4g'] ['A1g', 'E1u']
B2u 0 [0] ['z'] ['E4u'] ['A1g', 'E1u']
E1u 0 [0] ['z'] ['A1u', 'A2u', 'E2u'] ['A1g', 'E1u']
E1g 1 [0] ['z'] ['A1g', 'A2g', 'E2g'] ['A1g', 'E1u']
E2g 0 [0] ['z'] ['E1g', 'E3g'] ['A1g', 'E1u']
E2u 0 [0] ['z'] ['E1u', 'E3u'] ['A1g', 'E1u']
E3u 0 [0] ['z'] ['E2u', 'E4u'] ['A1g', 'E1u']
E3g 0 [0] ['z'] ['E2g', 'E4g'] ['A1g', 'E1u']
E4g 0 [0] ['z'] ['B1g', 'B2g', 'E3g'] ['A1g', 'E1u']
E4u 0 [0] ['z'] ['B1u', 'B2u', 'E3u'] ['A1g', 'E1u']
E1u A1g 1 [-1, 1, -1, 1] ['x', 'y'] ['A1g', 'A2g', 'E2g'] ['A1g', 'E1u']
A2g 1 [-1, 1, -1, 1] ['x', 'y'] ['A1g', 'A2g', 'E2g'] ['A1g', 'E1u']
A1u 0 [-1, 1, -1, 1] ['x', 'y'] ['A1u', 'A2u', 'E2u'] ['A1g', 'E1u']
A2u 0 [-1, 1, -1, 1] ['x', 'y'] ['A1u', 'A2u', 'E2u'] ['A1g', 'E1u']
B1g 0 [-1, 1, -1, 1] ['x', 'y'] ['B1g', 'B2g', 'E3g'] ['A1g', 'E1u']
B1u 0 [-1, 1, -1, 1] ['x', 'y'] ['B1u', 'B2u', 'E3u'] ['A1g', 'E1u']
B2g 0 [-1, 1, -1, 1] ['x', 'y'] ['B1g', 'B2g', 'E3g'] ['A1g', 'E1u']
B2u 0 [-1, 1, -1, 1] ['x', 'y'] ['B1u', 'B2u', 'E3u'] ['A1g', 'E1u']
E1u 0 [-1, 1, -1, 1] ['x', 'y'] ['E1u', 'E3u'] ['A1g', 'E1u']
E1g 0 [-1, 1, -1, 1] ['x', 'y'] ['E1g', 'E3g'] ['A1g', 'E1u']
E2g 1 [-1, 1, -1, 1] ['x', 'y'] ['A1g', 'A2g', 'E2g', 'E4g'] ['A1g', 'E1u']
E2u 0 [-1, 1, -1, 1] ['x', 'y'] ['A1u', 'A2u', 'E2u', 'E4u'] ['A1g', 'E1u']
E3u 0 [-1, 1, -1, 1] ['x', 'y'] ['B1u', 'B2u', 'E1u', 'E3u'] ['A1g', 'E1u']
E3g 0 [-1, 1, -1, 1] ['x', 'y'] ['B1g', 'B2g', 'E1g', 'E3g'] ['A1g', 'E1u']
E4g 0 [-1, 1, -1, 1] ['x', 'y'] ['E2g', 'E4g'] ['A1g', 'E1u']
E4u 0 [-1, 1, -1, 1] ['x', 'y'] ['E2u', 'E4u'] ['A1g', 'E1u'] 
[11]:

# Results are pushed to self.continuum, in dictionary and Pandas DataFrame formats, and can be manipulated using standard functionality.
# The subset of allowed values are also set to a separate DataFrame and list.
symBasis.continuum['allowed']['PD']

[11]:
allowed m pol result terms
Dipole Target
A2u E1g True [0] [z] [A1g, A2g, E2g] [A1g, E1u]
E1u A1g True [-1, 1, -1, 1] [x, y] [A1g, A2g, E2g] [A1g, E1u]
A2g True [-1, 1, -1, 1] [x, y] [A1g, A2g, E2g] [A1g, E1u]
E2g True [-1, 1, -1, 1] [x, y] [A1g, A2g, E2g, E4g] [A1g, E1u] 
[12]:

symBasis.continuum['allowed']['list']

[12]:

[('A2u', 'E1g'), ('E1u', 'A1g'), ('E1u', 'A2g'), ('E1u', 'E2g')]

In this case we find that the allowed symmetries are \(E1g\) only for the \(z\)-component, and \(A1g,A2g,E2g\) for the \((x,y)\)-component. In \(D_{\infty h}\) these correspond to \(\Pi_g\) and \((\Sigma_g^+,\Sigma_g^-,\Delta_g)\) continua respectively.

Checking allowed spherical harmonics

As detailed in the basic Symmetrized harmonics demo, the full basis can be displayed using the displayXlm method.

[13]:

# The full basis table
symBasis.displayXlm()   # Default case shows real harmonics

symBasis.displayXlm(YlmType='comp')  # Set for complex harmonics
b
l 0 1 2 3 4
Character ($\Gamma$) SALC (h) PFIX ($\mu$) m
A1g 0 0 0 1.0
1 0 0 1.0
2 0 0 1.0
A2u 0 0 0 1.0
1 0 0 1.0
E1g 0 0 -1 1.0
1 1 -1.0
1 0 -1 1.0
1 1 -1.0
E1u 0 0 1 1.0
1 -1 1.0
1 0 1 1.0
1 -1 1.0
E2g 0 0 2 1.0
1 -2 1.0
1 0 2 1.0
1 -2 1.0
E2u 0 0 -2 1.0
1 2 -1.0
E3g 0 0 -3 1.0
1 3 -1.0
E3u 0 0 3 1.0
1 -3 1.0
E4g 0 0 4 1.0
1 -4 1.0
b
l 0 1 2 3 4
Character ($\Gamma$) SALC (h) PFIX ($\mu$) m
A1g 0 0 0 (1+0j)
1 0 0 (1+0j)
2 0 0 (1+0j)
A2u 0 0 0 (1+0j)
1 0 0 (1+0j)
E1g 0 0 -1 (0.7071067811865475+0j)
1 (0.7071067811865475-0j)
1 -1 (-0.7071067811865475-0j)
1 (0.7071067811865475-0j)
1 0 -1 (0.7071067811865475+0j)
1 (0.7071067811865475-0j)
1 -1 (-0.7071067811865475-0j)
1 (0.7071067811865475-0j)
E1u 0 0 -1 (0.7071067811865475+0j)
1 (-0.7071067811865475-0j)
1 -1 (-0.7071067811865475+0j)
1 (-0.7071067811865475-0j)
1 0 -1 (0.7071067811865475+0j)
1 (-0.7071067811865475-0j)
1 -1 (-0.7071067811865475+0j)
1 (-0.7071067811865475-0j)
E2g 0 0 -2 (0.7071067811865475+0j)
2 (0.7071067811865475+0j)
1 -2 (-0.7071067811865475+0j)
2 (0.7071067811865475+0j)
1 0 -2 (0.7071067811865475+0j)
2 (0.7071067811865475+0j)
1 -2 (-0.7071067811865475+0j)
2 (0.7071067811865475+0j)
E2u 0 0 -2 (0.7071067811865475+0j)
2 (-0.7071067811865475+0j)
1 -2 (-0.7071067811865475-0j)
2 (-0.7071067811865475+0j)
E3g 0 0 -3 (0.7071067811865475+0j)
3 (0.7071067811865475-0j)
1 -3 (-0.7071067811865475-0j)
3 (0.7071067811865475-0j)
E3u 0 0 -3 (0.7071067811865475+0j)
3 (-0.7071067811865475-0j)
1 -3 (-0.7071067811865475+0j)
3 (-0.7071067811865475-0j)
E4g 0 0 -4 (0.7071067811865475+0j)
4 (0.7071067811865475+0j)
1 -4 (-0.7071067811865475+0j)
4 (0.7071067811865475+0j)

Additionally, if self.continuum is defined, these can be filtered. Currently this requires symFilter = True to be set, and will use the Target symmetry as a filter by default (or pass symFilterChannel = 'Channel' to use a different channel).

[14]:

# Basis table with the Character values limited to those defined in self.continuum['allowed']['PD'] Target column
symBasis.displayXlm(symFilter = 'Target')   # ['Character ($\Gamma$)','Target'])
b
l 0 1 2 3 4
Character ($\Gamma$) SALC (h) PFIX ($\mu$) m
A1g 0 0 0 1.0
1 0 0 1.0
2 0 0 1.0
E1g 0 0 -1 1.0
1 1 -1.0
1 0 -1 1.0
1 1 -1.0
E2g 0 0 2 1.0
1 -2 1.0
1 0 2 1.0
1 -2 1.0

In this case we find \(A1g\) corresponds to even-\(l\), \(m=0\), and the \(En\) components contain \(m=\pm n\) terms. Note, also, that there are no terms for \(A2g\) - even though it is allowed by the direct product, there are no liner combinations of spherical harmonics which span this irrep.

Defining allowed matrix elements

From the direct products the allowed matrix elements, in terms of spherical harmonics, can be defined.

Note that a copy of the basis can be assigned in ePSproc style matrix element format with the toEPSproc method (see Symmetrized harmonics demo pt. II for more info). Here the irreps for the point group are mapped to the continuum (free electron) symmetry label Cont, but no other symmetries are specified. For a specific ionization event, only subsets of these values are required, and also will be correlated with specific dipole terms, so some additional configuration is required.

[15]:

# Run conversion with a different dimMap & dataType - note this includes all symmetries, and both real and complex harmonic basis sets
dataType = 'matE'

# With custom dim mapping (optional)...
dimMap = {'C':'Cont', 'mu':'it'}   # Default dimMap = {'C':'Cont','h':'it', 'mu':'muX'}  # TBC - decide on default here.
# dimMap = {'C':'Cont','h':'it', 'mu':'muX'}  # Default case

symBasis.toePSproc(dataType=dataType, dimMap=dimMap)
symBasis.coeffs[dataType]

*** Mapping coeffs to ePSproc dataType = matE
Remapped dims: {'C': 'Cont', 'mu': 'it'}
Added dim Eke
Added dim Targ
Added dim Total
Added dim mu
Added dim Type

[15]:

<xarray.Dataset>
Dimensions:   (Type: 1, mu: 1, Eke: 1, h: 3, it: 2, LM: 45, Sym: 9)
Coordinates:
  * Type      (Type) <U1 'U'
  * mu        (mu) float64 nan
  * Eke       (Eke) int64 0
  * h         (h) int64 0 1 2
  * it        (it) int64 0 1
  * LM        (LM) MultiIndex
  - l         (LM) int64 0 0 0 0 0 0 0 0 0 1 1 1 1 ... 3 3 3 3 4 4 4 4 4 4 4 4 4
  - m         (LM) int64 -4 -3 -2 -1 0 1 2 3 4 -4 -3 ... 4 -4 -3 -2 -1 0 1 2 3 4
  * Sym       (Sym) MultiIndex
  - Cont      (Sym) object 'A1g' 'A2u' 'E1g' 'E1u' 'E2g' 'E2u' 'E3g' 'E3u' 'E4g'
  - Targ      (Sym) object 'U' 'U' 'U' 'U' 'U' 'U' 'U' 'U' 'U'
  - Total     (Sym) object 'U' 'U' 'U' 'U' 'U' 'U' 'U' 'U' 'U'
Data variables:
    b (real)  (Type, mu, Eke, h, it, LM, Sym) float64 nan nan nan ... nan nan
    b (comp)  (Type, mu, Eke, h, it, LM, Sym) complex128 (nan+nanj) ... (nan+...
Attributes:
    dataType:  matE
    name:      Symmetrized harmonics
    PG:        D10h
    lmax:      4
[16]:

# Tabulate results with Pandas conversion - complex components
import epsproc as ep
matEPD, _ = ep.multiDimXrToPD(symBasis.coeffs[dataType]['b (comp)'], colDims = 'Cont', thres=1e-4)
matEPD.fillna('')

[16]:
Cont A1g A2u E1g E1u E2g E2u E3g E3u E4g
Eke Targ Total Type h it l m mu
0 U U U 0 0 0 0 NaN (1+0j)
1 -1 NaN (0.7071067811865475+0j)
0 NaN (1+0j)
1 NaN (-0.7071067811865475-0j)
2 -2 NaN (0.7071067811865475+0j)
-1 NaN (0.7071067811865475+0j)
1 NaN (0.7071067811865475-0j)
2 NaN (0.7071067811865475+0j)
3 -3 NaN (0.7071067811865475+0j)
-2 NaN (0.7071067811865475+0j)
2 NaN (-0.7071067811865475+0j)
3 NaN (-0.7071067811865475-0j)
4 -4 NaN (0.7071067811865475+0j)
-3 NaN (0.7071067811865475+0j)
3 NaN (0.7071067811865475-0j)
4 NaN (0.7071067811865475+0j)
1 1 -1 NaN (-0.7071067811865475+0j)
1 NaN (-0.7071067811865475-0j)
2 -2 NaN (-0.7071067811865475+0j)
-1 NaN (-0.7071067811865475-0j)
1 NaN (0.7071067811865475-0j)
2 NaN (0.7071067811865475+0j)
3 -3 NaN (-0.7071067811865475+0j)
-2 NaN (-0.7071067811865475-0j)
2 NaN (-0.7071067811865475+0j)
3 NaN (-0.7071067811865475-0j)
4 -4 NaN (-0.7071067811865475+0j)
-3 NaN (-0.7071067811865475-0j)
3 NaN (0.7071067811865475-0j)
4 NaN (0.7071067811865475+0j)
1 0 2 0 NaN (1+0j)
3 -1 NaN (0.7071067811865475+0j)
0 NaN (1+0j)
1 NaN (-0.7071067811865475-0j)
4 -2 NaN (0.7071067811865475+0j)
-1 NaN (0.7071067811865475+0j)
1 NaN (0.7071067811865475-0j)
2 NaN (0.7071067811865475+0j)
1 3 -1 NaN (-0.7071067811865475+0j)
1 NaN (-0.7071067811865475-0j)
4 -2 NaN (-0.7071067811865475+0j)
-1 NaN (-0.7071067811865475-0j)
1 NaN (0.7071067811865475-0j)
2 NaN (0.7071067811865475+0j)
2 0 4 0 NaN (1+0j)

Running the method assignSymMuTerms will additionally attempt to subselect only the allowed symmetries, based on \(\Gamma_{\mathrm{electron}}\otimes\Gamma_{\mathrm{dipole}} \supseteq \Gamma_{\mathrm{symmetric}}\). This uses self.continuum['dict'] to define the allowed terms, and subselects accordingly. Polarization mu terms are also set according to symmetry.

[17]:

# Default assignment will set the allowed terms for dipole-allowed symmetries only, i.e. only uses (dipole X continuum) product terms, with mu values defined by symmetry.
symBasis.assignSymMuTerms()
symBasis.coeffs['symAllowed']['PD'].fillna('')

Assigned dipole-allowed terms for dim = 'Cont' to self.coeffs['symAllowed']

[17]:
Cont A1g E1g E2g
Eke Targ Total Type h it l m mu
0 U U U 0 0 0 0 -1 (1+0j)
1 (1+0j)
2 -2 -1 (0.7071067811865475+0j)
1 (0.7071067811865475+0j)
-1 0 (0.7071067811865475+0j)
1 0 (0.7071067811865475-0j)
2 -1 (0.7071067811865475+0j)
1 (0.7071067811865475+0j)
1 2 -2 -1 (-0.7071067811865475+0j)
1 (-0.7071067811865475+0j)
-1 0 (-0.7071067811865475-0j)
1 0 (0.7071067811865475-0j)
2 -1 (0.7071067811865475+0j)
1 (0.7071067811865475+0j)
1 0 2 0 -1 (1+0j)
1 (1+0j)
4 -2 -1 (0.7071067811865475+0j)
1 (0.7071067811865475+0j)
-1 0 (0.7071067811865475+0j)
1 0 (0.7071067811865475-0j)
2 -1 (0.7071067811865475+0j)
1 (0.7071067811865475+0j)
1 4 -2 -1 (-0.7071067811865475+0j)
1 (-0.7071067811865475+0j)
-1 0 (-0.7071067811865475-0j)
1 0 (0.7071067811865475-0j)
2 -1 (0.7071067811865475+0j)
1 (0.7071067811865475+0j)
2 0 4 0 -1 (1+0j)
1 (1+0j)

In some cases, the above may be sufficient. However, for a specific case, we can additionally assign all the symmetries matching ePolyScat definitions, and use the full direct products.

In this case:

  1. Cont is the continuum (free electron) symmetry, \(\Gamma_{\mathrm{e}}\).

  2. Targ is the target state symmetry, \(\Gamma_{\mathrm{ion}}\).

  3. Total is the overall symmetry of the scattering state, \(\Gamma_{\mathrm{scat}}=\Gamma_{\mathrm{ion}}\otimes\Gamma_{\mathrm{electron}}\).

[18]:

# To assign specific terms, use self.assignMissingSym
# Note this can take a single value, or a list which must match the size of the Sym multiindex defined in the Xarray dataset.
symBasis.assignMissingSym('Targ', sIon)

*** Updated self.coeffs['matE'] with new coords.

[19]:

# To define terms from produts, use self.assignMissingSymProd
symBasis.assignMissingSymProd()

# To update the outputs, rerun assignSymMuTerms
symBasis.assignSymMuTerms()
symBasis.coeffs['symAllowed']['PD'].fillna('')

Assigned 'Total' from A1g x E1u = ['E1u']
Assigned 'Total' from A2u x E1u = ['E1g']
Found multiple products for E1g x E1u = ['A1u', 'A2u', 'E2u']; assigning 'Total' as A1u.
Found multiple products for E1u x E1u = ['A1g', 'A2g', 'E2g']; assigning 'Total' as A1g.
Found multiple products for E2g x E1u = ['E1u', 'E3u']; assigning 'Total' as E1u.
Found multiple products for E2u x E1u = ['E1g', 'E3g']; assigning 'Total' as E1g.
Found multiple products for E3g x E1u = ['E2u', 'E4u']; assigning 'Total' as E2u.
Found multiple products for E3u x E1u = ['E2g', 'E4g']; assigning 'Total' as E2g.
Found multiple products for E4g x E1u = ['B1u', 'B2u', 'E3u']; assigning 'Total' as B1u.
*** Updated self.coeffs['matE'] with new coords.
Assigned dipole-allowed terms for dim = 'Cont' to self.coeffs['symAllowed']

[19]:
Cont A1g E1g E2g
Eke Targ Total Type h it l m mu
0 E1u A1u U 0 0 2 -1 0 (0.7071067811865475+0j)
1 0 (0.7071067811865475-0j)
1 2 -1 0 (-0.7071067811865475-0j)
1 0 (0.7071067811865475-0j)
1 0 4 -1 0 (0.7071067811865475+0j)
1 0 (0.7071067811865475-0j)
1 4 -1 0 (-0.7071067811865475-0j)
1 0 (0.7071067811865475-0j)
E1u U 0 0 0 0 -1 (1+0j)
1 (1+0j)
2 -2 -1 (0.7071067811865475+0j)
1 (0.7071067811865475+0j)
2 -1 (0.7071067811865475+0j)
1 (0.7071067811865475+0j)
1 2 -2 -1 (-0.7071067811865475+0j)
1 (-0.7071067811865475+0j)
2 -1 (0.7071067811865475+0j)
1 (0.7071067811865475+0j)
1 0 2 0 -1 (1+0j)
1 (1+0j)
4 -2 -1 (0.7071067811865475+0j)
1 (0.7071067811865475+0j)
2 -1 (0.7071067811865475+0j)
1 (0.7071067811865475+0j)
1 4 -2 -1 (-0.7071067811865475+0j)
1 (-0.7071067811865475+0j)
2 -1 (0.7071067811865475+0j)
1 (0.7071067811865475+0j)
2 0 4 0 -1 (1+0j)
1 (1+0j)

This contains the same terms as previously, but split by Total symmetry groupings - in this particular case \(\mu=0\) and \(\mu=\pm 1\) correspond to different total symmetries.

Note that, as compared to ePolyScat outputs, there may be some differences in these mappings. In particular there may be additional \(\mu\) or \(m\) phase terms, particularly for degenerate states, which may need to be addressed. These are currently split by muX in the default mappings, but usually these terms can be collapsed into a degeneracy and phase factor, or otherwise interpreted as equivalent degenerate components.

Comparison with ePS matrix elements

For comparison of a given symmetry with sample ePolyScat calculations, results can be computed locally or pulled from the web. Some sample/test datasets can be found as part of the ePSproc repo, which includes \(N_2\). For more datasets try ePSdata, and data can be pulled using the python interface.

[20]:

# Pull N2 data from ePSproc Github repo
import wget
from pathlib import Path

# URLs for test ePSproc datasets - n2
# For more datasets use ePSdata, see https://epsproc.readthedocs.io/en/dev/demos/ePSdata_download_demo_300720.html
urls = {'n2PU':"https://github.com/phockett/ePSproc/blob/master/data/photoionization/n2_multiorb/n2_1pu_0.1-50.1eV_A2.inp.out",
        'n2SU':"https://github.com/phockett/ePSproc/blob/master/data/photoionization/n2_multiorb/n2_3sg_0.1-50.1eV_A2.inp.out"}

# Set data dir
dataPath = Path(Path.cwd(), 'n2Data')

# Create and pull files if dir not present (NOTE doesn't check per file here)
if not dataPath.is_dir():
    dataPath.mkdir()

    # Pull files with wget
    for k,v in urls.items():
        wget.download(v+'?raw=true',out=dataPath.as_posix())  # For Github add '?raw=true' to URL


# List files
list(dataPath.glob('*.out'))

[20]:

[PosixPath('/home/jovyan/code-share/jupyter-shared/PEMtk_dev_2022/symmetry/n2Data/n2_1pu_0.1-50.1eV_A2.inp.out'),
 PosixPath('/home/jovyan/code-share/jupyter-shared/PEMtk_dev_2022/symmetry/n2Data/n2_3sg_0.1-50.1eV_A2.inp.out')]

[21]:

# # Import data
# # For more details on ePSproc usage see https://epsproc.readthedocs.io/en/dev/demos/ePSproc_class_demo_161020.html

# # from epsproc.classes.multiJob import ePSmultiJob
# from epsproc.classes.base import ePSbase

# # Instantiate class object.
# # Minimally this needs just the dataPath, if verbose = 1 is set then some useful output will also be printed.
# data = ePSbase(dataPath, verbose = 1)

# # ScanFiles() - this will look for data files on the path provided, and read from them.
# data.scanFiles()

[22]:

# Import data - PEMtk object
# For more details on ePSproc usage see https://epsproc.readthedocs.io/en/dev/demos/ePSproc_class_demo_161020.html

# from epsproc.classes.multiJob import ePSmultiJob
# from epsproc.classes.base import ePSbase
from pemtk.fit.fitClass import pemtkFit

# Instantiate class object.
# Minimally this needs just the dataPath, if verbose = 1 is set then some useful output will also be printed.
data = pemtkFit(fileBase=dataPath, verbose = 1)

# ScanFiles() - this will look for data files on the path provided, and read from them.
data.scanFiles()


*** Job subset details
Key: subset
No 'job' info set for self.data[subset].

*** Job orb6 details
Key: orb6
Dir /home/jovyan/code-share/jupyter-shared/PEMtk_dev_2022/symmetry/n2Data, 1 file(s).
{   'batch': 'ePS n2, batch n2_1pu_0.1-50.1eV, orbital A2',
    'event': ' N2 A-state (1piu-1)',
    'orbE': -17.09691397835426,
    'orbLabel': '1piu-1'}

*** Job orb5 details
Key: orb5
Dir /home/jovyan/code-share/jupyter-shared/PEMtk_dev_2022/symmetry/n2Data, 1 file(s).
{   'batch': 'ePS n2, batch n2_3sg_0.1-50.1eV, orbital A2',
    'event': ' N2 X-state (3sg-1)',
    'orbE': -17.34181645456815,
    'orbLabel': '3sg-1'}

[23]:

# Tabulate some demo matrix elements
singleEdata = ep.matEleSelector(data.data['orb6']['matE'], inds={'Eke':slice(0.001,1,1),'Type':'L'}, thres=1e-3)  #, drop=False, sq=False)   # CURRENTLY ALWAYS DROPPING TYPE HERE
# singleEdata = ep.matEleSelector(data.data['orb6']['matE'], inds={'Eke':slice(0.001,1,1)}, thres=1e-3)
# singleEdata = data.Esubset(key='orb6',dataType='matE',Erange=[0,1,1],Etype='Eke')

matEPD, _ = ep.multiDimXrToPD(singleEdata, colDims = 'Cont', thres=1e-4)  #, squeeze=False, dropna=False)
# matEPD.fillna('')

[24]:

# Compare results
# display(symBasis.coeffs['symAllowed']['PD'].fillna(''))
# display(matEPD.fillna(''))

[25]:

# Compare results
# For side-by-side tables, code adapted from https://stackoverflow.com/a/50899244

from IPython.display import display_html
# df1 = symBasis.coeffs['symAllowed']['PD'].droplevel('h').droplevel('Type').sort_index().fillna('')   # .sort_index(level='it', sort_remaining=False)
df1 = symBasis.coeffs['symAllowed']['PD'].droplevel('Type').sort_index().fillna('')   # .sort_index(level='it', sort_remaining=False)
df2 = matEPD.fillna('').sort_index().sort_index(level='Total', ascending=[False]).sort_index(axis=1, ascending=False)  # Swap sym label ordering to match symmetry case. Note ascenting=[False] for multindex case single level only.

df1_styler = df1.style.set_table_attributes("style='display:inline'").set_caption('<b>Symmetry basis</b>')
df2_styler = df2.style.set_table_attributes("style='display:inline'").set_caption('<b>ePS basis</b>')

display_html(df1_styler._repr_html_()+df2_styler._repr_html_(), raw=True)
Symmetry basis
              Cont A1g E1g E2g
Eke Targ Total h it l m mu      
0 E1u A1u 0 0 2 -1 0 0.707107+0.000000j
1 0 0.707107-0.000000j
1 2 -1 0 -0.707107-0.000000j
1 0 0.707107-0.000000j
1 0 4 -1 0 0.707107+0.000000j
1 0 0.707107-0.000000j
1 4 -1 0 -0.707107-0.000000j
1 0 0.707107-0.000000j
E1u 0 0 0 0 -1 1.000000+0.000000j
1 1.000000+0.000000j
2 -2 -1 0.707107+0.000000j
1 0.707107+0.000000j
2 -1 0.707107+0.000000j
1 0.707107+0.000000j
1 2 -2 -1 -0.707107+0.000000j
1 -0.707107+0.000000j
2 -1 0.707107+0.000000j
1 0.707107+0.000000j
1 0 2 0 -1 1.000000+0.000000j
1 1.000000+0.000000j
4 -2 -1 0.707107+0.000000j
1 0.707107+0.000000j
2 -1 0.707107+0.000000j
1 0.707107+0.000000j
1 4 -2 -1 -0.707107+0.000000j
1 -0.707107+0.000000j
2 -1 0.707107+0.000000j
1 0.707107+0.000000j
2 0 4 0 -1 1.000000+0.000000j
1 1.000000+0.000000j
ePS basis
            Cont SG PG DG
Eke Targ Total it l m mu      
0.100000 PU SU 1 2 -1 0 -2.896797+3.321073j
1 0 -2.896797+3.321073j
4 -1 0 -0.004197+0.035179j
1 0 -0.004197+0.035179j
2 2 -1 0 -3.321073-2.896797j
1 0 3.321073+2.896797j
4 -1 0 -0.035179-0.004197j
1 0 0.035179+0.004197j
PU 1 0 0 -1 1.455352-0.168378j
1 1.455352-0.168378j
2 -2 1 2.629738-0.365731j
0 -1 0.833527-0.752452j
1 0.833527-0.752452j
2 -1 2.629738-0.365731j
4 -2 1 0.011578-0.013731j
0 -1 0.005033-0.008819j
1 0.005033-0.008819j
2 -1 0.011578-0.013731j
2 0 0 -1 0.168378+1.455352j
1 -0.168378-1.455352j
2 -2 1 0.365731+2.629738j
0 -1 0.752452+0.833527j
1 -0.752452-0.833527j
2 -1 -0.365731-2.629738j
4 -2 1 0.013731+0.011578j
0 -1 0.008819+0.005033j
1 -0.008819-0.005033j
2 -1 -0.013731-0.011578j 
[26]:

# Compare results
# For side-by-side tables, code adapted from https://stackoverflow.com/a/50899244

from IPython.display import display_html

# Full case
# # df1 = symBasis.coeffs['symAllowed']['PD'].droplevel('h').droplevel('Type').sort_index().fillna('')   # .sort_index(level='it', sort_remaining=False)
# df1 = symBasis.coeffs['symAllowed']['PD'].droplevel('Type').sort_index().fillna('')   # .sort_index(level='it', sort_remaining=False)
# df2 = matEPD.fillna('').sort_index().sort_index(level='Total', ascending=[False]).sort_index(axis=1, ascending=False)  # Swap sym label ordering to match symmetry case. Note ascenting=[False] for multindex case single level only.

# df1_styler = df1.style.set_table_attributes("style='display:inline'").set_caption('<b>Symmetry basis</b>')
# df2_styler = df2.style.set_table_attributes("style='display:inline'").set_caption('<b>ePS basis</b>')

# display_html(df1_styler._repr_html_()+df2_styler._repr_html_(), raw=True)


# Check specific symmetry sets...

# df1['A1g'].dropna()
syms = ['A1g','SG']
# df1 = symBasis.coeffs['symAllowed']['PD'].droplevel('Type').sort_index()[syms[0]].dropna().to_frame()
df1 = symBasis.coeffs['symAllowed']['PD'].droplevel('h').droplevel('Type').sort_index()[syms[0]].dropna().to_frame()
df2 = matEPD.sort_index().sort_index(level='Total', ascending=[False]).sort_index(axis=1, ascending=False)[syms[1]].dropna().to_frame()

df1_styler = df1.style.set_table_attributes("style='display:inline'").set_caption('<b>Symmetry basis</b>')
df2_styler = df2.style.set_table_attributes("style='display:inline'").set_caption('<b>ePS basis</b>')

display_html(df1_styler._repr_html_()+df2_styler._repr_html_(), raw=True)
Symmetry basis
              A1g
Eke Targ Total it l m mu  
0 E1u E1u 0 0 0 -1 1.000000+0.000000j
1 1.000000+0.000000j
2 0 -1 1.000000+0.000000j
1 1.000000+0.000000j
4 0 -1 1.000000+0.000000j
1 1.000000+0.000000j
ePS basis
              SG
Eke Targ Total it l m mu  
0.100000 PU PU 1 0 0 -1 1.455352-0.168378j
1 1.455352-0.168378j
2 0 -1 0.833527-0.752452j
1 0.833527-0.752452j
4 0 -1 0.005033-0.008819j
1 0.005033-0.008819j
2 0 0 -1 0.168378+1.455352j
1 -0.168378-1.455352j
2 0 -1 0.752452+0.833527j
1 -0.752452-0.833527j
4 0 -1 0.008819+0.005033j
1 -0.008819-0.005033j

Here the ePS results have two degenerate continuua, it=1,2, with a phase rotation applied between them. The symmetry-defined harmonics here equate to the it=1 (in-phase) component only.

[27]:

# Check specific sets...

# df1['A1g'].dropna()
syms = ['E1g','PG']
# df1 = symBasis.coeffs['symAllowed']['PD'].droplevel('Type').sort_index()[syms[0]].dropna().to_frame()
df1 = symBasis.coeffs['symAllowed']['PD'].droplevel('h').droplevel('Type').sort_index()[syms[0]].dropna().to_frame()
df2 = matEPD.sort_index().sort_index(level='Total', ascending=[False]).sort_index(axis=1, ascending=False)[syms[1]].dropna().to_frame()

df1_styler = df1.style.set_table_attributes("style='display:inline'").set_caption('<b>Symmetry basis</b>')
df2_styler = df2.style.set_table_attributes("style='display:inline'").set_caption('<b>ePS basis</b>')

display_html(df1_styler._repr_html_()+df2_styler._repr_html_(), raw=True)
Symmetry basis
              E1g
Eke Targ Total it l m mu  
0 E1u A1u 0 2 -1 0 0.707107+0.000000j
1 0 0.707107-0.000000j
4 -1 0 0.707107+0.000000j
1 0 0.707107-0.000000j
1 2 -1 0 -0.707107-0.000000j
1 0 0.707107-0.000000j
4 -1 0 -0.707107-0.000000j
1 0 0.707107-0.000000j
ePS basis
              PG
Eke Targ Total it l m mu  
0.100000 PU SU 1 2 -1 0 -2.896797+3.321073j
1 0 -2.896797+3.321073j
4 -1 0 -0.004197+0.035179j
1 0 -0.004197+0.035179j
2 2 -1 0 -3.321073-2.896797j
1 0 3.321073+2.896797j
4 -1 0 -0.035179-0.004197j
1 0 0.035179+0.004197j

In this case the sets agree - both have two degenerate components, with in-phase and anti-phase (opposite sign) \(m\) pairs.

[28]:

# Check specific sets...

# df1['A1g'].dropna()
syms = ['E2g','DG']
df1 = symBasis.coeffs['symAllowed']['PD'].droplevel('h').droplevel('Type').sort_index()[syms[0]].dropna().to_frame()
# df1 = symBasis.coeffs['symAllowed']['PD'].droplevel('Type').sort_index()[syms[0]].dropna().to_frame()
df2 = matEPD.sort_index().sort_index(level='Total', ascending=[False]).sort_index(axis=1, ascending=False)[syms[1]].dropna().to_frame()

df1_styler = df1.style.set_table_attributes("style='display:inline'").set_caption('<b>Symmetry basis</b>')
df2_styler = df2.style.set_table_attributes("style='display:inline'").set_caption('<b>ePS basis</b>')

display_html(df1_styler._repr_html_()+df2_styler._repr_html_(), raw=True)
Symmetry basis
              E2g
Eke Targ Total it l m mu  
0 E1u E1u 0 2 -2 -1 0.707107+0.000000j
1 0.707107+0.000000j
2 -1 0.707107+0.000000j
1 0.707107+0.000000j
4 -2 -1 0.707107+0.000000j
1 0.707107+0.000000j
2 -1 0.707107+0.000000j
1 0.707107+0.000000j
1 2 -2 -1 -0.707107+0.000000j
1 -0.707107+0.000000j
2 -1 0.707107+0.000000j
1 0.707107+0.000000j
4 -2 -1 -0.707107+0.000000j
1 -0.707107+0.000000j
2 -1 0.707107+0.000000j
1 0.707107+0.000000j
ePS basis
              DG
Eke Targ Total it l m mu  
0.100000 PU PU 1 2 -2 1 2.629738-0.365731j
2 -1 2.629738-0.365731j
4 -2 1 0.011578-0.013731j
2 -1 0.011578-0.013731j
2 2 -2 1 0.365731+2.629738j
2 -1 -0.365731-2.629738j
4 -2 1 0.013731+0.011578j
2 -1 -0.013731-0.011578j

Here the ePS results are restricted to terms with \(m\) and \(mu\) of opposite sign. This is a restriction imposed by angular momentum selection rules, so is not included in the symmetry-defined case. Note, however, that the two it components are correctly accounted for in this case (for \(m\neq0\)).

In conclusion, the current mappings should be suitable for simulation and reconstruction, but care should be taken to:

  1. Confirm symmetry and angular momentum relations for a given case.

  2. Apply additional transformations for comparison with computational results.

  3. Add degeneracy factors if required (otherwise will be subsumed into matrix element values).

TODO: address some of these points in this notebook.

Versions

[29]:

import scooby
scooby.Report(additional=['epsproc', 'pemtk', 'holoviews', 'hvplot', 'xarray', 'matplotlib', 'bokeh', 'pyshtools'])

[29]:
Fri Apr 21 12:38:30 2023 EDT
OS Linux CPU(s) 32 Machine x86_64 Architecture 64bit
RAM 50.1 GiB Environment Jupyter
Python 3.9.10 | packaged by conda-forge | (main, Feb 1 2022, 21:24:11) [GCC 9.4.0]
epsproc 1.3.2-dev pemtk 0.0.1 holoviews 1.15.4 hvplot 0.8.3
xarray 2022.3.0 matplotlib 3.5.1 bokeh 2.4.2 pyshtools v4.10
numpy 1.21.5 scipy 1.8.0 IPython 8.1.1 scooby 0.7.1 
[30]:

# Check current Git commit for local ePSproc version
import pemtk
from pathlib import Path
!git -C {Path(pemtk.__file__).parent} branch
!git -C {Path(pemtk.__file__).parent} log --format="%H" -n 1

* master
  mfFittingDev
cfd7b4078b515a208f77d9ee4efe6b074831d5a5

[31]:

# Check current remote commits
!git ls-remote --heads https://github.com/phockett/pemtk

cfd7b4078b515a208f77d9ee4efe6b074831d5a5        refs/heads/master
3f4686dffdbb310f15692f978ba36d6a3d15e8d3        refs/heads/mfFittingDev