内坐标转换计算

发布时间:2024年01月04日

前言

化学这边的库太多了。
cs这边的库太少了。
去看化学的库太累了。
写一个简单的实现思路,让cs的人能看懂。

向量夹角的范围

[0, pi)
这是合理的。
因为两个向量只能构成一个平面系统,平面系统内的夹角不能超过pi。

二面角的范围

涉及二面角,说明坐标空间至少是E(3),可以更高维。

严格定义

严格意义上,二面角的定义是2个半平面的夹角。
在严格定义下,二面角的取值范围必然是 [0, pi)

但上面这个定义显然不是我们cs人喜欢的。
因为在E(3)中,如果使用[0,pi)的二面角定义坐标系统,则会有手性问题

例如:

下图中2个D都符合要求。
(图略)

本质上是因为,与某一半平面A的夹角为 [ 0 , π ) [0, \pi) [0,π) 的半平面有2个。
需要一个额外的sign来指示,point处于这两个半平面中的哪一个。

使用 [ ? π , π ) [-\pi, \pi) [?π,π)的广义二面角,可以唯一确定一个三维空间内的相对位置。

这个广义二面角事实上等价于 旋转角 。
旋转角就是采用四元数计算的,范围也是 [ 0 , 2 π ) [0, 2\pi) [0,2π)

两个向量的旋转角,是指从向量p1开始,逆时针旋转,转到向量p2时,所转过的角度, 范围是 0 ~ 360度

定义

给定一个有序列(ordered sequence) = [A,B,C,D]。
约定 A,B,C,D 构成的广义二面角为向量 B A ? , C B ? , D C ? \vec{BA}, \vec{CB}, \vec{DC} BA ,CB ,DC 构成的2个平面法向量
n ? 1 = n ? C B A \vec{n}_1=\vec{n}_{CBA} n 1?=n CBA? n ? 2 = n ? D C B \vec{n}_2=\vec{n}_{DCB} n 2?=n DCB?之间, n ? 2 \vec{n}_{2} n 2? 逆时针转到 n ? 1 \vec{n}_{1} n 1? 的旋转角。

(后一个面的法向量转到前一个面)

可以这样想,在E3空间中,给定2个向量,用右手定则可以得到他们叉乘的法向量。
从该法向量逆向往正向看去,就能建立一个2D平面坐标系。
在该2D平面坐标系上,v2逆时针旋转到v1的夹角是唯一确定的,取值[0,2pi)。

https://zhuanlan.zhihu.com/p/45404840

笛卡尔坐标转内坐标

Convention

在这里插入图片描述

示例输入:

import torch

cart_coord = torch.tensor([
	[1,0,0],[0,0,0],[0,2,0],[0,0,2]
	])

我们期望得到的内坐标

tensor([ 1.0000, 2.0000, 2.8284, 1.5708, 0.7854, -1.5708])

转成易读形式:

1.0000, 2.0000, 2.8284, pi/2, pi/4, -pi/2

键角[0, pi) 。

对最后一个元素,即广义二面角,需要做说明。

在该特殊例子中,
后一个面DCB的法向量是 B A ? \vec{BA} BA (右手定则),
前一个面CBA的法向量是 D B ? \vec{DB} DB (右手定则)。

后者逆时针转到前者的角度是 -pi/2。 ( 取值范围 [-pi, pi) )。

进一步地,我们可以证明,在general case中,
后一个面按右手定则得到的法向量,等于 B A ? \vec{BA} BA 垂直于 C B ? \vec{CB} CB 的分量。
前一个面按右手定则得到的法向量,等于 D C ? \vec{DC} DC 垂直于 C B ? \vec{CB} CB 的分量。

于是二面角 (D,C,B,A)本质上就是 B A ? \vec{BA} BA D C ? \vec{DC} DC C B ? \vec{CB} CB 垂直平面上的分量之间的夹角。
理解这个性质对下面的内容会有帮助。

内坐标转笛卡尔算法

前置知识: 向量绕任意轴旋转矩阵

https://zhuanlan.zhihu.com/p/380237903
https://blog.csdn.net/FreeSouthS/article/details/112576370
https://zhuanlan.zhihu.com/p/56587491

Rodrigues’旋转公式
设旋转轴 n ? = ( n x , n y , n z ) \vec{n}=(n_x,n_y,n_z) n =(nx?,ny?,nz?)

写出 n ? \vec{n} n 的叉乘矩阵:
N = [ 0 , ? n z , n y n z , 0 , ? n x ? n y , n x , 0 ] N=\left[\begin{array}{ccc}0, & -n_z, & n_y \\ n_z, & 0, & -n_x \\ -n_y, & n_x, & 0\end{array}\right] N= ?0,nz?,?ny?,??nz?,0,nx?,?ny??nx?0? ?

于是叉乘:
n × p = [ 0 , ? n z , n y n z , 0 , ? n x ? n y , n x , 0 ] p = N p \mathbf{n} \times \mathbf{p}=\left[\begin{array}{ccc}0, & -n_z, & n_y \\ n_z, & 0, & -n_x \\ -n_y, & n_x, & 0\end{array}\right] \mathbf{p}=\mathbf{N} \mathbf{p} n×p= ?0,nz?,?ny?,??nz?,0,nx?,?ny??nx?0? ?p=Np

于是绕 n ? \vec{n} n 的旋转矩阵为:
R = I + sin ? ( θ ) N + ( 1 ? cos ? ( θ ) ) N 2 \mathbf{R}=\mathbf{I}+\sin (\theta) \mathbf{N}+(1-\cos (\theta)) \mathbf{N}^2 R=I+sin(θ)N+(1?cos(θ))N2

旋转后向量
p ? ′ = R p ? \vec{p}'=\mathbf{R}\vec{p} p ?=Rp ?

或者用向量形式
p ′ = p ⊥ ′ + p ∥ ′ = cos ? ( θ ) p ⊥ + sin ? ( θ ) ( n × p ) + p ∥ = cos ? ( θ ) ( p ? p ∥ ) + p ∥ + sin ? ( θ ) ( n × p ) = cos ? ( θ ) p + ( 1 ? cos ? ( θ ) ) ( n ? p ) n + sin ? ( θ ) ( n × p ) \begin{gathered}\mathbf{p}^{\prime}=\mathbf{p}_{\perp}^{\prime}+\mathbf{p}_{\|}^{\prime}=\cos (\theta) \mathbf{p}_{\perp}+\sin (\theta)(\mathbf{n} \times \mathbf{p})+\mathbf{p}_{\|} \\ =\cos (\theta)\left(\mathbf{p}-\mathbf{p}_{\|}\right)+\mathbf{p}_{\|}+\sin (\theta)(\mathbf{n} \times \mathbf{p}) \\ =\cos (\theta) \mathbf{p}+(1-\cos (\theta))(\mathbf{n} \cdot \mathbf{p}) \mathbf{n}+\sin (\theta)(\mathbf{n} \times \mathbf{p})\end{gathered} p=p?+p?=cos(θ)p?+sin(θ)(n×p)+p?=cos(θ)(p?p?)+p?+sin(θ)(n×p)=cos(θ)p+(1?cos(θ))(n?p)n+sin(θ)(n×p)?

即:

p ′ = cos ? θ ( p ? ( n ? p ) n ) + ( n ? p ) n + sin ? θ ( n × p ) \mathbf{p}^{\prime}=\cos\theta (\mathbf{p}-(\mathbf{n} \cdot \mathbf{p}) \mathbf{n})+(\mathbf{n} \cdot \mathbf{p}) \mathbf{n}+\sin\theta(\mathbf{n} \times \mathbf{p}) p=cosθ(p?(n?p)n)+(n?p)n+sinθ(n×p)

写成这种形式,因为计算机里面算三角函数的开销大于坐标运算。

内坐标转笛卡尔坐标;二面角复原

欲求 D C ? \vec{DC} DC , 将其分解为 C B ? \vec{CB} CB 上的分量,和 垂直于 C B ? \vec{CB} CB 的分量。
(align with , orthogonal to C B ? \vec{CB} CB )

已知距离和bond angle,align with 的分量很好求
D C ? a l i g n = d cos ? α C B ? ∣ C B ? ∣ \vec{DC}_{align} = d\cos \alpha\frac{\vec{CB}}{|\vec{CB}|} DC align?=dcosαCB CB ?

复杂一点的是垂直分量。
先计算 B A ? \vec{BA} BA C B ? \vec{CB} CB 上的分量。
t ? = B A ? ? C B ? ∣ C B ? ∣ C B ? ∣ C B ? ∣ = ( B A ? ? u ? C B ) u ? C B \vec{t}= \frac{\vec{BA}\cdot \vec{CB}}{|\vec{CB}|} \frac{\vec{CB}}{|\vec{CB}|} = (\vec{BA}\cdot \vec{u}_{CB}) \vec{u}_{CB} t =CB BA ?CB ?CB CB ?=(BA ?u CB?)u CB? , u ? \vec{u} u 表示单位向量。
然后得到 B A ? \vec{BA} BA 垂直于 C B ? \vec{CB} CB 的分量
v ? = B A ? ? t ? \vec{v}=\vec{BA} - \vec{t} v =BA ?t

从后文可以知道,这一步如果把 B A ? \vec{BA} BA 归一化为单位向量也无妨,因为我们在乎的只有方向。

按上文所言,以 C B ? \vec{CB} CB 为视角,从逆向往正向看,建立平面坐标系。
v ? \vec{v} v 为该平面上的 x ′ x' x 轴,
约定该平面上 y ′ y' y 轴正向是 v ? × C B ? \vec{v}\times\vec{CB} v ×CB (右手定则)。
D C ? \vec{DC} DC 在该平面上的分量即为 D C ? \vec{DC} DC 垂直于 C B ? \vec{CB} CB 的分量。

根据上文的旋转角(广义二面角)定义, D C ? \vec{DC} DC 在该平面上的分量,
即为 v ? \vec{v} v 逆时针旋转{dihedral}度得到的向量。

即,问题变成了求 v ? \vec{v} v C B ? \vec{CB} CB 逆时针旋转 {dihedral} 度得到的向量。

应用上文所言的旋转公式
其中 旋转轴为 n ? = C B ? \vec{n}=\vec{CB} n =CB ,
w ? = cos ? θ ( v ? ? ( n ? ? v ? ) n ? ) + ( n ? ? v ? ) n ? + sin ? θ ( n ? × v ? ) \vec{w} = \cos \theta (\vec{v}-(\vec{n} \cdot \vec{v})\vec{n}) + (\vec{n}\cdot \vec{v})\vec{n} +\sin\theta (\vec{n}\times \vec{v}) w =cosθ(v ?(n ?v )n )+(n ?v )n +sinθ(n ×v )

由于按照定义 v ? \vec{v} v 是垂直于 C B ? \vec{CB} CB 的分量,故两者内积为0。
上式进一步化简为
w ? = cos ? θ v ? + sin ? θ ( n ? × v ? ) \vec{w}=\cos\theta \vec{v} +\sin\theta (\vec{n}\times \vec{v}) w =cosθv +sinθ(n ×v )

旋转后的方向得到了,再考虑长度。
由键长知
D C ? o r t h o = ∣ d sin ? α ∣ w ? ∣ w ? ∣ \vec{DC}_{ortho}=|d\sin\alpha| \frac{\vec{w} }{|\vec{w} |} DC ortho?=dsinαw w ?, 由于约定了键角[0,pi), sin ? α > 0 \sin \alpha >0 sinα>0
D C ? o r t h o = d sin ? α w ? ∣ w ? ∣ \vec{DC}_{ortho}=d\sin\alpha \frac{\vec{w} }{|\vec{w} |} DC ortho?=dsinαw w ?

最后
D C ? = D C ? a l i g n + D C ? o r t h o = d cos ? α u ? C B + d sin ? α u ? w \vec{DC}=\vec{DC}_{align} + \vec{DC}_{ortho} = d\cos \alpha \vec{u}_{CB} + d\sin\alpha \vec{u}_{w} DC =DC align?+DC ortho?=dcosαu CB?+dsinαu w?;

u ? \vec{u} u 表示单位向量。

于是D的坐标 = D C ? + C \vec{DC}+C DC +C

前三个点处理

第一个点,按习惯固定(0,0,0)
第二个点,我看大部分库都默认放到z-axis上。于是(0, 0, dst)。
第三个点,按照距离和键角可以获得一个圆锥,约定第三个点放在zoy平面的y轴正半面上。

或者可以用兼容第4个点的方式说,
在y轴正向有一个假想点y, 于是(C,B,A,y) 构成的二面角为0。

注意!

本文默认坐标是(x,y,z)顺序。

测试样例

下图为例,左边是输入坐标。
右边是我们希望还原得到的坐标。
A在原点,B在z轴,c在zoy平面。
这等价于对原输入做一次平移和一次90度旋转。
在这里插入图片描述

# test
#输入
input_coord=torch.tensor([
	[1,0,0],[0,0,0],[0,2,0],[0,0,2]
	])
# 内坐标
inner_coord = torch.tensor(
	[1.0000,  2.0000,  2.8284,  1.5708,  0.7854, -1.5708]
)
# 还原坐标
output_coord=torch.tensor([
    [0.0000e+00, 0.0000e+00, 0.0000e+00],
    [0.0000e+00, 0.0000e+00, 1.0000e+00],
    [0.0000e+00, 2.0000e+00, 1.0000e+00],
    [2.0000e+00, 1.1921e-07, 1.0000e+00]])

其他坐标系的兼容

主流的一些3D库,化学库,有些用的(z,x,y) 或者 (z, y, x)坐标顺序。
这个也简单。
我们返回的坐标交换一下顺序就能得到其他坐标系了。
output_coord[:, [1,2,0]] -> (z,x,y) 坐标。

speed test

cart2internal + internal2cart。

1w个样本约30.6s。
100w 样本约 1h。

文章来源:https://blog.csdn.net/w55100/article/details/134690602
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。