Structural alignment

align(x, y; mass = nothing)
align!(x, y; mass = nothing)

Aligns two structures (sets of points in 3D space). Solves the "Procrustes" problem, which is to find the best translation, and rotation, that aligns the two structures, minimizing the RMSD between them.

Structures are expected to be of the same size, and the correspondence is assumed from the vector indices.

align returns a new vector containing the coordinates of x aligned to y. align! modifies the input vector x in place.

align(x, y; mass = nothing)
align!(x, y; mass = nothing)

Aligns two structures (sets of points in 3D space). Solves the "Procrustes" problem, which is to find the best translation, and rotation, that aligns the two structures, minimizing the RMSD between them.

Structures are expected to be of the same size, and the correspondence is assumed from the vector indices.

align returns a new vector containing the coordinates of x aligned to y. align! modifies the input vector x in place.

rmsd(x::AbstractVector,y::AbstractVector)
rmsd(simulation::Simulation, indices::AbstractVector{Int}; mass = nothing, reference_frame = nothing, show_progress = true)

Computes the root mean square deviation (RMSD) between two sets of points in space, or along a trajectory.

The rmsd(x,y) method computes the RMSD between two sets of points x and y. The sets are expected to be of the same size, and the correspondence is assumed from the vector indices.

The rmsd(simulation::Simulation, indices) method computes the RMSD along a trajectory. With the following options:

• The indices vector contains the indices of the atoms to be considered.
• The mass argument can be used to provide the mass of the atoms if they are not the same.
• The reference_frame argument can be used to provide a reference frame to align the trajectory to:
  • If reference_frame == nothing, the first frame will be used (default behavior).
  • If reference_frame == :average, the average structure will be used.
  • If reference_frame is an integer, the frame at that index will be used as reference.

Examples

If the set is compared toi tself, the RMSD should be zero:

julia> using MolSimToolkit, MolSimToolkit.Testing

julia> using PDBTools

julia> ca = coor(readPDB(Testing.namd_pdb), "name CA");

julia> rmsd(ca, ca)
0.0

Computing the rmsd along a trajectory

julia> using MolSimToolkit, MolSimToolkit.Testing

julia> using PDBTools

julia> atoms = readPDB(Testing.namd_pdb);

julia> simulation = Simulation(Testing.namd_pdb, Testing.namd_traj);

julia> cas = findall(sel"name CA", atoms); # CA indices

julia> rmsd(simulation, cas; show_progress=false)
5-element Vector{Float64}:
 0.0
 2.8388710154609034
 2.9776998440690385
 2.4621444212469483
 3.8035683196100796

julia> rmsd(simulation, cas; reference_frame=:average, show_progress=false)
5-element Vector{Float64}:
 1.8995986972454748
 2.1512244220536973
 1.5081703191869376
 1.1651111324544219
 2.757039151265317

rmsd_matrix(
    simulation::Simulation, 
    indices::AbstractVector{Int}; 
    mass::Union{AbstractVector{Int}, Nothing} = nothing,
    align::Bool = true,
    show_progress = true,
)

Computes the RMSD matrix for a set of atoms along a trajectory.

The indices vector contains the indices of the atoms to be considered. The mass argument can be used to provide the mass of the atoms if they are not the same. The align argument can be used to align the frames before computing the RMSD.

The show_progress argument can be used to show a progress bar.

Example

julia> using MolSimToolkit, MolSimToolkit.Testing

julia> using PDBTools

julia> atoms = readPDB(Testing.namd_pdb);

julia> cas = findall(Select("name CA"), atoms); # CA indices

julia> simulation = Simulation(Testing.namd_pdb, Testing.namd_traj);

julia> rmsd_matrix(simulation, cas; show_progress=false)
5×5 Matrix{Float64}:
 0.0      2.83887  2.9777   2.46214  3.80357
 2.83887  0.0      2.35492  2.64463  4.68028
 2.9777   2.35492  0.0      2.08246  3.46149
 2.46214  2.64463  2.08246  0.0      2.97835
 3.80357  4.68028  3.46149  2.97835  0.0