Repository: biotite-dev/biotite

Issue:
Support computation of lDDT
The *lDDT* (https://doi.org/10.1093/bioinformatics/btt473) has become a popular measure for evaluation of computational protein models against a reference model.

Using Biotite's `CellList` this should be rather easy and fast to implement.

Discussion:
@nscorley: Hi @padix-key -- I'm happy to pick this one up if you haven't progressed too far already, as we need this functionality for a project anyhow!@padix-key: Hi, I currently already have a WIP branch. Expect a PR soon.

Solution patch:
diff --git a/benchmarks/structure/benchmark_compare.py b/benchmarks/structure/benchmark_compare.py
new file mode 100644
index 00000000..0a04223f
--- /dev/null
+++ b/benchmarks/structure/benchmark_compare.py
@@ -0,0 +1,32 @@
+import itertools
+from pathlib import Path
+import pytest
+import biotite.structure as struc
+import biotite.structure.io.pdbx as pdbx
+from tests.util import data_dir
+
+
+@pytest.fixture
+def atoms():
+    pdbx_file = pdbx.BinaryCIFFile.read(Path(data_dir("structure")) / "1gya.bcif")
+    atoms = pdbx.get_structure(pdbx_file)
+    # Reduce the number of atoms to speed up the benchmark
+    return atoms[..., atoms.element != "H"]
+
+
+@pytest.mark.benchmark
+@pytest.mark.parametrize(
+    "multi_model, aggregation",
+    itertools.product([False, True], ["all", "chain", "residue", "atom"]),
+)
+def benchmark_lddt(atoms, multi_model, aggregation):
+    """
+    Compute lDDT on different aggregation levels.
+    """
+    reference = atoms[0]
+    if multi_model:
+        subject = atoms
+    else:
+        subject = atoms[0]
+
+    struc.lddt(reference, subject, aggregation=aggregation)
diff --git a/doc/apidoc.json b/doc/apidoc.json
index 281e9285..d53fb514 100644
--- a/doc/apidoc.json
+++ b/doc/apidoc.json
@@ -321,6 +321,7 @@
             "rmsd",
 	        "rmspd",
             "rmsf",
+            "lddt",
             "tm_score"
         ],
         "General analysis" : [
diff --git a/src/biotite/structure/celllist.pyx b/src/biotite/structure/celllist.pyx
index 712729d5..a61609b5 100644
--- a/src/biotite/structure/celllist.pyx
+++ b/src/biotite/structure/celllist.pyx
@@ -28,10 +28,10 @@ ctypedef np.uint8_t uint8
 cdef class CellList:
     """
     __init__(atom_array, cell_size, periodic=False, box=None, selection=None)
-    
+
     This class enables the efficient search of atoms in vicinity of a
     defined location.
-    
+
     This class stores the indices of an atom array in virtual "cells",
     each corresponding to a specific coordinate interval.
     If the atoms in vicinity of a specific location are searched, only
@@ -41,7 +41,7 @@ cdef class CellList:
     after the :class:`CellList` has been created.
     Therefore a :class:`CellList` saves calculation time in those
     cases, where vicinity is checked for multiple locations.
-    
+
     Parameters
     ----------
     atom_array : AtomArray or ndarray, dtype=float, shape=(n,3)
@@ -64,14 +64,14 @@ cdef class CellList:
         If provided, only the atoms masked by this array are stored in
         the cell list. However, the indices stored in the cell list
         will still refer to the original unfiltered `atom_array`.
-            
+
     Examples
     --------
-    
+
     >>> cell_list = CellList(atom_array, cell_size=5)
     >>> near_atoms = atom_array[cell_list.get_atoms(np.array([1,2,3]), radius=7.0)]
     """
-    
+
     # The atom coordinates
     cdef float32[:,:] _coord
     # A boolean mask that covers the selected atoms
@@ -99,8 +99,8 @@ cdef class CellList:
     cdef int _orig_length
     cdef float32[:] _orig_min_coord
     cdef float32[:] _orig_max_coord
-    
-    
+
+
     @cython.initializedcheck(False)
     @cython.boundscheck(False)
     @cython.wraparound(False)
@@ -119,14 +119,10 @@ cdef class CellList:
         # if 'periodic' is true
         self._orig_length = coord.shape[0]
         self._box = None
-        if coord.ndim != 2:
-            raise ValueError("Coordinates must have shape (n,3)")
-        if coord.shape[0] == 0:
-            raise ValueError("Coordinates must not be empty")
-        if coord.shape[1] != 3:
-            raise ValueError("Coordinates must have form (x,y,z)")
-        if np.isnan(coord).any():
-            raise ValueError("Coordinates contain NaN values")
+        if selection is None:
+            _check_coord(coord)
+        else:
+            _check_coord(coord[selection])
 
         if periodic:
             if box is not None:
@@ -145,7 +141,7 @@ cdef class CellList:
                 raise ValueError("Box contains NaN values")
             coord = move_inside_box(coord, self._box)
             coord, indices = repeat_box_coord(coord, self._box)
-        
+
         if self._has_initialized_cells():
             raise Exception("Duplicate call of constructor")
         self._cells = None
@@ -155,23 +151,23 @@ cdef class CellList:
         self._coord = coord.astype(np.float32, copy=False)
         self._cellsize = cell_size
         # calculate how many cells are required for each dimension
-        min_coord = np.min(coord, axis=0).astype(np.float32)
-        max_coord = np.max(coord, axis=0).astype(np.float32)
+        min_coord = np.nanmin(coord, axis=0).astype(np.float32)
+        max_coord = np.nanmax(coord, axis=0).astype(np.float32)
         self._min_coord = min_coord
         self._max_coord = max_coord
         cell_count = (((max_coord - min_coord) / cell_size) +1).astype(int)
         if self._periodic:
-            self._orig_min_coord = np.min(coord[:self._orig_length], axis=0) \
+            self._orig_min_coord = np.nanmin(coord[:self._orig_length], axis=0) \
                                    .astype(np.float32)
-            self._orig_max_coord = np.max(coord[:self._orig_length], axis=0) \
+            self._orig_max_coord = np.nanmax(coord[:self._orig_length], axis=0) \
                                    .astype(np.float32)
-        
+
         # ndarray of pointers to C-arrays
         # containing indices to atom array
         self._cells = np.zeros(cell_count, dtype=np.uint64)
         # Stores the length of the C-arrays
         self._cell_length = np.zeros(cell_count, dtype=np.int32)
-        
+
         # Prepare selection
         if selection is not None:
             self._has_selection = True
@@ -183,7 +179,7 @@ cdef class CellList:
                 )
         else:
             self._has_selection = False
-        
+
         # Fill cells
         for atom_array_i in range(self._coord.shape[0]):
             # Only put selected atoms into cell list
@@ -208,34 +204,34 @@ cdef class CellList:
                     # Store new cell pointer and length
                     self._cell_length[i,j,k] = length
                     self._cells[i,j,k] = <ptr> cell_ptr
-            
-    
+
+
     def __dealloc__(self):
         if self._has_initialized_cells():
             deallocate_ptrs(self._cells)
-    
-    
+
+
     @cython.initializedcheck(False)
     @cython.boundscheck(False)
     @cython.wraparound(False)
     def create_adjacency_matrix(self, float32 threshold_distance):
         """
         create_adjacency_matrix(threshold_distance)
-        
+
         Create an adjacency matrix for the atoms in this cell list.
 
         An adjacency matrix depicts which atoms *i* and *j* have
         a distance lower than a given threshold distance.
         The values in the adjacency matrix ``m`` are
         ``m[i,j] = 1 if distance(i,j) <= threshold else 0``
-        
+
         Parameters
         ----------
         threshold_distance : float
             The threshold distance. All atom pairs that have a distance
             lower than this value are indicated by ``True`` values in
             the resulting matrix.
-        
+
         Returns
         -------
         matrix : ndarray, dtype=bool, shape=(n,n)
@@ -244,7 +240,7 @@ cdef class CellList:
             :class:`CellList`, the rows and columns corresponding to
             atoms, that are not masked by the selection, have all
             elements set to ``False``.
-        
+
         Notes
         -----
         The highest performance is achieved when the the cell size is
@@ -269,7 +265,7 @@ cdef class CellList:
         if threshold_distance < 0:
             raise ValueError("Threshold must be a positive value")
         cdef int i=0
-        
+
         # Get atom position for all original positions
         # (no periodic copies)
         coord = np.asarray(self._coord[:self._orig_length])
@@ -287,17 +283,17 @@ cdef class CellList:
             return matrix
         else:
             return self.get_atoms(coord, threshold_distance, as_mask=True)
-        
-    
+
+
     @cython.initializedcheck(False)
     @cython.boundscheck(False)
     @cython.wraparound(False)
     def get_atoms(self, np.ndarray coord, radius, bint as_mask=False):
         """
         get_atoms(coord, radius, as_mask=False)
-        
+
         Find atoms with a maximum distance from given coordinates.
-        
+
         Parameters
         ----------
         coord : ndarray, dtype=float, shape=(3,) or shape=(m,3)
@@ -321,7 +317,7 @@ cdef class CellList:
         as_mask : bool, optional
             If true, the result is returned as boolean mask, instead
             of an index array.
-        
+
         Returns
         -------
         indices : ndarray, dtype=int32, shape=(p,) or shape=(m,p)
@@ -335,7 +331,7 @@ cdef class CellList:
             The values are true for atoms in the atom array,
             that are in the defined vicinity.
             Only returned with `as_mask` set to true.
-            
+
         See Also
         --------
         get_atoms_in_cells
@@ -376,7 +372,7 @@ cdef class CellList:
             A      18  PRO HA     H         2.719    3.181    1.316
             A      18  PRO HB3    H         2.781    3.223    3.618
             A      18  PRO CB     C         3.035    4.190    3.187
-        
+
         Get adjacent atoms for mutliple positions:
 
         >>> cell_list = CellList(atom_array, 3)
@@ -403,14 +399,14 @@ cdef class CellList:
         cdef float32 sq_radius
         cdef float32[:] sq_radii
         cdef np.ndarray cell_radii
-        
+
         cdef int[:,:] all_indices
         cdef int[:,:] indices
         cdef float32[:,:] coord_v
 
         if len(coord) == 0:
             return _empty_result(as_mask)
-        
+
         # Handle periodicity for the input coordinates
         if self._periodic:
             coord = move_inside_box(coord, self._box)
@@ -437,7 +433,7 @@ cdef class CellList:
         )
         # These have to be narrowed down in the next step
         # using the Euclidian distance
-        
+
         # Filter all indices from all_indices
         # where squared distance is smaller than squared radius
         # Using the squared distance is computationally cheaper than
@@ -464,23 +460,23 @@ cdef class CellList:
                         array_i += 1
             if array_i > max_array_length:
                 max_array_length = array_i
-        
+
         return self._post_process(
             np.asarray(indices)[:, :max_array_length],
             as_mask, is_multi_coord
         )
-    
-    
+
+
     @cython.boundscheck(False)
     @cython.wraparound(False)
     def get_atoms_in_cells(self, np.ndarray coord,
                            cell_radius=1, bint as_mask=False):
         """
         get_atoms_in_cells(coord, cell_radius=1, as_mask=False)
-        
+
         Find atoms with a maximum cell distance from given
         coordinates.
-        
+
         Instead of using the radius as maximum euclidian distance to the
         given coordinates,
         the radius is measured as the amount of cells:
@@ -489,7 +485,7 @@ cdef class CellList:
         that the atoms indices from this cell and the 8 surrounding
         cells are returned and so forth.
         This is more efficient than `get_atoms()`.
-        
+
         Parameters
         ----------
         coord : ndarray, dtype=float, shape=(3,) or shape=(m,3)
@@ -518,7 +514,7 @@ cdef class CellList:
         as_mask : bool, optional
             If true, the result is returned as boolean mask, instead
             of an index array.
-        
+
         Returns
         -------
         indices : ndarray, dtype=int32, shape=(p,) or shape=(m,p)
@@ -563,8 +559,8 @@ cdef class CellList:
             coord, cell_radius, is_multi_radius
         )
         return self._post_process(array_indices, as_mask, is_multi_coord)
-    
-    
+
+
     @cython.boundscheck(False)
     @cython.wraparound(False)
     def _get_atoms_in_cells(self,
@@ -573,7 +569,7 @@ cdef class CellList:
                             bint is_multi_radius):
         """
         Get the indices of atoms in `cell_radii` adjacency of `coord`.
-        
+
         Parameters
         ----------
         coord : ndarray, dtype=float32, shape=(n,3)
@@ -583,7 +579,7 @@ cdef class CellList:
         is_multi_radius : bool
             True indicates, that all values in `cell_radii` are the
             same.
-        
+
         Returns
         -------
         array_indices : ndarray, dtype=int32, shape=(m,p)
@@ -607,8 +603,8 @@ cdef class CellList:
         cdef int max_array_length \
             = self._find_adjacent_atoms(coord, array_indices, cell_radii)
         return array_indices[:, :max_array_length]
-    
-    
+
+
     @cython.boundscheck(False)
     @cython.wraparound(False)
     cdef int _find_adjacent_atoms(self,
@@ -618,7 +614,7 @@ cdef class CellList:
         """
         This method fills the given empty index array
         with actual indices of adjacent atoms.
-    
+
         Since the length of 'indices' (second dimension) is
         the worst case assumption, this method returns the actual
         required length, i.e. the highest length of all arrays
@@ -632,7 +628,7 @@ cdef class CellList:
         cdef int pos_i, array_i, cell_i
         cdef int max_array_length = 0
         cdef int cell_r
-        
+
         cdef ptr[:,:,:] cells = self._cells
         cdef int[:,:,:] cell_length = self._cell_length
 
@@ -663,7 +659,7 @@ cdef class CellList:
             if array_i > max_array_length:
                 max_array_length = array_i
         return max_array_length
-    
+
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
@@ -695,8 +691,8 @@ cdef class CellList:
                 return indices
             else:
                 return indices[0]
-    
-    
+
+
     @cython.initializedcheck(False)
     @cython.boundscheck(False)
     @cython.wraparound(False)
@@ -706,7 +702,7 @@ cdef class CellList:
         i[0] = <int>((x - self._min_coord[0]) / self._cellsize)
         j[0] = <int>((y - self._min_coord[1]) / self._cellsize)
         k[0] = <int>((z - self._min_coord[2]) / self._cellsize)
-    
+
     @cython.initializedcheck(False)
     @cython.boundscheck(False)
     @cython.wraparound(False)
@@ -718,7 +714,7 @@ cdef class CellList:
         if z < self._min_coord[2] or z > self._max_coord[2]:
             return False
         return True
-    
+
     @cython.initializedcheck(False)
     @cython.boundscheck(False)
     @cython.wraparound(False)
@@ -737,7 +733,7 @@ cdef class CellList:
                     break
                 matrix[i, index] = True
         return np.asarray(matrix, dtype=bool)
-    
+
     cdef inline bint _has_initialized_cells(self):
         # Memoryviews are not initialized on class creation
         # This method checks if the _cells memoryview was initialized
@@ -751,6 +747,20 @@ cdef class CellList:
             return False
 
 
+def _check_coord(coord):
+    """
+    Perform checks on validity of coordinates.
+    """
+    if coord.ndim != 2:
+        raise ValueError("Coordinates must have shape (n,3)")
+    if coord.shape[0] == 0:
+        raise ValueError("Coordinates must not be empty")
+    if coord.shape[1] != 3:
+        raise ValueError("Coordinates must have form (x,y,z)")
+    if np.isnan(coord).any():
+        raise ValueError("Coordinates contain NaN values")
+
+
 def _empty_result(as_mask):
     """
     Create return value for :func:`get_atoms()` and
@@ -774,11 +784,11 @@ def _prepare_vectorization(np.ndarray coord, radius, radius_dtype):
     radii/coordinates were given.
 
     The shapes before and after conversion are:
-       
+
        - coord: (3, ), radius: scalar -> coord: (1,3), radius: (1,)
        - coord: (n,3), radius: scalar -> coord: (n,3), radius: (n,)
        - coord: (n,3), radius: (n,  ) -> coord: (n,3), radius: (n,)
-    
+
     Thes resulting values have the same dimensionality for all cases and
     can be handeled uniformly by `get_atoms()` and
     `get_atoms_in_cells()`.
@@ -798,7 +808,7 @@ def _prepare_vectorization(np.ndarray coord, radius, radius_dtype):
         raise ValueError(
             f"Invalid shape for input coordinates"
         )
-    
+
     if isinstance(radius, np.ndarray):
         # Multiple radii
         # Check whether amount of coordinates match amount of radii
diff --git a/src/biotite/structure/compare.py b/src/biotite/structure/compare.py
index bdce1d7a..60d9cab5 100644
--- a/src/biotite/structure/compare.py
+++ b/src/biotite/structure/compare.py
@@ -9,11 +9,16 @@ comparing multiple structures with each other.
 
 __name__ = "biotite.structure"
 __author__ = "Patrick Kunzmann"
-__all__ = ["rmsd", "rmspd", "rmsf", "average"]
+__all__ = ["rmsd", "rmspd", "rmsf", "average", "lddt"]
 
+import collections.abc
+import warnings
 import numpy as np
 from biotite.structure.atoms import AtomArrayStack, coord
+from biotite.structure.celllist import CellList
+from biotite.structure.chains import get_chain_count, get_chain_positions
 from biotite.structure.geometry import index_distance
+from biotite.structure.residues import get_residue_count, get_residue_positions
 from biotite.structure.util import vector_dot
 
 
@@ -242,6 +247,218 @@ def average(atoms):
         return mean_coords
 
 
+def lddt(
+    reference,
+    subject,
+    aggregation="all",
+    atom_mask=None,
+    partner_mask=None,
+    inclusion_radius=15,
+    distance_bins=(0.5, 1.0, 2.0, 4.0),
+    exclude_same_residue=True,
+    symmetric=False,
+):
+    """
+    Calculate the *local Distance Difference Test* (lDDT) score of a structure with
+    respect to its reference.
+    :footcite:`Mariani2013`
+
+    Parameters
+    ----------
+    reference : AtomArray
+        The reference structure.
+    subject : AtomArray or AtomArrayStack or ndarray, dtype=float, shape=(n,3) or shape=(m,n,3)
+        The structure(s) to evaluate with respect to `reference`.
+        The number of atoms must be the same as in `reference`.
+        Alternatively, coordinates can be provided directly as
+        :class:`ndarray`.
+    aggregation : {'all', 'chain', 'residue', 'atom'} or ndarray, shape=(n,), dtype=int, optional
+        Defines on which scale the lDDT score is calculated.
+
+        - `'all'`: The score is computed over all contacts.
+        - `'chain'`: The score is calculated for each chain separately.
+        - `'residue'`: The score is calculated for each residue separately.
+        - `'atom'`: The score is calculated for each atom separately.
+
+        Alternatively, an array of aggregation bins can be provided, i.e. each contact
+        is assigned to the corresponding bin.
+    atom_mask : ndarray, shape=(n,), dtype=bool, optional
+        If given, the contacts are only computed for the masked atoms.
+        Atoms excluded by the mask do not have any contacts and their *lDDT* would
+        be NaN in case of ``aggregation="atom"``.
+        Providing this mask can significantly speed up the computation, if
+        only for certain chains/residues/atoms the *lDDT* is of interest.
+    partner_mask : ndarray, shape=(n,), dtype=bool, optional
+        If given, only contacts **to** the masked atoms are considered.
+        While `atom_mask` does not alter the *lDDT* for the masked atoms,
+        `partner_mask` does, as for each atom only the masked atoms are considered
+        as potential contact partners.
+    inclusion_radius : float, optional
+        Pairwise atom distances are considered within this radius in `reference`.
+    distance_bins : list of float, optional
+        The distance bins for the score calculation, i.e if a distance deviation is
+        within the first bin, the score is 1, if it is outside all bins, the score is 0.
+    exclude_same_residue : bool, optional
+        If set to False, distances between atoms of the same residue are also
+        considered.
+        By default, only atom distances between different residues are considered.
+    symmetric : bool, optional
+        If set to true, the *lDDT* score is computed symmetrically.
+        This means both contacts found in the `reference` and `subject` structure are
+        considered.
+        Hence the score is independent of which structure is given as `reference` and
+        `subject`.
+        Note that in this case `subject` must be an :class:`AtomArray` as well.
+        By default, only contacts in the `reference` are considered.
+
+    Returns
+    -------
+    lddt : float or ndarray, dtype=float
+        The lDDT score for each model and aggregation bin.
+        The shape depends on `subject` and `aggregation`:
+        If `subject` is an :class:`AtomArrayStack` (or equivalent coordinate
+        :class:`ndarray`), a dimension depicting each model is added.
+        if `aggregation` is not ``'all'``, a second dimension with the length equal to
+        the number of aggregation bins is added (i.e. number of chains, residues, etc.).
+        If both, an :class:`AtomArray` as `subject` and ``aggregation='all'`` is passed,
+        a float is returned.
+
+    Notes
+    -----
+    The lDDT score measures how well the pairwise atom distances in a model match the
+    corresponding distances in a reference.
+    Hence, like :func:`rmspd()` it works superimposition-free, but instead of capturing
+    the global deviation, only the local environment within the `inclusion_radius` is
+    considered.
+
+    Note that by default, also hydrogen atoms are considered in the distance
+    calculation.
+    If this is undesired, the hydrogen atoms can be removed prior to the calculation.
+
+    References
+    ----------
+
+    .. footbibliography::
+
+    Examples
+    --------
+
+    Calculate the global lDDT of all models to the first model:
+
+    >>> reference = atom_array_stack[0]
+    >>> subject = atom_array_stack[1:]
+    >>> print(lddt(reference, subject))
+    [0.799 0.769 0.792 0.836 0.799 0.752 0.860 0.769 0.825 0.777 0.760 0.787
+     0.790 0.783 0.804 0.842 0.769 0.797 0.757 0.852 0.811 0.786 0.805 0.755
+     0.734 0.794 0.771 0.778 0.842 0.772 0.815 0.789 0.828 0.750 0.826 0.739
+     0.760]
+
+    Calculate the residue-wise lDDT for a single model:
+
+    >>> subject = atom_array_stack[1]
+    >>> print(lddt(reference, subject, aggregation="residue"))
+    [0.599 0.692 0.870 0.780 0.830 0.881 0.872 0.658 0.782 0.901 0.888 0.885
+     0.856 0.795 0.847 0.603 0.895 0.878 0.871 0.789]
+
+    As example for custom aggregation, calculate the lDDT for each chemical element:
+
+    >>> unique_elements = np.unique(reference.element)
+    >>> element_bins = np.array(
+    ...     [np.where(unique_elements == element)[0][0] for element in reference.element]
+    ... )
+    >>> element_lddt = lddt(reference, subject, aggregation=element_bins)
+    >>> for element, lddt_for_element in zip(unique_elements, element_lddt):
+    ...     print(f"{element}: {lddt_for_element:.3f}")
+    C: 0.837
+    H: 0.770
+    N: 0.811
+    O: 0.808
+
+    If the reference structure has more atoms resolved than the subject structure,
+    the missing atoms can be indicated with *NaN* values:
+
+    >>> reference = atom_array_stack[0]
+    >>> subject = atom_array_stack[1].copy()
+    >>> # Simulate the situation where the first residue is missing in the subject
+    >>> subject.coord[subject.res_id == 1] = np.nan
+    >>> global_lddt = lddt(reference, subject)
+    >>> print(f"{global_lddt:.3f}")
+    0.751
+    """
+    reference_coord = coord(reference)
+    subject_coord = coord(subject)
+    if subject_coord.shape[-2] != reference_coord.shape[-2]:
+        raise IndexError(
+            f"The given reference has {reference_coord.shape[-2]} atoms, but the "
+            f"subject has {subject_coord.shape[-2]} atoms"
+        )
+
+    contacts = _find_contacts(
+        reference, atom_mask, partner_mask, inclusion_radius, exclude_same_residue
+    )
+    if symmetric:
+        subject_contacts = _find_contacts(
+            subject, atom_mask, partner_mask, inclusion_radius, exclude_same_residue
+        )
+        contacts = np.concatenate((contacts, subject_contacts), axis=0)
+        # Adding additional contacts may introduce duplicates between the existing and
+        # new ones -> filter them out
+        contacts = np.unique(contacts, axis=0)
+    if (
+        isinstance(aggregation, str)
+        and aggregation == "all"
+        and atom_mask is None
+        and partner_mask is None
+    ):
+        # Remove duplicate pairs as each pair appears twice
+        # (if i is in threshold distance to j, j is also in threshold distance to i)
+        # keep only the pair where i < j
+        # This improves performance due to less distances that need to be computed
+        # The assumption also only works when no atoms are masked
+        contacts = contacts[contacts[:, 0] < contacts[:, 1]]
+
+    reference_distances = index_distance(reference_coord, contacts)
+    subject_distances = index_distance(subject_coord, contacts)
+    deviations = np.abs(subject_distances - reference_distances)
+    distance_bins = np.asarray(distance_bins)
+    fraction_preserved_bins = np.count_nonzero(
+        deviations[..., np.newaxis] <= distance_bins[np.newaxis, :], axis=-1
+    ) / len(distance_bins)
+
+    # Aggregate the fractions over the desired level
+    if isinstance(aggregation, str) and aggregation == "all":
+        # Average over all contacts
+        return np.mean(fraction_preserved_bins, axis=-1)
+    else:
+        # A string is also a 'Sequence'
+        # -> distinguish between string and array, list, etc.
+        if isinstance(
+            aggregation, (np.ndarray, collections.abc.Sequence)
+        ) and not isinstance(aggregation, str):
+            return _average_over_indices(
+                fraction_preserved_bins,
+                bins=np.asarray(aggregation)[contacts[:, 0]],
+            )
+        elif aggregation == "chain":
+            return _average_over_indices(
+                fraction_preserved_bins,
+                bins=get_chain_positions(reference, contacts[:, 0]),
+                n_bins=get_chain_count(reference),
+            )
+        elif aggregation == "residue":
+            return _average_over_indices(
+                fraction_preserved_bins,
+                bins=get_residue_positions(reference, contacts[:, 0]),
+                n_bins=get_residue_count(reference),
+            )
+        elif aggregation == "atom":
+            return _average_over_indices(
+                fraction_preserved_bins, contacts[:, 0], reference.array_length()
+            )
+        else:
+            raise ValueError(f"Invalid aggregation level '{aggregation}'")
+
+
 def _sq_euclidian(reference, subject):
     """
     Calculate squared euclidian distance between atoms in two
@@ -272,3 +489,149 @@ def _sq_euclidian(reference, subject):
         )
     dif = subject_coord - reference_coord
     return vector_dot(dif, dif)
+
+
+def _to_sparse_indices(all_contacts):
+    """
+    Create tuples of contact indices from the :meth:`CellList.get_atoms()` return value.
+
+    In other words, they would mark the non-zero elements in a dense contact matrix.
+
+    Parameters
+    ----------
+    all_contacts : ndarray, dtype=int, shape=(m,n)
+        The contact indices as returned by :meth:`CellList.get_atoms()`.
+        Padded with -1, in the second dimension.
+        Dimension *m* marks the query atoms, dimension *n* marks the contact atoms.
+
+    Returns
+    -------
+    combined_indices : ndarray, dtype=int, shape=(l,2)
+        The contact indices.
+        Each column contains the query and contact atom index.
+    """
+    # Find rows where a query atom has at least one contact
+    non_empty_indices = np.where(np.any(all_contacts != -1, axis=1))[0]
+    # Take those rows and flatten them
+    contact_indices = all_contacts[non_empty_indices].flatten()
+    # For each row the corresponding query atom is the same
+    # Hence in the flattened form the query atom index is simply repeated
+    query_indices = np.repeat(non_empty_indices, all_contacts.shape[1])
+    combined_indices = np.stack([query_indices, contact_indices], axis=1)
+    # Remove the padding values
+    return combined_indices[contact_indices != -1]
+
+
+def _find_contacts(
+    atoms=None,
+    atom_mask=None,
+    partner_mask=None,
+    inclusion_radius=15,
+    exclude_same_residue=True,
+):
+    """
+    Find contacts between the atoms in the given structure.
+
+    Parameters
+    ----------
+    atoms : AtomArray
+        The structure to find the contacts for.
+    atom_mask : ndarray, shape=(n,), dtype=bool, optional
+        If given, the contacts are only computed for the masked atoms.
+        Atoms excluded by the mask do not have any contacts and their *lDDT* would
+        be NaN in case of ``aggregation="atom"``.
+        Providing this mask can significantly speed up the computation, if
+        only for certain chains/residues/atoms the *lDDT* is of interest.
+    partner_mask : ndarray, shape=(n,), dtype=bool, optional
+        If given, only contacts **to** the masked atoms are considered.
+        While `atom_mask` does not alter the *lDDT* for the masked atoms,
+        `partner_mask` does, as for each atom only the masked atoms are considered
+        as potential contact partners.
+    inclusion_radius : float, optional
+        Pairwise atom distances are considered within this radius.
+    exclude_same_residue : bool, optional
+        If set to False, distances between atoms of the same residue are also
+        considered.
+        By default, only atom distances between different residues are considered.
+
+    Returns
+    -------
+    contacts : ndarray, shape=(n,2), dtype=int
+        The array of contacts.
+        Each element represents a pair of atom indices that are in contact.
+    """
+    coords = coord(atoms)
+    selection = ~np.isnan(coords).any(axis=-1)
+    if partner_mask is not None:
+        selection &= partner_mask
+    # Use a cell list to find atoms within inclusion radius in O(n) time complexity
+    cell_list = CellList(coords, inclusion_radius, selection=selection)
+    # Pairs of indices for atoms within the inclusion radius
+    if atom_mask is None:
+        all_contacts = cell_list.get_atoms(coords, inclusion_radius)
+    else:
+        filtered_contacts = cell_list.get_atoms(coords[atom_mask], inclusion_radius)
+        # Map the contacts for the masked atoms to the original coordinates
+        # Rows that were filtered out by the mask are fully padded with -1
+        # consistent with the padding of `get_atoms()`
+        all_contacts = np.full(
+            (coords.shape[0], filtered_contacts.shape[-1]),
+            -1,
+            dtype=filtered_contacts.dtype,
+        )
+        all_contacts[atom_mask] = filtered_contacts
+    # Convert into pairs of indices
+    contacts = _to_sparse_indices(all_contacts)
+
+    if exclude_same_residue:
+        # Find the index of the residue for each atom
+        residue_indices = get_residue_positions(atoms, contacts.flatten()).reshape(
+            contacts.shape
+        )
+        # Remove contacts between atoms of the same residue
+        contacts = contacts[residue_indices[:, 0] != residue_indices[:, 1]]
+    else:
+        # In any case self-contacts should not be considered
+        contacts = contacts[contacts[:, 0] != contacts[:, 1]]
+    return contacts
+
+
+def _average_over_indices(values, bins, n_bins=None):
+    """
+    For each unique index in `bins`, average the corresponding values in `values`.
+
+    Based on
+    https://stackoverflow.com/questions/79140661/how-to-sum-values-based-on-a-second-index-array-in-a-vectorized-manner
+
+    Parameters
+    ----------
+    values : ndarray, shape=(..., n)
+        The values to average.
+    bins : ndarray, shape=(n,) dtype=int
+        Associates each value from `values` with a bin.
+    n_bins : int
+        The total number of bins.
+        This is necessary as the some bin in `bins`may be empty.
+        By default the number of bins is determined from `bins`.
+
+    Returns
+    -------
+    averaged : ndarray, shape=(..., k)
+        The averaged values.
+        *k* is the maximum value in `bins` + 1.
+    """
+    if n_bins is None:
+        n_elements_per_bin = np.bincount(bins)
+        n_bins = len(n_elements_per_bin)
+    else:
+        n_elements_per_bin = np.bincount(bins, minlength=n_bins)
+    # The last dimension is replaced by the number of bins
+    # Broadcasting in 'np.add.at()' requires the replaced dimension to be the first
+    aggregated = np.zeros((n_bins, *values.shape[:-1]), dtype=values.dtype)
+    np.add.at(aggregated, bins, np.swapaxes(values, 0, -1))
+    # If an atom has no contacts, the corresponding value is NaN
+    # This result is expected, hence the warning is ignored
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        # Bring the bin dimension into the last dimension again
+        return np.swapaxes(aggregated, 0, -1) / n_elements_per_bin


Test patch:
diff --git a/tests/structure/test_compare.py b/tests/structure/test_compare.py
index 1895c8ba..70ea3351 100644
--- a/tests/structure/test_compare.py
+++ b/tests/structure/test_compare.py
@@ -2,11 +2,12 @@
 # under the 3-Clause BSD License. Please see 'LICENSE.rst' for further
 # information.
 
+import itertools
 from os.path import join
 import numpy as np
 import pytest
 import biotite.structure as struc
-import biotite.structure.io as strucio
+import biotite.structure.io.pdbx as pdbx
 from tests.util import data_dir
 
 
@@ -36,21 +37,31 @@ def test_rmsf(stack, as_coord):
     )
 
 
+@pytest.fixture(scope="module")
+def models():
+    """
+    Get a short multi-model structure
+    """
+    pdbx_file = pdbx.BinaryCIFFile.read(join(data_dir("structure"), "1l2y.bcif"))
+    return pdbx.get_structure(pdbx_file)
+
+
 @pytest.fixture
-def load_stack_superimpose():
-    stack = strucio.load_structure(join(data_dir("structure"), "1l2y.bcif"))
-    # Superimpose with first frame
-    bb_mask = struc.filter_peptide_backbone(stack[0])
-    supimp, _ = struc.superimpose(stack[0], stack, atom_mask=bb_mask)
-    return stack, supimp
+def superimposed_models(models):
+    """
+    Models superimposed on the first model
+    """
+    bb_mask = struc.filter_peptide_backbone(models[0])
+    supimp, _ = struc.superimpose(models[0], models, atom_mask=bb_mask)
+    return models, supimp
 
 
-def test_rmsd_gmx(load_stack_superimpose):
+def test_rmsd_gmx(superimposed_models):
     """
     Comparison of RMSD values computed with Biotite with results
     obtained from GROMACS 2021.5.
     """
-    stack, supimp = load_stack_superimpose
+    stack, supimp = superimposed_models
     rmsd = struc.rmsd(stack[0], supimp) / 10
 
     # Gromacs RMSDs -> Without mass-weighting:
@@ -102,12 +113,12 @@ def test_rmsd_gmx(load_stack_superimpose):
     assert np.allclose(rmsd, rmsd_gmx, atol=1e-03)
 
 
-def test_rmspd_gmx(load_stack_superimpose):
+def test_rmspd_gmx(superimposed_models):
     """
     Comparison of the RMSPD computed with Biotite with results
     obtained from GROMACS 2021.5.
     """
-    stack, _ = load_stack_superimpose
+    stack, _ = superimposed_models
     rmspd = struc.rmspd(stack[0], stack) / 10
 
     # Gromacs RMSDist:
@@ -159,12 +170,12 @@ def test_rmspd_gmx(load_stack_superimpose):
     assert np.allclose(rmspd, rmspd_gmx, atol=1e-03)
 
 
-def test_rmsf_gmx(load_stack_superimpose):
+def test_rmsf_gmx(superimposed_models):
     """
     Comparison of RMSF values computed with Biotite with results
     obtained from GROMACS 2021.5.
     """
-    stack, supimp = load_stack_superimpose
+    stack, supimp = superimposed_models
     ca_mask = (stack[0].atom_name == "CA") & (stack[0].element == "C")
     rmsf = struc.rmsf(struc.average(supimp[:, ca_mask]), supimp[:, ca_mask]) / 10
 
@@ -196,3 +207,199 @@ def test_rmsf_gmx(load_stack_superimpose):
     )
 
     assert np.allclose(rmsf, rmsf_gmx, atol=1e-02)
+
+
+@pytest.mark.parametrize(
+    "multi_model, as_coord", itertools.product([False, True], [False, True])
+)
+def test_lddt_perfect(models, multi_model, as_coord):
+    """
+    Check if the lDDT of a structure with itself as reference is 1.0
+    """
+    reference = models[0]
+    if multi_model:
+        subject = models[0:1]
+    else:
+        subject = models[0]
+
+    if as_coord:
+        reference = reference.coord
+        subject = subject.coord
+        # Coordinates can only be used if 'exclude_same_residue' is False
+        exclude_same_residue = False
+    else:
+        exclude_same_residue = True
+
+    lddt = struc.lddt(reference, subject, exclude_same_residue=exclude_same_residue)
+
+    if multi_model:
+        assert lddt.tolist() == [1.0]
+    else:
+        assert lddt == 1.0
+
+
+def test_lddt_consistency(models):
+    """
+    Check if the computed lDDT via :func:`lddt` is the same as the reference value from
+    the https://swissmodel.expasy.org/assess web server.
+    """
+    # As each model needs to be uploaded separately, only the first 8 models were run
+    REFERENCE_LDDT = [
+        1.00,
+        0.86,
+        0.83,
+        0.85,
+        0.88,
+        0.85,
+        0.81,
+        0.89,
+    ]
+
+    # The web server computes the lDDT without hydrogen atoms
+    models = models[:, models.element != "H"]
+
+    reference = models[0]
+    subject = models[: len(REFERENCE_LDDT)]
+    lddt = struc.lddt(reference, subject)
+
+    assert lddt.tolist() == pytest.approx(REFERENCE_LDDT, abs=1e-2)
+
+
+@pytest.mark.parametrize("multi_model", [False, True])
+def test_lddt_custom_aggregation(models, multi_model):
+    """
+    Check if custom `aggregation` in :func:`lddt` works by giving a custom aggregation
+    that simply aggregates all values, i.e. it should give the same result as the
+    `all` aggregation.
+    """
+    reference = models[0]
+    if multi_model:
+        subject = models
+    else:
+        subject = models[1]
+    # Every atom goes into the same aggregation bin
+    aggregation_bins = np.zeros(models.array_length(), dtype=int)
+
+    ref_lddt = struc.lddt(reference, subject, aggregation="all")
+
+    test_lddt = struc.lddt(reference, subject, aggregation=aggregation_bins)
+
+    # There is only one aggregation bin -> only one lDDT value (for each value)
+    assert test_lddt.shape[-1] == 1
+    assert test_lddt[..., 0].tolist() == ref_lddt.tolist()
+
+
+@pytest.mark.parametrize(
+    "multi_model, aggregation",
+    itertools.product([False, True], ["chain", "residue", "atom"]),
+)
+def test_lddt_aggregation_levels(models, multi_model, aggregation):
+    """
+    Check if the predefined aggregation levels :func:`lddt()` return the expected shape.
+    Furthermore, the average of each aggregated lDDT values should be similar to the
+    global lDDT value (not exactly, as the global lDDT is weighted by the number of
+    contacts in each residue/chain).
+    """
+    ABS_ERROR = 0.02
+
+    reference = models[0]
+    if multi_model:
+        subject = models
+    else:
+        subject = models[1]
+
+    expected_shape = ()
+    if multi_model:
+        expected_shape = expected_shape + (subject.stack_depth(),)
+    if aggregation == "chain":
+        expected_shape = expected_shape + (struc.get_chain_count(subject),)
+    elif aggregation == "residue":
+        expected_shape = expected_shape + (struc.get_residue_count(subject),)
+    elif aggregation == "atom":
+        expected_shape = expected_shape + (models.array_length(),)
+
+    lddt = struc.lddt(reference, subject, aggregation=aggregation)
+    all_lddt = struc.lddt(reference, subject, aggregation="all")
+
+    assert lddt.shape == expected_shape
+    if multi_model:
+        assert np.mean(lddt, axis=-1).tolist() == pytest.approx(
+            all_lddt.tolist(), abs=ABS_ERROR
+        )
+    else:
+        assert np.mean(lddt) == pytest.approx(all_lddt, abs=ABS_ERROR)
+
+
+# Dividing by zero contacts for masked atoms would normally raise a warning
+@pytest.mark.filterwarnings("error:.*invalid value encountered in divide")
+@pytest.mark.parametrize("seed", range(10))
+def test_lddt_mask(models, seed):
+    """
+    Using the same `subject_mask` and `partner_mask` should give the same score as
+    using the same mask on the structure directly, before it is passed to
+    :func:`lddt()`.
+    """
+    rng = np.random.default_rng(seed)
+    mask = rng.choice([False, True], size=models.array_length())
+
+    reference = models[0]
+    subject = models[1]
+
+    # Do not aggregate to make the test more strict
+    ref_lddt = struc.lddt(reference[mask], subject[mask], aggregation="atom")
+
+    test_lddt = struc.lddt(
+        reference, subject, aggregation="atom", atom_mask=mask, partner_mask=mask
+    )
+
+    # The subject mask defines for which atoms the contact partners should be identified
+    # In consequence, for all atoms that are masked out, the lDDT should be NaN
+    assert np.isnan(test_lddt[..., ~mask]).all()
+    test_lddt = test_lddt[mask]
+    assert test_lddt.tolist() == pytest.approx(ref_lddt.tolist())
+
+
+def test_custom_lddt_symmetric(models):
+    """
+    Check that in :func:`lddt()` with ``symmetric=True`` the *lDDT* score is independent
+    of the order of the structures.
+    """
+    model1 = models[0]
+    model2 = models[1]
+
+    # Do not aggregate to make the test more strict
+    model1_symmetric_lddt = struc.lddt(model1, model2, "atom", symmetric=True)
+    model2_symmetric_lddt = struc.lddt(model2, model1, "atom", symmetric=True)
+
+    # The order of the structures should not matter
+    assert model1_symmetric_lddt.tolist() == pytest.approx(
+        model2_symmetric_lddt.tolist()
+    )
+
+
+@pytest.mark.parametrize(
+    "what_is_nan, ref_lddt",
+    [
+        # No reference coordinates (except one atom)
+        # -> set of reference contacts is empty
+        ("reference", np.nan),
+        # No subject coordinates (except one atom)
+        # -> no recovered contacts
+        ("subject", 0.0),
+    ],
+)
+def test_lddt_with_missing_atoms(models, what_is_nan, ref_lddt):
+    """
+    Check if NaN values in the coordinates are handled correctly.
+    """
+    reference = models[0].copy()
+    subject = models[1].copy()
+    if what_is_nan == "reference":
+        reference.coord[1:, :] = np.nan
+    elif what_is_nan == "subject":
+        subject.coord[1:, :] = np.nan
+
+    test_lddt = struc.lddt(reference, subject)
+
+    # Use 'np.array_equal()' to handle equality of NaN values
+    assert np.array_equal(test_lddt, ref_lddt, equal_nan=True)