哲●科●文Jerkwin 类别 标签 链接 关于 订阅

GROMACS使用amber19sb力场

类别: 科   标签: gmx amber   阅读次数:   版权: (CC) BY-NC-SA ✕

amber力场一直在更新, 目前最新的是amber19sb, 但gmx自带的amber力场却一直没有更新, 至今还是amber99sb. 虽然网上也有gmx格式的amber14sb力场(存在一个小问题, 但容易解决, 见amber14SB力场报错“No default Improper Dih. types”的解决办法), 但gmx格式的amber19sb力场我暂时还没有看到. 如果要在gmx中使用这个力场, 我觉得可以先使用amber工具生成amber拓扑, 再使用acpype之类的工具将其转换为gmx格式的拓扑. 这种解决方法能否真正走通, 我没有测试, 无法作答. 但即便可行, 其中也还有一个很难解决的问题, 就是amber19sb中的蛋白残基使用了cmap二面角, 而且还基于残基名称匹配, 这违背了传统力场的常规作法, 导致很难兼容. 网上对此也有不少讨论, 可供参考.

既然如此, 那就干脆放下那些投机取巧, 兜兜转转的作法吧, 暴虎冯河, 从根本上解决. 本质问题无非有两个:

  1. 得到gmx格式的amber19sb力场;
  2. 让gmx支持amber19sb力场的cmap.

前一个问题容易解决, 无非是个格式转换问题; 但后一个问题, 在不修改gmx源码的情况下, 恕我愚钝, 能想到的可行方法是修改tpr文件. 这种作法很困难, 也很麻烦, 不易掌握, 通用性也不好. 所以我觉得虽然修改gmx源码也很困难, 但更值得一试, 只要处理好了, 一劳永逸. 代价就是这只能是个特制的gmx版本. 当然, 如果gmx官方版本支持这种作法并更新到新版本中, 那是最理想的.

鉴于此, 我就自告奋勇, 试着创建一个gmx格式的amber19sb力场, 同时修改gmx源码以支持新型的cmap二面角.

于个人而言, 目的如下:

于他人而言, 主要益处则是提供了一个gmx使用amber19sb力场的简单方法, 而无需学习amber工具的使用(得承认, amber工具我用得也不是很熟).

创建gmx格式的amber19sb力场文件

我们这里只转换蛋白的amber19sb力场文件. 为此, 先看看amber的原始力场文件.

如果安装了amber20, 打开E:\amber20\dat\leap\cmd\leaprc.protein.ff19SB(Windows下的路径, Linux下的类比可知), 这是amber19sb力场的调用文件, 内容如下:

leaprc.protein.ff19SB
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
logFile leap.log
#
# ----- leaprc for loading the ff19SB force field
# ----- NOTE: this is designed for PDB format 3!  用于 PDB 3 格式?
#
#	load atom type hybridizations  原子类型
addAtomTypes {
	{ "H"   "H" "sp3" }
	{ "HO"  "H" "sp3" }
	{ "HS"  "H" "sp3" }
	{ "H1"  "H" "sp3" }
	{ "H2"  "H" "sp3" }
	{ "H3"  "H" "sp3" }
	{ "H4"  "H" "sp3" }
	{ "H5"  "H" "sp3" }
	{ "HW"  "H" "sp3" }
	{ "HC"  "H" "sp3" }
	{ "HA"  "H" "sp3" }
	{ "HP"  "H" "sp3" }
	{ "HZ"  "H" "sp3" }
	{ "OH"  "O" "sp3" }
	{ "OS"  "O" "sp3" }
	{ "O"   "O" "sp2" }
	{ "O2"  "O" "sp2" }
	{ "OP"  "O" "sp2" }
	{ "OW"  "O" "sp3" }
	{ "CT"  "C" "sp3" }
	{ "CX"  "C" "sp3" }
	{ "XC"  "C" "sp3" }
	{ "C8"  "C" "sp3" }
	{ "2C"  "C" "sp3" }
	{ "3C"  "C" "sp3" }
	{ "CH"  "C" "sp3" }
	{ "CS"  "C" "sp2" }
	{ "C"   "C" "sp2" }
	{ "CO"  "C" "sp2" }
	{ "C*"  "C" "sp2" }
	{ "CA"  "C" "sp2" }
	{ "CB"  "C" "sp2" }
	{ "CC"  "C" "sp2" }
	{ "CN"  "C" "sp2" }
	{ "CM"  "C" "sp2" }
	{ "CK"  "C" "sp2" }
	{ "CQ"  "C" "sp2" }
	{ "CD"  "C" "sp2" }
	{ "C5"  "C" "sp2" }
	{ "C4"  "C" "sp2" }
	{ "CP"  "C" "sp2" }
	{ "CI"  "C" "sp3" }
	{ "CJ"  "C" "sp2" }
	{ "CW"  "C" "sp2" }
	{ "CV"  "C" "sp2" }
	{ "CR"  "C" "sp2" }
	{ "CA"  "C" "sp2" }
	{ "CY"  "C" "sp2" }
	{ "C0"  "Ca" "sp3" }
	{ "MG"  "Mg" "sp3" }
	{ "N"   "N" "sp2" }
	{ "NA"  "N" "sp2" }
	{ "N2"  "N" "sp2" }
	{ "N*"  "N" "sp2" }
	{ "NP"  "N" "sp2" }
	{ "NQ"  "N" "sp2" }
	{ "NB"  "N" "sp2" }
	{ "NC"  "N" "sp2" }
	{ "NT"  "N" "sp3" }
	{ "NY"  "N" "sp2" }
	{ "N3"  "N" "sp3" }
	{ "S"   "S" "sp3" }
	{ "SH"  "S" "sp3" }
	{ "P"   "P" "sp3" }
	{ "LP"  ""  "sp3" }
	{ "EP"  ""  "sp3" }
	{ "F"   "F" "sp3" }
	{ "Cl"  "Cl" "sp3" }
	{ "Br"  "Br" "sp3" }
	{ "I"   "I"  "sp3" }
}
#
#	Load the main parameter set.  主参数集
#
set default cmap on                         # 启用了cmap, 转换时需要特殊处理
parm19 = loadamberparams parm19.dat         # 基本参数
frcmod19SB = loadamberparams frcmod.ff19SB  # 修正参数, 覆盖基本参数
#
#	Load main chain and terminating amino acid libraries
#
loadOff amino19.lib                         # 主链残基数据库
loadOff aminoct12.lib                       # C端残基
loadOff aminont12.lib                       # N端残基

#
#	Define the PDB name map for the amino acids
#
addPdbResMap {
  { 0 "HYP" "NHYP" } { 1 "HYP" "CHYP" }
  { 0 "ALA" "NALA" } { 1 "ALA" "CALA" }
  { 0 "ARG" "NARG" } { 1 "ARG" "CARG" }
  { 0 "ASN" "NASN" } { 1 "ASN" "CASN" }
  { 0 "ASP" "NASP" } { 1 "ASP" "CASP" }
  { 0 "CYS" "NCYS" } { 1 "CYS" "CCYS" }
  { 0 "CYX" "NCYX" } { 1 "CYX" "CCYX" }
  { 0 "GLN" "NGLN" } { 1 "GLN" "CGLN" }
  { 0 "GLU" "NGLU" } { 1 "GLU" "CGLU" }
  { 0 "GLY" "NGLY" } { 1 "GLY" "CGLY" }
  { 0 "HID" "NHID" } { 1 "HID" "CHID" }
  { 0 "HIE" "NHIE" } { 1 "HIE" "CHIE" }
  { 0 "HIP" "NHIP" } { 1 "HIP" "CHIP" }
  { 0 "ILE" "NILE" } { 1 "ILE" "CILE" }
  { 0 "LEU" "NLEU" } { 1 "LEU" "CLEU" }
  { 0 "LYS" "NLYS" } { 1 "LYS" "CLYS" }
  { 0 "MET" "NMET" } { 1 "MET" "CMET" }
  { 0 "PHE" "NPHE" } { 1 "PHE" "CPHE" }
  { 0 "PRO" "NPRO" } { 1 "PRO" "CPRO" }
  { 0 "SER" "NSER" } { 1 "SER" "CSER" }
  { 0 "THR" "NTHR" } { 1 "THR" "CTHR" }
  { 0 "TRP" "NTRP" } { 1 "TRP" "CTRP" }
  { 0 "TYR" "NTYR" } { 1 "TYR" "CTYR" }
  { 0 "VAL" "NVAL" } { 1 "VAL" "CVAL" }
  { 0 "HIS" "NHIS" } { 1 "HIS" "CHIS" }
}

#
# assume that most often proteins use HIE
#
NHIS = NHIE
HIS = HIE
CHIS = CHIE

力场参数文件

可以看到, 19sb中用于蛋白的原子类型并不多, 与99sb相比, 新增了几个原子类型, 如CX, XC, C8, 2C, 3C, CO等.

根据这个文件中的调用说明, 我们很容易找到

这两个文件中给出了原子类型, 成键, 非键, cmap参数, 需要将它们转换成相应的gmx的格式, 放到

转换时注意单位, 公式系数, 二面角换算, cmap处理. 其中的二面角处理尤其需要仔细参考amber手册中的相关说明. 还有就是, 相同成键项存在多套参数时, 优先使用frcmod.ff19SB中的参数.

残基数据库

接下来, 处理氨基酸残基. 以主链残基为例吧, 末端残基处理方式相同. 根据调用说明, 查看E:\amber20\dat\leap\lib\amino19.lib:

amino19.lib
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
!!index array str 本文件中的所有残基, 28个
 "ALA"
 "ARG"
 "ASH"
 "ASN"
 "ASP"
 "CYM"
 "CYS"
 "CYX"
 "GLH"
 "GLN"
 "GLU"
 "GLY"
 "HID"
 "HIE"
 "HIP"
 "HYP"
 "ILE"
 "LEU"
 "LYN"
 "LYS"
 "MET"
 "PHE"
 "PRO"
 "SER"
 "THR"
 "TRP"
 "TYR"
 "VAL"
!entry.ALA.unit.atoms table  str name  str type  int typex  int resx  int flags  int seq  int elmnt  dbl chg
! ALA 残基, 原子名称, 原子类型, 电荷
 "N" "N" 0 1 131072 1 7 -0.415700
 "H" "H" 0 1 131072 2 1 0.271900
 "CA" "XC" 0 1 131072 3 6 0.033700
 "HA" "H1" 0 1 131072 4 1 0.082300
 "CB" "CT" 0 1 131072 5 6 -0.182500
 "HB1" "HC" 0 1 131072 6 1 0.060300
 "HB2" "HC" 0 1 131072 7 1 0.060300
 "HB3" "HC" 0 1 131072 8 1 0.060300
 "C" "C" 0 1 131072 9 6 0.597300
 "O" "O" 0 1 131072 10 8 -0.567900
!entry.ALA.unit.atomspertinfo table  str pname  str ptype  int ptypex  int pelmnt  dbl pchg
 "N" "N" 0 -1 0.0
 "H" "H" 0 -1 0.0
 "CA" "XC" 0 -1 0.0
 "HA" "H1" 0 -1 0.0
 "CB" "CT" 0 -1 0.0
 "HB1" "HC" 0 -1 0.0
 "HB2" "HC" 0 -1 0.0
 "HB3" "HC" 0 -1 0.0
 "C" "C" 0 -1 0.0
 "O" "O" 0 -1 0.0
!entry.ALA.unit.boundbox array dbl
 -1.000000
 0.0
 0.0
 0.0
 0.0
!entry.ALA.unit.childsequence single int
 2
!entry.ALA.unit.connect array int
 1
 9
!entry.ALA.unit.connectivity table  int atom1x  int atom2x  int flags
! 成键连接, 用于推断成键项
 1 2 1
 1 3 1
 3 4 1
 3 5 1
 3 9 1
 5 6 1
 5 7 1
 5 8 1
 9 10 1
!entry.ALA.unit.hierarchy table  str abovetype  int abovex  str belowtype  int belowx
 "U" 0 "R" 1
 "R" 1 "A" 1
 "R" 1 "A" 2
 "R" 1 "A" 3
 "R" 1 "A" 4
 "R" 1 "A" 5
 "R" 1 "A" 6
 "R" 1 "A" 7
 "R" 1 "A" 8
 "R" 1 "A" 9
 "R" 1 "A" 10
!entry.ALA.unit.name single str
 "ALA"
!entry.ALA.unit.positions table  dbl x  dbl y  dbl z
! 每个原子的坐标, gmx不需要, 但可用于补全蛋白, 构建多肽
 3.325770 1.547909 -1.607204E-06
 3.909407 0.723611 -2.739882E-06
 3.970048 2.845795 -1.311163E-07
 3.671663 3.400129 -0.889820
 3.576965 3.653838 1.232143
 3.877484 3.115795 2.131197
 4.075059 4.623017 1.205786
 2.496995 3.801075 1.241379
 5.485541 2.705207 -4.398755E-06
 6.008824 1.593175 -8.449768E-06
!entry.ALA.unit.residueconnect table  int c1x  int c2x  int c3x  int c4x  int c5x  int c6x
 1 9 0 0 0 0
!entry.ALA.unit.residues table  str name  int seq  int childseq  int startatomx  str restype  int imagingx
 "ALA" 1 11 1 "p" 0
!entry.ALA.unit.residuesPdbSequenceNumber array int
 0
!entry.ALA.unit.solventcap array dbl
 -1.000000
 0.0
 0.0
 0.0
 0.0
!entry.ALA.unit.velocities table  dbl x  dbl y  dbl z
 0.0 0.0 0.0
 0.0 0.0 0.0
 0.0 0.0 0.0
 0.0 0.0 0.0
 0.0 0.0 0.0
 0.0 0.0 0.0
 0.0 0.0 0.0
 0.0 0.0 0.0
 0.0 0.0 0.0
 0.0 0.0 0.0

! ALA 结束, 开始下一个 ARG, 格式类似
!entry.ARG.unit.atoms table  str name  str type  int typex  int resx  int flags  int seq  int elmnt  dbl chg
... ...

19sb中蛋白残基种类总数为28, 除20种天然氨基酸残基及其不同质子化状态外, 还有一个非标准残基HYP. 每种残基的所有信息都已列出, 只要将这些信息处理一下, 就可以将其转换成

中的信息.

将结果与gmx自带的99sb, 或网上的14sb对比一下, 可以看到大部分都是一样的, 只有少量不同.

至于力场目录中的其他文件, 因为关系不大, 直接使用旧的就好了.

19sb的推荐使用的水模型是4位点OPC, 所以最好把这个水模型也一起放到力场文件中以方便使用.

修改gmx源码以支持amber19sb的CMAP

I have made this longer than usual because I have not had time to make it shorter.

—- Blaise Pascal

可行的修改方法无穷无尽, 我遵循最小改动原则, 改动越小越好, 因为改动越多越容易出bug, 也越容易引起兼容性问题. 可惜的是改动越小, 所花时间越多.

For the sake of brevity, I use the term “rays” and more especially “X rays” in order to distinguish them from others.

—- Wilhelm Conrad Rontgen

我就以自己熟悉并常用的gmx2019.6版本进行修改, 为避免与原始gmx想混淆, 我就称修改后的版本为gromacsX, 简写gmxx, 倒也合适, 因为我不知道后面还会有什么改动, 所以就留个未知数. 当然, 这个名称也有碰瓷的嫌疑.

gmx对cmap的支持并不好, 如果够好的话, 我们也就不需要改源码了. 对于其他成键项, 拓扑处理时都是支持#define宏定义的, 但cmap什么都不支持, 只有一个函数类型, 默认为1, 处理拓扑时这个函数类型也没有用到. 不过这恰好给我们留下了可操作的空间. 我们可以将其改为支持任意数字, 并通过相应的数字对应到到不同的cmap数据就可以达到目的了. 这样修改的话, 对使用cmap的charmm力场文件也要做相应的调整.

这样看下来, 要修改的源代码有两处:

测试

准备好amber19sb.ff力场文件, 编译好gmxx, 之后用起来就和以前一样了. 当然, 将gmxx生成的tpr文件用于更高版本的gmx也是可行的.

更严格的测试需要针对各种蛋白体系, 分别用gmxx和amber计算能量和力, 检查在误差范围内是否一致. 这就先缓一缓吧. 因为我有点累了.

◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆
前一篇:　秋天的告别
后一篇:　gmx_MMPBSA脚本更新：处理不连续索引
访问人次(2015年7月 9日起): | 最后更新: 2024-11-01 02:53:58 UTC | 版权所有 © 2008 - 2024 Jerkwin 个