1# coding: utf-8 2# Copyright (c) Pymatgen Development Team. 3# Distributed under the terms of the MIT License. 4 5""" 6This module define a WulffShape class to generate the Wulff shape from 7a lattice, a list of indices and their corresponding surface energies, 8and the total area and volume of the wulff shape,the weighted surface energy, 9the anisotropy and shape_factor can also be calculated. 10In support of plotting from a given view in terms of miller index. 11 12The lattice is from the conventional unit cell, and (hkil) for hexagonal 13lattices. 14 15If you use this code extensively, consider citing the following: 16 17Tran, R.; Xu, Z.; Radhakrishnan, B.; Winston, D.; Persson, K. A.; Ong, S. P. 18(2016). Surface energies of elemental crystals. Scientific Data. 19""" 20 21import itertools 22import logging 23import warnings 24 25import numpy as np 26import plotly.graph_objs as go 27from scipy.spatial import ConvexHull 28 29from pymatgen.core.structure import Structure 30from pymatgen.util.coord import get_angle 31from pymatgen.util.string import unicodeify_spacegroup 32 33__author__ = "Zihan Xu, Richard Tran, Shyue Ping Ong" 34__copyright__ = "Copyright 2013, The Materials Virtual Lab" 35__version__ = "0.1" 36__maintainer__ = "Zihan Xu" 37__email__ = "zix009@eng.ucsd.edu" 38__date__ = "May 5 2016" 39 40logger = logging.getLogger(__name__) 41 42 43def hkl_tuple_to_str(hkl): 44 """ 45 Prepare for display on plots 46 "(hkl)" for surfaces 47 Agrs: 48 hkl: in the form of [h, k, l] or (h, k, l) 49 """ 50 str_format = "($" 51 for x in hkl: 52 if x < 0: 53 str_format += "\\overline{" + str(-x) + "}" 54 else: 55 str_format += str(x) 56 str_format += "$)" 57 return str_format 58 59 60def get_tri_area(pts): 61 """ 62 Given a list of coords for 3 points, 63 Compute the area of this triangle. 64 65 Args: 66 pts: [a, b, c] three points 67 """ 68 a, b, c = pts[0], pts[1], pts[2] 69 v1 = np.array(b) - np.array(a) 70 v2 = np.array(c) - np.array(a) 71 area_tri = abs(np.linalg.norm(np.cross(v1, v2)) / 2) 72 return area_tri 73 74 75class WulffFacet: 76 """ 77 Helper container for each Wulff plane. 78 """ 79 80 def __init__(self, normal, e_surf, normal_pt, dual_pt, index, m_ind_orig, miller): 81 """ 82 :param normal: 83 :param e_surf: 84 :param normal_pt: 85 :param dual_pt: 86 :param index: 87 :param m_ind_orig: 88 :param miller: 89 """ 90 self.normal = normal 91 self.e_surf = e_surf 92 self.normal_pt = normal_pt 93 self.dual_pt = dual_pt 94 self.index = index 95 self.m_ind_orig = m_ind_orig 96 self.miller = miller 97 self.points = [] 98 self.outer_lines = [] 99 100 101class WulffShape: 102 """ 103 Generate Wulff Shape from list of miller index and surface energies, 104 with given conventional unit cell. 105 surface energy (Jm^2) is the length of normal. 106 107 Wulff shape is the convex hull. 108 Based on: 109 http://scipy.github.io/devdocs/generated/scipy.spatial.ConvexHull.html 110 111 Process: 112 1. get wulff simplices 113 2. label with color 114 3. get wulff_area and other properties 115 116 .. attribute:: debug (bool) 117 118 .. attribute:: alpha 119 transparency 120 121 .. attribute:: color_set 122 123 .. attribute:: grid_off (bool) 124 125 .. attribute:: axis_off (bool) 126 127 .. attribute:: show_area 128 129 .. attribute:: off_color 130 color of facets off wulff 131 132 .. attribute:: structure 133 Structure object, input conventional unit cell (with H ) from lattice 134 135 .. attribute:: miller_list 136 list of input miller index, for hcp in the form of hkil 137 138 .. attribute:: hkl_list 139 modify hkill to hkl, in the same order with input_miller 140 141 .. attribute:: e_surf_list 142 list of input surface energies, in the same order with input_miller 143 144 .. attribute:: lattice 145 Lattice object, the input lattice for the conventional unit cell 146 147 .. attribute:: facets 148 [WulffFacet] for all facets considering symm 149 150 .. attribute:: dual_cv_simp 151 simplices from the dual convex hull (dual_pt) 152 153 .. attribute:: wulff_pt_list 154 155 .. attribute:: wulff_cv_simp 156 simplices from the convex hull of wulff_pt_list 157 158 .. attribute:: on_wulff 159 list for all input_miller, True is on wulff. 160 161 .. attribute:: color_area 162 list for all input_miller, total area on wulff, off_wulff = 0. 163 164 .. attribute:: miller_area 165 ($hkl$): area for all input_miller 166 167 """ 168 169 def __init__(self, lattice, miller_list, e_surf_list, symprec=1e-5): 170 """ 171 Args: 172 lattice: Lattice object of the conventional unit cell 173 miller_list ([(hkl), ...]: list of hkl or hkil for hcp 174 e_surf_list ([float]): list of corresponding surface energies 175 symprec (float): for recp_operation, default is 1e-5. 176 """ 177 178 if any(se < 0 for se in e_surf_list): 179 warnings.warn("Unphysical (negative) surface energy detected.") 180 181 self.color_ind = list(range(len(miller_list))) 182 183 self.input_miller_fig = [hkl_tuple_to_str(x) for x in miller_list] 184 # store input data 185 self.structure = Structure(lattice, ["H"], [[0, 0, 0]]) 186 self.miller_list = tuple(tuple(x) for x in miller_list) 187 self.hkl_list = tuple((x[0], x[1], x[-1]) for x in miller_list) 188 self.e_surf_list = tuple(e_surf_list) 189 self.lattice = lattice 190 self.symprec = symprec 191 192 # 2. get all the data for wulff construction 193 # get all the surface normal from get_all_miller_e() 194 self.facets = self._get_all_miller_e() 195 logger.debug(len(self.facets)) 196 197 # 3. consider the dual condition 198 dual_pts = [x.dual_pt for x in self.facets] 199 dual_convex = ConvexHull(dual_pts) 200 dual_cv_simp = dual_convex.simplices 201 # simplices (ndarray of ints, shape (nfacet, ndim)) 202 # list of [i, j, k] , ndim = 3 203 # i, j, k: ind for normal_e_m 204 # recalculate the dual of dual, get the wulff shape. 205 # conner <-> surface 206 # get cross point from the simplices of the dual convex hull 207 wulff_pt_list = [self._get_cross_pt_dual_simp(dual_simp) for dual_simp in dual_cv_simp] 208 209 wulff_convex = ConvexHull(wulff_pt_list) 210 wulff_cv_simp = wulff_convex.simplices 211 logger.debug(", ".join([str(len(x)) for x in wulff_cv_simp])) 212 213 # store simplices and convex 214 self.dual_cv_simp = dual_cv_simp 215 self.wulff_pt_list = wulff_pt_list 216 self.wulff_cv_simp = wulff_cv_simp 217 self.wulff_convex = wulff_convex 218 219 self.on_wulff, self.color_area = self._get_simpx_plane() 220 221 miller_area = [] 222 for m, in_mill_fig in enumerate(self.input_miller_fig): 223 miller_area.append(in_mill_fig + " : " + str(round(self.color_area[m], 4))) 224 self.miller_area = miller_area 225 226 def _get_all_miller_e(self): 227 """ 228 from self: 229 get miller_list(unique_miller), e_surf_list and symmetry 230 operations(symmops) according to lattice 231 apply symmops to get all the miller index, then get normal, 232 get all the facets functions for wulff shape calculation: 233 |normal| = 1, e_surf is plane's distance to (0, 0, 0), 234 normal[0]x + normal[1]y + normal[2]z = e_surf 235 236 return: 237 [WulffFacet] 238 """ 239 all_hkl = [] 240 color_ind = self.color_ind 241 planes = [] 242 recp = self.structure.lattice.reciprocal_lattice_crystallographic 243 recp_symmops = self.lattice.get_recp_symmetry_operation(self.symprec) 244 245 for i, (hkl, energy) in enumerate(zip(self.hkl_list, self.e_surf_list)): 246 for op in recp_symmops: 247 miller = tuple(int(x) for x in op.operate(hkl)) 248 if miller not in all_hkl: 249 all_hkl.append(miller) 250 normal = recp.get_cartesian_coords(miller) 251 normal /= np.linalg.norm(normal) 252 normal_pt = [x * energy for x in normal] 253 dual_pt = [x / energy for x in normal] 254 color_plane = color_ind[divmod(i, len(color_ind))[1]] 255 planes.append(WulffFacet(normal, energy, normal_pt, dual_pt, color_plane, i, hkl)) 256 257 # sort by e_surf 258 planes.sort(key=lambda x: x.e_surf) 259 return planes 260 261 def _get_cross_pt_dual_simp(self, dual_simp): 262 """ 263 |normal| = 1, e_surf is plane's distance to (0, 0, 0), 264 plane function: 265 normal[0]x + normal[1]y + normal[2]z = e_surf 266 267 from self: 268 normal_e_m to get the plane functions 269 dual_simp: (i, j, k) simplices from the dual convex hull 270 i, j, k: plane index(same order in normal_e_m) 271 """ 272 matrix_surfs = [self.facets[dual_simp[i]].normal for i in range(3)] 273 matrix_e = [self.facets[dual_simp[i]].e_surf for i in range(3)] 274 cross_pt = np.dot(np.linalg.inv(matrix_surfs), matrix_e) 275 return cross_pt 276 277 def _get_simpx_plane(self): 278 """ 279 Locate the plane for simpx of on wulff_cv, by comparing the center of 280 the simpx triangle with the plane functions. 281 """ 282 on_wulff = [False] * len(self.miller_list) 283 surface_area = [0.0] * len(self.miller_list) 284 for simpx in self.wulff_cv_simp: 285 pts = [self.wulff_pt_list[simpx[i]] for i in range(3)] 286 center = np.sum(pts, 0) / 3.0 287 # check whether the center of the simplices is on one plane 288 for plane in self.facets: 289 abs_diff = abs(np.dot(plane.normal, center) - plane.e_surf) 290 if abs_diff < 1e-5: 291 on_wulff[plane.index] = True 292 surface_area[plane.index] += get_tri_area(pts) 293 294 plane.points.append(pts) 295 plane.outer_lines.append([simpx[0], simpx[1]]) 296 plane.outer_lines.append([simpx[1], simpx[2]]) 297 plane.outer_lines.append([simpx[0], simpx[2]]) 298 # already find the plane, move to the next simplices 299 break 300 for plane in self.facets: 301 plane.outer_lines.sort() 302 plane.outer_lines = [line for line in plane.outer_lines if plane.outer_lines.count(line) != 2] 303 return on_wulff, surface_area 304 305 def _get_colors(self, color_set, alpha, off_color, custom_colors={}): 306 """ 307 assign colors according to the surface energies of on_wulff facets. 308 309 return: 310 (color_list, color_proxy, color_proxy_on_wulff, miller_on_wulff, 311 e_surf_on_wulff_list) 312 """ 313 import matplotlib as mpl 314 import matplotlib.pyplot as plt 315 316 color_list = [off_color] * len(self.hkl_list) 317 color_proxy_on_wulff = [] 318 miller_on_wulff = [] 319 e_surf_on_wulff = [(i, e_surf) for i, e_surf in enumerate(self.e_surf_list) if self.on_wulff[i]] 320 321 c_map = plt.get_cmap(color_set) 322 e_surf_on_wulff.sort(key=lambda x: x[1], reverse=False) 323 e_surf_on_wulff_list = [x[1] for x in e_surf_on_wulff] 324 if len(e_surf_on_wulff) > 1: 325 cnorm = mpl.colors.Normalize(vmin=min(e_surf_on_wulff_list), vmax=max(e_surf_on_wulff_list)) 326 else: 327 # if there is only one hkl on wulff, choose the color of the median 328 cnorm = mpl.colors.Normalize( 329 vmin=min(e_surf_on_wulff_list) - 0.1, 330 vmax=max(e_surf_on_wulff_list) + 0.1, 331 ) 332 scalar_map = mpl.cm.ScalarMappable(norm=cnorm, cmap=c_map) 333 334 for i, e_surf in e_surf_on_wulff: 335 color_list[i] = scalar_map.to_rgba(e_surf, alpha=alpha) 336 if tuple(self.miller_list[i]) in custom_colors.keys(): 337 color_list[i] = custom_colors[tuple(self.miller_list[i])] 338 color_proxy_on_wulff.append(plt.Rectangle((2, 2), 1, 1, fc=color_list[i], alpha=alpha)) 339 miller_on_wulff.append(self.input_miller_fig[i]) 340 scalar_map.set_array([x[1] for x in e_surf_on_wulff]) 341 color_proxy = [plt.Rectangle((2, 2), 1, 1, fc=x, alpha=alpha) for x in color_list] 342 343 return ( 344 color_list, 345 color_proxy, 346 color_proxy_on_wulff, 347 miller_on_wulff, 348 e_surf_on_wulff_list, 349 ) 350 351 def show(self, *args, **kwargs): 352 r""" 353 Show the Wulff plot. 354 355 Args: 356 *args: Passed to get_plot. 357 **kwargs: Passed to get_plot. 358 """ 359 self.get_plot(*args, **kwargs).show() 360 361 def get_line_in_facet(self, facet): 362 """ 363 Returns the sorted pts in a facet used to draw a line 364 """ 365 366 lines = list(facet.outer_lines) 367 pt = [] 368 prev = None 369 while len(lines) > 0: 370 if prev is None: 371 l = lines.pop(0) 372 else: 373 for i, l in enumerate(lines): 374 if prev in l: 375 l = lines.pop(i) 376 if l[1] == prev: 377 l.reverse() 378 break 379 # make sure the lines are connected one by one. 380 # find the way covering all pts and facets 381 pt.append(self.wulff_pt_list[l[0]].tolist()) 382 pt.append(self.wulff_pt_list[l[1]].tolist()) 383 prev = l[1] 384 385 return pt 386 387 def get_plot( 388 self, 389 color_set="PuBu", 390 grid_off=True, 391 axis_off=True, 392 show_area=False, 393 alpha=1, 394 off_color="red", 395 direction=None, 396 bar_pos=(0.75, 0.15, 0.05, 0.65), 397 bar_on=False, 398 units_in_JPERM2=True, 399 legend_on=True, 400 aspect_ratio=(8, 8), 401 custom_colors={}, 402 ): 403 """ 404 Get the Wulff shape plot. 405 406 Args: 407 color_set: default is 'PuBu' 408 grid_off (bool): default is True 409 axis_off (bool): default is Ture 410 show_area (bool): default is False 411 alpha (float): chosen from 0 to 1 (float), default is 1 412 off_color: Default color for facets not present on the Wulff shape. 413 direction: default is (1, 1, 1) 414 bar_pos: default is [0.75, 0.15, 0.05, 0.65] 415 bar_on (bool): default is False 416 legend_on (bool): default is True 417 aspect_ratio: default is (8, 8) 418 custom_colors ({(h,k,l}: [r,g,b,alpha}): Customize color of each 419 facet with a dictionary. The key is the corresponding Miller 420 index and value is the color. Undefined facets will use default 421 color site. Note: If you decide to set your own colors, it 422 probably won't make any sense to have the color bar on. 423 units_in_JPERM2 (bool): Units of surface energy, defaults to 424 Joules per square meter (True) 425 426 Return: 427 (matplotlib.pyplot) 428 """ 429 import matplotlib as mpl 430 import matplotlib.pyplot as plt 431 import mpl_toolkits.mplot3d as mpl3 432 433 ( 434 color_list, 435 color_proxy, 436 color_proxy_on_wulff, 437 miller_on_wulff, 438 e_surf_on_wulff, 439 ) = self._get_colors(color_set, alpha, off_color, custom_colors=custom_colors) 440 441 if not direction: 442 # If direction is not specified, use the miller indices of 443 # maximum area. 444 direction = max(self.area_fraction_dict.items(), key=lambda x: x[1])[0] 445 446 fig = plt.figure() 447 fig.set_size_inches(aspect_ratio[0], aspect_ratio[1]) 448 azim, elev = self._get_azimuth_elev([direction[0], direction[1], direction[-1]]) 449 450 wulff_pt_list = self.wulff_pt_list 451 452 ax = mpl3.Axes3D(fig, azim=azim, elev=elev) 453 454 for plane in self.facets: 455 # check whether [pts] is empty 456 if len(plane.points) < 1: 457 # empty, plane is not on_wulff. 458 continue 459 # assign the color for on_wulff facets according to its 460 # index and the color_list for on_wulff 461 plane_color = color_list[plane.index] 462 pt = self.get_line_in_facet(plane) 463 # plot from the sorted pts from [simpx] 464 tri = mpl3.art3d.Poly3DCollection([pt]) 465 tri.set_color(plane_color) 466 tri.set_edgecolor("#808080") 467 ax.add_collection3d(tri) 468 469 # set ranges of x, y, z 470 # find the largest distance between on_wulff pts and the origin, 471 # to ensure complete and consistent display for all directions 472 r_range = max([np.linalg.norm(x) for x in wulff_pt_list]) 473 ax.set_xlim([-r_range * 1.1, r_range * 1.1]) 474 ax.set_ylim([-r_range * 1.1, r_range * 1.1]) 475 ax.set_zlim([-r_range * 1.1, r_range * 1.1]) # pylint: disable=E1101 476 # add legend 477 if legend_on: 478 color_proxy = color_proxy 479 if show_area: 480 ax.legend( 481 color_proxy, 482 self.miller_area, 483 loc="upper left", 484 bbox_to_anchor=(0, 1), 485 fancybox=True, 486 shadow=False, 487 ) 488 else: 489 ax.legend( 490 color_proxy_on_wulff, 491 miller_on_wulff, 492 loc="upper center", 493 bbox_to_anchor=(0.5, 1), 494 ncol=3, 495 fancybox=True, 496 shadow=False, 497 ) 498 ax.set_xlabel("x") 499 ax.set_ylabel("y") 500 ax.set_zlabel("z") 501 502 # Add colorbar 503 if bar_on: 504 cmap = plt.get_cmap(color_set) 505 cmap.set_over("0.25") 506 cmap.set_under("0.75") 507 bounds = [round(e, 2) for e in e_surf_on_wulff] 508 bounds.append(1.2 * bounds[-1]) 509 norm = mpl.colors.BoundaryNorm(bounds, cmap.N) 510 # display surface energies 511 ax1 = fig.add_axes(bar_pos) 512 cbar = mpl.colorbar.ColorbarBase( 513 ax1, 514 cmap=cmap, 515 norm=norm, 516 boundaries=[0] + bounds + [10], 517 extend="both", 518 ticks=bounds[:-1], 519 spacing="proportional", 520 orientation="vertical", 521 ) 522 units = "$J/m^2$" if units_in_JPERM2 else r"$eV/\AA^2$" 523 cbar.set_label("Surface Energies (%s)" % (units), fontsize=25) 524 525 if grid_off: 526 ax.grid("off") 527 if axis_off: 528 ax.axis("off") 529 return plt 530 531 def get_plotly( 532 self, 533 color_set="PuBu", 534 off_color="red", 535 alpha=1, 536 custom_colors={}, 537 units_in_JPERM2=True, 538 ): 539 """ 540 Get the Wulff shape as a plotly Figure object. 541 542 Args: 543 color_set: default is 'PuBu' 544 alpha (float): chosen from 0 to 1 (float), default is 1 545 off_color: Default color for facets not present on the Wulff shape. 546 custom_colors ({(h,k,l}: [r,g,b,alpha}): Customize color of each 547 facet with a dictionary. The key is the corresponding Miller 548 index and value is the color. Undefined facets will use default 549 color site. Note: If you decide to set your own colors, it 550 probably won't make any sense to have the color bar on. 551 units_in_JPERM2 (bool): Units of surface energy, defaults to 552 Joules per square meter (True) 553 554 Return: 555 (plotly.graph_objs.Figure) 556 """ 557 558 units = "Jm⁻²" if units_in_JPERM2 else "eVÅ⁻²" 559 ( 560 color_list, 561 color_proxy, 562 color_proxy_on_wulff, 563 miller_on_wulff, 564 e_surf_on_wulff, 565 ) = self._get_colors(color_set, alpha, off_color, custom_colors=custom_colors) 566 567 planes_data, color_scale, ticktext, tickvals = [], [], [], [] 568 for plane in self.facets: 569 if len(plane.points) < 1: 570 # empty, plane is not on_wulff. 571 continue 572 573 plane_color = color_list[plane.index] 574 plane_color = (1, 0, 0, 1) if plane_color == off_color else plane_color # set to red for now 575 576 pt = self.get_line_in_facet(plane) 577 x_pts, y_pts, z_pts = [], [], [] 578 for p in pt: 579 x_pts.append(p[0]) 580 y_pts.append(p[1]) 581 z_pts.append(p[2]) 582 583 # remove duplicate x y z pts to save time 584 all_xyz = [] 585 # pylint: disable=E1133,E1136 586 [all_xyz.append(list(coord)) for coord in np.array([x_pts, y_pts, z_pts]).T if list(coord) not in all_xyz] 587 all_xyz = np.array(all_xyz).T 588 x_pts, y_pts, z_pts = all_xyz[0], all_xyz[1], all_xyz[2] 589 index_list = [int(i) for i in np.linspace(0, len(x_pts) - 1, len(x_pts))] 590 591 tri_indices = np.array(list(itertools.combinations(index_list, 3))).T 592 hkl = self.miller_list[plane.index] 593 hkl = unicodeify_spacegroup("(" + "%s" * len(hkl) % hkl + ")") 594 color = "rgba(%.5f, %.5f, %.5f, %.5f)" % tuple(np.array(plane_color) * 255) 595 596 # note hoverinfo is incompatible with latex, need unicode instead 597 planes_data.append( 598 go.Mesh3d( 599 x=x_pts, 600 y=y_pts, 601 z=z_pts, 602 i=tri_indices[0], 603 j=tri_indices[1], 604 k=tri_indices[2], 605 hovertemplate="<br>%{text}<br>" + "%s=%.3f %s<br>" % ("\u03b3", plane.e_surf, units), 606 color=color, 607 text=[r"Miller index: %s" % hkl] * len(x_pts), 608 hoverinfo="name", 609 name="", 610 ) 611 ) 612 613 # normalize surface energy from a scale of 0 to 1 for colorbar 614 norm_e = (plane.e_surf - min(e_surf_on_wulff)) / (max(e_surf_on_wulff) - min(e_surf_on_wulff)) 615 c = [norm_e, color] 616 if c not in color_scale: 617 color_scale.append(c) 618 ticktext.append("%.3f" % plane.e_surf) 619 tickvals.append(norm_e) 620 621 # Add colorbar 622 color_scale = sorted(color_scale, key=lambda c: c[0]) 623 colorbar = go.Mesh3d( 624 x=[0], 625 y=[0], 626 z=[0], 627 colorbar=go.ColorBar( 628 title={ 629 "text": r"Surface energy %s" % units, 630 "side": "right", 631 "font": {"size": 25}, 632 }, 633 ticktext=ticktext, 634 tickvals=tickvals, 635 ), 636 colorscale=[[0, "rgb(255,255,255, 255)"]] + color_scale, # fix the scale 637 intensity=[0, 0.33, 0.66, 1], 638 i=[0], 639 j=[0], 640 k=[0], 641 name="y", 642 showscale=True, 643 ) 644 planes_data.append(colorbar) 645 646 # Format aesthetics: background, axis, etc. 647 axis_dict = dict( 648 title="", 649 autorange=True, 650 showgrid=False, 651 zeroline=False, 652 ticks="", 653 showline=False, 654 showticklabels=False, 655 showbackground=False, 656 ) 657 fig = go.Figure(data=planes_data) 658 fig.update_layout( 659 dict( 660 showlegend=True, 661 scene=dict(xaxis=axis_dict, yaxis=axis_dict, zaxis=axis_dict), 662 ) 663 ) 664 665 return fig 666 667 def _get_azimuth_elev(self, miller_index): 668 """ 669 Args: 670 miller_index: viewing direction 671 672 Returns: 673 azim, elev for plotting 674 """ 675 if miller_index in [(0, 0, 1), (0, 0, 0, 1)]: 676 return 0, 90 677 cart = self.lattice.get_cartesian_coords(miller_index) 678 azim = get_angle([cart[0], cart[1], 0], (1, 0, 0)) 679 v = [cart[0], cart[1], 0] 680 elev = get_angle(cart, v) 681 return azim, elev 682 683 @property 684 def volume(self): 685 """ 686 Volume of the Wulff shape 687 """ 688 return self.wulff_convex.volume 689 690 @property 691 def miller_area_dict(self): 692 """ 693 Returns {hkl: area_hkl on wulff} 694 """ 695 return dict(zip(self.miller_list, self.color_area)) 696 697 @property 698 def miller_energy_dict(self): 699 """ 700 Returns {hkl: surface energy_hkl} 701 """ 702 return dict(zip(self.miller_list, self.e_surf_list)) 703 704 @property 705 def surface_area(self): 706 """ 707 Total surface area of Wulff shape. 708 """ 709 return sum(self.miller_area_dict.values()) 710 711 @property 712 def weighted_surface_energy(self): 713 """ 714 Returns: 715 sum(surface_energy_hkl * area_hkl)/ sum(area_hkl) 716 """ 717 return self.total_surface_energy / self.surface_area 718 719 @property 720 def area_fraction_dict(self): 721 """ 722 Returns: 723 (dict): {hkl: area_hkl/total area on wulff} 724 """ 725 return {hkl: area / self.surface_area for hkl, area in self.miller_area_dict.items()} 726 727 @property 728 def anisotropy(self): 729 """ 730 Returns: 731 (float) Coefficient of Variation from weighted surface energy 732 The ideal sphere is 0. 733 """ 734 square_diff_energy = 0 735 weighted_energy = self.weighted_surface_energy 736 area_frac_dict = self.area_fraction_dict 737 miller_energy_dict = self.miller_energy_dict 738 739 for hkl, energy in miller_energy_dict.items(): 740 square_diff_energy += (energy - weighted_energy) ** 2 * area_frac_dict[hkl] 741 return np.sqrt(square_diff_energy) / weighted_energy 742 743 @property 744 def shape_factor(self): 745 """ 746 This is useful for determining the critical nucleus size. 747 A large shape factor indicates great anisotropy. 748 See Ballufi, R. W., Allen, S. M. & Carter, W. C. Kinetics 749 of Materials. (John Wiley & Sons, 2005), p.461 750 751 Returns: 752 (float) Shape factor. 753 """ 754 return self.surface_area / (self.volume ** (2 / 3)) 755 756 @property 757 def effective_radius(self): 758 """ 759 Radius of the Wulffshape when the 760 Wulffshape is approximated as a sphere. 761 762 Returns: 763 (float) radius. 764 """ 765 return ((3 / 4) * (self.volume / np.pi)) ** (1 / 3) 766 767 @property 768 def total_surface_energy(self): 769 """ 770 Total surface energy of the Wulff shape. 771 772 Returns: 773 (float) sum(surface_energy_hkl * area_hkl) 774 """ 775 tot_surface_energy = 0 776 for hkl, energy in self.miller_energy_dict.items(): 777 tot_surface_energy += energy * self.miller_area_dict[hkl] 778 return tot_surface_energy 779 780 @property 781 def tot_corner_sites(self): 782 """ 783 Returns the number of vertices in the convex hull. 784 Useful for identifying catalytically active sites. 785 """ 786 return len(self.wulff_convex.vertices) 787 788 @property 789 def tot_edges(self): 790 """ 791 Returns the number of edges in the convex hull. 792 Useful for identifying catalytically active sites. 793 """ 794 all_edges = [] 795 for facet in self.facets: 796 edges = [] 797 pt = self.get_line_in_facet(facet) 798 799 lines = [] 800 for i, p in enumerate(pt): 801 if i == len(pt) / 2: 802 break 803 lines.append(tuple(sorted(tuple([tuple(pt[i * 2]), tuple(pt[i * 2 + 1])])))) 804 805 for i, p in enumerate(lines): 806 if p not in all_edges: 807 edges.append(p) 808 809 all_edges.extend(edges) 810 811 return len(all_edges) 812