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

sim = Simulation(Testing.namd_pdb, Testing.namd_traj; 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
using PDBTools
using Plots: plot

sim = Simulation(Testing.namd_pdb, Testing.namd_traj)

# 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.3506959154055
 32.32651999592368
 33.38388153244739
 34.739991307891934
 30.118300797521297

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])
Simulation(pdb_file::String, trajectory_file::String; frames=9:2:20)
Simulation(pdb_file::String, trajectory_file::String; first=1, last=nothing, step=1)
Simulation(atoms::AbstractVector{<:AtomType}, trajectory_file::String; first=1, last=nothing, step=1)

Creates a new Simulation object.

The first constructor creates a Simulation object from a PDB or mmCIF file and a trajectory file. It will use the PDBTools.Atom for the atom type, which will populate the atoms vector of the Simulation object. Currently, other atom types are supported, if the MolSimToolkit.atomic_mass(::AtomType) function is defined for the atom type.

With the second constructor, the atoms vector is passed as an argument. This is useful when the atoms are provided by a different source than the PDB file.

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.

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; frame=nothing): get vector of atoms of the simulation, frame=nothing, the positions of the current frame will be returned.

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

julia> simulation = Simulation(
           Testing.namd_pdb, Testing.namd_traj; 
           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)