Extract data on harmonic bonds

extractdata/bonding.jl

@@ -21,28 +21,33 @@ function parsestring(str)
     return String(str[2:end-1])
 end
 
-function parsexmlline(f, line, tag, keyname)
+function parsexmlline(f, line, tag, keynames...; skip=())
     @assert startswith(line, "<$tag ")
     @assert endswith(line, "/>")
     kv = split(strip(line[length(tag)+2:end-2]), ' ')
-    key = ""
+    key = String[]
     vals = Pair{Symbol,Any}[]
     for kvp in kv
         k, v = split(kvp, '=')
-        if k == keyname
-            key = parsestring(v)
+        k ∈ skip && continue
+        if k ∈ keynames
+            push!(key, parsestring(v))
         else
             push!(vals, Symbol(k) => f(k, v))
         end
     end
+    @assert length(key) == length(keynames)
+    key = length(key) == 1 ? only(key) : (key...,)
     return key => (; vals...)
 end
 
-atomtypes, residues = open("protein.ff14SB.xml", "r") do io
+atomtypes, residues, bondlengths, bondangles = open("protein.ff14SB.xml", "r") do io
     line = readline(io)
     @assert line == "<ForceField>"
     atomtypes = Dict{String, @NamedTuple{element::String, mass::Float32, name::String}}()
     residues = Dict{String, @NamedTuple{atoms::Dict{String, @NamedTuple{charge::Float32, type::String}}, bonds::Vector{Tuple{String,String}}, externalbonds::Vector{String}}}()
+    harmonicbonds = Dict{Tuple{String,String}, @NamedTuple{length::Float32}}()
+    harmonicangles = Dict{Tuple{String,String,String}, @NamedTuple{angle::Float32}}()
     parsexmblock(io, "</ForceField>") do line
         if line == "<AtomTypes>"
             parsexmblock(io, "</AtomTypes>") do line
@@ -92,9 +97,31 @@ atomtypes, residues = open("protein.ff14SB.xml", "r") do io
                     error("Unknown Residues line $line")
                 end
             end
+        elseif line == "<HarmonicBondForce>"
+            parsexmblock(io, "</HarmonicBondForce>") do line
+                push!(harmonicbonds, parsexmlline(line, "Bond", "type1", "type2"; skip=("k",)) do k, v
+                    if k == "length"
+                        return parse(Float32, v[2:end-1])
+                    else
+                        error("Unknown Bond key $k")
+                    end
+                end)
+            end
+        elseif line == "<HarmonicAngleForce>"
+            parsexmblock(io, "</HarmonicAngleForce>") do line
+                push!(harmonicangles, parsexmlline(line, "Angle", "type1", "type2", "type3"; skip=("k",)) do k, v
+                    if k == "k"
+                        return parse(Float32, v[2:end-1])  # strip the quotes
+                    elseif k == "angle"
+                        return parse(Float32, v[2:end-1])  # strip the quotes
+                    else
+                        error("Unknown Angle key $k")
+                    end
+                end)
+            end
         end
     end
-    atomtypes, residues
+    atomtypes, residues, Dict(k => 10 * v.length for (k, v) in harmonicbonds), Dict(k => v.angle for (k, v) in harmonicangles)
 end
 
 open(joinpath(dirname(@__DIR__), "src", "bonding.jl"), "w") do io
@@ -117,4 +144,18 @@ open(joinpath(dirname(@__DIR__), "src", "bonding.jl"), "w") do io
         println(io, "               externalbonds = ", v.externalbonds, "),")
     end
     println(io, ")")
+
+    println(io, "\nconst bondlengths = Dict{Tuple{String,String}, Float32}(")
+    hb = sort!(collect(bondlengths); by=first)
+    for pr in hb
+        println(io, "    ", pr, ',')
+    end
+    println(io, ")")
+
+    println(io, "\nconst bondangles = Dict{Tuple{String,String,String}, Float32}(")
+    ha = sort!(collect(bondangles); by=first)
+    for pr in ha
+        println(io, "    ", pr, ',')
+    end
+    println(io, ")")
 end