RMSD 계산하기 (Kabsch algorithm) - Python

3차원 좌표계에서 Minimum RMSD (Root Mean Squared Deviation)을 구하기 위한 두 가지 알고리즘에 대한 내용이다.

실제로 사용할 때에는 아래 RMSD library를 활용하자. (Kabsch, Quaternion 모두 가능)

rmsd

charnley

Kabsch Algorithm

Quaternion Algorithm (TBU)

RMSD

Kabsch Algorithm

두 구조

X, Y \in \mathbb{R}^{N \times 3}

간의 minimum RMSD를 구할 때, minimum RMSD를 갖도록

X

를 이동시키는 어떤 affine tranformation

RX + \vec{t}

이 있다고 하자.

이때, optimal rotation matrix (

R

)와 translation matrix (

\vec{t}

)를 deterministic하게 구하는 알고리즘이 Kabsch algorithm 이다.

참고로, 순서가 중요하다.

식

RX + \vec{t}

에서 볼 수 있듯이, 먼저 i) rotation matrix 를 구하고 ii) translation matrix를 구한다. 반대로 하면 안 된다!

그 이유는 rotation과 translation이 서로 교환법칙이 성립하지 않기 때문이다. (not commutative)

따라서, 원점(아무리 rotation을 해도 이동하지 않는 지점)으로 중심을 이동하여 rotation matrix를 찾고, rotation 한 뒤 translation matrix를 구한다.

Kabsch Algorithm은 크게 세 단계로 나뉜다.

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
복사

Find Optimal Rotation Matrix

Cross-covariance matrix 계산

H = X^\top Y

(

X, Y

모두

(N,3)

)

# compute cross covariance matrix
h = torch.einsum('ij, ik -> jk', a, b) # (3, 3)
Python
복사

Rotation matrix 계산 (with SVD)

참고) Singular Value Decomposition (SVD)의 기하학적인 의미

H = U\Sigma V^\top

d = \text{sign}(\text{det}(VU^\top))

→ proper rotation일 경우,

d=1

, improper rotation일 경우,

d=-1

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
복사

간단한 증명

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
복사