Merge pull request #28 from m3g/molsimtoolkit_shared

Molsimtoolkit shared

.github/workflows/ci.yml

@@ -2,6 +2,9 @@ name: CI
 on:
   - push
   - pull_request
+concurrency:
+  group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }}
+  cancel-in-progress: true
 env:
   JULIA_NUM_THREADS: 2
 jobs:

Project.toml

@@ -11,6 +11,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 EasyFit = "fde71243-0cda-4261-b7c7-4845bd106b21"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MolSimToolkitShared = "aea21201-7910-416e-b888-414710bebccb"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 PDBTools = "e29189f1-7114-4dbd-93d0-c5673a921a58"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
@@ -39,6 +40,7 @@ Documenter = "1.3"
 EasyFit = "0.6.8"
 LaTeXStrings = "1.3"
 LinearAlgebra = "1.9"
+MolSimToolkitShared = "1.1.2"
 OffsetArrays = "1.13"
 PDBTools = "1.8.1, 2"
 Plots = "1.39"
@@ -65,4 +67,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
 [targets]
-test = ["Aqua", "Test", "TestItemRunner", "BenchmarkTools", "DelimitedFiles", "Rotations", "Documenter", "Plots"]
+test = ["Aqua", "Test", "TestItemRunner", "BenchmarkTools", "DelimitedFiles", "Rotations", "Documenter", "Plots"]

docs/Project.toml

@@ -3,6 +3,7 @@ Chemfiles = "46823bd8-5fb3-5f92-9aa0-96921f3dd015"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 MolSimToolkit = "054db54f-6694-444d-9bbb-e9ecdbfe77be"
+MolSimToolkitShared = "aea21201-7910-416e-b888-414710bebccb"
 PDBTools = "e29189f1-7114-4dbd-93d0-c5673a921a58"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

docs/make.jl

@@ -4,8 +4,9 @@ using Documenter
 push!(LOAD_PATH, "../")
 push!(LOAD_PATH, "../src/")
 using MolSimToolkit
+using MolSimToolkitShared
 makedocs(
-    modules=[MolSimToolkit],
+    modules=[MolSimToolkit, MolSimToolkitShared],
     sitename="MolSimToolkit.jl",
     doctest=false,
     pages=[
@@ -14,7 +15,7 @@ makedocs(
             "Distances and misc." => "Structural_properties.md",
             "Dihedral angle analysis" => "Dihedrals.md",
             "Secondary structure" => "secondary_structures.md",
-            "Structural alignment" => "procrustes.md",
+            "Structural alignment" => "structural_alignment.md",
         ],
         "Simulation statistics" => Any[ 
             " Block averages" => "block_averages.md",

docs/src/Developer.md

@@ -31,7 +31,7 @@ Filter = (f) -> (nameof(f) === :unitcell)
 ## Wrap coordinates
 
 ```@autodocs
-Modules = [ MolSimToolkit ]
+Modules = [ MolSimToolkitShared ]
 Pages = [ "wrap.jl" ]
 ```
 

docs/src/procrustes.md → docs/src/structural_alignment.md

@@ -8,5 +8,9 @@ CollapsedDocStrings = true
 
 ```@autodocs
 Modules = [ MolSimToolkit ]
-Pages = [ "procrustes.jl" ]
+Pages = [ "structural_alignment.jl" ]
+```
+```@autodocs
+Modules = [ MolSimToolkitShared ]
+Pages = [ "structural_alignment.jl" ]
 ```

src/MolSimToolkit.jl

@@ -5,6 +5,17 @@ import PDBTools
 import OffsetArrays
 import LaTeXStrings # only because Aqua complains: used in the Plotting extensions
 
+# Names, shared between different packages
+import MolSimToolkitShared: center_of_mass, 
+                            distances,
+                            coordination_number,
+                            bulk_coordination,
+                            wrap,
+                            wrap_to_first,
+                            align, 
+                            align!, 
+                            rmsd 
+
 using TestItems: @testitem
 using StaticArrays: FieldVector, SMatrix, MVector
 using LinearAlgebra: norm, cross, dot
@@ -19,6 +30,7 @@ export align, align!, rmsd, rmsd_matrix
 export intermittent_correlation
 export bulk_coordination
 export coordination_number
+export center_of_mass
 
 # Reexported from ProteinSecondaryStructures for convenience
 using ProteinSecondaryStructures: dssp_run, stride_run, 
@@ -27,10 +39,6 @@ export dssp_run, stride_run, ss_code, ss_number, ss_name
 # SS trajectory functions
 export ss_map, ss_mean, ss_heatmap
 
-# New method added to this function, which is reexported
-import PDBTools.center_of_mass
-export center_of_mass
-
 # Version of the package: used for printing in some places
 const version = pkgversion(@__MODULE__)
 
@@ -44,9 +52,6 @@ include("../test/Testing.jl")
 include("./datastructures/Simulation.jl")
 include("./datastructures/Positions.jl")
 
-# Coordinate PBC wrapping functions
-include("./wrap.jl")
-
 # Structural properties
 include("./miscelaneous_functions/distances.jl")
 include("./miscelaneous_functions/dihedrals.jl")
@@ -62,7 +67,7 @@ include("./miscelaneous_functions/bulk_coordination.jl")
 include("./miscelaneous_functions/coordination_number.jl")
 
 #  Structural alignment
-include("./procrustes.jl")
+include("./structural_alignment.jl")
 
 # Analysis functions and modules
 include("./BlockAverages.jl")

src/miscelaneous_functions/center_of_mass.jl

@@ -38,7 +38,7 @@ julia> cm = center_of_mass(protein_indices, simulation, coor)
     The `iref=nothing` option was added in version 1.22.0.
 
 """
-function PDBTools.center_of_mass(
+function center_of_mass(
     indices::AbstractVector{<:Integer},
     simulation::Simulation,
     p::FramePositions;
@@ -61,7 +61,7 @@ function PDBTools.center_of_mass(
     return Point3D(cm /= totmass)
 end
 
-@testitem "centerofmass" begin
+@testitem "center_of_mass" begin
     using MolSimToolkit.Testing
     import PDBTools
     simulation = Simulation(Testing.namd_pdb, Testing.namd_traj)

src/procrustes.jl → src/structural_alignment.jl

@@ -1,146 +1,19 @@
-using LinearAlgebra: eigvecs
-using StaticArrays: SVector, MMatrix, SMatrix
-
-_center_of_mass(x::AbstractVector{<:AbstractVector}, ::Nothing) = sum(x) / length(x)
-_center_of_mass(x::AbstractVector{<:AbstractVector}, mass::AbstractVector) =
-    sum(x[i] * mass[i] for i in eachindex(x, mass)) / sum(mass)
-
 """
-    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.
-
-"""
-function align end
-@doc (@doc align) align!
-
-function align(
-    x::AbstractVector{<:AbstractVector},
-    y::AbstractVector{<:AbstractVector};
-    mass=nothing
-)
-    xnew = copy(x)
-    return align!(xnew, y; mass)
-end
-
-function align!(
-    x::AbstractVector{<:AbstractVector},
-    y::AbstractVector{<:AbstractVector};
-    mass=nothing,
-    # Auxiliary arrays that might be preallocated
-    xm=zeros(3, length(x)),
-    xp=zeros(3, length(x))
-)
-    length(x) == length(y) || throw(DimensionMismatch("x and y must have the same length"))
-    (length(x[1]) != 3 || length(x[2]) != 3) && throw(DimensionMismatch("x and y must be 3D vectors"))
-
-    cmx = _center_of_mass(x, mass)
-    cmy = _center_of_mass(y, mass)
-    for i in eachindex(x, y)
-        x[i] -= cmx
-        y[i] -= cmy
-    end
-
-    for i in eachindex(x, y)
-        xm[1:3, i] .= y[i] .- x[i]
-        xp[1:3, i] .= y[i] .+ x[i]
-    end
-
-    q = zeros(MMatrix{4,4,eltype(xm),16})
-    for i in eachindex(x)
-        q[1, 1] = q[1, 1] + sum(abs2, @view(xm[1:3, i]))
-        q[1, 2] = q[1, 2] + xp[2, i] * xm[3, i] - xm[2, i] * xp[3, i]
-        q[1, 3] = q[1, 3] + xm[1, i] * xp[3, i] - xp[1, i] * xm[3, i]
-        q[1, 4] = q[1, 4] + xp[1, i] * xm[2, i] - xm[1, i] * xp[2, i]
-        q[2, 2] = q[2, 2] + xp[2, i]^2 + xp[3, i]^2 + xm[1, i]^2
-        q[2, 3] = q[2, 3] + xm[1, i] * xm[2, i] - xp[1, i] * xp[2, i]
-        q[2, 4] = q[2, 4] + xm[1, i] * xm[3, i] - xp[1, i] * xp[3, i]
-        q[3, 3] = q[3, 3] + xp[1, i]^2 + xp[3, i]^2 + xm[2, i]^2
-        q[3, 4] = q[3, 4] + xm[2, i] * xm[3, i] - xp[2, i] * xp[3, i]
-        q[4, 4] = q[4, 4] + xp[1, i]^2 + xp[2, i]^2 + xm[3, i]^2
-    end
-    q[2, 1] = q[1, 2]
-    q[3, 1] = q[1, 3]
-    q[3, 2] = q[2, 3]
-    q[4, 1] = q[1, 4]
-    q[4, 2] = q[2, 4]
-    q[4, 3] = q[3, 4]
-    q = SMatrix(q)
-
-    # Computing the eigenvectors 'v' of the q matrix
-    v = eigvecs(q)
-
-    # Compute rotation matrix
-    u = zeros(MMatrix{3,3,Float64,9})
-    u[1, 1] = v[1, 1]^2 + v[2, 1]^2 - v[3, 1]^2 - v[4, 1]^2
-    u[1, 2] = 2.0 * (v[2, 1] * v[3, 1] + v[1, 1] * v[4, 1])
-    u[1, 3] = 2.0 * (v[2, 1] * v[4, 1] - v[1, 1] * v[3, 1])
-    u[2, 1] = 2.0 * (v[2, 1] * v[3, 1] - v[1, 1] * v[4, 1])
-    u[2, 2] = v[1, 1]^2 + v[3, 1]^2 - v[2, 1]^2 - v[4, 1]^2
-    u[2, 3] = 2.0 * (v[3, 1] * v[4, 1] + v[1, 1] * v[2, 1])
-    u[3, 1] = 2.0 * (v[2, 1] * v[4, 1] + v[1, 1] * v[3, 1])
-    u[3, 2] = 2.0 * (v[3, 1] * v[4, 1] - v[1, 1] * v[2, 1])
-    u[3, 3] = v[1, 1]^2 + v[4, 1]^2 - v[2, 1]^2 - v[3, 1]^2
-    u = SMatrix(u)
-
-    # Rotate to align x to y 
-    for i in eachindex(x)
-        x[i] = u * x[i]
-    end
-
-    # Move aligned x to the original center of mass of y
-    for i in eachindex(x, y)
-        x[i] += cmy
-        y[i] += cmy
-    end
-
-    return x
-end
-
-"""
-    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.
+Computes the root mean square deviation (RMSD) between two sets of points in 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.
+# Arguments
 
-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:
+- `indices` vector contains the indices of the atoms to be considered. 
+- `mass` argument can be used to provide the mass of the atoms if they are not the same.
+- `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:
-
-```jldoctest
-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
 
 ```jldoctest; filter = r"([0-9]+\\.[0-9]{2})[0-9]+" => s"\\1***"
@@ -172,14 +45,6 @@ julia> rmsd(simulation, cas; reference_frame=:average, show_progress=false)
 ```
 
 """
-function rmsd(x::AbstractVector, y::AbstractVector)
-    rmsd = 0.0
-    for i in eachindex(x, y)
-        rmsd += sum(abs2, x[i] .- y[i])
-    end
-    return sqrt(rmsd / length(x))
-end
-
 function rmsd(
     simulation::Simulation, indices::AbstractVector{Int};
     mass=nothing,
@@ -237,7 +102,7 @@ function rmsd(
     return rmsds
 end
 
-@testitem "procrustes" begin
+@testitem "rmsd" begin
     using MolSimToolkit
     using MolSimToolkit.Testing: namd_pdb, namd_traj
     using StaticArrays: SVector
@@ -319,6 +184,12 @@ 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.
 
+# Returns 
+
+A symetric matrix with the RMSD values between each pair of frames. For example, in 
+a trajectory with 5 frames, the matrix will be a 5x5 matrix with the RMSD values
+between the structures of each pair of frames.
+
 # Example
 
 ```jldoctest; filter = r"([0-9]+\\.[0-9]{2})[0-9]+" => s"\\1***"
@@ -391,6 +262,3 @@ end
 
     @test_throws ArgumentError rmsd_matrix(simulation, cas; mass=[1, 2, 3, 4, 5])
 end
-
-
-