技术背景
当我们打开一个用于表示分子构象的xyz文件或者pdb文件,很容易可以理解这种基于笛卡尔坐标的空间表征方法。但是除了笛卡尔坐标表示方法之外,其实也有很多其他的方法用于粗粒化或者其他目的的表征方法,比如前一篇文章中所介绍的在AlphaFold2中所使用的残基的刚体表示方法。而这种刚体坐标,在本质上来说也是一种特殊的分子内坐标表示方法,因为对于每一个残基而言只有旋转和平移的自由度,而残基内部是保持互相之间相对静止的。换句话说,每一个残基的内坐标是保持不变的,本文主要介绍分子的内坐标表示方法。
具体表示方法
代码实现
其实这个算法逻辑是很简单的,我们更多的注重一个原生算子的使用以及代码的复用。以下是几个相关的关注点:
- 在计算距离、角度和二面角的过程中,我们都会使用到序列原子之间的相对矢量(B, A-1, D),那么在计算过一次之后我们应该保存下来以供几个不同的函数使用。
- 在numpy或者是一些常用的深度学习框架中,我们最好在代码实现阶段就去避免
这种情况的出现,一般在遇到除法、反三角函数或者对数函数的时候,我们可以在对应的位置加一个小量
以避免出现nan的情况。
- 在计算相对矢量的时候我们一般使用的是错位相减,比如可以使用crd1:-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, 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) + 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, 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_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)。