|
|
本帖最后由 chieko 于 2023-11-29 17:35 编辑
请教一下大家，我想要统计体系中水形成的氢键取向的分布，现在用一个pbc的水盒子测试脚本，5 nm的OPC3水盒子，20 ns NPT：
脚本很简单，主要是先用MDAnalysis统计氢键，然后依次读取每个氢键中供体上H和受体O的坐标，用这两个坐标计算与z方向的角度，最后统计角度分布并输出；
理论上bulk水的氢键角度分布应该是各向同性的，显示一条直线（如J. Chem. Phys. 143, 114708 (2015)， https://doi.org/10.1063/1.4931414）；
但我计算出来的是这样的不均匀的分布↓；（左：我的结果；右：文献的，绿色为bulk）
脚本中用MDA分析氢键的部分：
- hbonds = HydrogenBondAnalysis(universe=u,
- donors_sel=None,
- hydrogens_sel=f'name HW1 HW2 and {center_region} and (prop z >= {z1}) and (prop z <= {z2})', # HBs in water.
- acceptors_sel=f'name OW and {center_region} and (prop z >= {z1}) and (prop z <= {z2})',
- d_a_cutoff=3.5,
- d_h_a_angle_cutoff=138.1, # set as 138.1 to fit gmx (30).
- update_selections=True)
- hbonds.run(start=fps1,stop=fps2,verbose=True)
- for hbond in tqdm(hbonds.results.hbonds, desc="Processing for all hydrogen bonds"):
- frame = int(hbond[0]) # Frame of this hb
- h_idx = int(hbond[2]) # Hydrogen atom index
- a_idx = int(hbond[3]) # Acceptor atom index
- u.trajectory[frame]
- vec = u.atoms.positions[h_idx] - u.atoms.positions[a_idx]
- angle = np.degrees(np.arccos(vec[2] / np.linalg.norm(vec)))
- angles.append(angle)
- angle_distribution = np.array(angles)
麻烦Sob老师和各位老师看看问题出在哪里，谢谢大家！
|
|
发表于 Post on 2023-11-29 17:11:20
收藏 Add to favorites
楼主 Author