Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

SpectrumWidget: Switch excited energies to atomic units #65

Draft
wants to merge 6 commits into
base: main
Choose a base branch
from
Draft

SpectrumWidget: Switch excited energies to atomic units #65

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
0ee24ca
Fix intensity
danielhollas Sep 2, 2022
dacd6f2
Fixes
danielhollas Sep 5, 2022
7445092
Adjust header and slider
danielhollas Sep 5, 2022
ecb8dec
Fix width unit
danielhollas Sep 5, 2022
5b56c78
Fix once more
danielhollas Sep 5, 2022
5700643
Width fix
danielhollas Sep 5, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
69 changes: 35 additions & 34 deletions aiidalab_ispg/spectrum.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,10 @@
XyData = DataFactory("array.xy")


# Atomic units to electronvolts
AUtoEV = 27.211386245


@unique
class EnergyUnit(Enum):
EV = "eV"
Expand Down Expand Up @@ -61,18 +65,16 @@ def update(self):


class Spectrum(object):
AUtoCm = 8.478354e-30
COEFF = (
constants.pi
* AUtoCm**2
* 8.478354e-30**2 # AUtoCm
* AUtoEV
* 1e4
/ (3 * constants.hbar * constants.epsilon_0 * constants.c)
/ (2 * constants.hbar * constants.epsilon_0 * constants.c)
)
# Transition Dipole to Osc. Strength in atomic units
COEFF_NEW = COEFF * 3 / 2

def __init__(self, transitions: dict, nsample: int):
# Excitation energies in eV
# Excitation energies in a.u.
self.excitation_energies = np.array(
[tr["energy"] for tr in transitions], dtype=float
)
Expand All @@ -83,52 +85,46 @@ def __init__(self, transitions: dict, nsample: int):
# Number of molecular geometries sampled from ground state distribution
self.nsample = nsample

def _get_energy_range_ev(self):
"""Get spectrum energy range in eV"""
def _get_energy_range_au(self):
"""Get spectrum energy range in a.u."""
# NOTE: We don't include zero to prevent
# division by zero when converting to wavelength
x_min = max(0.01, self.excitation_energies.min() - 2.0)
x_max = self.excitation_energies.max() + 2.0
x_min = max(0.01, self.excitation_energies.min() - 1.0 / AUtoEV)
x_max = self.excitation_energies.max() + 1.0 / AUtoEV
return x_min, x_max

@staticmethod
def get_energy_unit_factor(unit: EnergyUnit):
"""Returns a multiplication factor to go from eV to other energy units"""
"""Returns a multiplication factor to go from a.u. to other energy units"""

# https://physics.nist.gov/cgi-bin/cuu/Info/Constants/basis.html
# TODO: We should probably start from atomic units
if unit is EnergyUnit.EV:
return 1.0
return AUtoEV
# TODO: Construct these factors from scipy.constants or use pint
elif unit is EnergyUnit.NM:
return 1239.8
return 1239.8 * AUtoEV
elif unit is EnergyUnit.CM:
# https://physics.nist.gov/cgi-bin/cuu/Convert?exp=0&num=1&From=ev&To=minv&Action=Only+show+factor
return 8065.547937
return 8065.547937 * AUtoEV

def calc_lorentzian_spectrum(self, x, y, tau: float):
normalization_factor = tau / 2 / constants.pi / self.nsample
unit_factor = self.COEFF_NEW

for exc_energy, osc_strength in zip(
self.excitation_energies, self.osc_strengths
):
prefactor = normalization_factor * unit_factor * osc_strength
prefactor = normalization_factor * self.COEFF * osc_strength
y += prefactor / ((x - exc_energy) ** 2 + (tau**2) / 4)

def calc_gauss_spectrum(self, x, y, sigma: float):
normalization_factor = (
1 / np.sqrt(2 * constants.pi) / sigma / self.nsample
)
unit_factor = self.COEFF_NEW
normalization_factor = 1 / np.sqrt(2 * constants.pi) / sigma / self.nsample
for exc_energy, osc_strength in zip(
self.excitation_energies, self.osc_strengths
):
prefactor = normalization_factor * unit_factor * osc_strength
prefactor = normalization_factor * self.COEFF * osc_strength
y += prefactor * np.exp(-((x - exc_energy) ** 2) / 2 / sigma**2)

def get_spectrum(self, kernel: BroadeningKernel, width: float, x_unit: EnergyUnit):
x_min, x_max = self._get_energy_range_ev()
x_min, x_max = self._get_energy_range_au()

# TODO: How to determine this properly to cover a given interval?
n_sample = 500
Expand Down Expand Up @@ -184,7 +180,12 @@ def __init__(self, **kwargs):
)

self.width_slider = ipw.FloatSlider(
min=0.05, max=1, step=0.05, value=0.1, description="Width / eV"
min=0.01,
max=0.5,
step=0.01,
value=0.05,
description="Width / eV",
continuous_update=True,
)

self.kernel_selector = ipw.ToggleButtons(
Expand Down Expand Up @@ -286,11 +287,13 @@ def _prepare_payload(self):
delimiter = "\t"

fieldnames = [
f"Energy / {self.energy_unit_selector.value}",
f"Energy / {self.energy_unit_selector.value.value}",
f"Intensity / {self.intensity_unit}",
f"{self.kernel_selector.value.value} broadening, width = {self.width_slider.value} eV",
]
with SpooledTemporaryFile(mode="w+", newline="", max_size=10000000) as csvfile:
csvfile.write(f"# {fieldnames[0]}{delimiter}{fieldnames[1]}\n")
header = delimiter.join(fieldnames)
csvfile.write(f"# {header}\n")
writer = csv.writer(csvfile, delimiter=delimiter)
writer.writerows(zip(x, y))
csvfile.seek(0)
Expand All @@ -312,19 +315,17 @@ def _validate_transitions(self):

def _handle_width_update(self, change):
"""Redraw spectra when user changes broadening width via slider"""
width = change["new"]
self._plot_spectrum(
width=width,
width=change["new"] / AUtoEV,
kernel=self.kernel_selector.value,
energy_unit=self.energy_unit_selector.value,
)

def _handle_kernel_update(self, change):
"""Redraw spectra when user changes kernel for broadening"""
kernel = change["new"]
self._plot_spectrum(
width=self.width_slider.value,
kernel=kernel,
width=self.width_slider.value / AUtoEV,
kernel=change["new"],
energy_unit=self.energy_unit_selector.value,
)

Expand All @@ -336,7 +337,7 @@ def _handle_energy_unit_update(self, change):
self.figure.get_figure().xaxis.axis_label = xlabel

self._plot_spectrum(
width=self.width_slider.value,
width=self.width_slider.value / AUtoEV,
kernel=self.kernel_selector.value,
energy_unit=energy_unit,
)
Expand Down Expand Up @@ -443,7 +444,7 @@ def reset(self):
@traitlets.observe("transitions")
def _observe_transitions(self, change):
self._plot_spectrum(
width=self.width_slider.value,
width=self.width_slider.value / AUtoEV,
kernel=self.kernel_selector.value,
energy_unit=self.energy_unit_selector.value,
)
Expand Down
11 changes: 3 additions & 8 deletions aiidalab_ispg/steps.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@
except ImportError:
print("ERROR: Could not find aiidalab_atmospec_workchain module!")

from aiidalab_ispg.spectrum import SpectrumWidget
from aiidalab_ispg.spectrum import EnergyUnit, Spectrum, SpectrumWidget

StructureData = DataFactory("structure")
TrajectoryData = DataFactory("array.trajectory")
Expand Down Expand Up @@ -591,17 +591,12 @@ def reset(self):
self.process = None
self.spectrum.reset()

# TODO: Move this to the workflow
def _orca_output_to_transitions(self, output_dict, geom_index):
# TODO: Use atomic units both for energies and osc. strengths
CM2EV = 1 / 8065.547937
# TODO: Add error handling
AUtoCM = Spectrum.get_energy_unit_factor(EnergyUnit.CM)
en = output_dict["etenergies"]
osc = output_dict["etoscs"]
assert len(en) == len(osc)
# TODO: Use atomic units both for energies and osc. strengths
return [
{"energy": tr[0] * CM2EV, "osc_strength": tr[1], "geom_index": geom_index}
{"energy": tr[0] / AUtoCM, "osc_strength": tr[1], "geom_index": geom_index}
for tr in zip(en, osc)
]

Expand Down