From 880f321b8d4e4a9a33a3627e5360a3cfd70729cb Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Wed, 29 Jul 2026 07:22:46 -0500 Subject: [PATCH 1/2] Extract Amber LJ params This parses the Lennard-Jones parameters from Amber's `protein.ff14SB.xml` file, storing the result in `atomtypes`. --- NEWS.md | 5 +++ Project.toml | 2 +- docs/src/api.md | 4 ++ extractdata/bonding.jl | 73 ++++++++++++++++++++++++++++-- src/bonding.jl | 100 ++++++++++++++++++++++++++--------------- test/runtests.jl | 19 ++++++++ 6 files changed, 164 insertions(+), 39 deletions(-) diff --git a/NEWS.md b/NEWS.md index d288cdbc..62981779 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,10 @@ # BioStructures.jl release notes +## v4.8.0 - Jul 2026 + +* `atomtypes` entries gain the ff14SB Lennard-Jones parameters `sigma` (Å) and `epsilon` (kcal/mol). +* `ff14SB_scale14` holds the force field's 1-4 nonbonded scaling factors. + ## v4.7.3 - Jul 2026 * Edge cases for `bondangle` and `dihedralangle` are handled correctly. diff --git a/Project.toml b/Project.toml index d82de7a0..b926ee91 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "BioStructures" uuid = "de9282ab-8554-53be-b2d6-f6c222edabfc" -version = "4.7.3" +version = "4.8.0" authors = ["Joe G Greener "] [deps] diff --git a/docs/src/api.md b/docs/src/api.md index 3d308d0c..d9f3886c 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -21,6 +21,8 @@ Order = [:module, :type, :constant, :function, :macro] ## Non-exported names +- [`BioStructures.atomtypes`](@ref) +- [`BioStructures.ff14SB_scale14`](@ref) - [`BioStructures.x`](@ref) - [`BioStructures.x!`](@ref) - [`BioStructures.y`](@ref) @@ -36,6 +38,8 @@ Private = false Order = [:module, :type, :constant, :function, :macro] ``` ```@docs +atomtypes +ff14SB_scale14 x x! y diff --git a/extractdata/bonding.jl b/extractdata/bonding.jl index e9e8459d..4e2e1551 100644 --- a/extractdata/bonding.jl +++ b/extractdata/bonding.jl @@ -41,13 +41,16 @@ function parsexmlline(f, line, tag, keynames...; skip=()) return key => (; vals...) end -atomtypes, residues, bondlengths, bondangles = open("protein.ff14SB.xml", "r") do io +rawatomtypes, residues, bondlengths, bondangles, nonbonded, scale14 = open("protein.ff14SB.xml", "r") do io line = readline(io) @assert line == "" 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}}() + # The XML is OpenMM's, so lengths are nm and energies kJ/mol; convert to Å and kcal/mol. + nonbonded = Dict{String, @NamedTuple{epsilon::Float64, sigma::Float64}}() + scale14 = nothing parsexmblock(io, "") do line if line == "" parsexmblock(io, "") do line @@ -119,22 +122,86 @@ atomtypes, residues, bondlengths, bondangles = open("protein.ff14SB.xml", "r") d end end) end + elseif startswith(line, "$", line) + m === nothing && error("Unrecognized NonbondedForce header $line") + scale14 = (coulomb = parse(Float32, m.captures[1]), lj = parse(Float32, m.captures[2])) + parsexmblock(io, "") do line + # declares that charges come from + # the residue templates, which is where residuedata already takes them from. + startswith(line, " Å + elseif k == "epsilon" + return parse(Float64, v[2:end-1]) / 4.184 # kJ/mol -> kcal/mol + else + error("Unknown nonbonded Atom key $k") + end + end) + end end end - atomtypes, residues, Dict(k => 10 * v.length for (k, v) in harmonicbonds), Dict(k => v.angle for (k, v) in harmonicangles) + atomtypes, residues, Dict(k => 10 * v.length for (k, v) in harmonicbonds), Dict(k => v.angle for (k, v) in harmonicangles), nonbonded, scale14 end +scale14 === nothing && error("no NonbondedForce block found") +let extra = setdiff(keys(nonbonded), (v.name for v in values(rawatomtypes))) + isempty(extra) || error("nonbonded parameters given for unknown types: $(sort!(collect(extra)))") +end +atomtypes = Dict{String, @NamedTuple{element::String, mass::Float32, name::String, sigma::Float32, epsilon::Float32}}( + class => begin + lj = get(() -> error("no nonbonded parameters for atom type $class ($(v.name))"), nonbonded, v.name) + (; v.element, v.mass, v.name, sigma = Float32(lj.sigma), epsilon = Float32(lj.epsilon)) + end for (class, v) in rawatomtypes) + +const atomtypesdoc = raw""" + atomtypes + +The atom types of the Amber ff14SB force field, keyed by type class (e.g., `atomtypes["CT"]`). +Each value is a `NamedTuple` with fields: + +- `element`: the chemical element symbol. +- `mass`: the atomic mass in amu. +- `name`: the force field's full name for the type, e.g., `"protein-CT"`. This is the value + stored in the `type` field of `residuedata`'s per-atom entries. +- `sigma`: the Lennard-Jones distance parameter in Å, the separation at which the pair + potential crosses zero. Amber's published parameter tables instead list `Rmin/2`, where + `Rmin = 2^(1/6) * sigma` is the separation of minimum energy. +- `epsilon`: the Lennard-Jones well depth in kcal/mol. Type `"HO"` has `epsilon = 0` and a + placeholder `sigma`: it carries no Lennard-Jones interaction. + +Partial charges are not properties of the type, as they vary between residues; they are in +`residuedata`. + +See also `ff14SB_scale14` for the scaling of nonbonded interactions between close neighbors. +""" + +const scale14doc = raw""" + ff14SB_scale14 + +The factors by which the ff14SB force field scales nonbonded interactions between atoms +separated by exactly three bonds (1-4 pairs): `coulomb` for electrostatics and `lj` for the +Lennard-Jones terms of `atomtypes`. +""" + +printdoc(io, doc) = (println(io, "\"\"\""); print(io, doc); println(io, "\"\"\"")) + open(joinpath(dirname(@__DIR__), "src", "bonding.jl"), "w") do io println(io, "# This file is auto-generated by extractdata/bonding.jl; do not edit directly.") println(io, "# It defines only data tables (atomtypes, residuedata, bondlengths, bondangles);") println(io, "# code that uses them belongs in a hand-written file such as src/atombonds.jl.\n") - println(io, "const atomtypes = Dict{String, @NamedTuple{element::String, mass::Float32, name::String}}(") + printdoc(io, atomtypesdoc) + println(io, "const atomtypes = Dict{String, @NamedTuple{element::String, mass::Float32, name::String, sigma::Float32, epsilon::Float32}}(") at = sort!(collect(atomtypes); by=first) for pr in at println(io, " ", pr, ',') end println(io, ")\n") + printdoc(io, scale14doc) + println(io, "const ff14SB_scale14 = ", scale14, "\n") + println(io, "const RDADict = Dict{String, @NamedTuple{charge::Float32, type::String}}") println(io, "const residuedata = Dict{String, @NamedTuple{atoms::RDADict, bonds::Vector{Tuple{String,String}}, externalbonds::Vector{String}}}(") diff --git a/src/bonding.jl b/src/bonding.jl index e118c49e..614f18da 100644 --- a/src/bonding.jl +++ b/src/bonding.jl @@ -2,43 +2,73 @@ # It defines only data tables (atomtypes, residuedata, bondlengths, bondangles); # code that uses them belongs in a hand-written file such as src/atombonds.jl. -const atomtypes = Dict{String, @NamedTuple{element::String, mass::Float32, name::String}}( - "2C" => (element = "C", mass = 12.01f0, name = "protein-2C"), - "3C" => (element = "C", mass = 12.01f0, name = "protein-3C"), - "C" => (element = "C", mass = 12.01f0, name = "protein-C"), - "C*" => (element = "C", mass = 12.01f0, name = "protein-C*"), - "C8" => (element = "C", mass = 12.01f0, name = "protein-C8"), - "CA" => (element = "C", mass = 12.01f0, name = "protein-CA"), - "CB" => (element = "C", mass = 12.01f0, name = "protein-CB"), - "CC" => (element = "C", mass = 12.01f0, name = "protein-CC"), - "CN" => (element = "C", mass = 12.01f0, name = "protein-CN"), - "CO" => (element = "C", mass = 12.01f0, name = "protein-CO"), - "CR" => (element = "C", mass = 12.01f0, name = "protein-CR"), - "CT" => (element = "C", mass = 12.01f0, name = "protein-CT"), - "CV" => (element = "C", mass = 12.01f0, name = "protein-CV"), - "CW" => (element = "C", mass = 12.01f0, name = "protein-CW"), - "CX" => (element = "C", mass = 12.01f0, name = "protein-CX"), - "H" => (element = "H", mass = 1.008f0, name = "protein-H"), - "H1" => (element = "H", mass = 1.008f0, name = "protein-H1"), - "H4" => (element = "H", mass = 1.008f0, name = "protein-H4"), - "H5" => (element = "H", mass = 1.008f0, name = "protein-H5"), - "HA" => (element = "H", mass = 1.008f0, name = "protein-HA"), - "HC" => (element = "H", mass = 1.008f0, name = "protein-HC"), - "HO" => (element = "H", mass = 1.008f0, name = "protein-HO"), - "HP" => (element = "H", mass = 1.008f0, name = "protein-HP"), - "HS" => (element = "H", mass = 1.008f0, name = "protein-HS"), - "N" => (element = "N", mass = 14.01f0, name = "protein-N"), - "N2" => (element = "N", mass = 14.01f0, name = "protein-N2"), - "N3" => (element = "N", mass = 14.01f0, name = "protein-N3"), - "NA" => (element = "N", mass = 14.01f0, name = "protein-NA"), - "NB" => (element = "N", mass = 14.01f0, name = "protein-NB"), - "O" => (element = "O", mass = 16.0f0, name = "protein-O"), - "O2" => (element = "O", mass = 16.0f0, name = "protein-O2"), - "OH" => (element = "O", mass = 16.0f0, name = "protein-OH"), - "S" => (element = "S", mass = 32.06f0, name = "protein-S"), - "SH" => (element = "S", mass = 32.06f0, name = "protein-SH"), +""" + atomtypes + +The atom types of the Amber ff14SB force field, keyed by type class (e.g., `atomtypes["CT"]`). +Each value is a `NamedTuple` with fields: + +- `element`: the chemical element symbol. +- `mass`: the atomic mass in amu. +- `name`: the force field's full name for the type, e.g., `"protein-CT"`. This is the value + stored in the `type` field of `residuedata`'s per-atom entries. +- `sigma`: the Lennard-Jones distance parameter in Å, the separation at which the pair + potential crosses zero. Amber's published parameter tables instead list `Rmin/2`, where + `Rmin = 2^(1/6) * sigma` is the separation of minimum energy. +- `epsilon`: the Lennard-Jones well depth in kcal/mol. Type `"HO"` has `epsilon = 0` and a + placeholder `sigma`: it carries no Lennard-Jones interaction. + +Partial charges are not properties of the type, as they vary between residues; they are in +`residuedata`. + +See also `ff14SB_scale14` for the scaling of nonbonded interactions between close neighbors. +""" +const atomtypes = Dict{String, @NamedTuple{element::String, mass::Float32, name::String, sigma::Float32, epsilon::Float32}}( + "2C" => (element = "C", mass = 12.01f0, name = "protein-2C", sigma = 3.3996694f0, epsilon = 0.1094f0), + "3C" => (element = "C", mass = 12.01f0, name = "protein-3C", sigma = 3.3996694f0, epsilon = 0.1094f0), + "C" => (element = "C", mass = 12.01f0, name = "protein-C", sigma = 3.3996694f0, epsilon = 0.086f0), + "C*" => (element = "C", mass = 12.01f0, name = "protein-C*", sigma = 3.3996694f0, epsilon = 0.086f0), + "C8" => (element = "C", mass = 12.01f0, name = "protein-C8", sigma = 3.3996694f0, epsilon = 0.1094f0), + "CA" => (element = "C", mass = 12.01f0, name = "protein-CA", sigma = 3.3996694f0, epsilon = 0.086f0), + "CB" => (element = "C", mass = 12.01f0, name = "protein-CB", sigma = 3.3996694f0, epsilon = 0.086f0), + "CC" => (element = "C", mass = 12.01f0, name = "protein-CC", sigma = 3.3996694f0, epsilon = 0.086f0), + "CN" => (element = "C", mass = 12.01f0, name = "protein-CN", sigma = 3.3996694f0, epsilon = 0.086f0), + "CO" => (element = "C", mass = 12.01f0, name = "protein-CO", sigma = 3.3996694f0, epsilon = 0.086f0), + "CR" => (element = "C", mass = 12.01f0, name = "protein-CR", sigma = 3.3996694f0, epsilon = 0.086f0), + "CT" => (element = "C", mass = 12.01f0, name = "protein-CT", sigma = 3.3996694f0, epsilon = 0.1094f0), + "CV" => (element = "C", mass = 12.01f0, name = "protein-CV", sigma = 3.3996694f0, epsilon = 0.086f0), + "CW" => (element = "C", mass = 12.01f0, name = "protein-CW", sigma = 3.3996694f0, epsilon = 0.086f0), + "CX" => (element = "C", mass = 12.01f0, name = "protein-CX", sigma = 3.3996694f0, epsilon = 0.1094f0), + "H" => (element = "H", mass = 1.008f0, name = "protein-H", sigma = 1.0690784f0, epsilon = 0.0157f0), + "H1" => (element = "H", mass = 1.008f0, name = "protein-H1", sigma = 2.471353f0, epsilon = 0.0157f0), + "H4" => (element = "H", mass = 1.008f0, name = "protein-H4", sigma = 2.5105526f0, epsilon = 0.015f0), + "H5" => (element = "H", mass = 1.008f0, name = "protein-H5", sigma = 2.4214628f0, epsilon = 0.015f0), + "HA" => (element = "H", mass = 1.008f0, name = "protein-HA", sigma = 2.5996425f0, epsilon = 0.015f0), + "HC" => (element = "H", mass = 1.008f0, name = "protein-HC", sigma = 2.6495328f0, epsilon = 0.0157f0), + "HO" => (element = "H", mass = 1.008f0, name = "protein-HO", sigma = 10.0f0, epsilon = 0.0f0), + "HP" => (element = "H", mass = 1.008f0, name = "protein-HP", sigma = 1.9599771f0, epsilon = 0.0157f0), + "HS" => (element = "H", mass = 1.008f0, name = "protein-HS", sigma = 1.0690784f0, epsilon = 0.0157f0), + "N" => (element = "N", mass = 14.01f0, name = "protein-N", sigma = 3.2499986f0, epsilon = 0.17f0), + "N2" => (element = "N", mass = 14.01f0, name = "protein-N2", sigma = 3.2499986f0, epsilon = 0.17f0), + "N3" => (element = "N", mass = 14.01f0, name = "protein-N3", sigma = 3.2499986f0, epsilon = 0.17f0), + "NA" => (element = "N", mass = 14.01f0, name = "protein-NA", sigma = 3.2499986f0, epsilon = 0.17f0), + "NB" => (element = "N", mass = 14.01f0, name = "protein-NB", sigma = 3.2499986f0, epsilon = 0.17f0), + "O" => (element = "O", mass = 16.0f0, name = "protein-O", sigma = 2.9599218f0, epsilon = 0.21f0), + "O2" => (element = "O", mass = 16.0f0, name = "protein-O2", sigma = 2.9599218f0, epsilon = 0.21f0), + "OH" => (element = "O", mass = 16.0f0, name = "protein-OH", sigma = 3.0664735f0, epsilon = 0.2104f0), + "S" => (element = "S", mass = 32.06f0, name = "protein-S", sigma = 3.5635948f0, epsilon = 0.25f0), + "SH" => (element = "S", mass = 32.06f0, name = "protein-SH", sigma = 3.5635948f0, epsilon = 0.25f0), ) +""" + ff14SB_scale14 + +The factors by which the ff14SB force field scales nonbonded interactions between atoms +separated by exactly three bonds (1-4 pairs): `coulomb` for electrostatics and `lj` for the +Lennard-Jones terms of `atomtypes`. +""" +const ff14SB_scale14 = (coulomb = 0.8333333f0, lj = 0.5f0) + const RDADict = Dict{String, @NamedTuple{charge::Float32, type::String}} const residuedata = Dict{String, @NamedTuple{atoms::RDADict, bonds::Vector{Tuple{String,String}}, externalbonds::Vector{String}}}( "ACE" => (atoms = RDADict("CH3" => (charge = -0.3662, type = "protein-CT"), "HH31" => (charge = 0.1123, type = "protein-HC"), "C" => (charge = 0.5972, type = "protein-C"), "HH32" => (charge = 0.1123, type = "protein-HC"), "HH33" => (charge = 0.1123, type = "protein-HC"), "O" => (charge = -0.5679, type = "protein-O")), diff --git a/test/runtests.jl b/test/runtests.jl index cf0b75ba..af266d80 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3567,6 +3567,25 @@ end # Bond lengths and angles @test BioStructures.bondlengths[("protein-2C", "protein-2C")] == 1.526f0 @test BioStructures.bondangles[("protein-CT", "protein-C", "protein-N")] == 2.035054f0 + + # Lennard-Jones parameters + at = BioStructures.atomtypes + @test length(at) == 34 + # Every type referenced by a residue template has parameters + typenames = Set(t.name for t in values(at)) + @test all(a.type ∈ typenames for r in values(rd) for a in values(r.atoms)) + # Amber tabulates Rmin/2 = 1.9080 Å and ε = 0.0860 kcal/mol for type C + @test at["C"].epsilon ≈ 0.0860f0 atol=1e-5 + @test at["C"].sigma ≈ 3.39967f0 atol=1e-4 + @test 2^(1/6) * at["C"].sigma / 2 ≈ 1.9080f0 atol=1e-3 + for (class, t) in at + @test isfinite(t.sigma) && t.sigma > 0 + @test isfinite(t.epsilon) && t.epsilon ≥ 0 + # HO is a non-interacting dummy: zero well depth and a placeholder sigma + @test (t.epsilon > 0) == (class != "HO") + end + @test BioStructures.ff14SB_scale14.coulomb ≈ 5/6 + @test BioStructures.ff14SB_scale14.lj == 0.5f0 end @testset "Secondary structure" begin From caec19fe331608fc3cd59fd26686373460b520b3 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Wed, 29 Jul 2026 07:29:38 -0500 Subject: [PATCH 2/2] Fix a path inconsistency --- extractdata/bonding.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/extractdata/bonding.jl b/extractdata/bonding.jl index 4e2e1551..d9bf7150 100644 --- a/extractdata/bonding.jl +++ b/extractdata/bonding.jl @@ -3,8 +3,10 @@ using Downloads -if !isfile(joinpath(@__DIR__, "protein.ff14SB.xml")) - Downloads.download("https://raw.githubusercontent.com/openmm/openmm/refs/heads/master/wrappers/python/openmm/app/data/amber14/protein.ff14SB.xml", "protein.ff14SB.xml") +const xmlfile = joinpath(@__DIR__, "protein.ff14SB.xml") + +if !isfile(xmlfile) + Downloads.download("https://raw.githubusercontent.com/openmm/openmm/refs/heads/master/wrappers/python/openmm/app/data/amber14/protein.ff14SB.xml", xmlfile) end function parsexmblock(f, io::IO, key) @@ -41,7 +43,7 @@ function parsexmlline(f, line, tag, keynames...; skip=()) return key => (; vals...) end -rawatomtypes, residues, bondlengths, bondangles, nonbonded, scale14 = open("protein.ff14SB.xml", "r") do io +rawatomtypes, residues, bondlengths, bondangles, nonbonded, scale14 = open(xmlfile, "r") do io line = readline(io) @assert line == "" atomtypes = Dict{String, @NamedTuple{element::String, mass::Float32, name::String}}()