Create kabsch.jl

[cossio]

Set of functions to compute the least RMSD between two sets of points, up to rigid motions. If you have two sets of points P, Q, you should first center them both to the origin:

center!(P) center!(Q)

Next, rotate one into the other:

rotatedP = zeros(size(P)) kabsch!(rotatedP, P, Q)

Finally compute the RMSD:

rmsd(rotatedP, Q)

[diegozea]

Many Thanks @cossio ! This looks fine, I will merged when Travis and AppVeyor finish the tests. Could you please recommend me or add some tests for that functions?

[coveralls]

Coverage decreased (-0.5%) to 82.145% when pulling 085e1c198eb8b6e0f91d2717cab4475c7c5f3b3d on cossio:patch-1 into 522022eacb7c52ec2aad8eab27f044a02d759a30 on diegozea:master.

[diegozea]

I merged it! If you have special test cases, please tell me. I will try to add some tests later.

Muchas gracias por la contribución @cossio ! :D

[cossio]

I would use another package, like R's Bio3D, to calculate the RMSD of some random structures, and compare. I did something similar for a few structures, but didn't save it. If you do it, let me know if it passes the checks.

[cossio]

@diegozea por nada!

[diegozea]

Ok. I'm going to do that later today using Bio3D without adding the dependency. Thanks @cossio !

[diegozea]

@cossio How can I rotate the Calpha coordinates and use that rotation to rotate the rest of the structure atoms?

[cossio]

@diegozea You'd have to modify my code, so that it computes the rotation matrix using only the Ca coordinates. Then you apply the rotation matrix to the full structure

[diegozea]

@cossio i.e. I want to rotate csa, and use that rotation to transform ca

a = [ 1  1 0
      2  1 0
      2 .9 0 ] # 3 x 3
b = [  1 1 0 
    1+cos(pi/4) 1+sin(pi/4)  0 ]
center!(b)
sa = a[1:2,:] # 2 x 3
csa = sa .- mean(sa,1)
ca  = a  .- mean(sa,1)
r = Array(Float64,2,3)
kabsch!(r, csa, b)

[cossio]

@diegozea Just modify the kabsch! I wrote to return the rotation matrix, which is u * χ * v'

[diegozea]

I did it, but it doesn't work

[cossio]

function kabsch(X::Matrix{Float64}, Y::Matrix{Float64})
    @assert size(X) == size(Y)
    A::Matrix{Float64} = X' * Y
    χ = eye(A)
    χ[end,end] = sign(det(A))
    u::Matrix{Float64}, σ::Vector{Float64}, v::Matrix{Float64} = svd(A)
    return u * χ * v'
end

[diegozea]

The last point of csa is collapsed to get the same coordinates than the second point.

[diegozea]

julia> a = [ 1  1 0
             2  1 0
             2 .9 0 ]
3x3 Array{Float64,2}:
 1.0  1.0  0.0
 2.0  1.0  0.0
 2.0  0.9  0.0

julia> b = [  1 1 0 
           1+cos(pi/4) 1+sin(pi/4)  0 ]
2x3 Array{Float64,2}:
 1.0      1.0      0.0
 1.70711  1.70711  0.0

julia> center!(b)

julia> sa = a[1:2,:]
2x3 Array{Float64,2}:
 1.0  1.0  0.0
 2.0  1.0  0.0

julia> csa = sa .- mean(sa,1)
2x3 Array{Float64,2}:
 -0.5  0.0  0.0
  0.5  0.0  0.0

julia> ca  = a  .- mean(sa,1)
3x3 Array{Float64,2}:
 -0.5   0.0  0.0
  0.5   0.0  0.0
  0.5  -0.1  0.0

julia> function kabsch!(rotatedX::Matrix{Float64}, X::Matrix{Float64}, Y::Matrix{Float64})
         @assert size(rotatedX) == size(X) == size(Y)
         A::Matrix{Float64} = X' * Y
         χ = eye(A)
         χ[end,end] = sign(det(A))
         u::Matrix{Float64}, σ::Vector{Float64}, v::Matrix{Float64} = svd(A)
         A_mul_B!(rotatedX, X, u * χ * v')
         u * χ * v'
       end
kabsch! (generic function with 1 method)

julia> r = Array(Float64,2,3)
2x3 Array{Float64,2}:
 4.24399e-314  0.0  0.0
 1.061e-314    0.0  0.0

julia> f = kabsch!(r, csa, b)
3x3 Array{Float64,2}:
 0.707107  0.707107  0.0
 0.0       0.0       0.0
 0.0       0.0       1.0

julia> r
2x3 Array{Float64,2}:
 -0.353553  -0.353553  0.0
  0.353553   0.353553  0.0

julia> R = Array(Float64, 3, 3)
3x3 Array{Float64,2}:
 6.92975e-310  6.92975e-310  6.92975e-310
 2.97128e-316  2.97129e-316  2.97129e-316
 2.97129e-316  6.92975e-310  2.97129e-316

julia> A_mul_B!(R, ca, f)
3x3 Array{Float64,2}:
 -0.353553  -0.353553  0.0
  0.353553   0.353553  0.0
  0.353553   0.353553  0.0

[diegozea]

@cossio Do you know why are (0.5,0.0,0.0) and (0.5,-0.1,0.0) transformed into (0.353553,0.353553,0.0)... Is it because u * χ * v' second row is full of 0.0?

[diegozea]

@cossio I'm using the new version you post here, but I'm still getting the same error. Is it a problem of my example? Am I centering the matrix in the correct way? Thanks!

[diegozea]

Ok, I think I solve it :)

[diegozea]

@cossio It's solved! Thanks! :D

[cossio]

@diegozea Great!