技術背景

當我們打開一個用于表示分子構象的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]]]

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

總結概要

本文主要介紹了在numpy的框架下實作的分子內坐標的計算，類似的方法可以應用于MindSpore和Pytorch、Jax等深度學習相關的框架中，分子的內坐標，可以更加直觀的描述分子內的相對運動，通過鍵長鍵角和二面角這三個引數，

著作權宣告

本文首發鏈接為：https://www.cnblogs.com/dechinphy/p/inner_crd.html

作者ID：DechinPhy

更多原著文章請參考：https://www.cnblogs.com/dechinphy/

打賞專用鏈接：https://www.cnblogs.com/dechinphy/gallery/image/379634.html

騰訊云專欄同步：https://cloud.tencent.com/developer/column/91958

CSDN同步鏈接：https://blog.csdn.net/baidu_37157624?spm=1008.2028.3001.5343

51CTO同步鏈接：https://blog.51cto.com/u_15561675

參考鏈接

https://zhuanlan.zhihu.com/p/438712498#:~:text=而內坐標通過分子中各原子相對位置,來確定原子位置，因此只要選定參考原子，內坐標系下的分子坐標天生滿足旋轉平移不變性，分子中全部原子的內坐標總稱為“Z-矩陣” (Z-matrix)，

“留一手”加劇內卷，“講不清”浪費時間，

轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/554788.html

標籤：其他

上一篇：邊緣計算簡介