Skip to content

Geometry & Molecular Operations

mol_ops

Referenced by:

SEA dataclass 

SEA(label: str, subset: ndarray, axis: ndarray)

Symmetry equivalent atoms (SEA).

SEAs are atoms that can be swapped with no distinguishable change in the molecule.

Parameters:

  • label (str or None) –

    Rotor type of the SEA set (e.g. Single Atom, Linear, Spherical, Regular Polygon, Oblate Symmetric Top).

  • subset (ndarray) –

    Atom indices in the molecule that belong to this SEA set, shape (N,).

  • axis (ndarray or None) –

    Candidate rotational symmetry axis, shape (3,).

Returned by:

transform 

transform(positions: ndarray, M: ndarray) -> ndarray

Transform coordinates of molecule by matrix M and return new positions.

Parameters:

  • positions (ndarray) –

    Molecule positions.

  • M (ndarray) –

    Transformation matrix (e.g. rotation, reflection, etc.), shape (3,3)

Returns:

  • Atoms –

    Molecule with transformed atom coordinates

Source code in molsympy/core/mol_ops.py 
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
@njit(cache=True)
def transform(positions: "np.ndarray", M: "np.ndarray") -> "np.ndarray":
    """
    Transform coordinates of molecule by matrix M and return new positions.

    Parameters
    ----------
    positions: np.ndarray
        Molecule positions.
    M: np.ndarray
        Transformation matrix (e.g. rotation, reflection, etc.), shape (3,3)

    Returns
    -------
    ase.Atoms
        Molecule with transformed atom coordinates
    """
    return np.dot(positions, np.transpose(M))

distance_matrix 

distance_matrix(positions)

Calculate the interatomic distance matrix as all pairwise distances between atoms.

Parameters:

  • positions –

    Numpy array of positions of Molecule.

Returns:

  • ndarray –

    Interatomic distance matrix, shape(len(mol),len(mol))

Source code in molsympy/core/mol_ops.py 
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
@njit(cache=True)
def distance_matrix(positions):
    """
    Calculate the interatomic distance matrix as all pairwise distances between atoms.

    Parameters
    ----------
    positions: np.ndarray
        Numpy array of positions of Molecule.

    Returns
    -------
    np.ndarray
        Interatomic distance matrix, shape(len(mol),len(mol))
    """
    natoms = len(positions)
    dm = np.zeros((natoms, natoms))
    for i in range(natoms):
        for j in range(i, natoms):
            dx = positions[i,0] - positions[j,0]
            dy = positions[i,1] - positions[j,1]
            dz = positions[i,2] - positions[j,2]
            dm[i,j] = np.sqrt(dx*dx + dy*dy + dz*dz)
            dm[j,i] = dm[i,j]
    return dm

find_SEAs 

find_SEAs(mol)

Find sets of symmetry equivalent atoms. Permutations of the distance matrix reveal which atoms form symmetry equivalent sets.

Parameters:

  • mol –

    Molecule object.

Returns:

  • List[SEA] –

    List of symmetry equivalent atom sets

Source code in molsympy/core/mol_ops.py 
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
def find_SEAs(mol):
    """
    Find sets of symmetry equivalent atoms.
    Permutations of the distance matrix reveal which atoms form symmetry equivalent sets.

    Parameters
    ----------
    mol: ase.Atoms
        Molecule object.

    Returns
    -------
    List[SEA]
        List of symmetry equivalent atom sets
    """
    sea_id = _jit_find_SEAs(len(mol), mol.positions, mol.info["geom_tol"])
    n_seas = int(sea_id.max()) + 1
    SEAs = []
    for k in range(n_seas):
        subset = np.where(sea_id == k)[0].astype(np.int64)
        SEAs.append(SEA("", subset, np.zeros(3)))
    return SEAs

get_SEAs_from_atom_map 

get_SEAs_from_atom_map(atom_map)

Find sets of symmetry equivalent atoms based on a Symtext.atom_map. Form symmetry equivalent sets from who each atom maps to.

Parameters:

  • atom_map –

    Symtext.atom_map

Returns:

  • List[SEA] –

    List of symmetry equivalent atom sets

Source code in molsympy/core/mol_ops.py 
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
def get_SEAs_from_atom_map(atom_map):
    """
    Find sets of symmetry equivalent atoms based on a Symtext.atom_map.
    Form symmetry equivalent sets from who each atom maps to.

    Parameters
    ----------
    atom_map: np.ndarray
        Symtext.atom_map

    Returns
    -------
    List[SEA]
        List of symmetry equivalent atom sets
    """
    SEAs = []
    subset_list = _jit_get_SEAs_from_atom_map(atom_map)

    for subset in subset_list:
        SEAs.append(SEA("", subset, np.zeros(3)))

    return SEAs

calcmoit 

calcmoit(atoms)

Calculates the moment of inertia tensor for a list of atoms.

Parameters:

  • atoms –

    Set of atoms.

Returns:

  • ndarray –

    Cartesian moment of inertia tensor, shape(3,3).

Source code in molsympy/core/mol_ops.py 
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
def calcmoit(atoms):
    """
    Calculates the moment of inertia tensor for a list of atoms.

    Parameters
    ----------
    atoms: ase.Atoms
        Set of atoms.

    Returns
    -------
    np.ndarray
        Cartesian moment of inertia tensor, shape(3,3).
    """
    atoms.translate(-atoms.get_center_of_mass())
    masses = atoms.get_masses()
    positions = atoms.positions

    return _jit_calcmoit(positions, masses)

mol_orient

mol_orient.py — Molecular orientation utilities.

Provides rotate_mol_to_symels(), which aligns paxis with z and saxis with x so that the molecule sits in the canonical frame expected by pg_to_symels().

Referenced by:

_jit_rotate_mol_to_symels 

_jit_rotate_mol_to_symels(positions, paxis, saxis)

Align paxis to z-axis, saxis to x-axis.

Parameters:

  • positions (ndarray) –
  • paxis (ndarray) –
  • saxis (ndarray) –

Returns:

  • tuple(ndarray, ndarray, ndarray) –

    new_positions, rmat, rmat_inv

Source code in molsympy/core/mol_orient.py 
16
17
18
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
@njit(cache=True)
def _jit_rotate_mol_to_symels(positions, paxis, saxis):
    """
    Align paxis to z-axis, saxis to x-axis.

    Parameters
    ----------
    positions : np.ndarray
    paxis : np.ndarray
    saxis : np.ndarray

    Returns
    -------
    tuple(np.ndarray, np.ndarray, np.ndarray)
        new_positions, rmat, rmat_inv
    """
    if float_isclose(vec_norm(paxis), 0.0, atol=global_tol):
        rmat = rmat_inv = np.eye(3)
        return positions, rmat, rmat_inv
    z = paxis
    if float_isclose(vec_norm(saxis), 0.0, atol=global_tol):
        x = np.zeros(3)
        for trial_vec in np.eye(3):
            x = np.cross(trial_vec, z)
            if not float_isclose(vec_norm(x), 0.0, atol=global_tol):
                x = normalize(x)
                break
        y = normalize(np.cross(z, x))
    else:
        x = saxis
        y = np.cross(z, x)
    rmat = np.column_stack((x, y, z))
    rmat_inv = rmat.T
    new_positions = transform(positions, rmat_inv)
    return new_positions, rmat, rmat_inv

rotate_mol_to_symels 

rotate_mol_to_symels(mol, paxis, saxis)

Rotate molecule so that paxis aligns with z and saxis aligns with x.

Returns the rotated molecule and the forward/inverse rotation matrices so that computed properties can be rotated back to the original orientation.

Parameters:

  • mol (Atoms) –
  • paxis ((ndarray, shape(3))) –
  • saxis ((ndarray, shape(3))) –

Returns:

  • tuple(Atoms, ndarray, ndarray) –

    Rotated molecule, rotation matrix, inverse rotation matrix (shape 3×3).

Source code in molsympy/core/mol_orient.py 
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
def rotate_mol_to_symels(mol, paxis, saxis):
    """
    Rotate molecule so that paxis aligns with z and saxis aligns with x.

    Returns the rotated molecule and the forward/inverse rotation matrices so
    that computed properties can be rotated back to the original orientation.

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

    Returns
    -------
    tuple(ase.Atoms, np.ndarray, np.ndarray)
        Rotated molecule, rotation matrix, inverse rotation matrix (shape 3×3).
    """
    new_mol = mol.copy()
    paxis = np.asarray(paxis, dtype=np.float64)
    saxis = np.asarray(saxis, dtype=np.float64)
    positions, rmat, rmat_inv = _jit_rotate_mol_to_symels(mol.positions, paxis, saxis)
    new_mol.positions = positions
    return new_mol, rmat, rmat_inv

atom_mapping

atom_mapping.py — Atom permutation map builders.

Provides _get_atom_mapping() and _get_linear_atom_mapping(), which build the (n_atoms × n_symels) integer array that records where each atom goes under each symmetry operation.

Referenced by:

_where_you_go 

_where_you_go(positions, geom_tol, atom, rrep) -> int

Return the index of the atom that atom maps to under rrep. Return -1 for failure.

Source code in molsympy/core/atom_mapping.py 
16
17
18
19
20
21
22
23
24
25
26
27
28
@njit(cache=True)
def _where_you_go(positions, geom_tol, atom, rrep) -> int:
    """
    Return the index of the atom that *atom* maps to under *rrep*.
    Return -1 for failure.
    """
    ratom = np.dot(rrep, positions[atom, :].T)
    tol2 = geom_tol*geom_tol
    for i in range(len(positions)):
        dist = positions[i,:] - ratom
        if (dist[0]*dist[0] + dist[1]*dist[1] + dist[2]*dist[2]) < tol2:
            return i
    return np.int64(-1)

_jit_get_atom_mapping 

_jit_get_atom_mapping(positions, geom_tol, rreps)

Build the full (n_atoms × n_symels) permutation map.

Source code in molsympy/core/atom_mapping.py 
30
31
32
33
34
35
36
37
38
39
@njit(cache=True)
def _jit_get_atom_mapping(positions, geom_tol, rreps):
    """Build the full (n_atoms × n_symels) permutation map."""
    natoms = len(positions)
    nsymels = rreps.shape[0]
    amap = np.empty((natoms, nsymels), dtype=np.int64)
    for atom in range(natoms):
        for s in range(nsymels):
            amap[atom, s] = _where_you_go(positions, geom_tol, atom, rreps[s])
    return amap

_get_atom_mapping 

_get_atom_mapping(mol, symels)

Build the (n_atoms × n_symels) atom permutation map for non-linear groups.

Parameters:

  • mol (Atoms) –
  • symels (List[Symel]) –

Returns:

  • ndarray –

    shape (n_atoms, n_symels)

Raises:

  • Exception –

    If any atom fails to map under a symmetry operation.

Source code in molsympy/core/atom_mapping.py 
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
def _get_atom_mapping(mol, symels):
    """
    Build the (n_atoms × n_symels) atom permutation map for non-linear groups.

    Parameters
    ----------
    mol : ase.Atoms
    symels : List[Symel]

    Returns
    -------
    : np.ndarray
        shape (n_atoms, n_symels)

    Raises
    ------
    Exception
        If any atom fails to map under a symmetry operation.
    """
    rreps = np.array([s.rrep for s in symels])
    amap = _jit_get_atom_mapping(mol.positions, mol.info["geom_tol"]*1.1, rreps)
    failed = np.argwhere(amap == -1)
    if len(failed) > 0:
        atom, s = int(failed[0, 0]), int(failed[0, 1])
        raise Exception(
            f"Atom {atom} not mapped to another atom "
            f"under symel {symels[s]}\nPositions: {mol.positions}\nRreps: {rreps[s]}"
        )
    return amap

_get_linear_atom_mapping 

_get_linear_atom_mapping(mol, pg)

Build the permutation map for linear point groups (C0v / D0h).

For C0v: identity only (each atom maps to itself). For D0h: identity + inversion.

Parameters:

  • mol (Atoms) –
  • pg (PointGroup) –

Returns:

  • ndarray –
Source code in molsympy/core/atom_mapping.py 
 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
def _get_linear_atom_mapping(mol, pg):
    """
    Build the permutation map for linear point groups (C0v / D0h).

    For C0v: identity only (each atom maps to itself).
    For D0h: identity + inversion.

    Parameters
    ----------
    mol : ase.Atoms
    pg : PointGroup

    Returns
    -------
    np.ndarray
    """
    natoms = len(mol)
    positions = mol.positions
    geom_tol = mol.info["geom_tol"]
    amap = np.arange(natoms, dtype=int).reshape((natoms, 1))
    if pg.family == "D":
        ungerade_map = np.zeros(natoms, dtype=int)
        inversion_symel = Symel("i", None, -np.eye(3), None, None, None)
        for atom in range(natoms):
            w = _where_you_go(positions, geom_tol, atom, inversion_symel.rrep)
            if w != -1:
                ungerade_map[atom] = w
            else:
                raise Exception(f"Atom {atom} not mapped to another atom under symel i")
        return np.column_stack((amap, ungerade_map))
    return amap