Skip to content

Utilities

sym_ops

Referenced by:

rotation_matrix 

rotation_matrix(axis, theta)

Create rotation matrix about an axis by theta in radians.

Parameters:

  • axis –

    Cartesian vector defining rotation axis, shape(3,).

  • theta –

    Angle of rotation in radians.

Returns:

  • array –

    Matrix defining rotation on column vector, shape (3,3).

Source code in minimalsym/core/sym_ops.py 
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
@njit
def rotation_matrix(axis, theta):
    """
    Create rotation matrix about an axis by theta in radians.

    Parameters
    ----------
    axis: np.array
        Cartesian vector defining rotation axis, shape(3,).
    theta: float
        Angle of rotation in radians.

    Returns
    -------
    np.array
        Matrix defining rotation on column vector, shape (3,3).
    """
    kmat1 = np.array([[0.0, -axis[2], axis[1]],
                      [axis[2], 0.0, -axis[0]],
                      [-axis[1], axis[0], 0.0]])
    kmat2 = kmat1 @ kmat1
    rodriguesrm = np.eye(3) + np.sin(theta)*kmat1 + (1.0 - np.cos(theta))*kmat2
    return rodriguesrm

reflection_matrix 

reflection_matrix(axis)

Create reflection matrix about a plane defined by its normal vector.

Parameters:

  • axis –

    Cartesian vector defining the plane normal vector, shape(3,).

Returns:

  • array –

    Matrix defining reflection on column vector, shape (3,3).

Source code in minimalsym/core/sym_ops.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
@njit
def reflection_matrix(axis):
    """
    Create reflection matrix about a plane defined by its normal vector.

    Parameters
    ----------
    axis: np.array
        Cartesian vector defining the plane normal vector, shape(3,).

    Returns
    -------
    np.array
        Matrix defining reflection on column vector, shape (3,3).
    """
    M = np.zeros((3, 3))
    for i in range(3):
        for j in range(i, 3):
            if i == j:
                M[i,i] = 1 - 2*(axis[i]**2)
            else:
                M[i,j] = -2 * axis[i] * axis[j]
                M[j,i] = M[i,j]
    return M

inversion_matrix 

inversion_matrix()

Create cartesian inversion matrix.

Returns:

  • array –

    Matrix defining inversion, shape(3,3).

Source code in minimalsym/core/sym_ops.py 
54
55
56
57
58
59
60
61
62
63
64
@njit
def inversion_matrix():
    """
    Create cartesian inversion matrix.

    Returns
    -------
    np.array
        Matrix defining inversion, shape(3,3).
    """
    return -1*np.eye(3)

Cn 

Cn(axis, n)

Wrapper around rotation_matrix for producing a C_n rotation about axis.

Parameters:

  • axis –

    Cartesian vector defining rotation axis, shape(3,).

  • n –

    Defines rotation angle by theta = 2 pi / n.

Returns:

  • array –

    Matrix defining proper rotation on column vector, shape(3,3).

Source code in minimalsym/core/sym_ops.py 
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
@njit
def Cn(axis, n):
    """
    Wrapper around rotation_matrix for producing a C_n rotation about axis.

    Parameters
    ----------
    axis: np.array
        Cartesian vector defining rotation axis, shape(3,).
    n: int
        Defines rotation angle by theta = 2 pi / n.

    Returns
    -------
    np.array
        Matrix defining proper rotation on column vector, shape(3,3).
    """
    theta = 2*np.pi/n
    return rotation_matrix(axis, theta)

Sn 

Sn(axis, n)

Improper rotation S_n about an axis.

Parameters:

  • axis –

    Cartesian vector defining rotation axis, shape(3,).

  • n –

    Defines rotation angle by theta = 2 pi / n.

Returns:

  • array –

    Matrix defining improper rotation on column vector, shape (3,3).

Source code in minimalsym/core/sym_ops.py 
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
@njit
def Sn(axis, n):
    """
    Improper rotation S_n about an axis.

    Parameters
    ----------
    axis: np.array
        Cartesian vector defining rotation axis, shape(3,).
    n: int
        Defines rotation angle by theta = 2 pi / n.

    Returns
    -------
    np.array
        Matrix defining improper rotation on column vector, shape (3,3).
    """
    return np.dot(reflection_matrix(axis), Cn(axis, n))

vec_norm_axis 

vec_norm_axis(v)

Compute Euclidean norm (Numba-compatible).

Source code in minimalsym/core/sym_ops.py 
105
106
107
108
109
@njit
def vec_norm_axis(v):
    """Compute Euclidean norm (Numba-compatible)."""
    # Vectorize the norm calculation over rows
    return np.sqrt(np.sum(v**2, axis=1))

vec_norm 

vec_norm(v)

Compute Euclidean norm (Numba-compatible).

Source code in minimalsym/core/sym_ops.py 
111
112
113
114
115
@njit
def vec_norm(v):
    """Compute Euclidean norm (Numba-compatible)."""
    #return np.sqrt(np.sum(v**2))
    return np.sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2])

normalize 

normalize(a)

Normalize vector a to unit length, return zero vector shape(n, ) if the input vector is of zero length.

Parameters:

  • a –

    Vector of arbitrary magnitude, shape(n,).

Returns:

  • array –

    Normalized vector shape(n,) or zero vector shape(n, ) if the magnitude of a is less than the global tolerance.

Source code in minimalsym/core/sym_ops.py 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
@njit
def normalize(a):
    """
    Normalize vector a to unit length, return zero vector shape(n, ) if the input vector is of zero length.

    Parameters
    ----------
    a: np.array
        Vector of arbitrary magnitude, shape(n,).

    Returns
    -------
    np.array
        Normalized vector shape(n,) or zero vector shape(n, ) if the magnitude of ``a`` is less than the global tolerance.
    """
    n = vec_norm(a)
    if n <= global_tol:
        return np.zeros(3)
    return a / n

float_isclose 

float_isclose(a: float, b: float, rtol: float = 1e-05, atol: float = 1e-08) -> bool

Compare two floats (Numba-compatible).

Source code in minimalsym/core/sym_ops.py 
137
138
139
140
@njit
def float_isclose(a: float, b: float, rtol:float=1e-05, atol:float=1e-08) -> bool:
    """Compare two floats (Numba-compatible)."""
    return np.abs(a - b) <= (atol + rtol * max(np.abs(b), np.abs(a)))

inertia_isclose 

inertia_isclose(a: float, b: float, rtol: float = 1e-05, atol: float = 1e-08) -> bool

Compare two inertia (Numba-compatible).

Source code in minimalsym/core/sym_ops.py 
142
143
144
145
146
147
148
@njit
def inertia_isclose(a: float, b: float, rtol:float=1e-05, atol:float=1e-08) -> bool:
    """Compare two inertia (Numba-compatible)."""
    if a == 0.0 or b == 0.0:
        return float_isclose(a, b, rtol=0.0, atol=atol)
    else:
        return float_isclose(a, b, rtol=rtol, atol=0.0)

issame_axis 

issame_axis(a, b, tol=NUMERICAL_TOL)

Return True if vectors a and b are colinear within the global tolerance.

Paremeters

a: np.array Vector a, shape(n,). b: np.array Vector b, shape(n,). tol: float Tolerance for error, default is global_tol.

Returns:

  • bool –

    True if vectors are collinear, False if not collinear or if either vector has zero length.

Source code in minimalsym/core/sym_ops.py 
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
@njit
def issame_axis(a, b, tol=global_tol):
    """
    Return True if vectors a and b are colinear within the global tolerance.

    Paremeters
    ----------
    a: np.array
        Vector a, shape(n,).
    b: np.array
        Vector b, shape(n,).
    tol: float
        Tolerance for error, default is ``global_tol``.

    Returns
    -------
    bool
        True if vectors are collinear, False if not collinear or if either vector has zero length.
    """
    A_vector = normalize(a)
    B_vector = normalize(b)
    if A_vector is None or B_vector is None:
        return False
    if (A_vector == np.zeros(3)).all() or (B_vector == np.zeros(3)).all():
        return False
    d = np.abs(np.dot(A_vector, B_vector))
    return float_isclose(d, 1.0, atol=tol)

isfactor 

isfactor(n, a)

Return True if a divides n.

Parameters:

  • n –

    Dividend.

  • a –

    Divisor.

Returns:

  • bool –

    True if a divides n with remainder 0.

Source code in minimalsym/core/sym_ops.py 
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
@njit
def isfactor(n, a):
    """
    Return True if a divides n.

    Parameters
    ----------
    n: int
        Dividend.
    a: int
        Divisor.

    Returns
    -------
    bool
        True if ``a`` divides ``n`` with remainder 0.
    """
    return n % a == 0

reduce 

reduce(n, i)

Divide n and i by their greatest common divisor g.

Parameters:

  • n –
  • i –

Returns:

  • tuple –

    Tuple of n/g and i/g.

Source code in minimalsym/core/sym_ops.py 
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
@njit
def reduce(n, i):
    """
    Divide n and i by their greatest common divisor g.

    Parameters
    ----------
    n: int
    i: int

    Returns
    -------
    tuple
        Tuple of n/g and i/g.
    """
    g = _gcd(n, i)
    return n//g, i//g

_gcd 

_gcd(A, B)

Euclid algorithm for finding the greatest common divisor between A and B.

Parameters:

  • A –
  • B –

Returns:

  • int –

    Greatest common divisor between A and B.

Source code in minimalsym/core/sym_ops.py 
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
@njit
def _gcd(A, B):
    """
    Euclid algorithm for finding the greatest common divisor between A and B.

    Parameters
    ----------
    A: int
    B: int

    Returns
    -------
    int
        Greatest common divisor between A and B.
    """
    a = max(A, B)
    b = min(A, B)
    if a == 0:
        return b
    elif b == 0:
        return a
    else:
        r = a % b
        return _gcd(b, r)

unique_sorted 

unique_sorted(arr)

Return the sorted unique elements of a 1D array.

Parameters:

  • arr (ndarray) –

    Input 1D array of comparable elements.

Returns:

  • ndarray –

    Sorted array containing the unique elements of arr.

Notes

This is a Numba-compatible alternative to np.unique, implemented using sorting followed by linear duplicate removal.

Equality is tested using exact comparison (!=), so this function is best suited for integer or discretized floating-point data.

Source code in minimalsym/core/sym_ops.py 
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
@njit
def unique_sorted(arr):
    """
    Return the sorted unique elements of a 1D array.

    Parameters
    ----------
    arr : np.ndarray
        Input 1D array of comparable elements.

    Returns
    -------
    np.ndarray
        Sorted array containing the unique elements of `arr`.

    Notes
    -----
    This is a Numba-compatible alternative to `np.unique`, implemented
    using sorting followed by linear duplicate removal.

    Equality is tested using exact comparison (`!=`), so this function
    is best suited for integer or discretized floating-point data.
    """
    if len(arr) == 0:
        return arr

    sorted_arr = np.sort(arr)

    # allocate output (max possible size)
    out = np.empty_like(sorted_arr)
    count = 0

    out[count] = sorted_arr[0]
    count += 1

    for i in range(1, len(sorted_arr)):
        if sorted_arr[i] != out[count-1]:
            out[count] = sorted_arr[i]
            count += 1

    return out[:count]

canonical 

canonical(v)

Returns canonical form of vector. Ensures that the first non-zero element is positive.

Parameters:

  • v ((ndarray, shape(3))) –

    Vector to canonicalize

Returns:

  • (ndarray, shape(3)) –

    Canonical vector

Source code in minimalsym/core/sym_ops.py 
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
@njit
def canonical(v):
    """
    Returns canonical form of vector.
    Ensures that the first non-zero element is positive.

    Parameters
    ----------
    v : np.ndarray, shape (3,)
        Vector to canonicalize

    Returns
    -------
    np.ndarray, shape (3,)
        Canonical vector
    """
    v = normalize(v)
    for i in range(len(v)):
        if abs(v[i]) < 1e-08:
            continue
        if v[i] < 0:
            v = -v
        break
    return v

generate_cyclic_axes 

generate_cyclic_axes(static_axis: array, saxis: array, n: int, num_elem_generate: int = -1)

Generate a set of axes by rotating a reference axis around a fixed axis.

Parameters:

  • static_axis ((ndarray, shape(3))) –

    Unit vector defining the rotation axis (e.g., principal symmetry axis).

  • saxis ((ndarray, shape(3))) –

    Reference axis to be rotated about static_axis.

  • n (int) –

    Order of the rotation symmetry (Cn). Defines the angular step 2π/n. If n is even, it is internally doubled to ensure full coverage of distinct orientations.

  • num_elem_generate (int, default: -1 ) –

    Number of rotated axes to generate. If <= 0, defaults to n.

Returns:

  • (ndarray, shape(num_elem_generate, 3)) –

    Array of rotated axes. The first row is saxis, and subsequent rows are obtained by successive rotations about static_axis.

Notes

This function is typically used to generate symmetry-equivalent axes (e.g., C2 axes perpendicular to a principal axis in dihedral groups).

Source code in minimalsym/core/sym_ops.py 
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
@njit
def generate_cyclic_axes(static_axis:np.array, saxis:np.array, n:int, num_elem_generate:int = -1):
    """
    Generate a set of axes by rotating a reference axis around a fixed axis.

    Parameters
    ----------
    static_axis : np.ndarray, shape (3,)
        Unit vector defining the rotation axis (e.g., principal symmetry axis).
    saxis : np.ndarray, shape (3,)
        Reference axis to be rotated about `static_axis`.
    n : int
        Order of the rotation symmetry (Cn). Defines the angular step 2π/n.
        If `n` is even, it is internally doubled to ensure full coverage of
        distinct orientations.
    num_elem_generate : int, optional
        Number of rotated axes to generate. If <= 0, defaults to `n`.

    Returns
    -------
    np.ndarray, shape (num_elem_generate, 3)
        Array of rotated axes. The first row is `saxis`, and subsequent rows
        are obtained by successive rotations about `static_axis`.

    Notes
    -----
    This function is typically used to generate symmetry-equivalent axes
    (e.g., C2 axes perpendicular to a principal axis in dihedral groups).
    """
    if num_elem_generate <= 0:
        num_elem_generate = n
    if n % 2 == 0:
        n *= 2
    rotated_axes = np.empty((num_elem_generate, 3), dtype=np.float64)
    rotated_axes[0, :] = saxis

    for ii in range(1, num_elem_generate):
        theta = 2 * np.pi * ii / n
        R = rotation_matrix(static_axis, theta)
        rotated_axis = R @ saxis
        #rotated_axis = rotated_axis - np.dot(rotated_axis, static_axis) * static_axis
        rotated_axes[ii, :] = normalize(rotated_axis)

    return rotated_axes

constants

Shared numerical constants for the core package.

Referenced by:

NUMERICAL_TOL module-attribute 

NUMERICAL_TOL: float = 1e-08

Internal floating-point comparison tolerance used across mol_ops and sym_ops. Governs axis collinearity checks, zero-vector detection, etc.

Used by:

IH_C2_C3_ANGLE module-attribute 

IH_C2_C3_ANGLE: float = 0.36486382647383764

Angle (radians) between a C2 axis and the nearest C3 axis in an icosahedron. Derived from the golden ratio phi = (1+sqrt(5))/2: theta = arccos(phi^2 / sqrt(1 + phi^4)) = arctan(1/phi^2) ~= 0.3649 rad ~= 20.9 deg Used in find_point_group to identify the secondary axis for Ih/I molecules.

IH_ANGLE_TOL module-attribute 

IH_ANGLE_TOL: float = 0.0001

Angular tolerance for the IH_C2_C3_ANGLE comparison.

PRINT_WARNINGS module-attribute 

PRINT_WARNINGS: bool = False

If True warnings and debug messages are printed.