计算化学公社

标题: 统计体相水氢键z轴取向的分布，得到的分布却显示各向异性 [打印本页]

作者
Author: chieko 时间: 2023-11-29 17:11
标题: 统计体相水氢键z轴取向的分布，得到的分布却显示各向异性
本帖最后由 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）
(, 下载次数 Times of downloads: 12)
脚本中用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老师和各位老师看看问题出在哪里，谢谢大家！

作者
Author: Lacrimosa 时间: 2023-11-29 18:05
本帖最后由 Lacrimosa 于 2023-11-29 18:14 编辑

可能是归一化的问题，把结果除以一个sin(theta)看看

作者
Author: chieko 时间: 2023-11-29 21:12
本帖最后由 chieko 于 2023-11-30 10:13 编辑

Lacrimosa 发表于 2023-11-29 18:05
可能是归一化的问题，把结果除以一个sin(theta)看看

非常感谢您的建议！
我尝试了一下，数据确实平了，但是还想请教两个问题：
1、这样得到的值大约在0.88%左右，理论上似乎应该是1/180=0.00556即0.55%左右，如文献中显示那样；
2、为什么您会想到除以sin(theta)呢？是只因为形状相似（笑），还是有什么地方我没有注意到呢？或许这也有助于我解决问题1；
再次感谢您！

作者
Author: Lacrimosa 时间: 2023-11-30 12:21

chieko 发表于 2023-11-29 21:12
非常感谢您的建议！
我尝试了一下，数据确实平了，但是还想请教两个问题：
1、这样得到的值大约在0.88% ...

第一个问题我没太看懂，
第二个问题简单来说就是当用你的方法统计角度分布的时候，不同角度都乘了个权重值，这个权重正比于该角度下圆环区域的大小（见链接中的图），要把这个权重除掉才能得到真实的角度分布，详细的内容可以参考下面的链接
http://keatonb.github.io/archivers/uniforminclination

作者
Author: chieko 时间: 2023-11-30 16:11
本帖最后由 chieko 于 2023-11-30 16:13 编辑

Lacrimosa 发表于 2023-11-30 12:21
第一个问题我没太看懂，
第二个问题简单来说就是当用你的方法统计角度分布的时候，不同角度都乘了个权重 ...

感谢您提供的思路！这确实是一个看起来很简单，但实际上很容易被忽略的点。
虽然我现在还是没有完全想明白——因为在MD里面计算角度，似乎并不涉及“观察球壳上某点相对于观测点的角度”？这样也就不涉及不同角度下球壳面积不同，导致的取样概率不一样的问题了。HBs与z轴的角度值并不依赖于某个观察点，而是直接用原子坐标计算的。
其实我前两天也做了个小测试，随机在（20，30）的范围内生成50万行，每行6个随机数，每行的6个数代表随机的两个原子坐标。然后计算每对坐标的连线与z轴的角度分布。
我本觉得这够随机了，统计的角度分布应该得到一条直线——但意外的是结果近似sin曲线（还不完全是）。现在看来，可能正是您提到的原因。

另外，上面我第一个问题的意思就是说，因为统计角度时是在0-180°上每间隔1°统计的，因此对于各向同性的体系，这个直线的分布值应该是1/180=0.00556即0.55%。
至于为什么我的结果不是这个值，可能也是归一化的问题吧？这倒是好解决（笑，就是不知道会不会不严谨）。

再次感谢您的解释和帮助！

作者
Author: Lacrimosa 时间: 2023-11-30 17:07

chieko 发表于 2023-11-30 16:11
感谢您提供的思路！这确实是一个看起来很简单，但实际上很容易被忽略的点。
虽然我现在还是没有完全想明 ...

对于第一个问题，举个例子，如果你在程序里给每个氢键形成的向量进行归一化，最终这些氢键向量就组成了一个半径为1的球面了。这个时候统计向量出现的次数和直接数球面上的点是同样的操作。尽管你在程序中转化成了角度，实际上也还是在数点而已，这个时候2*Pi*sinθdθ的权重就自然被考虑进去了。

第二个问题的信息不太够，不好判断呢。如果你有180个点，分布是一条y=0.88的直线，那求和以后不等于1啊。

作者
Author: chieko 时间: 2023-11-30 22:47

Lacrimosa 发表于 2023-11-30 17:07
对于第一个问题，举个例子，如果你在程序里给每个氢键形成的向量进行归一化，最终这些氢键向量就组成了一 ...

我好像有点明白了，我再去自己研究的体系上试试，看结果和预期是否吻合；
至于结果除以sin(theta)后积分不为1的问题，可能暗示这其中还是有点什么地方没有弄对，我暂时也想不出什么原因了，可能只好先“强行”让它归一化了。
非常感谢您的帮助！！

作者
Author: naoki 时间: 2024-9-9 00:18
本帖最后由 naoki 于 2024-9-9 07:41 编辑

Lacrimosa 发表于 2023-11-29 19:05
可能是归一化的问题，把结果除以一个sin(theta)看看

感觉这样的话0°和180°附近会有过大的权重啊

作者
Author: Lacrimosa 时间: 2024-9-9 13:18

naoki 发表于 2024-9-9 00:18
感觉这样的话0°和180°附近会有过大的权重啊

实际操作中好像不会遇到这个问题，毕竟对角度划分bins的时候不可能取到sin(0)

作者
Author: naoki 时间: 2024-9-9 15:36

Lacrimosa 发表于 2024-9-9 14:18
实际操作中好像不会遇到这个问题，毕竟对角度划分bins的时候不可能取到sin(0)

我取到了0.5，然后一柱擎天

作者
Author: Lacrimosa 时间: 2024-9-9 21:44

naoki 发表于 2024-9-9 15:36
我取到了0.5，然后一柱擎天

如果把样本量继续增大说不定能降下来？

作者
Author: yxzinc 时间: 2024-9-10 09:59
bulk water相关角度分布本来就该是你第一张图那样的啊，见Nat Commun 11, 5809 (2020).的fig5黑线，你第二张图才是做了什么数据处理吧

作者
Author: 云非侠 时间: 2024-9-25 16:43
我的观点：向量与Z轴夹角的概率分布就是应该是抛物线而非概率相等的直线。理由如下：我也遇到了这个问题，在盒子里随机放入线性分子，统计线性分子与X/Y/Z轴的夹角分布概率，也是得到了从0-180度，概率密度先增大后减小，在90度附近概率最大的抛物线形状。按2楼的方法将结果除以sin(theta)之后的确得到了之前预想的直线，但是概率加和也是不等于1。于是我又对同一个体系求了与X轴、Y轴的夹角概率分布，发现也是相同的抛物线，查找资料得到这样一句话"当分子与某一轴的夹角为 90 度时，它们的取向是最随机的。此时，分子在两个平面内的运动自由度最大，从而导致这一位置的分布概率最大",我试着在纸上画了三维坐标图，想象一个球体，发现当分子与Z轴的夹角为90度的时候，自由度的确是最大的，分子可以在XY平面内任意分布，当与Z的夹角从90度缩小到60度、30度、0度的时候，自由度逐渐减小，最后只能与Z轴平行，根本就不自由了。可能这样还是很难理解，那就再极端一点，直接比较与Z轴的夹角为0度还是90度，分子分布的概率是不是相同的，在0度时，分子或者说向量只能与Z轴平行，在XY平面没有任何变化的空间，而当夹角为90度时，向量在XY平面既可以选择与Y轴或者X轴平行，也可以自由地在XY轴平面内选择一个方向朝向。由此可见，90度的概率就是应该比0度或者180度大才对。以上是我的一点观点。

作者
Author: chieko 时间: 2024-11-5 21:52
本帖最后由 chieko 于 2024-11-5 22:13 编辑

搞明白了，说到底还是自己数学基础不好hhh；其实关于此早有严格的数学描述，见：
【空间中两随机向量间夹角的概率密度分布】 https://jerkwin.github.io/2013/0 ... %E5%88%86%E5%B8%83/
例如对于3维空间，随机向量夹角θ的分布的概率为：P(θ) = 1/2 sin(θ)
感谢大家的讨论！

欢迎光临计算化学公社 (http://bbs.keinsci.com/)