diff --git a/aiida_fleur/calculation/fleur.py b/aiida_fleur/calculation/fleur.py index f4cc50041..e27a265d4 100644 --- a/aiida_fleur/calculation/fleur.py +++ b/aiida_fleur/calculation/fleur.py @@ -104,6 +104,9 @@ class FleurCalculation(CalcJob): #files for greensfunctions _GREENSF_HDF5_FILE_NAME = 'greensf.hdf' + #Name of the folder containing the hubbard1 results + _HUBBARD1_FOLDER_NAME = 'Hubbard1' + # files for hybrid functionals _COULOMB1_FILE_NAME = 'coulomb1' _MIXBAS_FILE_NAME = 'mixbas' @@ -388,11 +391,13 @@ def prepare_for_submission(self, folder): if modes['relax']: # if l_f="T" retrieve relax.xml mode_retrieved_filelist.append(self._RELAX_FILE_NAME) - if modes['ldau']: + if modes['ldau'] or modes['ldahia']: if with_hdf5: mode_retrieved_filelist.append(self._NMMPMAT_HDF5_FILE_NAME) else: mode_retrieved_filelist.append(self._NMMPMAT_FILE_NAME) + if modes['ldahia']: + mode_retrieved_filelist.append(self._HUBBARD1_FOLDER_NAME) if with_hdf5 and modes['greensf']: mode_retrieved_filelist.append(self._GREENSF_HDF5_FILE_NAME) if modes['cf_coeff']: diff --git a/aiida_fleur/workflows/hubbard1.py b/aiida_fleur/workflows/hubbard1.py new file mode 100644 index 000000000..d6542e61e --- /dev/null +++ b/aiida_fleur/workflows/hubbard1.py @@ -0,0 +1,372 @@ +""" +This module contains the FleurHubbard1WorkChain, which is used to perform +a complete hubbard 1 calculation including calculating model parameters beforehand +""" +from collections import defaultdict +import copy + +from aiida.common import AttributeDict +from aiida.common.exceptions import NotExistent +from aiida.engine import WorkChain, ToContext, if_, calcfunction +from aiida import orm + +from aiida_fleur.workflows.scf import FleurScfWorkChain +from .cfcoeff import FleurCFCoeffWorkChain + + +class FleurHubbard1WorkChain(WorkChain): + """ + WorkChain for performing Hubbard-1 calculations + """ + + _workflowversion = '0.0.1' + + _default_wf_para = { + 'ldahia_dict': None, + 'soc': 'auto', + 'exchange_constants': None, + 'exchange_field': None, + 'cf_coefficients': None, + 'occupation_converged': 0.01, + 'matrix_elements_converged': 0.001, + 'itmax_hubbard1': 5, + 'energy_contours': { + 'shape': 'semircircle', + 'eb': -1.0, + 'n': 128 + }, + 'energy_grid': { + 'ellow': -1.0, + 'elup': 1.0, + 'ne': 5400 + }, + 'inpxml_changes': [], + } + + @classmethod + def define(cls, spec): + super().define(spec) + spec.expose_inputs(FleurCFCoeffWorkChain, + namespace='cfcoeff', + namespace_options={ + 'required': False, + 'populate_defaults': False + }) + spec.expose_inputs(FleurScfWorkChain, namespace='hubbard1_scf', namespace_options={ + 'required': True, + }) + spec.input('wf_parameters', valid_type=orm.Dict, required=False) + + #TODO: hubbard1 calculation in loop until converged? + spec.outline( + cls.start, cls.validate_inputs, + if_(cls.preliminary_calcs_needed)(cls.run_preliminary_calculations, cls.inspect_preliminary_calculations), + cls.run_hubbard1_calculation, cls.inspect_hubbard1_calculation, cls.return_results) + + spec.output('output_hubbard1_wc_para', valid_type=orm.Dict) + spec.expose_outputs(FleurScfWorkChain) + + spec.exit_code(230, 'ERROR_INVALID_INPUT_PARAM', message='Invalid workchain parameters.') + spec.exit_code(231, 'ERROR_INVALID_INPUT_CONFIG', message='Invalid input configuration.') + spec.exit_code(377, 'ERROR_CFCOEFF_CALCULATION_FAILED', message='Crystal field coefficient calculation failed') + spec.exit_code(378, 'ERROR_HUBBARD1_CALCULATION_FAILED', message='Hubbard 1 calculation failed') + + def start(self): + """ + init context and some parameters + """ + self.report(f'INFO: started hubbard1 workflow version {self._workflowversion}') + + self.ctx.run_preliminary = False + self.ctx.fleurinp_hubbard1 = None + self.ctx.successful = True + self.ctx.info = [] + self.ctx.errors = [] + self.ctx.warnings = [] + + wf_default = copy.deepcopy(self._default_wf_para) + if 'wf_parameters' in self.inputs: + wf_dict = self.inputs.wf_parameters.get_dict() + else: + wf_dict = wf_default + + for key, val in wf_default.items(): + wf_dict[key] = wf_dict.get(key, val) + self.ctx.wf_dict = wf_dict + + def validate_inputs(self): + """ + Validate the input configuration + """ + + extra_keys = {key for key in self.ctx.wf_dict if key not in self._default_wf_para} + if extra_keys: + error = f'ERROR: input wf_parameters for Orbcontrol contains extra keys: {extra_keys}' + self.report(error) + return self.exit_codes.ERROR_INVALID_INPUT_PARAM + + inputs = self.inputs + if 'cfcoeff' in inputs: + if 'wf_parameters' in inputs: + if inputs.wf_parameters.get_dict().get('convert_to_stevens') is True: + error = 'ERROR: you gave cfcoeff input with explcitly setting convert_to_setvens to True. This needs to be False' + self.report(error) + return self.exit_codes.ERROR_INVALID_INPUT_PARAM + self.ctx.run_preliminary = True + if self.ctx.wf_dict['cf_coefficients'] is not None: + error = 'ERROR: you gave cfcoeff input + explicit cf_coefficients in wf_parameters' + self.report(error) + return self.exit_codes.ERROR_INVALID_INPUT_CONFIG + + if self.ctx.wf_dict['exchange_constants'] is not None and \ + self.ctx.wf_dict['exchange_field'] is not None: + error = 'ERROR: you gave exchange_constants input + exchange_field in wf_parameters' + self.report(error) + return self.exit_codes.ERROR_INVALID_INPUT_CONFIG + + def preliminary_calcs_needed(self): + """ + Return whether calculations need to be run before the hubbard 1 calculation + """ + return self.ctx.run_preliminary + + def run_preliminary_calculations(self): + """ + Run calculations needed before hubbard 1 + """ + + self.report('INFO: Starting Preliminary calculations') + + calcs = {} + if 'cfcoeff' in self.inputs: + self.report('INFO: Starting Crystal field calculation') + + inputs = AttributeDict(self.exposed_inputs(FleurCFCoeffWorkChain, namespace='cfcoeff')) + + if 'wf_parameters' in inputs: + if 'convert_to_stevens' not in inputs: + inputs.wf_parameters = orm.Dict(dict={ + **inputs.wf_parameters.get_dict(), 'convert_to_stevens': False + }) + else: + inputs.wf_parameters = orm.Dict(dict={'convert_to_stevens': False}) + + calcs['cfcoeff'] = self.submit(FleurCFCoeffWorkChain, **inputs) + + return ToContext(**calcs) + + def inspect_preliminary_calculations(self): + """ + Check that the preliminary calculations finished sucessfully + """ + if not self.ctx.cfcoeff.is_finished_ok: + error = ('ERROR: CFCoeff workflow was not successful') + self.report(error) + return self.exit_codes.ERROR_CFCOEFF_CALCULATION_FAILED + + try: + self.ctx.cfcoeff.outputs.output_cfcoeff_wc_para + except NotExistent: + message = ('ERROR: CFCoeff workflow failed, no output node') + self.ctx.errors.append(message) + return self.exit_codes.ERROR_CFCOEFF_CALCULATION_FAILED + + def create_hubbard1_input(self): + """ + Create the inpxml_changes input to create the hubbard 1 input + """ + inputs_hubbard1 = AttributeDict(self.exposed_inputs(FleurScfWorkChain, namespace='hubbard1_scf')) + + if 'wf_parameters' in inputs_hubbard1: + scf_wf_para = inputs_hubbard1.wf_parameters.get_dict() + else: + scf_wf_para = {} + + scf_wf_para.setdefault('inpxml_changes', []).append(('set_inpchanges', { + 'change_dict': { + 'minoccdistance': self.ctx.wf_dict['occupation_converged'], + 'minmatdistance': self.ctx.wf_dict['matrix_elements_converged'], + 'itmaxhubbard1': self.ctx.wf_dict['itmax_hubbard1'] + } + })) + + scf_wf_para['inpxml_changes'].append(('set_simple_tag', { + 'tag_name': 'realAxis', + 'changes': self.ctx.wf_dict['energy_grid'], + 'create_parents': True + })) + + complex_contours = self.ctx.wf_dict['energy_contours'] + if not isinstance(complex_contours, list): + complex_contours = [complex_contours] + + #Build the list of changes (we need to provide them in one go to not override each other) + contour_tags = defaultdict(list) + for contour in complex_contours: + shape = contour.pop('shape') + + if shape.lower() == 'semircircle': + contour_tags['contourSemicircle'].append(contour) + elif shape.lower() == 'dos': + contour_tags['contourDOS'].append(contour) + elif shape.lower() == 'rectangle': + contour_tags['contourRectangle'].append(contour) + + scf_wf_para['inpxml_changes'].append(('set_complex_tag', { + 'tag_name': 'greensFunction', + 'changes': dict(contour_tags) + })) + + if 'cfcoeff' in self.ctx: + #Take the output crystal field cofficients from the CFCoeffWorkchain + #TODO: Select spin channel + cf_coefficients = self.ctx.cfcoeff.outputs.output_cfcoeff_wc_para.dict.cf_coefficients_spin_up + if all('/' not in key for key in cf_coefficients): + self.report('WARNING: Crystal field coefficients from multiple atomtypes not supported fully') + _, cf_coefficients = cf_coefficients.popitem( + ) #Just take one (To support multiple we need to split up species) + + #Drop coefficients with negative m + cf_coefficients = {key: coeff for key, coeff in cf_coefficients.items() if '-' not in key} + else: + cf_coefficients = self.ctx.wf_dict.get('cf_coefficients', {}) + + coefficient_tags = [] + cf_coefficients = cf_coefficients or {} + for key, coeff in cf_coefficients.items(): + l, m = key.split('/') + coefficient_tags.append({'l': l, 'm': m, 'value': coeff}) + + exc_constant_tags = [] + exc_constant = self.ctx.wf_dict.get('exchange_constants', {}) + exc_constant = exc_constant or {} + for l, exc_dict in exc_constant.items(): + exc_constant_tags.append({'l': l, **exc_dict}) + + addarg_tags = [] + if self.ctx.wf_dict['soc'] != 'auto': + addarg_tags.append({'key': 'xiSOC', 'value': self.ctx.wf_dict['soc']}) + + if self.ctx.wf_dict['exchange_field'] is not None: + addarg_tags.append({'key': 'Exc', 'value': self.ctx.wf_dict['exchange_field']}) + + for species_name, ldahia_dict in self.ctx.wf_dict['ldahia_dict'].items(): + if coefficient_tags: + ldahia_dict = {**ldahia_dict, 'cfCoeff': coefficient_tags} + + if exc_constant_tags: + ldahia_dict = {**ldahia_dict, 'exc': exc_constant_tags} + + if addarg_tags: + ldahia_dict.setdefault('addarg', []).extend(addarg_tags) + + scf_wf_para['inpxml_changes'].append(('set_species', { + 'species_name': species_name, + 'attributedict': { + 'ldahia': ldahia_dict + } + })) + scf_wf_para['inpxml_changes'].extend(self.ctx.wf_dict.get('inpxml_changes', [])) + + inputs_hubbard1.wf_parameters = orm.Dict(dict=scf_wf_para) + + return inputs_hubbard1 + + def run_hubbard1_calculation(self): + """ + Start the Hubbard1 calculation by submitting the SCF workchain with the correct inputs + """ + self.report('INFO: Starting Hubbard 1 calculation') + + inputs = self.create_hubbard1_input() + result = self.submit(FleurScfWorkChain, **inputs) + return ToContext(hubbard1=result) + + def inspect_hubbard1_calculation(self): + """ + Inspect the finished Hubbard 1 calculation + """ + if not self.ctx.hubbard1.is_finished_ok: + error = ('ERROR: Hubbard1 SCF workflow was not successful') + self.ctx.successful = False + self.control_end_wc(error) + return self.exit_codes.ERROR_HUBBARD1_CALCULATION_FAILED + + try: + self.ctx.hubbard1.outputs.output_scf_wc_para + except NotExistent: + error = ('ERROR: Hubbard1 SCF workflow failed, no output node') + self.ctx.successful = False + self.control_end_wc(error) + return self.exit_codes.ERROR_HUBBARD1_CALCULATION_FAILED + + def return_results(self): + """ + Return results of the hubbard1 workchain + """ + + try: # if something failed, we still might be able to retrieve something + last_calc_out = self.ctx.hubbard1.outputs.last_calc.output_parameters + retrieved = self.ctx.hubbard1.outputs.last_calc.retrieved + except (NotExistent, AttributeError): + last_calc_out = None + retrieved = None + + outputnode_dict = {} + + outputnode_dict['workflow_name'] = self.__class__.__name__ + outputnode_dict['workflow_version'] = self._workflowversion + outputnode_dict['warnings'] = self.ctx.warnings + outputnode_dict['info'] = self.ctx.info + outputnode_dict['errors'] = self.ctx.errors + outputnode_dict['successful'] = self.ctx.successful + #TODO: What should be in here + # - magnetic moments + # - orbital moments + # What else + + outputnode = orm.Dict(dict=outputnode_dict) + + outdict = create_hubbard1_result_node(outpara=outputnode, + last_calc_out=last_calc_out, + last_calc_retrieved=retrieved) + + if self.ctx.hubbard1: + self.out_many(self.exposed_outputs(self.ctx.hubbard1, FleurScfWorkChain)) + + for link_name, node in outdict.items(): + self.out(link_name, node) + + if not self.ctx.successful: + return self.exit_codes.ERROR_HUBBARD1_CALCULATION_FAILED + + def control_end_wc(self, errormsg): + """ + Controlled way to shutdown the workchain. will initialize the output nodes + The shutdown of the workchain will has to be done afterwards + """ + self.ctx.successful = False + self.ctx.abort = True + self.report(errormsg) # because return_results still fails somewhen + self.ctx.errors.append(errormsg) + self.return_results() + + +@calcfunction +def create_hubbard1_result_node(**kwargs): + """ + This is a pseudo wf, to create the right graph structure of AiiDA. + This wokfunction will create the output node in the database. + It also connects the output_node to all nodes the information commes from. + So far it is just also parsed in as argument, because so far we are to lazy + to put most of the code overworked from return_results in here. + """ + outpara = kwargs['outpara'] #Should always be there + outdict = {} + outputnode = outpara.clone() + outputnode.label = 'output_hubbard1_wc_para' + outputnode.description = ( + 'Contains self-consistency/magnetic information results and information of an FleurHubbard1WorkChain run.') + + outdict['output_hubbard1_wc_para'] = outputnode + return outdict diff --git a/aiida_fleur/workflows/scf.py b/aiida_fleur/workflows/scf.py index c2ba8ba3b..9b3d2d081 100644 --- a/aiida_fleur/workflows/scf.py +++ b/aiida_fleur/workflows/scf.py @@ -71,6 +71,8 @@ class FleurScfWorkChain(WorkChain): 'kpoints_force_even': False, 'kpoints_force_gamma': False, 'nmmp_converged': 0.002, + 'hubbard1_occ_converged': 0.01, + 'hubbard1_elem_converged': 0.001, 'mode': 'density', # 'density', 'energy', 'force' or 'gw' 'add_comp_para': { 'only_even_MPI': False, @@ -216,10 +218,14 @@ def start(self): self.ctx.all_forces = [] self.ctx.total_energy = [] self.ctx.nmmp_distance = [] + self.ctx.hubbard1_occ_distance = [] + self.ctx.hubbard1_elem_distance = [] self.ctx.energydiff = 10000 self.ctx.forcediff = 10000 self.ctx.last_charge_density = 10000 self.ctx.last_nmmp_distance = -10000 + self.ctx.last_hubbard1_occ_distance = -10000 + self.ctx.last_hubbard1_elem_distance = -10000 self.ctx.warnings = [] # "debug": {}, self.ctx.errors = [] @@ -669,6 +675,15 @@ def get_res(self): if nmmp_distances is not None: self.ctx.nmmp_distance.extend(nmmp_distances) + if 'ldahia_info' in output_dict: + occ_distances = output_dict['ldahia_info'].get('occupation_distance', []) + elem_distances = output_dict['ldahia_info'].get('element_distance', []) + + if occ_distances: + #Assume both worked + self.ctx.hubbard1_occ_distance.extend(occ_distances) + self.ctx.hubbard1_elem_distance.extend(elem_distances) + if mode == 'force': forces = output_dict.get('force_atoms', []) if forces is not None: @@ -702,6 +717,10 @@ def get_res(self): 'Assuming that the calculatin should be continued') self.ctx.last_nmmp_distance = self.ctx.wf_dict['nmmp_converged'] + 1 + if self.ctx.hubbard1_occ_distance: + self.ctx.last_hubbard1_occ_distance = self.ctx.hubbard1_occ_distance[-1] + self.ctx.last_hubbard1_elem_distance = self.ctx.hubbard1_elem_distance[-1] + def condition(self): """ check convergence condition @@ -709,6 +728,7 @@ def condition(self): self.report('INFO: checking condition FLEUR') mode = self.ctx.wf_dict['mode'] ldau_notconverged = False + ldahia_notconverged = False energy = self.ctx.total_energy if len(energy) >= 2: @@ -726,13 +746,22 @@ def condition(self): self.ctx.last_nmmp_distance >= self.ctx.wf_dict['nmmp_converged']: ldau_notconverged = True + if self.ctx.last_hubbard1_occ_distance > 0.0 and \ + self.ctx.last_hubbard1_occ_distance >= self.ctx.wf_dict['hubbard1_occ_converged']: + ldahia_notconverged = True + + if self.ctx.last_hubbard1_elem_distance > 0.0 and \ + self.ctx.last_hubbard1_elem_distance >= self.ctx.wf_dict['hubbard1_elem_converged']: + ldahia_notconverged = True + + ignore_convergence_criteria = ldau_notconverged or ldahia_notconverged if mode == 'density': if self.ctx.wf_dict['density_converged'] >= self.ctx.last_charge_density: - if not ldau_notconverged: + if not ignore_convergence_criteria: return False elif mode in ('energy', 'spex'): if self.ctx.wf_dict['energy_converged'] >= self.ctx.energydiff: - if not ldau_notconverged: + if not ignore_convergence_criteria: return False elif mode == 'force': if self.ctx.last_charge_density is None: @@ -741,7 +770,7 @@ def condition(self): except NotExistent: pass else: - if not ldau_notconverged: + if not ignore_convergence_criteria: return False elif self.ctx.wf_dict['density_converged'] >= self.ctx.last_charge_density: @@ -750,7 +779,7 @@ def condition(self): except NotExistent: pass else: - if not ldau_notconverged: + if not ignore_convergence_criteria: return False if self.ctx.loop_count >= self.ctx.max_number_runs: @@ -786,6 +815,14 @@ def return_results(self): if self.ctx.last_nmmp_distance > 0.0: last_nmmp_distance = self.ctx.last_nmmp_distance + last_hubbard1_occ_distance = None + if self.ctx.last_hubbard1_occ_distance > 0.0: + last_hubbard1_occ_distance = self.ctx.last_hubbard1_occ_distance + + last_hubbard1_elem_distance = None + if self.ctx.last_hubbard1_elem_distance > 0.0: + last_hubbard1_elem_distance = self.ctx.last_hubbard1_elem_distance + outputnode_dict = {} outputnode_dict['workflow_name'] = self.__class__.__name__ outputnode_dict['workflow_version'] = self._workflowversion @@ -803,6 +840,10 @@ def return_results(self): outputnode_dict['total_energy_units'] = 'Htr' outputnode_dict['nmmp_distance'] = last_nmmp_distance outputnode_dict['nmmp_distance_all'] = self.ctx.nmmp_distance + outputnode_dict['hubbard1_occ_distance'] = last_hubbard1_occ_distance + outputnode_dict['hubbard1_occ_distance_all'] = self.ctx.hubbard1_occ_distance + outputnode_dict['hubbard1_elem_distance'] = last_hubbard1_elem_distance + outputnode_dict['hubbard1_elem_distance_all'] = self.ctx.hubbard1_elem_distance outputnode_dict['last_calc_uuid'] = last_calc_uuid outputnode_dict['total_wall_time'] = self.ctx.total_wall_time outputnode_dict['total_wall_time_units'] = 's' @@ -857,6 +898,11 @@ def return_results(self): self.report(f'INFO: The LDA+U density matrix is converged to {self.ctx.last_nmmp_distance} change ' 'of all matrix elements') + if self.ctx.last_hubbard1_occ_distance > 0.0: + self.report( + f'INFO: The LDA+Hubbard 1 density matrix is converged to {self.ctx.last_hubbard1_occ_distance} change ' + f'of occupation and {self.ctx.last_hubbard1_elem_distance} change of all matrix elements') + outputnode_t = Dict(dict=outputnode_dict) # this is unsafe so far, because last_calc_out could not exist... if last_calc_out: diff --git a/setup.json b/setup.json index 6eb09da3b..5af9802a1 100644 --- a/setup.json +++ b/setup.json @@ -85,6 +85,7 @@ "fleur.strain = aiida_fleur.workflows.strain:FleurStrainWorkChain", "fleur.eos = aiida_fleur.workflows.eos:FleurEosWorkChain", "fleur.cfcoeff = aiida_fleur.workflows.cfcoeff:FleurCFCoeffWorkChain", + "fleur.hubbard1 = aiida_fleur.workflows.hubbard1:FleurHubbard1WorkChain", "fleur.init_cls = aiida_fleur.workflows.initial_cls:FleurInitialCLSWorkChain", "fleur.corehole = aiida_fleur.workflows.corehole:FleurCoreholeWorkChain", "fleur.mae = aiida_fleur.workflows.mae:FleurMaeWorkChain", diff --git a/tests/test_entrypoints.py b/tests/test_entrypoints.py index 637077e22..2a51a5abb 100644 --- a/tests/test_entrypoints.py +++ b/tests/test_entrypoints.py @@ -178,3 +178,10 @@ def test_fleur_cfcoeff_wc_entry_point(self): workflow = WorkflowFactory('fleur.cfcoeff') assert workflow == FleurCFCoeffWorkChain + + def test_fleur_hubbard1_wc_entry_point(self): + from aiida.plugins import WorkflowFactory + from aiida_fleur.workflows.hubbard1 import FleurHubbard1WorkChain + + workflow = WorkflowFactory('fleur.hubbard1') + assert workflow == FleurHubbard1WorkChain diff --git a/tests/test_workflows_builder_init.py b/tests/test_workflows_builder_init.py index 270e7c402..fcdc2626a 100644 --- a/tests/test_workflows_builder_init.py +++ b/tests/test_workflows_builder_init.py @@ -164,3 +164,11 @@ def test_fleur_cfcoeff_wc_init(self): from aiida_fleur.workflows.cfcoeff import FleurCFCoeffWorkChain builder = FleurCFCoeffWorkChain.get_builder() + + def test_fleur_hubbard1_wc_init(self): + """ + Test the interface of the cfcoeff workchain + """ + from aiida_fleur.workflows.hubbard1 import FleurHubbard1WorkChain + + builder = FleurHubbard1WorkChain.get_builder()