Source code for aiida_quantumespresso.calculations.functions.xspectra.get_xps_spectra

# -*- coding: utf-8 -*-
"""CalcFunction to compute the spectrum from ``XpsWorkchain``."""
from aiida import orm
from aiida.engine import calcfunction
import numpy as np


@calcfunction
def get_spectra_by_element(elements_list, equivalent_sites_data, voight_gamma, voight_sigma, **kwargs):  # pylint: disable=too-many-statements
    """Generate the XPS spectra for each element.

    Calculate the core level shift and binding energy for each element.
    Generate the final spectra using the Voigt profile.

    :param elements_list: a List object defining the list of elements to consider
            when producing spectrum.
    :param equivalent_sites_data: an Dict object containing symmetry data.
    :param voight_gamma: a Float node for the gamma parameter of the voigt profile.
    :param voight_sigma: a Float node for the sigma parameter of the voigt profile.
    :param structure: the StructureData object to be analysed
    :returns: Dict objects for all generated spectra and associated binding energy
            and core level shift.

    """
    from scipy.special import voigt_profile  # pylint: disable=no-name-in-module

    ground_state_node = kwargs.pop('ground_state', None)
    correction_energies = kwargs.pop('correction_energies', orm.Dict()).get_dict()
    incoming_param_nodes = {key: value for key, value in kwargs.items() if key != 'metadata'}
    group_state_energy = None
    if ground_state_node is not None:
        group_state_energy = ground_state_node.get_dict()['energy']
    elements = elements_list.get_list()
    sigma = voight_sigma.value
    gamma = voight_gamma.value
    equivalency_data = equivalent_sites_data.get_dict()

    data_dict = {element: {} for element in elements}
    for key in incoming_param_nodes:
        xspectra_out_params = incoming_param_nodes[key].get_dict()
        multiplicity = equivalency_data[key]['multiplicity']
        element = equivalency_data[key]['symbol']
        total_energy = xspectra_out_params['energy']
        data_dict[element][key] = {'element': element, 'multiplicity': multiplicity, 'total_energy': total_energy}

    result = {}
    core_level_shifts = {}
    binding_energies = {}
    for element in elements:
        spectra_list = []
        for key in data_dict[element]:
            site_multiplicity = data_dict[element][key]['multiplicity']
            spectra_list.append((site_multiplicity, float(data_dict[element][key]['total_energy']), key))
        spectra_list.sort(key=lambda entry: entry[1])
        lowest_total_energy = spectra_list[0][1]
        core_level_shift = [(entry[0], entry[1] - lowest_total_energy, entry[2]) for entry in spectra_list]
        core_level_shifts[element] = core_level_shift
        result[f'{element}_cls'] = orm.Dict(dict={entry[2]: entry[1] for entry in core_level_shift})

        if group_state_energy is not None:
            binding_energy = [(entry[0], entry[1] - group_state_energy + correction_energies[element], entry[2])
                              for entry in spectra_list]
            binding_energies[element] = binding_energy
            result[f'{element}_be'] = orm.Dict(dict={entry[2]: entry[1] for entry in binding_energy})

    fwhm_voight = gamma / 2 + np.sqrt(gamma**2 / 4 + sigma**2)

    def spectra_broadening(points, label='cls_spectra'):
        """Broadening base on the binding energy."""
        result_spectra = {}
        for element in elements:
            final_spectra_y_arrays = []
            final_spectra_y_labels = []
            final_spectra_y_units = []

            total_multiplicity = sum(i[0] for i in points[element])

            final_spectra = orm.XyData()
            max_core_level_shift = points[element][-1][1]
            min_core_level_shift = points[element][0][1]
            # 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 atoms, index in zip(points[element], range(len(points[element]))):
                # Weight for the spectra of every atom
                intensity = atoms[0]
                relative_peak_position = atoms[1]
                final_spectra_y_labels.append(f'{element}{index}_xps')
                final_spectra_y_units.append('sigma')
                final_spectra_y_arrays.append(
                    intensity * voigt_profile(x_energy_range - relative_peak_position, sigma, gamma) /
                    total_multiplicity
                )

            final_spectra_y_labels.append(f'{element}_total_xps')
            final_spectra_y_units.append('sigma')
            final_spectra_y_arrays.append(sum(final_spectra_y_arrays))

            final_spectra_x_label = 'energy'
            final_spectra_x_units = 'eV'
            final_spectra_x_array = x_energy_range
            final_spectra.set_x(final_spectra_x_array, final_spectra_x_label, final_spectra_x_units)
            final_spectra.set_y(final_spectra_y_arrays, final_spectra_y_labels, final_spectra_y_units)
            result_spectra[f'{element}_{label}'] = final_spectra
        return result_spectra

    result.update(spectra_broadening(core_level_shifts))
    if ground_state_node is not None:
        result.update(spectra_broadening(binding_energies, label='be_spectra'))
    return result