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(mol, positions, masses, geom_tol)

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

Source code in minimalsym/core/pg_detect.py 
61
62
63
64
65
def _classify_linear(mol, positions, masses, geom_tol):
    """Classify a linear molecule (Ia ~= 0). Returns PointGroupResult."""
    paxis = _linear_mol_axis(mol)
    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 minimalsym/core/pg_detect.py 
 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
130
131
132
133
134
135
136
137
138
139
140
141
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)
    num_C2 = _num_C2(mol, seas)
    if num_C2 is None:
        if PRINT_WARNINGS:
            warnings.warn("Molecule was wrongly classified as a spherical top, (num_C2 is None), probably due to high eigen_tol. " \
                "Process will continue as general symmetry.")
        return _classify_general(mol, positions, masses, geom_tol)
    n, axes = num_C2
    invertable = transform_isequivalent(positions, masses, geom_tol, inversion_matrix())
    is_spherical = True

    if n >= 15:
        # Icosahedral: paxis = C5 axis, saxis = C2 axis from golden-ratio geometry.
        try:
            c2_axis = axes[0]
            c3s = _find_C3s_for_Ih(mol)
            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
            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)
            pg = "Ih" if invertable else "I"
        except RuntimeError:
            is_spherical = False
    elif n == 9:
        # Octahedral: paxis and saxis are two orthogonal C4 axes.
        try:
            c4s = _find_C4s_for_Oh(mol)
            paxis, saxis = c4s[0], c4s[1]
            pg = "Oh" if invertable else "O"
        except RuntimeError:
            is_spherical = False

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

            # Detect reflection symmetry (any sigma plane)
            sigmav_chk, _ = _is_there_sigmav(mol, seas, paxis)
            sigmah_chk = _is_there_sigmah(mol, 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"
        except RuntimeError:
            is_spherical = False
    else:
        is_spherical = False
    if not is_spherical:
        if PRINT_WARNINGS:
            warnings.warn(f"Molecule was wrongly classified as a spherical top, (n is {n}), probably due to high eigen_tol or geom_tol. " \
                "Process will continue as general symmetry.")
        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 minimalsym/core/pg_detect.py 
143
144
145
146
147
148
149
150
151
152
153
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 minimalsym/core/pg_detect.py 
155
156
157
158
159
160
161
162
163
164
165
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(mol, 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 minimalsym/core/pg_detect.py 
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
def _classify_subfamily(mol, 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)
    ortho_c2_chk, c2_ortho = _is_there_ortho_c2(mol, seas, paxis)
    sigmav_chk, sigmav = _is_there_sigmav(mol, seas, paxis)
    sigmah_chk = _is_there_sigmah(mol, paxis)

    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:
        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:
            pg = "D" + str(Cn_order) + "h"
        elif sigmav_chk:
            S2n = Sn(paxis, Cn_order * 2)
            if transform_isequivalent(positions, masses, geom_tol, S2n):
                pg = "D" + str(Cn_order) + "d"
            else:
                pg = "D" + str(Cn_order)
        else:
            pg = "D" + str(Cn_order)
    elif sigmah_chk:
        pg = "C" + str(Cn_order) + "h"
    elif sigmav_chk:
        pg = "C" + str(Cn_order) + "v"
        if mol_is_planar(mol):
            saxis = _planar_mol_axis(mol)
        elif sigmav is not None and hasattr(sigmav, '__len__') and any(sigmav):
            saxis = normalize(np.cross(paxis, sigmav))
    else:
        S2n = Sn(paxis, Cn_order * 2)
        if transform_isequivalent(positions, masses, geom_tol, S2n):
            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 minimalsym/core/pg_detect.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
239
240
241
242
243
244
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 = _highest_order_axis(rots)
        paxis = rots[0].axis
    else:
        c2 = _find_a_c2(mol, seas)
        if c2 is None:
            # 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(mol, seas, np.zeros(3))
            if sigmav_chk:
                if sigmav is not None:
                    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(mol, 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 minimalsym/core/pg_detect.py 
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
281
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
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(mol, 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 minimalsym/core/rotation_detection.py 
22
23
24
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 minimalsym/core/rotation_detection.py 
36
37
38
39
40
41
42
43
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 = unique_rotation_elements(out)
    return out

unique_rotation_elements 

unique_rotation_elements(elements)

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

Source code in minimalsym/core/rotation_detection.py 
47
48
49
50
51
52
53
def unique_rotation_elements(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

_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 minimalsym/core/rotation_detection.py 
 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
 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
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 minimalsym/core/rotation_detection.py 
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
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]):
            if normalize(positions[i, :]) is not None and not np.allclose(positions[i, :], np.zeros(3)):
                axis = normalize(positions[i, :])
                break
        re = RotationElement(axis, 0)
        return [re]
    rsu = _rotation_set_union(rotation_set)
    out = []
    for i in rsu:
        rotation_is_valid = True
        for suborder in range(2, i.order+1):
            if i.order % suborder == 0:
                rmat = Cn(i.axis, suborder)
                if not transform_isequivalent(positions, masses, geom_tol, rmat):
                    rotation_is_valid = False
                    break
        if rotation_is_valid:
            out.append(i)
    return out

_linear_mol_axis 

_linear_mol_axis(mol)

Return the axis that best aligns with a linear molecule.

Parameters:

  • mol (Atoms) –

Returns:

  • array –

    shape (3,)

Source code in minimalsym/core/rotation_detection.py 
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
def _linear_mol_axis(mol):
    """
    Return the axis that best aligns with a linear molecule.

    Parameters
    ----------
    mol : ase.Atoms

    Returns
    -------
    : np.array
        shape (3,)
    """
    coords = mol.positions - mol.positions.mean(axis=0)
    _, _, 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 minimalsym/core/rotation_detection.py 
189
190
191
192
193
194
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 minimalsym/core/rotation_detection.py 
199
200
201
202
203
204
205
@njit
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 minimalsym/core/rotation_detection.py 
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
@njit
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 is None or (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), return_all=False)

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

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

Source code in minimalsym/core/rotation_detection.py 
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
@njit
def _c2a(positions, masses, geom_tol, sea_subset,
         exclude_axis=np.zeros(3), return_all=False):
    """
    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 = []
    for i in range(length):
        for j in range(i + 1, length):
            raw = positions[sea_subset[i], :] + positions[sea_subset[j], :]
            result = _validate_c2_candidate(raw, positions, masses, geom_tol, exclude_axis)
            if not (result == np.zeros(3)).all():
                if return_all:
                    out.append(result)
                else:
                    return [result]
    return out

_c2b 

_c2b(positions, masses, geom_tol, sea_subset, exclude_axis=zeros(3), return_all=False)

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

Candidate = normalised pos[i].

Source code in minimalsym/core/rotation_detection.py 
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
@njit
def _c2b(positions, masses, geom_tol, sea_subset,
         exclude_axis=np.zeros(3), return_all=False):
    """
    Find C_2 axes from individual atom position vectors within a SEA.

    Candidate = normalised pos[i].
    """
    length = len(sea_subset)
    out = []
    for i in range(length):
        result = _validate_c2_candidate(
            positions[sea_subset[i], :], positions, masses, geom_tol, exclude_axis
        )
        if not (result == np.zeros(3)).all():
            if return_all:
                out.append(result)
            else:
                return [result]
    return out

_c2c 

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

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 minimalsym/core/rotation_detection.py 
283
284
285
286
287
288
289
290
291
292
293
294
295
296
@njit
def _c2c(positions, masses, geom_tol, sea1_subset, sea2_subset,
         exclude_axis=np.zeros(3)):
    """
    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], :]
    return _validate_c2_candidate(
        np.cross(rij, rkl), positions, masses, geom_tol, exclude_axis
    )

_find_a_c2 

_find_a_c2(mol, SEAs)

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

Parameters:

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

Returns:

  • (array, shape(3) or None) –
Source code in minimalsym/core/rotation_detection.py 
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
def _find_a_c2(mol, SEAs):
    """
    Search for any C_2 axis; return the first one found.

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

    Returns
    -------
    np.array, shape (3,) or None
    """
    positions = mol.positions
    masses = mol.get_masses()
    geom_tol = mol.info["geom_tol"]
    for sea in SEAs:
        a = _c2a(positions, masses, geom_tol, sea.subset)
        if len(a) != 0:
            return a[0]
        b = _c2b(positions, masses, geom_tol, sea.subset)
        if len(b) != 0:
            return b[0]
        if sea.label == "Linear":
            for sea2 in SEAs:
                if sea == sea2:
                    continue
                elif sea2.label == "Linear":
                    c = _c2c(positions, masses, geom_tol, sea.subset, sea2.subset)
                    if not (c == np.zeros(3)).all():
                        return c
    return None

_is_there_ortho_c2 

_is_there_ortho_c2(mol, SEAs, paxis)

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

Parameters:

  • mol (Atoms) –
  • SEAs (List[SEA]) –
  • paxis ((array, shape(3))) –

Returns:

  • tuple(bool, array or None) –
Source code in minimalsym/core/rotation_detection.py 
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
def _is_there_ortho_c2(mol, SEAs, paxis):
    """
    Search for a C_2 axis orthogonal to paxis; return the first one found.

    Parameters
    ----------
    mol : ase.Atoms
    SEAs : List[SEA]
    paxis : np.array, shape (3,)

    Returns
    -------
    tuple(bool, np.array or None)
    """
    ortho_tol = mol.info["geom_tol"] / _compute_R_max(mol.positions, paxis) * 1.10
    positions = mol.positions
    masses = mol.get_masses()
    geom_tol = mol.info["geom_tol"]
    for sea in SEAs:
        b = _c2b(positions, masses, geom_tol, sea.subset, exclude_axis=paxis)
        if len(b) != 0:
            b = b[0]
            if abs(np.dot(b, paxis)) <= ortho_tol:
                return True, b
        else:
            a = _c2a(positions, masses, geom_tol, sea.subset, exclude_axis=paxis)
            if len(a) != 0:
                a = a[0]
                if abs(np.dot(a, paxis)) <= ortho_tol:
                    return True, a
            else:
                if sea.label == "Linear":
                    for sea2 in SEAs:
                        if sea == sea2:
                            continue
                        elif sea2.label == "Linear":
                            c = _c2c(positions, masses, geom_tol, sea.subset, sea2.subset,
                                     exclude_axis=paxis)
                            if not (c == np.zeros(3)).all() and abs(np.dot(c, paxis)) <= ortho_tol:
                                return True, c
    return False, None

_num_C2 

_num_C2(mol, SEAs)

Count distinct C_2 axes and return them.

Parameters:

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

Returns:

  • tuple(int, List[array]) or None –
Source code in minimalsym/core/rotation_detection.py 
377
378
379
380
381
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
def _num_C2(mol, SEAs):
    """
    Count distinct C_2 axes and return them.

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

    Returns
    -------
    tuple(int, List[np.array]) or None
    """
    axes = []
    positions = mol.positions
    masses = mol.get_masses()
    geom_tol = mol.info["geom_tol"]
    for sea in SEAs:
        a = _c2a(positions, masses, geom_tol, sea.subset, return_all=True)
        if len(a) != 0:
            axes.extend(a)
        b = _c2b(positions, masses, geom_tol, sea.subset, return_all=True)
        if len(b) != 0:
            axes.extend(b)
    if len(axes) < 1:
        return None
    unique_axes = [axes[0]]
    for i in axes:
        check = True
        for j in unique_axes:
            if issame_axis(i, j):
                check = False
                break
        if check:
            unique_axes.append(i)
    return len(unique_axes), 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(mol, paxis)

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

Parameters:

  • mol (Atoms) –
  • paxis ((array, shape(3))) –

Returns:

  • bool –
Source code in minimalsym/core/reflection_detection.py 
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
def _is_there_sigmah(mol, paxis):
    """
    Check for a horizontal reflection plane (normal = paxis).

    Parameters
    ----------
    mol : ase.Atoms
    paxis : np.array, shape (3,)

    Returns
    -------
    bool
    """
    sigmah = reflection_matrix(paxis)
    return transform_isequivalent(mol.positions, mol.get_masses(), mol.info["geom_tol"], sigmah)

_is_there_sigmav 

_is_there_sigmav(mol, SEAs, paxis)

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

Parameters:

  • mol (Atoms) –
  • SEAs (List[SEA]) –
  • paxis ((array, shape(3))) –

Returns:

  • tuple(bool, array or None) –
Source code in minimalsym/core/reflection_detection.py 
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
def _is_there_sigmav(mol, SEAs, paxis):
    """
    Check for vertical reflection planes (normal perpendicular to paxis).

    Parameters
    ----------
    mol : ase.Atoms
    SEAs : List[SEA]
    paxis : np.array, shape (3,)

    Returns
    -------
    tuple(bool, np.array or None)
    """
    axes = []
    positions = mol.positions
    masses = mol.get_masses()
    geom_tol = mol.info["geom_tol"]
    for sea in SEAs:
        length = len(sea.subset)
        if length < 2:
            continue
        A = sea.subset[0]
        for i in range(1, length):
            B = sea.subset[i]
            n = normalize(positions[A, :] - positions[B, :])
            if n is not None and 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(mol):
            return True, _planar_mol_axis(mol)
        else:
            return False, None
    unique_axes = [axes[0]]
    for i in axes:
        check = True
        for j in unique_axes:
            if issame_axis(i, j):
                check = False
                break
        if check:
            unique_axes.append(i)
    for i in unique_axes:
        if not issame_axis(i, paxis):
            return True, i
    return False, None

mol_is_planar 

mol_is_planar(mol)

Check if all atoms lie in a common plane.

Parameters:

  • mol (Atoms) –

Returns:

  • bool –
Source code in minimalsym/core/reflection_detection.py 
 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
def mol_is_planar(mol):
    """
    Check if all atoms lie in a common plane.

    Parameters
    ----------
    mol : ase.Atoms

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

_planar_mol_axis 

_planar_mol_axis(mol)

Return the normal to the plane of a planar molecule.

Parameters:

  • mol (Atoms) –

Returns:

  • array –

    shape (3,) or None

Source code in minimalsym/core/reflection_detection.py 
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
def _planar_mol_axis(mol):
    """
    Return the normal to the plane of a planar molecule.

    Parameters
    ----------
    mol : ase.Atoms

    Returns
    -------
    : np.array
        shape (3,) or None
    """
    coords = mol.positions - mol.positions.mean(axis=0)
    _, _, 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 minimalsym/core/special_geometry.py 
19
20
21
22
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
82
83
84
85
86
87
88
89
90
91
92
93
@njit
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 = []
    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:
        unique_axes = [c3_axes[0]]
        for i in c3_axes:
            check = True
            for j in unique_axes:
                if issame_axis(i, j):
                    check = False
                    break
            if check:
                unique_axes.append(i)
        chk = len(unique_axes)
    else:
        chk = 0
    if chk != 10:
        if PRINT_WARNINGS:
            print("DEBUG: C3 axes count:", chk)
        raise RuntimeError(
            "Unexpected number of C3 axes for Ih point group, expected 10."
        )
    return unique_axes

_find_C3s_for_Ih 

_find_C3s_for_Ih(mol)

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

Parameters:

  • mol (Atoms) –

Returns:

  • (List[array], shape(3)) –
Source code in minimalsym/core/special_geometry.py 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
def _find_C3s_for_Ih(mol):
    """
    Find the 10 unique C3 axes for an Ih/I molecule so paxis and saxis can be defined.

    Parameters
    ----------
    mol : ase.Atoms

    Returns
    -------
    List[np.array], shape (3,)        
    """
    return _jit_find_C3s_for_Ih(len(mol), mol.positions, mol.get_masses(), mol.info["geom_tol"])

_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 minimalsym/core/special_geometry.py 
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
@njit
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 minimalsym/core/special_geometry.py 
147
148
149
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
218
219
220
221
222
223
224
225
226
227
228
229
230
231
@njit
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 = []
    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:
        unique_axes = [c4_axes[0]]
        for i in c4_axes:
            check = True
            for j in unique_axes:
                if issame_axis(i, j):
                    check = False
                    break
            if check:
                unique_axes.append(i)
        chk = len(unique_axes)
    else:
        chk = 0
    if chk != 3:
        if PRINT_WARNINGS:
            print("DEBUG: c4 axes count", chk)
        raise RuntimeError(
            "Unexpected number of C4 axes for Oh point group, expected 3."
        )
    return unique_axes

_find_C4s_for_Oh 

_find_C4s_for_Oh(mol)

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

Parameters:

  • mol (Atoms) –

Returns:

  • (List[array], shape(3)) –
Source code in minimalsym/core/special_geometry.py 
234
235
236
237
238
239
240
241
242
243
244
245
246
def _find_C4s_for_Oh(mol):
    """
    Find the 3 C4 axes for an Oh/O molecule so paxis and saxis can be defined.

    Parameters
    ----------
    mol : ase.Atoms

    Returns
    -------
    List[np.array], shape (3,)
    """
    return _jit_find_C4s_for_Oh(len(mol), mol.positions, mol.get_masses(), mol.info["geom_tol"])