|||
笔者编写了一个统计水分子间氢键个数的小脚本,与大家分享一下。该脚本基于 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)
参考网址:
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
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-10-19 21:29
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社