OCS ePS results processing + AF-BLMs (Jan 2022 update) + test ADMs from VM 14/01/22¶
Quick attempt at AF-BLMs for OCS, based on guessed/reasonable ADMs - action starts below after setup stage). Status: basically working, may be some remaining bugs with ADMs, cos^2 and renormalisation, but trends with alignment should be correct.
For methods: https://epsproc.readthedocs.io/en/dev/demos/ePSproc_class_demo_161020.html
Additional notes below
Setup¶
In [1]:
!hostname
In [2]:
!conda env list
In [3]:
import sys
import os
from pathlib import Path
import numpy as np
# import epsproc as ep
import xarray as xr
import matplotlib.pyplot as plt
from datetime import datetime as dt
timeString = dt.now()
import epsproc as ep
# Plotters
import hvplot.xarray
from epsproc.plot import hvPlotters
# Multijob class dev code
from epsproc.classes.multiJob import ePSmultiJob
plotWidth = 700
hvPlotters.setPlotters(width = plotWidth, snsStyle='whitegrid')
hvPlotters.opts.defaults(hvPlotters.opts.Curve(height=plotWidth)) # Fix plot height! Currently not set by default.
In [4]:
# For class, above settings don't take, not sure why, something to do with namespaces/calling sequence?
# Overriding snsStyle does work however... although NOT CONSISTENTLY????
# AH, looks like ordering matters - set_style LAST (.set seems to override)
import seaborn as sns
sns.set(rc={'figure.figsize':(10,6)}) # Set figure size in inches
sns.set_context("paper")
sns.set_style("whitegrid") # Set plot style
sns.set_palette("Paired") # Set colour mapping
# Try direct fig type setting for PDF output figs
from IPython.display import set_matplotlib_formats
# set_matplotlib_formats('png', 'pdf')
set_matplotlib_formats('svg', 'pdf')
In [5]:
# xr.set_options(display_style='html')
Load data¶
In [6]:
import warnings
# warnings.filterwarnings('once') # Skip repeated numpy deprecation warnings in current build (xr15 env)
warnings.filterwarnings('ignore') # Skip repeated numpy deprecation warnings in current build (xr15 env)
In [7]:
# # Scan for subdirs, based on existing routine in getFiles()
fileBase = Path('/home/paul/ePS/OCS/OCS_survey') # Data dir on Stimpy
In [8]:
# TODO: fix orb label here, currently relies on (different) fixed format
data = ePSmultiJob(fileBase, verbose = 0)
data.scanFiles()
data.jobsSummary()
In [9]:
# Fix labels for plots - currently have some mis-named files!
# Set state labels
v = ['C','B','A','X']
sDict = {n:v[k] for k,n in enumerate(range(9,13))}
for key in data.data.keys():
data.data[key]['jobNotes']['orbKey'] = key # Existing orb key, from filename
data.data[key]['jobNotes']['orbGroup'] = int(key.strip('orb')) + 1 # Corrected orb group
# data.data[key]['jobNotes']['orbLabel'] = f"HOMO - {12 - int(key.strip('orb')) - 1}"
data.data[key]['jobNotes']['orbLabel'] = sDict[data.data[key]['jobNotes']['orbGroup']]
# print(data.data[key]['jobNotes']['orbLabel'])
System properties¶
Note orbital numbering in table below:
Orbis Gamess file output orbital numbering, but energy but not grouped by degeneracy.OrbGrpis grouped numbering by degeneracy, also used by ePolyScat, and will be used for labels etc. below.
BUT - file names are 0-indexed (oops), so give OrbGrp-1 here. Sorry.
In [10]:
data.molSummary()
In [11]:
# Fix orb names! NOW FIXED ABOVE
AF-$\beta_{LM}$ (UPDATES 13/12/21)¶
14/01/22 PH
With additional realistic ADMs from VM.
13/12/21 PH
Revisiting AF cases with (hopefully) better match to experimental alignment.
- Add in cos^2(theta) metric.
- Ballpark ADMs assuming Gaussian case.
- Calculate & plot AF-BLMs for these cases.
Alignment metrics¶
Here we'll start with the molecular alignment defined in terms of axis distribution moments (ADMs), given by parameters $A^K_{Q,S}(t)$:
$ P(\Omega,t) = \sum_{K,Q,S} A^K_{Q,S}(t)D^K_{Q,S}(\Omega)$
Where $P(\Omega,t)$ is the full axis distribution probability, expanded for Euler angles $\Omega$.
This reduces to the 2D case if $S=0$:
$ P(\theta,\phi,t) = \sum_{K,Q} A^K_{Q,0}(t)D^K_{Q,0}(\Omega) = \sum_{K,Q} A^K_{Q}(t)Y_{K,Q}(\Omega)$
In the current code (v1.3.0, 16/08/21), the ADMs can be set directly, and expanded trivially for the 2D case (3D case to do). (Note that the 3D case is supported for AF-BLM calculations, but not yet for direct expansion & visualisation.)
Expectation values of $\langle cos^n(\theta, t)\rangle$ can be calculated directly from the $P(\theta,t)$ distributions:
$\langle cos^n(\theta, t)\rangle = \int_{\theta} cos^n(\theta)P(\theta,t) d\theta$
See, e.g.:
[1]
P. Hockett, “General phenomenology of ionization from aligned molecular ensembles,” New Journal of Physics, vol. 17, no. 2, p. 023069, Feb. 2015, doi: 10.1088/1367-2630/17/2/023069.
[2] Reid, K.L. (2018) ‘Accessing the molecular frame through strong-field alignment of distributions of gas phase molecules’, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 376(2115), p. 20170158. doi:10.1098/rsta.2017.0158.
Set test ADMs¶
Roughly eyeballed for a Gaussian distribution with $0<\sigma^{2}<0.2$ as per Figure 1 in Reid 2018 [2].
NOTE: not totally sure on the normalization here. Assuming othonorm sph, may be an additional sqrt(4pi) required for Y00? https://shtools.github.io/SHTOOLS/complex-spherical-harmonics.html#orthonormalized
UPDATE: setting for 4$pi$- normalised & converting to orthonorm (for numerical functions) looks good for comparison with ref. [2]. Note also additional factor of sqrt(2) in metrics at the moment - due to 0:$pi$ default domain and/or renorm(?)... should also confirm there is no additional ^2 in distribution generation.
In [12]:
# Load test ADMs from .mat files.
# Code adapted from N2 case, https://epsproc.readthedocs.io/en/latest/methods/geometric_method_dev_pt3_AFBLM_090620_010920_dev_bk100920.html#Test-compared-to-experimental-N2-AF-results...
from scipy.io import loadmat
ADMdataDir = Path('~/ePS/OCS/OCS_alignment_ADMs_VM_140122')
fList = ep.getFiles(fileBase = ADMdataDir.expanduser(), fType='.mat')
ADMs = []
ADMLabels = []
for f in fList:
item = Path(f).name.rstrip('ocs.mat') # Get & set variable name
if item == 'time':
item = 'timeocs'
t = loadmat(f)[item][0,:]
else:
K = int(item.strip('D'))
item = item + 'ave'
# ADMs.append(np.asarray([K,0,0,loadmat(f)[item][:,0]]))
# ADMs.append([K,0,0,loadmat(f)[item][:,0]])
ADMs.append(loadmat(f)[item][:,0])
ADMLabels.append([K,0,0])
# Add K = 0 (population) term?
# NOTE - may need additional renorm?
addPop = True
if addPop:
ADMs.append(np.ones(t.size) * np.sqrt(4*np.pi))
ADMLabels.append([0,0,0])
ADMLabels = np.array(ADMLabels)
In [13]:
# Set ADMs for increasing alignment...
tPoints = 10
# inputADMs = np.asarray([[0,0,0, *np.ones(10)], [2,0,0, *np.linspace(0,1,tPoints)], [4,0,0, *np.linspace(0,0.5,tPoints)]])
# inputADMs = [[0,0,0, *np.ones(10) * np.sqrt(4*np.pi)],
# [2,0,0, *np.linspace(1,2,tPoints)],
# [4,0,0, *np.linspace(0.3,2.8,tPoints)],
# [6,0,0, *np.linspace(0,3,tPoints)],
# [8,0,0, *np.linspace(0,3,tPoints)]]
# Ballpark ADMs from Fig. 1 in [2]
# inputADMs = [[0,0,0, *np.ones(10) * np.sqrt(4*np.pi)],
# [2,0,0, *np.linspace(0.3,2.2,tPoints)],
# [4,0,0, *np.linspace(0,1.5,tPoints)],
# [6,0,0, *np.linspace(0,1.5,tPoints)],
# [8,0,0, *np.linspace(0,1.5,tPoints)]]
warnings.filterwarnings('ignore') # Skip repeated numpy deprecation warnings in current build (xr15 env)
# For manual sets
# ADMs = ep.setADMs(ADMs = inputADMs) # TODO WRAP TO CLASS (IF NOT ALREADY!)
# For ADMs from VM
ADMs = ep.setADMs(ADMs = ADMs, KQSLabels = ADMLabels, t=t) # TODO WRAP TO CLASS (IF NOT ALREADY!)
attrCpy = ADMs.attrs
ADMs = ADMs * (2*ADMs.K+1)/(8*np.pi**2) # Additional renormalisation by (2*K+1)/8*pi^2
# Set to separate dataset for testing
ADMs.attrs = attrCpy # Propagate attrs
ADMs.attrs['jobLabel'] = 'Alignment data testing'
data.data['ADM'] = {'ADM': ADMs}
ADMs
# ep.lmPlot(ADMs, xDim='t');
Out[13]:
In [14]:
ADMs.unstack().squeeze().real.hvplot.line(x='t').overlay('K')
Out[14]:
In [15]:
# Selection & downsampling - adapted from https://epsproc.readthedocs.io/en/latest/methods/geometric_method_dev_pt3_AFBLM_090620_010920_dev_bk100920.html#Test-compared-to-experimental-N2-AF-results...
# See PEMtk for updated routines, https://pemtk.readthedocs.io/en/latest/fitting/PEMtk_fitting_basic_demo_030621-full.html#Subselect-data
# See Xarray docs for basics https://xarray.pydata.org/en/stable/user-guide/indexing.html#indexing-with-dimension-names
trange=[38, 44] # Set range in ps for calc
tStep=2 # Set tStep for downsampling
data.data['ADM'] = {'ADM': ADMs.sel(t=slice(trange[0],trange[1], tStep))} # Set and update
# tMask = (ADMs.t>trange[0]) * (ADMs.t<trange[1])
# ind = np.nonzero(tMask) #[0::tStep]
# At = ADMs['time'][:,ind].squeeze()
# ADMin = ADMs['ADM'][:,ind]
print(f"Selecting {data.data['ADM']['ADM'].t.size} points")
# plt.plot(At, np.real(ADMin.T))
# plt.legend(ADMs['ADMlist'])
# plt.xlabel('t/ps')
# plt.ylabel('ADM')
# plt.show()
# # Set in Xarray
# ADMX = ep.setADMs(ADMs = ADMs['ADM'][:,ind], t=At, KQSLabels = ADMs['ADMlist'], addS = True)
# # ADMX
data.data['ADM']['ADM'].unstack().squeeze().real.hvplot.line(x='t').overlay('K')
Out[15]:
Compute $P(\theta)$ distribtuions + metrics¶
In [16]:
# Axis distributions P(theta) vs. t (summed over phi)
# Pt = data.padPlot(keys = 'ADM', dataType='ADM', Etype = 't', pStyle='grid', reducePhi='sum', returnFlag = True) #, facetDims=['t', 'Eke'])
Pt = data.padPlot(keys = 'ADM', dataType='ADM', Etype = 't', pType='r', pStyle='grid', reducePhi='sum', returnFlag = True) #, facetDims=['t', 'Eke']) # Real plot to check values OK.
In [17]:
# The returned data array contains the plotted values
Pt
Out[17]:
In [18]:
# Check no negative values...
print(Pt.min() > 0)
In [19]:
def computeCNs(Pt, N = None, Nrange = [2,10], Nstep = 2):
"""
Function to compute <cos^N> expectation values from P(theta) distributions.
NOTE: currently assumes 1D case. Should generalise!
"""
if N is None:
N = np.arange(Nrange[0], Nrange[1], Nstep)
# Generate cos^N basis
# cn = xr.ones_like(Pt) * (np.cos(Pt.Theta)**N) # Use multiplication to keep t coord - this fails for array-like N with broadcast errors
cn = xr.ones_like(Pt) * np.cos(Pt.Theta).expand_dims({'N':N}).T.pipe(np.power,N) # OK with transpose
Jt = np.sin(Pt.Theta) # Jacobian
norm = (Pt * Jt).sum('Theta')
cnMetric = (Pt * cn * Jt).sum('Theta')/norm
# cnMetric = Pt.dot(cn, dims='Theta') # Multiply and sum over Theta - dot version (no Jacobian or renorm)
# cnMetric = (Pt.dot(cn * Jt, dims='Theta')/Pt.dot(Jt, dims='Theta')) # Dot version with renorm
return cnMetric
In [20]:
# Compute and plot for a range of N
PtNorm = Pt #/Pt.sum('Theta') # Norm integral to 1 - now set in function
cnMetric = computeCNs(PtNorm)
# cnMetric.plot() # Basic plot
sf = np.sqrt(2) # Seem to need additional sqrt(2) normalisation? Probably SPH definition and/or integration domain and/or distribution generation params.
# sf = 1
ep.lmPlot(cnMetric*sf, xDim='t'); # With lmPlot()
In [21]:
# Seem to need additional sqrt(2) normalisation? Probably SPH definition?
# (cnMetric*np.sqrt(2)).sel(N=2).plot()
# (cnMetric*np.sqrt(2)).plot.line(x='t')
# Testing exact values
# sf = np.sqrt(2)
# sf = 2
# (cnMetric*sf).plot.line(x='t')
# (ADMs).plot.line(x='t')
(cnMetric*sf).fillna(0).hvplot.line(x='t').overlay('N')
Out[21]:
AFBLMs for aligned case¶
New plotting code to update current package versions!
In [22]:
# Compute AFBLMs for aligned case
keys = data._keysCheck(None)
keys = list(set(keys) - set(['ADM']))
# data.AFBLM(keys = keys, AKQS = ADMs, selDims = {'Type':'L'}) # NOTE THIS USES independently set ADMs!
data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'], selDims = {'Type':'L'}, thres=1e-3) # thres=1e-4) BREAKS KERNEL #, thresDims=['Eke','t']) # TODO: currently thersDims must be in matE only for AFBLM routine.
# data.AFBLM(keys = keys, selDims = {'Type':'L'}) # Default case... sets/uses ADMs per key
In [23]:
# Quick stack to dict & XR dataset
dataType = 'AFBLM'
thres = 1e-2
selDims = None
Etype = 'Eke'
pDict = []
for key in keys:
# subset = ep.matEleSelector(data.data[key][dataType], thres=thres, inds = selDims, dims = [Etype, 't'], sq = True)
subset = ep.matEleSelector(data.data[key][dataType], thres=thres, inds = selDims, dims = 't', sq = True)
subset.name = 'BLM'
subset = ep.plotTypeSelector(subset, pType = 'r')
# pDict.append(subset.expand_dims({'Orb':[key]}))
# pDict.append(subset.expand_dims({'Orb':[data.data[key][dataType].attrs['jobLabel']]}))
pDict.append(subset.expand_dims({'Orb':[data.data[key]['jobNotes']['orbLabel']]}))
# hvDS = hv.Dataset(subset.unstack('BLM'))
# display(hvDS.to(hv.Curve, kdims=Etype).overlay(['l','m']))
# pDict[key] = hvDS.to(hv.Curve, kdims=Etype).overlay(['l','m'])
# hv.HoloMap(pDict, kdims = ['Orb']) # Individual plots OK, but won't stack?
xrAF = xr.concat(pDict, 'Orb')
In [24]:
xrAF.unstack().dropna(dim = 't', how = 'all').drop('m') #.fillna(0)
Out[24]:
In [25]:
# SUPER SLOW WITH LARGER DATASETS - must be a way to improve on this? Seems a bit better with subsets and if NaNs dropped properly, also if high-dim vars are overlaid rather than Holomapped
# hvPlotters.opts.defaults(hvPlotters.opts.Curve(height=700))
# lSel = [2,4,6]
# hvDS = hvPlotters.hv.Dataset(xrAF.unstack('BLM'))
# hvDS.select(l=lSel).to(hvPlotters.hv.Curve, kdims=Etype).overlay(['l','m'])
# hvDS.to(hvPlotters.hv.Curve, kdims=Etype).overlay(['l','m'])
# Try with subsets to speed up plotters
hvDS = hvPlotters.hv.Dataset(xrAF.unstack().dropna(dim = 't', how = 'all').drop('m').fillna(0))
# hvDS.to(hvPlotters.hv.Curve, kdims=Etype).overlay(['l','m'])
In [31]:
hvDS.to(hvPlotters.hv.Curve, kdims=Etype).overlay(['l'])
Out[31]:
In [26]:
# This is slow, but not so bad. NOW DROPPING TOOLTIPS and LEGEND??????
# With additional overlays
# hvDS.select(l=lSel, m=0).to(hvPlotters.hv.Curve, kdims=Etype).overlay(['t']).layout('Orb').cols(1)
hvDS.to(hvPlotters.hv.Curve, kdims=Etype).overlay(['t']).layout('Orb').cols(1)
Out[26]:
In [27]:
# Getting non-contiguous plots here... looks funky...? Might be automatic thresholding in AFBLM routine (==NaNs/zeros) - should look at more closely
hvDS.to(hvPlotters.hv.Curve, kdims='t').overlay([Etype]).layout('Orb').cols(1)
Out[27]:
In [30]:
# With additional overlays
# This seems to break legend?
# hvDS.select(l=[2,4,6]).to(hvPlotters.hv.Curve, kdims=Etype).overlay(['l','m','t']).layout('Orb').cols(1)
hvDS.to(hvPlotters.hv.Curve, kdims=Etype).overlay(['l','t']).layout('Orb').cols(1)
Out[30]:
In [ ]:
# Check also cross-secions
XS = np.abs(xrAF.XSrescaled) # Should be in MB
# XS = np.abs(xrAF.XSraw) # Unscaled/uncorrect XS
XS.name = 'XS/MB'
XS.hvplot.line(x='Eke').overlay('t').layout('Orb').cols(1)
Versions¶
Code¶
In [ ]:
# print(data.jobInfo['ePolyScat'][0])
# print('Run: ' + jobInfo['Starting'][0].split('at')[1])
In [ ]:
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter', 'holoviews'])
In [ ]:
# Check current Git commit for local ePSproc version
!git -C {Path(ep.__file__).parent} branch
!git -C {Path(ep.__file__).parent} log --format="%H" -n 1
In [ ]:
# Check current remote commits
!git ls-remote --heads git://github.com/phockett/ePSproc
Notes¶
14/01/22 PH
Rerunning with some additional test ADMs from Varun Makhija.
07/01/22 PH
- Fixed/tested cos^2 values. (Maybe - now matching Reid 2018, but may not be correct in all cases?)
- Updated plotting with HV for BLMs (in dev stage).
- Checked XS values too - may have some -ve cases/phase issues to sort? (Cf. general AF code testing notes elsewhere.)
Versions
- TIDY version for distribution: https://phockett.github.io/ePSdata/OCS-preliminary/OCS_orb11_proc_131221-JAKE.ipynb
- Previous HOMO results: https://phockett.github.io/ePSdata/OCS-preliminary/OCS_orb11_proc_090621-JAKE.html
- Previous all-orb restuls: https://phockett.github.io/ePSdata/OCS-preliminary/OCS_orbs8-11_AFPAD_preliminary_110621.html
13/12/21 PH
STATUS: basics in place, need to trim code and fix a few things still.
Revisiting AF cases with (hopefully) better match to experimental alignment.
- Add in cos^2(theta) metric.
- Ballpark ADMs assuming Gaussian case.
TODO:
- Fix orb numering issue! (See note below.)
FIXED:
- Note 'it' and 'labels' handling changed in recent versions of ePSproc, and degenerate dim hanling now correct!
09/06/21 PH
Quick look at initial ePS results for orb 11 (HOMO), 2.5eV step size.
TODO:
- Currently AF code fails for Xarray = 0.17, issue with phaseCons setting to coord, sigh.
- AF normalisation not correct here, due to doubly-degenerate state? Summing over 'it' not working correctly either, needs a debug (or missing some settings...?).
- lmPlot() slice and cmapping to fix, in cases of missing/dropped dims, and for clims if possible.
For methods: https://epsproc.readthedocs.io/en/dev/demos/ePSproc_class_demo_161020.html
In [ ]: