笔者编写了一个统计水分子间氢键个数的小脚本，与大家分享一下。该脚本基于 python 3 编写而成， 可以直接应用于 pymol 软件。脚本的使用方法与在 pymol 中使用其他脚本的方法一样，即把 *.py 文件放到与 PyMOL.exe 相同的目录下，接着 run *.py ，此时 *.py 文件中扩展的命令生效，最后执行该命令。这里，小脚本的使用方法是：1) run sum_hb.py; 2) sum_hb 。需要注意的是，该脚本仅能计算水分子之间的氢键，事实上它能很方便地扩展到含S，N等原子的体系的计算。

       sum_hb.py 文件在附件中( sum_hb.py )，这里也把代码贴出来供大家参考：

# -*- coding: utf-8 -*-

"""

Created on Tue Nov 21 14:54:28 2017

Usage: 1) run sum_hb.py; 2) sum_hb

Author: Mu Liu-Hua, 2017

"""

#声明函数包

import numpy as np

import math

from pymol import cmd, stored

def sum_hb():

#点乘函数

   def dot_product(xyz_1,xyz_2):

       return xyz_1[0] * xyz_2[0] + xyz_1[1] *xyz_2[1] + xyz_1[2] * xyz_2[2]  

   def distant(xyz_O_1,xyz_O_2):

       return math.sqrt(math.pow(xyz_O_1[0] - xyz_O_2[0],2) + math.pow(xyz_O_1[1] \

                        - xyz_O_2[1],2) + math.pow(xyz_O_1[2] - xyz_O_2[2],2))

#计算长度的函数

   def length(xyz_O_1):

       return math.sqrt(math.pow(xyz_O_1[0],2) + math.pow(xyz_O_1[1],2) + math.pow(xyz_O_1[2],2))      

   n_frames = 2

   arguments = ([3.2, 45])#氢键的判别标准

   N_O=[]

   N_H=[]

   stored.list=[]

   cmd.iterate("all","stored.list.append((name))")

   N_element = stored.list

   for index_0, element_name in enumerate(N_element):

       if element_name[0:1] == "O":

           N_O.append(index_0+1)

       if element_name[0:1] =="H":

           N_H.append(index_0+1)

#统计氢键的个数    

   hydrogen_number = np.zeros([n_frames-1])

   frames_index = np.zeros([n_frames-1])    

   for n in range(n_frames-1):

       frames_index[n]= n+1

       N_xyz = cmd.get_coords('all', 1)

       possible_bonds = []

       for index_1 in range(len(N_O)-1):

           xyz_O_1 = np.zeros([1,3])

           O_1 = N_O[index_1]

           xyz_O_1 = N_xyz[O_1-1]

           for index_2 in range(index_1 + 1, len(N_O)):

               xyz_O_2 = np.zeros([1,3])

               O_2 = N_O[index_2]

               xyz_O_2 = N_xyz[O_2-1]

               dist = distant(xyz_O_1,xyz_O_2)

               if dist < arguments[0]:

                   possible_bonds.append((O_1, O_2))

#利用氢原子与氧原子的距离作进一步判断                                      

       possible_bonds_refined = []

       for possible_bond in possible_bonds:

           O1_index = possible_bond[0]

           O2_index = possible_bond[1]

           O1_xyz = N_xyz[O1_index-1]

           O2_xyz = N_xyz[O2_index-1]

           for H_index_1 in range(len(N_H)):

               H_index = N_H[H_index_1]

               H_xyz = N_xyz[H_index-1]

               dist1 = distant(O1_xyz,H_xyz)

               if dist1 < arguments[0]:

                   dist2 = distant(O2_xyz,H_xyz)

                   if dist2 < arguments[0]:

                       if dist1 < 1.3 or dist2 < 1.3:

                           if dist1 < 1.3:

                               possible_bonds_refined.append((H_index, O1_index, O2_index))

                           else:

                               possible_bonds_refined.append((H_index, O2_index, O1_index))

#利用角度作最终判断

       hydrogen_bonds = []

       for possible_bond_refined in possible_bonds_refined:

           hydrogen_index = possible_bond_refined[0]

           bond_donor_index = possible_bond_refined[1]

           bond_acceptor_index = possible_bond_refined[2]

           bond_donor = N_xyz[bond_donor_index-1]

           hydrogen = N_xyz[hydrogen_index-1]

           bond_acceptor = N_xyz[bond_acceptor_index-1]

           vector1 = ([hydrogen[0] - bond_donor[0], hydrogen[1] - \

                       bond_donor[1], hydrogen[2] - bond_donor[2]])

           vector2 = ([bond_acceptor[0] - bond_donor[0], bond_acceptor[1]\

                       - bond_donor[1], bond_acceptor[2] - bond_donor[2]])

           angle = math.acos((dot_product(vector1,vector2))/(length(vector1) * \

                              length(vector2))) * 180 / math.pi # at most 180

           if angle <= arguments[1]:

               hydrogen_bonds.append(possible_bond_refined)

       hydrogen_number[n]=len(hydrogen_bonds)

   print ("hydrogen_bonds are:",hydrogen_bonds)

   print ("hydrogen_number is:",int(hydrogen_number))

cmd.extend("sum_hb", sum_hb)

sum_hb.py

参考网址：

HBonanza: a computer algorithm for molecular-dynamics-trajectory hydrogen-bond analysis. http://rocce-vm0.ucsd.edu/data/sw/hosted/hbonanza/ . 2017.11.21

PyMOL Command Reference. http://pymol.org/pymol-command-ref.html 2017.11.21