计算模拟轨迹中水的方位四面体序参数Orientational Tetrahedral Order的脚本

对抗路达摩

水的四面体序参数（Orientational Tetrahedral Order）是一种用来表征水的结构特征的参数，其定义如图所示:



更多有关四面体序参数的资料参考：
Relationship between structural order and the anomalies of liquid water | Nature
5 Minutes Tutorial — Order 0.0.1 documentation


在分子模拟中，我们可以计算模拟轨迹中每个水分子在每一个时刻的序参数，进一步得到序参数的平均值、分布、随时间的变化、随环境的变化来帮助我们研究各种过程中水的结构性质。
那么这里我们使用一个python脚本来实现。代码主要基于MDAnalysis库，有基础的python和MDAnalysis使用经验的朋友可以轻松地修改代码应用到你的体系。

在这个脚本及show-case中，我们实现了模拟轨迹内距离蛋白质残基GLN16附近5A以内的水分子的序参数的计算，并绘制了序参数的分布。
根据定义，水的OTO由其自己的氧原子和离其最近的四个水分子的氧原子的相对位置关系决定。所以代码需要实现的两个事情：1.找到离每个水分子最近的4个水分子。2.根据这五个水分子的坐标计算序参数。

程序中最重要的函数是cal_order_frame，这个函数有两个输入waters和water_buffer。两个都是MDAnalysis的atomgroup对象。其中waters就是所有你需要计算序参数的水分子的氧原子，在这里就
是name OW and around 5 (resid 15), 也就是GLN16(因为蛋白从GLY2开始，所以GLN16的resid是15）周围5A所有名字叫做OW的原子。而water_buffer这里则定义为'name OW and around 8.5 (resid 15)'，这是因为离waters里面一个water最近的四个water，不一定都在waters里面。所以设置了3.5A的缓冲区域，对于每一个waters中的water，我们将会在water_buffer里面寻找离他最近的氧原子。实际上是为了计算效率的一种考虑(不用遍历所有的水)。这里3.5A是一个比较安全比较大的值，理论上可以设计的更小，但我因为不是研究水的，在这里为了正确性，取了一个比较保守的值。
有了这两个atomgroups之后，cal_order_frame就会遍历waters中的每一个water，并在water_buffer中寻找离它最近的四个氧原子。之后会将这个氧和离它最近的四个氧传入函数cal_order计算该水分子的OTO，最后我们会将每一个水的序参数存在一个list里面并返回。那么我们只需要遍历轨迹，每一帧调用cal_order_frame，并把每一帧返回的list加起来，就可以得到所有蛋白质残基GLN16附近5A以内的水分子的序参数的值，然后就可以根据此画出序参数的分布图。需要注意的是在这个过程中选择的atomgroup要使用动态选取，即updating=True，来保证每一帧选的都是GLN16附近的水。

该脚本最会输出两个文件，一个是用matplotlib画的分布图Orientational Tetrahedral Order.png，一个是记录了每一个OTO的值的result.txt。分布图如下所示：（使用我提供的文件绘制的结果会更不平滑很多，因为我为了附件大小skip了很多frame）




代码的缺陷：在这个代码中我们没有考虑periodic boundary conditions（PBC），在我们的show case中，由于蛋白质的平动和转动都被消除了，且残基离盒子的边界远远大于8.5A，所以不会带来任何的问题。在实际使用中，如果你需要考虑的水分子在盒子边界旁边且你使用了PBC，那么你需要修改代码来显式地考虑PBC。主要要修改两个代码：1.寻找离水分子最近的4个水分子的部分的代码。 2.根据三个点的坐标计算cos的代码。

最后一个有意思的点：在这个计算残基附近水的例子中，代码的计算复杂度随着盒子大小是O（1）的，也就是常数复杂度（随着轨迹frame数自然是O（n））。因为每一帧的计算量只和5A小球中水的个数n以及8.5A大球中水的个数m有关，而m和n显然与盒子大小。但这个代码又确实说不上快，这就告诉我们常数复杂度和快不快也没必然联系。

牧生

【水分子F3/F4序参数计算编程-哔哩哔哩】 https://b23.tv/3v8sGTW

看一下这个视频呢

对抗路达摩

牧生 发表于 2023-9-2 21:26
【水分子F3/F4序参数计算编程-哔哩哔哩】 https://b23.tv/3v8sGTW

看一下这个视频呢

“序参数”有很多种，这里的F3F4和我这里用的OPO形式不一样，描述的事情也有区别。本质上可以看做一个映射，是高维的多体空间坐标投影到一维坐标轴上；这样的映射有很多种，每种的值有不同的含义。对一个序参数来说，比较重要的就是这个值能不能比较好地来反应水的结构，在水的不同结构下有明显的特征取值。

Lacrimosa

目前好像确实没什么现成的程序能算这个，大多都是自己写程序分析

王五

作者您好，我看MDAnalysis官网里说，选择原子方式“prop”是默认考虑周期性边界的，这样是不是就不需要修改您所提到的两块内容，如果是其他要考虑周期性边界的选择方式的话，是不是将periodic=True设置在select_atoms就可以了，好像也不用修改角度和水分子选择的代码？

对抗路达摩

王五 发表于 2023-12-28 14:51
作者您好，我看MDAnalysis官网里说，选择原子方式“prop”是默认考虑周期性边界的，这样是不是就不需要修改 ...

很有可能！
可以找一个是否正确考虑PBC有影响的例子，然后测试一下periodic=True和False的结果有没有区别

chieko

按照原文描述（https://doi.org/10.1038/35053024），这个参数q应该介于0到1之间呀，会不会哪里有什么问题

——————————————————————————————
刚搞明白，是<q>应该在0到1之间