Structures and unit cells

read_pdb, read_mmcif, read_xyz and wget return a Structure: a vector of atoms that also carries the metadata of the system: its unit cell, title, and connectivity. A Structure behaves as a vector of Atoms, such that all functions that operate on vectors of atoms also operate on structures.

Structure(atoms::AbstractVector{<:Atom}=Atom{Nothing}[]; unitcell=nothing)
Structure(atoms::AbstractVector{<:Atom}, metadata::MetaData)

julia> using PDBTools

julia> s = Structure(read_pdb(PDBTools.test_file("small.pdb")); unitcell=UnitCell(10.0, 10.0, 10.0));

julia> s.unitcell
UnitCell(a=10.000, b=10.000, c=10.000, α=90.00, β=90.00, γ=90.00)

julia> s.unitcell = UnitCell(a=10.0, b=20.0, c=30.0, γ=120.0)
UnitCell(a=10.000, b=20.000, c=30.000, α=90.00, β=90.00, γ=120.00)

julia> length(s[1:10])
10

julia> s[1:10].unitcell == s.unitcell
true

MetaData

Metadata

The metadata fields are accessible as properties of the structure, and are nothing when not defined:

The metadata can be set in the Structure constructor (e.g. Structure(atoms; title="My system")), or by assigning the properties (e.g. str.title = "My system").

Reading structures

The unit cell is read from the CRYST1 field of PDB files, or from the _cell fields of mmCIF files:

using PDBTools, LazyArtifacts
str = read_pdb(PDBTools.test_file("pbc.pdb"))

   Structure{Atom{Nothing}} with 130241 atoms - with unit cell
   index name resname chain   resnum  residue        x        y        z occup  beta model segname index_pdb
       1    N     ALA     P        1        1  103.848   97.435   20.875  1.00  0.00     1    PROT         1
       2  HT1     ALA     P        1        1  103.498   96.484   21.110  1.00  0.00     1    PROT         2
⋮
  130240   O3    GLYC     D     4302    26788   26.526   49.111   13.033  1.00  0.00     1    GLYC    130240
  130241   H3    GLYC     D     4302    26788   27.053   48.661   12.369  1.00  0.00     1    GLYC    130241

str.unitcell

UnitCell(a=107.845, b=107.845, c=107.845, α=90.00, β=90.00, γ=90.00)

If the file does not define a unit cell (or defines the 1Å unit cube, which indicates that the structure is not crystallographic, as in NMR or cryo-EM structures), str.unitcell is nothing.

Title

The title is read from the TITLE records of PDB files, the _struct.title field of mmCIF files, or the comment line of xyz files:

ethanol = read_pdb(PDBTools.test_file("ethanol.pdb"))
ethanol.title

"ETHANOL MOLECULE WITH EXPLICIT HYDROGENS"

Connectivity

The bonds of the atoms are read from the CONECT records of PDB files, or from the _geom_bond fields of mmCIF files. The connectivity is a Vector{Vector{Int32}}, where connectivity[i] contains the indices of the atoms bonded to the atom with index i:

ethanol.connectivity

9-element Vector{Vector{Int32}}:
 [2, 4, 5, 6]
 [1, 3, 7, 8]
 [2, 9]
 ⋮
 [2]
 [3]

The indices are the sequential indices of the atoms in the file, index(atom), which are not changed by selections, filtering, or reordering of the atoms. Thus, the connectivity remains valid for substructures. For example, the atoms bonded to the oxygen atom of ethanol are:

heavy = select(ethanol, "not element H")
O = heavy[findfirst(at -> name(at) == "O1", heavy)]
[ name(ethanol[j]) for j in ethanol.connectivity[index(O)] ]

2-element Vector{InlineStrings.String15}:
 "C2"
 "H6"

When writing, write_pdb and write_mmcif number the atoms sequentially in the file (1, 2, ...) and convert the connectivity to these numbers. Bonds to atoms that are not written are dropped:

tmpfile = tempname() * ".pdb"
write_pdb(tmpfile, heavy)
read_pdb(tmpfile).connectivity

3-element Vector{Vector{Int32}}:
 [2]
 [1, 3]
 [2]

Indexing and iteration

A Structure can be indexed and iterated as a vector of atoms:

str[1]

   index name resname chain   resnum  residue        x        y        z occup  beta model segname index_pdb
       1    N     ALA     P        1        1  103.848   97.435   20.875  1.00  0.00     1    PROT         1

collect(eachchain(str))

8-element Vector{Chain}[ 
    Chain(P-3801 atoms)
    Chain(A-29997 atoms)
    ⋮
    Chain(E-106 atoms)
    Chain(D-60228 atoms)
]

Indexing with ranges, selections, filter, copy, and other functions that create new vectors of atoms return new Structures, with a copy of the metadata. Thus, the unit cell is preserved:

protein = select(str, "protein")
protein.unitcell

UnitCell(a=107.845, b=107.845, c=107.845, α=90.00, β=90.00, γ=90.00)

The underlying vector of atoms is available as str.atoms.

Unit cells

The unit cell is set with a UnitCell, which can be constructed from the lattice parameters, from the sides of an orthorhombic box, or from a matrix with the lattice vectors as columns:

str.unitcell = UnitCell(85.0, 85.0, 85.0)

UnitCell(a=85.000, b=85.000, c=85.000, α=90.00, β=90.00, γ=90.00)

str.unitcell = UnitCell(a=85.0, b=85.0, c=85.0, γ=120.0)

UnitCell(a=85.000, b=85.000, c=85.000, α=90.00, β=90.00, γ=120.00)

UnitCell(matrix::AbstractMatrix)
UnitCell(sides::AbstractVector)
UnitCell(a, b, c, α=90, β=90, γ=90)
UnitCell(; a, b, c, α=90, β=90, γ=90)
UnitCell(lattice::NamedTuple)

julia> using PDBTools

julia> uc = UnitCell(10.0, 20.0, 30.0)
UnitCell(a=10.000, b=20.000, c=30.000, α=90.00, β=90.00, γ=90.00)

julia> UnitCell([10.0, 20.0, 30.0]) == uc
true

julia> UnitCell(a=10.0, b=20.0, c=30.0, γ=120.0)
UnitCell(a=10.000, b=20.000, c=30.000, α=90.00, β=90.00, γ=120.00)

julia> matrix_to_lattice(uc)
(a = 10.0f0, b = 20.0f0, c = 30.0f0, α = 90.0f0, β = 90.0f0, γ = 90.0f0)

Unit cells are written by write_pdb and write_mmcif.

Periodic boundary conditions

The functions that support periodic boundary conditions (for example, contact_map, hydrogen_bonds, sasa_particles, and visualize) accept the unitcell keyword argument as:

• unitcell=false (default): no periodic boundary conditions.
• unitcell=true: use the unit cell of the Structure.
• unitcell=UnitCell(...): use the given unit cell.

str = read_pdb(PDBTools.test_file("pbc.pdb"), "protein")
sasa(sasa_particles(str; unitcell=true))

10838.204f0