Search

RMSD 계산하기 (Kabsch algorithm) - Python

3차원 좌표계에서 Minimum RMSD (Root Mean Squared Deviation)을 구하기 위한 두 가지 알고리즘에 대한 내용이다.
실제로 사용할 때에는 아래 RMSD library를 활용하자. (Kabsch, Quaternion 모두 가능)
rmsd
charnley

Kabsch Algorithm

두 구조 X,Y∈RN×3X, Y \in \mathbb{R}^{N \times 3} 간의 minimum RMSD를 구할 때, minimum RMSD를 갖도록 XX를 이동시키는 어떤 affine tranformation RX+t⃗RX + \vec{t} 이 있다고 하자.
이때, optimal rotation matrix (RR)와 translation matrix (t⃗\vec{t})를 deterministic하게 구하는 알고리즘이 Kabsch algorithm 이다.
참고로, 순서가 중요하다.
식 RX+t⃗RX + \vec{t} 에서 볼 수 있듯이, 먼저 i) rotation matrix 를 구하고 ii) translation matrix를 구한다. 반대로 하면 안 된다!
그 이유는 rotation과 translation이 서로 교환법칙이 성립하지 않기 때문이다. (not commutative)
따라서, 원점(아무리 rotation을 해도 이동하지 않는 지점)으로 중심을 이동하여 rotation matrix를 찾고, rotation 한 뒤 translation matrix를 구한다.
Kabsch Algorithm은 크게 세 단계로 나뉜다.
1.
Move to Origin
예시를 위해 random point 정의
import torch x = torch.randn(10, 3) # (N, 3) y = torch.randn(10, 3) # (N, 3)
Python
복사
mean_x = x.mean(dim=-2, keepdim=True) mean_y = y.mean(dim=-2, keepdim=True) # move to origin x = x - mean_x y = y - mean_y
Python
복사
2.
Find Optimal Rotation Matrix
a.
Cross-covariance matrix 계산
H=X⊤YH = X^\top Y (X,YX, Y 모두 (N,3)(N,3))
# compute cross covariance matrix h = torch.einsum('ij, ik -> jk', a, b) # (3, 3)
Python
복사
b.
Rotation matrix 계산 (with SVD)
참고) Singular Value Decomposition (SVD)의 기하학적인 의미
H=UΣV⊤H = U\Sigma V^\top
d=sign(det(VU⊤))d = \text{sign}(\text{det}(VU^\top)) → proper rotation일 경우, d=1d=1, improper rotation일 경우, d=−1d=-1
R=V(10001000d)U⊤R = V \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & d\end{pmatrix}U^\top
# compute optimal rotation matrix u, _, vt = torch.svd(h) ut, v = u.transpose(-1, -2), vt.transpose(-1, -2) d = torch.sign(torch.linalg.det(torch.einsum('ij, jk -> ik', v, ut))) diag = torch.eye(3) diag[2, 2] = d rot_matrix = torch.einsum('ij, jk, kl -> il', v, diag, ut)
Python
복사
간단한 증명
3.
Find Optimal Translation Vector
# compute optimal translation vector trans_vec = mean_y.squeeze(-2) - torch.einsum('ij, j -> i', rot_matrix, mean_x.squeeze(-2))
Python
복사

Quaternion Algorithm (TBU)

onlinelibrary.wiley.com
https://onlinelibrary.wiley.com/doi/epdf/10.1002/jcc.20110

RMSD

Kabsch나 Quaternion을 활용해서 align 한 뒤에는, 간단히 point-by-point로 RMSD를 구하면 된다.
def rmsd(X: torch.Tensor, Y:torch.Tensor) -> float: rot_matrix, trans_vec = kabsch(X, Y) X = torch.einsum('ij, jk -> ik', X, rot_matrix) X = X + trans_vec diff = X - Y return torch.sqrt((diff*diff).sum() / X.shape[0])
Python
복사