而这里寻找角平分线的算法也并不难，在python中只需要把ij→ij→和ik→ik→这两个矢量归一化之后，相加再取一个相反反向，再归一化一次，就可以得到需要添加氢原子的位置了。需要注意的是，这里长度单位使用的是埃而不是纳米，因此找到长度为1埃的氢键即可。

if type == 'c6': left_arrow = crd[j] - crd[i] left_arrow /= np.linalg.norm(left_arrow) right_arrow = crd[k] - crd[i] right_arrow /= np.linalg.norm(right_arrow) h_arrow = -1 * (left_arrow right_arrow) h_arrow /= np.linalg.norm(h_arrow) 2. dihedral--二面角补氢

在上一个环结构中氢原子所连接的重原子，同时与其他另外两个原子相连接，而本章节中二面角形式的氢原子所连接的重原子，只与一个另外的重原子相连。同样的我们需要保障补充的氢原子跟这三个重原子处在同一个平面内。然后保障二面角的中心旋转对称性，就可以找到需要添加氢原子的位置。

因为依然是在同一个平面内进行处理，因此也有比较简单的操作可以实现，相应的python代码如下：

if type == 'dihedral': h_arrow = crd[j] - crd[k] h_arrow /= np.linalg.norm(h_arrow)

这里使用的方法就是把向量jk→jk→平移了一下，保持一样的角度即可。跟GROMACS中的方法略有点不同的是，在这里我们不追求一定达到109度的角度，更多的是使用系统本身自带的向量来进行操作。当然，这样实现的前提是，我们所得到的蛋白质的重原子结构，已经是一个相对比较合理的构象，否则按照本算法实施也会出现一些问题。

同样的还是一个平面内的操作，类似于乙烯的结构，我们在得到i、j、k、l这4个位置的原子坐标之后就可以通过平移来得到需要补充的两个氢原子的位置。

这里为了跟其他几种补氢的算法保持输入维度的一致，我们并没有直接将4个重原子的位置作为输入，而是取了其中的3个原子，再通过这3个原子的位置去推导氢原子的位置，相应python代码如下：

if type == 'c2h4': h_arrow_1 = crd[j] - crd[k] h1 = (h_arrow_1/np.linalg.norm(h_arrow_1) crd[i]) middle_arrow = (crd[i] - crd[j]) middle_arrow /= np.linalg.norm(middle_arrow) middle_arrow *= np.linalg.norm(h_arrow_1) h_arrow_2 = -h_arrow_1 middle_arrow h2 = (h_arrow_2/np.linalg.norm(h_arrow_2) crd[i])

在具体算法中，第一个氢的位置可以用跟上一个章节中二面角一样的算法推导出来，而第二个氢的位置，我们是假设这个乙烯结构中每一条边的长度都大致相等，这样根据等边三角形的矢量闭环关系，可以推导出来第二个氢原子的位置。

一个sp3杂化的碳原子可以连接其他的4个原子，而甲基−CH3−CH3是很常见的一种基团，我们需要在这个碳上面补3个氢原子。那么此时除了碳原子的位置，和碳原子直接相连的一个原子的位置之外，我们还需要一个额外的原子，就是碳原子的次近邻原子，这样我们就有了三个原子可以去构造一个平面。

因为需要补氢的数量有3个，因此整体上算法会相对复杂一些。首先，补第一个氢原子位置时，可以参考二面角的补法，直接补上一个氢原子。补第二个和第三个氢原子时，需要用到一个平移旋转矩阵，其中又分为三个步骤：平移ij→ij→到Z轴上->分别旋转第一个补的氢120度和240度->平移回原始的位置。相关python代码如下所示：

if type == 'ch3': upper_arrow = crd[k] - crd[j] upper_arrow /= np.linalg.norm(upper_arrow) h1 = -upper_arrow crd[i] axes = crd[j] - crd[i] rotate_matrix = rotate_by_axis(axes, 2 * np.pi / 3) h2 = np.dot(rotate_matrix, h1-crd[i]) h2 /= np.linalg.norm(h2) h2 = crd[i] rotate_matrix = rotate_by_axis(axes, 4 * np.pi / 3) h3 = np.dot(rotate_matrix, h1-crd[i]) h3 /= np.linalg.norm(h3) h3 = crd[i]

这里还用到了一个旋转矩阵，其形式为：

x′=(vx⋅vx⋅(1−cosθ) cosθ)⋅x (vx⋅vy⋅(1−cosθ)−vz⋅sinθ)⋅y (vx⋅vz⋅(1−cosθ) vy⋅sinθ)⋅zy′=(vx⋅vy⋅(1−cosθ) vz⋅sinθ)⋅x (vy⋅vy⋅(1−cosθ) cosθ)⋅y (vy⋅vz⋅(1−cosθ)−vx⋅sinθ)⋅zz′=(vx⋅vz⋅(1−cosθ)−vy⋅sinθ)⋅x (vy⋅vz⋅(1−cosθ) vx⋅sinθ)⋅y (vz⋅vz⋅(1−cosθ) cosθ)⋅zx′=(vx⋅vx⋅(1−cosθ) cosθ)⋅x (vx⋅vy⋅(1−cosθ)−vz⋅sinθ)⋅y (vx⋅vz⋅(1−cosθ) vy⋅sinθ)⋅zy′=(vx⋅vy⋅(1−cosθ) vz⋅sinθ)⋅x (vy⋅vy⋅(1−cosθ) cosθ)⋅y (vy⋅vz⋅(1−cosθ)−vx⋅sinθ)⋅zz′=(vx⋅vz⋅(1−cos⁡θ)−vy⋅sin⁡θ)⋅x (vy⋅vz⋅(1−cos⁡θ) vx⋅sin⁡θ)⋅y (vz⋅vz⋅(1−cos⁡θ) cos⁡θ)⋅z

相应的python代码实现为：

def rotate_by_axis(axis, theta): """Rotate an atom by a given axis with angle theta. Args: axis: The rotate axis. theta: The rotate angle. Returns: The rotate matrix. """ vx, vy, vz = axis[0], axis[1], axis[2] return np.array([[vx*vx*(1-np.cos(theta)) np.cos(theta), vx*vy*(1-np.cos(theta))-vz*np.sin(theta), vx*vz*(1-np.cos(theta)) vy*np.sin(theta)], [vx*vy*(1-np.cos(theta)) vz*np.sin(theta), vy*vy*(1-np.cos(theta)) np.cos(theta), vy*vz*(1-np.cos(theta))-vx*np.sin(theta)], [vx*vz*(1-np.cos(theta))-vy*np.sin(theta), vy*vz*(1-np.cos(theta)) vx*np.sin(theta), vz*vz*(1-np.cos(theta)) np.cos(theta)]])

类似的，这种算法对于初始构象有要求，如果是一个无序混乱的原子系统，则是没有办法通过这种算法来加氢的。

还是sp3杂化的碳原子，但是此时该碳原子已经跟其他三个重原子成键，因此有一个多余的键可以跟氢原子结合生成氢键。由于sp3杂化的特殊性，形成的结构会是一个接近于正四面体的形状。此时同样为了保持接口一致，我们选取包含sp3碳原子在内的一共3个原子，来定位氢原子的位置。