diff --git a/.bumpversion.cfg b/.bumpversion.cfg index fe3ab7c1..9351d207 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,5 +1,5 @@ [bumpversion] -current_version = 1.1.12 +current_version = 1.1.13 commit = True tag = True parse = (?P\d+)\.(?P\d+)\.(?P\d+)(\-(?P[a-z]+))? diff --git a/aiida_kkr/__init__.py b/aiida_kkr/__init__.py index 6a9c1669..04ae409b 100644 --- a/aiida_kkr/__init__.py +++ b/aiida_kkr/__init__.py @@ -2,4 +2,4 @@ AiiDA KKR """ -__version__ = '1.1.12' +__version__ = '1.1.13' diff --git a/aiida_kkr/calculations/kkr.py b/aiida_kkr/calculations/kkr.py index 45101dc0..bbcc9f35 100644 --- a/aiida_kkr/calculations/kkr.py +++ b/aiida_kkr/calculations/kkr.py @@ -7,7 +7,7 @@ import os import numpy as np from aiida.engine import CalcJob -from aiida.orm import CalcJobNode, load_node, RemoteData, Dict, StructureData, KpointsData, Bool +from aiida.orm import CalcJobNode, load_node, RemoteData, Dict, StructureData, KpointsData, Bool, FolderData from .voro import VoronoiCalculation from ..tools.common_workfunctions import get_natyp from aiida.common.utils import classproperty @@ -27,7 +27,7 @@ __copyright__ = (u'Copyright (c), 2017, Forschungszentrum Jülich GmbH, ' 'IAS-1/PGI-1, Germany. All rights reserved.') __license__ = 'MIT license, see LICENSE.txt file' -__version__ = '0.12.2' +__version__ = '0.12.5' __contributors__ = ('Jens Bröder', 'Philipp Rüßmann') verbose = False @@ -55,6 +55,7 @@ class KkrCalculation(CalcJob): # List of optional input files (may be mandatory for some settings in inputcard) _SHAPEFUN = 'shapefun' # mandatory if nonspherical calculation _SCOEF = 'scoef' # mandatory for KKRFLEX calculation and some functionalities + _BFIELD = 'bfield.dat' # mandatory if = True _NONCO_ANGLES = 'nonco_angle.dat' # mandatory if noncollinear directions are used that are not (theta, phi)= (0,0) for all atoms _NONCO_ANGLES_IMP = 'nonco_angle_imp.dat' # mandatory for GREENIMP option (scattering code) _SHAPEFUN_IMP = 'shapefun_imp' # mandatory for GREENIMP option (scattering code) @@ -101,6 +102,7 @@ class KkrCalculation(CalcJob): _DECIFILE = 'decifile' # BdG mode _BDG_POT = 'den_lm_ir.%0.3i.%i.txt' + _BDG_CHI_NS = 'den_lm_ns_*.npy' # template.product entry point defined in setup.json _default_parser = 'kkr.kkrparser' @@ -215,6 +217,25 @@ def define(cls, spec): equal to the number of atoms. """ ) + spec.input( + 'bfield', + valid_type=Dict, + required=False, + help="""Non-collinear exteral B-field used for constraint calculations. + + The Dict node should be of the form + initial_noco_angles = Dict(dict={ + 'theta': [theta_at1, theta_at2, ..., theta_atN], + # list theta values in degrees (0..180) + 'phi': [phi_at1, phi_at2, ..., phi_atN], + # list phi values in degrees (0..360) + 'magnitude': [magnitude at_1, ..., magnitude at_N] + # list of magnitude of the applied fields in Ry units + }) + Note: The length of the theta, phi and magnitude lists have to be + equal to the number of atoms. + """ + ) spec.input( 'deciout_parent', valid_type=RemoteData, @@ -230,6 +251,15 @@ def define(cls, spec): the kkrflex_* files are copied to the retrieved (can clutter the database) or are ony left in the remote folder.""" ) + spec.input( + 'anomalous_density', + valid_type=FolderData, + required=False, + help="""FolderData that contains anomalous density input files for + the KKRhost BdG calculation. If these are not give the code looks + for them in the retrieved of the parent calculation and takes them + from there.""" + ) # define outputs spec.output( @@ -295,6 +325,10 @@ def prepare_for_submission(self, tempfolder): if 'initial_noco_angles' in self.inputs: parameters = self._use_initial_noco_angles(parameters, structure, tempfolder) + # write bfield.dat file and add '= True' to input parameters + if 'bfield' in self.inputs: + parameters = self._use_nonco_bfield(parameters, structure, tempfolder) + # activate decimation mode and copy decifile from deciout parent if 'deciout_parent' in self.inputs: parameters = self._use_decimation(parameters, tempfolder) @@ -802,10 +836,13 @@ def _get_BdG_filelist(self, parameters, natom, nspin): self.report(f'retrieve BdG? {retrieve_BdG_files}') if retrieve_BdG_files: + # anomalous density files for all atoms if present for iatom in range(natom): for ispin in range(nspin): print('adding files for BdG mode') add_files += [self._BDG_POT % (iatom + 1, ispin + 1)] + # radially-averaged anomalous density matrix (for triplet components etc.) + add_files.append(self._BDG_CHI_NS) #also retrieve BdG DOS files for anomalous density and hole part for BdGadd in ['_eh', '_he', '_hole']: @@ -1015,6 +1052,51 @@ def _use_initial_noco_angles(self, parameters, structure, tempfolder): return parameters + def _use_nonco_bfield(self, parameters, structure, tempfolder): + """ + Set external non-collinear bfield (writes bfield.dat to tempfolder) used in constraint calculations. + """ + self.report('Found `bfield` input node, writing nonco_angle.dat file') + + # extract number of atoms for length comparison + natom = get_natyp(structure) + + change_values = [['', True]] + parameters = _update_params(parameters, change_values) + + # extract magnitude, theta and phi values from input node + mags = self.inputs.bfield['magnitude'] + if len(mags) != natom: + raise InputValidationError( + 'Error: `magnitude` list in `bfield` input node needs to have the same length as number of atoms!' + ) + thetas = self.inputs.bfield['theta'] + if len(thetas) != natom: + raise InputValidationError( + 'Error: `theta` list in `bfield` input node needs to have the same length as number of atoms!' + ) + phis = self.inputs.bfield['phi'] + if len(phis) != natom: + raise InputValidationError( + 'Error: `phi` list in `bfield` input node needs to have the same length as number of atoms!' + ) + + # now write kkrflex_angle file + with tempfolder.open(self._BFIELD, 'w') as bfield_file: + bfield_file.write('# theta [deg] phi [deg] magnitude [Ry]\n') + for iatom in range(natom): + theta, phi = thetas[iatom], phis[iatom] + magnitude = mags[iatom] + # check consistency + if theta < 0. or theta > 180.: + raise InputValidationError( + f'Error: theta value out of range (0..180): iatom={iatom}, theta={theta}' + ) + # write line + bfield_file.write(f' {theta} {phi} {magnitude}\n') + + return parameters + def _use_decimation(self, parameters, tempfolder): """ Activate decimation mode and copy decifile from output of deciout_parent calculation @@ -1063,13 +1145,26 @@ def _copy_BdG_pot(self, retrieved, tempfolder): Activate BdG mode and copy den_lm_ir files of the previous output to the input of this calculation. """ - BDG_POT_FILES = [i for i in retrieved.list_object_names() if self._BDG_POT.split('.')[0] in i] + if 'anomalous_density' in self.inputs: + # this means we have a FolderData input that contains the + # anomalous density files + adens_folder = self.inputs.anomalous_density + else: + # if no anomalous density is given as an input node we check + # if there are any anomalous density files in the parent retrieved + # and take them from there if present + adens_folder = retrieved + + # list of den_lm_ir files (anomalous density per atom) + BDG_POT_FILES = [i for i in adens_folder.list_object_names() if self._BDG_POT.split('.')[0] in i] # add 'den-lm_ir' files to input for BdG_pot in BDG_POT_FILES: - self.report(f'Copy BdG potential {BdG_pot}') - with retrieved.open(BdG_pot, 'r') as file_handle: + self.report(f'Copy BdG potential {BdG_pot} from {adens_folder.uuid}') + # read from parent + with adens_folder.open(BdG_pot, 'r') as file_handle: file_txt = file_handle.readlines() + # write to tempfolder with tempfolder.open(BdG_pot, 'w') as file_handle: file_handle.writelines(file_txt) diff --git a/aiida_kkr/tools/__init__.py b/aiida_kkr/tools/__init__.py index a906ffa6..d80a0d0f 100644 --- a/aiida_kkr/tools/__init__.py +++ b/aiida_kkr/tools/__init__.py @@ -18,6 +18,7 @@ from .kick_out_core_states import * from .neworder_potential import * from .find_parent import get_calc_from_remote, get_remote, get_parent +from .bdg_tools import get_anomalous_density_data # expose structure finder from VoronoiCalculation diff --git a/aiida_kkr/tools/bdg_tools.py b/aiida_kkr/tools/bdg_tools.py new file mode 100644 index 00000000..9703d318 --- /dev/null +++ b/aiida_kkr/tools/bdg_tools.py @@ -0,0 +1,48 @@ +""" +Helper tools that deal with the anomalous density of the BdG formalism in KKR +""" + +from aiida.engine import calcfunction +from aiida.orm import FolderData + + +@calcfunction +def get_anomalous_density_data(retrieved, rename_files=None): + """ + Extract anomalous density files from a retrieved folder of a KkrCalculation + and copy into a new FolderData. This FolderData is then returned and can be + used as the anomalous_density FolderData input to a KkrCalculation. + + :param retrieved: retrieved FolderData of a parent KKRhost BdG calculation where anomalous density files are stored (called `den_lm_ir.AAA.1.txt` where AAA is an atom index) + + :param rename_files: Optional Dict node where mappings of file names are defined. This is helpful if the atom index in the new structure is different from the original calculation (e.g. atom N of the original structure corresponds to atom M of the new structure, then the renaming dict should be {N: M}). Indices that are not found or that map to None are skipped and will not appear in the returned FolderData! + + :returns: anomalous_density FolderData which contains the (possibly renamed) anomalous density files + """ + BdG_files = [i for i in retrieved.list_object_names() if 'den_lm_ir' in i] + + anomalous_density = FolderData() + for fname in BdG_files: + with retrieved.open(fname, 'r') as _fin: + # default is to use the same name as in the input + fname_out = fname + + # rename file, if natom is found in the renaming dict + if rename_files is not None: + # find atom index (should be integer) + natom = int(fname.split('.')[1]) + # we pick the new atom index and recreate the filename + natom_new = {int(k): v for k, v in rename_files.get_dict().items()}.get(natom) + if natom_new is None: + # if no renaming is mapping is found + # we do not copy this anomalous density but leave it out + fname_out = None + else: + # construct changed name + fname_out = 'den_lm_ir.%0.3i.1.txt' % natom_new + + if fname_out is not None: + # copy input file to FolderData node, maybe with changed name + anomalous_density.put_object_from_filelike(_fin, fname_out) + + return anomalous_density diff --git a/aiida_kkr/tools/plot_kkr.py b/aiida_kkr/tools/plot_kkr.py index e2aca2e8..a368e508 100644 --- a/aiida_kkr/tools/plot_kkr.py +++ b/aiida_kkr/tools/plot_kkr.py @@ -16,7 +16,7 @@ __copyright__ = (u'Copyright (c), 2018, Forschungszentrum Jülich GmbH, ' 'IAS-1/PGI-1, Germany. All rights reserved.') __license__ = 'MIT license, see LICENSE.txt file' -__version__ = '0.7.0' +__version__ = '0.7.1' __contributors__ = ('Philipp Rüßmann') @@ -1896,7 +1896,7 @@ def get_ef_from_parent(node): except: with node.outputs.retrieved.open('output.0.txt', mode='r') as file_handle: txt = file_handle.readlines() - iline = search_string('Fermi energy', txt) + iline = search_string('Fermi energy =', txt) if iline >= 0: ef = txt[iline].split('=')[1] ef = float(ef.split()[0]) diff --git a/aiida_kkr/workflows/_decimation.py b/aiida_kkr/workflows/_decimation.py index a45c21ec..1e57ecd1 100644 --- a/aiida_kkr/workflows/_decimation.py +++ b/aiida_kkr/workflows/_decimation.py @@ -10,14 +10,16 @@ from aiida.engine import WorkChain, ToContext, calcfunction from aiida.orm import Code, Dict, Int, Float, RemoteData, KpointsData, XyData, StructureData, FolderData from aiida_kkr.tools.common_workfunctions import test_and_get_codenode -from aiida_kkr.tools import kkrparams +from aiida_kkr.tools import kkrparams, get_anomalous_density_data +from aiida_kkr.workflows.bs import set_energy_params + import numpy as np from masci_tools.io.common_functions import get_Ry2eV __copyright__ = (u'Copyright (c), 2020, Forschungszentrum Jülich GmbH, ' 'IAS-1/PGI-1, Germany. All rights reserved.') __license__ = 'MIT license, see LICENSE.txt file' -__version__ = '0.1.0' +__version__ = '0.2.1' __contributors__ = u'Philipp Rüßmann' _eV2Ry = 1.0 / get_Ry2eV() @@ -48,15 +50,22 @@ class kkr_decimation_wc(WorkChain): 'kmesh': [50, 50, 50]}, # k-mesh used in dos calculation } - :param wf_parameters: (Dict); Workchain specifications - :param options: (Dict); specifications for the computer (used in decimation step only) - :param remote_data: (RemoteData), mandatory; either parent slab or previous decimation calculation - :param kkr: (Code), mandatory; KKR code for running deci-out and decimation steps - :param voronoi: (Code), mandatory if starting from slab calculation; voronoi code for auxiliary calculations - - :return results: (Dict), Information of workflow results - like Success, last result node, list with convergence behavior - :return deci_calc: (RemoteData), Remote data of decimation calculation, used to reuse decimation setup if calculation is continued (e.g. for DOS after scf) + :param wf_parameters: Dict node with workchain parameters (see kkr_decimation_wc.get_wf_defaults()) + :param options: Dict node with specifications for the computer (used in decimation step only) + :param remote_data: mandatory RemoteData node of either a parent slab or previous decimation calculation + :param kkr: mandatory Code node with KKR code for running deci-out and decimation steps + :param voronoi: Code node that is mandatory if starting from slab calculation. Is the voronoi code for auxiliary calculations + :param kpoints: KpointsData node that triggers a band structure calculation. The kpoints specify the k-point path along which the bandstructure is computed with the qdos mode of KKRhost. + :param calc_parameters: Dict node that contains KKR parameters which overwrites settings from the slab parent. + + :returns structure_decimate: StructureData node of the structure of the decimation region. + :returns structure_substrate: StructureData node of the structure of thesubstrate lattice continuation. + :returns out_params_calc_deci_out: Dict node of the output parameters of the deci-out calculation. + :returns out_params_calc_decimate: Dict node of the output parameters of the decimation calculation. + :returns out_remote_calc_decimate: RemoteData node of the decimation calculation. + :returns out_retrieved_calc_decimate: retrieved FolderData node of the decimation calculation. + :returns dos_data: XyData node with the DOS data at finite imaginary part in the energy contour. Only present in DOS mode. + :returns dos_data_interpol: XyData node with the interpolated DOS data onto the real axis. Only present in DOS mode. """ _workflowversion = __version__ @@ -215,6 +224,14 @@ def define(cls, spec): spec.exit_code( 302, 'ERROR_VORONOICODE_NOT_CORRECT', 'The code you provided for voronoi does not use the plugin kkr.voro' ) + spec.exit_code( + 303, 'ERROR_VORO_SUBST_FAILED', 'The voronoi step for the starting potential of the substrate failed.' + ) + spec.exit_code( + 304, 'ERROR_VORO_DECI_FAILED', + 'The voronoi step for the starting potential of the decimation region failed.' + ) + spec.exit_code(305, 'ERROR_DECIOUT_FAILED', 'The deci-out step (writeout of continuation GF) failed.') ################################################################################################################### # functions from outline @@ -225,7 +242,7 @@ def start(self): """ self.report(f'INFO: started KKR decimation workflow version {self._workflowversion}') - ####### init ####### + ####### init ####### # input para wf_dict = self.inputs.wf_parameters.get_dict() @@ -309,10 +326,12 @@ def prepare_deci_from_slab(self): deci_calc = deci_remote.get_incoming(node_class=KkrCalculation).first().node struc_substrate, voro_substrate = VoronoiCalculation.find_parent_structure(deci_calc.inputs.deciout_parent) - self.ctx.struc_substrate, self.ctx.voroaux_substrate = struc_substrate, voro_substrate + self.ctx.struc_substrate = struc_substrate + self.ctx.voroaux_substrate = voro_substrate struc_deci, voro_deci = VoronoiCalculation.find_parent_structure(deci_remote) - self.ctx.struc_decimation, self.ctx.voroaux_decimation = struc_deci, voro_deci + self.ctx.struc_decimation = struc_deci + self.ctx.voroaux_decimation = voro_deci self.ctx.startpot_substrate = voro_substrate.inputs.potential_overwrite self.ctx.startpot_decimate = voro_deci.inputs.potential_overwrite @@ -327,36 +346,30 @@ def prepare_deci_from_slab(self): self.ctx.params_overwrite = self.inputs.calc_parameters # take care of qdos or dos modes - qdos_mode = 0 + qdos_mode = False if 'kpoints' in self.inputs: - qdos_mode = 1 - if qdos_mode == 1 or self.ctx.dosmode: + qdos_mode = True + if qdos_mode or self.ctx.dosmode: if self.ctx.params_overwrite is None: params_new = {} else: params_new = self.ctx.params_overwrite.get_dict() # set dos contour from dos_params - params_new['NPOL'] = 0 - params_new['NPT1'] = 0 - params_new['NPT2'] = self.ctx.dos_params['nepts'] - params_new['IEMXD'] = self.ctx.dos_params['nepts'] - params_new['NPT3'] = 0 - params_new['TEMPR'] = self.ctx.dos_params['tempr'] ef_ry = self.ctx.slab_calc.outputs.output_parameters['fermi_energy'] - params_new['EMIN'] = ef_ry + self.ctx.dos_params['emin_EF'] * _eV2Ry - params_new['EMAX'] = ef_ry + self.ctx.dos_params['emax_EF'] * _eV2Ry - params_new['BZDIVIDE'] = self.ctx.dos_params['kmesh'] - if qdos_mode == 1: - if params_new['TEMPR'] > 100.: - params_new['TEMPR'] = 50. - params_new['BZDIVIDE'] = [1, 1, 1] + _rename = {'emin_EF': 'emin', 'emax_EF': 'emax'} + econt_new = {_rename.get(k, k.lower()): v for k, v in self.ctx.dos_params.items()} + if qdos_mode: + if econt_new['tempr'] > 100.: + econt_new['tempr'] = 50. + econt_new['kmesh'] = [1, 1, 1] + params_new = set_energy_params(econt_new, ef_ry, kkrparams(**params_new)) self.ctx.params_overwrite = Dict(dict=params_new) alat_slab = self.ctx.slab_calc.outputs.output_parameters['alat_internal'] out = make_decimation_param_nodes( - self.ctx.slab_calc.inputs.parameters, Float(alat_slab), self.ctx.struc_decimation, Int(self.ctx.nkz), - self.ctx.params_overwrite + self.ctx.slab_calc.inputs.parameters, Float(alat_slab), self.ctx.struc_decimation, self.ctx.struc_substrate, + Int(self.ctx.nkz), self.ctx.params_overwrite ) self.ctx.dsubstrate = out['dsubstrate'] @@ -376,56 +389,67 @@ def run_voroaux(self): run auxiliary voronoi steps if needed """ # only needed if parent is not already a decimation calculation - if not self.ctx.parent_is_deci: - # set up voronoi calculation for substrate - builder = VoronoiCalculation.get_builder() - builder.code = self.inputs.voronoi - builder.parameters = self.ctx.dsubstrate - builder.structure = self.ctx.struc_substrate - builder.metadata.label = 'auxiliary_voronoi_substrate' # pylint: disable=no-member - builder.metadata.options = self.ctx.options # pylint: disable=no-member - builder.metadata.options['resources'] = {'tot_num_mpiprocs': 1, 'num_machines': 1} # pylint: disable=no-member - builder.potential_overwrite = self.ctx.startpot_substrate - - # submit voroaux for substrate calculation - future_substrate = self.submit(builder) - self.report(f'INFO: running voroaux for substrate (pk: {future_substrate.pk})') - - # set up voronoi calculation for decimation - builder = VoronoiCalculation.get_builder() - builder.code = self.inputs.voronoi - builder.parameters = self.ctx.ddecimation - builder.structure = self.ctx.struc_decimation - builder.metadata.label = 'auxiliary_voronoi_decimation' # pylint: disable=no-member - builder.metadata.options = self.ctx.options # pylint: disable=no-member - builder.metadata.options['resources'] = {'tot_num_mpiprocs': 1, 'num_machines': 1} # pylint: disable=no-member - builder.potential_overwrite = self.ctx.startpot_decimation - - # submit voroaux for substrate calculation - future_decimation = self.submit(builder) - self.report(f'INFO: running voroaux for decimation region (pk: {future_decimation.pk})') - - return ToContext(voroaux_substrate=future_substrate, voroaux_decimation=future_decimation) - - else: + if self.ctx.parent_is_deci: self.report('Skip voroaux steps due to previous decimation input') + return + + # set up voronoi calculation for substrate + builder = VoronoiCalculation.get_builder() + builder.code = self.inputs.voronoi + builder.parameters = self.ctx.dsubstrate + builder.structure = self.ctx.struc_substrate + builder.metadata.label = 'auxiliary_voronoi_substrate' # pylint: disable=no-member + builder.metadata.options = self.ctx.options # pylint: disable=no-member + builder.metadata.options['resources'] = {'tot_num_mpiprocs': 1, 'num_machines': 1} # pylint: disable=no-member + builder.potential_overwrite = self.ctx.startpot_substrate + + # submit voroaux for substrate calculation + future_substrate = self.submit(builder) + self.report(f'INFO: running voroaux for substrate (pk: {future_substrate.pk})') + + # set up voronoi calculation for decimation + builder = VoronoiCalculation.get_builder() + builder.code = self.inputs.voronoi + builder.parameters = self.ctx.ddecimation + builder.structure = self.ctx.struc_decimation + builder.metadata.label = 'auxiliary_voronoi_decimation' # pylint: disable=no-member + builder.metadata.options = self.ctx.options # pylint: disable=no-member + builder.metadata.options['resources'] = {'tot_num_mpiprocs': 1, 'num_machines': 1} # pylint: disable=no-member + builder.potential_overwrite = self.ctx.startpot_decimation + + # submit voroaux for substrate calculation + future_decimation = self.submit(builder) + self.report(f'INFO: running voroaux for decimation region (pk: {future_decimation.pk})') + + return ToContext(voroaux_substrate=future_substrate, voroaux_decimation=future_decimation) def run_deciout(self): """ run KKR calculation for deci-out step """ + if not self.ctx.voroaux_substrate.is_finished_ok: + return self.exit_codes.ERROR_VORO_SUBST_FAILED # pylint: disable=no-member + builder = KkrCalculation.get_builder() builder.code = self.inputs.kkr builder.parameters = self.ctx.dsubstrate builder.metadata.options = self.ctx.options # pylint: disable=no-member - # force serial run: + + # force serial run, otherwise KKRhost code does not work: builder.metadata.options['resources'] = {'tot_num_mpiprocs': 1, 'num_machines': 1} # pylint: disable=no-member builder.metadata.label = 'deci-out' # pylint: disable=no-member builder.parent_folder = self.ctx.voroaux_substrate.outputs.remote_folder + # create and set initial nonco_angles if needed if 'initial_noco_angles' in self.ctx.slab_calc.inputs: builder.initial_noco_angles = self.ctx.noco_angles_substrate + # for BdG mode we have to set the correct anomalous density + is_BdG, adens = self._get_adens_substrate() + if is_BdG: + builder.anomalous_density = adens + + # finally submit calculation future = self.submit(builder) self.report(f'INFO: running deci-out step (pk: {future.pk})') @@ -435,6 +459,11 @@ def run_decimation(self): """ run KKR calculation for decimation step """ + if not self.ctx.deciout_calc.is_finished_ok: + return self.exit_codes.ERROR_DECIOUT_FAILED # pylint: disable=no-member + if not self.ctx.voroaux_decimation.is_finished_ok: + return self.exit_codes.ERROR_VORO_DECI_FAILED # pylint: disable=no-member + builder = KkrCalculation.get_builder() builder.code = self.inputs.kkr builder.parameters = self.ctx.ddecimation @@ -445,12 +474,19 @@ def run_decimation(self): if 'kpoints' in self.inputs: builder.kpoints = self.inputs.kpoints self.report('INFO: detected kpoints input: run qdos calculation') + # create and set initial nonco_angles if needed if 'initial_noco_angles' in self.ctx.slab_calc.inputs: builder.initial_noco_angles = self.ctx.noco_angles_decimation + # for BdG mode we have to set the correct anomalous density + is_BdG, adens = self._get_adens_decimate() + if is_BdG: + builder.anomalous_density = adens + + # finally submit calculation future = self.submit(builder) - self.report(f'INFO: running deci-out step (pk: {future.pk})') + self.report(f'INFO: running decimation step (pk: {future.pk})') return ToContext(decimation_calc=future) @@ -521,6 +557,64 @@ def _create_startpots(self): self.ctx.startpot_substrate = startpot_substrate self.ctx.startpot_decimation = startpot_deci + def _get_adens_substrate(self): + """ + Extract the correct anomalous density which is used for the substrate + """ + adens = None + + # first check if calculation is in BdG mode + dd = self.ctx.dsubstrate.get_dict() + is_BdG = dd.get('USE_BDG', dd.get('', False)) + self.report(f'is BdG? {is_BdG}') # debug output + + if is_BdG: + # now get the atom indices from the slab parent and relabel + # then starting with atom index 1 (remember that Fortran starts counting at 1 and not 0) + retrieved = self.ctx.slab_calc.outputs.retrieved + params = self.ctx.slab_calc.inputs.parameters.get_dict() + nrbasis = params.get('', params.get('NRBASIS')) + nplayer = self.ctx.nplayer + nprinc = self.ctx.nprinc + rename_files = Dict( + dict=dict( + # next nrbasis atoms after nplayer*nprinc are the substrate atoms + # format: (index in slab, index in substrate) + # Note: AiiDA needs the key in the dict to be a string instead of an integer + [(str(nplayer * nprinc + i + 1), i + 1) for i in range(nrbasis)] + ) + ) + # copy and relabel the anomalous density files + adens = get_anomalous_density_data(retrieved, rename_files) + + return is_BdG, adens + + def _get_adens_decimate(self): + """ + Extract the correct anomalous density which is used for the decimation region + """ + adens = None + + # first check if calculation is in BdG mode + dd = self.ctx.ddecimation.get_dict() + is_BdG = dd.get('USE_BDG', dd.get('', False)) + + if is_BdG: + # now copy the anomalous density files to a FolderData that can be the input to the KkrCalculation + retrieved = self.ctx.slab_calc.outputs.retrieved + nplayer = self.ctx.nplayer + nprinc = self.ctx.nprinc + rename_files = Dict( + dict=dict( + # the first nplayer*nprinc are the decimation region + [(str(i + 1), i + 1) for i in range(nplayer * nprinc)] + ) + ) + # copy only the files in the renaming list, others are ignored and not copied + adens = get_anomalous_density_data(retrieved, rename_files) + + return is_BdG, adens + ################################################################################################################### @@ -587,21 +681,17 @@ def get_substrate_structure(nprinc, nplayer, slab_calc): return struc_substrate -@calcfunction -def make_decimation_param_nodes(slab_calc_params, slab_alat, struc_deci, nkz, params_overwrite=None): +def _make_d_substrate(d, params_overwrite, slab_alat, nkz): """ - Create parameter nodes for deci-out and + Create the input parameters for the substrate calculation """ - # prepare decimation params - d = {k: v for k, v in slab_calc_params.get_dict().items() if v is not None} - - # make kkr params for substrate calculation (i.e. deci-out mode, needs to run in serial!) dsubstrate = d.copy() for k in kkr_decimation_wc._keys2d: if k in dsubstrate: dsubstrate.pop(k) dsubstrate = kkrparams(**dsubstrate) dsubstrate.set_value('NSTEPS', 1) + # add deci-out option runopts = dsubstrate.get_value('RUNOPT') if runopts is None: @@ -609,29 +699,38 @@ def make_decimation_param_nodes(slab_calc_params, slab_alat, struc_deci, nkz, pa runopts = [i for i in runopts if i != 'DECIMATE'] # remove decimate flag runopts += ['deci-out'] dsubstrate.set_value('RUNOPT', runopts) + # increase BZDIVIDE[2] if needed bzdiv = dsubstrate.get_value('BZDIVIDE') if bzdiv is not None: bzdiv[2] = nkz.value dsubstrate.set_value('BZDIVIDE', bzdiv) + # overwrite params from input node if params_overwrite is not None: for k, v in params_overwrite.get_dict().items(): dsubstrate.set_value(k, v) + # use alat from slab calculation dsubstrate = { k: v for k, v in dsubstrate.get_dict().items() if v is not None } # clean up removing None values from dict dsubstrate['ALATBASIS'] = slab_alat.value dsubstrate['use_input_alat'] = True - # now create Dict node with params - dsubstrate = Dict(dict=dsubstrate) + return dsubstrate + + +def _make_d_deci(d, struc_deci, params_overwrite, slab_alat): + """ + Create the input parameters for the substrate calculation + """ # set kkr params for substrate writeout ddeci = kkrparams(**d) ddeci.set_value('NSTEPS', 1) ddeci.set_value('DECIFILES', ['vacuum', 'decifile']) ddeci.set_multiple_values(**struc_deci.extras['kkr_settings']) + # add decimate runopt runopts = ddeci.get_value('RUNOPT') if runopts is None: @@ -642,11 +741,55 @@ def make_decimation_param_nodes(slab_calc_params, slab_alat, struc_deci, nkz, pa if params_overwrite is not None: for k, v in params_overwrite.get_dict().items(): ddeci.set_value(k, v) + # use alat from slab calculation ddeci = {k: v for k, v in ddeci.get_dict().items() if v is not None} # clean up removing None values from dict ddeci['ALATBASIS'] = slab_alat.value ddeci['use_input_alat'] = True - # now create Dict node with params + + return ddeci + + +def _adapt_array_sizes(params_dict, pick_layers): + """ + Check the params dict for array entries which should be changed to the correct size. + + This is needed because the decimation region is smaller than the parent slab structure + and the substrate only uses a part from the middle of the arrays. + """ + # this list of keywords was copied from kkrparams of masci-tools + for key in [ + '', '', '', '', 'XINIPOL', '', '', '', + '' + ]: + val = params_dict.get(key) + if val is not None: + # update size by picking layers + val = np.array(val)[pick_layers] + params_dict[key] = val + + +@calcfunction +def make_decimation_param_nodes(slab_calc_params, slab_alat, struc_deci, struc_substrate, nkz, params_overwrite=None): + """ + Create parameter nodes for deci-out and decimation steps + """ + # prepare decimation params + d = {k: v for k, v in slab_calc_params.get_dict().items() if v is not None} + + # make kkr params for substrate calculation (i.e. deci-out mode, needs to run in serial!) + + dsubstrate = _make_d_substrate(d, params_overwrite, slab_alat, nkz) + ddeci = _make_d_deci(d, struc_deci, params_overwrite, slab_alat) + + # modify array inputs to the right size + Ndeci = len(struc_deci.sites) + _adapt_array_sizes(ddeci, pick_layers=[i for i in range(Ndeci)]) + Nsubstrate = len(struc_substrate.sites) + _adapt_array_sizes(dsubstrate, pick_layers=[Ndeci + i for i in range(Nsubstrate)]) + + # now create Dict nodes with params + dsubstrate = Dict(dict=dsubstrate) ddeci = Dict(dict=ddeci) return {'dsubstrate': dsubstrate, 'ddeci': ddeci} diff --git a/aiida_kkr/workflows/bs.py b/aiida_kkr/workflows/bs.py index ee525bde..44b9adf5 100644 --- a/aiida_kkr/workflows/bs.py +++ b/aiida_kkr/workflows/bs.py @@ -23,7 +23,7 @@ __copyright__ = (u'Copyright (c), 2020, Forschungszentrum Jülich GmbH, ' 'IAS-1/PGI-1, Germany. All rights reserved.') __license__ = 'MIT license, see LICENSE.txt file' -__version__ = '0.1.2' +__version__ = '0.1.5' __contributors__ = (u'Rubel Mozumder', u'Philipp Rüßmann') @@ -61,6 +61,7 @@ class kkr_bs_wc(WorkChain): 'nepts': 96, # number of energy points 'RCLUSTZ': None, # can be used to increase the cluster radius if a value is set here 'tempr': 50., # smearing temperature in K + 'kmesh': None, # k-point integration mesh, only useful for CPA calculation } _options_default = { @@ -69,7 +70,10 @@ class kkr_bs_wc(WorkChain): 'num_machines': 1 }, 'withmpi': True, - 'queue_name': '' + 'queue_name': '', + 'prepend_text': '', + 'append_text': '', + 'additional_retrieve_list': None } @classmethod @@ -125,6 +129,13 @@ def define(cls, spec): help="""Initial non-collinear angles for the magnetic moments. See KkrCalculation for details. If this is found in the input potentially extracted nonco angles from the parent calulation are overwritten!""" ) + # maybe overwrite some settings from the KKRhost convergence run + spec.input( + 'params_kkr_overwrite', + valid_type=Dict, + required=False, + help='Overwrite some input parameters of the parent KKR calculation.' + ) # Here outputs are defined spec.output('results_wf', valid_type=Dict, required=True) @@ -181,6 +192,11 @@ def start(self): if options_dict == {}: self.report('INFO: Using default wf Options') options_dict = self._options_default + self.ctx.append_text = options_dict.get('append_text', self._options_default['append_text']) + self.ctx.prepend_text = options_dict.get('prepend_text', self._options_default['prepend_text']) + self.ctx.additional_retrieve_list = options_dict.get( + 'additional_retrieve_list', self._options_default['additional_retrieve_list'] + ) self.ctx.withmpi = options_dict.get('withmpi', self._options_default['withmpi']) self.ctx.resources = options_dict.get('resources', self._options_default['resources']) self.ctx.max_wallclock_seconds = options_dict.get( @@ -243,19 +259,31 @@ def validate_input(self): output_remote = last_calc.outputs.remote_folder self.inputs.remote_data = output_remote + + # extract structure + struc_kkr, _ = VoronoiCalculation.find_parent_structure(self.inputs.remote_data) + # save if structure is an alloy + self.ctx.struc_is_alloy = struc_kkr.is_alloy + # To validate for kpoints if 'kpoints' in inputs: self.ctx.BS_kpoints = inputs.kpoints input_ok = True self.ctx.structure_data = 'None (kpoints taken from input)' else: - struc_kkr, remote_voro = VoronoiCalculation.find_parent_structure(self.inputs.remote_data) #create an auxiliary structure with unique kind_names, this leads to using the input structure in the seekpath method instead of finding the primitive one - saux = StructureData(cell=struc_kkr.cell) + cell = np.array(struc_kkr.cell) + if not struc_kkr.pbc[2]: + # 2D structure, make sure the third bravais vector points along z + cell[2] = np.cross(cell[0], cell[1]) + saux = StructureData(cell=cell) for isite, site in enumerate(struc_kkr.sites): kind = struc_kkr.get_kind(site.kind_name) saux.append_atom( - name='atom' + str(isite) + ':' + site.kind_name, symbols=kind.symbol, position=site.position + name='atom' + str(isite) + ':' + site.kind_name, + symbols=kind.symbols, + weights=kind.weights, + position=site.position ) # use auxiliary structure inside k-point generator output = get_explicit_kpoints_path(saux) @@ -307,6 +335,13 @@ def set_params_BS(self): """ params = self.ctx.input_params_KKR + # maybe overwrite some inputs + if 'params_kkr_overwrite' in self.inputs: + self.report(f'found params_kkr_overwrite: {self.inputs.params_kkr_overwrite.get_dict()}') + updatenode = self.inputs.params_kkr_overwrite + updatenode.label = 'params overwrite' + params = update_params_wf(params, updatenode) + input_dict = params.get_dict() para_check = kkrparams() try: @@ -336,7 +371,9 @@ def set_params_BS(self): descr = 'added missing default keys, ' ##+++ Starts to add the NTP2, EMAX and EMIN from the econt_new = self.ctx.BS_params_dict - econt_new['kmesh'] = [1, 1, 1] # overwrite kmesh since the kpoints are used from the input + if self.ctx.struc_is_alloy: + if econt_new['kmesh'] is None: + econt_new['kmesh'] = [1, 1, 1] # overwrite kmesh since the kpoints are used from the input kkr_calc = self.inputs.remote_data.get_incoming(node_class=KkrCalculation).first().node ef = kkr_calc.outputs.output_parameters.get_dict()['fermi_energy'] # unit in Ry self.ctx.fermi_energy = ef ## in Ry unit @@ -382,6 +419,12 @@ def get_BS(self): } if self.ctx.custom_scheduler_commands: options['custom_scheduler_commands'] = self.ctx.custom_scheduler_commands + if self.ctx.append_text: + options['append_text'] = self.ctx.append_text + if self.ctx.prepend_text: + options['prepend_text'] = self.ctx.prepend_text + if self.ctx.additional_retrieve_list: + options['additional_retrieve_list'] = self.ctx.additional_retrieve_list # get inputs for band structure calculation inputs = get_inputs_kkr( @@ -513,8 +556,7 @@ def set_energy_params(econt_new, ef, para_check): elif key in ['emax', 'EMAX']: key = 'EMAX' val = (ef + val / evscal) # Converting to the Ry (unit of the energy) - elif key in ['tempr' - 'TEMPR']: + elif key in ['tempr', 'TEMPR']: key = 'TEMPR' elif key in ['RCLUSTZ', 'rclustz']: key = 'RCLUSTZ' diff --git a/aiida_kkr/workflows/dos.py b/aiida_kkr/workflows/dos.py index d7664e97..faae55c3 100644 --- a/aiida_kkr/workflows/dos.py +++ b/aiida_kkr/workflows/dos.py @@ -29,7 +29,7 @@ __copyright__ = (u'Copyright (c), 2017, Forschungszentrum Jülich GmbH, ' 'IAS-1/PGI-1, Germany. All rights reserved.') __license__ = 'MIT license, see LICENSE.txt file' -__version__ = '0.8.2' +__version__ = '0.8.3' __contributors__ = u'Philipp Rüßmann' @@ -68,6 +68,9 @@ class kkr_dos_wc(WorkChain): 'max_wallclock_seconds': 60 * 60, 'withmpi': True, # execute KKR with mpi or without 'custom_scheduler_commands': '', # some additional scheduler commands + 'prepend_text': '', + 'append_text': '', + 'additional_retrieve_list': None } # intended to guide user interactively in setting up a valid wf_params node @@ -115,6 +118,13 @@ def define(cls, spec): help="""Initial non-collinear angles for the magnetic moments. See KkrCalculation for details. If this is found in the input potentially extracted nonco angles from the parent calulation are overwritten!""" ) + # maybe overwrite some settings from the KKRhost convergence run + spec.input( + 'params_kkr_overwrite', + valid_type=orm.Dict, + required=False, + help='Overwrite some input parameters of the parent KKR calculation.' + ) # define outputs spec.output('results_wf', valid_type=orm.Dict, required=True, help='Results collected by the workflow.') @@ -193,6 +203,11 @@ def start(self): if options_dict == {}: options_dict = self._options_default self.report('INFO: using default options') + self.ctx.append_text = options_dict.get('append_text', self._options_default['append_text']) + self.ctx.prepend_text = options_dict.get('prepend_text', self._options_default['prepend_text']) + self.ctx.additional_retrieve_list = options_dict.get( + 'additional_retrieve_list', self._options_default['additional_retrieve_list'] + ) # set values, or defaults self.ctx.withmpi = options_dict.get( @@ -354,6 +369,14 @@ def set_params_dos(self): 'KKR parameter node extracted from parent parameters and wf_parameter input node.' paranode_dos = update_params_wf(self.ctx.input_params_KKR, updatenode) + + # maybe overwrite some inputs + if 'params_kkr_overwrite' in self.inputs: + self.report(f'found params_kkr_overwrite: {self.inputs.params_kkr_overwrite.get_dict()}') + updatenode = self.inputs.params_kkr_overwrite + updatenode.label = 'params overwrite' + paranode_dos = update_params_wf(paranode_dos, updatenode) + self.ctx.dos_kkrparams = paranode_dos def get_dos(self): @@ -376,6 +399,12 @@ def get_dos(self): } if self.ctx.custom_scheduler_commands: options['custom_scheduler_commands'] = self.ctx.custom_scheduler_commands + if self.ctx.append_text: + options['append_text'] = self.ctx.append_text + if self.ctx.prepend_text: + options['prepend_text'] = self.ctx.prepend_text + if self.ctx.additional_retrieve_list: + options['additional_retrieve_list'] = self.ctx.additional_retrieve_list inputs = get_inputs_kkr( code, remote, options, label, description, parameters=params, serial=(not self.ctx.withmpi) diff --git a/aiida_kkr/workflows/jijs.py b/aiida_kkr/workflows/jijs.py index 73cdcdff..14094025 100644 --- a/aiida_kkr/workflows/jijs.py +++ b/aiida_kkr/workflows/jijs.py @@ -19,7 +19,7 @@ __copyright__ = (u'Copyright (c), 2022, Forschungszentrum Jülich GmbH, ' 'IAS-1/PGI-1, Germany. All rights reserved.') __license__ = 'MIT license, see LICENSE.txt file' -__version__ = '0.1.0' +__version__ = '0.1.1' __contributors__ = (u'Philipp Rüßmann') @@ -359,25 +359,25 @@ def _get_para_jij(self, params): para_jij.set_value('NSTEPS', 1) # one-shot run # accuracy settings - tempr = self.ctx.wf_dict['accuracy'].get('TEMPR') + tempr = self.ctx.wf_dict.get('accuracy', {}).get('TEMPR') if tempr is not None: para_jij.set_value('TEMPR', tempr) # slightly reduce temperature - rclustz = self.ctx.wf_dict['accuracy'].get('RCLUSTZ') + rclustz = self.ctx.wf_dict.get('accuracy', {}).get('RCLUSTZ') if rclustz is not None: para_jij.set_value('RCLUSTZ', rclustz) # increase cluster radius - kmesh = self.ctx.wf_dict['accuracy'].get('kmesh') + kmesh = self.ctx.wf_dict.get('accuracy', {}).get('kmesh') if kmesh is not None: para_jij.set_value('BZDIVIDE', kmesh) # increase k-points para_jij.set_value('KPOIBZ', np.product(kmesh)) # array dimension # array dimensions - NATOMIMPD = self.ctx.wf_dict['accuracy'].get('NATOMIMPD') + NATOMIMPD = self.ctx.wf_dict.get('accuracy', {}).get('NATOMIMPD') if NATOMIMPD is not None: para_jij.set_value('NATOMIMPD', NATOMIMPD) # array dimension - NSHELD = self.ctx.wf_dict['accuracy'].get('NSHELD') + NSHELD = self.ctx.wf_dict.get('accuracy', {}).get('NSHELD') if NSHELD is not None: para_jij.set_value('NSHELD', NSHELD) # array dimension @@ -389,8 +389,8 @@ def _get_para_jij(self, params): # set optional Jij parameters # i and j index for Jij calculation in internal units # uses site index (i.e. needs to be <=10) - JIJSITEI = self.ctx.wf_dict['jijsite_i'] - JIJSITEJ = self.ctx.wf_dict['jijsite_j'] + JIJSITEI = self.ctx.wf_dict.get('jijsite_i') + JIJSITEJ = self.ctx.wf_dict.get('jijsite_j') if JIJSITEI is not None: if JIJSITEJ is None: JIJSITEJ = JIJSITEI @@ -406,7 +406,7 @@ def _get_jijrad(self): """ # Jij radius in Ang. - jijrad_ang = self.ctx.wf_dict['jijrad_ang'] + jijrad_ang = self.ctx.wf_dict.get('jijrad_ang') # find structure from calculation struc, _ = find_parent_structure(self.ctx.parent_calc) diff --git a/aiida_kkr/workflows/kkr_imp.py b/aiida_kkr/workflows/kkr_imp.py index 13268950..3560212e 100644 --- a/aiida_kkr/workflows/kkr_imp.py +++ b/aiida_kkr/workflows/kkr_imp.py @@ -542,6 +542,7 @@ def run_voroaux(self): kkrcode = self.inputs.kkr imp_info = self.inputs.impurity_info voro_params = self.ctx.voro_params_dict + if self.ctx.do_gf_calc: self.report('INFO: get converged host remote from inputs to extract structure for Voronoi calculation') converged_host_remote = self.inputs.remote_data_host @@ -553,6 +554,9 @@ def run_voroaux(self): GF_host_calc = remote_data_gf_node.get_incoming(link_label_filter=u'remote_folder').first().node converged_host_remote = GF_host_calc.inputs.parent_folder + # find host structure + structure_host, voro_calc = VoronoiCalculation.find_parent_structure(converged_host_remote) + # get previous kkr parameters following remote_folder->calc->parameters links prev_kkrparams = converged_host_remote.get_incoming(link_label_filter='remote_folder' ).first().node.get_incoming(link_label_filter='parameters' @@ -570,10 +574,24 @@ def run_voroaux(self): # add or overwrite some parameters (e.g. things that are only used by voronoi) calc_params_dict = calc_params.get_dict() # add some voronoi specific parameters automatically if found (RMTREF should also set RMTCORE to the same value) - if '' in list(calc_params_dict.keys()): + if calc_params_dict.get('', None) is not None: self.report('INFO: add rmtcore to voro params') self.ctx.change_voro_params[''] = calc_params_dict[''] self.report(self.ctx.change_voro_params) + + # add some voronoi-specific settings starting from host's voronoi run (might have been overwritten in between) + # this is necessary to make sure that voronoi creates the same radial mesh for the impurity potential, otherwise KKRimp will fail + if voro_calc.inputs.parameters.get_dict().get('RUNOPT', None) is not None: + self.report("INFO: copy runopt from host's voronoi run") + runopt = self.ctx.change_voro_params.get('RUNOPT', None) + if runopt is None: + runopt = [] + runopt += voro_calc.inputs.parameters['RUNOPT'] + self.ctx.change_voro_params['RUNOPT'] = runopt + if voro_calc.inputs.parameters.get_dict().get('RCLUSTZ', None) is not None: + self.report("INFO: copy RCLUSTZ from host's voronoi run") + self.ctx.change_voro_params['RCLUSTZ'] = voro_calc.inputs.parameters['RCLUSTZ'] + changed_params = False for key, val in self.ctx.change_voro_params.items(): if key in ['RUNOPT', 'TESTOPT']: @@ -589,9 +607,6 @@ def run_voroaux(self): updatenode.description = 'Overwritten voronoi input parameter from kkr_imp_wc input.' calc_params = update_params_wf(calc_params, updatenode) - # find host structure - structure_host, voro_calc = VoronoiCalculation.find_parent_structure(converged_host_remote) - # for every impurity, generate a structure and launch the voronoi workflow # to get the auxiliary impurity startpotentials self.ctx.voro_calcs = {} diff --git a/aiida_kkr/workflows/kkr_imp_sub.py b/aiida_kkr/workflows/kkr_imp_sub.py index 72e0af42..753002cb 100644 --- a/aiida_kkr/workflows/kkr_imp_sub.py +++ b/aiida_kkr/workflows/kkr_imp_sub.py @@ -290,6 +290,7 @@ def start(self): self.ctx.nsteps = wf_dict.get('nsteps', self._wf_default['nsteps']) self.ctx.broyden_num = wf_dict.get('broyden-number', self._wf_default['broyden-number']) self.ctx.nsimplemixfirst = wf_dict.get('nsimplemixfirst', self._wf_default['nsimplemixfirst']) + self.ctx.mesh_params = wf_dict.get('accuracy_params', {}) # initial magnetization self.ctx.mag_init = wf_dict.get('mag_init', self._wf_default['mag_init']) diff --git a/pyproject.toml b/pyproject.toml index 62ef5dc8..2fa345b8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -3,7 +3,7 @@ requires = ["setuptools", "wheel"] [tool.poetry] name = "aiida-kkr" -version = "1.1.12" +version = "1.1.13" description = "AiiDA plugin for the KKR code" authors = ["Philipp Rüssmann ", "Jens Bröder ", diff --git a/setup.json b/setup.json index d5e6a06a..79582dcc 100644 --- a/setup.json +++ b/setup.json @@ -20,7 +20,7 @@ "Natural Language :: English", "Framework :: AiiDA" ], - "version": "1.1.12", + "version": "1.1.13", "reentry_register": true, "install_requires": [ "aiida-core >= 1.0.0b6,<3.0.0",