技術背景

當我們開啟一個用於表示分子構象的xyz檔案或者pdb檔案，很容易可以理解這種基於笛卡爾座標的空間表徵方法。但是除了笛卡爾座標表示方法之外，其實也有很多其他的方法用於粗粒化或者其他目的的表徵方法，比如前一篇文章中所介紹的在AlphaFold2中所使用的殘基的剛體表示方法。而這種剛體座標，在本質上來說也是一種特殊的分子內座標表示方法，因為對於每一個殘基而言只有旋轉和平移的自由度，而殘基內部是保持互相之間相對靜止的。換句話說，每一個殘基的內座標是保持不變的，本文主要介紹分子的內座標表示方法。

具體表示方法

在笛卡爾座標系中，我們使用絕對座標來表示每一個原子的空間位置，雖然也可以用於計算分子之間的相對位置，但是如果每一次更新之後都需要重新計算一遍這個相對位置的話，在演化效率上來說會比較低。因此，我們考慮有沒有一種方法，可以直接對分子的「相對位置」進行演化，直到演化結束之後，再轉化回笛卡爾座標進行視覺化。其實，市面上已經有一些軟體可以直接視覺化這種基於「相對位置」的內座標了，這裡我們主要探討從絕對座標到相對座標的演演算法。

絕對座標描述的是B個體系，每個體系A個原子的三維空間座標\(\textbf{r}_i\)，在python程式設計中我們可以用張量的維度來描述這個三維座標(B, A, D)。
相對座標的第一個參量是原子之間的距離\(l_{i,i+1}\)，其張量維度為(B, A-1, 1)。
相對座標的第二個參量是原子之間的夾角\(a_{i,i+1,i+2}\)，其張量維度為(B, A-2, 1)。
相對座標的第三個參量是原子之間的二面角\(d_{i,i+1,i+2,i+3}\)，其張量維度為(B, A-3, 1)。

最後在計算得到所有的內座標參量之後，我們可以用concat的方法把它們按照最後一個維度進行拼接，並且在原子數維度進行擴充套件，最終得到一個(B, A, 3)的張量，也就是我們所需要的最終的內座標。

程式碼實現

其實這個演演算法邏輯是很簡單的，我們更多的注重一個原生運算元的使用以及程式碼的複用。以下是幾個相關的關注點：

在計算距離、角度和二面角的過程中，我們都會使用到序列原子之間的相對向量(B, A-1, D)，那麼在計算過一次之後我們應該儲存下來以供幾個不同的函數使用。
在numpy或者是一些常用的深度學習框架中，我們最好在程式碼實現階段就去避免\(\frac{x}{0}\)這種情況的出現，一般在遇到除法、反三角函數或者對數函數的時候，我們可以在對應的位置加一個小量\(\epsilon=1e-08\)以避免出現nan的情況。
在計算相對向量的時候我們一般使用的是錯位相減，比如可以使用crd[1:]-crd[:-1]，但是這裡我們在計算過程中使用的是numpy.roll對陣列進行卷動之後做減法，最後再去掉一個結果。事實上用前面的這種演演算法會更加簡單高效一些。

# inner_crd.py
import numpy as np
np.random.seed(1)
EPSILON = 1e-08

def get_vec(crd):
    """ Get the vector of the sequential coordinate.
    """
    # (B, A, D)
    crd_ = np.roll(crd, -1, axis=-2)
    vec = crd_ - crd
    # (B, A-1, D)
    return vec[:, :-1, :]

def get_dis(crd):
    """ Get the distance of the sequential coordinate.
    """
    # (B, A-1, D)
    vec = get_vec(crd)
    # (B, A-1, 1)
    dis = np.linalg.norm(vec, axis=-1, keepdims=True)
    return dis, vec

def get_angle(crd):
    """ Get the bond angle of the sequential coordinate.
    """
    # (B, A-1, 1), (B, A-1, D)
    dis, vec = get_dis(crd)
    vec_ = np.roll(vec, -1, axis=-2)
    dis_ = np.roll(dis, -1, axis=-2)
    # (B, A-1, 1)
    angle = np.einsum('ijk,ijk->ij', vec, vec_)[..., None] / (dis * dis_ + EPSILON)
    # (B, A-2, 1), (B, A-1, 1), (B, A-1, D)
    return np.arccos(angle[:, :-1, :]), dis, vec

def get_dihedral(crd):
    """ Get the dihedrals of the sequential coordinate.
    """
    # (B, A-2, 1), (B, A-1, 1), (B, A-1, D)
    angle, dis, vec_0 = get_angle(crd)
    # (B, A-1, D)
    vec_1 = np.roll(vec_0, -1, axis=-2)
    vec_2 = np.roll(vec_1, -1, axis=-2)
    vec_01 = np.cross(vec_0, vec_1)
    vec_12 = np.cross(vec_1, vec_2)
    vec_01 /= np.linalg.norm(vec_01, axis=-1, keepdims=True) + EPSILON
    vec_12 /= np.linalg.norm(vec_12, axis=-1, keepdims=True) + EPSILON
    # (B, A-1, 1)
    dihedral = np.einsum('ijk,ijk->ij', vec_01, vec_12)[..., None]
    # (B, A-3, 1), (B, A-2, 1), (B, A-1, 1)
    return np.arccos(dihedral[:, :-2, :]), angle, dis

def get_inner_crd(crd):
    """ Concat the distance, angles and dihedrals to get the inner coordinate.
    """
    # (B, A-3, 1), (B, A-2, 1), (B, A-1, 1)
    dihedral, angle, dis = get_dihedral(crd)
    # (B, A, 1)
    dihedral_ = np.pad(dihedral, ((0, 0), (3, 0), (0, 0)), mode='constant', constant_values=0)
    angle_ = np.pad(angle, ((0, 0), (2, 0), (0, 0)), mode='constant', constant_values=0)
    dis_ = np.pad(dis, ((0, 0), (1, 0), (0, 0)), mode='constant', constant_values=0)
    # (B, A, 3)
    inner_crd = np.concatenate((dis_, angle_, dihedral_), axis=-1)
    return inner_crd

if __name__ == '__main__':
    B = 1
    A = 6
    D = 3
    # (B, A, D)
    origin_crd = np.random.random((B, A, D))
    # (B, A, 3)
    icrd = get_inner_crd(origin_crd)
    print (icrd)

上述程式碼執行的輸出結果如下所示：

[[[0.         0.         0.        ]
  [0.59214856 0.         0.        ]
  [0.38167145 1.89801242 0.        ]
  [0.46143538 1.2138982  1.46589893]
  [0.86899521 2.32255675 1.61009033]
  [0.84368274 2.92999231 1.97853456]]]

這個結果就是我們所需要的分子內座標。