added center of faces to brillouin zone plot

src/solid_state_tools.py

@@ -810,6 +810,61 @@ def replace_greek(x):
     return res
 
 
+def find_center_of_faces(coordinates, triangles):
+    faces = construct_faces(coordinates, triangles)
+    centers = np.zeros((len(faces), 3))
+
+    for i, vertex in enumerate(faces):
+        centers[i, :] = np.mean(coordinates[list(vertex), :], axis=0)
+    return centers
+
+
+def is_in_plane(base_triangle, test_triangle, coordinates):
+    normal = construct_normal_from_triangle(coordinates, base_triangle)
+    normal_alt = construct_normal_from_triangle(coordinates, test_triangle)
+
+    if base_triangle[0] == test_triangle[0]:
+        test_vec = coordinates[base_triangle[0], :] - coordinates[test_triangle[1], :]
+    else:
+        test_vec = coordinates[base_triangle[0], :] - coordinates[test_triangle[0], :]
+
+    test_vec = test_vec/np.linalg.norm(test_vec)
+    return (abs(np.dot(normal, normal_alt)) > 0.999) and (abs(np.dot(normal, test_vec)) < 0.001)
+
+
+def construct_normal_from_triangle(coordinates, triangle):
+    v1 = coordinates[triangle[1], :] - coordinates[triangle[0], :]
+    v2 = coordinates[triangle[2], :] - coordinates[triangle[0], :]
+    normal = np.cross(v1, v2)
+    return normal/np.linalg.norm(normal)
+
+
+def is_list_in_list_of_lists(small_list, large_list):
+    for el in large_list:
+        if sorted(small_list) == sorted(el):
+            return True
+    return False
+
+
+def construct_faces(points, triangles):
+    used_triangles = []
+    faces = []
+    for base_triangle in triangles:
+        if is_list_in_list_of_lists(base_triangle, used_triangles):
+            continue
+        used_triangles.append(base_triangle)
+        face = set(base_triangle)
+        for test_triangle in triangles:
+            if is_list_in_list_of_lists(test_triangle, used_triangles):
+                continue
+            if is_in_plane(base_triangle, test_triangle, points):
+                face.update(test_triangle)
+                used_triangles.append(test_triangle)
+        faces.append(list(face))
+
+    return faces
+
+
 def get_materials_dos_from_id(id):
     with MPRester(API_KEY) as m:
         dos = m.get_dos_by_material_id(id)
@@ -945,3 +1000,5 @@ def get_materials_band_structure_from_id(id):
     # out = parser.parse_cif_file('/home/jannick/OpenDFT_projects/LiBH4/1504402.cif')
     # general_handler = GeneralHandler()
     # print(general_handler.available_handlers)
+
+

src/visualization.py

@@ -109,6 +109,8 @@ def __init__(self, parent):
         self.text_plots = []
         self.glyph_points = None
         self.picker = None
+        self.w_points = None
+        self.centers = None
 
     def clear_plot(self):
         self.scene.mlab.clf(figure=self.scene.mayavi_scene)
@@ -117,6 +119,7 @@ def set_crystal_structure(self, crystal_structure):
         self.crystal_structure = crystal_structure
         self.w_points = sst.construct_brillouin_vertices(crystal_structure)
         self.brillouin_edges = sst.construct_convex_hull(self.w_points)
+        self.centers = sst.find_center_of_faces(self.w_points, self.brillouin_edges)
 
     def set_path(self, k_path):
         self.k_path = k_path
@@ -171,12 +174,8 @@ def plot_path(self):
     def plot_brillouin_zone(self, plot_connections=True):
 
         self.wpoints_plot = np.append(self.w_points, np.array([[0, 0, 0]]), axis=0)
-        # TODO find general way to make the center of facets points
-        # for i in [-1,1]:
-        #     for j in range(3):
-        #         self.wpoints_plot = np.append(self.wpoints_plot, np.array([i*0.5*self.crystal_structure.inv_lattice_vectors[j,:]]), axis=0)
+        self.wpoints_plot = np.append(self.wpoints_plot, self.centers, axis=0)
 
-        # self.scene.mlab.points3d(0.0, 0.0, 0.0, color=(0.7, 0.7, 0.7), scale_factor=.1, figure=self.scene.mayavi_scene)
         self.plot_of_vertices = self.scene.mlab.points3d(self.wpoints_plot[:, 0], self.wpoints_plot[:, 1],
                                                          self.wpoints_plot[:, 2], color=(0.7, 0.7, 0.7),
                                                          scale_factor=.1, figure=self.scene.mayavi_scene)
@@ -186,12 +185,6 @@ def plot_brillouin_zone(self, plot_connections=True):
                                         self.brillouin_edges, opacity=0.3, color=(0.5, 0.5, 0.5), tube_radius=2,
                                         figure=self.scene.mayavi_scene)
 
-        # if plot_connections:
-        #     for i, connection in enumerate(self.brillouin_edges):
-        #         for con in connection:
-        #             bond = [i, con]
-        #             self.scene.mlab.plot3d(self.w_points[bond, 0], self.w_points[bond, 1], self.w_points[bond, 2],figure=self.scene.mayavi_scene,tube_radius=0.01)
-
         self.plot_unit_vectors()
 
         self.outline = self.scene.mlab.outline(line_width=3, figure=self.scene.mayavi_scene)
@@ -1724,8 +1717,8 @@ def set_dark_mode_matplotlib(f, ax, color):
 if __name__ == "__main__":
     from solid_state_tools import CrystalStructure, PhononEigenvectors,StructureParser
 
-    atoms = np.array([[0, 0, 0, 6], [1/3, 1/3, 0, 6]])
-    unit_cell = 4.65 * np.array([[0.5,  0.8660254040,  0.0], [-0.5,0.8660254040,0.0], [0, 0, 6.0]])
+    atoms = np.array([[0, 0, 0, 6], [1/4, 1/4, 1/4, 6]])
+    unit_cell = 6.6 * np.array([[0,  0.5,  0.5], [0.5,0.0,0.5], [0.5, 0.5, 0.0]])
 
     # unit_cell = 20*np.eye(3,3)
     # N = 50
@@ -1741,20 +1734,22 @@ def set_dark_mode_matplotlib(f, ax, color):
 
     freqs = np.ones((1, 1))
     # mode = np.array([np.random.randn(3*N_atoms)])[None,:,:]
-    mode = np.array([[0,0,1,0,0,1]])[None,:,:]
-
-    k_vector = np.array([[0.1,0.1,0]])
-
-    eig = PhononEigenvectors(freqs,mode,k_vector)
+    # mode = np.array([[0,0,1,0,0,1]])[None,:,:]
+    #
+    # k_vector = np.array([[0.1,0.1,0]])
+    #
+    # eig = PhononEigenvectors(freqs,mode,k_vector)
+    #
+    # vis = PhononVisualization(crystal_structure)
+    # vis.phonon_eigenvectors = eig
+    # vis.plot_phonons(0,0)
+    # vis.configure_traits()
 
-    vis = PhononVisualization(crystal_structure)
-    vis.phonon_eigenvectors = eig
-    vis.plot_phonons(0,0)
+    vis = BrillouinVisualization(None)
+    vis.set_crystal_structure(crystal_structure)
     vis.configure_traits()
 
 
-
-
     # x, y, z = np.ogrid[-5:5:64j, -5:5:64j, -5:5:64j]
     # data = np.sin(3 * x) / x + 0.05 * z ** 2 + np.cos(3 * y)
     #