diff --git a/src/biotite/structure/tm.py b/src/biotite/structure/tm.py
index 508af70a3..48fcd2132 100644
--- a/src/biotite/structure/tm.py
+++ b/src/biotite/structure/tm.py
@@ -533,8 +533,7 @@ def _tm_score(distance, reference_length):
     tm_score : float or ndarray
         The TM score for the given distances.
     """
-    d0 = max(_d0(reference_length), _D0_MIN)
-    return 1 / (1 + (distance / d0) ** 2)
+    return 1 / (1 + (distance / _d0(reference_length)) ** 2)
 
 
 def _d0(reference_length):
@@ -552,7 +551,12 @@ def _d0(reference_length):
         The :math:`d_0` threshold.
     """
     # Constants taken from Zhang2004
-    return 1.24 * (reference_length - 15) ** (1 / 3) - 1.8
+    return max(
+        # Avoid complex solutions -> clip to positive values
+        # For short sequence lengths _D0_MIN takes precedence anyway
+        1.24 * max((reference_length - 15), 0) ** (1 / 3) - 1.8,
+        _D0_MIN,
+    )
 
 
 def _get_reference_length(user_parameter, reference_length, subject_length):
diff --git a/tests/structure/test_tm.py b/tests/structure/test_tm.py
index c02b7b486..e39bddcef 100644
--- a/tests/structure/test_tm.py
+++ b/tests/structure/test_tm.py
@@ -158,6 +158,23 @@ def test_superimpose_consistency(fixed_pdb_id, mobile_pdb_id, ref_tm_score):
     )
 
 
+def test_short_structure():
+    """
+    Check that the :func:`tm_score()` and :func:`superimpose_structural_homologs()` work
+    for structures shorter than 15 residues, because in that case we got
+
+    ``TypeError: '>' not supported between instances of 'float' and 'complex'``
+
+    before.
+    """
+    pdbx_file = pdbx.BinaryCIFFile.read(Path(data_dir("structure")) / "1l2y.bcif")
+    atoms = pdbx.get_structure(pdbx_file, model=1)
+    atoms = atoms[atoms.res_id < 15]
+
+    _, _, fixed_i, mobile_i = struc.superimpose_structural_homologs(atoms, atoms)
+    struc.tm_score(atoms, atoms, fixed_i, mobile_i)
+
+
 def _transform_random_affine(coord, rng):
     coord = struc.translate(coord, rng.uniform(low=0, high=10, size=3))
     coord = struc.rotate(coord, rng.uniform(low=0, high=2 * np.pi, size=3))