Simulation

The Simulation object is the central data structure of MolSimToolkit. It holds the atom data, the trajectory file, and the current frame. Most analysis functions in MolSimToolkit accept a Simulation as their primary argument.

Iterating over frames

The standard way to process a trajectory is to iterate over all frames with a for loop:

using MolSimToolkit, MolSimToolkit.Testing, LazyArtifacts

sim = Simulation(test_file("popc_membrane.pdb"), test_file("popc_membrane.dcd"); first=2, step=2, last=4)

for frame in sim
    @show frame_index(sim)
    @show positions(frame)[1]
end

frame_index(sim) = 2
(positions(frame))[1] = [2.001203775405884, 7.035768985748291, 43.27251434326172]
frame_index(sim) = 4
(positions(frame))[1] = [-0.39684221148490906, 12.047886848449707, 37.44126892089844]

The pairs iterator yields (frame_index, frame) tuples — useful when you need the actual frame number in the trajectory:

for (i, frame) in pairs(sim)
    @show i, frame_index(sim)
end

(i, frame_index(sim)) = (2, 2)
(i, frame_index(sim)) = (4, 4)

enumerate yields (counter, frame) tuples — useful when you just need a sequential count:

for (i, frame) in enumerate(sim)
    @show i, frame_index(sim)
end

(i, frame_index(sim)) = (1, 2)
(i, frame_index(sim)) = (2, 4)

Example: end-to-end distance

As a simple illustration, we compute the end-to-end distance of the protein (distance between the C$_\alpha$ atoms of the first and last residues) across all frames, applying periodic boundary conditions when wrapping the second atom position.

using MolSimToolkit, MolSimToolkit.Testing, LazyArtifacts
using PDBTools
using Plots: plot

sim = Simulation(test_file("popc_membrane.pdb"), test_file("popc_membrane.dcd"))

# Indices of the Cα atoms of the first and last residues
ats = get_atoms(sim)
protein_ca = findall(sel"protein and name CA", ats)
i1, i2 = first(protein_ca), last(protein_ca)

# Iterate over frames and compute end-to-end distances
end_to_end = Float64[]
for frame in sim
    p = positions(frame)
    p1 = p[i1]
    p2 = wrap(p[i2], p1, unitcell(frame))  # apply PBC
    push!(end_to_end, distance(p2,  p1))
end

end_to_end

5-element Vector{Float64}:
 42.350696430918255
 32.32651865283743
 33.38388339753901
 34.739990657172726
 30.118298361705577

We can plot the result as a function of the frame number:

plot(MolSimStyle,
    frame_range(sim), end_to_end;
    xlabel="frame",
    ylabel="end-to-end distance / Å",
    label=nothing,
)

Creating a Simulation

MolSimToolkit.Simulation — Type

Simulation(pdb_file::String, trajectory_file::String; frames=[1,2,3,5], unitcell=nothing)
Simulation(pdb_file::String, trajectory_file::String; frames=9:2:20, unitcell=nothing)
Simulation(pdb_file::String, trajectory_file::String; first=1, last=nothing, step=1, unitcell=nothing)
Simulation(pdb_file::String, trajectory_file1::String, trajectory_file2::String, ...; kargs...)
Simulation(pdb_file::String, trajectory_files::AbstractVector{<:AbstractString}; kargs...)
Simulation(atoms::AbstractVector{<:AtomType}, trajectory_file::String; kargs...)
Simulation(atoms::AbstractVector{<:AtomType}, trajectory_file1::String, trajectory_file2::String, ...; kargs...)
Simulation(atoms::AbstractVector{<:AtomType}, trajectory_files::AbstractVector{<:AbstractString}; kargs...)

Creates a new Simulation object.

The first argument defines the atoms of the simulation, and can be:

pdb_file: the name of a PDB or mmCIF file. The atoms are read with PDBTools.read_pdb, and will be of type PDBTools.Atom.
atoms: a vector of atoms. This is useful when the atoms are provided by a different source than a PDB file, or were modified after being read. The atom type is not restricted to PDBTools.Atom: other atom types are supported, if the MolSimToolkit.atomic_mass(::AtomType) function is defined for the atom type. In this case path_pdb(simulation) returns nothing.

The following arguments are the trajectory files. More than one trajectory file can be provided, either as additional arguments or as a vector of file names. The files are considered sequential parts of the same simulation: frames are numbered continuously across files, in the order the files were given, and all analyses can be performed on the concatenated trajectory. All files must contain the same number of atoms.

The keyword arguments kargs... are the same for all constructors.

The frames or first, last, and step arguments can be used to specify the frames to be iterated over:

frames can be a vector of frame indices, e. g., frames=[1,2,3,5] or frames=9:2:20.
first, last, and step are Integers that specify the frames to be iterated over. If last is not specified, the last frame in the trajectory will be used.

By default, the unit cell of each frame is read from the trajectory file. The unitcell keyword sets a fixed unit cell, which will be used for all frames, ignoring the unit cell of the trajectory (useful for trajectories without unit cell information, or of constant-volume simulations):

unitcell=true: use the unit cell of the PDB file (or of the atoms, if they are a PDBTools.Structure with a unit cell).
unitcell=uc: use the unit cell uc, a UnitCell object, which can be created from the sides of an orthorhombic box, the matrix of lattice vectors, or the lattice parameters (see the UnitCell documentation of PDBTools). Other values are converted with UnitCell(uc), e. g. unitcell=[10.0, 10.0, 10.0].

A Simulation object contains a trajectory file and a PDB data of the atoms. It can be iterated over to obtain the frames in the trajectory. The Simulation object is a mutable struct that contains the following data, that can be retrieved by the corresponding functions:

frame_range(::Simulation): the list of frames to be iterated over
frame_index(::Simulation): the index of the current frame in the trajectory
length(::Simulation): the number of frames to be iterated over in the trajectory file, considering the current range
raw_length(::Simulation): the number of frames in the trajectory file
get_atoms(::Simulation): get the atoms of the simulation, with the positions and unit cell of the current frame, as a PDBTools.Structure.
unitcell(::Frame): the unit cell of a frame (unitcell(current_frame(simulation)) for the current frame).

The Simulation object can also be manipulated by the following functions:

close(::Simulation): closes the trajectory file
restart!(::Simulation): restarts the iteration over the trajectory file, current frame is the first frame.
first_frame!(::Simulation): restarts the iteration over the trajectory file and places the current frame at the first frame in the trajectory
current_frame(::Simulation): returns the current frame in the trajectory
next_frame!(::Simulation): reads the next frame in the trajectory file and returns it. Moves the current frame to the next one.
set_frame_range!(::Simulation; first, last, step): resets the range of frames to be iterated over.
get_frame!(::Simulation, iframe): returns the frame at the given index in the trajectory.
goto_frame!(::Simulation, iframe): similar to get_frame!, but returns the Simulation object placed in the desired frame.

One important feature of the Simulation object is that it can be iterated over, frame by frame.

The pairs iterator can also be used to iterate over the frames, returning a tuple with the frame index and the frame itself.

The enumerate iterator can also be used to iterate over the frames, returning a tuple with the frame counter and the frame itself.

Examples

julia> using MolSimToolkit, MolSimToolkit.Testing, LazyArtifacts

julia> simulation = Simulation(
           test_file("popc_membrane.pdb"), test_file("popc_membrane.dcd"); 
           first = 2, step = 2, last = 4,
           # or frames = [2,4]
       );

julia> for frame in simulation 
           @show frame_index(simulation)
           # show x coordinate of first atom 
           @show positions(frame)[1].x
       end
frame_index(simulation) = 2
((positions(frame))[1]).x = 2.001203775405884
frame_index(simulation) = 4
((positions(frame))[1]).x = -0.39684221148490906 

julia> for (i, frame) in pairs(simulation)
           @show i, frame_index(simulation)
       end
(i, frame_index(simulation)) = (2, 2)
(i, frame_index(simulation)) = (4, 4)  

julia> for (i, frame) in enumerate(simulation)
           @show i, frame_index(simulation)
       end
(i, frame_index(simulation)) = (1, 2)
(i, frame_index(simulation)) = (2, 4)