成键项的能量与力 —— FPGA 实现公式记录
1. 键长 — 2 个原子 \((i,j)\)
几何量
\[\vec r_{ij} = \vec r_i - \vec r_j, \qquad r = \|\vec r_{ij}\|, \qquad \hat r_{ij} = \frac{\vec r_{ij}}{r}\]
能量
\[V_{\text{bond}} = \tfrac12 k_b\,(r - r_0)^2\]
力
\[\boxed{\ \vec F_i = -k_b\,(r - r_0)\,\hat r_{ij}\ } \qquad \vec F_j = -\vec F_i\]
分量形式：
\[F_{i,x} = -k_b\,(r-r_0)\,\frac{x_i - x_j}{r} \qquad (\text{y, z 同理})\]
2. 键角 — 3 个原子 \((i,j,k)\)，\(j\) 为顶点
几何量
\[\vec a = \vec r_i - \vec r_j, \qquad \vec b = \vec r_k - \vec r_j, \qquad \vec n = \vec a \times \vec b\]
\[\boxed{\ \theta = \operatorname{atan2}\big(\|\vec n\|,\ \ \vec a\cdot\vec b\big)\ } \qquad \theta \in [0,\pi]\]
能量
\[V_{\text{angle}} = \tfrac12 k_\theta\,(\theta - \theta_0)^2\]
力 —— 教科书形式
\[\vec F_i = -\frac{k_\theta(\theta-\theta_0)}{r_{ij}\sin\theta}\big(\cos\theta\,\hat r_{ij} - \hat r_{kj}\big), \qquad
\vec F_k = -\frac{k_\theta(\theta-\theta_0)}{r_{kj}\sin\theta}\big(\cos\theta\,\hat r_{kj} - \hat r_{ij}\big)\]
\[\vec F_j = -(\vec F_i + \vec F_k)\]
力 —— 硬件实际采用的形式（严格恒等）
\[\boxed{\ \vec F_i = k_\theta(\theta-\theta_0)\,\frac{\vec n \times \vec a}{\|\vec n\|\,(\vec a\cdot\vec a)}\ }
\qquad
\boxed{\ \vec F_k = k_\theta(\theta-\theta_0)\,\frac{\vec b \times \vec n}{\|\vec n\|\,(\vec b\cdot\vec b)}\ }\]
\[\vec F_j = -(\vec F_i + \vec F_k)\]
等价性证明：展开双重叉积 \((\vec a\times\vec b)\times\vec a = \vec b(\vec a\cdot\vec a) - \vec a(\vec a\cdot\vec b)\)，并用 \(\|\vec n\| = \|a\|\|b\|\sin\theta\)、\(\vec a\cdot\vec b = \|a\|\|b\|\cos\theta\) 代入即得。
3. 二面角 — 4 个原子 \((i,j,k,l)\)
几何量
\[\vec b_1 = \vec r_j - \vec r_i, \qquad \vec b_2 = \vec r_k - \vec r_j, \qquad \vec b_3 = \vec r_l - \vec r_k\]
\[\vec n_1 = \vec b_1 \times \vec b_2, \qquad \vec n_2 = \vec b_2 \times \vec b_3\]
\[\boxed{\ \varphi = \operatorname{atan2}\big(\|\vec b_2\|\,(\vec b_1\cdot\vec n_2),\ \ \vec n_1\cdot\vec n_2\big)\ } \qquad \varphi \in (-\pi,\pi]\]
关于 \(y\) 的化简：原式 \(y = (\vec n_1\times\vec n_2)\cdot\hat b_2\)。展开双重叉积
\(\vec n_1\times\vec n_2 = (\vec b_1\times\vec b_2)\times(\vec b_2\times\vec b_3) = [\vec b_1\cdot(\vec b_2\times\vec b_3)]\,\vec b_2\)（因 \(\vec b_2 \perp \vec b_2\times\vec b_3\)），
于是 \(y = \|\vec b_2\|\,(\vec b_1\cdot\vec n_2)\)。严格恒等，省掉一个 100 bit 除法器、六个 49×49 乘法器和三个 99 bit 减法器。
二面角是有符号的（覆盖整个圆周），所以用 atan2 而非 arccos —— 这里 \(y\) 可正可负，需要象限折叠。这是与键角的关键区别。
能量（周期形式）
\[V_{\text{torsion}} = k_\varphi\,\big[\,1 + \cos(n\varphi - \delta)\,\big]\]
\(n\) 是多重度（C–C 单键取 3，双键取 2），\(\delta\) 是相位，\(k_\varphi\) 是势垒高度。
用余弦而非平方，是因为扭转有多个等价极小，不像键角只有一个 \(\theta_0\)。”\(1+\)” 保证 \(V\ge 0\)，输出符号位可以直接接地。
力
\[\frac{\partial V}{\partial \varphi} = -n\,k_\varphi \sin(n\varphi - \delta)\]
记
\[\vec g_1 = \frac{\|\vec b_2\|}{\|\vec n_1\|^2}\,\vec n_1, \qquad
\vec g_3 = \frac{\|\vec b_2\|}{\|\vec n_2\|^2}\,\vec n_2, \qquad
c_1 = \frac{\vec b_1\cdot\vec b_2}{\|\vec b_2\|^2}, \qquad
c_3 = \frac{\vec b_3\cdot\vec b_2}{\|\vec b_2\|^2}\]
\[\vec g_2 = -c_1\,\vec g_1 - c_3\,\vec g_3\]
则梯度为
\[\frac{\partial\varphi}{\partial \vec r_i} = -\vec g_1, \qquad
\frac{\partial\varphi}{\partial \vec r_j} = \vec g_1 - \vec g_2, \qquad
\frac{\partial\varphi}{\partial \vec r_k} = \vec g_2 - \vec g_3, \qquad
\frac{\partial\varphi}{\partial \vec r_l} = \vec g_3\]
\[\boxed{\ \vec F_x = -\frac{\partial V}{\partial\varphi}\cdot\frac{\partial\varphi}{\partial \vec r_x} = n\,k_\varphi \sin(n\varphi-\delta)\cdot\frac{\partial\varphi}{\partial \vec r_x}\ }\]
\(\vec g_1,\vec g_2,\vec g_3\) 分别是 \(\partial\varphi/\partial\vec b_1,\partial\vec b_2,\partial\vec b_3\)。
由 \(\vec b_1 = \vec r_j-\vec r_i\) 等，链式法则直接给出上面四个式子，
且 \(\sum_x \partial\varphi/\partial\vec r_x = \vec 0\) 自动成立。
⚠️ \(\vec g_2\) 的两个符号都是负的。凭记忆很容易写错（我第一次就写成 \(+c_1\vec g_1 - c_3\vec g_3\)，
有限差分相对误差 45，而力的总和照样为零 —— 平移不变不能用来验证分量是否正确）。