new functions in tools_STM for visualization

aiida_kkr/tools/tools_STM_scan.py

@@ -13,13 +13,12 @@
 __copyright__ = (u'Copyright (c), 2023, 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')
+__version__ = '0.1.1'
+__contributors__ = (u'Philipp Rüßmann, Raffaele Aliberti')
 
 ##############################################################################
 # combine impurty clusters
 
-
 def convert_to_imp_cls(host_structure, imp_info):
     """
     convert imp info to rcls form
@@ -228,3 +227,227 @@ def create_combined_potential_node_cf(add_position, host_remote, imp_potential_n
     pot_combined_node = create_combined_potential_node(add_position, host_remote, imp_potential_node)
 
     return pot_combined_node
+
+
+#####################################################################
+# STM pathfinder 
+
+
+def STM_pathfinder(host_structure):
+    from aiida_kkr.tools import find_parent_structure
+    from ase.spacegroup import Spacegroup
+
+    """This function is used to help visualize the scanned positions 
+       and the symmetries that are present in the system            """
+
+    """
+    inputs: 
+    host_struture : RemoteData : The Remote data contains all the information needed to create the path to scan
+    
+    outputs: 
+    struc_info : Dict  : Dictionary containing the structural information of the film
+    matrices   : Array : Array containing the matrices that generate the symmetries of the system
+    
+    """
+
+    def info_creation(structure):
+        from ase.spacegroup import get_spacegroup
+            # List of the Bravais vectors
+        vec_list = structure.cell.tolist()
+
+        # Find the Bravais vectors that are in plane vectors
+        plane_vectors = {'plane_vectors':[], 'space_group':''}
+        for vec in vec_list:
+            # Is this sufficient to find all the in-plane vectors?
+            if vec[2] == 0:
+                plane_vectors['plane_vectors'].append(vec[:2])
+
+        space_symmetry = get_spacegroup(ase_struc)
+        plane_vectors['space_group'] = space_symmetry.no
+
+        return plane_vectors
+
+    def symmetry_finder(struc_info):
+        # Here we get the symmetry operations that are possible
+        symmetry_matrices = Spacegroup(struc_info['space_group'])
+
+        # Reduce the dimensionality, we only want the 2D matrices
+        matrices = []
+        for element in symmetry_matrices.get_rotations():
+            matrices.append(element[:2, :2])
+
+        # Uniqueness of the elements must be preserved
+        unique_matrices = []
+        for matrix in matrices:
+            if not any(np.array_equal(matrix, m) for m in unique_matrices):
+                unique_matrices.append(matrix)
+
+        return unique_matrices
+
+    struc = find_parent_structure(host_structure)
+    # clone the structure since it has already been saved in AiiDA and cannot be modified
+    supp_struc = struc.clone()
+
+    # If the structure is not periodic in every direction we force it to be.
+    supp_struc.pbc = (True, True, True)
+
+    # ASE struc 
+    ase_struc = supp_struc.get_ase()
+
+    # Structural informations are stored here
+    struc_info = info_creation(ase_struc)
+
+    # The structural informations are then used to find the symmetries of the system
+    symm_matrices = symmetry_finder(struc_info)
+
+    return struc_info, symm_matrices 
+
+@engine.calcfunction
+def STM_pathfinder_cf(host_structure):
+    """
+    Calcfunction that gives back the structural information of the film, and the symmetries of the system
+    """
+
+    struc_info, symm_matrices = STM_pathfinder(host_structure)
+
+    return struc_info, symm_matrices
+
+
+###################################################################
+# lattice generation (function of lattice plot)
+
+
+def lattice_generation(x_len, y_len, rot, vec):
+    import math    
+
+    """
+    x_len : int  : value to create points between - x and x.
+    y_len : int  : value to create points between - y and y.
+    rot   : list : list of the rotation matrices given by the symmetry of the system.
+    vec   : list : list containing the two Bravais vectors.
+    """
+
+    # Here we create a grid in made of points whic are the linear combination of the lattice vectors
+    lattice_points = []
+
+    for i in range(-x_len, x_len+1):
+        lattice_points_col = []
+        for j in range(-y_len, y_len+1):
+            p = [i*x + j*y for x,y in zip(vec[0], vec[1])]
+            lattice_points_col.append(p)
+        lattice_points.append(lattice_points_col)
+
+    # Eliminiatio of the symmetrical sites
+    points_to_eliminate = []
+
+    for i in range(-x_len, x_len+1):
+        for j in range(-y_len, y_len+1):
+            if lattice_points[i][j][0] >= 0 and lattice_points[i][j][1] >= 0:
+                for element in rot[1:]:
+                    point = np.dot(element, lattice_points[i][j])
+                    if point[0] >= 0 and point[1] >=0:
+                        continue
+                    else:
+                        points_to_eliminate.append(point)
+
+    points_to_scan = []
+
+    for i in range(-x_len, x_len + 1):
+        for j in range(-y_len, y_len + 1):
+            eliminate = False
+            for elem in points_to_eliminate:
+                # Since there can be some numerical error in the dot product we use the isclose function
+                if (math.isclose(elem[0], lattice_points[i][j][0], abs_tol=1e-4) and math.isclose(elem[1], lattice_points[i][j][1], abs_tol=1e-4)):
+                    eliminate = True
+            if not eliminate:
+                points_to_scan.append(lattice_points[i][j])
+
+    return points_to_eliminate, points_to_scan
+
+
+###############################################################################################
+# lattice plot
+
+
+def lattice_plot(plane_vectors, symm_vec, symm_matrices, grid_length_x, grid_length_y):
+    from aiida_kkr.tools import lattice_generation
+    import matplotlib.pyplot as plt
+    from matplotlib.lines import Line2D
+
+    origin = np.array([[0, 0],[0, 0]])
+    # Generation of the points to plot
+    unused , used = lattice_generation(grid_length_x, grid_length_y, symm_matrices, plane_vectors)
+
+    # Plotting of the points
+    for element in unused:
+        plt.scatter(element[0], element[1], marker='s', s=130, c='#33638DFF')
+
+    for element in used:
+        plt.scatter(element[0], element[1], marker='D', s=130, c='#FDE725FF')       
+
+    # Plot of the crystal symmetry directions, tag must be activated. 
+    if symm_vec:
+        import numpy.linalg as lin
+
+        for element in symm_matrices:
+            eig_val, eig_vec = lin.eig(element)
+
+        for element in eig_vec:
+            plt.quiver(*origin, element[0], element[1],  alpha=1, color='#B8DE29FF', angles='xy', scale_units='xy', scale=1)
+
+    # Plot of the Bravais lattice
+    for element in plane_vectors:
+        plt.quiver(*origin, element[0], element[1], color='#3CBB75FF', angles='xy', scale_units='xy', scale=1)
+
+    legend_elements = [
+    Line2D([0], [0], color='#33638DFF', lw=2, label='Unscanned Sites', marker='s'),
+    Line2D([0], [0], color='#FDE725FF', lw=2, label='Scanned Sites', marker='D'),
+    Line2D([0], [0], color='#3CBB75FF', lw=2, label='Bravais lattice'),
+    ]
+    plt.legend(handles=legend_elements, bbox_to_anchor=(0.75, -0.15))
+
+    plt.title('Lattice plot and symmetry directions')
+    plt.ylabel('y direction')
+    plt.xlabel('x direction')
+    #plt.xticks(np.arange(-grid_length, grid_length, float(plane_vectors[0][0])))
+    #plt.set_cmap(cmap)
+    plt.grid(linestyle='--')    
+    plt.show()
+
+
+##########################################################
+# find linear combination coefficients
+
+
+def find_linear_combination_coefficients(plane_vectors, vectors):
+    from operator import itemgetter
+    """This helper function takes the planar vectors and a list of vectors
+       and return the coefficients in the base of the planar vectors"""
+
+    # Formulate the system of equations Ax = b
+    A = np.vstack((plane_vectors[0], plane_vectors[1])).T
+
+    # Solve the system of equations using least squares method
+    data = []
+    for element in vectors:
+        b = element 
+        # We use the least square mean error procedure to evaulate the units of da and db
+        # lstsq returns: coeff, residue, rank, and singular value
+        # We only need the coefficient.
+        data.append(np.linalg.lstsq(A, b, rcond=None))
+
+    indices = []
+    for element in data:
+        supp = []
+        for elem in element[0]:
+            # Here we round to an integer, this is because of the numerical error 
+            # which is present inside the calculation. 
+            supp.append(round(elem))
+        indices.append(supp)
+
+    # Before returning the indices, we reorder them first from the lowest to the highest valued
+    # on the x axis and then from the lowest to the highest on the y axis. 
+
+    indices = sorted(indices, key=itemgetter(0,1))
+
+    return indices

aiida_kkr/workflows/bs.py

@@ -107,7 +107,7 @@ def define(cls, spec):
             'remote_data',
             valid_type=RemoteData,
             required=True,
-            help='Parent folder of previoously converged KkrCalculation'
+            help='Parent folder of previously converged KkrCalculation'
         )
         spec.input('kkr', valid_type=Code, required=True, help='KKRhost code, needed to run the qdos KkrCalculation')
         spec.input(

aiida_kkr/workflows/kkr_imp_dos.py

@@ -161,7 +161,9 @@ def define(cls, spec):
         spec.output('dos_data_interpol', valid_type=XyData, required=False)
         spec.output('dos_data_lm', valid_type=XyData, required=False)
         spec.output('dos_data_interpol_lm', valid_type=XyData, required=False)
-        spec.output('gf_dos_remote', valid_type=RemoteData, required=False, help='RemoteData node of the computed host GF.')
+        spec.output(
+            'gf_dos_remote', valid_type=RemoteData, required=False, help='RemoteData node of the computed host GF.'
+        )
 
         # Here the structure of the workflow is defined
         spec.outline(
@@ -468,7 +470,7 @@ def run_imp_dos(self):
             gf_writeout_remote = self.inputs.gf_dos_remote
             gf_writeout_calc = gf_writeout_remote.get_incoming(node_class=CalcJobNode).first().node
             self.ctx.pk_flexcalc = gf_writeout_calc.pk
-        
+
         self.ctx.gf_writeout_remote = gf_writeout_remote
 
         options = self.ctx.options_params_dict
@@ -737,18 +739,18 @@ def parse_impdosfiles(folder, natom, nspin, ef, use_lmdos):
                 with folder.open(name0 + '.atom=%0.2i_spin%i.dat' % (iatom, ispin)) as dosfile:
                     tmp = loadtxt(dosfile)
                     dos.append(tmp)
-                
+
                 with folder.open(name0 + '.interpol.atom=%0.2i_spin%i.dat' % (iatom, ispin)) as dosfile:
                     tmp = loadtxt(dosfile)
                     if len(tmp) > 0:
                         dos_int.append(tmp)
-                                     
+
             except:
                 # new file names with 3 digits for atom numbers
                 with folder.open(name0 + '.atom=%0.3i_spin%i.dat' % (iatom, ispin)) as dosfile:
                     tmp = loadtxt(dosfile)
                     dos.append(tmp)
-                
+
                 with folder.open(name0 + '.interpol.atom=%0.3i_spin%i.dat' % (iatom, ispin)) as dosfile:
                     try:
                         tmp = loadtxt(dosfile)
@@ -757,7 +759,6 @@ def parse_impdosfiles(folder, natom, nspin, ef, use_lmdos):
                     # can happen if there are too few points to interpolate on
                     if len(tmp) > 0:
                         dos_int.append(tmp)
-
 
     dos = array(dos)
     if len(dos_int) > 0:
@@ -770,10 +771,13 @@ def parse_impdosfiles(folder, natom, nspin, ef, use_lmdos):
     dos[:, :, 0] = (dos[:, :, 0] - ef.value) * eVscale
     dos[:, :, 1:] = dos[:, :, 1:] / eVscale
     if len(dos_int) > 0:
-        dos_int[:, :, 0] = (dos_int[:, :, 0] - ef.value) * eVscale
-        dos_int[:, :, 1:] = dos_int[:, :, 1:] / eVscale
-    else :
-        pass 
+        try:
+            dos_int[:, :, 0] = (dos_int[:, :, 0] - ef.value) * eVscale
+            dos_int[:, :, 1:] = dos_int[:, :, 1:] / eVscale
+        except:
+            message = 'Not enough data were present in the interpolated dos, due to a lack of energy points to interpolate on'
+    else:
+        pass
 
     # create output nodes
     dosnode = XyData()