diff --git a/aiida_kkr/calculations/kkr.py b/aiida_kkr/calculations/kkr.py index 490d7e66..8d7724f4 100644 --- a/aiida_kkr/calculations/kkr.py +++ b/aiida_kkr/calculations/kkr.py @@ -19,6 +19,7 @@ ) from masci_tools.io.common_functions import get_alat_from_bravais, get_Ang2aBohr from aiida_kkr.tools.tools_kkrimp import make_scoef, write_scoef_full_imp_cls +from aiida_kkr.tools.find_parent import find_parent_structure from masci_tools.io.kkr_params import __kkr_default_params__, kkrparams import six from six.moves import range @@ -443,7 +444,7 @@ def _get_structure_inputs(self, parent_calc, parameters): structure = None self.logger.info('KkrCalculation: Get structure node from voronoi parent') try: - structure, voro_parent = VoronoiCalculation.find_parent_structure(parent_calc) + structure, voro_parent = find_parent_structure(parent_calc) except: self.logger.error(f'KkrCalculation: Could not get structure from Voronoi parent ({parent_calc}).') raise ValidationError(f'Cound not find structure node from parent {parent_calc}') @@ -626,7 +627,7 @@ def _get_local_copy_list(self, parent_calc, voro_parent, tempfolder, parameters) # add shapefun file from voronoi parent if needed if self._SHAPEFUN not in copylist: try: - struc, voro_parent = VoronoiCalculation.find_parent_structure(parent_calc) + struc, voro_parent = find_parent_structure(parent_calc) except ValueError: return self.exit_codes.ERROR_NO_SHAPEFUN_FOUND # pylint: disable=no-member # copy shapefun from retrieved of voro calc diff --git a/aiida_kkr/calculations/kkrnano.py b/aiida_kkr/calculations/kkrnano.py index 2880ea21..2764ad74 100644 --- a/aiida_kkr/calculations/kkrnano.py +++ b/aiida_kkr/calculations/kkrnano.py @@ -2,22 +2,21 @@ """ Input plug-in for a KKRnano calculation. """ -from aiida.engine import CalcJob -from aiida.orm import CalcJobNode, Dict, Bool, Float, RemoteData, StructureData -from aiida_kkr.data.strucwithpot import StrucWithPotData -from aiida.common.datastructures import (CalcInfo, CodeInfo) + +import numpy as np from masci_tools.io.common_functions import get_Ang2aBohr +from aiida.orm import CalcJobNode, Dict, Bool, Float, RemoteData, StructureData +from aiida.engine import CalcJob from aiida.common import NotExistent +from aiida.common.datastructures import (CalcInfo, CodeInfo) from aiida.common.exceptions import InputValidationError, UniquenessError - -from aiida_kkr.calculations.voro import VoronoiCalculation - -import numpy as np +from aiida_kkr.tools.find_parent import get_remote, get_parent +from aiida_kkr.data.strucwithpot import StrucWithPotData __copyright__ = (u'Copyright (c), 2021, Forschungszentrum Jülich GmbH, ' 'IAS-1/PGI-1, Germany. All rights reserved.') __license__ = 'MIT license, see LICENSE.txt file' -__version__ = '0.0.1' +__version__ = '0.0.2' __contributors__ = ('Markus Struckmann', 'Philipp Rüßmann') @@ -704,9 +703,9 @@ def find_parent_struc_from_voro_or_stwpd(self, parent_folder): """ iiter = 0 Nmaxiter = 1000 - parent_folder_tmp = VoronoiCalculation.get_remote(parent_folder) + parent_folder_tmp = get_remote(parent_folder) while not self._has_struc(parent_folder_tmp) and iiter < Nmaxiter: - parent_folder_tmp = VoronoiCalculation.get_remote(VoronoiCalculation.get_parent(parent_folder_tmp)) + parent_folder_tmp = get_remote(get_parent(parent_folder_tmp)) iiter += 1 if iiter % 200 == 0: print( diff --git a/aiida_kkr/calculations/voro.py b/aiida_kkr/calculations/voro.py index ec5c7037..5ec5520c 100644 --- a/aiida_kkr/calculations/voro.py +++ b/aiida_kkr/calculations/voro.py @@ -10,7 +10,7 @@ from aiida.common.exceptions import (InputValidationError, ValidationError) from aiida.common.datastructures import (CalcInfo, CodeInfo) from aiida_kkr.tools.common_workfunctions import generate_inputcard_from_structure, check_2Dinput_consistency, vca_check -from aiida.common.exceptions import UniquenessError, NotExistent +from aiida.common.exceptions import UniquenessError import os import six @@ -275,86 +275,11 @@ def _is_KkrCalc(self, calc): return is_KKR - @classmethod - def _get_struc(self, parent_calc): - """ - Get structure from a parent_folder (result of a calculation, typically a remote folder) - """ - return parent_calc.inputs.structure - - @classmethod - def _has_struc(self, parent_folder): - """ - Check if parent_folder has structure information in its input - """ - success = True - if 'structure' not in parent_folder.get_incoming().all_link_labels(): - success = False - return success - - @classmethod - def get_remote(self, parent_folder): - """ - get remote_folder from input if parent_folder is not already a remote folder - """ - parent_folder_tmp0 = parent_folder - try: - parent_folder_tmp = parent_folder_tmp0.get_incoming().get_node_by_label('remote_folder') - except NotExistent: - parent_folder_tmp = parent_folder_tmp0 - return parent_folder_tmp - - @classmethod - def get_parent(self, input_folder): - """ - get the parent folder of the calculation. If not parent was found return input folder - """ - input_folder_tmp0 = input_folder - - # first option: parent_calc_folder (KkrimpCalculation) - try: - parent_folder_tmp = input_folder_tmp0.get_incoming().get_node_by_label('parent_calc_folder') - return_input = False - except NotExistent: - return_input = True - - # second option: parent_folder (KkrCalculation) - try: - parent_folder_tmp = input_folder_tmp0.get_incoming().get_node_by_label('parent_folder') - return_input = False - except NotExistent: - return_input = return_input & True - - # third option: parent_KKR option (special mode of VoronoiCalculation) - try: - parent_folder_tmp = input_folder_tmp0.get_incoming().get_node_by_label('parent_KKR') - return_input = False - except NotExistent: - return_input = return_input & True - - if return_input: - parent_folder_tmp = input_folder_tmp0 - - return parent_folder_tmp - @classmethod def find_parent_structure(self, parent_folder): """ Find the Structure node recuresively in chain of parent calculations (structure node is input to voronoi calculation) + This is a copy of the find_parent_structure that moved to tools.find_parent to keep backwards compatibility. """ - iiter = 0 - Nmaxiter = 1000 - parent_folder_tmp = self.get_remote(parent_folder) - while not self._has_struc(parent_folder_tmp) and iiter < Nmaxiter: - parent_folder_tmp = self.get_remote(self.get_parent(parent_folder_tmp)) - iiter += 1 - if iiter % 200 == 0: - print( - 'Warning: find_parent_structure takes quite long (already searched {} ancestors). Stop after {}'. - format(iiter, Nmaxiter) - ) - if self._has_struc(parent_folder_tmp): - struc = self._get_struc(parent_folder_tmp) - return struc, parent_folder_tmp - else: - raise ValueError('structure not found') + from aiida_kkr.tools.find_parent import find_parent_structure + return find_parent_structure(parent_folder) diff --git a/aiida_kkr/data/strucwithpot.py b/aiida_kkr/data/strucwithpot.py index 702d5a90..9686f443 100644 --- a/aiida_kkr/data/strucwithpot.py +++ b/aiida_kkr/data/strucwithpot.py @@ -1,8 +1,7 @@ from aiida.orm import CalcJobNode, Dict, StructureData, SinglefileData, load_node, Data from aiida.common.exceptions import InputValidationError from aiida.common import NotExistent -#from aiida_kkr.tools import find_parent_structure #deprecated -from aiida_kkr.calculations.voro import VoronoiCalculation +from aiida_kkr.tools.find_parent import find_parent_structure, get_parent, get_remote import os @@ -173,7 +172,7 @@ def get_strucwithpot_from_Voronoi(self, calcnode): """ #TODO: Check if input is valid - structure = VoronoiCalculation.find_parent_structure(calcnode)[0] + structure = find_parent_structure(calcnode)[0] cwd = os.getcwd() #shapefun @@ -275,7 +274,7 @@ def get_strucwithpot_from_KKRnanoConvert(self, calcnode): raise InputValidationError( 'Only the convert step output can be processed! If this has not been done, yet, the parameter `builder.convert=Bool(False)` can be used and the process be run with 1 MPI, to obtain ASCII-potential files' ) - structure = VoronoiCalculation.find_parent_structure(calcnode)[0] + structure = find_parent_structure(calcnode)[0] pot_list = [] for item in calcnode.outputs.retrieved.list_object_names(): if item.find('vpot') == 0: @@ -299,32 +298,6 @@ def get_strucwithpot_from_KKRnanoConvert(self, calcnode): return structure, shape_list, pot_list - def _get_remote(self, parent_folder): - """ - get remote_folder from input if parent_folder is not already a remote folder - """ - parent_folder_tmp0 = parent_folder - try: - parent_folder_tmp = parent_folder_tmp0.get_incoming().get_node_by_label('remote_folder') - except NotExistent: - parent_folder_tmp = parent_folder_tmp0 - return parent_folder_tmp - - #TODO remove redundancy aiida_kkr/calculations/voro.py -> There a private file - def _get_parent(self, input_folder): - """ - get the parent folder of the calculation. If not parent was found return input folder - """ - input_folder_tmp0 = input_folder - try: - parent_folder_tmp = input_folder_tmp0.get_incoming().get_node_by_label('parent_calc_folder') - except NotExistent: - try: - parent_folder_tmp = input_folder_tmp0.get_incoming().get_node_by_label('parent_folder') - except NotExistent: - parent_folder_tmp = input_folder_tmp0 - return parent_folder_tmp - def find_parent_shapefun(self, parent_folder): """ Find the shape files recursively in chain of parent calculations, either to be extracted from "shapefun" file or "shapes" files @@ -332,15 +305,13 @@ def find_parent_shapefun(self, parent_folder): iiter = 0 Nmaxiter = 1000 - parent_folder_tmp = self._get_parent(parent_folder) + parent_folder_tmp = get_parent(parent_folder) print(parent_folder_tmp) parent_folder_tmp_listdir = parent_folder_tmp.listdir( ) #requires remote ssh connection, therefore much quicker this way while not 'shape.0000001' in parent_folder_tmp_listdir and not 'shapefun' in parent_folder_tmp_listdir and iiter < Nmaxiter: - parent_folder_tmp = self._get_remote( - self._get_parent(parent_folder_tmp) - ) #at this point the result is a CalcNode! + parent_folder_tmp = get_remote(get_parent(parent_folder_tmp)) #at this point the result is a CalcNode! parent_folder_tmp_listdir = parent_folder_tmp.outputs.remote_folder.listdir() iiter += 1 print(parent_folder_tmp) diff --git a/aiida_kkr/parsers/kkrnano.py b/aiida_kkr/parsers/kkrnano.py index 80fd4db5..b0fd3abc 100644 --- a/aiida_kkr/parsers/kkrnano.py +++ b/aiida_kkr/parsers/kkrnano.py @@ -12,7 +12,6 @@ from aiida.parsers.parser import Parser from aiida.orm import Dict, CalcJobNode from aiida_kkr.calculations.kkrnano import KKRnanoCalculation -#from aiida_kkr.tools import find_parent_structure from masci_tools.io.common_functions import (search_string, open_general) import numpy as np from io import StringIO diff --git a/aiida_kkr/tools/__init__.py b/aiida_kkr/tools/__init__.py index 6dbef3a1..a906ffa6 100644 --- a/aiida_kkr/tools/__init__.py +++ b/aiida_kkr/tools/__init__.py @@ -17,6 +17,7 @@ from .multi_imps_data_extract import MultiImpuritiesData from .kick_out_core_states import * from .neworder_potential import * +from .find_parent import get_calc_from_remote, get_remote, get_parent # expose structure finder from VoronoiCalculation @@ -29,8 +30,11 @@ def find_parent_structure(calc_or_remote, return_voro=False): :param calc_or_remote: CalcJobNode or RemoteData node of VoronoiCalculation or KkrCalculation :return struc: parent StructureData node """ - from aiida_kkr.calculations.voro import VoronoiCalculation + from .find_parent import find_parent_structure - struc, voro_calc = VoronoiCalculation.find_parent_structure(calc_or_remote) + struc, voro_calc = find_parent_structure(calc_or_remote) - return struc + if return_voro: + return struc, voro_calc + else: + return struc diff --git a/aiida_kkr/tools/find_parent.py b/aiida_kkr/tools/find_parent.py new file mode 100644 index 00000000..2b7431eb --- /dev/null +++ b/aiida_kkr/tools/find_parent.py @@ -0,0 +1,107 @@ +""" +Helper tools to find parent calculations or structures. +""" + +from aiida.common.exceptions import NotExistent + + +def get_struc(parent_calc): + """ + Get structure from a parent_folder (result of a calculation, typically a remote folder) + """ + return parent_calc.inputs.structure + + +def has_struc(parent_folder): + """ + Check if parent_folder has structure information in its input + """ + success = True + if 'structure' not in parent_folder.get_incoming().all_link_labels(): + success = False + return success + + +def get_remote(parent_folder): + """ + get remote_folder from input if parent_folder is not already a remote folder + """ + parent_folder_tmp0 = parent_folder + try: + parent_folder_tmp = parent_folder_tmp0.get_incoming().get_node_by_label('remote_folder') + except NotExistent: + parent_folder_tmp = parent_folder_tmp0 + return parent_folder_tmp + + +def get_parent(input_folder): + """ + Get the parent folder of the calculation. If no parent was found return input folder + """ + input_folder_tmp0 = input_folder + + # first option: parent_calc_folder (KkrimpCalculation) + try: + parent_folder_tmp = input_folder_tmp0.get_incoming().get_node_by_label('parent_calc_folder') + return_input = False + except NotExistent: + return_input = True + + # second option: parent_folder (KkrCalculation) + try: + parent_folder_tmp = input_folder_tmp0.get_incoming().get_node_by_label('parent_folder') + return_input = False + except NotExistent: + return_input = return_input & True + + # third option: parent_KKR option (special mode of VoronoiCalculation) + try: + parent_folder_tmp = input_folder_tmp0.get_incoming().get_node_by_label('parent_KKR') + return_input = False + except NotExistent: + return_input = return_input & True + + if return_input: + parent_folder_tmp = input_folder_tmp0 + + return parent_folder_tmp + + +def find_parent_structure(parent_folder): + """ + Find the Structure node recuresively in chain of parent calculations (structure node is input to voronoi calculation) + """ + iiter = 0 + Nmaxiter = 1000 + parent_folder_tmp = get_remote(parent_folder) + while not has_struc(parent_folder_tmp) and iiter < Nmaxiter: + parent_folder_tmp = get_remote(get_parent(parent_folder_tmp)) + iiter += 1 + if iiter % 200 == 0: + print( + 'Warning: find_parent_structure takes quite long (already searched {} ancestors). Stop after {}'.format( + iiter, Nmaxiter + ) + ) + if has_struc(parent_folder_tmp): + struc = get_struc(parent_folder_tmp) + return struc, parent_folder_tmp + else: + raise ValueError('structure not found') + + +def get_calc_from_remote(calc_remote): + """ + Get the parent calculation from a RemoteData + """ + from aiida.orm import CalcJobNode + from aiida.orm import RemoteData + + if not isinstance(calc_remote, RemoteData): + raise ValueError('input node is not a RemoteData folder') + + parents = calc_remote.get_incoming(node_class=CalcJobNode).all() + if len(parents) != 1: + raise ValueError('Parent is not unique!') + + return parents[0].node diff --git a/aiida_kkr/tools/jij_tools.py b/aiida_kkr/tools/jij_tools.py new file mode 100644 index 00000000..1f02951e --- /dev/null +++ b/aiida_kkr/tools/jij_tools.py @@ -0,0 +1,567 @@ +""" +Tools for Jij calculations and parsing +""" + +# parsing jij output + +import numpy as np +from masci_tools.io.common_functions import search_string, get_aBohr2Ang, get_Ry2eV +from aiida.orm import ArrayData, StructureData, Dict, load_node, Bool +from aiida.engine import calcfunction, submit +from aiida_kkr.calculations import KkrCalculation +from aiida_kkr.tools import find_parent_structure + + +def get_sites(structure): + """Get all sites also for a CPA structure""" + sites = [] # for CPA + for site in structure.sites: + sitekind = structure.get_kind(site.kind_name) + for ikind in range(len(sitekind.symbols)): + sites.append(site) + return sites + + +def get_jijs_shells(jij_calc, verbose=False): + """read the jij.atom files from the retrieved and determine if the file has DMI or not""" + + # read the jij.atom files + jijs_shells_all = [] + for jijfile in [i for i in jij_calc.outputs.retrieved.list_object_names() if 'Jij' in i]: + i_index = int(jijfile.split('.atom')[1]) + if verbose: + print('load jij.atom', i_index) + with jij_calc.outputs.retrieved.open(jijfile) as f: + try: + tmp = np.loadtxt(f) + except: + f.seek(0) + txt = f.readlines() + tmp = np.loadtxt(txt[:-1]) + # sort by radius + jij_atom = tmp[tmp[:, 0].argsort()] + if len(jij_atom[0]) > 4: + # add i index, needed for mapping to spirit data + # only needed for DMI mode + jij_atom[:, -1] = i_index + if verbose: + print('shape jij.atom', jij_atom.shape) + jijs_shells_all += jij_atom.tolist() + + if verbose: + print('jijs_shells.shape', len(jijs_shells_all)) + + if len(jijs_shells_all) == 0: + raise ValueError('Did not extract any Jij shells') + + # format of jijs_shells differs in old and new solver: + if len(jijs_shells_all[0]) == 4: + # - isotropic exchange constants only (old solver). The columns refer to: + #  * [0] |Rij| in units of the lattice constant + # * [1] Jij in Ryd + # * [2] shell-number + # * [3] atom type of atom j + # * [4] atom type of atom i + dmimode = False + else: + # - full exchange tensor (new solver). The columns refer to: + # * [0] ∣Rij∣ in units of the lattice constant + # * [1] Jij (isotropic part) in Ryd + # * [2] Dij (anti-symmetric DMI part) in Ryd + # * [3] Sij (diagonal traceless part) in Ryd + # * [4] Aij (off-diagonal symmetric part) in Ryd + # * [5-7] Rj− Ri (3 component vector) + # * [8] atom type of atom j + # * [9] atom type of atom i + dmimode = True + + return jijs_shells_all, dmimode + + +def expand_jijs_iso(jij_calc, jijs_shells_all, shells, alat): + """expand jijs from shells to all pairs which is easier to use in spirit + This is the roune that deals with isotropic interactions (i.e. using the old solver)""" + + # extract basis vectors ot the structure, needed to convert from cartesian coordinates (x,y,z) + # to relative coordinates (da, db, dc) such that (x, y, z)^T = da*a + db*b + dc*c + # with a,b,c being the three unit vectors + structure = find_parent_structure(jij_calc) + cell = np.array(structure.cell) + cell_positions = np.array([i.position for i in structure.sites]) + + # transformation matrix from absolute to realtive coordinates + U2rel = np.linalg.inv(cell.transpose()) + + # expand shells data to all structure that spirit can understand (i.e. take all pairs) + jijs_expanded, positions_expanded = [], [] + for jshell in jijs_shells_all: + for s in shells[int(jshell[3] - 1)]: + # i and j indices in the unit cell + iatom, jatom = int(s[0] - 1), int(s[1] - 1) + + # absolute distance + # Rx, Ry, Rz = alat*(s[7]-s[4]), alat*(s[8]-s[5]), alat*(s[9]-s[6]) + Rx, Ry, Rz = alat * s[7], alat * s[8], alat * s[9] + + # substract sub-lattice distance to find correct unit cell scalings + R = np.array([Rx, Ry, Rz]) + R -= cell_positions[jatom] + # R -= (cell_positions[jatom] - cell_positions[i]) + x, y, z = R + + # calculate da, db, dc ( i.e. the mutiplicities of the unit vectors) + da, db, dc = np.array(np.round(np.dot(U2rel, [x, y, z])), dtype=int) + + # i, j, da, db, dc, Jij (meV) + jijs_expanded.append([iatom, jatom, da, db, dc, jshell[1] * get_Ry2eV() * 1000]) + + positions_expanded.append(np.array([Rx, Ry, Rz]) - cell_positions[iatom]) + + # convert to numpy array + jijs_expanded = np.array(jijs_expanded) + positions_expanded = np.array(positions_expanded) + + return jijs_expanded, positions_expanded + + +def expand_jijs_dmi(jijs_shells_x, jijs_shells_y, jijs_shells_z, shells, cell, sites, alat, verbose=False): + """Bring output Jijs of the new solver into the right form + (KKR output already is expanded in shells but it needs to be brought into the right form + by combining the x,y,z calculations and the Rvec information) + This will give the full Jij tensor and prodiuce the columns that spirit can understand + """ + + # conversion factor from Ry to meV + Ry2meV = 1000 * get_Ry2eV() + + # transformation matrix from absolute to realtive coordinates + U2rel = np.linalg.inv(cell.transpose()) + + jijs_expanded, positions_expanded = [], [] + for ipos, jz in enumerate(jijs_shells_z): + # also load the information from m||x and m||y + jx = jijs_shells_x[ipos] + jy = jijs_shells_y[ipos] + + if verbose: + print(ipos, jx, jy, jz) + + # get i and j indices of the atoms + # -1 to convert from fortran output to python indexing (starting at 0) + iatom = int(jz[-1]) - 1 + jatom = int(jz[-2]) - 1 + + # Rvec = pos_i - pos_j + R_ij[j], with R_ij[j] given in the Jij.atom files (in alat units) + dx = alat * jz[5] + dy = alat * jz[6] + dz = alat * jz[7] + x = sites[iatom].position[0] - sites[jatom].position[0] + dx + y = sites[iatom].position[1] - sites[jatom].position[1] + dy + z = sites[iatom].position[2] - sites[jatom].position[2] + dz + + # calculate da, db, dc ( i.e. the mutiplicities of the unit vectors) + da, db, dc = np.array(np.round(np.dot(U2rel, [x, y, z])), dtype=int) + + if verbose: + if ipos == 0: + print(cell) + print(iatom, jatom, dx, dy, dz, x, y, z, da, db, dc) + + #collect full Jij tensor + # 1: J, 2: D, 3: S, 4: A + Jxx = (jz[1] + jz[3] + jy[1] + jy[3]) / 2. * Ry2meV + Jxy = jz[4] * Ry2meV + Jxz = jy[4] * Ry2meV + Jyx = jz[4] * Ry2meV + Jyy = (jz[1] - jz[3] + jx[1] + jx[3]) / 2. * Ry2meV + Jyz = jx[4] * Ry2meV + Jzx = jy[4] * Ry2meV + Jzy = jx[4] * Ry2meV + Jzz = (jx[1] - jx[3] + jy[1] - jy[3]) / 2. * Ry2meV + # calculate DMI vector + Dx = jx[2] * Ry2meV + # for y component no minus sign: ERROR in documentation! + #Dy = -jy[2] * Ry2meV + Dy = jy[2] * Ry2meV + Dz = jz[2] * Ry2meV + + # collect data in big arrays + # i, j, da, db, dc, Jij (meV), Dij vector(meV), full Jij tensor (meV) + jijs_expanded.append([ + iatom, jatom, da, db, dc, (Jxx + Jyy + Jzz) / 3., Dx, Dy, Dz, Jxx, Jxy, Jxz, Jyx, Jyy, Jyz, Jzx, Jzy, Jzz + ]) + positions_expanded.append([dx, dy, dz]) + + # convert collected data to numpy array + jijs_expanded = np.array(jijs_expanded) + positions_expanded = np.array(positions_expanded) + + return jijs_expanded, positions_expanded + + +def get_jij_structure(structure, jijs_expanded, jij_calc): + """make auxiliary structure that has only the ij sites from the Jij calculation""" + struc_jij_sites = StructureData(cell=structure.cell) + struc_jij_sites.pbc = structure.pbc + all_sites_jij = set(list(jijs_expanded[:, 0]) + list(jijs_expanded[:, 1])) + isite, icount, mappings = -1, 0, [] # for mapping to the sites of the reduced structure + mappings_back = {} + for site in structure.sites: + sitekind = structure.get_kind(site.kind_name) + for ikind, symbol in enumerate(sitekind.symbols): + isite += 1 + # take only structues for which Jij couplings are extracted + if isite in all_sites_jij: + mappings.append([icount, isite]) + mappings_back[isite] = icount + icount += 1 + if ikind == 0: + struc_jij_sites.append_atom( + position=site.position, symbols=sitekind.symbols, weights=sitekind.weights + ) + + # add spin moments extra (used by spirit plugin) + calc_conv = jij_calc.inputs.parent_folder.get_incoming(node_class=KkrCalculation).first().node + mu_s_all = calc_conv.outputs.output_parameters['magnetism_group']['spin_moment_per_atom'] + mu_s = [] + for ij in mappings: + mu_s.append(mu_s_all[ij[1]]) + struc_jij_sites.set_extra('spin_moments', mu_s) + + return struc_jij_sites, mappings_back, mu_s + + +def get_shells(jij_calc, verbose=False): + # read the shells information that is needed to map to the complete list of pairs + with jij_calc.outputs.retrieved.open('shells.dat') as f: + txt = f.readlines() + nshell = int(txt.pop(0).split()[0]) + naez = len(find_parent_structure(jij_calc).sites) + ioffset = 0 + shells = [] + for ishell in range(nshell): + nat = int(txt[ioffset].split()[1]) + shell = [] + for iline in range(nat): + shell.append(txt[ioffset + 1 + iline].split()) + shells.append(np.array(shell, dtype=float)) + ioffset += 1 + nat + if verbose: + print('found shells:', shells) + return shells + + +@calcfunction +def parse_jij_calc_data( + jij_calc_retrieved, jij_calc_x_retrieved=None, jij_calc_y_retrieved=None, verbose=lambda: Bool(False) +): + """ + Parse the output of a Jij calculation from the retreived folder of a KkrCalculation + + :params jij_calc_retrieved: retrieved folder output of a KkrCalculation which ran with the Jij inputs (m||z is assumed) + :params jij_calc_x_retrieved: like jij_calc_retrieved but for m||x (only needed for new solver) + :params jij_calc_y_retrieved: like jij_calc_retrieved but for m||y (only needed for new solver) + :params verbose: True/False can be used to print debugging output + + :returns: + {'jij_data': jij_data, 'structure_jij_sites': struc_jij_sites} + where + * jij_data is and ArrayData that contains the expanded Jij's (see 'array_descriptions' extra for more details) + * struc_jij_sites is the reduced structure that contains onlt the atoms which have Jij couplings + (comes from the input of the KkrCalculation). The mappings to the original structure (and their i,j indices) is given as an extra + """ + verbose = verbose.value + + # extract kkrCalculation from retreived child + jij_calc = jij_calc_retrieved.get_incoming(node_class=KkrCalculation).first().node + if verbose: + print('jij(z) calculation:', jij_calc.uuid) + + # extract basis vectors ot the structure, needed to convert from cartesian coordinates (x,y,z) + # to relative coordinates (da, db, dc) such that (x, y, z)^T = da*a + db*b + dc*c + # with a,b,c being the three unit vectors + structure = find_parent_structure(jij_calc) + cell = np.array(structure.cell) + natyp = len(get_sites(structure)) + if verbose: + print('found structure with natyp:', natyp) + + # in KKR everything is scaled with the lattice constant + alat = jij_calc.outputs.output_parameters['alat_internal'] * get_aBohr2Ang() + + # read jij.atom files + jijs_shells_z, dmimode = get_jijs_shells(jij_calc, verbose) + if verbose: + print('shape jij_shells_z:', len(jijs_shells_z), dmimode) + + # read jij.atom files if a calculation in x and y are given in the input + dmimode_x, dmimode_y = False, False + if dmimode and jij_calc_x_retrieved is not None: + jij_calc_x = jij_calc_x_retrieved.get_incoming(node_class=KkrCalculation).first().node + if verbose: + print('jij(x) calculation:', jij_calc_x.uuid) + jijs_shells_x, dmimode_x = get_jijs_shells(jij_calc_x) + if verbose: + print('shape jij_shells_x:', len(jijs_shells_x), dmimode_x) + if not dmimode_x: + raise ValueError('jij_calc_x is not a DMI calculation (i.e. used old solver)') + if dmimode and jij_calc_y_retrieved is not None: + jij_calc_y = jij_calc_y_retrieved.get_incoming(node_class=KkrCalculation).first().node + if verbose: + print('jij(y) calculation:', jij_calc_y.uuid) + jijs_shells_y, dmimode_y = get_jijs_shells(jij_calc_y) + if verbose: + print('shape jij_shells_y:', len(jijs_shells_y), dmimode_y) + if not dmimode_y: + raise ValueError('jij_calc_y is not a DMI calculation (i.e. used old solver)') + # consistency check + if dmimode and not (dmimode_x and dmimode_y): + raise ValueError('Found dmimode but not all calculations (x,y and z) were given correctly') + + # read the shells information + shells = get_shells(jij_calc, verbose) + # print('shells', shells) + + # expand shells data to a structure that spirit can understand (i.e. take all pairs) + if not dmimode: + if verbose: + print('expand isotropic Jijs') + jijs_expanded, positions_expanded = expand_jijs_iso(jij_calc, jijs_shells_z, shells, alat) + else: + if verbose: + print('expand anisotropic Jijs, Dij etc') + jijs_expanded, positions_expanded = expand_jijs_dmi( + jijs_shells_x, jijs_shells_y, jijs_shells_z, shells, cell, get_sites(structure), alat, verbose=verbose + ) + + # create an auxiliary structure that contains only the sites which are used in the Jij step + # (i.e. we drop all sites where we don't have couplings) + struc_jij_sites, mappings_back, mu_s = get_jij_structure(structure, jijs_expanded, jij_calc) + + if verbose: + print(f'reduced structure has {len(struc_jij_sites.sites)} sites') + + # change i,j indices to match smaller structure + # and correct sign if mu_i and mu_j are oriented antiparallel + for ii, iatom in enumerate(jijs_expanded[:, 0]): + jatom = jijs_expanded[ii, 1] + jijs_expanded[ii, 0] = mappings_back[int(iatom)] + jijs_expanded[ii, 1] = mappings_back[int(jatom)] + # correct sign + smu_i = np.sign(mu_s[int(jijs_expanded[ii, 0])]) + smu_j = np.sign(mu_s[int(jijs_expanded[ii, 1])]) + sign_factor = smu_i * smu_j + jijs_expanded[ii, 5] *= sign_factor # change sign of Jij + jijs_expanded[ii, 6:9] *= sign_factor # change sign of Dij + + # now collect the outputs in AiiDA Array objects + jij_data = ArrayData() + # jij_data.set_array('Jij_shells', jijs_shells_z) + # if dmimode: + # # also add x and y shells output + # jij_data.set_array('Jij_shells_x', jijs_shells_x) + # jij_data.set_array('Jij_shells_y', jijs_shells_y) + jij_data.set_array('Jij_expanded', jijs_expanded) + jij_data.set_array('positions_expanded', positions_expanded) + # add description to extras + jij_data.extras['array_descriptions'] = { + # 'Jij_shells': """Jij output in the shells that KKR found. + # The format differs for the type of calculation: + # - isotropic exchange constants only (old solver). The columns refer to: + # * [0] |Rij| in units of the lattice constant + # * [1] Jij in Ryd + # * [2] shell-number + # * [3] atom type of atom j + # - full exchange tensor (new solver). The columns refer to: + # * [0] ∣Rij∣ in units of the lattice constant + # * [1] Jij (isotropic part) in Ryd + # * [2] Dij (anti-symmetric DMI part) in Ryd + # * [3] Sij (diagonal traceless part) in Ryd + # * [4] Aij (off-diagonal symmetric part) in Ryd + # * [5-7] Rj - Ri (3 component vector) + # * [8] atom type of atom j + + # if the full exchange tensor is calculated the Jij_shells_x and Jij_shells_y also exist""", + 'Jij_expanded': + 'i, j, da, db, dc, Jij (meV) [, Dij vector (x, y, z in meV), full Jij tensor (xx, xy, xz, yx, yy, yz, zx, zy, zz in meV)]', + 'positions_expanded': 'x, y, z (Ang.) positions of all pairs in Jij_expanded', + } + + # add extras to generated structure for quick access + struc_jij_sites.extras['mappings_ij'] = mappings_back + struc_jij_sites.extras['uuid_struc_orig'] = structure.uuid + struc_jij_sites.extras['uuid_jij_data'] = jij_data.uuid + + # return dict (link_label: node) + return {'jij_data': jij_data, 'structure_jij_sites': struc_jij_sites} + + +def parse_jij_calc(jij_calc, jij_calc_x=None, jij_calc_y=None, verbose=False): + """ + Parse outcome of Jij calculation and return jij_data and structure_jij_sites. + Calculate not only Jij but also Dij vector if 3 calculations for m ||z, m||x and m||y are given. + + :param jij_calc: KkrCalculation with Jij inputs for m || z + :param jij_calc_x: optional, KkrCalculation with Jij inputs for m || x + :param jij_calc_y: optional, KkrCalculation with Jij inputs for m || y + :param verbose: boolean that activate verbose writeout during parsing + + :return jij_data: ArrayData that contains the expanded Jij's (see 'array_descriptions' extra of the node for more details) + :return struc_jij_sites: the reduced structure that contains onlt the atoms which have Jij couplings + (comes from the input of the KkrCalculation). The mappings to the original structure (and their i,j indices) is given as an extra. + """ + + if not jij_calc.is_finished_ok: + raise ValueError('Jij (z) calculation not finished ok (yet)!') + if jij_calc_x is not None and not jij_calc_x.is_finished_ok: + raise ValueError('Jij (x) calculation not finished ok (yet)!') + if jij_calc_y is not None and not jij_calc_y.is_finished_ok: + raise ValueError('Jij (y) calculation not finished ok (yet)!') + + if 'parsed_jij_data' in jij_calc.extras: + # load from extras + if verbose: + print('found parsed jij data, load it now') + jij_data_uuid = jij_calc.extras['parsed_jij_data']['jij_data_uuid'] + struc_jij_sites_uuid = jij_calc.extras['parsed_jij_data']['struc_jij_sites_uuid'] + jij_data = load_node(jij_data_uuid) + structure_jij_sites = load_node(struc_jij_sites_uuid) + else: + # no old data found, do parsing here + # the following calcfunction keeps the data provenance by taking the retrieved folders as inputs + if jij_calc_x is None: + jij_calc_x_retrieved = None + else: + jij_calc_x_retrieved = jij_calc_x.outputs.retrieved + if jij_calc_y is None: + jij_calc_y_retrieved = None + else: + jij_calc_y_retrieved = jij_calc_y.outputs.retrieved + jij_parsed = parse_jij_calc_data( + jij_calc.outputs.retrieved, jij_calc_x_retrieved, jij_calc_y_retrieved, verbose=Bool(verbose) + ) + jij_data, structure_jij_sites = jij_parsed['jij_data'], jij_parsed['structure_jij_sites'] + + # then save as extra to be able to load it in the next run + if verbose: + print('finished parsing jij data, save it as extra now') + jij_calc.set_extra( + 'parsed_jij_data', { + 'jij_data_uuid': jij_data.uuid, + 'struc_jij_sites_uuid': structure_jij_sites.uuid + } + ) + + # finally return result + return jij_data, structure_jij_sites + + +def load_or_submit_cpa_jijcalcs( + scf_cpa_wf, + JIJSITEI=None, + JIJSITEJ=None, + tempr=400., + rclustz=0.9, + kmesh=[100, 100, 100], + NATOMIMPD=500, + NSHELD=2000, + JIJRAD=2.0, + options=None, + uuid_x=None, + uuid_y=None, + uuid_z=None +): + """Load or submit Jij calculation for CPA""" + + scf_remote = scf_cpa_wf.outputs.last_RemoteData + last_calc = scf_remote.get_incoming(node_class=KkrCalculation).first().node + + builder = last_calc.get_builder_restart() + builder.parent_folder = scf_remote + + # set Jij parameters + para_Jij = {k: v for k, v in scf_cpa_wf.outputs.last_InputParameters.get_dict().items() if v} + para_Jij['TEMPR'] = tempr # slightly reduce temperature + para_Jij['RCLUSTZ'] = rclustz # increase cluster radius + para_Jij['BZDIVIDE'] = kmesh # increase k-points + para_Jij['NSTEPS'] = 1 # one-shot + para_Jij['NATOMIMPD'] = NATOMIMPD # array dimension + para_Jij['NSHELD'] = NSHELD # array dimension + para_Jij['KPOIBZ'] = np.product(kmesh) # array dimension + # add 'XCPL' runopt to list of runopts (activates Jij calculation) + runopts = para_Jij.get('RUNOPT', []) + runopts.append('XCPL ') + para_Jij['RUNOPT'] = runopts + # set Jij parameters + # i and j index for Jij calculation in internal units + # uses site index (i.e. needs to be <=10) + if JIJSITEI is not None: + para_Jij['JIJSITEI'] = JIJSITEI + if JIJSITEJ is None: + JIJSITEJ = JIJSITEI + if JIJSITEJ is not None: + para_Jij['JIJSITEJ'] = JIJSITEJ + para_Jij['JIJRAD'] = JIJRAD # radius in lattice constants up to which the Jijs are calculated + + builder.parameters = Dict(dict=para_Jij) + + # starting angles for 3 directions, needed to extract full Jij tensor + + Nsites = len(get_sites(find_parent_structure(last_calc))) + + init_angles_x = Dict( + dict={ + 'fix_dir': [True for i in range(Nsites)], + 'theta': [90.0 for i in range(Nsites)], + 'phi': [0.0 for i in range(Nsites)], + } + ) + + init_angles_y = Dict( + dict={ + 'fix_dir': [True for i in range(Nsites)], + 'theta': [90.0 for i in range(Nsites)], + 'phi': [90.0 for i in range(Nsites)], + } + ) + + init_angles_z = Dict( + dict={ + 'fix_dir': [True for i in range(Nsites)], + 'theta': [0.0 for i in range(Nsites)], + 'phi': [0.0 for i in range(Nsites)], + } + ) + + # submit m||z calculation + builder.initial_noco_angles = init_angles_z + builder.metadata.label = 'Jij_z' + if options is not None: + builder.metadata.options = options + if uuid_z is not None: + calc_Jij_z = load_node(uuid_z) + else: + calc_Jij_z = submit(builder) + print('submitted z', calc_Jij_z) + + builder.initial_noco_angles = init_angles_y + builder.metadata.label = 'Jij_y' + if options is not None: + builder.metadata.options = options + if uuid_y is not None: + calc_Jij_y = load_node(uuid_y) + else: + calc_Jij_y = submit(builder) + print('submitted y', calc_Jij_y) + + builder.initial_noco_angles = init_angles_x + builder.metadata.label = 'Jij_x' + if options is not None: + builder.metadata.options = options + if uuid_x is not None: + calc_Jij_x = load_node(uuid_x) + else: + calc_Jij_x = submit(builder) + print('submitted x', calc_Jij_x) + + return calc_Jij_x, calc_Jij_y, calc_Jij_z diff --git a/aiida_kkr/workflows/__init__.py b/aiida_kkr/workflows/__init__.py index b7004378..0d8659b7 100644 --- a/aiida_kkr/workflows/__init__.py +++ b/aiida_kkr/workflows/__init__.py @@ -14,3 +14,4 @@ from .kkr_imp_dos import kkr_imp_dos_wc from ._combine_imps import combine_imps_wc from ._decimation import kkr_decimation_wc +from .jijs import kkr_jij_wc diff --git a/aiida_kkr/workflows/jijs.py b/aiida_kkr/workflows/jijs.py new file mode 100644 index 00000000..73cdcdff --- /dev/null +++ b/aiida_kkr/workflows/jijs.py @@ -0,0 +1,465 @@ +#!/usr/bin/env python +# coding: utf-8 +""" +This module contains the workflow that can be used to calculate the exchange coupling constants +""" + +import numpy as np +from aiida.engine import WorkChain, ToContext, calcfunction +from aiida.orm import Dict, RemoteData, StructureData, ArrayData, CalcJobNode, Code +from aiida_kkr.tools.find_parent import get_calc_from_remote, find_parent_structure +from aiida_kkr.tools.jij_tools import parse_jij_calc, get_sites +from aiida_kkr.tools.common_workfunctions import test_and_get_codenode, get_parent_paranode, update_params_wf, get_inputs_kkr +from aiida_kkr.calculations.kkr import KkrCalculation +from aiida_kkr.tools.save_output_nodes import create_out_dict_node +from masci_tools.io.common_functions import get_alat_from_bravais +from masci_tools.io.kkr_params import kkrparams +from masci_tools.util.constants import BOHR_A + +__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' +__contributors__ = (u'Philipp Rüßmann') + + +class kkr_jij_wc(WorkChain): + """ + Workchain for calculation of exchange coupling constants Jij and Dij if parent calculation used the SOC solver. + + inputs:: + + :param wf_parameters: optional Dict node of workchain specifications, contains settings like Jij radius cutoff, + selection of sites for i and j and numerical cutoffs. None values in the accuracy sub-dict + means that values from parent calculation are coptied. + :param remote_data: mandatory RemoteData node of parent (i.e. converged) KkrCalculation + :param kkr: optional Code for KKRhost executable (if not given the same as in the parent calculation is used) + :param options: optional Dict computer options like scheduler command or parallelization + + returns:: + + :return jij_data: ArrayData with the arrays 'Jij_expanded' (Table of all Jij and Dij pairs) and 'positions_expanded' (positions of all ij pairs) + :return structure_jij_sites: StructureData + """ + + _wf_version = __version__ + _wf_default = { + 'jijrad_ang': 5.0, # default cutoff radius + 'jijsite_i': None, # use all sites by default + 'jijsite_j': None, # use all sites by default + 'accuracy': { # accuracy settings, typically set to larger values than in scf run + 'NATOMIMPD': 500, + 'NSHELD': 2000, + 'TEMPR': None, + 'RCLUSTZ': None, + 'kmesh': None, + }, + } + _options_default = { + 'max_wallclock_seconds': 36000, + 'resources': { + 'num_machines': 1 + }, + 'withmpi': True, + 'queue_name': '' + } + + @classmethod + def get_wf_defaults(self, silent=False): + """ + Return the default values of the workflow parameters (wf_parameters input node) + """ + if not silent: + print(f'Version of the kkr_jij_wc workflow: {self._wf_version}') + return self._wf_default.copy() + + @classmethod + def define(cls, spec): + """ + Layout of the workflow, defines the input nodes and the outline of the workchain + """ + super(kkr_jij_wc, cls).define(spec) + + # here inputs are defined + spec.input( + 'wf_parameters', + valid_type=Dict, + required=False, + default=lambda: Dict(dict=cls._wf_default), + help='Parameters of the bandstructure workflow (see output of kkr_bs_wc.get_wf_default() for more details).' + ) + + spec.input( + 'options', + valid_type=Dict, + required=False, + default=lambda: Dict(dict=cls._options_default), + help= + 'Computer options (walltime etc.) passed onto KkrCalculation, fall back to settings from parent calculation if not given' + ) + + spec.input( + 'remote_data', + valid_type=RemoteData, + required=True, + help='Parent folder of previously converged KkrCalculation' + ) + + spec.input('kkr', valid_type=Code, required=True, help='KKRhost code, needed to run the Jij KkrCalculation') + + 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) + spec.output('jij_data', valid_type=ArrayData, required=True) + spec.output('structure_jij_sites', valid_type=StructureData, required=True) + + # Here outlines are being specified + spec.outline( + # For initialiging workflow + cls.start, + cls.validate_input, + cls.set_jij_params, + cls.submit_Jij_calcs, + cls.return_results + ) + # definition of exit code in case something goes wrong in this workflow + spec.exit_code( + 160, 'ERROR_KKRCODE_NOT_CORRECT', 'The code you provided for kkr does not use the plugin kkr.kkr' + ) + spec.exit_code(161, 'ERROR_INVALID_PARENT', 'Parent calculation is not valid') + spec.exit_code(162, 'ERROR_CALC_FAILED', 'KKR Band Structure calculation failed') + spec.exit_code(163, 'ERROR_PARSING_FAILED', 'Parsing of Jij calculations failed') + + def start(self): + """ + set up context of the workflow + """ + self.report(f'INFO: started KKR Jij workflow version {self._wf_version}') + if 'wf_parameters' in self.inputs: + wf_dict = self.inputs.wf_parameters.get_dict() + else: + wf_dict = {} + # add missing default valuesi + for key, val in self._wf_default.items(): + if ((key not in wf_dict.keys()) and (key.swapcase() not in wf_dict.keys()) and (val is not None)): + self.report(f'INFO: Using default wf parameter {key}: {val}') + wf_dict[key] = val + + if 'options' in self.inputs: + options_dict = self.inputs.options.get_dict() + else: + options_dict = self._options_default + self.ctx.options = options_dict + + 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( + 'max_wallclock_seconds', self._options_default['max_wallclock_seconds'] + ) + self.ctx.queue = options_dict.get('queue_name', self._options_default['queue_name']) + self.ctx.custom_scheduler_commands = options_dict.get('custom_scheduler_commands', '') + self.ctx.wf_dict = wf_dict + self.report( + 'INFO: use the following parameter:\n' + 'withmpi: {}\n' + 'Resources: {}\n' + 'Walltime (s): {}\n' + 'queue name: {}\n' + 'scheduler command: {}\n' + 'Workflow parameters: {}\n'.format( + self.ctx.withmpi, self.ctx.resources, self.ctx.max_wallclock_seconds, self.ctx.queue, + self.ctx.custom_scheduler_commands, wf_dict + ) + ) + + def validate_input(self): + """ + validate inputs + """ + + # save parent calculation + input_remote = self.inputs.remote_data + parents = input_remote.get_incoming(node_class=CalcJobNode).all() + if len(parents) != 1: + # check if parent is unique + return self.exit_codes.ERROR_INVALID_PARENT # pylint: disable=no-member + self.ctx.parent_calc = get_calc_from_remote(input_remote) + + # validate for kkrcode + try: + test_and_get_codenode(self.inputs.kkr, 'kkr.kkr', use_exceptions=True) + except ValueError: + return self.exit_codes.ERROR_KKRCODE_NOT_CORRECT # pylint: disable=no-member + + # get parent input parameter + self.ctx.input_params_KKR = get_parent_paranode(self.inputs.remote_data) + + def set_jij_params(self): + """ + set kkr parameters for the Jij calculation + """ + + # input parameters from parent + 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) + + # set Jij parameters + para_jij, runopts = self._get_para_jij(params) + updatenode = Dict(dict=para_jij.get_dict()) + updatenode.label = 'Jij params' + paranode_jij = update_params_wf(params, updatenode) + self.ctx.jij_params = paranode_jij + + # find out if we have a calculation with or without SOC (then no DMI is calculated) + self.ctx.noSOC = True + if 'NEWSOSOL' in runopts or para_jij.get_value(''): + self.ctx.noSOC = False + + def submit_Jij_calcs(self): + """ + submit the KkrCalcultion with the Jij settings + """ + + # get inputs for band structure calculation + inputs = get_inputs_kkr( + self.inputs.kkr, + self.inputs.remote_data, + self.ctx.options, + label='Jij_calc', + description='', + parameters=self.ctx.jij_params, + serial=(not self.ctx.withmpi) + ) + + # noSOC, only m||z + if self.ctx.noSOC: + inputs.metadata.label = 'jij_calc_z' # pylint: disable=no-member + jij_calc_z = self.submit(KkrCalculation, **inputs) + self.ctx.jij_calc_z = jij_calc_z + self.ctx.jij_calc_x = None + self.ctx.jij_calc_y = None + else: + # create nonco angles for the three calculations + nonco_angles = _make_nonco_angles(parent_remote=self.inputs.remote_data) + init_angles_x, init_angles_y, init_angles_z = nonco_angles['init_angles_x'], nonco_angles[ + 'init_angles_y'], nonco_angles['init_angles_z'] + + # submit m||z calculation + inputs.initial_noco_angles = init_angles_z + inputs.metadata.label = 'jij_calc_z' # pylint: disable=no-member + jij_calc_z = self.submit(KkrCalculation, **inputs) + self.ctx.jij_calc_z = jij_calc_z + + # submit m||x calculation + inputs.initial_noco_angles = init_angles_x + inputs.metadata.label = 'jij_calc_x' # pylint: disable=no-member + jij_calc_x = self.submit(KkrCalculation, **inputs) + self.ctx.jij_calc_x = jij_calc_x + + # submit m||y calculation + inputs.initial_noco_angles = init_angles_y + inputs.metadata.label = 'jij_calc_y' # pylint: disable=no-member + jij_calc_y = self.submit(KkrCalculation, **inputs) + self.ctx.jij_calc_y = jij_calc_y + + # add to context (needed to tell aiida to wait for processes to finish) + futures = {'jij_calc_z': self.ctx.jij_calc_z} + if self.ctx.jij_calc_x is not None: + futures['jij_calc_x'] = self.ctx.jij_calc_x + if self.ctx.jij_calc_y is not None: + futures['jij_calc_y'] = self.ctx.jij_calc_y + return ToContext(**futures) + + def return_results(self): + """ + Collect results, parse Jij output and link output nodes to workflow node + """ + + # check if calculations finished ok + success = True + if self.ctx.jij_calc_z is not None and not self.ctx.jij_calc_z.is_finished_ok: + success = False + if self.ctx.jij_calc_x is not None and not self.ctx.jij_calc_x.is_finished_ok: + success = False + if self.ctx.jij_calc_y is not None and not self.ctx.jij_calc_y.is_finished_ok: + success = False + if not success: + self.ctx.successful = False + error = f'ERROR Jij calculation failed somehow it is in state {self.ctx.jij_calc_z.process_state}' + if self.ctx.jij_calc_x is not None: + error += f'; {self.ctx.jij_calc_x.process_state} (x)' + if self.ctx.jij_calc_y is not None: + error += f'; {self.ctx.jij_calc_y.process_state} (y)' + self.report(error) + return self.exit_codes.ERROR_CALC_FAILED # pylint: disable=no-member + + # now parse calculation output + try: + jij_data, structure_jij_sites = parse_jij_calc( + self.ctx.jij_calc_z, jij_calc_x=self.ctx.jij_calc_x, jij_calc_y=self.ctx.jij_calc_y, verbose=False + ) + except Exception as err: + self.report(f'Caught error when trying to parse Jij output:{err}') + return self.exit_codes.ERROR_PARSING_FAILED # pylint: disable=no-member + + # collect output nodes + outdict = {'jij_data': jij_data, 'structure_jij_sites': structure_jij_sites} + + # create dict to store results of workflow output + outputnode_dict = {} + outputnode_dict['workflow_name'] = self.__class__.__name__ + outputnode_dict['workflow_version'] = self._wf_version + outputnode_dict['successful'] = success + + # create output node with data-provenance + outputnode = Dict(dict=outputnode_dict) + + # link to the output nodes + link_nodes = outdict.copy() + + outdict['results_wf'] = create_out_dict_node(outputnode, **link_nodes) + + # create links to output nodes + for link_name, node in outdict.items(): + self.out(link_name, node) + + self.report('INFO: done with Jij workflow!') + + # Helper functions + + def _get_para_jij(self, params): + """ + Set the Jij parameters from the input. + Returns a kkrparams instance with the set values + """ + + # get input parameters + input_dict = params.get_dict() + para_jij = kkrparams(**input_dict) + + # set Jij parameters + # add 'XCPL' runopt to list of runopts (activates Jij calculation) + runopts = input_dict.get('RUNOPT') + if runopts is None: + runopts = [] + runopts.append('XCPL ') + para_jij.set_value('RUNOPT', runopts) + para_jij.set_value('NSTEPS', 1) # one-shot run + + # accuracy settings + tempr = self.ctx.wf_dict['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') + if rclustz is not None: + para_jij.set_value('RCLUSTZ', rclustz) # increase cluster radius + + kmesh = self.ctx.wf_dict['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') + if NATOMIMPD is not None: + para_jij.set_value('NATOMIMPD', NATOMIMPD) # array dimension + + NSHELD = self.ctx.wf_dict['accuracy'].get('NSHELD') + if NSHELD is not None: + para_jij.set_value('NSHELD', NSHELD) # array dimension + + # Jij settings + jijrad = self._get_jijrad() + if jijrad is not None: + para_jij.set_value('JIJRAD', jijrad) # radius in lattice constants up to which the Jijs are calculated + + # 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'] + if JIJSITEI is not None: + if JIJSITEJ is None: + JIJSITEJ = JIJSITEI + para_jij.set_value('JIJSITEI', [len(JIJSITEI)] + [i + 1 for i in JIJSITEI]) + if JIJSITEJ is not None: + para_jij.set_value('JIJSITEJ', [len(JIJSITEJ)] + [i + 1 for i in JIJSITEJ]) + + return para_jij, runopts + + def _get_jijrad(self): + """ + get Jij radius convert from Ang to internal alat units + """ + + # Jij radius in Ang. + jijrad_ang = self.ctx.wf_dict['jijrad_ang'] + + # find structure from calculation + struc, _ = find_parent_structure(self.ctx.parent_calc) + self.ctx.structure = struc + + # get alat from structure + alat_ang = get_alat_from_bravais(np.array(struc.cell), struc.pbc[2]) + + # maybe use value provided in input instead + para = {k.lower(): v for k, v in self.ctx.parent_calc.inputs.parameters.get_dict().items()} + if 'use_alat_input' in para: + alat_ang = para.get('ALATBASIS') * BOHR_A + + # now have Jij radius in alat units + jijrad = jijrad_ang / alat_ang + + return jijrad + + +@calcfunction +def _make_nonco_angles(parent_remote): + """ + Create nonco angles for the 3 directions (x y, z) + """ + # find structure to count number of sites + structure = find_parent_structure(parent_remote)[0] + Nsites = len(get_sites(structure)) + + # create nonco angles for m||z + init_angles_z = Dict( + dict={ + 'fix_dir': [True for i in range(Nsites)], + 'theta': [0.0 for i in range(Nsites)], + 'phi': [0.0 for i in range(Nsites)], + } + ) + + # create nonco angles for m||x + init_angles_x = Dict( + dict={ + 'fix_dir': [True for i in range(Nsites)], + 'theta': [90.0 for i in range(Nsites)], + 'phi': [0.0 for i in range(Nsites)], + } + ) + + # create nonco angles for m||y + init_angles_y = Dict( + dict={ + 'fix_dir': [True for i in range(Nsites)], + 'theta': [90.0 for i in range(Nsites)], + 'phi': [90.0 for i in range(Nsites)], + } + ) + + return {'init_angles_x': init_angles_x, 'init_angles_y': init_angles_y, 'init_angles_z': init_angles_z} diff --git a/aiida_kkr/workflows/kkr_scf.py b/aiida_kkr/workflows/kkr_scf.py index e3f29943..602a4c47 100644 --- a/aiida_kkr/workflows/kkr_scf.py +++ b/aiida_kkr/workflows/kkr_scf.py @@ -426,7 +426,6 @@ def start(self): self.ctx.calcs = [] self.ctx.abort = False # flags used internally to check whether the individual steps were successful - self.ctx.kkr_converged = False self.ctx.dos_ok = False self.ctx.voro_step_success = False self.ctx.kkr_step_success = False @@ -1283,8 +1282,6 @@ def inspect_kkr(self): self.report(f'INFO: last_remote: {self.ctx.last_remote}') if self.ctx.kkr_step_success and found_last_calc_output: - # check convergence - self.ctx.kkr_converged = last_calc_output['convergence_group']['calculation_converged'] # check rms, compare spin and charge values and take bigger one rms_charge = last_calc_output['convergence_group']['rms'] # returning 0 if not found allows to reuse older verisons (e.g. in caching) @@ -1312,7 +1309,14 @@ def inspect_kkr(self): if self.ctx.kkr_step_success and self.convergence_on_track(): self.ctx.rms_all_steps += rms_all_iter_last_calc self.ctx.neutr_all_steps += neutr_all_iter_last_calc + + # check if calculation converged + self.ctx.kkr_converged = True + if rms_max > self.ctx.convergence_criterion: + self.ctx.kkr_converged = False + else: + # if last step did not succeed we know the calculation did not converge self.ctx.kkr_converged = False self.report(f'INFO: kkr_converged: {self.ctx.kkr_converged}') diff --git a/docs/source/module_guide/workflows.rst b/docs/source/module_guide/workflows.rst index e9e11362..23e5a817 100644 --- a/docs/source/module_guide/workflows.rst +++ b/docs/source/module_guide/workflows.rst @@ -48,6 +48,14 @@ Find Green Function writeout for KKRimp :private-members: :special-members: +KKRhost Jij calculation +----------------------- +.. automodule:: aiida_kkr.workflows.jijs + :members: + :private-members: + :special-members: + + KKRimp self-consistency ----------------------- .. automodule:: aiida_kkr.workflows.kkr_imp_sub diff --git a/docs/source/user_guide/workflows.rst b/docs/source/user_guide/workflows.rst index 8b78d8c8..4eb5571c 100644 --- a/docs/source/user_guide/workflows.rst +++ b/docs/source/user_guide/workflows.rst @@ -454,6 +454,71 @@ Finally we run the workflow:: from aiida.work import run run(kkr_flex_wc, label='test_gf_writeout', description='My test KKRflex calculation.', kkr=kkrcode, remote_data=kkr_remote_folder, options=options, wf_parameters=wf_params) + + +Exchange coupling constants ++++++++++++++++++++++++++++ + +Calculation of the exchange coupling constants (Jij's and Dij's) can be done with the ``kkr_jij_wc`` workchain starting from the remote folder of a parent calculation. + +.. note:: + Use ``kkr_jij_wc.get_wf_defaults()`` to get the default values for the ``wf_parameters`` input. + +Inputs: + + * ``wf_parameters`` (Dict, optional): Workchain settings where the Jij radius in Angstroem units is given (defaults to 5 Ang.). With ``jijsite_i`` and ``jijsite_j`` one can specify to calculate the Jij's only for some pairs i,j. + + * ``options`` (Dict, optional): Computer options (scheduler command, parallelization, walltime etc.) + + * ``remote_data`` (RemoteData, mandatory): Parent folder of a converged KkrCalculation. + + * ``kkr`` (Code, mandatory): KKRhost code (i.e. using ``kkr.kkr`` plugin). + + * ``params_kkr_overwrite`` (Dict, optional): Optional set of KKR parameters that overwrite the settings extracted from the parent calculation. + +Returns nodes: + * ``jij_data`` (ArrayData): Jij data with the arrays ``Jij_expanded`` that contains all (i, j, da, db, dc, Jij [, Dij]) pairs and ``positions_expanded`` that contains the corresponding positions (i.e. the offset of j vs i). + + * ``structure_jij_sites`` (StructureData): Structure with the Jij sites that match the mapping in ``jij_data`` + +Example Usage: +^^^^^^^^^^^^^^ + +To start the Band Structure calculation the steps: + +:: + + from aiida_kkr.workflows import kkr_jij_wc + + # converged KKR calculation + kkr_calc_converged = load_node() + + # create process builder + builder = kkr_jij_wc.get_builder() + + builder.parent_folder = kkr_calc_converged.outputs.remote_folder + + builder.kkr = parent.inputs.code + builder.options = Dict(dict={...}) # settings for the computer that we use + + wfd = kkr_jij_wc.get_wf_defaults() + wfd['jijrad_ang'] = 5.0 # set at least the Jij cutoff radius in Ang units + builder.wf_parameters = Dict(dict=wfd) + + # maybe overwrite some input parameters + # here we switch on the SOC mode starting from a no SOC calculation + builder.params_kkr_overwrite = Dict(dict={ + 'NPAN_LOG': 5, + 'NPAN_EQ': 15, + 'NCHEB': 12, + 'R_LOG': 0.6, + 'USE_CHEBYCHEV_SOLVER': True, + 'SET_CHEBY_NOSOC': True + }) + + # submit calculation + jij_wf = submit(builder) + KKR impurity self consistency @@ -723,7 +788,7 @@ Equation of states Workflow: ``aiida_kkr.workflows.eos`` -.. warning:: Not implemented yet! +.. warning:: Not documented yet! Check KKR parameter convergence @@ -753,6 +818,3 @@ the ground state or not. Then the unit cell could be doubled to compute the antiferromagnetic state. In case of noncollinear magnetism the full Jij tensor should be analyzed. - - - diff --git a/setup.json b/setup.json index 9aabb790..9be8e234 100644 --- a/setup.json +++ b/setup.json @@ -86,7 +86,8 @@ "kkr.imp = aiida_kkr.workflows.kkr_imp:kkr_imp_wc", "kkr.imp_sub = aiida_kkr.workflows.kkr_imp_sub:kkr_imp_sub_wc", "kkr.imp_dos = aiida_kkr.workflows.kkr_imp_dos:kkr_imp_dos_wc", - "kkr.decimation = aiida_kkr.workflows._decimation:kkr_decimation_wc" + "kkr.decimation = aiida_kkr.workflows._decimation:kkr_decimation_wc", + "kkr.jij = aiida_kkr.workflows.jijs:kkr_jij_wc" ], "console_scripts": [ "aiida-kkr = aiida_kkr.cmdline:cmd_root" diff --git a/tests/data_dir/kkr_jij_wc-nodes-7316580c3630215166fd2681063d792e.tar.gz b/tests/data_dir/kkr_jij_wc-nodes-7316580c3630215166fd2681063d792e.tar.gz new file mode 100644 index 00000000..5a363a12 Binary files /dev/null and b/tests/data_dir/kkr_jij_wc-nodes-7316580c3630215166fd2681063d792e.tar.gz differ diff --git a/tests/data_dir/kkr_jij_wc-nodes-b8b421acfdbc243b6c9a1728b88e7612.tar.gz b/tests/data_dir/kkr_jij_wc-nodes-b8b421acfdbc243b6c9a1728b88e7612.tar.gz new file mode 100644 index 00000000..f2bcb087 Binary files /dev/null and b/tests/data_dir/kkr_jij_wc-nodes-b8b421acfdbc243b6c9a1728b88e7612.tar.gz differ diff --git a/tests/files/jij_input.aiida b/tests/files/jij_input.aiida new file mode 100644 index 00000000..e6142540 Binary files /dev/null and b/tests/files/jij_input.aiida differ diff --git a/tests/run_all.sh b/tests/run_all.sh index a4e8d1f6..620ddba5 100755 --- a/tests/run_all.sh +++ b/tests/run_all.sh @@ -115,6 +115,7 @@ elif [[ ! -z "$GITHUB_SUITE" ]]; then ./workflows/test_dos_wc.py \ ./workflows/test_gf_writeout_wc.py \ ./workflows/test_scf_wc_simple.py \ + ./workflows/test_jij_wc.py \ ./workflows/test_eos.py \ ./workflows/test_kkrimp_sub_wc.py \ $addopt diff --git a/tests/test_entrypoints.py b/tests/test_entrypoints.py index 441d0ca8..7608b5f7 100644 --- a/tests/test_entrypoints.py +++ b/tests/test_entrypoints.py @@ -179,6 +179,13 @@ def test_decimation_entry_point(self): wf = WorkflowFactory('kkr.decimation') assert wf == kkr_decimation_wc + def test_jij_entry_point(self): + from aiida_kkr.workflows.jijs import kkr_jij_wc + from aiida.plugins import WorkflowFactory + + wf = WorkflowFactory('kkr.jij') + assert wf == kkr_jij_wc + def test_kkrimp_workchain_entry_point(self): from aiida_kkr.workflows.kkr_imp import kkr_imp_wc from aiida.plugins import WorkflowFactory diff --git a/tests/tools/test_find_parent.py b/tests/tools/test_find_parent.py new file mode 100755 index 00000000..d6c3c6b5 --- /dev/null +++ b/tests/tools/test_find_parent.py @@ -0,0 +1,39 @@ +#!/usr/bin/env python + +import pytest +from aiida.orm import load_node +from aiida_kkr.tools.find_parent import find_parent_structure, get_calc_from_remote +from ..conftest import import_with_migration + + +def test_find_parent_structure(): + """ + find parent structure and Voronoi parent calculation from KKR calculation + """ + + # load necessary files from db_dump files + import_with_migration('files/db_dump_kkrcalc.tar.gz') + kkr_calc = load_node('3058bd6c-de0b-400e-aff5-2331a5f5d566') + + # now find voronoi parent and structure + struc, voro_parent = find_parent_structure(kkr_calc) + + # check result + assert struc.uuid == 'e51ee6a1-bd27-4901-9612-7bac256bf117' + assert voro_parent.uuid == '559b9d9b-3525-402e-9b24-ecd8b801853c' + + +def test_get_calc_from_remote(): + """ + find parent calc from remote + """ + + # load necessary files from db_dump files + import_with_migration('files/db_dump_kkrcalc.tar.gz') + kkr_calc = load_node('3058bd6c-de0b-400e-aff5-2331a5f5d566') + + # now find voronoi parent and structure + parent_from_remote = get_calc_from_remote(kkr_calc.outputs.remote_folder) + + # check result + assert parent_from_remote.uuid == kkr_calc.uuid diff --git a/tests/workflows/test_jij_wc.py b/tests/workflows/test_jij_wc.py new file mode 100644 index 00000000..7b2ff6f3 --- /dev/null +++ b/tests/workflows/test_jij_wc.py @@ -0,0 +1,143 @@ +#!/usr/bin/env python +# coding: utf-8 +""" +Tests for Jij workflow +""" + +import pytest +from ..dbsetup import * +from aiida_testing.export_cache._fixtures import run_with_cache, export_cache, load_cache, hash_code_by_entrypoint +from aiida.manage.tests.pytest_fixtures import ( + aiida_local_code_factory, aiida_localhost, temp_dir, aiida_profile, clear_database, clear_database_after_test, + clear_database_before_test +) +from ..conftest import kkrhost_local_code, data_dir, import_with_migration + + +@pytest.mark.timeout(240, method='thread') +def test_jij(clear_database_before_test, kkrhost_local_code, run_with_cache, ndarrays_regression): + """ + Jij test with SOC + """ + + from aiida.orm import load_node, Dict + from aiida_kkr.workflows.jijs import kkr_jij_wc + + # for runing in local computer + options = Dict( + dict={ + 'queue_name': queuename, + 'resources': { + 'num_machines': 1 + }, + 'max_wallclock_seconds': 5 * 60, + 'withmpi': False, + 'custom_scheduler_commands': '' + } + ) + + # import calculation which is used as parent calculation + import_with_migration('files/jij_input.aiida') + parent = load_node('d3e62972-2f38-470b-8bbf-300c21847a6a') + + # now set up process builder + builder = kkr_jij_wc.get_builder() + wfd = kkr_jij_wc.get_wf_defaults() + wfd['jijrad_ang'] = 5.0 + builder.wf_parameters = Dict(dict=wfd) + builder.remote_data = parent.outputs.remote_folder + builder.kkr = kkrhost_local_code + builder.options = options + + # make test run faster + builder.params_kkr_overwrite = Dict(dict={ + 'LMAX': 2, + 'BZDIVIDE': [10, 10, 10], + 'NPT1': 3, + 'NPT2': 10, + 'NPT3': 3, + }) + + # run the calculation using cached data is available + out, _ = run_with_cache(builder, data_dir=data_dir) + + # check results + print('check results', out) + n = out['results_wf'] + assert n.get_dict().get('successful') + # check Jij data arrays + check_dict = { + 'Jij_expanded': out['jij_data'].get_array('Jij_expanded'), + 'positions_expanded': out['jij_data'].get_array('positions_expanded') + } + print(check_dict['Jij_expanded'][:10]) + ndarrays_regression.check(check_dict) + + +@pytest.mark.timeout(240, method='thread') +def test_jij_soc(clear_database_before_test, kkrhost_local_code, run_with_cache, ndarrays_regression): + """ + Jij test with SOC + """ + + from aiida.orm import load_node, Dict + from aiida_kkr.workflows.jijs import kkr_jij_wc + + # for runing in local computer + options = Dict( + dict={ + 'queue_name': queuename, + 'resources': { + 'num_machines': 1 + }, + 'max_wallclock_seconds': 12 * 60, + 'withmpi': False, + 'custom_scheduler_commands': '' + } + ) + + # import calculation which is used as parent calculation + import_with_migration('files/jij_input.aiida') + parent = load_node('d3e62972-2f38-470b-8bbf-300c21847a6a') + + # now set up process builder + builder = kkr_jij_wc.get_builder() + wfd = kkr_jij_wc.get_wf_defaults() + wfd['jijrad_ang'] = 5.0 + builder.wf_parameters = Dict(dict=wfd) + builder.remote_data = parent.outputs.remote_folder + builder.kkr = kkrhost_local_code + builder.options = options + + # make test run faster + builder.params_kkr_overwrite = Dict( + dict={ + 'RUNOPT': ['NEWSOSOL'], + 'LMAX': 2, + 'BZDIVIDE': [10, 10, 10], + 'NPT1': 3, + 'NPT2': 10, + 'NPT3': 3, + 'NPAN_LOG': 5, + 'NPAN_EQ': 12, + 'NCHEB': 12, + 'R_LOG': 0.6, + '': True, + '': True + } + ) + + # run the calculation using cached data is available + out, _ = run_with_cache(builder, data_dir=data_dir) + + # check results + print('check results', out) + n = out['results_wf'] + assert n.get_dict().get('successful') + # check Jij data arrays + check_dict = { + 'Jij_expanded': out['jij_data'].get_array('Jij_expanded'), + 'positions_expanded': out['jij_data'].get_array('positions_expanded') + } + print(check_dict['Jij_expanded'][:10]) + ndarrays_regression.check(check_dict) diff --git a/tests/workflows/test_jij_wc/test_jij.npz b/tests/workflows/test_jij_wc/test_jij.npz new file mode 100644 index 00000000..a526dac9 Binary files /dev/null and b/tests/workflows/test_jij_wc/test_jij.npz differ diff --git a/tests/workflows/test_jij_wc/test_jij_soc.npz b/tests/workflows/test_jij_wc/test_jij_soc.npz new file mode 100644 index 00000000..6a84528c Binary files /dev/null and b/tests/workflows/test_jij_wc/test_jij_soc.npz differ