Add plugin: X-ray Photoelectron Spectroscopy (XPS) (#518)

This PR adds a plugin for calculating X-ray Photoelectron Spectroscopy (XPS) spectra using the XpsWorkChain of aiida-quantumespresso package. The Setting panel allows users to select the core-level (element and orbital) to be calculated. The availability of core-evel in the panel depends on the corresponding pseudopotentials.

The Result panel displays the final XPS spectrum
- Broadened spectrum using a Voight profile (combined Lorenzian and Gaussian).
- The chemical shift values (differences in total energy relative to the lowest value)
- Absolute binding energy
- The Offset Energies (δ)

To make the setting as simple as possible for user. The following default core_hole_treatment are used:
- `full` for molecule
- `xch_smear` for metal
- `xch_fixed` for insulator

The supercell_min_parameter is updated based on protocol.

---------


Co-authored-by: Peter Gillespie <55498719+PNOGillespie@users.noreply.github.com>

setup.cfg

@@ -55,6 +55,7 @@ aiidalab_qe.plugins.xas = pseudo_toc.yaml
 aiidalab_qe.properties =
     bands = aiidalab_qe.plugins.bands:bands
     pdos = aiidalab_qe.plugins.pdos:pdos
+    xps = aiidalab_qe.plugins.xps:xps
     electronic_structure = aiidalab_qe.plugins.electronic_structure:electronic_structure
     xas = aiidalab_qe.plugins.xas:xas
 

src/aiidalab_qe/app/structure/examples/SiO2.xyz

@@ -1,8 +1,8 @@
 6
-Lattice="4.1801 0.0 0.0 0.0 4.1801 0.0 0.0 0.0 2.6678" Properties=species:S:1:pos:R:3:tags:I:1 spacegroup="P 42/m n m" unit_cell=conventional pbc="T T T"
-Si       0.00000000       0.00000000       0.00000000        0
-Si       2.09005000       2.09005000       1.33390000        0
-O       1.28203667       1.28203667       0.00000000        1
-O       2.89806333       2.89806333       0.00000000        1
-O       3.37208667       0.80801333       1.33390000        1
-O       0.80801333       3.37208667       1.33390000        1
+Lattice="4.1801 0.0 0.0 0.0 4.1801 0.0 0.0 0.0 2.6678" Properties=species:S:1:pos:R:3 spacegroup="P 42/m n m" unit_cell=conventional pbc="T T T"
+Si       0.00000000       0.00000000       0.00000000
+Si       2.09005000       2.09005000       1.33390000
+O       1.28203667       1.28203667       0.00000000
+O       2.89806333       2.89806333       0.00000000
+O       3.37208667       0.80801333       1.33390000
+O       0.80801333       3.37208667       1.33390000

src/aiidalab_qe/plugins/xps/__init__.py

@@ -0,0 +1,18 @@
+from aiidalab_qe.common.panel import OutlinePanel
+
+from .result import Result
+from .setting import Setting
+from .workchain import workchain_and_builder
+
+
+class XpsOutline(OutlinePanel):
+    title = "X-ray photoelectron spectroscopy (XPS)"
+    help = """"""
+
+
+xps = {
+    "outline": XpsOutline,
+    "setting": Setting,
+    "result": Result,
+    "workchain": workchain_and_builder,
+}

src/aiidalab_qe/plugins/xps/result.py

@@ -0,0 +1,246 @@
+"""XPS results view widgets
+
+"""
+import ipywidgets as ipw
+
+from aiidalab_qe.common.panel import ResultPanel
+
+
+def export_xps_data(outputs):
+    """Export the data from the XPS workchain"""
+
+    chemical_shifts = {}
+    symmetry_analysis_data = outputs.symmetry_analysis_data.get_dict()
+    equivalent_sites_data = symmetry_analysis_data["equivalent_sites_data"]
+    if "chemical_shifts" in outputs:
+        for key, data in outputs.chemical_shifts.items():
+            ele = key[:-4]
+            chemical_shifts[ele] = data.get_dict()
+    binding_energies = {}
+    if "binding_energies" in outputs:
+        for key, data in outputs.binding_energies.items():
+            ele = key[:-3]
+            binding_energies[ele] = data.get_dict()
+
+    return (
+        chemical_shifts,
+        binding_energies,
+        equivalent_sites_data,
+    )
+
+
+def xps_spectra_broadening(
+    points, equivalent_sites_data, gamma=0.3, sigma=0.3, label=""
+):
+    """Broadening the XPS spectra with Voigt function and return the spectra data"""
+
+    import numpy as np
+    from scipy.special import voigt_profile  # pylint: disable=no-name-in-module
+
+    result_spectra = {}
+    fwhm_voight = gamma / 2 + np.sqrt(gamma**2 / 4 + sigma**2)
+    for element, point in points.items():
+        result_spectra[element] = {}
+        final_spectra_y_arrays = []
+        total_multiplicity = sum(
+            [equivalent_sites_data[site]["multiplicity"] for site in point]
+        )
+        max_core_level_shift = max(point.values())
+        min_core_level_shift = min(point.values())
+        # Energy range for the Broadening function
+        x_energy_range = np.linspace(
+            min_core_level_shift - fwhm_voight - 1.5,
+            max_core_level_shift + fwhm_voight + 1.5,
+            500,
+        )
+        for site in point:
+            # Weight for the spectra of every atom
+            intensity = equivalent_sites_data[site]["multiplicity"]
+            relative_core_level_position = point[site]
+            y = (
+                intensity
+                * voigt_profile(
+                    x_energy_range - relative_core_level_position, sigma, gamma
+                )
+                / total_multiplicity
+            )
+            result_spectra[element][site] = [x_energy_range, y]
+            final_spectra_y_arrays.append(y)
+        total = sum(final_spectra_y_arrays)
+        result_spectra[element]["total"] = [x_energy_range, total]
+    return result_spectra
+
+
+class Result(ResultPanel):
+    title = "XPS"
+    workchain_labels = ["xps"]
+
+    def __init__(self, node=None, **kwargs):
+        super().__init__(node=node, **kwargs)
+
+    def _update_view(self):
+        import plotly.graph_objects as go
+
+        spectrum_select_prompt = ipw.HTML(
+            """
+            <div style="line-height: 140%; padding-top: 10px; padding-right: 10px; padding-bottom: 0px;"><b>
+            Select spectrum to plot</b></div>"""
+        )
+
+        voigt_profile_help = ipw.HTML(
+            """<div style="line-height: 140%; padding-top: 10px; padding-bottom: 10px">
+                Set the <a href="https://en.wikipedia.org/wiki/Voigt_profile" target="_blank">Voigt profile</a> to broaden the XPS spectra:
+            </div>"""
+        )
+
+        spectra_type = ipw.ToggleButtons(
+            options=[
+                ("Chemical shift", "chemical_shift"),
+                ("Binding energy", "binding_energy"),
+            ],
+            value="chemical_shift",
+        )
+        gamma = ipw.FloatSlider(
+            value=0.3,
+            min=0.1,
+            max=1,
+            description="Lorentzian profile ($\gamma$)",
+            disabled=False,
+            style={"description_width": "initial"},
+        )
+        sigma = ipw.FloatSlider(
+            value=0.3,
+            min=0.1,
+            max=1,
+            description="Gaussian profile ($\sigma$)",
+            disabled=False,
+            style={"description_width": "initial"},
+        )
+        fill = ipw.Checkbox(
+            description="Fill",
+            value=True,
+            disabled=False,
+            style={"description_width": "initial"},
+        )
+        paras = ipw.HBox(
+            children=[
+                gamma,
+                sigma,
+            ]
+        )
+        # get data
+        (
+            chemical_shifts,
+            binding_energies,
+            equivalent_sites_data,
+        ) = export_xps_data(self.outputs.xps)
+        self.spectrum_select_options = [
+            key.split("_")[0] for key in chemical_shifts.keys()
+        ]
+        spectrum_select = ipw.Dropdown(
+            description="",
+            disabled=False,
+            value=self.spectrum_select_options[0],
+            options=self.spectrum_select_options,
+            layout=ipw.Layout(width="20%"),
+        )
+        # init figure
+        g = go.FigureWidget(
+            layout=go.Layout(
+                title=dict(text="XPS"),
+                barmode="overlay",
+            )
+        )
+        g.layout.xaxis.title = "Chemical shift (eV)"
+        g.layout.xaxis.autorange = "reversed"
+        #
+        spectra = xps_spectra_broadening(
+            chemical_shifts, equivalent_sites_data, gamma=gamma.value, sigma=sigma.value
+        )
+        # only plot the selected spectrum
+        for site, d in spectra[spectrum_select.value].items():
+            g.add_scatter(x=d[0], y=d[1], fill="tozeroy", name=site.replace("_", " "))
+
+        def response(change):
+            data = []
+            if spectra_type.value == "chemical_shift":
+                points = chemical_shifts
+                xaxis = "Chemical Shift (eV)"
+            else:
+                points = binding_energies
+                xaxis = "Binding Energy (eV)"
+            #
+            spectra = xps_spectra_broadening(
+                points, equivalent_sites_data, gamma=gamma.value, sigma=sigma.value
+            )
+
+            for site, d in spectra[spectrum_select.value].items():
+                data.append(
+                    {
+                        "x": d[0],
+                        "y": d[1],
+                        "site": site,
+                    }
+                )
+            fill_type = "tozeroy" if fill.value else None
+            with g.batch_update():
+                if len(g.data) == len(data):
+                    for i in range(len(data)):
+                        g.data[i].x = data[i]["x"]
+                        g.data[i].y = data[i]["y"]
+                        g.data[i].fill = fill_type
+                        g.data[i].name = data[i]["site"].replace("_", " ")
+
+                else:
+                    g.data = []
+                    for d in data:
+                        g.add_scatter(
+                            x=d["x"], y=d["y"], fill=fill_type, name=d["site"]
+                        )
+                g.layout.barmode = "overlay"
+                g.layout.xaxis.title = xaxis
+
+        spectra_type.observe(response, names="value")
+        spectrum_select.observe(response, names="value")
+        gamma.observe(response, names="value")
+        sigma.observe(response, names="value")
+        fill.observe(response, names="value")
+        correction_energies_table = self._get_corrections()
+        self.children = [
+            spectra_type,
+            ipw.HBox(
+                children=[
+                    spectrum_select_prompt,
+                    spectrum_select,
+                ]
+            ),
+            voigt_profile_help,
+            paras,
+            fill,
+            g,
+            correction_energies_table,
+        ]
+
+    def _get_corrections(self):
+        from aiida.orm.utils.serialize import deserialize_unsafe
+
+        ui_parameters = self.node.base.extras.get("ui_parameters", {})
+        ui_parameters = deserialize_unsafe(ui_parameters)
+        correction_energies = ui_parameters["xps"]["correction_energies"]
+        # create a table for the correction energies using ipywidgets
+        # These offsets mainly depend on chemical element and core level, and are determined by comparing the calculated core electron binding energy to the experimental one.
+        correction_energies_table = ipw.HTML(
+            """<h4>Offset Energies (δ)</h4>
+            <div>When comparing the calculated binding energies to the experimental data, these should be corrected by the given offset listed below.</div>
+            <table><tr><th>Core-level</th><th>Value (eV)</th></tr>"""
+        )
+        for core_level, value in correction_energies.items():
+            element = core_level.split("_")[0]
+            if element not in self.spectrum_select_options:
+                continue
+            exp = -value["exp"]
+            sign = "" if exp < 0 else "+"
+            correction_energies_table.value += (
+                f"<tr><td>{core_level}</td><td>{sign}{exp}</td></tr>"
+            )
+        return correction_energies_table