Skip to content

Symmetry Detection

pg_detect

pg_detect.py — Point-group detection.

This module exposes a single public function, find_point_group(), which implements the Beruski & Vidal (2013) decision-tree algorithm.

The megafunction is decomposed into three private classifiers: _classify_linear() — Ia ~= 0 (linear molecules) _classify_spherical_top() — Ia ~= Ib ~= Ic (cubic/icosahedral) _classify_symmetric_top() — two equal MOIT eigenvalues

Internal implementation is split across three focused sub-modules: rotation_detection.py — Cn axis searches reflection_detection.py — sigma plane searches special_geometry.py — icosahedral and octahedral geometry

Referenced by:

PointGroupResult dataclass 

PointGroupResult(pg: str, paxis: ndarray, saxis: ndarray)

Return value of find_point_group().

Attributes:

  • pg (str) –

    Schoenflies point-group symbol (e.g. "C2v", "D3h", "Oh").

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

    Principal symmetry axis in the original molecule frame. Zero vector when not applicable (e.g. C1).

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

    Secondary axis defining the canonical orientation. Zero vector when not applicable.

Returned by: Used by:

_classify_linear 

_classify_linear(positions, masses, geom_tol)

Classify a linear molecule (Ia ~= 0). Returns PointGroupResult.

Source code in molsympy/core/pg_detect.py 
65
66
67
68
69
def _classify_linear(positions, masses, geom_tol):
    """Classify a linear molecule (Ia ~= 0). Returns PointGroupResult."""
    paxis = _linear_mol_axis(positions)
    pg = "D0h" if transform_isequivalent(positions, masses, geom_tol, inversion_matrix()) else "C0v"
    return PointGroupResult(pg=pg, paxis=paxis, saxis=np.zeros(3))

_classify_spherical_top 

_classify_spherical_top(mol, positions, masses, geom_tol)

Classify a spherical top (Ia ~= Ib ~= Ic). Discriminates T/O/I families by the number of distinct C2 axes. Returns PointGroupResult.

Source code in molsympy/core/pg_detect.py 
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 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
158
159
160
161
162
163
164
165
def _classify_spherical_top(mol, positions, masses, geom_tol):
    """
    Classify a spherical top (Ia ~= Ib ~= Ic).
    Discriminates T/O/I families by the number of distinct C2 axes.
    Returns PointGroupResult.
    """
    seas = find_SEAs(mol)
    list_sea_subset = to_typed_list([sea.subset for sea in seas])
    c2_axes = _num_C2(positions, masses, geom_tol, list_sea_subset, 15)
    c2_axes = np.array(c2_axes)
    n = len(c2_axes)
    if n == 0:
        logger.warning("Molecule was wrongly classified as a spherical top, (num_C2 couldn't find any C2 axis), probably due to high eigen_tol. " \
            "Process will continue as general symmetry.")
        return _classify_general(mol, positions, masses, geom_tol)
    invertable = transform_isequivalent(positions, masses, geom_tol, inversion_matrix())
    found_is_spherical = False

    if n >= 15:
        # Icosahedral: paxis = C5 axis, saxis = C2 axis from golden-ratio geometry.
        try:
            c2_axis = c2_axes[0]
            c3s = _find_C3s_for_Ih(mol, seas)
            saxis = np.zeros(3)
            for c3 in c3s:
                if np.isclose(np.arccos(abs(np.dot(c3, c2_axis))), IH_C2_C3_ANGLE, atol=IH_ANGLE_TOL):
                    taxis = normalize(np.cross(c3, c2_axis))
                    saxis = normalize(np.cross(taxis, c2_axis))
                    break
            if (saxis != np.zeros(3)).all():
                phi = (1 + np.sqrt(5.0)) / 2
                theta = np.arccos(phi / np.sqrt(1 + phi**2))
                paxis = np.dot(rotation_matrix(saxis, theta), c2_axis)
                if validate_I(positions, masses, geom_tol, paxis, saxis):
                    pg = "Ih" if invertable else "I"
                    found_is_spherical = True
        except RuntimeError:
            pass
    if n >= 9 and found_is_spherical == False:
        # Octahedral: paxis and saxis are two orthogonal C4 axes.
        try:
            #c4s = #_find_C4s_for_Oh(mol, seas)
            all_c4s = []
            for c2_axis in c2_axes:
                if validate_cn_subrotations(4, c2_axis, positions, masses, geom_tol):
                    all_c4s.append(c2_axis)
            all_c4s = np.array(all_c4s)
            ortho_c4_groups = get_orthogonal_groups(all_c4s, 1e-3)
            for c4s  in ortho_c4_groups:
                c4s = orthonalize_3(c4s[0], c4s[1], c4s[2])

                if validate_O(positions, masses, geom_tol, c4s):
                    paxis, saxis = c4s[0], c4s[1]
                    pg = "Oh" if invertable else "O"
                found_is_spherical = True
        except RuntimeError:
            pass

    if n >= 3 and found_is_spherical == False:
        ortho_c2_groups = get_orthogonal_groups(c2_axes, 1e-3)
        for c2_group in ortho_c2_groups:
            try:
                ortho_c2_axes = orthonalize_3(c2_group[0], c2_group[1], c2_group[2])

                if validate_T(positions, masses, geom_tol, ortho_c2_axes) == False:
                    continue

                # Tetrahedral (n == 3): use two of the three C2 axes.
                paxis, saxis, taxis = ortho_c2_axes[0], ortho_c2_axes[1], ortho_c2_axes[2]

                # Detect reflection symmetry (any sigma plane)
                sigmav_chk, _ = _is_there_sigmav(positions, masses, geom_tol, list_sea_subset, paxis)
                sigmah_chk = _is_there_sigmah(positions, masses, geom_tol, paxis)

                # Detect improper rotation S4 (characteristic of Td/Th)
                S4 = Sn(paxis, 4)
                has_S4 = transform_isequivalent(positions, masses, geom_tol, S4)

                if invertable:
                    pg = "Th"
                # elif sigmav_chk or sigmah_chk or has_S4:
                elif has_S4: # Must have S4
                    pg = "Td"
                else:
                    pg = "T"
                found_is_spherical = True
                break
            except RuntimeError:
                pass
    if found_is_spherical == False:
        logger.warning("Molecule was wrongly classified as a spherical top, (number of c2 axes is %d), probably due to high eigen_tol or geom_tol. " \
            "Process will continue as general symmetry.", n)
        return _classify_general(mol, positions, masses, geom_tol)

    return PointGroupResult(pg=pg, paxis=paxis, saxis=saxis)

_validate_all_c2_ortho 

_validate_all_c2_ortho(positions, masses, geom_tol, paxis, c2_ortho, Cn_order)

Validate the c2 orthogonal rotations for a molecule. Generates all the c2 orthogonal axis from c2_ortho, generates their rotation matrix and validates it.

Source code in molsympy/core/pg_detect.py 
167
168
169
170
171
172
173
174
175
176
177
def _validate_all_c2_ortho(positions, masses, geom_tol, paxis, c2_ortho, Cn_order):
    """
    Validate the c2 orthogonal rotations for a molecule.
    Generates all the c2 orthogonal axis from c2_ortho,
    generates their rotation matrix and validates it.
    """
    c2_ortho_axes = generate_cyclic_axes(paxis, c2_ortho, Cn_order)
    for c2_ortho_axis in c2_ortho_axes:
        if not transform_isequivalent(positions, masses, geom_tol, Cn(c2_ortho_axis, 2)):
            return False
    return True

_validate_all_sigmav 

_validate_all_sigmav(positions, masses, geom_tol, paxis, sigmav, Cn_order)

Validate the vertical mirror planes for a molecule. Generates all the norm axis of the mirror planes from sigmav, generates their reflection matrix and validates it.

Source code in molsympy/core/pg_detect.py 
179
180
181
182
183
184
185
186
187
188
189
def _validate_all_sigmav(positions, masses, geom_tol, paxis, sigmav, Cn_order):
    """
    Validate the vertical mirror planes for a molecule.
    Generates all the norm axis of the mirror planes from sigmav,
    generates their reflection matrix and validates it.
    """
    sigmav_axes = generate_cyclic_axes(paxis, sigmav, Cn_order)
    for sigmav_axis in sigmav_axes:
        if not transform_isequivalent(positions, masses, geom_tol, reflection_matrix(normalize(sigmav_axis))):
            return False
    return True

_classify_subfamily 

_classify_subfamily(seas, positions, masses, geom_tol, paxis, Cn_order)

Determine the point-group subfamily (h/v/d/S2n/pure) once paxis and Cn_order are known. Returns the full Schoenflies symbol and updated saxis.

Source code in molsympy/core/pg_detect.py 
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
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
def _classify_subfamily(seas, positions, masses, geom_tol, paxis, Cn_order):
    """
    Determine the point-group subfamily (h/v/d/S2n/pure) once paxis and
    Cn_order are known. Returns the full Schoenflies symbol and updated saxis.
    """
    saxis = np.zeros(3)
    list_sea_subset = to_typed_list([sea.subset for sea in seas])
    ortho_c2_chk, c2_ortho = _is_there_ortho_c2(positions, masses, geom_tol, list_sea_subset, paxis)
    sigmav_chk, sigmav = _is_there_sigmav(positions, masses, geom_tol, list_sea_subset, paxis)
    sigmah_chk = _is_there_sigmah(positions, masses, geom_tol, paxis)
    inversion_chk = transform_isequivalent(positions, masses, geom_tol, inversion_matrix())
    if sigmah_chk:
        sn_chk = validate_sn_subrotations(Cn_order, paxis, positions, masses, geom_tol, sigmah_chk)

    if ortho_c2_chk:
        c2_ortho = normalize(c2_ortho - np.dot(c2_ortho, paxis) * paxis)
        ortho_c2_chk = _validate_all_c2_ortho(positions, masses, geom_tol, paxis, c2_ortho, Cn_order)

    if sigmav_chk:
        if ortho_c2_chk:
            if sigmah_chk:
                # Dnh
                sigmav = normalize(np.cross(c2_ortho, paxis))
                sigmav_chk = _validate_all_sigmav(positions, masses, geom_tol, paxis, sigmav, Cn_order)
            else:
                # Dnd
                sigmav = normalize(sum(generate_cyclic_axes(paxis, c2_ortho, Cn_order, 2)))
                sigmav_chk = _validate_all_sigmav(positions, masses, geom_tol, paxis, sigmav, Cn_order)
        else:
            # Cnv
            sigmav = normalize(sigmav - np.dot(sigmav, paxis) * paxis)
            sigmav_chk = _validate_all_sigmav(positions, masses, geom_tol, paxis, sigmav, Cn_order)

    if ortho_c2_chk:
        saxis = c2_ortho
        if sigmah_chk and (Cn_order % 2 or inversion_chk) and sigmav_chk and sn_chk:
            pg = "D" + str(Cn_order) + "h"
        elif sigmav_chk:
            s2n_chk = validate_sn_subrotations(Cn_order * 2, paxis, positions, masses, geom_tol)
            if s2n_chk:
                pg = "D" + str(Cn_order) + "d"
            else:
                pg = "D" + str(Cn_order)
        else:
            pg = "D" + str(Cn_order)
    elif sigmah_chk and (Cn_order % 2 or inversion_chk) and sn_chk:
        pg = "C" + str(Cn_order) + "h"
    elif sigmav_chk:
        pg = "C" + str(Cn_order) + "v"
        if mol_is_planar(positions, geom_tol):
            saxis = _planar_mol_axis(positions)
        elif not (sigmav == np.zeros(3)).all() and hasattr(sigmav, '__len__') and any(sigmav):
            saxis = normalize(np.cross(paxis, sigmav))
    else:
        s2n_chk = validate_sn_subrotations(Cn_order * 2, paxis, positions, masses, geom_tol)
        if s2n_chk:
            pg = "S" + str(2 * Cn_order)
        else:
            pg = "C" + str(Cn_order)

    return pg, saxis

_classify_general 

_classify_general(mol, positions, masses, geom_tol)

Classify a general symmetry. Returns PointGroupResult.

Source code in molsympy/core/pg_detect.py 
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
281
282
283
def _classify_general(mol, positions, masses, geom_tol):
    """
    Classify a general symmetry.
    Returns PointGroupResult.
    """
    paxis = np.zeros(3)
    seas = find_SEAs(mol)
    rot_set = _find_rotation_sets(mol, seas)
    rots = _find_rotations(mol, rot_set)

    if len(rots) >= 1:
        Cn_order = rots[0].order
        paxis = rots[0].axis
    else:
        list_sea_subset = to_typed_list([sea.subset for sea in seas])
        c2 = _find_a_c2(positions, masses, geom_tol, list_sea_subset)
        if (c2 == np.zeros(3)).all():
            # No proper rotation -> Ci, Cs, or C1.
            if transform_isequivalent(positions, masses, geom_tol, inversion_matrix()):
                return PointGroupResult(pg="Ci", paxis=paxis, saxis=np.zeros(3))
            sigmav_chk, sigmav = _is_there_sigmav(positions, masses, geom_tol, list_sea_subset, np.zeros(3))
            if sigmav_chk:
                if not (sigmav == np.zeros(3)).all():
                    paxis = sigmav
                return PointGroupResult(pg="Cs", paxis=paxis, saxis=np.zeros(3))
            return PointGroupResult(pg="C1", paxis=paxis, saxis=np.zeros(3))
        paxis = c2
        Cn_order = 2

    pg, saxis = _classify_subfamily(seas, positions, masses, geom_tol, paxis, Cn_order)
    return PointGroupResult(pg=pg, paxis=paxis, saxis=saxis)

find_point_group 

find_point_group(mol)

Find the point group of a molecule.

Returns the Schoenflies symbol and the primary/secondary axes that define the canonical orientation with respect to the generated symmetry elements.

Algorithm summary (Beruski & Vidal 2013)
  1. Compute the moment-of-inertia tensor (MOIT) and its eigenvalues Ia <= Ib <= Ic.
  2. Classify the rotor type from the eigenvalues:
  3. Linear (Ia ~= 0): C0v or D0h.
  4. Spherical top (Ia ~= Ib ~= Ic): T, Td, Th, O, Oh, I, Ih.
  5. Symmetric top and Asymmetric rotor: Cn, Cnv, Cnh, Dn, Dnh, Dnd, Sn, C1, Cs, Ci.
  6. Subfamily (h/v/d) is determined from sigma_h, sigma_v, and ortho-C2.

Based on: Beruski, Otávio; Vidal, Luciano N. J. Comp. Chem. 2013. doi:10.1002/jcc.23493

Parameters:

  • mol (Atoms) –

    Molecule with mol.info["geom_tol"] set to the geometric tolerance.

Returns:

Raises:

  • Exception –

    If mol.info["geom_tol"] is not set.

Source code in molsympy/core/pg_detect.py 
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
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
def find_point_group(mol):
    """
    Find the point group of a molecule.

    Returns the Schoenflies symbol and the primary/secondary axes that define
    the canonical orientation with respect to the generated symmetry elements.

    Algorithm summary (Beruski & Vidal 2013)
    ----------------------------------------
    1. Compute the moment-of-inertia tensor (MOIT) and its eigenvalues
       Ia <= Ib <= Ic.
    2. Classify the rotor type from the eigenvalues:
       - Linear (Ia ~= 0): C0v or D0h.
       - Spherical top (Ia ~= Ib ~= Ic): T, Td, Th, O, Oh, I, Ih.
       - Symmetric top and Asymmetric rotor: Cn, Cnv, Cnh, Dn, Dnh, Dnd, Sn, C1, Cs, Ci.
    3. Subfamily (h/v/d) is determined from sigma_h, sigma_v, and ortho-C2.

    Based on:
        Beruski, Otávio; Vidal, Luciano N. J. Comp. Chem. 2013.
        doi:10.1002/jcc.23493

    Parameters
    ----------
    mol : ase.Atoms
        Molecule with mol.info["geom_tol"] set to the geometric tolerance.

    Returns
    -------
    : PointGroupResult

    Raises
    ------
    Exception
        If mol.info["geom_tol"] is not set.
    """
    try:
        geom_tol = mol.info["geom_tol"]
    except KeyError:
        raise Exception("Atoms object geometric tolerance hasn't been set. Set it with Atoms.info[\"geom_tol\"]=.")

    try:
        eigen_tol = mol.info['eigen_tol']
    except KeyError:
        raise Exception("Atoms object eigen tolerance hasn't been set. Set it with Atoms.info[\"eigen_tol\"]=.")

    mol = mol.copy()
    positions = mol.positions
    masses = mol.get_masses()

    # Diagonalize MOIT; eigenvectors become candidate symmetry axes.
    moit = calcmoit(mol)
    evals_mol, evecs_mol = np.linalg.eigh(moit)
    _idx = evals_mol.argsort()
    Ia_mol, Ib_mol, Ic_mol = evals_mol[_idx]

    if inertia_isclose(Ia_mol, 0.0, atol=geom_tol, rtol=eigen_tol):
        return _classify_linear(positions, masses, geom_tol)

    elif inertia_isclose(Ia_mol, Ib_mol, atol=geom_tol, rtol=eigen_tol) and inertia_isclose(Ia_mol, Ic_mol, atol=geom_tol, rtol=eigen_tol):
        return _classify_spherical_top(mol, positions, masses, geom_tol)

    return _classify_general(mol, positions, masses, geom_tol)

rotation_detection

rotation_detection.py — Search for proper rotation axes (Cn).

Public names consumed by pg_detect.py: RotationElement, _find_rotation_sets, _find_rotations, _linear_mol_axis, _find_a_c2, _is_there_ortho_c2, _num_C2, _highest_order_axis

Referenced by:

RotationElement 

RotationElement(axis, order)

Data structure holding a candidate rotation axis and its order.

Returned by: Used by:
Source code in molsympy/core/rotation_detection.py 
23
24
25
def __init__(self, axis, order):
    self.axis = axis
    self.order = order

_rotation_set_union 

_rotation_set_union(rotation_set: RotationElement) -> list[RotationElement]

Return the union of all per-SEA rotation sets.

Source code in molsympy/core/rotation_detection.py 
37
38
39
40
41
42
43
44
def _rotation_set_union(rotation_set:RotationElement) -> list[RotationElement]:
    """Return the union of all per-SEA rotation sets."""
    out = rotation_set[0]
    if len(rotation_set) > 1:
        for i in range(1, len(rotation_set)):
            out.extend(rotation_set[i][:])
    out = sort_descending_RotationElements(unique_RotationElements(out))
    return out

unique_RotationElements 

unique_RotationElements(elements)

Return a list of RotationElement objects with duplicates removed, using eq.

Source code in molsympy/core/rotation_detection.py 
48
49
50
51
52
53
54
def unique_RotationElements(elements):
    """Return a list of RotationElement objects with duplicates removed, using __eq__."""
    unique = []
    for el in elements:
        if not any(el == u for u in unique):  # uses __eq__
            unique.append(el)
    return unique

sort_descending_RotationElements 

sort_descending_RotationElements(elements)

Return a list of RotationElement objects sorted by RotationElement.order in descending order.

Source code in molsympy/core/rotation_detection.py 
56
57
58
59
def sort_descending_RotationElements(elements):
    """Return a list of RotationElement objects sorted by
    RotationElement.order in descending order."""
    return sorted(elements, key=lambda x: x.order, reverse=True)

_find_rotation_sets 

_find_rotation_sets(mol, SEAs) -> list[list[RotationElement]]

For each SEA, find the set of possible RotationElements.

Parameters:

  • mol (Atoms) –
  • SEAs (List[SEA]) –

Returns:

Source code in molsympy/core/rotation_detection.py 
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 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
def _find_rotation_sets(mol, SEAs) -> list[list[RotationElement]]:
    """
    For each SEA, find the set of possible RotationElements.

    Parameters
    ----------
    mol : ase.Atoms
    SEAs : List[SEA]

    Returns
    -------
    List[List[RotationElement]]
    """
    geom_tol = mol.info["geom_tol"]
    eigen_tol = mol.info['eigen_tol']
    out_all_SEAs = []
    for sea in SEAs:
        length = len(sea.subset)
        out_per_SEA = []
        if length < 2:
            sea.label = "Single Atom"
        elif length == 2:
            sea.label = "Linear"
            sea.axis = normalize(mol[sea.subset[0]].position)
        else:
            sea_mol = mol[sea.subset]
            sea_mol.translate(-sea_mol.get_center_of_mass())
            evals, evecs = np.linalg.eigh(calcmoit(mol[sea.subset]))
            idx = evals.argsort()
            Ia, Ib, Ic = evals[idx]
            Iav, Ibv, Icv = [evecs[:, i] for i in idx]
            if inertia_isclose(Ia, Ib, atol=geom_tol, rtol=eigen_tol) and inertia_isclose(Ia, Ic, atol=geom_tol, rtol=eigen_tol):
                sea.label = "Spherical"
            elif inertia_isclose(Ia + Ib, Ic, atol=geom_tol, rtol=eigen_tol):
                axis = Icv
                sea.axis = axis
                if inertia_isclose(Ia, Ib, atol=geom_tol, rtol=eigen_tol):
                    sea.label = "Regular Polygon"
                    for i in range(2, length + 1):
                        if isfactor(length, i):
                            out_per_SEA.append(RotationElement(axis, i))
                else:
                    sea.label = "Irregular Polygon"
                    for i in range(2, length):
                        if isfactor(length, i):
                            out_per_SEA.append(RotationElement(axis, i))
            else:
                if not (inertia_isclose(Ia, Ib, atol=geom_tol, rtol=eigen_tol) or inertia_isclose(Ib, Ic, atol=geom_tol, rtol=eigen_tol)):
                    sea.label = "Asymmetric Rotor"
                    for ax in [Iav, Ibv, Icv]:
                        out_per_SEA.append(RotationElement(ax, 2))
                else:
                    if inertia_isclose(Ia, Ib, atol=geom_tol, rtol=eigen_tol):
                        sea.label = "Oblate Symmetric Top"
                        axis = Icv
                        sea.axis = Icv
                    else:
                        sea.label = "Prolate Symmetric Top"
                        axis = Iav
                        sea.axis = Iav
                    k = length // 2
                    for i in range(2, k + 1):
                        if isfactor(k, i):
                            out_per_SEA.append(RotationElement(axis, i))
            if len(out_per_SEA) > 0:
                out_all_SEAs.append(out_per_SEA)
    return out_all_SEAs

_find_rotations 

_find_rotations(mol, rotation_set)

Find RotationElements in rotation_set that leave the molecule invariant.

Parameters:

Returns:

Source code in molsympy/core/rotation_detection.py 
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
161
162
163
164
165
166
def _find_rotations(mol, rotation_set):
    """
    Find RotationElements in rotation_set that leave the molecule invariant.

    Parameters
    ----------
    mol : ase.Atoms
    rotation_set : List[List[RotationElement]]

    Returns
    -------
    List[RotationElement]
    """
    positions = mol.positions
    masses = mol.get_masses()
    geom_tol = mol.info["geom_tol"]
    eigen_tol = mol.info['eigen_tol']
    if len(rotation_set) < 1:
        return []
    molmoit = calcmoit(mol)
    evals = np.sort(np.linalg.eigh(molmoit)[0])
    if evals[0] == 0.0 and inertia_isclose(evals[1], evals[2], atol=geom_tol, rtol=eigen_tol):
        for i in range(np.shape(positions)[0]):
            norm = normalize(positions[i, :])
            if not (norm == np.zeros(3)).all():
                axis = norm
                break
        re = RotationElement(axis, 0)
        return [re]
    rsu = _rotation_set_union(rotation_set)
    for i in rsu:
        rotation_is_valid = validate_cn_subrotations(i.order, i.axis, positions, masses, geom_tol)
        if rotation_is_valid:
            return [i]
    return []

validate_cn_subrotations 

validate_cn_subrotations(order: int, axis: ndarray, positions: ndarray, masses: ndarray, geom_tol: float) -> bool

Validate subrotations

Parameters:

  • order (int) –

    Order of rotation

  • axis (ndarray) –

    Axis of rotation

  • positions (ndarray) –

    Atom positions, shape(n, 3)

  • masses (ndarray) –

    Atom masses, shape(n, )

  • geom_tol (float) –

    Geometric tolerance

Returns:

  • bool –

    True if all subrotations and rotation are valid

Source code in molsympy/core/rotation_detection.py 
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
@njit(cache=True)
def validate_cn_subrotations(order:int, axis:np.ndarray, positions:np.ndarray, masses:np.ndarray, geom_tol:float) -> bool:
    """
    Validate subrotations

    Parameters
    ----------
    order: int
        Order of rotation
    axis: np.ndarray
        Axis of rotation
    positions: np.ndarray
        Atom positions, shape(n, 3)
    masses: np.ndarray
        Atom masses, shape(n, )
    geom_tol: float
        Geometric tolerance

    Returns
    -------
    : bool
        True if all subrotations and rotation are valid
    """
    for subrotation in range(1, order):
        rmat = Cn(axis, (order/subrotation))
        if not transform_isequivalent(positions, masses, geom_tol, rmat):
            return False

    return True

validate_sn_subrotations 

validate_sn_subrotations(order: int, axis: ndarray, positions: ndarray, masses: ndarray, geom_tol: float, has_sigmah: bool = False) -> bool

Validate improper subrotations, Sn

Parameters:

  • order (int) –

    Order of rotation

  • axis (ndarray) –

    Axis of rotation

  • positions (ndarray) –

    Atom positions, shape(n, 3)

  • masses (ndarray) –

    Atom masses, shape(n, )

  • geom_tol (float) –

    Geometric tolerance

  • has_sigmah (bool, default: False ) –

    If True improper rotation is done for each rotation in order. Else it is only done for reachable states of pure improper rotations.

Returns:

  • bool –

    True if all subrotations and rotation are valid

Source code in molsympy/core/rotation_detection.py 
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
@njit(cache=True)
def validate_sn_subrotations(order:int, axis:np.ndarray, positions:np.ndarray, masses:np.ndarray, geom_tol:float, has_sigmah:bool=False) -> bool:
    """
    Validate improper subrotations, Sn

    Parameters
    ----------
    order: int
        Order of rotation
    axis: np.ndarray
        Axis of rotation
    positions: np.ndarray
        Atom positions, shape(n, 3)
    masses: np.ndarray
        Atom masses, shape(n, )
    geom_tol: float
        Geometric tolerance
    has_sigmah: bool
        If True improper rotation is done for each rotation in order.
        Else it is only done for reachable states of pure improper rotations.

    Returns
    -------
    : bool
        True if all subrotations and rotation are valid
    """
    if has_sigmah:
        for subrotation in range(1, order):
            rmat = Sn(axis, (order/subrotation))
            if not transform_isequivalent(positions, masses, geom_tol, rmat):
                return False
    else:
        for subrotation in range(1, order * 2, 2):
            subrotation = subrotation % order
            if subrotation == 0:
                continue
            rmat = Sn(axis, (order/subrotation))
            if not transform_isequivalent(positions, masses, geom_tol, rmat):
                return False

    return True

_linear_mol_axis 

_linear_mol_axis(positions)

Return the axis that best aligns with a linear molecule.

Parameters:

  • positions (ndarray) –

Returns:

  • ndarray –

    shape (3,)

Source code in molsympy/core/rotation_detection.py 
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
@njit(cache=True)
def _linear_mol_axis(positions):
    """
    Return the axis that best aligns with a linear molecule.

    Parameters
    ----------
    positions : np.ndarray

    Returns
    -------
    : np.ndarray
        shape (3,)
    """
    coords = positions - mean_axis_0_numba(positions)
    _, _, vh = np.linalg.svd(coords, full_matrices=False)
    return normalize(vh[0])

_highest_order_axis 

_highest_order_axis(rotations)

Return the order of the highest-order rotation in the list.

Source code in molsympy/core/rotation_detection.py 
263
264
265
266
267
268
def _highest_order_axis(rotations):
    """Return the order of the highest-order rotation in the list."""
    m = rotations[0].order
    for elem_rot in rotations:
        m = max(m, elem_rot.order)
    return m

_compute_R_max 

_compute_R_max(positions, axis)

Return the maximum perpendicular distance of any atom from the given axis.

Source code in molsympy/core/rotation_detection.py 
273
274
275
276
277
278
279
@njit(cache=True)
def _compute_R_max(positions, axis):
    """Return the maximum perpendicular distance of any atom from the given axis."""
    axis = normalize(axis)
    proj = np.dot(positions, axis)[:, np.newaxis] * axis[np.newaxis, :]
    dists = vec_norm_axis(positions - proj)
    return np.max(dists)

_validate_c2_candidate 

_validate_c2_candidate(raw_axis, positions, masses, geom_tol, exclude_axis)

Shared C2 validation kernel: normalize → exclude-filter → invariance test.

Normalises raw_axis, rejects it if it coincides with exclude_axis (pass np.zeros(3) to disable the filter), and tests whether the corresponding C2 rotation leaves the molecule invariant.

Parameters:

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

    candidate vector (need not be unit)

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

    axis to reject (zeros(3) = no filter)

Returns:

  • (ndarray, shape(3)) –

    Normalised unit axis if the candidate passes all checks; zeros(3) otherwise.

Source code in molsympy/core/rotation_detection.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
306
307
308
309
310
@njit(cache=True)
def _validate_c2_candidate(raw_axis, positions, masses, geom_tol, exclude_axis):
    """
    Shared C2 validation kernel: normalize → exclude-filter → invariance test.

    Normalises *raw_axis*, rejects it if it coincides with *exclude_axis*
    (pass ``np.zeros(3)`` to disable the filter), and tests whether the
    corresponding C2 rotation leaves the molecule invariant.

    Parameters
    ----------
    raw_axis : np.ndarray, shape (3,)
        candidate vector (need not be unit)
    exclude_axis : np.ndarray, shape (3,)
        axis to reject (zeros(3) = no filter)

    Returns
    -------
    : np.ndarray, shape (3,)
        Normalised unit axis if the candidate passes all checks; zeros(3) otherwise.
    """
    c2_axis = normalize(raw_axis)
    if (c2_axis == np.zeros(3)).all():
        return np.zeros(3)
    if not (exclude_axis == np.zeros(3)).all() and issame_axis(c2_axis, exclude_axis):
        return np.zeros(3)
    if transform_isequivalent(positions, masses, geom_tol, Cn(c2_axis, 2)):
        return c2_axis
    return np.zeros(3)

_c2a 

_c2a(positions, masses, geom_tol, sea_subset, exclude_axis=zeros(3), max_num_c2=1, ortho_tol=-1)

Find C_2 axes from origin-to-midpoint vectors of atom pairs within a SEA.

Candidate = normalised (pos[i] + pos[j]).

Source code in molsympy/core/rotation_detection.py 
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
@njit(cache=True)
def _c2a(positions, masses, geom_tol, sea_subset,
         exclude_axis=np.zeros(3), max_num_c2=1, ortho_tol = -1):
    """
    Find C_2 axes from origin-to-midpoint vectors of atom pairs within a SEA.

    Candidate = normalised (pos[i] + pos[j]).
    """
    length = len(sea_subset)
    out = List.empty_list(types.float64[:])
    for ii in range(length):
        for jj in range(ii + 1, length):
            origin_midpoint = normalize(positions[sea_subset[ii], :] + positions[sea_subset[jj], :])
            if ortho_tol >= 0 and abs(np.dot(origin_midpoint, exclude_axis)) > ortho_tol:
                continue
            result = _validate_c2_candidate(origin_midpoint, positions, masses, geom_tol, exclude_axis)
            if not (result == np.zeros(3)).all():
                out.append(result)
                if len(out) >= max_num_c2:
                    out = get_unique_axes(out)
                    if len(out) >= max_num_c2:
                        return out
    return out

_c2b 

_c2b(positions, masses, geom_tol, sea_subset, exclude_axis=zeros(3), max_num_c2=1, ortho_tol=-1)

Find C_2 axes from individual atom position vectors within a SEA.

Candidate = normalised pos[i].

Source code in molsympy/core/rotation_detection.py 
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
@njit(cache=True)
def _c2b(positions, masses, geom_tol, sea_subset,
         exclude_axis=np.zeros(3), max_num_c2=1, ortho_tol = -1):
    """
    Find C_2 axes from individual atom position vectors within a SEA.

    Candidate = normalised pos[i].
    """
    length = len(sea_subset)
    out = List.empty_list(types.float64[:])
    for ii in range(length):
        from_origin = normalize(positions[sea_subset[ii]])
        if ortho_tol >= 0 and abs(np.dot(from_origin, exclude_axis)) > ortho_tol:
            continue
        result = _validate_c2_candidate(
            from_origin, positions, masses, geom_tol, exclude_axis
        )
        if not (result == np.zeros(3)).all():
            out.append(result)
            if len(out) >= max_num_c2:
                out = get_unique_axes(out)
                if len(out) >= max_num_c2:
                    return out
    return out

_c2c 

_c2c(positions, masses, geom_tol, sea1_subset, sea2_subset, exclude_axis=zeros(3), ortho_tol=-1)

Find a C_2 axis from the cross-product of two linear-SEA bond vectors.

Candidate = normalised cross(r_SEA1, r_SEA2). Returns the unit axis or zeros(3) if the candidate is rejected.

Source code in molsympy/core/rotation_detection.py 
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
@njit(cache=True)
def _c2c(positions, masses, geom_tol, sea1_subset, sea2_subset,
         exclude_axis=np.zeros(3), ortho_tol = -1):
    """
    Find a C_2 axis from the cross-product of two linear-SEA bond vectors.

    Candidate = normalised cross(r_SEA1, r_SEA2).
    Returns the unit axis or zeros(3) if the candidate is rejected.
    """
    rij = positions[sea1_subset[0], :] - positions[sea1_subset[1], :]
    rkl = positions[sea2_subset[0], :] - positions[sea2_subset[1], :]
    cross_prod = normalize(np.cross(rij, rkl))
    if ortho_tol >= 0 and abs(np.dot(cross_prod, exclude_axis)) > ortho_tol:
        return np.zeros(3)
    return _validate_c2_candidate(
        cross_prod, positions, masses, geom_tol, exclude_axis
    )

_find_a_c2 

_find_a_c2(positions: ndarray, masses: ndarray, geom_tol: float, list_sea_subset: list[ndarray]) -> ndarray

Search for any C_2 axis; return the first one found.

Parameters:

  • positions (ndarray) –
  • masses (ndarray) –
  • geom_tol (float) –
  • list_sea_subset (list[ndarray]) –

Returns:

  • (ndarray, shape(3)) –
Source code in molsympy/core/rotation_detection.py 
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
@njit(cache=True)
def _find_a_c2(positions: np.ndarray, masses: np.ndarray, geom_tol: float, list_sea_subset: list[np.ndarray]) -> np.ndarray:
    """
    Search for any C_2 axis; return the first one found.

    Parameters
    ----------
    positions: np.ndarray
    masses: np.ndarray
    geom_tol: float
    list_sea_subset: list[np.ndarray]

    Returns
    -------
    np.ndarray, shape (3,)
    """
    num_seas = len(list_sea_subset)
    for id_sea1 in range(num_seas):
        sea1 = list_sea_subset[id_sea1]
        a = _c2a(positions, masses, geom_tol, sea1)
        if len(a) != 0:
            return a[0]
        b = _c2b(positions, masses, geom_tol, sea1)
        if len(b) != 0:
            return b[0]
        if len(sea1) == 2:
            for id_sea2 in range(id_sea1, num_seas):
                sea2 = list_sea_subset[id_sea2]
                if len(sea2) == 2:
                    c = _c2c(positions, masses, geom_tol, sea1, sea2)
                    if not (c == np.zeros(3)).all():
                        return c
    return np.zeros(3)

_is_there_ortho_c2 

_is_there_ortho_c2(positions, masses, geom_tol, list_sea_subset, paxis)

Search for a C_2 axis orthogonal to paxis; return the first one found.

Parameters:

  • positions (ndarray) –
  • masses (ndarray) –
  • geom_tol (float) –
  • list_sea_subset –
  • paxis ((ndarray, shape(3))) –

Returns:

  • tuple(bool, ndarray) –
Source code in molsympy/core/rotation_detection.py 
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
@njit(cache=True)
def _is_there_ortho_c2(positions, masses, geom_tol, list_sea_subset, paxis):
    """
    Search for a C_2 axis orthogonal to paxis; return the first one found.

    Parameters
    ----------
    positions : np.ndarray
    masses : np.ndarray
    geom_tol : float
    list_sea_subset: list[np.ndarray]
    paxis : np.ndarray, shape (3,)

    Returns
    -------
    tuple(bool, np.ndarray)
    """
    ortho_tol = geom_tol / _compute_R_max(positions, paxis) * 1.10
    paxis = np.ascontiguousarray(paxis)
    positions = np.ascontiguousarray(positions)
    num_seas = len(list_sea_subset)
    for id_sea1 in range(num_seas):
        sea1 = list_sea_subset[id_sea1]
        b = _c2b(positions, masses, geom_tol, sea1, exclude_axis=paxis, ortho_tol=ortho_tol)
        if len(b) != 0:
            b = b[0]
            return True, b
        else:
            a = _c2a(positions, masses, geom_tol, sea1, exclude_axis=paxis, ortho_tol=ortho_tol)
            if len(a) != 0:
                a = a[0]
                return True, a
            else:
                if len(sea1) == 2:
                    for id_sea2 in range(id_sea1, num_seas):
                        sea2 = list_sea_subset[id_sea2]
                        if len(sea2) == 2:
                            c = _c2c(positions, masses, geom_tol, sea1, sea2,
                                     exclude_axis=paxis, ortho_tol=ortho_tol)
                            if not (c == np.zeros(3)).all():
                                return True, c
    return False, np.zeros(3)

_num_C2 

_num_C2(positions: ndarray, masses: ndarray, geom_tol: float, list_sea_subset: list[ndarray], max_num_c2: int) -> list[ndarray]

Count distinct C_2 axes and return them.

Parameters:

  • positions (ndarray) –
  • masses (ndarray) –
  • geom_tol (float) –
  • list_sea_subset (list[ndarray]) –
  • max_num_c2 (int) –

Returns:

  • list[ndarray] –
Source code in molsympy/core/rotation_detection.py 
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
@njit(cache=True)
def _num_C2(positions: np.ndarray, masses: np.ndarray, geom_tol: float, list_sea_subset: list[np.ndarray], max_num_c2: int) -> list[np.ndarray]:
    """
    Count distinct C_2 axes and return them.

    Parameters
    ----------
    positions: np.ndarray
    masses: np.ndarray
    geom_tol: float
    list_sea_subset: list[np.ndarray]
    max_num_c2: int

    Returns
    -------
    list[np.ndarray]
    """
    axes = List.empty_list(types.float64[:])
    for subset in list_sea_subset:
        a = _c2a(positions, masses, geom_tol, subset, max_num_c2=max_num_c2)
        if len(a) != 0:
            axes.extend(a)
        b = _c2b(positions, masses, geom_tol, subset, max_num_c2=max_num_c2)
        if len(b) != 0:
            axes.extend(b)
        if len(axes) > max_num_c2:
            axes = get_unique_axes(axes)
            if len(axes) >= max_num_c2:
                break
    if len(axes) < 1:
        return axes
    unique_axes = get_unique_axes(axes)
    return unique_axes

reflection_detection

reflection_detection.py — Search for reflection planes (sigma).

Public names consumed by pg_detect.py: _is_there_sigmah, _is_there_sigmav, mol_is_planar, _planar_mol_axis

Referenced by:

_is_there_sigmah 

_is_there_sigmah(positions, masses, geom_tol, paxis)

Check for a horizontal reflection plane (normal = paxis).

Parameters:

  • positions (ndarray) –
  • masses (ndarray) –
  • geom_tol (float) –
  • paxis ((ndarray, shape(3))) –

Returns:

  • bool –
Source code in molsympy/core/reflection_detection.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
@njit(cache=True)
def _is_there_sigmah(positions, masses, geom_tol, paxis):
    """
    Check for a horizontal reflection plane (normal = paxis).

    Parameters
    ----------
    positions : np.ndarray
    masses : np.ndarray
    geom_tol : float
    paxis : np.ndarray, shape (3,)

    Returns
    -------
    bool
    """
    sigmah = reflection_matrix(paxis)
    return transform_isequivalent(positions, masses, geom_tol, sigmah)

_is_there_sigmav 

_is_there_sigmav(positions, masses, geom_tol, list_sea_subset, paxis)

Check for vertical reflection planes (normal perpendicular to paxis).

Parameters:

  • positions (ndarray) –
  • masses (ndarray) –
  • geom_tol (float) –
  • list_sea_subset –
  • paxis ((ndarray, shape(3))) –

Returns:

  • tuple(bool, ndarray) –
Source code in molsympy/core/reflection_detection.py 
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
74
75
76
77
78
@njit(cache=True)
def _is_there_sigmav(positions, masses, geom_tol, list_sea_subset, paxis):
    """
    Check for vertical reflection planes (normal perpendicular to paxis).

    Parameters
    ----------
    positions : np.ndarray
    masses : np.ndarray
    geom_tol : float
    list_sea_subset: list[np.ndarray]
    paxis : np.ndarray, shape (3,)

    Returns
    -------
    tuple(bool, np.ndarray)
    """
    axes = List.empty_list(types.float64[:])
    for sea in list_sea_subset:
        length = len(sea)
        if length < 2:
            continue
        A = sea[0]
        for i in range(1, length):
            B = sea[i]
            n = normalize(positions[A, :] - positions[B, :])
            if not (n == np.zeros(3)).all():
                sigma = reflection_matrix(n)
                if transform_isequivalent(positions, masses, geom_tol, sigma):
                    axes.append(n)
    if len(axes) < 1:
        if mol_is_planar(positions, geom_tol):
            return True, _planar_mol_axis(positions)
        else:
            return False, np.zeros(3)
    unique_axes = get_unique_axes(axes)
    for i in unique_axes:
        if not issame_axis(i, paxis):
            return True, i
    return False, np.zeros(3)

mol_is_planar 

mol_is_planar(positions, geom_tol)

Check if all atoms lie in a common plane.

Parameters:

  • positions (Atoms) –

Returns:

  • bool –
Source code in molsympy/core/reflection_detection.py 
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
@njit(cache=True)
def mol_is_planar(positions, geom_tol):
    """
    Check if all atoms lie in a common plane.

    Parameters
    ----------
    positions : ase.Atoms

    Returns
    -------
    bool
    """
    rank = np.linalg.matrix_rank(positions, tol=geom_tol)
    if rank < 3:
        axis = _planar_mol_axis(positions)
        new_positions, _, _ = _jit_rotate_mol_to_symels(positions, axis, np.array([0.0, 0.0, 0.0]))
        matrix = np.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, -1.0]])
        positions_B = transform(new_positions, matrix)
        for i in range(len(positions)):
            if np.linalg.norm(new_positions[i,:] - positions_B[i,:]) >= geom_tol:
                return False
        return True
    return False

_planar_mol_axis 

_planar_mol_axis(positions)

Return the normal to the plane of a planar molecule.

Parameters:

  • positions (ndarray) –

Returns:

  • ndarray –

    shape (3,)

Source code in molsympy/core/reflection_detection.py 
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
@njit(cache=True)
def _planar_mol_axis(positions):
    """
    Return the normal to the plane of a planar molecule.

    Parameters
    ----------
    positions : np.ndarray

    Returns
    -------
    : np.ndarray
        shape (3,)
    """
    coords = positions - mean_axis_0_numba(positions)
    _, _, vh = np.linalg.svd(coords, full_matrices=False)
    return normalize(vh[-1])

special_geometry

special_geometry.py — Special-geometry axis detection for icosahedral (Ih/I) and octahedral (Oh/O) point groups.

Public names consumed by pg_detect.py: _find_C3s_for_Ih, _find_C4s_for_Oh

Referenced by:

_jit_find_C3s_for_Ih 

_jit_find_C3s_for_Ih(size, positions, masses, geom_tol)

Find the 10 unique C3 axes of an icosahedral (Ih/I) molecule.

Geometric basis

The icosahedron has 20 triangular faces; each face defines one C3 axis (through its centroid), but opposite faces share an axis, giving 10 distinct axes. Each C3 axis passes through a pair of antipodal triangular face-centers.

Strategy

Enumerate all triples (i, j, k) of atoms. A triple forms an equilateral triangle when all three pairwise squared distances are equal within geom_tol. The C3 axis candidate is the normal to the plane of the triangle. The candidate is accepted only if the full C3 rotation leaves the molecule invariant.

Complexity: O(n^3) atom triples.

Parameters:

  • size (int) –
  • positions ((ndarray, shape(n, 3))) –
  • masses ((ndarray, shape(n))) –
  • geom_tol (float) –

Returns:

  • (list[ndarray], shape(3)) –

    exactly 10 unique unit C3-axis vectors.

Raises:

  • Exception –

    If the number of unique axes found is not 10.

Source code in molsympy/core/special_geometry.py 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
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
74
75
76
77
78
79
80
81
@njit(cache=True)
def _jit_find_C3s_for_Ih(size, positions, masses, geom_tol):
    """
    Find the 10 unique C3 axes of an icosahedral (Ih/I) molecule.

    Geometric basis
    ---------------
    The icosahedron has 20 triangular faces; each face defines one C3 axis
    (through its centroid), but opposite faces share an axis, giving 10
    distinct axes. Each C3 axis passes through a pair of antipodal triangular
    face-centers.

    Strategy
    --------
    Enumerate all triples (i, j, k) of atoms. A triple forms an equilateral
    triangle when all three pairwise squared distances are equal within geom_tol.
    The C3 axis candidate is the normal to the plane of the triangle.
    The candidate is accepted only if the full C3 rotation leaves the molecule
    invariant.

    Complexity: O(n^3) atom triples.

    Parameters
    ----------
    size : int
    positions : np.ndarray, shape (n, 3)
    masses : np.ndarray, shape (n,)
    geom_tol : float

    Returns
    -------
    list[np.ndarray], shape (3,)
        exactly 10 unique unit C3-axis vectors.

    Raises
    ------
    Exception
        If the number of unique axes found is not 10.
    """
    c3_axes = List.empty_list(types.float64[:])
    for i in range(size):
        for j in range(i + 1, size):
            for k in range(j + 1, size):
                rij = positions[i, :] - positions[j, :]
                rjk = positions[j, :] - positions[k, :]
                rik = positions[i, :] - positions[k, :]
                nij2 = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2]
                njk2 = rjk[0]*rjk[0] + rjk[1]*rjk[1] + rjk[2]*rjk[2]
                nik2 = rik[0]*rik[0] + rik[1]*rik[1] + rik[2]*rik[2]
                if float_isclose(nij2, njk2, atol=geom_tol) and float_isclose(nij2, nik2, atol=geom_tol):
                    c3_axis = normalize(np.cross(rij, rjk))
                    if not (c3_axis == np.zeros(3)).all():
                        c3 = Cn(c3_axis, 3)
                        if transform_isequivalent(positions, masses, geom_tol, c3):
                            c3_axes.append(c3_axis)
    if len(c3_axes) >= 1:
        return get_unique_axes(c3_axes)
    else:
        return List.empty_list(types.float64[:])

_find_C3s_for_Ih 

_find_C3s_for_Ih(mol, seas)

Find the 10 unique C3 axes for an Ih/I molecule so paxis and saxis can be defined.

Parameters:

  • mol (Atoms) –
  • seas –

Returns:

  • (List[ndarray], shape(3)) –
Source code in molsympy/core/special_geometry.py 
 84
 85
 86
 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
def _find_C3s_for_Ih(mol, seas):
    """
    Find the 10 unique C3 axes for an Ih/I molecule so paxis and saxis can be defined.

    Parameters
    ----------
    mol : ase.Atoms
    seas: list[SEA]

    Returns
    -------
    List[np.ndarray], shape (3,)
    """

    c3_axes = List.empty_list(types.float64[:])
    positions = mol.positions
    masses = mol.get_masses()
    geom_tol = mol.info["geom_tol"]
    for sea in seas:
        c3_axes.extend(_jit_find_C3s_for_Ih(len(sea.subset), positions[sea.subset], masses[sea.subset], geom_tol))
        c3_axes = get_unique_axes(c3_axes)
        chk = len(c3_axes)
        if chk == 10:
            return c3_axes
    logger.debug("C3 axes count: %d", chk)
    raise RuntimeError(
        "Unexpected number of C3 axes for Ih point group, expected 10."
    )

_check_square 

_check_square(va, vb, a, b, c, d, positions, masses, geom_tol)

Test whether four edge lengths form a square and, if so, return the C4 axis.

A square has four equal sides (a==b==c==d within geom_tol). The C4 axis is the normal to the plane of the square: cross(edge1, edge2).

Parameters:

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

    Two adjacent edge vectors of the candidate quadrilateral.

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

    Two adjacent edge vectors of the candidate quadrilateral.

  • a (float) –

    Squared lengths of the four sides to compare.

  • b (float) –

    Squared lengths of the four sides to compare.

  • c (float) –

    Squared lengths of the four sides to compare.

  • d (float) –

    Squared lengths of the four sides to compare.

  • positions (ndarray) –
  • masses (ndarray) –
  • geom_tol (float) –

Returns:

  • (ndarray, shape(3)) –

    Unit C4-axis if valid square and rotation leaves molecule invariant; otherwise a zero vector.

Source code in molsympy/core/special_geometry.py 
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
@njit(cache=True)
def _check_square(va, vb, a, b, c, d, positions, masses, geom_tol):
    """
    Test whether four edge lengths form a square and, if so, return the C4 axis.

    A square has four equal sides (a==b==c==d within geom_tol). The C4 axis
    is the normal to the plane of the square: cross(edge1, edge2).

    Parameters
    ----------
    va, vb : np.ndarray, shape (3,)
        Two adjacent edge vectors of the candidate quadrilateral.
    a, b, c, d : float
        Squared lengths of the four sides to compare.
    positions, masses : np.ndarray
    geom_tol : float

    Returns
    -------
    np.ndarray, shape (3,)
        Unit C4-axis if valid square and rotation leaves molecule invariant;
        otherwise a zero vector.
    """
    if (float_isclose(a, b, atol=geom_tol) and
            float_isclose(c, d, atol=geom_tol) and
            float_isclose(a, c, atol=geom_tol)):
        c4_axis = normalize(np.cross(va, vb))
        if not (c4_axis == np.zeros(3)).all():
            c4 = Cn(c4_axis, 4)
            if transform_isequivalent(positions, masses, geom_tol, c4):
                return c4_axis
    return np.zeros(3)

_jit_find_C4s_for_Oh 

_jit_find_C4s_for_Oh(size, positions, masses, geom_tol)

Find the 3 unique C4 axes of an octahedral (Oh/O) molecule.

Geometric basis

The octahedron has 6 vertices in 3 orthogonal pairs; each pair defines one C4 axis, giving 3 perpendicular axes.

Strategy

Enumerate all quadruples (i, j, k, l). Each is tested in three cyclic orderings to check whether any forms a square. Complexity: O(n^4).

Parameters:

  • size (int) –
  • positions ((ndarray, shape(n, 3))) –
  • masses ((ndarray, shape(n))) –
  • geom_tol (float) –

Returns:

  • list of np.ndarray, shape (3,) –

    exactly 3 unique unit C4-axis vectors.

Raises:

  • Exception –

    If the number of unique axes found is not 3.

Source code in molsympy/core/special_geometry.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
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
@njit(cache=True)
def _jit_find_C4s_for_Oh(size, positions, masses, geom_tol):
    """
    Find the 3 unique C4 axes of an octahedral (Oh/O) molecule.

    Geometric basis
    ---------------
    The octahedron has 6 vertices in 3 orthogonal pairs; each pair defines one
    C4 axis, giving 3 perpendicular axes.

    Strategy
    --------
    Enumerate all quadruples (i, j, k, l). Each is tested in three cyclic
    orderings to check whether any forms a square. Complexity: O(n^4).

    Parameters
    ----------
    size : int
    positions : np.ndarray, shape (n, 3)
    masses : np.ndarray, shape (n,)
    geom_tol : float

    Returns
    -------
    list of np.ndarray, shape (3,)
        exactly 3 unique unit C4-axis vectors.

    Raises
    ------
    Exception
        If the number of unique axes found is not 3.
    """
    c4_axes = List.empty_list(types.float64[:])
    for i in range(size):
        for j in range(i + 1, size):
            for k in range(j + 1, size):
                for l in range(k + 1, size):
                    if i != j and k != l and i != k:
                        rij = positions[i, :] - positions[j, :]
                        rjk = positions[j, :] - positions[k, :]
                        rkl = positions[k, :] - positions[l, :]
                        ril = positions[i, :] - positions[l, :]
                        rik = positions[i, :] - positions[k, :]
                        rjl = positions[j, :] - positions[l, :]

                        nij2 = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2]
                        njk2 = rjk[0]*rjk[0] + rjk[1]*rjk[1] + rjk[2]*rjk[2]
                        nkl2 = rkl[0]*rkl[0] + rkl[1]*rkl[1] + rkl[2]*rkl[2]
                        nil2 = ril[0]*ril[0] + ril[1]*ril[1] + ril[2]*ril[2]
                        nik2 = rik[0]*rik[0] + rik[1]*rik[1] + rik[2]*rik[2]
                        njl2 = rjl[0]*rjl[0] + rjl[1]*rjl[1] + rjl[2]*rjl[2]

                        c4_axis = _check_square(rij, rjk, nij2, njk2, nkl2, nil2,
                                                positions, masses, geom_tol)
                        if not (c4_axis == np.zeros(3)).all():
                            c4_axes.append(c4_axis)
                        c4_axis = _check_square(rij, rjl, nij2, njl2, nkl2, nik2,
                                                positions, masses, geom_tol)
                        if not (c4_axis == np.zeros(3)).all():
                            c4_axes.append(c4_axis)
                        c4_axis = _check_square(rik, rjk, nik2, njk2, njl2, nil2,
                                                positions, masses, geom_tol)
                        if not (c4_axis == np.zeros(3)).all():
                            c4_axes.append(c4_axis)
    if len(c4_axes) >= 1:
        return get_unique_axes(c4_axes)
    else:
        return List.empty_list(types.float64[:])

_find_C4s_for_Oh 

_find_C4s_for_Oh(mol, seas)

Find the 3 C4 axes for an Oh/O molecule so paxis and saxis can be defined.

Parameters:

  • mol (Atoms) –
  • seas (list[SEA]) –

Returns:

  • (List[ndarray], shape(3)) –
Source code in molsympy/core/special_geometry.py 
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
def _find_C4s_for_Oh(mol, seas):
    """
    Find the 3 C4 axes for an Oh/O molecule so paxis and saxis can be defined.

    Parameters
    ----------
    mol : ase.Atoms
    seas : list[SEA]

    Returns
    -------
    List[np.ndarray], shape (3,)
    """
    c4_axes = List.empty_list(types.float64[:])
    positions = mol.positions
    masses = mol.get_masses()
    geom_tol = mol.info["geom_tol"]
    for sea in seas:
        c4_axes.extend(_jit_find_C4s_for_Oh(len(sea.subset), positions[sea.subset], masses[sea.subset], geom_tol))
        c4_axes = get_unique_axes(c4_axes)
        chk = len(c4_axes)
        if chk == 3:
            return c4_axes
    logger.debug("c4 axes count %d", chk)
    raise RuntimeError(
            "Unexpected number of C4 axes for Oh point group, expected 3."
        )