仓库源文，站点原文

layout: post title: GROMACS多糖模拟简单示例 categories:

• 科 tags:
• gmx
• amber chem: true

• 2017年04月10日 20:42:32

以前, 使用MD模拟糖分子并不容易，因为针对蛋白质和核酸设计出来的力场用在糖分子上不大合适, 因此需要对常用力场进行修正. 可参考下列文献(文献稍老, 请查阅更新的)

• CHARMM: Carbohydrate Solution Simulations: Producing A Force Field With Experimentally Consistent Primary Alcohol Rotational Frequencies And Populations. Michelle Kuttel, J. W. Brady, Kevin J. Naidoo; J. Comput. Chem. 23(13):1236-1243, 2002; 10.1002/jcc.10119

• OPLSAA: An Improved Opls-aa Force Field For Carbohydrates. D. Kony, W. Damm, S. Stoll, W. F. Van Gunsteren; J. Comput. Chem. 23(15):1416-1429, 2002; 10.1002/jcc.10139

• GROMOS: A New Gromos Force Field For Hexopyranose-based Carbohydrates. Roberto D. Lins, Philippe H. Hünenberger; J. Comput. Chem. 26(13):1400-1412, 2005; 10.1002/jcc.20275

目前, 使用比较广泛的多糖力场是GLYCAM:

GLYCAM06: A Generalizable Biomolecular Force Field. Carbohydrates

Karl N. Kirschner, Austin B. Yongye, Sarah M. Tschampel, Jorge González-outeiriño, Charlisa R. Daniels, B. Lachele Foley, Robert J. Woods; J. Comput. Chem. 29(4):622-655, 2007; 10.1002/jcc.20820

这里简单示例下如何利用AmberTools和acpype创建多糖的GROMACS输入文件并进行模拟.

所需程序

• AmberTools: 用于生成AMBER的输入文件. 可以下载编译安装官方版本, 经验不足的话, 比较麻烦. 建议使用我整合的Windows版, 参考 Windows下的AmberTools+RESP+ACPYPE.

• acpype脚本: 用于将AMBER输入文件转换为GROMACS输入文件. 可以在网上下载, 但要注意版本. 我整合的Windows版AmberTools中已经自带了acpype脚本及其可执行程序, 无须另外安装. 使用方法可参考 使用AmberTools+ACPYPE+Gaussian创建小分子GAFF力场的拓扑文件.

注意, 下面的示例基于我整合的Windows版AmberTools, 如果你使用其他版本, 可能需要对命令进行相应修改.

具体步骤

第一步: 了解基础知识

了解多糖的基本信息, 阅读AmberTools手册的相关内容.

相关资料可参考多糖模拟系列博文:

• 糖及其GLYCAM力场中的命名约定
• Amber构建寡糖的方法
• Amber构建糖蛋白的方法

第二步: 生成AMBER输入文件

以手册中的示例来熟悉多糖的构建方法.

新建文件gly

<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span style="color: #008800; font-style: italic"># 加载参数</span> <span style="color: #AA22FF">source</span> leaprc.GLYCAM_06

<span style="color: #008800; font-style: italic"># 定义多糖序列</span> <span style="color: #B8860B">gly</span> <span style="color: #666666">=</span> sequence<span style="color: #666666">{</span> ROH 4YB 4YB 3MB 0MA <span style="color: #666666">}</span>

<span style="color: #008800; font-style: italic"># 设定 4YB-4YB, 3MB-4YB 之间的扭转角 psi</span> impose gly <span style="color: #666666">{3</span> 2<span style="color: #666666">}</span> <span style="color: #666666">{</span> <span style="color: #666666">{</span>C1 O4 C4 H4 0.0<span style="color: #666666">}</span> <span style="color: #666666">}</span> impose gly <span style="color: #666666">{4</span> 3<span style="color: #666666">}</span> <span style="color: #666666">{</span> <span style="color: #666666">{</span>C1 O4 C4 H4 0.0<span style="color: #666666">}</span> <span style="color: #666666">}</span>

<span style="color: #008800; font-style: italic"># 保存 top, crd, mol2文件</span> saveamberparm gly gly.top gly.crd savepdb gly gly.pdb savemol2 gly gly.mol2

<span style="color: #008800; font-style: italic"># 退出</span> quit </pre></div> 这里建议将构型同时保存为pdb格式和mol2格式. 与pdb格式相比, mol2格式保存了更多与力场相关的信息, 如原子类型, 电荷等, 查看修改更方便, 也可以直接用于AmberTools的处理. 使用mol2格式的缺点在于有些可视化程序可能无法识别其中的原子类别, 导致显示异常(但VMD显示正常). 这种情况下使用pdb格式查看就可以了.

使用sleap命令来生成AMBER输入文件

sleap -f gly.in

这样我们就得到了多糖的Amber输入文件: 拓扑gly.top, 坐标gly.crd, 构型gly.mol2.

其构型如下

第三步: 转换为GROMACS输入文件

使用acpype将AMBER输入文件转化为GROMACS输入文件

acpype -p gly.top -x gly.crd -d

会得到GROMACS坐标文件gly_GMX.gro和拓扑文件gly_GMX.top. 此外acpype还会自动给出两个mdp文件: em.mdp和md.mdp, 但基本没什么用, 因为里面的设置太少.

第四步: 运行GROMACS模拟

得到坐标文件和拓扑文件后, 再准备好mdp文件和索引文件，遵照通常的流程进行模拟就可以了.

由于GLYCAM和AMBER力场属于一个系列, 水模型都推荐用TIP3P.

下面是gly的一段模拟轨迹

存在的问题

AMBER力场中, 分子内1-4非键相互作用的校正因子为0.5和0.8333, 而GLYCAM力场中二者都是1.0. 所以, 如果体系中只有多糖的话, 需要将上面第三步所得拓扑文件中[ defaults ]下面的fudgeLJ和fudgeQQ都改成1.0.

对于同时使用AMBER和GLYCAM力场的体系, 1-4非键作用的校正因子如何设置呢? 最简单的方法是不修改, 直接使用AMBER力场的设置0.5和0.8333. 这样的话GLYCAM的准确度可能会降低些, 但如果你更关心蛋白质的构型, 这种方法也未尝不可. 实际上, Amber和GROMACS都可以实现使用混合的1-4非键作用因子. 在GROMACS中只要把[ pairs ]部分的函数类型改成2, 并提供fudgeQQ, q_i, q_j, sigma_ij, epsilon_ij五个参数就可以了. 详细说明请参考GROMACS手册 5.3.4 分子内的对相互作用, 表 5.5: moleculetype 指令详解.

扩展应用

如果你理解了上面的作法, 请思考并解决下面两个问题:

• 如何构建环状多糖, 或周期性多糖?
• 如果多糖结构文件已知, 如何修改上面的流程获得其GROMACS输入文件?

评论

• 2017-04-17 15:37:30 Teddy 老师你好，如果多糖结构文件已知该怎么弄啊？

• 2017-04-17 22:34:17 Jerkwin 这是个作业, 你学了上面的, 就应该会做了.

• 2017-05-15 13:11:03 Teddy 不好意思，李老师，我学了上面的，还是弄不好，我是从英国剑桥晶体数据库下载的纤维素的单晶，然后转化成mol2文件，直接写

    source leaprc.GLYCAM_06
    cell = loadmol2 cell.mol2
    saveamberparm cell cell.top cell.crd
    savepdb  cell cell.pdb
    savemol2 cell cell.mol2
    quit
  

  出现很多错误 额。。应该不是这样弄吧？ 请李老师教一下，谢谢