利用软势修正低精度大型蛋白模型的原子重叠问题一例

Graphite

利用软势修正低精度大型蛋白模型的原子重叠问题一例

石墨 首发计算化学公社

    在预处理一些大型、柔性模型时，会遇到表现为”bad contact”、”inf(infinity) energy / force”、”can not be settled”、”step XXX.pdb 刷屏”这些情况，往往是由结构不合理或力场不合适，导致瞬间的力过高、体系崩溃。如果是一些小意外，比如盒子留的边再大1 A就好了、初始建模时自己粗心了，那倒还好。

    然而，有时候，在做蛋白、聚合物、某些大型体系时，我们拿到一些结构复杂的低质量模型，其中可能有很多处不合理的接触，这时候可以考虑一些特别的技巧，包括利用软势。

    这里以某个球蛋白为例，这个蛋白拿到手就是一个建模精度很差的结构，pdb文件里不含氢、也没有经过优化或者复查。它有约1000个氨基酸，加上氢之后有约3.3万个原子。它的二、三级结构没大问题，但细节上精度很不好。

    当使用gmx pdb2gmx和一些MD相关的可视化软件、工具包加氢后，运行正常的能量最小化，很快报错，提示13639号原子具有“无限大”的受力崩溃。

    打开VMD，选取"same residue as (within 3 of index 13639)"以显示周围的情况。果然，这里的两个苯丙氨酸残基有问题。注意画红圈的部位，很明显H原子重叠了，而且上方的N原子也离苯环太近，VMD甚至默认以为C-N之间成键了，而这两个残基并不是直接相连的，蛋白质里面也不可能存在这种乱连的情况。这固然一部分原因是gromacs并没有考虑这种情况，只是简单的按照一般键长、键角加了氢，但本质原因还是这个结构本身精度不好。

    如果是小分子，也许用gaussview还是chemcraft，稍微拉开它们，再做优化也就罢了。然而这个蛋白质有3万多个原子，里面这些“打架”的情况不可能人工处理。如果编程处理，比如将间距<1 A的原子在附近重新指定一个还算合理的位置，也并非不可，但要做到以后也能用，设计算法、调试，恐怕一两天是不够的。

    不过，即使VMD以为他们胡乱成键了，gmx pdb2gmx给出的top文件是严格按照序列来的。键、键角、二面角都严格按照残基指定，因此有了拓扑信息，对于MD程序来说，两个苯丙氨酸内部的键、键角、二面角互不干扰，只是LJ和静电部分因为接触太近而给出了过高的值，使得瞬间力和速度过大。

    这种时候不妨将“LJ和静电”换掉，也就是换一个pair_style。在LAMMPS中，可以使用pair_style soft。

    LJ势和用于计算静电的Columb势都是距离越近，受力急剧升高的函数。如果原子间距离仅有sigma的1/2，那么它的LJ相互作用强度接近平衡距离时强度的4000倍（2^12-2^6）。可以利用软势（soft）暂且代替它，“温和”地推开距离太近的原子。

    具体的做法是：

    1、将已经生成的gro文件和top文件输入Intermol，这是一个转换不同MD软件间格式的软件。生成Lammps适用的in文件和data文件，Intermol将它们命名为*_converted.input和*_convert.lmp。我们将其重命名为soft_min.in和soft_min.data。in文件包含模拟设置；data文件包含原子坐标和连接性，力场参数可以写在in文件里，也可以写在data文件里。

    2、将in文件里以下内容删除：

1. special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 0.8333
   
2. 
   
3. pair_style lj/cut/coul/long 9.0 9.0
   
4. pair_modify tail yes
   
5. kspace_style pppm 1e-8
   
6. 
   
7. pair_coeff 1 1   0.1700000   3.2500000
   
8. pair_coeff 2 2   0.0157000   1.0690800

复制代码

    这表示用经典的LJ/Coulumb作用描述原子对间的作用，并且设定了一些如cut-off、长程如何处理的细节。
    将他们替换为：

1. pair_style soft 1.0
   
2. pair_coeff * * 1.0

复制代码

    也就是采用soft势描述原子对间的作用，强度为1.0、参考距离和cut-off为1.0。考虑soft势的函数形式，这意味着小于1.0 A距离的原子间会受到一个不太大的斥力，而大于1.0 A的不受影响。
    在文件的最后再将run 0删掉，替换成以下内容。这表示进行能量最小化，并输出轨迹文件。

1. dump 1 all atom soft_min.lammpstrj   
   
2. minimize 1e-6 1e-6 10000 10000

复制代码

    最终的文件：

1. units real
   
2. atom_style full
   
3. boundary p p p
   
4. 
   
5. bond_style hybrid harmonic
   
6. angle_style hybrid harmonic
   
7. dihedral_style hybrid multi/harmonic
   
8. 
   
9. read_data soft_min.data
   
10. 
    
11. pair_style soft 1.0
    
12. pair_coeff * * 1.0
    
13. 
    
14. thermo 50
    
15. thermo_style custom ebond eangle edihed eimp epair pe
    
16. dump 1 all atom 200 soft_min.lammpstrj
    
17. 
    
18. min_style cg
    
19. minimize 1e-6 1e-6 10000 10000
    
20. 
    
21. write_data soft_min.data

复制代码

    单核执行lammps，不到一分钟即完成。这个蛋白本身的键长、键角等是没有多少阻力的，而体系的总势能却直接下降了3个数量级，可见初始结构的原子重叠问题非常严重。

    最后，原先位置的两个苯丙氨酸残基的结构变成了这样：

    而体系总体的二、三级结构不变，其余部分也得到了修正。最后，再利用VMD将LAMMPS的轨迹转为pdb文件，即可以为gromacs正常使用。

几个可能会关注的问题：

1、如果只是因为LJ和Coulumb在近距离过高，在gmx里将它们的参数（sigma、epsilon、电荷）降低、或直接关闭是否可行？

不可行，这不会改变函数在近距离直线上升的形式。另外，如果直接关闭，那么像残基外围的氢之间就完全没有了斥力，仍然会意外重叠。

2、用软势推开重叠原子之后，还需要换回经典LJ和Coulumb能量最小化吗？

需要，之后就正常做模拟，这一步只相当于代替我们手动精修初始结构。最终还是在真正要用的分子力场中跑。

3、gmx也有soft-core函数（本意用来做自由能计算的），是否也能用？

实测对于重叠比较严重、偏离比较明显的结构，在修正初始结构的情景下，其形式和效果都不如LAMMPS包含的这个cosine函数好用。

4、能用软势做正常模拟吗？它还有什么作用

它没有明确的物理意义，不能够用来做正常模拟。还可以用来优化缠绕的聚合物、将物质“蓬松地”填到空间里（巧妙设定强度和cutoff）、用来操纵结构制作一些特定的示意图等。

5、lammps和gromacs之间结构文件互相转化还有哪些办法常用？

lammpstrj文件本质上就和xyz文件差不多，只含有类型、坐标、电荷，很多软件都能直接读写。如果要包含连接性信息，则需要读写data文件。VMD可以通过topo readlammpsdata和topo writelammpsdata指令读写；OVITO可以读取其他MD程序的结构，并制作data文件；lammps本身也可以输出pdb和xtc等类型的文件。