Merge remote-tracking branch 'origin/master'

Author: DanPorter

Dans_Diffraction/classes_crystal.py

@@ -365,6 +365,96 @@ def start_gui(self):
         except ImportError:
             print('Sorry, you need to install tkinter!')
 
+    def search_distances(self, min_d=0.65, max_d=3.20, c_ele=None, elems=None,
+                         labels=None, simple=True):
+        """
+        Calculated atoms interatomic distances form each label.
+        :param c_ele (list,string): only sites with noted elements
+                                    if None all site
+        :param elems (list,string): only distances with noted elements
+                                    if None all site
+        :param min_d: minimum distance 
+        :param max_d: maximum distance
+        :return dictionary: 
+
+        """
+        xyz = self.Cell.calculateR(self.Atoms.uvw())
+        UVstar = self.Cell.UVstar()
+        UV = self.Cell.UV()
+        Lpar = self.Cell.lp()[:3]
+
+        Lran = np.array([0, 0, 0, 0, 0, 0])
+
+        for i in range(3):
+            Lran[i] = min(fg.distance2plane([0, 0, 0], UVstar[i], xyz))
+        for i in range(3):
+            Lran[3 + i] = min(fg.distance2plane(UV[i], UVstar[i] + UV[i], xyz))
+        for i, L in enumerate(Lran):
+            Lran[i] = np.ceil(max_d / Lpar[i % 3]) if L < max_d else 0
+
+        IntRes = self.Structure.generate_lattice(U=Lran[0] + Lran[3],
+                                                 V=Lran[1] + Lran[4],
+                                                 W=Lran[2] + Lran[5],
+                                                 centred=False)[:3]
+        B_uvw, B_Ele, B_label = IntRes
+        B_label, B_Ele = np.asarray(B_label), np.asarray(B_Ele)
+
+        B_uvw -= np.array(Lran[:3])
+        B_xyz = self.Cell.calculateR(B_uvw)
+
+        if c_ele is None:
+            c_ele = set(self.Atoms.type)
+        elif isinstance(c_ele, str):
+            c_ele = [c_ele]
+
+        if elems is None:
+            elems = set(self.Atoms.type)
+        elif isinstance(elems, str):
+            elems = [elems]
+
+        if labels is None:
+            labels = self.Atoms.label
+        elif isinstance(labels, str):
+            labels = [labels]
+
+        distances = {}
+        for i, atom in enumerate(xyz):
+            if not self.Atoms.type[i] in c_ele:
+                continue
+            if not self.Atoms.label[i] in labels:
+                continue
+            s_dist = {}
+            vdist = (B_xyz - atom) ** 2
+            vdist = np.sqrt(np.sum(vdist, axis=1))
+            cond1 = (vdist > min_d) * (vdist < max_d)
+            vlabel = B_label[cond1]
+            vele = B_Ele[cond1]
+            vdist = vdist[cond1]
+
+            Ord = np.argsort(vdist)
+            cond2 = [ele in elems for ele in vele[Ord]]
+            s_dist = {'dist': vdist[Ord][cond2],
+                      'label': list(vlabel[Ord][cond2]),
+                      'type': list(vele[Ord][cond2])}
+
+            distances[self.Atoms.label[i]] = s_dist
+
+        if simple:
+            # reduce the output for mixed site occupation
+            lab = [i for i,j in zip(self.Atoms.label, self.Atoms.type) if j in c_ele]
+            for i, site_i in enumerate(lab):
+                for site_j in lab[i + 1:]:
+                    uvw_i = self.Atoms[site_i].uvw()
+                    uvw_j = self.Atoms[site_j].uvw()
+                    if all(uvw_i == uvw_j):
+                        lab.remove(site_j)
+                        try:
+                            del distances[site_j]
+                        except KeyError:
+                            pass
+
+        return distances
+
     def __add__(self, other):
         return MultiCrystal([self, other])
 
@@ -564,7 +654,7 @@ def indexQ(self, Q):
         Convert coordinates [x,y,z], in an orthogonal basis, to
         coordinates [h,k,l], in the basis of the reciprocal lattice
                     H(h,k,l) = Q(x,y,z) / [A*,B*,C*]
-        
+
         E.G.
             HKL = indexQ([2.2046264, 1.2728417, 0.0000000]) # for a hexagonal system, a = 2.85
             > HKL = [1,0,0]
@@ -577,7 +667,6 @@ def calculateR(self, UVW):
         Convert coordinates [u,v,w], in the basis of the unit cell, to
         coordinates [x,y,z], in an orthogonal basis, in units of A
                     R(x,y,z) = uA + vB + wC
-        
         E.G.
             R = Cell.calculateR([0.1,0,0]) # for a hexagonal system, a = 2.85
             > R = array([[0.285, 0, 0]])
@@ -1030,6 +1119,8 @@ def __call__(self, u=[0], v=[0], w=[0], type=None,
         self.__init__(u, v, w, type, label, occupancy, uiso, mxmymz=None)
 
     def __getitem__(self, idx):
+        if isinstance(idx, str):
+            idx = self.label.index(idx)
         return self.atom(idx)
 
     def fromcif(self, cifvals):

Dans_Diffraction/classes_plotting.py

@@ -143,7 +143,43 @@ def plot_crystal(self, show_labels=False):
         plt.legend(fontsize=24, frameon=False)
 
         plt.title(self.xtl.name, fontsize=28, fontweight='bold')
-    
+
+    def plot_distance(self, min_d=0.65, max_d=3.20, labels=None,
+                      c_ele=None, elems=None, ranges=None, step=0.04):
+        """
+        Plot atoms interatomic distances form each label.
+        :param c_ele (list,string): only sites with noted elements
+                                    if None all site
+        :param elems (list,string): only distances with noted elements
+                                    if None all site
+        :param min_d: minimum distance
+        :param max_d: maximum distance
+        :return:
+        """
+        dist = self.xtl.search_distances(c_ele=c_ele, elems=elems,
+                                         labels=labels, min_d=min_d,
+                                         max_d=max_d)
+
+        if ranges is None:
+            all_d = np.hstack([i['dist'] for i in dist.values()])
+            ranges = (np.floor(min(all_d)), np.ceil(max(all_d)))
+        # Create plot
+        if len(dist) == 0:
+            print('no distance present')
+            return
+        fig, axs = plt.subplots(len(dist), constrained_layout=True)
+        if len(dist) == 1:
+            axs = [axs]
+        for i, site in enumerate(dist):
+            axs[i].set_title(site)
+            axs[i].hist(dist[site]['dist'], range=ranges,
+                        bins=int((ranges[1] - ranges[0]) / step))
+            axs[i].set(ylabel='n. Atoms')
+        axs[-1].set(xlabel='$\AA$')
+        if len(dist) > 1:
+            for ax in axs.flat:
+                ax.label_outer()
+
     def plot_layers(self, layers=None, layer_axis=2, layer_width=0.05, show_labels=False):
         """
         Separate the structure into layers along the chosen axis

Dans_Diffraction/functions_general.py

@@ -497,6 +497,25 @@ def distance2line(line_start, line_end, point):
     return np.sqrt(np.sum(vec_arb ** 2))
 
 
+def distance2plane(line_start, line_end, point):
+    """
+    Calculate distance from a plane (defined by a perpendicular vector) to an arbitary point in space
+    :param vector_start: array, position lying on plane and vector start
+    :param line_end:  array, position of the end of the vector
+    :param point: array, arbitary position in space
+    :return: float
+    """
+    line_start = np.asarray(line_start)
+    line_end = np.asarray(line_end)
+    point = np.asarray(point)
+
+    line_diff = line_end - line_start
+    unit_line = line_diff / np.sqrt(np.sum(line_diff ** 2))
+
+    vec_arb = np.dot((point - line_start), unit_line)
+    return abs(vec_arb)
+
+
 def vector_intersection(point1, direction1, point2, direction2):
     """
     Calculate the point in 2D where two lines cross.