内坐标转换计算

前言

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

向量夹角的范围

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

二面角的范围

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

严格定义

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

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

例如：

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

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

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

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

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

定义

给定一个有序列(ordered sequence) = [A,B,C,D]。
约定 A,B,C,D 构成的广义二面角为向量$\overrightarrow{BA},\overrightarrow{CB},\overrightarrow{DC}$ 构成的2个平面法向量
${\vec{n}}_1={\vec{n}}_{CBA}$ 与 ${\vec{n}}_2={\vec{n}}_{DCB}$之间， ${\vec{n}}_2$ 逆时针转到 ${\vec{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的法向量是$\overrightarrow{BA}$ （右手定则），
前一个面CBA的法向量是$\overrightarrow{DB}$ （右手定则）。

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

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

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

内坐标转笛卡尔算法

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

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

Rodrigues’旋转公式
设旋转轴$\vec{n}=(n_x,n_y,n_z)$ 。

写出$\vec{n}$的叉乘矩阵:
$N=\left[\begin{array}{ccc}0,& -n_z,& n_y\\ n_z,& 0,& -n_x\\ -n_y,& n_x,& 0\end{array}\right]$。

于是叉乘：
$\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{Np}$

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

旋转后向量
${\vec{p}}^{\prime }=\mathbf{R}\vec{p}$。

或者用向量形式
$\begin{array}{c}{\mathbf{p}}^{\prime }={\mathbf{p}}_{\perp }^{\prime }+{\mathbf{p}}_{\Vert }^{\prime }=\cos (\theta ){\mathbf{p}}_{\perp }+\sin (\theta )(\mathbf{n}\times \mathbf{p})+{\mathbf{p}}_{\Vert }\\ =\cos (\theta )\left(\mathbf{p}-{\mathbf{p}}_{\Vert }\right)+{\mathbf{p}}_{\Vert }+\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{array}$

即:

${\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})$

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

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

欲求 $\overrightarrow{DC}$， 将其分解为 $\overrightarrow{CB}$ 上的分量，和 垂直于 $\overrightarrow{CB}$ 的分量。
(align with , orthogonal to $\overrightarrow{CB}$)

已知距离和bond angle，align with 的分量很好求
${\overrightarrow{DC}}_{align}=d\cos \alpha \frac{\overrightarrow{CB}}{|\overrightarrow{CB}|}$。

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

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

按上文所言，以$\overrightarrow{CB}$为视角，从逆向往正向看，建立平面坐标系。
以$\vec{v}$为该平面上的 $x^{\prime }$ 轴，
约定该平面上 $y^{\prime }$ 轴正向是 $\vec{v}\times \overrightarrow{CB}$ (右手定则)。
$\overrightarrow{DC}$在该平面上的分量即为$\overrightarrow{DC}$垂直于$\overrightarrow{CB}$的分量。

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

即，问题变成了求 $\vec{v}$ 绕 $\overrightarrow{CB}$ 逆时针旋转 {dihedral} 度得到的向量。

应用上文所言的旋转公式
其中 旋转轴为 $\vec{n}=\overrightarrow{CB}$,
$\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})$。

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

旋转后的方向得到了，再考虑长度。
由键长知
${\overrightarrow{DC}}_{ortho}=|d\sin \alpha |\frac{\vec{w}}{|\vec{w}|}$， 由于约定了键角[0,pi), $\sin \alpha >0$。
${\overrightarrow{DC}}_{ortho}=d\sin \alpha \frac{\vec{w}}{|\vec{w}|}$。

最后
$\overrightarrow{DC}={\overrightarrow{DC}}_{align}+{\overrightarrow{DC}}_{ortho}=d\cos \alpha {\vec{u}}_{CB}+d\sin \alpha {\vec{u}}_w$;

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

于是D的坐标 = $\overrightarrow{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。