seekpath_structure_analysis: HubbardStructureData compatibility

The current implementation of the `seekpath_structure_analysis` calculation function
is not compatible with the `HubbardStructureData` data class. Although the function runs
without errors, the returned primitive structure is now simply a `StructureData`
instance, and hence we lose all information related to the Hubbard parameters.

Here we adapt the `seekpath_structure_analysis` calculation function so that it checks
if the provided input `structure` is a `HubbardStructureData` instance. If so, the
resulting primitive and conventional structures are converted into
`HubbardStructureData` instances with the same _onsite_ Hubbard parameters as the input
structure.

Intersite parameters are much more challenging to preserve correctly when primitivizing
the structure, since they cannot simply be defined in terms of the kinds of the
structure. Rather, they are defined on pairs of site indices and transations, which the
`get_explicit_kpoints_path` does not consider when primitivizing the structure. Hence,
the `seekpath_structure_analysis` function currently does not support intersite
parameters and will raise if the input `structure` has any defined.

Co-authored-by: AndresOrtegaGuerrero <[email protected]>

src/aiida_quantumespresso/calculations/functions/seekpath_structure_analysis.py

@@ -3,6 +3,8 @@
 from aiida.engine import calcfunction
 from aiida.orm import Data
 
+from aiida_quantumespresso.data.hubbard_structure import HubbardStructureData
+
 
 @calcfunction
 def seekpath_structure_analysis(structure, **kwargs):
@@ -28,4 +30,31 @@ def seekpath_structure_analysis(structure, **kwargs):
     # All keyword arugments should be `Data` node instances of base type and so should have the `.value` attribute
     unwrapped_kwargs = {key: node.value for key, node in kwargs.items() if isinstance(node, Data)}
 
-    return get_explicit_kpoints_path(structure, **unwrapped_kwargs)
+    result = get_explicit_kpoints_path(structure, **unwrapped_kwargs)
+
+    if isinstance(structure, HubbardStructureData):
+        result['primitive_structure'] = update_structure_with_hubbard(result['primitive_structure'], structure)
+        result['conv_structure'] = update_structure_with_hubbard(result['conv_structure'], structure)
+
+    return result
+
+
+def update_structure_with_hubbard(structure, orig_structure):
+    """Update the structure based on Hubbard parameters if the input structure is a HubbardStructureData."""
+    from aiida_quantumespresso.utils.hubbard import is_intersite_hubbard
+
+    hubbard_structure = HubbardStructureData.from_structure(structure)
+
+    if is_intersite_hubbard(orig_structure.hubbard):
+        raise NotImplementedError('Intersite Hubbard parameters are not yet supported.')
+
+    for parameter in orig_structure.hubbard.parameters:
+        hubbard_structure.initialize_onsites_hubbard(
+            atom_name=orig_structure.sites[parameter.atom_index].kind_name,
+            atom_manifold=parameter.atom_manifold,
+            value=parameter.value,
+            hubbard_type=parameter.hubbard_type,
+            use_kinds=True,
+        )
+
+    return hubbard_structure

tests/calculations/functions/test_seekpath_analysis.py

@@ -0,0 +1,69 @@
+# -*- coding: utf-8 -*-
+"""Tests for the `seekpath_structure_analysis` function for HubbbardStructureData."""
+import pytest
+
+from aiida_quantumespresso.calculations.functions.seekpath_structure_analysis import seekpath_structure_analysis
+from aiida_quantumespresso.data.hubbard_structure import HubbardStructureData
+
+
+# pylint: disable=W0621
+@pytest.mark.usefixtures('aiida_profile')
+def test_seekpath_analysis(data_regression):
+    """Test the `seekpath_structure_analysis` calculation function for HubbardStructureData."""
+    cell = [[5.43, 0.0, 0.0], [0.0, 5.43, 0.0], [0.0, 0.0, 5.43]]
+    sites = (
+        ('Co', 'Co0', (0.0, 0.0, 0.0)),
+        ('Co', 'Co0', (0.0, 2.715, 2.715)),
+        ('Co', 'Co1', (2.715, 0.0, 2.715)),
+        ('Co', 'Co1', (2.715, 2.715, 0.0)),
+    )
+    orig_structure = HubbardStructureData(cell=cell, sites=sites)
+
+    for hubbard_parameter in [
+        ('Co0', '3d', 4.0, 'Ueff', True),
+        ('Co1', '3d', 3.0, 'U', True),
+    ]:
+        orig_structure.initialize_onsites_hubbard(*hubbard_parameter)
+
+    result = seekpath_structure_analysis(orig_structure)
+
+    prim_structure = result['primitive_structure']
+    conv_structure = result['conv_structure']
+
+    assert isinstance(prim_structure, HubbardStructureData), 'Primitive structure should be a HubbardStructureData'
+    assert isinstance(conv_structure, HubbardStructureData), 'Conventional structure should be a HubbardStructureData'
+
+    assert prim_structure.hubbard.parameters != orig_structure.hubbard.parameters, \
+        'Primitive parameters should be different'
+    assert len(prim_structure.hubbard.parameters) == len(orig_structure.hubbard.parameters), \
+        'Primitive parameters should have the same length as original parameters'
+    assert all(
+        prim_param.atom_manifold == orig_param.atom_manifold
+        for prim_param, orig_param in zip(prim_structure.hubbard.parameters, orig_structure.hubbard.parameters)
+    ), 'Primitive cell parameter atom manifolds should match the original'
+
+    data_regression.check({
+        'primitive': {
+            'cell': prim_structure.cell,
+            'kinds': prim_structure.get_site_kindnames(),
+            'positions': [site.position for site in prim_structure.sites],
+            'hubbard': prim_structure.hubbard.to_list(),
+        },
+        'conventional': {
+            'cell': conv_structure.cell,
+            'kinds': conv_structure.get_site_kindnames(),
+            'positions': [site.position for site in conv_structure.sites],
+            'hubbard': conv_structure.hubbard.to_list(),
+        }
+    })
+
+
+# pylint: disable=W0621
+@pytest.mark.usefixtures('aiida_profile')
+def test_seekpath_analysis_intersite(generate_structure):
+    """Test that the `seekpath_structure_analysis` with intersite hubbard corrections fails."""
+    orig_structure = HubbardStructureData.from_structure(generate_structure('silicon-kinds'))
+    orig_structure.initialize_intersites_hubbard('Si0', '2p', 'Si1', '2p', 4.0, 'V', True)
+
+    with pytest.raises(NotImplementedError, match='Intersite Hubbard parameters are not yet supported.'):
+        seekpath_structure_analysis(orig_structure)

tests/calculations/functions/test_seekpath_analysis/test_seekpath_analysis.yml

@@ -0,0 +1,80 @@
+conventional:
+  cell:
+  - - 3.839589821843
+    - 0.0
+    - 0.0
+  - - 0.0
+    - 3.839589821843
+    - 0.0
+  - - 0.0
+    - 0.0
+    - 5.43
+  hubbard:
+  - - 0
+    - 3d
+    - 0
+    - 3d
+    - 4.0
+    - - 0
+      - 0
+      - 0
+    - Ueff
+  - - 1
+    - 3d
+    - 1
+    - 3d
+    - 3.0
+    - - 0
+      - 0
+      - 0
+    - U
+  kinds:
+  - Co0
+  - Co1
+  positions:
+  - - 0.0
+    - 0.0
+    - 0.0
+  - - 1.9197949109215
+    - 1.9197949109215
+    - 2.715
+primitive:
+  cell:
+  - - 3.839589821843
+    - 0.0
+    - 0.0
+  - - 0.0
+    - 3.839589821843
+    - 0.0
+  - - 0.0
+    - 0.0
+    - 5.43
+  hubbard:
+  - - 0
+    - 3d
+    - 0
+    - 3d
+    - 4.0
+    - - 0
+      - 0
+      - 0
+    - Ueff
+  - - 1
+    - 3d
+    - 1
+    - 3d
+    - 3.0
+    - - 0
+      - 0
+      - 0
+    - U
+  kinds:
+  - Co0
+  - Co1
+  positions:
+  - - 0.0
+    - 0.0
+    - 0.0
+  - - 1.9197949109215
+    - 1.9197949109215
+    - 2.715