自己做网站的优势,大网站开发语言,广州白云区网站建设公司,运城做网站公司前言
化学这边的库太多了。 cs这边的库太少了。 去看化学的库太累了。 写一个简单的实现思路#xff0c;让cs的人能看懂。
向量夹角的范围
[0, pi) 这是合理的。 因为两个向量只能构成一个平面系统#xff0c;平面系统内的夹角不能超过pi。
二面角的范围
涉及二面角让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 1n CBA 与 n ⃗ 2 n ⃗ D C B \vec{n}_2\vec{n}_{DCB} n 2n 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 torchcart_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−nx0 。
于是叉乘 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−nx0 pNp
于是绕 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 RIsin(θ)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)nsin(θ)(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)nsinθ(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 anglealign 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 aligndcosα∣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 orthodsinα∣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 alignDC orthodcosαu CBdsinα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_coordtorch.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_coordtorch.tensor([[0.0000e00, 0.0000e00, 0.0000e00],[0.0000e00, 0.0000e00, 1.0000e00],[0.0000e00, 2.0000e00, 1.0000e00],[2.0000e00, 1.1921e-07, 1.0000e00]])其他坐标系的兼容
主流的一些3D库化学库有些用的(z,x,y) 或者 (z, y, x)坐标顺序。 这个也简单。 我们返回的坐标交换一下顺序就能得到其他坐标系了。 output_coord[:, [1,2,0]] - (z,x,y) 坐标。
speed test
cart2internal internal2cart。
1w个样本约30.6s。 100w 样本约 1h。
