Add the XspectraCoreWorkChain (#872)

The `XspectraCoreWorkChain` takes a given input structure in which one
atomic site is designated as the site for the absorbing atom and run all
`XspectraCalculation`s as requested. The results are post-processed into
a powder spectrum. This WorkChain makes use of as many features already
implemented in `aiida-quantumespresso` as possible, however, it does
require the `aiida-shell` plugin to run `upf2plotcore.sh` to generate
the `core_wfc_data` input.

pyproject.toml

@@ -115,6 +115,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.xspectra.core' = 'aiida_quantumespresso.workflows.xspectra.core:XspectraCoreWorkChain'
 
 [tool.flit.module]
 name = 'aiida_quantumespresso'

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

@@ -0,0 +1,102 @@
+# -*- coding: utf-8 -*-
+"""CalcFunction to compute the powder spectrum of a set of XANES spectra from ``XspectraCalculation``."""
+from aiida.common import ValidationError
+from aiida.engine import calcfunction
+from aiida.orm import XyData
+
+
+@calcfunction
+def get_powder_spectrum(**kwargs): # pylint: disable=too-many-statements
+ """Combine the given spectra into a single "Powder" spectrum, representing the XAS of a powder sample.
+
+ The function expects between 1 and 3 XyData nodes from ``XspectraCalculation`` whose
+ polarisation vectors are the basis vectors of the original crystal structure (100, 010, 001).
+ """
+ spectra = [node for node in kwargs.values() if isinstance(node, XyData)]
+ vectors = [node.creator.res['xepsilon'] for node in spectra]
+
+ if len(vectors) > 3:
+ raise ValidationError(f'Expected between 1 and 3 XyData nodes as input, but {len(spectra)} were given.')
+
+ # If the system is isochoric (e.g. a cubic system) then the three basis vectors are
+ # equal to each other, thus we simply return the
+ if len(vectors) == 1:
+ vector = vectors[0]
+ vector_string = f'{float(vector[0])} {float(vector[1])} {float(vector[2])}'
+ if vector not in [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]:
+ raise ValidationError(
+ f'Polarisation vector ({vector_string}) does not correspond to a crystal basis vector (100, 010, 001)'
+ )
+
+ powder_spectrum = spectra[0]
+ powder_x = powder_spectrum.get_x()[1]
+ powder_y = powder_spectrum.get_y()[0][1]
+
+ powder_data = XyData()
+ powder_data.set_x(powder_x, 'energy', 'eV')
+ powder_data.set_y(powder_y, 'sigma', 'n/a')
+
+ # if the system is dichoric (e.g. a hexagonal system) then the A and B periodic
+ # dimensions are equal to each other by symmetry, thus the powder spectrum is simply
+ # the average of 2x the 1 0 0 eps vector and 1x the 0 0 1 eps vector
+ if len(vectors) == 2:
+ # Since the individual vectors are labelled, we can extract just the spectra needed
+ # to produce the powder and leave the rest
+
+ spectrum_a = None
+ spectrum_c = None
+ for vector, spectrum in zip(vectors, spectra):
+ vector_string = f'{float(vector[0])} {float(vector[1])} {float(vector[2])}'
+ if vector in [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]:
+ spectrum_a = spectrum
+ elif vector == [0.0, 0.0, 1.0]:
+ spectrum_c = spectrum
+ else:
+ raise ValidationError(
+ f'Polarisation vector ({vector_string}) does not correspond to a crystal basis vector'
+ ' (100, 010, 001)'
+ )
+
+ if spectrum_a and not spectrum_c:
+ raise ValidationError(f'Found no polarisation vector for the C-axis ([0, 0, 1]), found instead: {vectors}.')
+ powder_x = spectrum_a.get_x()[1]
+ yvals_a = spectrum_a.get_y()[0][1]
+ yvals_c = spectrum_c.get_y()[0][1]
+
+ powder_y = ((yvals_a * 2) + yvals_c) / 3
+ powder_data = XyData()
+ powder_data.set_x(powder_x, 'energy', 'eV')
+ powder_data.set_y(powder_y, 'sigma', 'n/a')
+
+ # if the system is trichoric (e.g. a monoclinic system) then no periodic dimensions
+ # are equal by symmetry, thus the powder spectrum is the average of the three basis
+ # dipole vectors (1.0 0.0 0.0, 0.0 1.0 0.0, 0.0 0.0 1.0)
+ if len(vectors) == 3:
+ # Since the individual vectors are labelled, we can extract just the spectra needed to
+ # produce the powder and leave the rest
+ for vector, spectra in zip(vectors, spectra):
+ vector_string = f'{float(vector[0])} {float(vector[1])} {float(vector[2])}'
+ if vector == [1.0, 0.0, 0.0]:
+ spectrum_a = spectra
+ elif vector == [0.0, 1.0, 0.0]:
+ spectrum_b = spectra
+ elif vector == [0.0, 0.0, 1.0]:
+ spectrum_c = spectra
+ else:
+ raise ValidationError(
+ f'Polarisation vector ({vector_string}) does not correspond to a crystal basis vector'
+ ' (100, 010, 001)'
+ )
+
+ powder_x = spectrum_a.get_x()[1]
+ yvals_a = spectrum_a.get_y()[0][1]
+ yvals_b = spectrum_b.get_y()[0][1]
+ yvals_c = spectrum_c.get_y()[0][1]
+
+ powder_y = (yvals_a + yvals_b + yvals_c) / 3
+
+ powder_data = XyData()
+ powder_data.set_x(powder_x, 'energy', 'eV')
+ powder_data.set_y(powder_y, 'sigma', 'n/a')
+
+ return powder_data

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

@@ -0,0 +1,66 @@
+# -*- coding: utf-8 -*-
+"""CalcFunction to merge multiple ``XyData`` nodes of calculated XANES spectra into a new ``XyData`` node."""
+from aiida.engine import calcfunction
+from aiida.orm import XyData
+
+
+@calcfunction
+def merge_spectra(**kwargs):
+ """Compile all calculated spectra into a single ``XyData`` node for easier plotting.
+
+ The keyword arguments must be an arbitrary number of ``XyData`` nodes from
+ the `output_spectra` of `XspectraCalculation`s, all other `kwargs` will be discarded at
+ runtime.
+
+ Returns a single ``XyData`` node where each set of y values is labelled
+ according to the polarisation vector used for the `XspectraCalculation`.
+ """
+ output_spectra = XyData()
+ y_arrays_list = []
+ y_units_list = []
+ y_labels_list = []
+
+ spectra = [node for label, node in kwargs.items() if isinstance(node, XyData)]
+
+ for spectrum_node in spectra:
+ calc_node = spectrum_node.creator
+ calc_out_params = calc_node.res
+ eps_vector = calc_out_params['xepsilon']
+
+ old_y_component = spectrum_node.get_y()
+ if len(old_y_component) == 1:
+ y_array = old_y_component[0][1]
+ y_units = old_y_component[0][2]
+ y_arrays_list.append(y_array)
+ y_units_list.append(y_units)
+ y_labels_list.append(f'sigma_{eps_vector[0]}_{eps_vector[1]}_{eps_vector[2]}')
+ elif len(old_y_component) == 3:
+ y_tot = old_y_component[0][1]
+ y_tot_units = old_y_component[0][2]
+ y_tot_label = f'sigma_tot_{eps_vector[0]}_{eps_vector[1]}_{eps_vector[2]}'
+ y_arrays_list.append(y_tot)
+ y_units_list.append(y_tot_units)
+ y_labels_list.append(y_tot_label)
+
+ y_up = old_y_component[1][1]
+ y_up_units = old_y_component[1][2]
+ y_up_label = f'sigma_up_{eps_vector[0]}_{eps_vector[1]}_{eps_vector[2]}'
+ y_arrays_list.append(y_up)
+ y_units_list.append(y_up_units)
+ y_labels_list.append(y_up_label)
+
+ y_down = old_y_component[2][1]
+ y_down_units = old_y_component[2][2]
+ y_down_label = f'sigma_down_{eps_vector[0]}_{eps_vector[1]}_{eps_vector[2]}'
+ y_arrays_list.append(y_down)
+ y_units_list.append(y_down_units)
+ y_labels_list.append(y_down_label)
+
+ x_array = spectrum_node.get_x()[1]
+ x_label = spectrum_node.get_x()[0]
+ x_units = spectrum_node.get_x()[2]
+
+ output_spectra.set_x(x_array, x_label, x_units)
+ output_spectra.set_y(y_arrays_list, y_labels_list, y_units_list)
+
+ return output_spectra

src/aiida_quantumespresso/workflows/protocols/core_hole_treatments.yaml

@@ -0,0 +1,29 @@
+default_inputs:
+ SYSTEM:
+ tot_charge: 1
+default_treatment: full
+treatments:
+ full:
+ description: 'Core-hole treatment using a formal countercharge of +1, equivalent to removing one electron from the system.'
+ half:
+ description: 'Core-hole treatment using a formal countercharge of +0.5, equivalent to removing half an electron from the system.'
+ SYSTEM:
+ tot_charge: 0.5
+ xch_fixed:
+ description: 'Core-hole treatment which places the excited electron into the conduction band (fixed occupations).'
+ SYSTEM:
+ occupations: fixed
+ tot_charge: 0
+ nspin: 2
+ tot_magnetization: 1
+ xch_smear:
+ description: 'Core-hole treatment which places the excited electron into the conduction band (smeared occupations).'
+ SYSTEM:
+ occupations: smearing
+ tot_charge: 0
+ nspin: 2
+ starting_magnetization(1): 0
+ none:
+ description: 'Applies no core-hole treatment (overrides the default tot_charge and changes it to 0).'
+ SYSTEM:
+ tot_charge: 0

src/aiida_quantumespresso/workflows/protocols/xspectra/core.yaml

@@ -0,0 +1,17 @@
+default_inputs:
+ abs_atom_marker: X
+ clean_workdir: True
+ get_powder_spectrum: False
+ run_replot: False
+ eps_vectors:
+ - [1., 0., 0.]
+ - [0., 1., 0.]
+ - [0., 0., 1.]
+default_protocol: moderate
+protocols:
+ moderate:
+ description: 'Protocol to perform XANES dipole calculations at normal precision and moderate computational cost.'
+ precise:
+ description: 'Protocol to perform XANES dipole calculations at high precision and higher computational cost.'
+ fast:
+ description: 'Protocol to perform XANES dipole calculations at low precision and minimal computational cost for testing purposes.'