ic

Compute internal coordinate representations of molecules.

Notes

• This module is heavily based off of the internal module of geomeTRIC. See the LICENSE-3RD-PARTY for license information.

Classes:

• ICType –

  The supported primitive internal coordinate types.

• IC –

  A base class for internal coordinate representations.

• RIC –

  A redundant internal coordinate representation.

• DLC –

  A delocalized internal coordinates representation

Attributes:

• ICDict –

  A dictionary of tensors partitioned by internal coordinate type.

• IC_ORDER –

  The order that internal coordinate types appear in a flattened representation.

• ConstraintDict –

  A dictionary of constraints for each type of internal coordinate. It is of the form

ICDict module-attribute

ICDict = dict[ICType, Tensor]

A dictionary of tensors partitioned by internal coordinate type.

IC_ORDER module-attribute

IC_ORDER = (DISTANCE, ANGLE, LINEAR, OUT_OF_PLANE, DIHEDRAL)

The order that internal coordinate types appear in a flattened representation.

ConstraintDict module-attribute

ConstraintDict = dict[ICType, tuple[Tensor, Tensor]]

A dictionary of constraints for each type of internal coordinate. It is of the form {ic_type: (idxs, values)} where idxs is a tensor of the indices of the atoms involved in the internal coordinates to constrain and values is a tensor of the values of that the internal coordinate should be constrained at.

ICType

Bases: Enum

The supported primitive internal coordinate types.

IC dataclass

Bases: ABC

A base class for internal coordinate representations.

Methods:

• compute_b –

  Computes jacobian of the internal coordinates with respect to the cartesian

• compute_q –

  Maps a set of cartesian coordinates to a set of internal coordinates.

• compute_dq –

  Compute the difference between two sets of coordinates in terms of internal

• guess_hess_q –

  Build an approximate Hessian that roughly follows Schlegel's guidelines.

• dq_to_x –

  Perturb a set of cartesian coordinates by a set of internal coordinate

Attributes:

• idxs (ICDict) –

  The indices of the atoms involved in each type of internal coordinate.

idxs instance-attribute

idxs: ICDict

The indices of the atoms involved in each type of internal coordinate.

compute_b abstractmethod

compute_b(coords_x: Tensor) -> Tensor

Computes jacobian of the internal coordinates with respect to the cartesian coordinates.

This is the B matrix in the Pulay-Fogarasi-Pulay approach.

Parameters:

• coords_x (Tensor) –

  The coordinates with shape=(n_atoms, 3).

Returns:

• Tensor –

  The jacobian, with shape=(n_ic, n_atoms * 3).

Source code in tico/ic.py

@abc.abstractmethod
def compute_b(self, coords_x: torch.Tensor) -> torch.Tensor:
    """Computes jacobian of the internal coordinates with respect to the cartesian
    coordinates.

    This is the B matrix in the Pulay-Fogarasi-Pulay approach.

    Args:
        coords_x: The coordinates with ``shape=(n_atoms, 3)``.

    Returns:
        The jacobian, with ``shape=(n_ic, n_atoms * 3)``.
    """

compute_q abstractmethod

compute_q(coords_x: Tensor) -> Tensor

Maps a set of cartesian coordinates to a set of internal coordinates.

Parameters:

• coords_x (Tensor) –

  The coordinates with shape=(n_atoms, 3).

Returns:

• Tensor –

  The flattened internal coordinate values with shape=(n_ic,).

Source code in tico/ic.py

@abc.abstractmethod
def compute_q(self, coords_x: torch.Tensor) -> torch.Tensor:
    """Maps a set of cartesian coordinates to a set of internal coordinates.

    Args:
        coords_x: The coordinates with ``shape=(n_atoms, 3)``.

    Returns:
        The flattened internal coordinate values with ``shape=(n_ic,)``.
    """

compute_dq abstractmethod

compute_dq(
    coords_x_a: Tensor, coords_x_b: Tensor
) -> Tensor

Compute the difference between two sets of coordinates in terms of internal coordinates.

Parameters:

• coords_x_a (Tensor) –

  The first coordinates with shape=(n_atoms, 3).

• coords_x_b (Tensor) –

  The second coordinates with shape=(n_atoms, 3).

Returns:

• Tensor –

  The internal coordinate displacements with shape=(n_ic,).

Source code in tico/ic.py

@abc.abstractmethod
def compute_dq(
    self, coords_x_a: torch.Tensor, coords_x_b: torch.Tensor
) -> torch.Tensor:
    """Compute the difference between two sets of coordinates in terms of internal
    coordinates.

    Args:
        coords_x_a: The first coordinates with ``shape=(n_atoms, 3)``.
        coords_x_b: The second coordinates with ``shape=(n_atoms, 3)``.

    Returns:
        The internal coordinate displacements with ``shape=(n_ic,)``.
    """

guess_hess_q abstractmethod

guess_hess_q(
    coords_x: Tensor, atomic_nums: Tensor
) -> Tensor

Build an approximate Hessian that roughly follows Schlegel's guidelines.

Parameters:

• coords_x (Tensor) –

  The coordinates with shape=(n_atoms, 3).

• atomic_nums (Tensor) –

  The atomic number of each atom with shape=(n_atoms,).

Returns:

• Tensor –

  The approximate Hessian with shape=(n_ic, n_ic).

Source code in tico/ic.py

@abc.abstractmethod
def guess_hess_q(
    self, coords_x: torch.Tensor, atomic_nums: torch.Tensor
) -> torch.Tensor:
    """Build an approximate Hessian that roughly follows Schlegel's guidelines.

    Args:
        coords_x: The coordinates with ``shape=(n_atoms, 3)``.
        atomic_nums: The atomic number of each atom with ``shape=(n_atoms,)``.

    Returns:
        The approximate Hessian with ``shape=(n_ic, n_ic)``.
    """

dq_to_x abstractmethod

dq_to_x(
    coords_x: Tensor, dq: Tensor
) -> tuple[Tensor, Tensor, bool]

Perturb a set of cartesian coordinates by a set of internal coordinate displacements.

Parameters:

• coords_x (Tensor) –

  The coordinates with shape=(n_atoms, 3).

• dq (Tensor) –

  The internal coordinate displacements with shape=(n_ic,).

Returns:

• tuple[Tensor, Tensor, bool] –

  The perturbed cartesian (with shape=(n_atoms, 3)) and internal coordinates (with shape=(n_ic,)), and a boolean indicating whether the perturbation was successful.

Source code in tico/ic.py

@abc.abstractmethod
def dq_to_x(
    self, coords_x: torch.Tensor, dq: torch.Tensor
) -> tuple[torch.Tensor, torch.Tensor, bool]:
    """Perturb a set of cartesian coordinates by a set of internal coordinate
    displacements.

    Args:
        coords_x: The coordinates with ``shape=(n_atoms, 3)``.
        dq: The internal coordinate displacements with ``shape=(n_ic,)``.

    Returns:
        The perturbed cartesian (with ``shape=(n_atoms, 3)``) and internal
        coordinates (with ``shape=(n_ic,)``), and a boolean indicating whether the
        perturbation was successful.
    """

RIC dataclass

Bases: IC

A redundant internal coordinate representation.

Methods:

• from_coords –

  Projects a set of cartesian coordinates onto a reduced set of delocalized

Attributes:

• idxs (ICDict) –

  The indices of the atoms involved in each type of internal coordinate.

idxs instance-attribute

idxs: ICDict

The indices of the atoms involved in each type of internal coordinate.

from_coords classmethod

from_coords(coords_x: Tensor, bond_idxs: Tensor) -> RIC

Projects a set of cartesian coordinates onto a reduced set of delocalized internal coordinates.

Parameters:

• coords_x (Tensor) –

  The coordinates with shape=(n_atoms, 3).

• bond_idxs (Tensor) –

  The atoms involved in each bond with shape=(n_bonds, 2).

Returns:

• RIC –

  The internal coordinate representation.

Source code in tico/ic.py

@classmethod
def from_coords(cls, coords_x: torch.Tensor, bond_idxs: torch.Tensor) -> "RIC":
    """Projects a set of cartesian coordinates onto a reduced set of delocalized
    internal coordinates.

    Args:
        coords_x: The coordinates with ``shape=(n_atoms, 3)``.
        bond_idxs: The atoms involved in each bond with ``shape=(n_bonds, 2)``.

    Returns:
        The internal coordinate representation.
    """
    coords_x = coords_x.double()

    graph = networkx.Graph(bond_idxs.detach().tolist())

    angle_idxs, linear_angle_idxs = _detect_angles(coords_x, graph)
    angle_idxs, out_of_plane_idxs = _detect_out_of_plane_angles(
        coords_x, graph, angle_idxs
    )

    dihedral_idxs = _detect_dihedrals(coords_x, graph, linear_angle_idxs)

    return_value = {
        ICType.DISTANCE: bond_idxs,
        ICType.ANGLE: angle_idxs,
        ICType.LINEAR: linear_angle_idxs,
        ICType.OUT_OF_PLANE: out_of_plane_idxs,
        ICType.DIHEDRAL: dihedral_idxs,
    }
    return RIC(
        {key: value for key, value in return_value.items() if len(value) > 0}
    )

DLC dataclass

Bases: IC

A delocalized internal coordinates representation

Methods:

• from_coords –

  Projects a set of cartesian coordinates onto a reduced set of delocalized

• compute_constr_delta –

  Compute the difference between the current internal coordinate values and the

• augment_hess_q –

  Augment a hessian with extra dimensions corresponding to the constrained

• project_grad_x –

  Project out the components of the gradient along the constrained degrees of

Attributes:

• idxs (ICDict) –

  The indices of the atoms involved in each type of internal coordinate.

• v (Tensor) –

  The projection of RIC onto the non-redundant basis.

• constr (ConstraintDict | None) –

  Constraints on the internal coordinates.

idxs instance-attribute

idxs: ICDict

The indices of the atoms involved in each type of internal coordinate.

v instance-attribute

v: Tensor

The projection of RIC onto the non-redundant basis.

constr class-attribute instance-attribute

constr: ConstraintDict | None = None

Constraints on the internal coordinates.

from_coords classmethod

from_coords(
    coords_x: Tensor,
    bond_idxs: Tensor,
    constr: ConstraintDict | None = None,
) -> DLC

Projects a set of cartesian coordinates onto a reduced set of delocalized internal coordinates.

Parameters:

• coords_x (Tensor) –

  The coordinates with shape=(n_atoms, 3).

• bond_idxs (Tensor) –

  The atoms involved in each bond with shape=(n_bonds, 2).

• constr (ConstraintDict | None, default: None ) –

  A dictionary of constraints on the internal coordinates. This dictionary should be of the form constr[ic_type] = (idxs, values) where idxs is a tensor of the indices of the atoms involved in the internal coordinate and values is a tensor of the values to constrain the internal coordinate at.

Returns:

• DLC –

  The internal coordinate representation.

Source code in tico/ic.py

@classmethod
def from_coords(
    cls,
    coords_x: torch.Tensor,
    bond_idxs: torch.Tensor,
    constr: ConstraintDict | None = None,
) -> "DLC":
    """Projects a set of cartesian coordinates onto a reduced set of delocalized
    internal coordinates.

    Args:
        coords_x: The coordinates with shape=(n_atoms, 3).
        bond_idxs: The atoms involved in each bond with shape=(n_bonds, 2).
        constr: A dictionary of constraints on the internal coordinates. This
            dictionary should be of the form ``constr[ic_type] = (idxs, values)``
            where ``idxs`` is a tensor of the indices of the atoms involved in the
            internal coordinate and ``values`` is a tensor of the values to
            constrain the internal coordinate at.

    Returns:
        The internal coordinate representation.
    """
    coords_x = coords_x.double()

    ric = RIC.from_coords(coords_x, bond_idxs)

    b_matrix = ric.compute_b(coords_x)
    g_matrix = b_matrix @ b_matrix.T

    eig_vals, eig_vecs = torch.linalg.eigh(g_matrix)
    eig_mask = torch.abs(eig_vals) > 1e-6

    v = eig_vecs[:, eig_mask]

    if constr is not None and len(constr) == 0:
        constr = None
    if constr is not None:
        v = cls._project_constr(v, constr, ric)

    return cls(ric.idxs, v, constr)

compute_constr_delta

compute_constr_delta(coords_x: Tensor) -> Tensor

Compute the difference between the current internal coordinate values and the constrained values.

Source code in tico/ic.py

def compute_constr_delta(self, coords_x: torch.Tensor) -> torch.Tensor:
    """Compute the difference between the current internal coordinate values and the
    constrained values."""
    constr_idxs = {ic_type: idxs for ic_type, (idxs, _) in self.constr.items()}
    constr_vals = {ic_type: vals for ic_type, (_, vals) in self.constr.items()}

    coords_q = _compute_q(coords_x, constr_idxs)

    return _compute_dq(constr_vals, coords_q)

augment_hess_q

augment_hess_q(
    coords_x: Tensor, grad_q: Tensor, hess_q: Tensor
)

Augment a hessian with extra dimensions corresponding to the constrained DoF.

Parameters:

• coords_x (Tensor) –

  The coordinates with shape=(n_atoms, 3).

• grad_q (Tensor) –

  The gradient of the internal coordinates with shape=(n_ic,).

• hess_q (Tensor) –

  The hessian of the internal coordinates with shape=(n_ic, n_ic).

Source code in tico/ic.py

def augment_hess_q(
    self, coords_x: torch.Tensor, grad_q: torch.Tensor, hess_q: torch.Tensor
):
    """Augment a hessian with extra dimensions corresponding to the constrained
    DoF.

    Args:
        coords_x: The coordinates with ``shape=(n_atoms, 3)``.
        grad_q: The gradient of the internal coordinates with ``shape=(n_ic,)``.
        hess_q: The hessian of the internal coordinates with ``shape=(n_ic, n_ic)``.
    """

    constr_to_ic_idx = _match_constr(self.constr, self.idxs)

    ni = len(grad_q)
    nc = len(constr_to_ic_idx)

    diag_idxs = list(range(len(constr_to_ic_idx)))

    ct = smee.utils.zeros_like((nc, ni), coords_x)
    ct[diag_idxs, diag_idxs] = 1.0 / self.v[constr_to_ic_idx, diag_idxs]

    nt = ni + nc

    constr_diff = -self.compute_constr_delta(coords_x)
    # au/rad (about 0.16 A / 17 deg)
    # constr_diff[:nc] = torch.clamp(constr_diff[:nc], -0.3, 0.3)

    hess_q_constr = smee.utils.zeros_like((nt, nt), coords_x)
    hess_q_constr[0:ni, 0:ni] = hess_q[:, :]
    hess_q_constr[ni:nt, 0:ni] = ct[:, :]
    hess_q_constr[0:ni, ni:nt] = ct.T[:, :]

    grad_q_constr = smee.utils.zeros_like(nt, coords_x)
    grad_q_constr[0:ni] = grad_q[:]
    grad_q_constr[ni:nt] = -constr_diff[:]

    return grad_q_constr, hess_q_constr

project_grad_x

project_grad_x(coords_x: Tensor, grad_x: Tensor) -> Tensor

Project out the components of the gradient along the constrained degrees of freedom.

Source code in tico/ic.py

def project_grad_x(
    self, coords_x: torch.Tensor, grad_x: torch.Tensor
) -> torch.Tensor:
    """Project out the components of the gradient along the constrained degrees of
    freedom."""

    if self.constr is None:
        return grad_x

    b_matrix = self.compute_b(coords_x)
    g_matrix_inv = tico.utils.pinv(b_matrix @ b_matrix.T)

    n_constr = sum(len(v) for i, v in self.constr.values())

    grad_q = torch.linalg.multi_dot([g_matrix_inv, b_matrix, grad_x])
    grad_q[:n_constr] = 0.0

    return torch.linalg.multi_dot([b_matrix.T, grad_q])