diff --git a/aiidalab_ispg/spectrum.py b/aiidalab_ispg/spectrum.py index 8340b24..45658fe 100644 --- a/aiidalab_ispg/spectrum.py +++ b/aiidalab_ispg/spectrum.py @@ -22,6 +22,10 @@ XyData = DataFactory("array.xy") +# Atomic units to electronvolts +AUtoEV = 27.211386245 + + @unique class EnergyUnit(Enum): EV = "eV" @@ -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 ) @@ -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 @@ -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( @@ -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) @@ -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, ) @@ -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, ) @@ -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, ) diff --git a/aiidalab_ispg/steps.py b/aiidalab_ispg/steps.py index c9e54b9..bd95f54 100644 --- a/aiidalab_ispg/steps.py +++ b/aiidalab_ispg/steps.py @@ -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") @@ -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) ]