diff --git a/pysisyphus/intcoords/setup.py b/pysisyphus/intcoords/setup.py
index 98091b51c..90eff0cfa 100644
--- a/pysisyphus/intcoords/setup.py
+++ b/pysisyphus/intcoords/setup.py
@@ -435,6 +435,7 @@ def get_hydrogen_bond_inds(atoms, coords3d, bond_inds, logger=None):
     hydrogen_inds = [i for i, a in enumerate(atoms) if a.lower() == "h"]
     x_inds = [i for i, a in enumerate(atoms) if a.lower() in "n o f p s cl".split()]
     hydrogen_bond_inds = list()
+    bond_sets = [set(map(int, bond)) for bond in bond_inds]
     for h_ind, x_ind in it.product(hydrogen_inds, x_inds):
         as_set = set((h_ind, x_ind))
         if as_set not in tmp_sets:
@@ -450,8 +451,13 @@ def get_hydrogen_bond_inds(atoms, coords3d, bond_inds, logger=None):
             distance = Stretch._calculate(coords3d, (h_ind, y_ind))
             vdw_rad_sum = VDW_RADII["h"] + VDW_RADII[y_atom]
             angle = Bend._calculate(coords3d, (x_ind, h_ind, y_ind))
-            if (1.1 * cov_rad_sum < distance < 0.9 * vdw_rad_sum) and (
-                angle > np.pi / 2
+            if (
+                (1.1 * cov_rad_sum < distance < 0.9 * vdw_rad_sum)
+                and (angle > np.pi / 2)
+                # Avoid adding hydrogen bonds that are already in bond_inds.
+                # This can happen when the user manually defined a bond, that
+                # is actually a hydrogen bond.
+                and {h_ind, y_ind} not in bond_sets
             ):
                 hydrogen_bond_inds.append((h_ind, y_ind))
                 log(
diff --git a/tests/test_intcoords/test_intcoords.py b/tests/test_intcoords/test_intcoords.py
index 8d2b8dd06..c9a44876d 100644
--- a/tests/test_intcoords/test_intcoords.py
+++ b/tests/test_intcoords/test_intcoords.py
@@ -292,7 +292,9 @@ def test_hydrogen_bonds():
     geom = geom_loader("lib:hydrogen_bond_test.xyz", coord_type="redund")
     typed_prims = geom.internal.typed_prims
     h_bonds = [tp for tp in typed_prims if tp[0] == PrimTypes.HYDROGEN_BOND]
-    h_bond_ref = set([frozenset(inds) for inds in ((6, 27), (43, 56), (15, 42))])
+    # (6, 27) was included below earlier, but this is already designated as a real bond
+    # and is not added again as hydrogen bond.
+    h_bond_ref = set([frozenset(inds) for inds in ((43, 56), (15, 42))])
     assert set([frozenset(inds) for _, *inds in h_bonds]) == h_bond_ref