使用numpy計(jì)算分子內(nèi)坐標(biāo)
技術(shù)背景
當(dāng)我們打開(kāi)一個(gè)用于表示分子構(gòu)象的xyz文件或者pdb文件，很容易可以理解這種基于笛卡爾坐標(biāo)的空間表征方法。但是除了笛卡爾坐標(biāo)表示方法之外，其實(shí)也有很多其他的方法用于粗?；蛘咂渌康牡谋碚鞣椒ǎ热缜耙黄恼轮兴榻B的在AlphaFold2中所使用的殘基的剛體表示方法。而這種剛體坐標(biāo)，在本質(zhì)上來(lái)說(shuō)也是一種特殊的分子內(nèi)坐標(biāo)表示方法，因?yàn)閷?duì)于每一個(gè)殘基而言只有旋轉(zhuǎn)和平移的自由度，而殘基內(nèi)部是保持互相之間相對(duì)靜止的。換句話說(shuō)，每一個(gè)殘基的內(nèi)坐標(biāo)是保持不變的，本文主要介紹分子的內(nèi)坐標(biāo)表示方法。
具體表示方法
在笛卡爾坐標(biāo)系中，我們使用絕對(duì)坐標(biāo)來(lái)表示每一個(gè)原子的空間位置，雖然也可以用于計(jì)算分子之間的相對(duì)位置，但是如果每一次更新之后都需要重新計(jì)算一遍這個(gè)相對(duì)位置的話，在演化效率上來(lái)說(shuō)會(huì)比較低。因此，我們考慮有沒(méi)有一種方法，可以直接對(duì)分子的“相對(duì)位置”進(jìn)行演化，直到演化結(jié)束之后，再轉(zhuǎn)化回笛卡爾坐標(biāo)進(jìn)行可視化。其實(shí)，市面上已經(jīng)有一些軟件可以直接可視化這種基于“相對(duì)位置”的內(nèi)坐標(biāo)了，這里我們主要探討從絕對(duì)坐標(biāo)到相對(duì)坐標(biāo)的算法。
絕對(duì)坐標(biāo)描述的是B個(gè)體系，每個(gè)體系A(chǔ)個(gè)原子的三維空間坐標(biāo)riri，在python編程中我們可以用張量的維度來(lái)描述這個(gè)三維坐標(biāo)(B, A, D)。
相對(duì)坐標(biāo)的第一個(gè)參量是原子之間的距離li,i+1li,i+1，其張量維度為(B, A-1, 1)。
相對(duì)坐標(biāo)的第二個(gè)參量是原子之間的夾角ai,i+1,i+2ai,i+1,i+2，其張量維度為(B, A-2, 1)。
相對(duì)坐標(biāo)的第三個(gè)參量是原子之間的二面角di,i+1,i+2,i+3di,i+1,i+2,i+3，其張量維度為(B, A-3, 1)。
最后在計(jì)算得到所有的內(nèi)坐標(biāo)參量之后，我們可以用concat的方法把它們按照最后一個(gè)維度進(jìn)行拼接，并且在原子數(shù)維度進(jìn)行擴(kuò)展，最終得到一個(gè)(B, A, 3)的張量，也就是我們所需要的最終的內(nèi)坐標(biāo)。
代碼實(shí)現(xiàn)
其實(shí)這個(gè)算法邏輯是很簡(jiǎn)單的，我們更多的注重一個(gè)原生算子的使用以及代碼的復(fù)用。以下是幾個(gè)相關(guān)的關(guān)注點(diǎn)：
在計(jì)算距離、角度和二面角的過(guò)程中，我們都會(huì)使用到序列原子之間的相對(duì)矢量(B, A-1, D)，那么在計(jì)算過(guò)一次之后我們應(yīng)該保存下來(lái)以供幾個(gè)不同的函數(shù)使用。
在numpy或者是一些常用的深度學(xué)習(xí)框架中，我們最好在代碼實(shí)現(xiàn)階段就去避免x0x0這種情況的出現(xiàn)，一般在遇到除法、反三角函數(shù)或者對(duì)數(shù)函數(shù)的時(shí)候，我們可以在對(duì)應(yīng)的位置加一個(gè)小量?=1e?08?=1e?08以避免出現(xiàn)nan的情況。
在計(jì)算相對(duì)矢量的時(shí)候我們一般使用的是錯(cuò)位相減，比如可以使用crd[1:]-crd[:-1]，但是這里我們?cè)谟?jì)算過(guò)程中使用的是numpy.roll對(duì)數(shù)組進(jìn)行滾動(dòng)之后做減法，最后再去掉一個(gè)結(jié)果。事實(shí)上用前面的這種算法會(huì)更加簡(jiǎn)單高效一些。
# inner_crd.pyimport numpy as npnp.random.seed(1)EPSILON = 1e-08def 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, vecdef 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, vecdef 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) + EPSILONvec_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, disdef 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_crdif __name__ == '__main__':B = 1A = 6D = 3# (B, A, D)origin_crd = np.random.random((B, A, D))# (B, A, 3)icrd = get_inner_crd(origin_crd)print (icrd)
上述代碼執(zhí)行的輸出結(jié)果如下所示：
[[[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]]]
這個(gè)結(jié)果就是我們所需要的分子內(nèi)坐標(biāo)。
總結(jié)概要
本文主要介紹了在numpy的框架下實(shí)現(xiàn)的分子內(nèi)坐標(biāo)的計(jì)算，類(lèi)似的方法可以應(yīng)用于MindSpore和Pytorch、Jax等深度學(xué)習(xí)相關(guān)的框架中。分子的內(nèi)坐標(biāo)，可以更加直觀的描述分子內(nèi)的相對(duì)運(yùn)動(dòng)，通過(guò)鍵長(zhǎng)鍵角和二面角這三個(gè)參數(shù)。
