Add Lennard-Jones parameters to atomtypes

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.

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 <jgreener@hotmail.co.uk>"]
 
 [deps]

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

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,13 +43,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(xmlfile, "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}}()
+    # 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, "</ForceField>") do line
         if line == "<AtomTypes>"
             parsexmblock(io, "</AtomTypes>") do line
@@ -119,22 +124,86 @@ atomtypes, residues, bondlengths, bondangles = open("protein.ff14SB.xml", "r") d
                     end
                 end)
             end
+        elseif startswith(line, "<NonbondedForce ")
+            m = match(r"^<NonbondedForce coulomb14scale=\"([^\"]*)\" lj14scale=\"([^\"]*)\">$", line)
+            m === nothing && error("Unrecognized NonbondedForce header $line")
+            scale14 = (coulomb = parse(Float32, m.captures[1]), lj = parse(Float32, m.captures[2]))
+            parsexmblock(io, "</NonbondedForce>") do line
+                # <UseAttributeFromResidue name="charge"/> declares that charges come from
+                # the residue templates, which is where residuedata already takes them from.
+                startswith(line, "<UseAttributeFromResidue ") && return nothing
+                push!(nonbonded, parsexmlline(line, "Atom", "type") do k, v
+                    if k == "sigma"
+                        return 10 * parse(Float64, v[2:end-1])  # nm -> Å
+                    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}}}(")

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")),

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