Add Feature: XpsWorkChain (#879)

Adds a WorkChain to compute XPS spectra using core-hole pseudopotentials within pw.x. Performs the following steps:
- Relaxes the input structure (optional)
- Calls get_xspectra_structures() to get all symmetrically-nonequivalent sites for each element and prepare structures suitable for SCF calculations (ground state and core-hole)
- Runs all SCF calculations
- Calls get_spectra_by_element() to return an XPS spectrum on a relative energy axis for each structure
- Returns the broadened spectrum using a Voight profile, the chemical shift values (differences in total energy relative to the lowest value) and binding energies for each element.

---------

Co-authored-by: superstar54 <xingwang1991@gmail.com>

pyproject.toml

@@ -117,6 +117,7 @@ aiida-quantumespresso = 'aiida_quantumespresso.cli:cmd_root'
 'quantumespresso.matdyn.base' = 'aiida_quantumespresso.workflows.matdyn.base:MatdynBaseWorkChain'
 'quantumespresso.pdos' = 'aiida_quantumespresso.workflows.pdos:PdosWorkChain'
 'quantumespresso.xspectra.base' = 'aiida_quantumespresso.workflows.xspectra.base:XspectraBaseWorkChain'
+'quantumespresso.xps' = 'aiida_quantumespresso.workflows.xps:XpsWorkChain'
 'quantumespresso.xspectra.core' = 'aiida_quantumespresso.workflows.xspectra.core:XspectraCoreWorkChain'
 'quantumespresso.xspectra.crystal' = 'aiida_quantumespresso.workflows.xspectra.crystal:XspectraCrystalWorkChain'
 

src/aiida_quantumespresso/calculations/functions/xspectra/get_xps_spectra.py

@@ -0,0 +1,112 @@
+# -*- 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

src/aiida_quantumespresso/workflows/protocols/xps.yaml

@@ -0,0 +1,13 @@
+default_inputs:
+ clean_workdir: True
+ abs_atom_marker: X
+ voight_gamma: 0.3
+ voight_sigma: 0.3
+default_protocol: moderate
+protocols:
+ moderate:
+ description: 'Protocol to perform an XPS calculation at normal precision at moderate computational cost.'
+ precise:
+ description: 'Protocol to perform an XPS calculation at high precision at higher computational cost.'
+ fast:
+ description: 'Protocol to perform an XPS calculation at low precision at minimal computational cost for testing purposes.'