Skip to content

geometry #

Compute internal coordinates (e.g. bond lengths).

Modules:

  • smee

    Differentiably evaluate energies of molecules using SMIRNOFF force fields

Functions:

  • compute_bond_vectors

    Computes the vectors between each atom pair specified by the atom_indices as

  • compute_angles

    Computes the angles [rad] between each atom triplet specified by the

  • compute_dihedrals

    Computes the dihedral angles [rad] between each atom quartet specified by the

  • polar_to_cartesian_coords

    Converts a set of polar coordinates into cartesian coordinates.

  • compute_v_site_coords

    Computes the positions of a set of virtual sites relative to a specified

  • add_v_site_coords

    Appends the coordinates of any virtual sites to a conformer (or batch of

compute_bond_vectors #

compute_bond_vectors(
    conformer: Tensor, atom_indices: Tensor
) -> tuple[Tensor, Tensor]

Computes the vectors between each atom pair specified by the atom_indices as well as their norms.

Parameters:

  • conformer (Tensor) –

    The conformer [Å] to compute the bond vectors for with shape=(n_atoms, 3) or shape=(n_confs, n_atoms, 3).

  • atom_indices (Tensor) –

    The indices of the atoms involved in each bond with shape=(n_bonds, 2)

Returns:

  • tuple[Tensor, Tensor]

    The bond vectors and their norms [Å].

Source code in smee/geometry.py
def compute_bond_vectors(
    conformer: torch.Tensor, atom_indices: torch.Tensor
) -> tuple[torch.Tensor, torch.Tensor]:
    """Computes the vectors between each atom pair specified by the ``atom_indices`` as
    well as their norms.

    Args:
        conformer: The conformer [Å] to compute the bond vectors for with
            ``shape=(n_atoms, 3)`` or ``shape=(n_confs, n_atoms, 3)``.
        atom_indices: The indices of the atoms involved in each bond with
            ``shape=(n_bonds, 2)``

    Returns:
        The bond vectors and their norms [Å].
    """

    if len(atom_indices) == 0:
        return (
            smee.utils.tensor_like([], other=conformer),
            smee.utils.tensor_like([], other=conformer),
        )

    is_batched = conformer.ndim == 3

    if not is_batched:
        conformer = torch.unsqueeze(conformer, 0)

    directions = conformer[:, atom_indices[:, 1]] - conformer[:, atom_indices[:, 0]]
    distances = torch.norm(directions, dim=-1)

    if not is_batched:
        directions = torch.squeeze(directions, dim=0)
        distances = torch.squeeze(distances, dim=0)

    return directions, distances

compute_angles #

compute_angles(
    conformer: Tensor, atom_indices: Tensor
) -> Tensor

Computes the angles [rad] between each atom triplet specified by the atom_indices.

Parameters:

  • conformer (Tensor) –

    The conformer [Å] to compute the angles for with shape=(n_atoms, 3) or shape=(n_confs, n_atoms, 3).

  • atom_indices (Tensor) –

    The indices of the atoms involved in each angle with shape=(n_angles, 3).

Returns:

  • Tensor

    The valence angles [rad].

Source code in smee/geometry.py
def compute_angles(conformer: torch.Tensor, atom_indices: torch.Tensor) -> torch.Tensor:
    """Computes the angles [rad] between each atom triplet specified by the
    ``atom_indices``.

    Args:
        conformer: The conformer [Å] to compute the angles for with
            ``shape=(n_atoms, 3)`` or ``shape=(n_confs, n_atoms, 3)``.
        atom_indices: The indices of the atoms involved in each angle with
            ``shape=(n_angles, 3)``.

    Returns:
        The valence angles [rad].
    """

    if len(atom_indices) == 0:
        return smee.utils.tensor_like([], other=conformer)

    is_batched = conformer.ndim == 3

    if not is_batched:
        conformer = torch.unsqueeze(conformer, 0)

    vector_ab = conformer[:, atom_indices[:, 1]] - conformer[:, atom_indices[:, 0]]
    vector_ac = conformer[:, atom_indices[:, 1]] - conformer[:, atom_indices[:, 2]]

    # tan theta = sin theta / cos theta
    #
    # ||a x b|| = ||a|| ||b|| sin theta
    #   a . b   = ||a|| ||b|| cos theta
    #
    # => tan theta = (a x b) / (a . b)
    angles = torch.atan2(
        torch.norm(torch.cross(vector_ab, vector_ac, dim=-1), dim=-1),
        (vector_ab * vector_ac).sum(dim=-1),
    )

    if not is_batched:
        angles = torch.squeeze(angles, dim=0)

    return angles

compute_dihedrals #

compute_dihedrals(
    conformer: Tensor, atom_indices: Tensor
) -> Tensor

Computes the dihedral angles [rad] between each atom quartet specified by the atom_indices.

Parameters:

  • conformer (Tensor) –

    The conformer [Å] to compute the dihedral angles for with shape=(n_atoms, 3) or shape=(n_confs, n_atoms, 3).

  • atom_indices (Tensor) –

    The indices of the atoms involved in each dihedral angle with shape=(n_dihedrals, 4).

Returns:

  • Tensor

    The dihedral angles [rad].

Source code in smee/geometry.py
def compute_dihedrals(
    conformer: torch.Tensor, atom_indices: torch.Tensor
) -> torch.Tensor:
    """Computes the dihedral angles [rad] between each atom quartet specified by the
    ``atom_indices``.

    Args:
        conformer: The conformer [Å] to compute the dihedral angles for with
            ``shape=(n_atoms, 3)`` or ``shape=(n_confs, n_atoms, 3)``.
        atom_indices: The indices of the atoms involved in each dihedral angle with
            ``shape=(n_dihedrals, 4)``.

    Returns:
        The dihedral angles [rad].
    """

    if len(atom_indices) == 0:
        return smee.utils.tensor_like([], other=conformer)

    is_batched = conformer.ndim == 3

    if not is_batched:
        conformer = torch.unsqueeze(conformer, 0)

    # Based on the OpenMM formalism.
    vector_ab = conformer[:, atom_indices[:, 0]] - conformer[:, atom_indices[:, 1]]
    vector_cb = conformer[:, atom_indices[:, 2]] - conformer[:, atom_indices[:, 1]]
    vector_cd = conformer[:, atom_indices[:, 2]] - conformer[:, atom_indices[:, 3]]

    vector_ab_cross_cb = torch.cross(vector_ab, vector_cb, dim=-1)
    vector_cb_cross_cd = torch.cross(vector_cb, vector_cd, dim=-1)

    vector_cb_norm = torch.norm(vector_cb, dim=-1).unsqueeze(-1)

    y = (
        torch.cross(vector_ab_cross_cb, vector_cb_cross_cd, dim=-1)
        * vector_cb
        / vector_cb_norm
    ).sum(axis=-1)

    x = (vector_ab_cross_cb * vector_cb_cross_cd).sum(axis=-1)

    phi = torch.atan2(y, x)

    if not is_batched:
        phi = torch.squeeze(phi, dim=0)

    return phi

polar_to_cartesian_coords #

polar_to_cartesian_coords(polar_coords: Tensor) -> Tensor

Converts a set of polar coordinates into cartesian coordinates.

Parameters:

  • polar_coords (Tensor) –

    The polar coordinates with shape=(n_coords, 3) and with columns of distance [Å], 'in plane angle' [rad] and 'out of plane' angle [rad].

Returns:

  • Tensor

    An array of the cartesian coordinates with shape=(n_coords, 3) and units of [Å].

Source code in smee/geometry.py
def polar_to_cartesian_coords(polar_coords: torch.Tensor) -> torch.Tensor:
    """Converts a set of polar coordinates into cartesian coordinates.

    Args:
        polar_coords: The polar coordinates with ``shape=(n_coords, 3)`` and with
            columns of distance [Å], 'in plane angle' [rad] and 'out of plane'
            angle [rad].

    Returns:
        An array of the cartesian coordinates with ``shape=(n_coords, 3)`` and units
        of [Å].
    """

    d, theta, phi = polar_coords[:, 0], polar_coords[:, 1], polar_coords[:, 2]

    cos_theta = torch.cos(theta)
    sin_theta = torch.sin(theta)

    cos_phi = torch.cos(phi)
    sin_phi = torch.sin(phi)

    # Here we use cos(phi) in place of sin(phi) and sin(phi) in place of cos(phi)
    # this is because we want phi=0 to represent a 0 degree angle from the x-y plane
    # rather than 0 degrees from the z-axis.
    coords = torch.stack(
        [d * cos_theta * cos_phi, d * sin_theta * cos_phi, d * sin_phi], dim=-1
    )
    return coords

compute_v_site_coords #

compute_v_site_coords(
    v_sites: VSiteMap,
    conformer: Tensor,
    force_field: TensorForceField,
) -> Tensor

Computes the positions of a set of virtual sites relative to a specified conformer or batch of conformers.

Parameters:

  • v_sites (VSiteMap) –

    A mapping between the virtual sites to add and their corresponding force field parameters.

  • conformer (Tensor) –

    The conformer(s) to add the virtual sites to with shape=(n_atoms, 3) or shape=(n_batches, n_atoms, 3) and units of [Å].

  • force_field (TensorForceField) –

    The force field containing the virtual site parameters.

Returns:

  • Tensor

    A tensor of virtual site positions [Å] with shape=(n_v_sites, 3) or shape=(n_batches, n_v_sites, 3).

Source code in smee/geometry.py
def compute_v_site_coords(
    v_sites: "smee.VSiteMap",
    conformer: torch.Tensor,
    force_field: "smee.TensorForceField",
) -> torch.Tensor:
    """Computes the positions of a set of virtual sites relative to a specified
    conformer or batch of conformers.

    Args:
        v_sites: A mapping between the virtual sites to add and their corresponding
            force field parameters.
        conformer: The conformer(s) to add the virtual sites to with
            ``shape=(n_atoms, 3)`` or ``shape=(n_batches, n_atoms, 3)`` and units of
            [Å].
        force_field: The force field containing the virtual site parameters.

    Returns:
        A tensor of virtual site positions [Å] with ``shape=(n_v_sites, 3)`` or
        ``shape=(n_batches, n_v_sites, 3)``.
    """

    is_batched = conformer.ndim == 3

    if not is_batched:
        conformer = torch.unsqueeze(conformer, 0)

    if len(v_sites.parameter_idxs) > 0:
        local_frame_coords = force_field.v_sites.parameters[v_sites.parameter_idxs]
        local_coord_frames = _build_v_site_coord_frames(v_sites, conformer, force_field)

        v_site_coords = _convert_v_site_coords(local_frame_coords, local_coord_frames)
    else:
        v_site_coords = smee.utils.zeros_like((len(conformer), 0, 3), other=conformer)

    if not is_batched:
        v_site_coords = torch.squeeze(v_site_coords, 0)

    return v_site_coords

add_v_site_coords #

add_v_site_coords(
    v_sites: VSiteMap,
    conformer: Tensor,
    force_field: TensorForceField,
) -> Tensor

Appends the coordinates of any virtual sites to a conformer (or batch of conformers) containing only atomic coordinates.

Notes
  • This function only supports appending v-sites to the end of the list of coordinates, and not interleaving them between existing atomic coordinates.

Parameters:

  • v_sites (VSiteMap) –

    A mapping between the virtual sites to add and their corresponding force field parameters.

  • conformer (Tensor) –

    The conformer(s) to add the virtual sites to with shape=(n_atoms, 3) or shape=(n_batches, n_atoms, 3) and units of [Å].

  • force_field (TensorForceField) –

    The force field containing the virtual site parameters.

Returns:

  • Tensor

    The full conformer(s) with both atomic and virtual site coordinates [Å] with shape=(n_atoms+n_v_sites, 3) or shape=(n_batches, n_atoms+n_v_sites, 3).

Source code in smee/geometry.py
def add_v_site_coords(
    v_sites: "smee.VSiteMap",
    conformer: torch.Tensor,
    force_field: "smee.TensorForceField",
) -> torch.Tensor:
    """Appends the coordinates of any virtual sites to a conformer (or batch of
    conformers) containing only atomic coordinates.

    Notes:
        * This function only supports appending v-sites to the end of the list of
          coordinates, and not interleaving them between existing atomic coordinates.

    Args:
        v_sites: A mapping between the virtual sites to add and their corresponding
            force field parameters.
        conformer: The conformer(s) to add the virtual sites to with
            ``shape=(n_atoms, 3)`` or ``shape=(n_batches, n_atoms, 3)`` and units of
            [Å].
        force_field: The force field containing the virtual site parameters.

    Returns:
        The full conformer(s) with both atomic and virtual site coordinates [Å] with
        ``shape=(n_atoms+n_v_sites, 3)`` or ``shape=(n_batches, n_atoms+n_v_sites, 3)``.
    """

    v_site_coords = compute_v_site_coords(v_sites, conformer, force_field)

    return torch.cat([conformer, v_site_coords], dim=(1 if conformer.ndim == 3 else 0))