diff --git a/src/clustering.jl b/src/clustering.jl
index 5459145..4fc1103 100644
--- a/src/clustering.jl
+++ b/src/clustering.jl
@@ -1235,6 +1235,15 @@ function _split_this_sbu!(toremove, graph, k, types, stopiftype, sbus)
             length(othersbu) == 1 && types[othersbu[1].v] === stopiftype && return nothing
         end
     end
+    if stopiftype === :O
+        # skip splitting if O has too many H-bonds (like in clathrate hydrates)
+        numhbonds = 0
+        neighsO = neighbors(graph, k)
+        for (x, _) in neighsO
+            numhbonds += types[x] === :H
+        end
+        length(neighsO) - numhbonds ≤ 2 && return nothing
+    end
     push!(toremove, k)
     neighs = reverse(neighbors(graph, k))
     n = length(neighs)
diff --git a/src/input.jl b/src/input.jl
index 63dc859..ed3367f 100644
--- a/src/input.jl
+++ b/src/input.jl
@@ -1274,6 +1274,59 @@ function remove_homoatomic_bonds!(graph::PeriodicGraph{D}, types, targets, reduc
     reduce_homometallic && _remove_homometallic_bonds!(graph, types, metallics)
 end
 
+function add_H_bonds!(graph::PeriodicGraph3D, types, cell::Cell, pos, opts::Options)
+    bonded_H = Tuple{Int,PeriodicVertex3D}[]
+    acceptors = Int[]
+    for (i, tH) in enumerate(types)
+        if tH in (:O, :N, :F, :Cl, :S)
+            push!(acceptors, i)
+        elseif tH === :H
+            for x in neighbors(graph, i)
+                v = first(x)
+                t = types[v]
+                if t in (:O, :N, :F, :Cl, :S)
+                    push!(bonded_H, (i, x))
+                    break
+                end
+            end
+        end
+    end
+    if isempty(bonded_H)
+        @ifwarn @warn "Could not find any H bonded to O, N, F, Cl or S although option H_bonds is set"
+        return
+    end
+    mat = Float64.(cell.mat)
+    pd2 = PeriodicDistance2(mat)
+    ofs = MVector{3,Int}(undef)
+    dmax2 = opts.Hbonds_dmax^2
+    θmax = opts.Hbonds_θmax
+    nmax = opts.Hbonds_nmax
+    accepted = Tuple{PeriodicVertex3D,Float64}[]
+    for (uH, (uD, ofsD)) in bonded_H
+        pH = pos[uH]
+        pD = pos[uD] .+ ofsD
+        empty!(accepted)
+        for uA in acceptors
+            uA == uD && continue
+            pA = pos[uA]
+            d = pd2(pH, pA; ofs)
+            1.0^2 ≤ d ≤ dmax2 || continue
+            θ = angle(mat*(pD .- (pA .+ ofs)), mat*(pD .- pH))
+            if θ < θmax
+                push!(accepted, (PeriodicVertex3D(uA, SVector{3,Int}(ofs)), d))
+            end
+        end
+        if length(accepted) > nmax
+            sort!(accepted; by=last)
+            resize!(accepted, nmax)
+        end
+        for (x, _) in accepted
+            add_edge!(graph, uH, x)
+        end
+    end
+    nothing
+end
+
 
 function parse_as_cif(cif::CIF, options::Options, name::String)
     guessed_bonds = false
@@ -1454,6 +1507,8 @@ function finalize_checks(cell::Cell, pos::Vector{SVector{3,Float64}}, types::Vec
         remove_homoatomic_bonds!(graph, types, options.ignore_homoatomic_bonds, options.reduce_homometallic_bonds)
     end
 
+    options.Hbonds && add_H_bonds!(graph, types, cell, pos, options)
+
     if isempty(attributions)
         crystalnothing = Crystal{Nothing}(cell, types, pos, graph, options)
         export_default(crystalnothing, "input", name, options.export_input)
diff --git a/src/options.jl b/src/options.jl
index 754151e..8443673 100644
--- a/src/options.jl
+++ b/src/options.jl
@@ -342,6 +342,13 @@ string is equivalent to `false`.
 - `max_polyhedron_radius`: an integer specifying the maximum number of bonds between two
   corners of the coordination polyhedron built for the [`Clustering.PE`](@ref Clustering) option.
   Default is 4.
+- `Hbonds`: set to `true` to include hydrogen bonds. Default is false.
+- `Hbonds_dmax`: the maximum length of a hydrogen bond. Only used if `Hbonds` is set.
+  Default is 2.5 Å.
+- `Hbonds_θmax`: the maximum angle of a hydrogen bond. Only used if `Hbonds` is set.
+  Default is 30°.
+- `Hbonds_nmax`: the maximum number of hydrogen bond per hydrogen. Only used if `Hbonds` is
+  set. Default is 1.
 
 ## Miscellaneous options
 These boolean options have a default value that may be determined by [`Bonding`](@ref),
@@ -371,7 +378,7 @@ These boolean options have a default value that may be determined by [`Bonding`]
 - `premerge_metalbonds`: when a periodic metallic SBU is detected, cluster together bonded
   metal atoms of the same kind before splitting the SBU.
 - `split_O_vertex`: if a vertex is composed of a single O, remove it and bond together all of
-  its neighbors. Default is true.
+  its neighbors, unless removing its hydrogen bonds would make it bivalent. Default is true.
 - `unify_sbu_decomposition`: apply the same rule to decompose both periodic and finite SBUs.
   Default is false.
 - `force_warn`: force printing warning and information even during `..._dataset` function
@@ -409,6 +416,10 @@ struct Options
     export_trimmed::String
     force_warn::Bool
     label_for_type::Bool
+    Hbonds::Bool
+    Hbonds_dmax::Float64
+    Hbonds_θmax::Float64
+    Hbonds_nmax::Int
 
     # Clustering options
     clusterings::Vector{_Clustering}
@@ -456,6 +467,10 @@ struct Options
                        export_trimmed=false,
                        force_warn=false,
                        label_for_type=false,
+                       Hbonds=false,
+                       Hbonds_dmax=2.5,
+                       Hbonds_θmax=30.0,
+                       Hbonds_nmax=1,
                        clusterings=_Clustering[Clustering.Auto],
                        bond_adjacent_sbus=false,
                        ignore_metal_cluster_bonds=nothing,
@@ -518,6 +533,10 @@ struct Options
             _export_trimmed,
             force_warn,
             label_for_type,
+            Hbonds,
+            Hbonds_dmax,
+            Hbonds_θmax,
+            Hbonds_nmax,
             clusterings,
             bond_adjacent_sbus,
             ignore_metal_cluster_bonds,