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:

  • ndarray –

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

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

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

    Returns
    -------
    np.ndarray
        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:

  • ndarray –

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

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

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

    Returns
    -------
    np.ndarray
        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:

  • ndarray –

    Matrix defining inversion, shape(3,3).

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

    Returns
    -------
    np.ndarray
        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:

  • ndarray –

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

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

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

    Returns
    -------
    np.ndarray
        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:

  • ndarray –

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

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

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

    Returns
    -------
    np.ndarray
        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 molsympy/core/sym_ops.py 
106
107
108
109
110
@njit(cache=True)
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 molsympy/core/sym_ops.py 
112
113
114
115
116
@njit(cache=True)
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:

  • ndarray –

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

Source code in molsympy/core/sym_ops.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
@njit(cache=True)
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.ndarray
        Vector of arbitrary magnitude, shape(n,).

    Returns
    -------
    np.ndarray
        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 molsympy/core/sym_ops.py 
138
139
140
141
@njit(cache=True)
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 molsympy/core/sym_ops.py 
143
144
145
146
147
148
149
@njit(cache=True)
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.ndarray Vector a, shape(n,). b: np.ndarray 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 molsympy/core/sym_ops.py 
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
@njit(cache=True)
def issame_axis(a, b, tol=global_tol):
    """
    Return True if vectors a and b are colinear within the global tolerance.

    Paremeters
    ----------
    a: np.ndarray
        Vector a, shape(n,).
    b: np.ndarray
        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 == 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 molsympy/core/sym_ops.py 
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
@njit(cache=True)
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 molsympy/core/sym_ops.py 
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
@njit(cache=True)
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 molsympy/core/sym_ops.py 
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
@njit(cache=True)
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 molsympy/core/sym_ops.py 
239
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
@njit(cache=True)
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 molsympy/core/sym_ops.py 
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
@njit(cache=True)
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: ndarray, saxis: ndarray, 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 molsympy/core/sym_ops.py 
306
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
@njit(cache=True)
def generate_cyclic_axes(static_axis:np.ndarray, saxis:np.ndarray, 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

axes_are_orthogonal 

axes_are_orthogonal(axes: ndarray, tol: float = NUMERICAL_TOL) -> bool

Return True if axes are all orthogonal to each other.

Source code in molsympy/core/sym_ops.py 
391
392
393
394
395
396
397
398
399
400
401
402
@njit(cache=True)
def axes_are_orthogonal(axes: np.ndarray, tol:float=global_tol) -> bool:
    """Return True if axes are all orthogonal to each other."""
    naxes = len(axes)
    axes = np.ascontiguousarray(axes)
    for ii in range(naxes):
        ii_axis = axes[ii]
        for jj in range(ii+1, naxes):
            d = np.abs(np.dot(ii_axis, axes[jj]))
            if float_isclose(d, 0.0, atol=tol) == False:
                return False
    return True

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.