Skip to content

shape_descriptors

make_N_invariants(coefficients)

Construct the N type invariants from SHT coefficients. If coefficients is of length n, the size of the result will be sqrt(n)

Parameters:

Name Type Description Default
coefficients ndarray

the set of spherical harmonic coefficients

 required

Returns:

Type Description
ndarray

np.ndarray the N type rotational invariants based on these coefficients
Source code in chmpy/shape/shape_descriptors.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
def make_N_invariants(coefficients) -> np.ndarray:
    """
    Construct the `N` type invariants from SHT coefficients.
    If coefficients is of length n, the size of the result will be sqrt(n)

    Arguments:
        coefficients (np.ndarray): the set of spherical harmonic coefficients

    Returns:
        np.ndarray the `N` type rotational invariants based on these coefficients
    """
    size = int(np.sqrt(len(coefficients)))
    invariants = np.empty(shape=(size), dtype=np.float64)
    for i in range(0, size):
        lower, upper = i ** 2, (i + 1) ** 2
        invariants[i] = np.sum(
            coefficients[lower : upper + 1]
            * np.conj(coefficients[lower : upper + 1])
        ).real
    return np.sqrt(invariants)

make_invariants(l_max, coefficients, kinds='NP')

Construct the N and/or P type invariants from SHT coefficients.

Parameters:

Name Type Description Default
l_max int

the maximum angular momentum of the coefficients

 required
coefficients ndarray

the set of spherical harmonic coefficients

 required
kinds str

which kinds of invariants to include

 'NP'

Returns:

Type Description
ndarray

np.ndarray the N and/or P type rotational invariants based on these coefficients
Source code in chmpy/shape/shape_descriptors.py 
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
def make_invariants(l_max, coefficients, kinds="NP") -> np.ndarray:
    """
    Construct the `N` and/or `P` type invariants from SHT coefficients.

    Arguments:
        l_max (int): the maximum angular momentum of the coefficients
        coefficients (np.ndarray): the set of spherical harmonic coefficients
        kinds (str, optional): which kinds of invariants to include

    Returns:
        np.ndarray the `N` and/or `P` type rotational invariants based on these coefficients
    """

    global _HAVE_WARNED_ABOUT_LMAX_P
    invariants = []
    if "N" in kinds:
        invariants.append(make_N_invariants(coefficients))
    if "P" in kinds:
        # Because we only have factorial precision (double precision)
        # in our clebsch implementation up to 70! l_max for P type
        # invariants is restricted to <= 23
        # TODO use a better clebsch gordan coefficients implementation
        # e.g. that in https://github.com/GXhelsinki/Clebsch-Gordan-Coefficients-
        pfunc = p_invariants_c
        MAX_L_MAX = 23
        if l_max > MAX_L_MAX:
            if not _HAVE_WARNED_ABOUT_LMAX_P:
                LOG.warn(
                    f"P type invariants only supported up to l_max = {MAX_L_MAX}: "
                    "will only using N type invariants beyond that."
                )
                _HAVE_WARNED_ABOUT_LMAX_P = True
            c = coefficients[: MAX_L_MAX * MAX_L_MAX]
            invariants.append(pfunc(c))
        else:
            invariants.append(pfunc(coefficients))
    return np.hstack(invariants)

promolecule_density_descriptor(sht, n_i, p_i, **kwargs)

Calculate the shape description of the promolecule density isosurface.

Parameters:

Name Type Description Default
sht SHT

the spherical harmonic transform object handle

 required
n_i ndarray

atomic numbers of the atoms

 required
p_i ndarray

Cartesian coordinates of the atoms

 required
**kwargs

keyword arguments for optional settings.

 {}

Keyword Args: isovalue (float): change the Hirshfeld weight value (default 0.5) with_property (str): calculate the combined shape + surface property descriptor using the specified property on the surface (e.g. d_norm, esp) bounds (Tuple): modify the lower/upper bounds on the search for the isovalue (default 0.1, 20.0) coefficients (bool): also return the coefficients of the SHT origin (np.ndarray): specify the center of the surface (default is the geometric centroid of the atoms) kinds (str): the kinds of invariants to calculate (default 'NP') Returns: np.ndarray: the rotation invariant descriptors of the promolecule surface shape

Source code in chmpy/shape/shape_descriptors.py 
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
def promolecule_density_descriptor(sht, n_i, p_i, **kwargs):
    """
    Calculate the shape description of the promolecule density isosurface.

    Args:
        sht (SHT): the spherical harmonic transform object handle
        n_i (np.ndarray): atomic numbers of the atoms
        p_i (np.ndarray): Cartesian coordinates of the atoms
        **kwargs: keyword arguments for optional settings.
    Keyword Args:
        isovalue (float): change the Hirshfeld weight value (default 0.5)
        with_property (str): calculate the combined shape + surface
            property descriptor using the specified property on the
            surface (e.g. d_norm, esp)
        bounds (Tuple): modify the lower/upper bounds on the search for
            the isovalue (default 0.1, 20.0)
        coefficients (bool): also return the coefficients of the SHT
        origin (np.ndarray): specify the center of the surface
            (default is the geometric centroid of the atoms)
        kinds (str): the kinds of invariants to calculate (default 'NP')
    Returns:
        np.ndarray: the rotation invariant descriptors of the promolecule surface shape
    """
    isovalue = kwargs.get("isovalue", 0.0002)
    property_function = kwargs.get("with_property", None)
    r_min, r_max = kwargs.get("bounds", (0.4, 20.0))
    pro = PromoleculeDensity((n_i, p_i))
    g = np.empty((sht.grid[0].size, 3), dtype=np.float32)
    x, y, z = sht.grid_cartesian
    g[:, 0] = x.flatten()
    g[:, 1] = y.flatten()
    g[:, 2] = z.flatten()

    o = kwargs.get("origin", np.mean(p_i, axis=0, dtype=np.float32))
    r = sphere_promolecule_radii(pro.dens, o, g, r_min, r_max, 1e-12, 30, isovalue).reshape(sht.grid[0].shape)
    if np.any(r < 0):
        raise ValueError(
            f"Unable to find isovalue {isovalue:.2f} in all directions for bounds ({r_min:.2f}, {r_max:.2f})"
        )
    real = True
    if property_function is not None:
        if property_function == "d_norm":
            property_function = lambda x: pro.d_norm(x)[1]
        elif property_function == "esp":
            from chmpy import Molecule

            els = pro.elements
            pos = pro.positions
            property_function = Molecule.from_arrays(els, pos).electrostatic_potential
        r = _compute_property_in_j_channel(sht, r, property_function)
        real = False
    l_max = sht.lmax
    coeffs = sht.analysis(r)
    coeff4inv = expand_coeffs_to_full(l_max, coeffs) if real else coeffs
    invariants = make_invariants(
        l_max, coeff4inv, kinds=kwargs.get("kinds", "NP")
    )

    if kwargs.get("coefficients", False):
        return coeffs, invariants
    return invariants

stockholder_weight_descriptor(sht, n_i, p_i, n_e, p_e, **kwargs)

Calculate the 'stockholder weight' shape descriptors based on the Hirshfeld weight i.e. ratio of electron density from the 'interior' to the total electron density.

Parameters:

Name Type Description Default
sht SHT

the spherical harmonic transform object handle

 required
n_i ndarray

atomic numbers of the interior atoms

 required
p_i ndarray

Cartesian coordinates of the interior atoms

 required
n_e ndarray

atomic numbers of the exterior atoms

 required
p_e ndarray

Cartesian coordinates of the exterior atoms

 required

Other Parameters:

Name Type Description
isovalue float

change the Hirshfeld weight value (default 0.5)
background float

include an optional 'background' electron density (default 0.0)
with_property str

calculate the combined shape + surface property descriptor using the specified property on the surface (e.g. d_norm, esp)
bounds Tuple

modify the lower/upper bounds on the search for the isovalue (default 0.1, 20.0)
coefficients bool

also return the coefficients of the SHT
origin ndarray

specify the center of the surface (default is the geometric centroid of the interior atoms)
kinds str

the kinds of invariants to calculate (default 'NP')

Returns: np.ndarray: the rotation invariant descriptors of the Hirshfeld surface shape

Source code in chmpy/shape/shape_descriptors.py 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 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
def stockholder_weight_descriptor(sht, n_i, p_i, n_e, p_e, **kwargs):
    """
    Calculate the 'stockholder weight' shape descriptors based on the
    Hirshfeld weight i.e. ratio of electron density from the 'interior'
    to the total electron density.

    Args:
        sht (SHT): the spherical harmonic transform object handle
        n_i (np.ndarray): atomic numbers of the interior atoms
        p_i (np.ndarray): Cartesian coordinates of the interior atoms
        n_e (np.ndarray): atomic numbers of the exterior atoms
        p_e (np.ndarray): Cartesian coordinates of the exterior atoms

    Keyword Args:
        isovalue (float): change the Hirshfeld weight value (default 0.5)
        background (float): include an optional 'background' electron
            density (default 0.0)
        with_property (str): calculate the combined shape + surface
            property descriptor using the specified property on the
            surface (e.g. d_norm, esp)
        bounds (Tuple): modify the lower/upper bounds on the search for
            the isovalue (default 0.1, 20.0)
        coefficients (bool): also return the coefficients of the SHT
        origin (np.ndarray): specify the center of the surface
            (default is the geometric centroid of the interior atoms)
        kinds (str): the kinds of invariants to calculate (default 'NP')
    Returns:
        np.ndarray: the rotation invariant descriptors of the Hirshfeld surface shape
    """
    isovalue = kwargs.get("isovalue", 0.5)
    background = kwargs.get("background", 0.0)
    property_function = kwargs.get("with_property", None)
    r_min, r_max = kwargs.get("bounds", (0.1, 20.0))
    o = kwargs.get("origin", np.mean(p_i, axis=0, dtype=np.float32))
    s = StockholderWeight.from_arrays(n_i, p_i, n_e, p_e, background=background)
    g = np.empty((sht.grid[0].size, 3), dtype=np.float32)
    x, y, z = sht.grid_cartesian
    g[:, 0] = x.flatten()
    g[:, 1] = y.flatten()
    g[:, 2] = z.flatten()

    r = sphere_stockholder_radii(s.s, o, g, r_min, r_max, 1e-7, 30, isovalue).reshape(sht.grid[0].shape)
    if np.any(r < 0):
        raise ValueError(
            f"Unable to find isovalue {isovalue:.2f} in all directions for bounds ({r_min:.2f}, {r_max:.2f})"
        )
    real = True
    if property_function is not None:
        if property_function == "d_norm":
            property_function = lambda x: s.d_norm(x)[3]
        elif property_function == "esp":
            from chmpy import Molecule

            els = s.dens_a.elements
            pos = s.dens_a.positions
            property_function = Molecule.from_arrays(
                s.dens_a.elements, s.dens_a.positions
            ).electrostatic_potential
        r = _compute_property_in_j_channel(sht, r, property_function, origin=o)
        real = False
    l_max = sht.lmax
    coeffs = sht.analysis(r)

    coeff4inv = expand_coeffs_to_full(l_max, coeffs) if real else coeffs
    invariants = make_invariants(
        l_max, coeff4inv, kinds=kwargs.get("kinds", "NP")
    )

    if kwargs.get("coefficients", False):
        return coeffs, invariants
    return invariants