本帖最后由 Acee 于 2023-4-13 14:01 编辑
GROMACS中伞型采样算法应用详解 原文链接(https://mp.weixin.qq.com/s/PG8niVc2r6aIY9pn_WbprA)
1. 背景介绍 伞型采样作为一种增强采样算法，被广泛应用于分子动力学模拟的各个领域，例如多肽构象牵引、分子穿越细胞膜、囊泡融合、囊泡细胞膜融合、分子二聚自由能计算、分子溶剂化自由能计算、蛋白质-配体解离牵引、丙氨酸二肽 Phi/Psi 旋转等一系列非常有意思的生物物理现象。部分例子如下图所示： 本文将详细讲解伞型采样算法在GROMACS软件中的使用方法，并且提供相应的mdp文件以供大家参考。本文以及公众号之后的相关文章都将以“精准，高效，简约，可靠”为撰写主旨。 2. 理论基础 我们通常认为的无偏分子动力学模拟（unbias molecular dynamics simulation）和物理化学里面描述的孤立体系有所类似，即系统与环境之间无物质交换也无能量的传递。就一般的凝聚态体系来说，倘若确定了体系中的温度，压力和物质的量，则系统的状态就被确定了。常规的无偏动力学模拟由于其本身底层势能函数和模拟时间的限制，对于高能态，过渡态或者稀有构象的采样率很低。因此，伞型采样，副本交换动力学等增强采样方式应运而生。这其实是一种从外部引入能量从而帮助体系越过本身需要花很长时间或者终身无法越过的能垒。其原理简述于下图： 3. 伞型采样在GROMACS中的应用 我们通常所说的伞型采样一般分为三步完成。 第一步是确定反应坐标，通过合理配置mdp文件中的pull code获取牵引构象。 第二步，合理离散反应坐标，确保关键的路径有足够的窗口来描述。 第三步，通过WHAM计算反应路径的平均力势（Potential of Mean Force, PMF）。其中第一步是最关键的一步，下面的内容会对反应坐标的确定和选择进行极为详细的描述。
code如下图所示（以distance作为坐标几何变化类型为例）： - ;Pull code
- pull = yes ; 有无质心牵引
- pull_ngroups = 2 ; 牵引组数目
- pull_ncoords = 1 ; 牵引坐标数目
- pull_group1_name = SOLU ; 牵引组名称，对应index文件里面的名称
- pull_group2_name = POT ;
- pull_coord1_type = umbrella ; 使用参考组与一个或多个组之间的伞型势来牵引质心
- pull_coord1_geometry = distance ; 连接两个组的向量进行牵引（distance/direction-periodic）
- pull_coord1_groups = 1 2 ; 与此牵引坐标作用的组的索引
- pull_coord1_dim = N N Y ; 选择此牵引坐标作用的维度，分别对应X,Y,Z三个分量（N表示关闭，Y表示开启）
- pull_coord1_rate = 0.005 ; 牵引速率（nm/ps,度/ps）
- pull_coord1_k = 1500 ; 弹簧刚度（kJ mol^-1^ nm^-2^）
- pull_coord1_start = yes ; 将初始构型的质心距离加到pull-coord1-init
很多GROMACS使用者都苦于不知如何设置pull code来完成理想牵引过程。在此，我将对该code里面的几个关键点进行介绍。以便于之后的研究者一上来就知道如何设置自己的体系的code。 第一点，geometry的设定。GROMACS中定义了八种坐标反应的类型，此处我们对常用的前三种反应类型的优缺点和适用情况进行讨论。分别是distance、direction和direction-periodic。 特点：
distance是通过距离增加从而达到牵引的目的。注意，这里是距离增加是指，两个反应组的距离在distance条件下的距离是不断增加的，直到模拟结束或者牵引的距离超过盒子牵引方向所在轴长度的一半（在想考察的牵引距离内，倘若沿着Z轴运动，则需要保证盒子Z轴的长度要大于牵引距离的二分之一，其他轴亦然）。这种反应类型的用处特别广泛，通过合理的设置牵引对象，可以用于多肽的牵引、小分子过膜、分子二聚结合自由能计算，以及囊泡融合等体系，是伞型采样中的万金油方法。 direction是沿着固定方向（pull-coord-vec）进行牵引（如x，y，z），直到模拟结束或者牵引的距离超过盒子牵引方向所在轴长度的一半。这种情况同样可以用于多肽的牵引，小分子过膜，以及囊泡融合等过程。也是万金油方法。 direction-periodic和direction类似，但不使用周期性边界条件对距离进行校正， 允许距离超过盒子长度的一半。(只有)通过使用牵引速率连续改变参考位置，将两个组推离开的距离超过盒子长度一半时，才需要这个选项。
缺点： distance无法通过pull-coord-vec来控制牵引方向，故牵引方向仅取决于两个组距离增加的方向。经常有人问，为什么自己用distance的时候，分子的运动方向不是朝我想要的方向运行，那么我想现在你应该明白了。此外，需要保证牵引距离小于盒子牵引方向所在轴长度的二分之一 direction无明显缺点，方向可控，速率可控。但同样需要保证牵引距离小于盒子牵引方向所在轴长度的二分之一。 direction-periodic无法使用NPT，即在模拟过程中盒子的大小不能发生变化。
4. 三种geometry在分子穿越细胞膜体系中的应用与讨论 （1）distance类型在分子穿越细胞膜体系中的应用和讨论： 实际上在常规的牵引体系中，我们通常将被牵引的分子设为一个组（SOLU），再将细胞膜设为一个组（MEMB）。但根据distance原理，我们这样设置后，SOLU和MEMB的距离会逐渐增加，这样就不可能实现穿越细胞膜的过程（当然也也许有的同学会说，我可以把牵引速率设为负值（pull_coord1_rate为负），这样SOLU不就可以向反方向运动了吗）。很不幸，这样的操作确实可以让SOLU短暂的向负方向运动从而靠近细胞膜，但是SOLU和MEMB的距离只会无限趋近于0。所以，SOLU是不可能穿过细胞膜的。想要使用distance作为geometry来实现分子穿越细胞膜的过程，我们主要注意两个关键点（并列关系）。第一种方法，将体系内原本的原子作为冻结组或者位置限制组（例如一个水分子），例如选择一个与膜的距离比被牵引组离膜的距离稍微远一点（1nm以内都可）的分子冻结或者限制起来，将该分子替代MEMB组。第二种方法，引入虚拟粒子作为冻结组或者位置限制组，将虚拟粒子放置在与膜的距离比被牵引组离膜的距离稍微远的位置。这两种方法引入的组我们暂且称为POT组。如下图所示： 通过这样的方式，随着模拟的运行，SOLU和POT的距离逐渐增加，直到完全穿过细胞膜。
（2）direction类型在分子穿越细胞膜体系中的应用和讨论：
这种情况就要比distance的类型简单多了，不需要引入冻结组，限制组或者虚原子。此处我们将被牵引分子称为SOLU，与MEMB连接成为反应坐标。确定好SOLU离膜的初始距离和牵引方向（pull-coord1-vec）后即可开始运行。如下图所示：
（3）direction-periodic类型在分子穿越细胞膜体系中的应用和讨论： 其设置思路和direction类似，但是盒子可以不用设那么大，但需要确保increasing distance要比Mirror distance 1 和2的和要小。原因是这种情况下GROMACS会自动计算两个索引组的质心距离，并按照最短的距离进行牵引。一般来说不需要太过担心，因为当前SOLU和MEMB的质心距离只有increasing distance的一半，是相当安全的。如下图所示：
通过几点的详细论述，我相信大家对GROMACS伞型采样的应用有了比较清晰的理解。 5. 三种geometry的pull code示例（分子穿膜） （1）distance（此处注意要对pull_group2_name施加一个冻结效果，或者限制效果，确保虚原子在整个模拟过程中位置不发生变化。个人推荐使用位置限制）
- ;Pull code
- pull = yes ; 有无质心牵引
- pull_ngroups = 2 ; 牵引组数目
- pull_ncoords = 1 ; 牵引坐标数目
- pull_group1_name = SOLU ; 牵引组名称，对应index文件里面的名称
- pull_group2_name = POT ;
- pull_coord1_type = umbrella ; 使用参考组与一个或多个组之间的伞型势来牵引质心
- pull_coord1_geometry = distance ; 连接两个组的向量进行牵引（distance/direction-periodic）
- pull_coord1_groups = 1 2 ; 与此牵引坐标作用的组的索引
- pull_coord1_dim = N N Y ; 选择此牵引坐标作用的维度，分别对应X,Y,Z三个分量（N表示关闭，Y表示开启）
- pull_coord1_rate = 0.005 ; 牵引速率（nm/ps,度/ps）
- pull_coord1_k = 1500 ; 弹簧刚度（kJ mol^-1^ nm^-2^）
- pull_coord1_start = yes ; 将初始构型的质心距离加到pull-coord1-init
（2）direction
- ;Pull code
- pull = yes
- pull_ngroups = 2
- pull_ncoords = 1
- pull_group1_name = SOLU
- pull_group2_name = MEMB ;从虚原子变为膜
- pull_coord1_type = umbrella
- pull_coord1_geometry = direction ; 从distance变为direction
- pull_coord1_groups = 1 2
- pull_coord1_dim = N N Y
- pull_coord1_rate = 0.005
- pull_coord1_k = 1500
- pull-coord1-vec = 0.0 0.0 1.0 ; 调控牵引组运动方向
- pull_coord1_start = yes
（3）direction-periodic
- ;Pull code
- pull = yes
- pull_ngroups = 2
- pull_ncoords = 1
- pull_group1_name = SOLU
- pull_group2_name = MEMB
- pull-group1-pbcatom = xxxxx ; 处理PBC时参考原子的编号
- pull-group2-pbcatom = xxxxx ;
- pull-pbc-ref-prev-step-com = yes ; 使用上一步的质心作为参考来处理周期性边界条件
- pull_coord1_type = umbrella
- pull_coord1_geometry = direction-periodic
- pull_coord1_groups = 1 2
- pull_coord1_dim = N N Y
- pull_coord1_rate = 0.005
- pull_coord1_k = 1500
- pull-coord1-vec = 0.0 0.0 1.0
- pull_coord1_start = yes
6. 几点注意 盒子大小：注意在使用distance和direction的时候，盒子在牵引轴的大小需要比牵引距离的一半大。 NPT/NVT选择：direction-periodic只能使用NVT。其他情况根据体系特点选择。但是要明白，不是什么体系都能无脑使用NPT。 弹簧刚度k：k就是弹簧系数，值越大，弹簧刚度越强。则引入偏置势的影响越大。一般小体系来说，1000-1500足够，如果是要牵引囊泡融合这样的体系，则需要更大的弹簧系数，具体值视体系大小而定。 伞型采样极为灵活，通常要实现理想的采样过程，需要多加尝试，反复调整。
后记 第一次写这样的文章，颇费心力。但由于笔者本身水平有限，细节地方还需要读者多加推敲。本期内容到此结束。下期会继续讲解伞型采样中后面两步的计算过程。另，欢迎加入我的交流群进行讨论。谢谢大家。 qq群：477872340；微信群需要暂时先加鄙人微信，我将大家拉入群聊，微信号：zys1220098477 参考 GROMACS official tutorial In-house data In-house data Ou, L., Chen, H., Yuan, B., & Yang, K. (2022). Membrane-Specific Binding of 4 nm Lipid Nanoparticles Mediated by an Entropy-Driven Interaction Mechanism. ACS nano, 16(11), 18090-18100. Amber official tutorial
致谢 特别感谢He Furui，Huang Ming等人对本文的帮助和支持。
| 
[复制链接 Copy URL]
发表于 Post on 2023-4-13 13:54:16
|
只看该作者 Only view this author
|只看大图 View big picture
|正序浏览
|阅读模式 Reading model
收藏 Add to favorites
楼主 Author