Feature: add radial analysis during cycle

Some structures may have onsite atoms of the same
specie but with different nearest neighbours (e.g. Fe3O4).
The hp.x code can handle this only specifying a maximum
radius, within which one can find the nearest neighbours.
As such radius may change during relaxation of the atomic
structure, an analysis is introduced to define at each iteration
such maximum radius.

src/aiida_quantumespresso_hp/workflows/hubbard.py

@@ -120,6 +120,8 @@ def define(cls, spec):
         spec.input('skip_relax_iterations', valid_type=orm.Int, required=False, validator=validate_positive,
             help=('The number of iterations for skipping the `relax` '
                   'step without performing check on parameters convergence.'))
+        spec.input('radial_analysis', valid_type=orm.Dict, required=False,
+            help='If specified, it performs a nearest neighbour analysis and feed the radius to hp.x')
         spec.input('relax_frequency', valid_type=orm.Int, required=False, validator=validate_positive,
             help='Integer value referring to the number of iterations to wait before performing the `relax` step.')
         spec.expose_inputs(PwRelaxWorkChain, namespace='relax',
@@ -255,6 +257,7 @@ def get_builder_from_protocol(
         builder.hubbard = hubbard
         builder.tolerance_onsite = orm.Float(inputs['tolerance_onsite'])
         builder.tolerance_intersite = orm.Float(inputs['tolerance_intersite'])
+        builder.radial_analysis = orm.Dict(inputs['radial_analysis'])
         builder.max_iterations = orm.Int(inputs['max_iterations'])
         builder.meta_convergence = orm.Bool(inputs['meta_convergence'])
         builder.clean_workdir = orm.Bool(inputs['clean_workdir'])
@@ -550,6 +553,20 @@ def run_hp(self):
         workchain = self.ctx.workchains_scf[-1]
 
         inputs = AttributeDict(self.exposed_inputs(HpWorkChain, namespace='hubbard'))
+
+        if 'radial_analysis' in self.inputs:
+            from qe_tools import CONSTANTS
+
+            kwargs = self.inputs.radial_analysis.get_dict()
+            hubbard_utils = HubbardUtils(self.ctx.current_hubbard_structure)
+            radius = hubbard_utils.get_intersites_radius(**kwargs)  # in Angstrom
+
+            parameters = inputs.hp.parameters.get_dict()
+            parameters['INPUTHP'].pop('num_neigh', None)
+            parameters['INPUTHP']['rmax'] = radius / CONSTANTS.bohr_to_ang
+
+            inputs.hp.parameters = orm.Dict(parameters)
+
         inputs.clean_workdir = self.inputs.clean_workdir
         inputs.hp.parent_scf = workchain.outputs.remote_folder
         inputs.hp.hubbard_structure = self.ctx.current_hubbard_structure
@@ -558,7 +575,7 @@ def run_hp(self):
         running = self.submit(HpWorkChain, **inputs)
 
         self.report(f'launching HpWorkChain<{running.pk}> iteration #{self.ctx.iteration}')
-        return ToContext(workchains_hp=append_(running))
+        self.to_context(**{'workchains_hp': append_(running)})
 
     def inspect_hp(self):
         """Analyze the last completed HpWorkChain.

src/aiida_quantumespresso_hp/workflows/protocols/hubbard.yaml

@@ -4,6 +4,10 @@ default_inputs:
     meta_convergence: True
     tolerance_onsite: 0.1
     tolerance_intersite: 0.01
+    radial_analysis:
+        radius_max: 10.0 # in Angstrom
+        thr: 0.01
+        nn_finder: 'crystal'
     scf:
         kpoints_distance: 0.4
 

tests/conftest.py

@@ -58,7 +58,7 @@ def _fixture_code(entry_point_name):
 
         try:
             return load_code(label=label)
-        except exceptions.NotExistent:
+        except (exceptions.NotExistent, exceptions.MultipleObjectsError):
             return InstalledCode(
                 label=label,
                 computer=fixture_localhost,

tests/workflows/protocols/test_hubbard/test_default.yml

@@ -19,6 +19,10 @@ hubbard:
 hubbard_structure: CoLiO2
 max_iterations: 10
 meta_convergence: true
+radial_analysis:
+  nn_finder: 'crystal'
+  radius_max: 10.0
+  thr: 0.01
 relax:
   base:
     kpoints_distance: 0.15

tests/workflows/test_hubbard.py

@@ -247,6 +247,28 @@ def test_relax_frequency(generate_workchain_hubbard, generate_inputs_hubbard):
     assert not process.should_run_relax()  # skip
 
 
+@pytest.mark.usefixtures('aiida_profile')
+def test_radial_analysis(
+    generate_workchain_hubbard,
+    generate_inputs_hubbard,
+    generate_scf_workchain_node,
+):
+    """Test `SelfConsistentHubbardWorkChain` outline when radial analysis is activated.
+
+    We want to make sure `rmax` is in `hp.parameters`.
+    """
+    inputs = generate_inputs_hubbard()
+    inputs['radial_analysis'] = Dict({})  # no need to specify inputs, it will use the defaults
+    process = generate_workchain_hubbard(inputs=inputs)
+
+    process.setup()
+    process.ctx.workchains_scf = [generate_scf_workchain_node(remote_folder=True)]
+    process.run_hp()
+
+    # parameters = process.ctx['workchains_hp'][-1].inputs['hp']['parameters'].get_dict()
+    # assert 'rmax' in parameters
+
+
 @pytest.mark.usefixtures('aiida_profile')
 def test_should_check_convergence(generate_workchain_hubbard, generate_inputs_hubbard):
     """Test `SelfConsistentHubbardWorkChain.should_check_convergence`."""