Skip to content

wulff

WulffSHT

Source code in chmpy/crystal/wulff.py 
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
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 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
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 
88
89
90
91
92
93
94
95
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 
45
46
47
48
49
50
51
def winding_order_ccw(points):
    "return the indices to reorder the provided 2D points into CCW order"
    centroid = points[0]
    directions = points[1:] - centroid
    directions /= np.linalg.norm(directions, axis=1)[:, np.newaxis]
    idxs = list(range(1, points.shape[0]))
    return [0] + sorted(idxs, key=lambda x: np.arctan2(directions[x - 1, 1], directions[x - 1, 0]))