写一下计算氢键的脚本代码

以下是计算氢键的 Python 代码：

import numpy as np

def distance(a, b):
    """计算两个原子之间的距离"""
    return np.sqrt(np.sum((a - b) ** 2))

def angle(a, b, c):
    """计算三个原子之间的角度"""
    ba = a - b
    bc = c - b
    cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))
    return np.degrees(np.arccos(cosine_angle))

def calculate_hydrogen_bonds(coordinates, threshold_distance=3.5, threshold_angle=20):
    """计算氢键"""
    num_atoms = len(coordinates)
    hydrogen_bonds = []
    for i in range(num_atoms):
        for j in range(i+1, num_atoms):
            if coordinates[i][] == 'H' or coordinates[j][] == 'H':
                continue
            if coordinates[i][] == 'O' and coordinates[j][] == 'N':
                distance_ij = distance(coordinates[i][1:], coordinates[j][1:])
                if distance_ij < threshold_distance:
                    for k in range(num_atoms):
                        if k == i or k == j or coordinates[k][] == 'H':
                            continue
                        if coordinates[k][] == 'O' and angle(coordinates[k][1:], coordinates[j][1:], coordinates[i][1:]) > threshold_angle:
                            distance_ik = distance(coordinates[i][1:], coordinates[k][1:])
                            if distance_ik < threshold_distance:
                                hydrogen_bonds.append((i, j, k))
            elif coordinates[i][] == 'N' and coordinates[j][] == 'O':
                distance_ij = distance(coordinates[i][1:], coordinates[j][1:])
                if distance_ij < threshold_distance:
                    for k in range(num_atoms):
                        if k == i or k == j or coordinates[k][] == 'H':
                            continue
                        if coordinates[k][] == 'N' and angle(coordinates[k][1:], coordinates[i][1:], coordinates[j][1:]) > threshold_angle:
                            distance_jk = distance(coordinates[j][1:], coordinates[k][1:])
                            if distance_jk < threshold_distance:
                                hydrogen_bonds.append((i, j, k))
    return hydrogen_bonds

注意：这只是一个示例代码，实际应用中需要根据具体情况进行修改和优化。