From 1950c345c3f82c3d06c3efbeea8252c4891f3666 Mon Sep 17 00:00:00 2001 From: Peter Gillespie <55498719+PNOGillespie@users.noreply.github.com> Date: Mon, 29 Jan 2024 22:33:25 +0100 Subject: [PATCH] Add Feature: XAS Plugin (#580) Adds a plugin to calculate XANES spectra using the XspectraCrystalWorkChain of AiiDA-QE. The Setting panel allows users to simply select which element to calculate XAS for, as well as tuning the supercell size to be used in the calculation. This panel also provides the option to change the type of core-hole approximation between FCH (removing the excited electron with tot_charge = 1) and XCH (allowing the excited electron to occupy the conduction band. The result will display both the complete powder XAS for the input structure and the contributions from each symmetrically non-equivalent site. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Xing Wang --- setup.cfg | 1 + src/aiidalab_qe/plugins/xas/__init__.py | 31 ++ src/aiidalab_qe/plugins/xas/pseudo_toc.yaml | 26 + src/aiidalab_qe/plugins/xas/result.py | 535 ++++++++++++++++++++ src/aiidalab_qe/plugins/xas/setting.py | 341 +++++++++++++ src/aiidalab_qe/plugins/xas/workchain.py | 111 ++++ tests/test_plugins_xas.py | 37 ++ 7 files changed, 1082 insertions(+) create mode 100644 src/aiidalab_qe/plugins/xas/__init__.py create mode 100644 src/aiidalab_qe/plugins/xas/pseudo_toc.yaml create mode 100644 src/aiidalab_qe/plugins/xas/result.py create mode 100644 src/aiidalab_qe/plugins/xas/setting.py create mode 100644 src/aiidalab_qe/plugins/xas/workchain.py create mode 100644 tests/test_plugins_xas.py diff --git a/setup.cfg b/setup.cfg index c16c2cf8f..8965d9585 100644 --- a/setup.cfg +++ b/setup.cfg @@ -56,6 +56,7 @@ aiidalab_qe.properties = bands = aiidalab_qe.plugins.bands:bands pdos = aiidalab_qe.plugins.pdos:pdos electronic_structure = aiidalab_qe.plugins.electronic_structure:electronic_structure + xas = aiidalab_qe.plugins.xas:xas [aiidalab] title = Quantum ESPRESSO diff --git a/src/aiidalab_qe/plugins/xas/__init__.py b/src/aiidalab_qe/plugins/xas/__init__.py new file mode 100644 index 000000000..1ca6f4473 --- /dev/null +++ b/src/aiidalab_qe/plugins/xas/__init__.py @@ -0,0 +1,31 @@ +from importlib import resources + +import yaml +from aiidalab_widgets_base import ComputationalResourcesWidget + +from aiidalab_qe.common.panel import OutlinePanel +from aiidalab_qe.plugins import xas as xas_folder + +from .result import Result +from .setting import Setting +from .workchain import workchain_and_builder + +PSEUDO_TOC = yaml.safe_load(resources.read_text(xas_folder, "pseudo_toc.yaml")) + + +class XasOutline(OutlinePanel): + title = "X-ray absorption spectroscopy (XAS)" + help = """""" + + +xs_code = ComputationalResourcesWidget( + description="xspectra.x", default_calc_job_plugin="quantumespresso.xspectra" +) + +xas = { + "outline": XasOutline, + "code": {"xspectra": xs_code}, + "setting": Setting, + "result": Result, + "workchain": workchain_and_builder, +} diff --git a/src/aiidalab_qe/plugins/xas/pseudo_toc.yaml b/src/aiidalab_qe/plugins/xas/pseudo_toc.yaml new file mode 100644 index 000000000..8c3299225 --- /dev/null +++ b/src/aiidalab_qe/plugins/xas/pseudo_toc.yaml @@ -0,0 +1,26 @@ +--- +xas_xch_elements: [C, Li] +pseudos: + pbe: + gipaw_pseudos: + C: C.pbe-n-kjgipaw_psl.1.0.0.UPF + Cu: Cu.pbe-n-van_gipaw.UPF + F: F.pbe-gipaw_kj_no_hole.UPF + Li: Li.pbe-s-rrkjus-gipaw.UPF + O: O.pbe-n-kjpaw_gipaw.UPF + Si: Si.pbe-van_gipaw.UPF + core_wavefunction_data: + C: C.pbe-n-kjgipaw_psl.1.0.0.dat + Cu: Cu.pbe-n-van_gipaw.dat + F: F.pbe-gipaw_kj_no_hole.dat + Li: Li.pbe-s-rrkjus-gipaw.dat + O: O.pbe-n-kjpaw_gipaw.dat + Si: Si.pbe-van_gipaw.dat + core_hole_pseudos: + 1s: + C: C.star1s.pbe-n-kjgipaw_psl.1.0.0.UPF + Cu: Cu.star1s-pbe-n-van_gipaw.UPF + F: F.star1s-pbe-gipaw_kj.UPF + Li: Li.star1s-pbe-s-rrkjus-gipaw-test_2.UPF + O: O.star1s.pbe-n-kjpaw_gipaw.UPF + Si: Si.star1s-pbe-van_gipaw.UPF diff --git a/src/aiidalab_qe/plugins/xas/result.py b/src/aiidalab_qe/plugins/xas/result.py new file mode 100644 index 000000000..c6874c09d --- /dev/null +++ b/src/aiidalab_qe/plugins/xas/result.py @@ -0,0 +1,535 @@ +"""XAS results view widgets + +""" +import base64 +import hashlib +from typing import Callable + +import ipywidgets as ipw +import numpy as np +from IPython.display import HTML, display +from scipy.interpolate import make_interp_spline + +from aiidalab_qe.common.panel import ResultPanel + + +class SpectrumDownloadButton(ipw.Button): + """Download button with dynamic content + The content is generated using a callback when the button is clicked. + Modified from responses to https://stackoverflow.com/questions/61708701/how-to-download-a-file-using-ipywidget-button#62641240 + """ + + def __init__(self, filename: str, contents: Callable[[], str], **kwargs): + super(SpectrumDownloadButton, self).__init__(**kwargs) + self.filename = filename + self.contents = contents + self.on_click(self.__on_click) + + def __on_click(self, b): + if self.contents is None: + pass # to avoid a crash because NoneType obviously can't be processed here + else: + contents: bytes = self.contents().encode("utf-8") + b64 = base64.b64encode(contents) + payload = b64.decode() + digest = hashlib.md5(contents).hexdigest() # bypass browser cache + id = f"dl_{digest}" + + display( + HTML( + f""" + + + + + + + + """ + ) + ) + + +def write_csv(dataset): + from pandas import DataFrame + + x_vals = dataset[0]["x"] + df_data = {"energy_ev": x_vals} + for entry in dataset: + if "site" in entry["name"]: + if entry["weighting"] != 1: + df_data[ + f'{entry["name"].capitalize().replace("_", " ")} (Weighted)' + ] = entry["y"] + df_data[ + f'{entry["name"].capitalize().replace("_", " ")} (Unweighted)' + ] = entry["y"] / entry["weighting"] + else: + df_data[entry["name"].capitalize().replace("_", " ")] = entry["y"] + else: + df_data[entry["name"]] = entry["y"] + + df = DataFrame(data=df_data) + df_energy_indexed = df.set_index("energy_ev") + + return df_energy_indexed.to_csv(header=True) + + +def export_xas_data(outputs): + if "final_spectra" in outputs.xas: + final_spectra = outputs.xas.final_spectra + symmetry_analysis_data = outputs.xas.symmetry_analysis_data.get_dict() + equivalent_sites_data = symmetry_analysis_data["equivalent_sites_data"] + + return ( + final_spectra, + equivalent_sites_data, + ) + else: + return None + + +def broaden_xas( + input_array, variable=False, gamma_hole=0.01, gamma_max=5, center_energy=15 +): + """Take an input spectrum and return a broadened spectrum as output using either a constant or variable parameter. + + :param input_array: The 2D array of x/y values to be broadened. Should be plotted with + little or no broadening before using the function. + :param gamma_hole: The broadening parameter for the Lorenzian broadening function. In constant mode (variable=False), + this value is applied to the entire spectrum. In variable mode (variable=True), this value defines + the starting broadening parameter of the arctangent function. Refers to the natural linewidth of + the element/XAS edge combination in question and (for elements Z > 10) should be based on reference + values from X-ray spectroscopy. + :param variable: Request a variable-energy broadening of the spectrum instead of the defaultconstant broadening. + Uses the functional form defined in Calandra and Bunau, PRB, 87, 205105 (2013). + :param gamma_max: The maximum lorenzian broadening to be applied in variable energy broadening mode. Refers to the + broadening effects at infinite energy above the main edge. + :param center_energy: The inflection point of the variable broadening function. Does not relate to experimental data + and must be tuned manually. + """ + + if variable: + if not all([gamma_hole, gamma_max, center_energy]): + missing = [ + i[0] + for i in zip( + ["gamma_hole", "gamma_max", "center_energy"], + [gamma_hole, gamma_max, center_energy], + ) + if i[1] is None + ] + raise ValueError( + f"The following variables were not defined {missing} and are required for variable-energy broadening" + ) + + x_vals = input_array[:, 0] + y_vals = input_array[:, 1] + + lorenz_y = np.zeros(len(x_vals)) + + if variable: + for x, y in zip(x_vals, y_vals): + if x < 0: # the function is bounded between gamma_hole and gamma_max + gamma_var = gamma_hole + else: + e = x / center_energy + + gamma_var = gamma_hole + gamma_max * ( + 0.5 + np.arctan((e - 1) / (e**2)) / np.pi + ) + + if y <= 1.0e-6: # do this to skip the calculation for very small values + lorenz_y = y + else: + lorenz_y += ( + gamma_var + / 2.0 + / np.pi + / ((x_vals - x) ** 2 + 0.25 * gamma_var**2) + * y + ) + else: + for x, y in zip(x_vals, y_vals): + lorenz_y += ( + gamma_hole + / 2.0 + / np.pi + / ((x_vals - x) ** 2 + 0.25 * gamma_hole**2) + * y + ) + + return np.column_stack((x_vals, lorenz_y)) + + +def get_aligned_spectra(core_wc_dict, equivalent_sites_dict): + """Return a set of spectra aligned according to the chemical shift (difference in Fermi level). + + Primarily this is a copy of ``get_spectra_by_element`` from AiiDA-QE which operates on only one + element. + """ + data_dict = {} + spectrum_dict = { + site: node.outputs.powder_spectrum for site, node in core_wc_dict.items() + } + for key, value in core_wc_dict.items(): + xspectra_out_params = value.outputs.parameters_xspectra__xas_0.get_dict() + energy_zero = xspectra_out_params["energy_zero"] + multiplicity = equivalent_sites_dict[key]["multiplicity"] + + if "total_multiplicity" not in data_dict: + data_dict["total_multiplicity"] = multiplicity + else: + data_dict["total_multiplicity"] += multiplicity + + data_dict[key] = { + "spectrum_node": spectrum_dict[key], + "multiplicity": multiplicity, + "energy_zero": energy_zero, + } + + spectra_list = [] + total_multiplicity = data_dict.pop("total_multiplicity") + for key in data_dict: + spectrum_node = data_dict[key]["spectrum_node"] + site_multiplicity = data_dict[key]["multiplicity"] + weighting = site_multiplicity / total_multiplicity + weighting_string = f"{site_multiplicity}/{total_multiplicity}" + spectrum_x = spectrum_node.get_x()[1] + spectrum_y = spectrum_node.get_y()[0][1] + spline = make_interp_spline(spectrum_x, spectrum_y) + norm_y = spline(spectrum_x) / np.trapz(spline(spectrum_x), spectrum_x) + weighted_spectrum = np.column_stack( + (spectrum_x, norm_y * (site_multiplicity / total_multiplicity)) + ) + spectra_list.append( + ( + weighted_spectrum, + key, + weighting, + weighting_string, + float(data_dict[key]["energy_zero"]), + ) + ) + + # Sort according to Fermi level, then correct to align all spectra to the + # highest value. Note that this is needed because XSpectra automatically aligns the + # final spectrum such that the system's Fermi level is at 0 eV. + spectra_list.sort(key=lambda entry: entry[-1]) + highest_level = spectra_list[0][-1] + energy_zero_corrections = [ + (entry[0], entry[1], entry[2], entry[3], entry[-1] - highest_level) + for entry in spectra_list + ] + aligned_spectra = [ + ( + entry[1], + entry[2], + entry[3], + np.column_stack((entry[0][:, 0] - entry[-1], entry[0][:, 1])), + ) + for entry in energy_zero_corrections + ] + + return aligned_spectra + + +class Result(ResultPanel): + title = "XAS" + workchain_labels = ["xas"] + + def __init__(self, node=None, **kwargs): + super().__init__(node=node, identifier="xas", **kwargs) + + def _update_view(self): + import plotly.graph_objects as go + + gamma_select_prompt = ipw.HTML( + """ +
+ Select parameters for spectrum broadening
""" + ) + + # PNOG: If someone knows a way to format certain words differently in HTML without causing a line-break, hit me up. + # For now, (17/10/23) we'll just have the function terms be in italics. + # Alternatively, if it's possible to format mathematical terms in HTML without a line-break, let me know + variable_broad_select_help = ipw.HTML( + """ +
+ Broadening parameters: + +

Γhole - Defines a constant Lorenzian broadening width for the whole spectrum. In "variable" mode, defines the initial broadening width of the ArcTangent function.

+

Γmax - Maximum Lorenzian broadening parameter at infinte energy in "variable" mode.

+

Ecenter - Defines the inflection point of the variable-energy broadening function.

+
+
+ Note that setting Γhole to 0 eV will simply plot the raw spectrum. +
+ """ + ) + spectrum_select_prompt = ipw.HTML( + """ +
+ Select spectrum to plot
""" + ) + final_spectra, equivalent_sites_data = export_xas_data(self.outputs) + xas_wc = [ + n for n in self.node.called if n.process_label == "XspectraCrystalWorkChain" + ][0] + core_wcs = { + n.get_metadata_inputs()["metadata"]["call_link_label"]: n + for n in xas_wc.called + if n.process_label == "XspectraCoreWorkChain" + } + core_wc_dict = { + key.replace("_xspectra", ""): value for key, value in core_wcs.items() + } + + spectrum_select_options = [key.split("_")[0] for key in final_spectra.keys()] + + spectrum_select = ipw.Dropdown( + description="", + disabled=False, + value=spectrum_select_options[0], + options=spectrum_select_options, + layout=ipw.Layout(width="20%"), + ) + + variable_broad_select = ipw.Checkbox( + value=False, + disabled=False, + description="Use variable energy broadening.", + style={"description_width": "initial", "opacity": 0.5}, + ) + + gamma_hole_select = ipw.FloatSlider( + value=0.0, + min=0.0, + max=5, + step=0.1, + description="$\Gamma_{hole}$", # noqa: W605 + disabled=False, + continuous_update=False, + orientation="horizontal", + readout=True, + ) + + gamma_max_select = ipw.FloatSlider( + value=5.0, + min=2.0, + max=10, + step=0.5, + continuous_update=False, + description="$\Gamma_{max}$", # noqa: W605 + disabled=True, + orientation="horizontal", + readout=True, + ) + + center_e_select = ipw.FloatSlider( + value=15.0, + min=5, + max=30, + step=0.5, + continuous_update=False, + description="$E_{center}$", + disabled=True, + orientation="horizontal", + readout=True, + ) + download_data = SpectrumDownloadButton( + filename=f"{spectrum_select.value}_XAS_Spectra.csv", + contents=None, + description="Download CSV", + icon="download", + ) + # # get data + # # init figure + g = go.FigureWidget( + layout=go.Layout( + title=dict(text="XAS"), + barmode="overlay", + ) + ) + + g.layout.xaxis.title = "Relative Photon Energy (eV)" + + chosen_spectrum = spectrum_select.value + chosen_spectrum_label = f"{chosen_spectrum}_xas" + spectra = final_spectra[chosen_spectrum_label] + + raw_spectrum = np.column_stack((spectra.get_x()[1], spectra.get_y()[0][1])) + + x = raw_spectrum[:, 0] + y = raw_spectrum[:, 1] + spline = make_interp_spline(x, y) + norm_y = spline(x) / np.trapz(spline(x), x) + element = chosen_spectrum_label.split("_")[0] + element_sites = [ + key + for key in equivalent_sites_data + if equivalent_sites_data[key]["symbol"] == element + ] + element_core_wcs = {} + total_multiplicity = 0 + for site in element_sites: + site_multiplicity = equivalent_sites_data[site]["multiplicity"] + total_multiplicity += site_multiplicity + element_core_wcs[site] = core_wc_dict[site] + + g.add_scatter(x=x, y=norm_y, name=f"{element} K-edge") + for entry in get_aligned_spectra( + core_wc_dict=element_core_wcs, equivalent_sites_dict=equivalent_sites_data + ): + g.add_scatter( + x=entry[-1][:, 0], + y=entry[-1][:, 1], + name=entry[0].capitalize().replace("_", " "), + ) + + def _update_download_selection(dataset, element): + download_data.contents = lambda: write_csv(dataset) + download_data.filename = f"{element}_XAS_Spectra.csv" + + def response(change): + chosen_spectrum = spectrum_select.value + chosen_spectrum_label = f"{chosen_spectrum}_xas" + element_sites = [ + key + for key in equivalent_sites_data + if equivalent_sites_data[key]["symbol"] == chosen_spectrum + ] + element_core_wcs = { + key: value + for key, value in core_wc_dict.items() + if key in element_sites + } + spectra = [] + final_spectrum_node = final_spectra[chosen_spectrum_label] + final_spectrum = np.column_stack( + (final_spectrum_node.get_x()[1], final_spectrum_node.get_y()[0][1]) + ) + final_x_vals = final_spectrum[:, 0] + final_y_vals = final_spectrum[:, 1] + final_spectrum_spline = make_interp_spline(final_x_vals, final_y_vals) + final_norm_y = final_spectrum_spline(final_x_vals) / np.trapz( + final_spectrum_spline(final_x_vals), final_x_vals + ) + spectra.append( + ( + f"{chosen_spectrum} K-edge", + 1, + "1", + np.column_stack((final_x_vals, final_norm_y)), + ) + ) + datasets = [] + for entry in get_aligned_spectra( + core_wc_dict=element_core_wcs, + equivalent_sites_dict=equivalent_sites_data, + ): + spectra.append(entry) + + for entry in spectra: + label = entry[0] + weighting = entry[1] + weighting_string = entry[2] + raw_spectrum = entry[-1] + x = raw_spectrum[:, 0] + y = raw_spectrum[:, 1] + if not variable_broad_select: + gamma_max_select.disabled = True + center_e_select.disabled = True + else: + gamma_max_select.disabled = False + center_e_select.disabled = False + + if gamma_hole_select.value == 0.0: + x = raw_spectrum[:, 0] + y = raw_spectrum[:, 1] + else: + broad_spectrum = broaden_xas( + raw_spectrum, + gamma_hole=gamma_hole_select.value, + gamma_max=gamma_max_select.value, + center_energy=center_e_select.value, + variable=variable_broad_select.value, + ) + x = broad_spectrum[:, 0] + y = broad_spectrum[:, 1] + + final_spline = make_interp_spline(x, y) + final_y_vals = final_spline(final_x_vals) + datasets.append( + { + "x": final_x_vals, + "y": final_y_vals, + "name": label, + "weighting": weighting, + "weighting_string": weighting_string, + } + ) + _update_download_selection(datasets, chosen_spectrum) + + with g.batch_update(): + # If the number of datasets is different from one update to the next, + # then we need to reset the data already in the Widget. Otherwise, we can + # simply override the data. This also helps since then changing the + # broadening is much smoother. + if len(datasets) == len( + g.data + ): # if the number of entries is the same, just update + for index, entry in enumerate(datasets): + g.data[index].x = entry["x"] + g.data[index].y = entry["y"] + if "site_" in entry["name"]: + g.data[index].name = ( + entry["name"].capitalize().replace("_", " ") + ) + else: + g.data[index].name = entry["name"] + else: # otherwise, reset the figure + g.data = () + for entry in datasets: + if "site_" in entry["name"]: + name = entry["name"].capitalize().replace("_", " ") + else: + name = entry["name"] + g.add_scatter(x=entry["x"], y=entry["y"], name=name) + + spectrum_select.observe(response, names="value") + gamma_hole_select.observe(response, names="value") + gamma_max_select.observe(response, names="value") + center_e_select.observe(response, names="value") + variable_broad_select.observe(response, names="value") + download_data.observe(response, names=["contents", "filename"]) + self.children = [ + ipw.HBox( + [ + ipw.VBox( + [ + spectrum_select_prompt, + spectrum_select, + gamma_select_prompt, + gamma_hole_select, + gamma_max_select, + center_e_select, + ], + layout=ipw.Layout(width="40%"), + ), + ipw.VBox( + [ + variable_broad_select, + variable_broad_select_help, + ], + layout=ipw.Layout(width="60%"), + ), + ] + ), + download_data, + g, + ] diff --git a/src/aiidalab_qe/plugins/xas/setting.py b/src/aiidalab_qe/plugins/xas/setting.py new file mode 100644 index 000000000..55b6ccb36 --- /dev/null +++ b/src/aiidalab_qe/plugins/xas/setting.py @@ -0,0 +1,341 @@ +# -*- coding: utf-8 -*- +"""Panel for XAS plugin. + +""" +import os +import tarfile +from importlib import resources +from pathlib import Path + +import ipywidgets as ipw +import requests +import traitlets as tl +import yaml +from aiida import orm + +from aiidalab_qe.common.panel import Panel +from aiidalab_qe.plugins import xas as xas_folder + +PSEUDO_TOC = yaml.safe_load(resources.read_text(xas_folder, "pseudo_toc.yaml")) +pseudo_data_dict = PSEUDO_TOC["pseudos"] +xch_elements = PSEUDO_TOC["xas_xch_elements"] + +base_url = "https://github.com/PNOGillespie/Core_Level_Spectra_Pseudos/raw/main" +head_path = f"{Path.home()}/.local/lib" +dir_header = "cls_pseudos" +functionals = ["pbe"] +core_wfc_dir = "core_wfc_data" +gipaw_dir = "gipaw_pseudos" +ch_pseudo_dir = "ch_pseudos/star1s" + + +def _load_or_import_nodes_from_filenames(in_dict, path, core_wfc_data=False): + for filename in in_dict.values(): + try: + orm.load_node(filename) + except BaseException: + if not core_wfc_data: + new_upf = orm.UpfData(f"{path}/{filename}", filename=filename) + new_upf.label = filename + new_upf.store() + else: + new_singlefile = orm.SinglefileData( + f"{path}/{filename}", filename="stdout" + ) + new_singlefile.label = filename + new_singlefile.store() + + +def _download_extract_pseudo_archive(func): + dir = f"{head_path}/{dir_header}/{func}" + archive_filename = f"{func}_ch_pseudos.tgz" + remote_archive_filename = f"{base_url}/{func}/{archive_filename}" + local_archive_filename = f"{dir}/{archive_filename}" + + env = os.environ.copy() + env["PATH"] = f"{env['PATH']}:{Path.home() / '.local' / 'lib'}" + + response = requests.get(remote_archive_filename, timeout=30) + response.raise_for_status() + with open(local_archive_filename, "wb") as handle: + handle.write(response.content) + handle.flush() + response.close() + + with tarfile.open(local_archive_filename, "r:gz") as tarfil: + tarfil.extractall(dir) + + +url = f"{base_url}" +for func in functionals: + dir = f"{head_path}/{dir_header}/{func}" + os.makedirs(dir, exist_ok=True) + archive_filename = f"{func}_ch_pseudos.tgz" + archive_found = False + for entry in os.listdir(dir): + if entry == archive_filename: + archive_found = True + if not archive_found: + _download_extract_pseudo_archive(func) + + +# Check all the pseudos/core-wfc data files in the TOC dictionary +# and load/check all of them before proceeding. Note that this +# approach relies on there not being multiple instances of nodes +# with the same label. +for func in functionals: + gipaw_pseudo_dict = pseudo_data_dict[func]["gipaw_pseudos"] + core_wfc_dict = pseudo_data_dict[func]["core_wavefunction_data"] + core_hole_pseudo_dict = pseudo_data_dict[func]["core_hole_pseudos"] + main_path = f"{head_path}/{dir_header}/{func}" + core_wfc_dir = f"{main_path}/core_wfc_data" + gipaw_dir = f"{main_path}/gipaw_pseudos" + ch_pseudo_dir = f"{main_path}/ch_pseudos/star1s" + # First, check that the local directories contain what's in the pseudo_toc + for pseudo_dir, pseudo_dict in zip( + [gipaw_dir, core_wfc_dir, ch_pseudo_dir], + [gipaw_pseudo_dict, core_wfc_dict, core_hole_pseudo_dict], + ): + pseudo_toc_mismatch = os.listdir(pseudo_dir) != pseudo_dict.values() + + # Re-download the relevant archive if there is a mismatch + if pseudo_toc_mismatch: + _download_extract_pseudo_archive(func) + + _load_or_import_nodes_from_filenames( + in_dict=gipaw_pseudo_dict, + path=gipaw_dir, + ) + _load_or_import_nodes_from_filenames( + in_dict=core_wfc_dict, path=core_wfc_dir, core_wfc_data=True + ) + _load_or_import_nodes_from_filenames( + in_dict=core_hole_pseudo_dict["1s"], path=ch_pseudo_dir + ) + + +class Setting(Panel): + title = "XAS Settings" + identifier = "xas" + input_structure = tl.Instance(orm.StructureData, allow_none=True) + protocol = tl.Unicode(allow_none=True) + + element_selection_title = ipw.HTML( + """
+

Element and Core-Hole Treatment Setting.

""" + ) + + # TODO: The element selection should lock the "Confirm" button if no elements have been + # selected for XAS calculation. + + element_selection_help = ipw.HTML( + """
+ To select elements for calculation of K-edge spectra:
+ (1) Tick the checkbox for each element symbol to select the element for calculation.
+ (2) Select the core-hole treatment scheme from the dropdown box.
+
+ There are three supported options for core-hole treatment:
+ - FCH: Remove one electron from the system (any occupations scheme).
+ - XCH (Smearing): places the excited electron into the conduction band (smeared occupations).
+ - XCH (Fixed): places the excited electron into the conduction band (fixed occupations).
+
+ For XAS calculations of most elements, the FCH treatment is recommended, however in some cases the XCH treatment should be used instead.
+ The recommended setting will be shown for each available element. + Note that only elements for which core-hole pseudopotential sets are available + will be shown.
+
""" + ) + # I will leave these objects here for now (15/11/23), but since the calculation of molecular + # systems is not really supported (neither in terms of XAS nor the main App itself) we should + # not present this option that essentially does nothing. + # structure_title = ipw.HTML( + # """
+ #

Structure

""" + # ) + # structure_help = ipw.HTML( + # """
+ # Below you can indicate if the material should be treated as a molecule + # or a crystal. + #
""" + # ) + supercell_title = ipw.HTML( + """
+

Cell size

""" + ) + supercell_help = ipw.HTML( + """
+ Define the minimum cell length in angstrom for the resulting supercell, and thus all output + structures. The default value of 8.0 angstrom will be used + if no input is given. Setting this value to 0.0 will + instruct the CF to not scale up the input structure. +
""" + ) + + def __init__(self, **kwargs): + self.gipaw_pseudos = pseudo_data_dict["pbe"]["gipaw_pseudos"] + self.core_hole_pseudos = pseudo_data_dict["pbe"]["core_hole_pseudos"]["1s"] + self.core_wfc_data_dict = pseudo_data_dict["pbe"]["core_wavefunction_data"] + + self.element_and_ch_treatment = ipw.VBox(layout=ipw.Layout(width="100%")) + + # self.structure_type = ipw.ToggleButtons( + # options=[ + # ("Molecule", "molecule"), + # ("Crystal", "crystal"), + # ], + # value="crystal", + # ) + self.supercell_min_parameter = ipw.FloatText( + value=8.0, + description="The minimum cell length (Å):", + disabled=False, + style={"description_width": "initial"}, + ) + + self.children = [ + # self.structure_title, + # self.structure_help, + # ipw.HBox( + # [self.structure_type], + # ), + self.element_selection_title, + self.element_selection_help, + ipw.HBox([self.element_and_ch_treatment], layout=ipw.Layout(width="95%")), + self.supercell_title, + self.supercell_help, + ipw.HBox( + [self.supercell_min_parameter], + ), + ] + + super().__init__(**kwargs) + + def get_panel_value(self): + elements_list = [] + core_hole_treatments = {} + for entry in self.element_and_ch_treatment.children: + if entry.children[0].value is True: + element = entry.children[0].description + ch_treatment = entry.children[1].value + elements_list.append(element) + core_hole_treatments[element] = ch_treatment + + pseudo_labels = {} + core_wfc_data_labels = {} + for element in elements_list: + pseudo_labels[element] = { + "gipaw": self.gipaw_pseudos[element], + "core_hole": self.core_hole_pseudos[element], + } + core_wfc_data_labels[element] = self.core_wfc_data_dict[element] + + parameters = { + "core_hole_treatments": core_hole_treatments, + "elements_list": elements_list, + # "structure_type": self.structure_type.value, + "pseudo_labels": pseudo_labels, + "core_wfc_data_labels": core_wfc_data_labels, + "supercell_min_parameter": self.supercell_min_parameter.value, + } + return parameters + + def set_panel_value(self, input_dict): + """Load a dictionary with the input parameters for the plugin.""" + + # set selected elements and core-hole treatments + elements_list = input_dict.get("elements_list", []) + for entry in self.element_and_ch_treatment.children: + element = entry.children[0].description + if element in elements_list: + entry.children[0].value = True + entry.children[1].value = input_dict["core_hole_treatments"][element] + else: + entry.children[0].value = False + entry.children[1].value = "full" + # set supercell min parameter + self.supercell_min_parameter.value = input_dict.get( + "supercell_min_parameter", 8.0 + ) + # self.structure_type.value = input_dict.get("structure_type", "crystal") + + @tl.observe("input_structure") + def _update_structure(self, _=None): + self._update_element_select_panel() + + def _update_element_select_panel(self): + if self.input_structure is None: + return + + starting_treatment_mapping = {"FCH": "full", "XCH": "xch_smear"} + ch_treatment_options = [ + ("FCH", "full"), + ("XCH (Smearing)", "xch_smear"), + ("XCH (Fixed)", "xch_fixed"), + ] + ch_pseudos = self.core_hole_pseudos + structure = self.input_structure + available_elements = [k for k in ch_pseudos] + elements_to_select = sorted( + [ + kind.symbol + for kind in structure.kinds + if kind.symbol in available_elements + ] + ) + treatment_options = () + + for element in elements_to_select: + if element in xch_elements: + recommended_treatment = "XCH" + else: + recommended_treatment = "FCH" + + treatment_options += ( + ipw.HBox( + [ + ipw.Checkbox( + description=element, + value=False, + disabled=False, + style={"description_width": "initial"}, + layout=ipw.Layout(width="7%"), + ), + ipw.Dropdown( + options=ch_treatment_options, + value=starting_treatment_mapping[recommended_treatment], + disabled=False, + layout=ipw.Layout(width="15%"), + ), + ipw.HTML( + f"Recommended treatment: {recommended_treatment} (PBE Core-Hole Pseudopotential)", + layout=ipw.Layout(width="78%"), + ), + ], + layout=ipw.Layout( + width="100%", + ), + ), + ) + + self.element_and_ch_treatment.children = treatment_options + + # For reference: + # This is the whole widget: + # print(f"{self.element_and_ch_treatment}\n") + + # This is the tuple of selected element and core-hole treatment: + # print(f"{self.element_and_ch_treatment.children[0]}\n") + + # This is the checkbox for the element, giving element name and whether to add it to the elements list + # print(f"{self.element_and_ch_treatment.children[0].children[0]}\n") + # print(f"{self.element_and_ch_treatment.children[0].children[0].value}\n") + # print(f"{self.element_and_ch_treatment.children[0].children[0].description}\n") + + # This is the dropdown for the core-hole treatment option: + # print(f"{self.element_and_ch_treatment.children[0].children[1]}\n") + # print(f"{self.element_and_ch_treatment.children[0].children[1].value}\n") + + def reset(self): + """Reset the panel to its initial state.""" + self.input_structure = None + # self.structure_type.value = "crystal" diff --git a/src/aiidalab_qe/plugins/xas/workchain.py b/src/aiidalab_qe/plugins/xas/workchain.py new file mode 100644 index 000000000..bfd175053 --- /dev/null +++ b/src/aiidalab_qe/plugins/xas/workchain.py @@ -0,0 +1,111 @@ +from importlib import resources + +import yaml +from aiida import orm +from aiida.plugins import WorkflowFactory +from aiida_quantumespresso.common.types import ElectronicType, SpinType + +from aiidalab_qe.plugins import xas as xas_folder + +XspectraCrystalWorkChain = WorkflowFactory("quantumespresso.xspectra.crystal") +PSEUDO_TOC = yaml.safe_load(resources.read_text(xas_folder, "pseudo_toc.yaml")) +pseudo_data_dict = PSEUDO_TOC["pseudos"] +xch_elements = PSEUDO_TOC["xas_xch_elements"] + + +def get_builder(codes, structure, parameters, **kwargs): + from copy import deepcopy + + protocol = parameters["workchain"]["protocol"] + xas_parameters = parameters["xas"] + core_hole_treatments = xas_parameters["core_hole_treatments"] + elements_list = xas_parameters["elements_list"] + supercell_min_parameter = xas_parameters["supercell_min_parameter"] + pseudo_labels = xas_parameters["pseudo_labels"] + core_wfc_data_labels = xas_parameters["core_wfc_data_labels"] + pseudos = {} + # Convert the pseudo and core_wfc_data node labels into nodes: + core_wfc_data = {k: orm.load_node(v) for k, v in core_wfc_data_labels.items()} + for element in elements_list: + pseudos[element] = { + k: orm.load_node(v) for k, v in pseudo_labels[element].items() + } + + # TODO should we override the cutoff_wfc, cutoff_rho by the new pseudo? + # In principle we should, if we know what that value is, but that would + # require testing them first... + + # (13/10/23) I'm keeping the part about molecules in for future reference, + # but we need to establish the protocol & backend code for XAS of molecules + # before thinking about a workflow. + # (22/01/24) Commented out the code for molecules, just so the option doesn't + # appear in the UI and confuse the user. + # is_molecule_input = ( + # True if xas_parameters.get("structure_type") == "molecule" else False + # ) + + structure_preparation_settings = { + "supercell_min_parameter": orm.Float(supercell_min_parameter), + # "is_molecule_input": orm.Bool(is_molecule_input), + } + spglib_settings = orm.Dict({"symprec": 1.0e-3}) + + pw_code = codes["pw"] + xs_code = codes["xspectra"] + overrides = { + "core": { + "scf": deepcopy(parameters["advanced"]), + # PG: Here, we set a "variable" broadening scheme, which actually defines a constant broadening + # The reason for this is that in "gamma_mode = constant", the Lorenzian broadening parameter + # is defined by "xgamma" (in "PLOT"), but this parameter *also* controls the broadening value + # used in the Lanczos algorithm to enhance the convergence rate. In order to give the user a + # final spectrum with minimal broadening, we use "gamma_mode = variable", which uses a different + # parameter set ("gamma_energy(1-2)", "gamma_value(1-2)") and thus allows us to decouple spectrum + # broadening from Lanczos broadening and avoid having to re-plot the final spectrum. + "xs_prod": { + "xspectra": { + "parameters": { + "PLOT": { + "gamma_mode": "variable", + "gamma_energy(1)": 0, + "gamma_energy(2)": 1, + "gamma_value(1)": 0.1, + "gamma_value(2)": 0.1, + } + } + } + }, + } + } + + builder = XspectraCrystalWorkChain.get_builder_from_protocol( + pw_code=pw_code, + xs_code=xs_code, + structure=structure, + protocol=protocol, + pseudos=pseudos, + elements_list=elements_list, + core_hole_treatments=core_hole_treatments, + core_wfc_data=core_wfc_data, + electronic_type=ElectronicType(parameters["workchain"]["electronic_type"]), + spin_type=SpinType(parameters["workchain"]["spin_type"]), + # TODO: We will need to merge the changes in AiiDA-QE PR#969 in order + # to better handle magnetic and Hubbard data. For now, we can probably + # leave it as it is. + initial_magnetic_moments=parameters["advanced"]["initial_magnetic_moments"], + overrides=overrides, + **kwargs, + ) + builder.pop("relax") + builder.pop("clean_workdir", None) + builder.spglib_settings = spglib_settings + builder.structure_preparation_settings = structure_preparation_settings + + return builder + + +workchain_and_builder = { + "workchain": XspectraCrystalWorkChain, + "exclude": ("clean_workdir", "structure", "relax"), + "get_builder": get_builder, +} diff --git a/tests/test_plugins_xas.py b/tests/test_plugins_xas.py new file mode 100644 index 000000000..4b45fc43a --- /dev/null +++ b/tests/test_plugins_xas.py @@ -0,0 +1,37 @@ +import pytest + + +@pytest.mark.usefixtures("sssp") +def test_settings(submit_app_generator): + """Test the settings of the xas app.""" + app = submit_app_generator(properties=["xas"]) + configure_step = app.configure_step + # test get_panel_value + # select the first elmement + configure_step.settings["xas"].element_and_ch_treatment.children[0].children[ + 0 + ].value = True + configure_step.settings["xas"].supercell_min_parameter.value = 4.0 + parameters = configure_step.settings["xas"].get_panel_value() + assert parameters["core_hole_treatments"] == {"Si": "full"} + assert parameters["pseudo_labels"] == { + "Si": { + "gipaw": "Si.pbe-van_gipaw.UPF", + "core_hole": "Si.star1s-pbe-van_gipaw.UPF", + } + } + assert parameters["core_wfc_data_labels"] == {"Si": "Si.pbe-van_gipaw.dat"} + assert parameters["supercell_min_parameter"] == 4.0 + # test set_panel_value + # update the parameters + parameters["supercell_min_parameter"] = 5.0 + parameters["core_hole_treatments"] = {"Si": "xch_smear"} + configure_step.settings["xas"].set_panel_value(parameters) + assert configure_step.settings["xas"].supercell_min_parameter.value == 5.0 + assert ( + configure_step.settings["xas"] + .element_and_ch_treatment.children[0] + .children[1] + .value + == "xch_smear" + )