Skip to content

wulff

WulffSHT

Source code in chmpy/crystal/wulff.py
class WulffSHT:
    def __init__(self, facet_normals, facet_energies, labels=None,
                 l_max=10, sht_object=None, scale=1):
        from chmpy.shape.sht import SHT


        self.facet_normals = facet_normals
        self.facet_energies = facet_energies
        self.labels = labels
        self.l_max = l_max
        if self.labels is None:
            self.labels = np.arange(len(facet_normals))

        if sht_object is None:
            sht = SHT(l_max)
        else:
            sht = sht_object

        self.grid = np.array(sht.grid_cartesian)
        dps = np.einsum('ijk,li->jkl', self.grid, self.facet_normals)
        mask = (dps > 0)
        intersections = np.inf * np.ones_like(dps)
        intersections[mask] = 1.0 / dps[mask]
        intersections *= scale * self.facet_energies

        facet_idx = np.argmin(intersections, axis=2, keepdims=True)
        self.distances = np.take_along_axis(intersections, facet_idx, axis=2).squeeze()
        self.facet_idx = facet_idx.squeeze()
        self.coeffs = sht.analysis(self.distances)

    def invariants(self, **kwargs):
        from chmpy.shape.shape_descriptors import expand_coeffs_to_full, make_invariants
        coeff4inv = expand_coeffs_to_full(self.l_max, self.coeffs)
        invariants = make_invariants(
            self.l_max, coeff4inv, kinds=kwargs.get("kinds", "NP")
        )
        return invariants

    def invariants_kazhdan(self, sht_object=None, **kwargs):
        from chmpy.shape.sht import SHT
        if sht_object is None:
            sht = SHT(self.l_max)
        else:
            sht = sht_object
        return sht.invariants_pure_python(self.coeffs)

    def power_spectrum(self, sht_object=None, **kwargs):
        from chmpy.shape.sht import SHT
        if sht_object is None:
            sht = SHT(self.l_max)
        else:
            sht = sht_object
        return sht.power_spectrum(self.coeffs)

    @classmethod
    def from_gmf_and_crystal(cls, gmf, crystal, **kwargs):
        hkl, energies = expand_symmetry_related_planes(gmf.hkl, gmf.energies, crystal.space_group)
        facet_normals = hkl @ crystal.uc.reciprocal_lattice
        facet_normals /= np.linalg.norm(facet_normals, axis=1)[:, None]
        return cls(facet_normals, energies, labels=energies, **kwargs)

    def __repr__(self):
        return f"WulffSHT<l={self.l_max}>"

__init__(facet_normals, facet_energies, labels=None, l_max=10, sht_object=None, scale=1)

Source code in chmpy/crystal/wulff.py
def __init__(self, facet_normals, facet_energies, labels=None,
             l_max=10, sht_object=None, scale=1):
    from chmpy.shape.sht import SHT


    self.facet_normals = facet_normals
    self.facet_energies = facet_energies
    self.labels = labels
    self.l_max = l_max
    if self.labels is None:
        self.labels = np.arange(len(facet_normals))

    if sht_object is None:
        sht = SHT(l_max)
    else:
        sht = sht_object

    self.grid = np.array(sht.grid_cartesian)
    dps = np.einsum('ijk,li->jkl', self.grid, self.facet_normals)
    mask = (dps > 0)
    intersections = np.inf * np.ones_like(dps)
    intersections[mask] = 1.0 / dps[mask]
    intersections *= scale * self.facet_energies

    facet_idx = np.argmin(intersections, axis=2, keepdims=True)
    self.distances = np.take_along_axis(intersections, facet_idx, axis=2).squeeze()
    self.facet_idx = facet_idx.squeeze()
    self.coeffs = sht.analysis(self.distances)

perpendicular_vector(p, q, r)

a unit vector perpendicular to the triangle p q r

Source code in chmpy/crystal/wulff.py
def perpendicular_vector(p, q, r):
    "a unit vector perpendicular to the triangle p q r"
    perp_vector = np.cross(q - p, r - p)
    dp = np.dot(perp_vector, p)
    if dp > 0:
        return perp_vector / np.linalg.norm(perp_vector)
    else:
        return - perp_vector / np.linalg.norm(perp_vector)

winding_order_ccw(points)

return the indices to reorder the provided 2D points into CCW order

Source code in chmpy/crystal/wulff.py
def winding_order_ccw(points):
    "return the indices to reorder the provided 2D points into CCW order"
    centroid = np.mean(points, axis=0)
    directions = points - centroid
    directions /= np.linalg.norm(directions, axis=1)[:, np.newaxis]
    idxs = list(range(points.shape[0]))
    return sorted(idxs, key=lambda x: np.arctan2(directions[x, 1], directions[x, 0]))