| 查看: 10933 | 回复: 16 | ||||||||||||||||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||||||||||||||||
御剑江湖荣誉版主 (著名写手)
天池冶海
|
[交流]
【转帖】用Amber+Gaussian做小分子力场 已有14人参与
|
|||||||||||||||
|
用Amber+Gaussian做小分子力场，不是我的研究领域，不过阅读此篇文章，发现不错，转载给此研究领域的虫子们。。。 第一步：生成小分子模板 蛋白质中各氨基酸残基的力参数是预先存在的，但是很多模拟过程会涉及配体分子，这些有机小分子有很高的多样性，他们的力参数和静电信息不可能预存在库文件中，需要根据需要自己计算生成模板。amber中的antechamber 程序就是生成小分子模板的。生成模板要进行量子化学计算，这一步可以由antechamber中附带的mopac完成，也可以由gaussian完成，这里介绍用gaussian计算的过程。 建议在计算前用sybyl软件将小分子预先优化，不要用gaussian优化，大基组从头计算进行几何优化花费的计算时间太长。gaussian计算的输入文件可以用antechamber程序直接生成，生成后去掉其中关于几何优化的参数即可将小分子优化后的结构存储为mol2各式，上传到工作目录，用antechamber程序生成gaussian输入文件，命令如下： antechamber -i 49.mol2 -fi mol2 -o 49.in -fo gzmat 这样可以生成49.in文件，下载到windows环境，运行gaussian计算这个文件，如果发现计算时间过长或者内存不足计算中断，可以修改文件选择小一些的基组。 获得输出文件49.out之后将它上传到工作目录，再用antechamber生成模板，命令如下： antechamber -i 49.out -fi gout -o 49mod.mol2 -fo mol2 -c resp 运行之后就会生成一个新的mol2文件，如果用看图软件打开这个文件会发现，原子的颜色很怪异，这是因为mol2的原子名称不是标准的原子名称，看图软件无法识别。下面一步是检查参数，因为可能会有一些特殊的参数在gaff中不存在需要程序注入， 命令如下： parmchk -i 49mod.mol2 -f mol2 -o 49mod 这样那些特殊的参数就存在49mod这个文件中了 第二步：处理蛋白质文件 amber自带的leap程序是处理蛋白质文件的，他可以读入PDB各式的蛋白质文件，根据已有的力场模板为蛋白质赋予键参数和静电参数。 PDB格式的文件有时会带有氢原子和孤对电子的信息，但是在这种格式下氢原子和孤对电子的命名不是标准命名，力场模板无法识别这种不标准的命名，因此需要将两者的信息删除 ATOM 12 1H ARG A 82 12.412 8.891 34.128 1.00 0.00 H 在PDB各式里面，氢原子的信息会在第13或者14列出现H字符，可以应用grep命令删除在13或者14列出现H的行命令如下： grep -v ‘^………….H’ 1t4j.pdb > x grep -v ‘^…………H’ x > 1t4j_noh.pdb 除了删除氢和孤对电子，还应该把文件中的结晶水、乙酸等分子删除，这些分子的信息常常集中在文件的尾部，可以直接删除。 处理过之后的蛋白质文件，只包括各氨基酸残基和小分子配体的重原子信息，模拟需要的氢原子和水分子将在leap中添加接下来需要进一步整理蛋白质文件，主要关注氨基酸的不同存在形态和小分子的原子名称。半胱氨酸有两种存在形态，有的是两个半胱氨酸通过二硫键相连，有的则是自由的，两种半胱氨酸在力场文件中是用不同的unit来表示的，这相当于是两个完全不同的氨基酸，需要手动更改蛋白质文件中半胱氨酸的名字，桥连的要用CYX，自由的用CYS 组氨酸有若干种质子态，和半胱氨酸一样，也需要查阅文献确定这些质子态，并更改残基名称 最后需要修改的是配体分子的原子名，这是工作量最大的事情，仔细观察可以发现，一般PDB文件中配体的各个原子名称，和我们上面通过antechamber 生成的49mod.mol2中原子名称并不一致，这会造成leap无法识别这些原子，无法为这些原子赋予力参数和静电参数， 因此需要手动更改蛋白质文件中配体分子的原子名称。 进行这一步可以同时用看图软件打开49mod.mol2和蛋白质文件，隐藏蛋白质文件中除配体分子以外的所有结构，旋转两个文件，让他们姿态相近，以方便观察，并且在各自均标识原子名称，然后用文字编辑软件打开蛋白质文件，核对看图软件中两个分子对应的原子名称，手动逐一修改。 修改原子名称需要极大的耐心和细心，如果发生错误下一步的操作会无法继续。我现在想到的也只有这个笨办法，如果谁还有别的好法子，欢迎告诉我。 现在文件的准备工作都已经完成，该开始正式的模拟了 第三步：生成拓扑文件和坐标文件 用amber进行分子动力学模拟需要坐标和拓扑文件，坐标文件记录了各个质点所座落的坐标，拓扑文件记录了整个体系各质点之间的链接状况、力参数电荷等信息。这两个文件是由leap 程序生成的amber中有两个leap程序，一个是纯文字界面的tleap，一个是带有图形界面的Xleap。但是amber的图形界面做得很差，用远程登录使用图形界面不仅麻烦而且慢，所以我比较偏爱使用tleap，两个leap的命令是完全一样的，其实用哪一个都无所谓。 启动tleap，在shell里输入命令tleap，leap就启动了，shell里会显示 -I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/prep to search path. -I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/lib to search path. -I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/parm to search path. -I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/cmd to search path. Welcome to LEaP! (no leaprc in search path) > 这个>是leap的提示符 下面要调入库文件。amber是模拟生物分子的好手，主要就是依靠专门为蛋白质多糖核酸量身订做的amber力场，力场的所有参数都存 储在库文件里，所以打开leap第一件事便是调入库文件。 amber提供了很多种库，这里我们只用到两个库，gaff和02库，输入命令： >source leaprc.gaff >source leaprc.ff02 之后两个库就调入进来了 这时可以用list命令看看库里都有什么： 这里面罗列的就是库里面的unit，包括20种氨基酸、糖以及核酸还有一些常见离子的参数 下面一步是调入配体分子的模板，首先补全参数，输入命令： >loadamberparams 49mod 然后读入模板文件,输入命令： >MOL = loadmol2 49mod.mol2 其中MOL是unit的名字，要保证这个名字和pdb文件中配体的残基名完全一致，否则系统仍然无法识别pdb文件中的小分子 现在再输入list命令会发现库里面多了一个unit: 那个就是配体分子的模板 下面读入pdb文件，输入命令： >comp = loadpdb 1t4j_noh.pdb 如果输入这个命令之后，屏幕上出现了大量的创建unit或者atom的信息，如下所示，则说明上面一步的pdb文件处理没有做好， 还需要重新处理，通常这种情况都发生在配体分子上面，有时则是因为蛋白质中存在特殊残基。 Creating new UNIT for residue: FRJ sequence: 1 Created a new atom named: O36 within residue: .R Created a new atom named: S33 within residue: .R Created a new atom named: O35 within residue: .R Created a new atom named: N34 within residue: .R 如果屏幕仅仅显示添加原子，这说明输入的PDB文件中缺失了部分重原子，leap根据模板自动补全了这些缺失的原子，这种情况 不会影响进一步的计算 Added missing heavy atom: .R.A Added missing heavy atom: .R.A Added missing heavy atom: .R.A Added missing heavy atom: .R.A 根据体系的具体情况，还可能要将成对的CYX残基中的二硫键相连，有时候还会链接其他原子，比如将糖基链接在氨基酸残基上 ，用bond命令可以完成，命令用法如下： >bond comp.35.SG comp.179.SG 其中comp是刚才读入的分子名称，35和179是残基序号，SG是CYX残基模板中硫原子的名称，用comp.35.SG这样的语法就可以定 位一个原子 如果希望进行考虑溶剂效应的动力学模拟，可能还需要为体系加上水，加水有很多种命令，这里只列举一个： >solvatebox comp TIP3PBOX 10.0 solvatebox命令是说要加上一个方形的周期水箱，comp指要加水的分子，TIP3PBOX是选择的水模板名称，10.0是水箱子的半径 有的体系总电荷不为0，为了模拟稳定，需要加入抗衡离子，系统会自动计算体系的静电场分布，在合适的位置加上离子，命令如下： >addions comp Na+ 0 意思是用钠离子把体系总电荷补平，还可以选择其他库里面有的离子。 完成这一步之后就可以输出拓扑文件和坐标文件了，但是为了方便起见，在运行最后一步之前要先把leap里加工好的分子单独存 成一个库文件，以后可以直接调用这个库文件，免得重复上面的操作： >saveoff comp 1taj.off 这样就生成了一个off文件在那里，下面输出拓扑文件和坐标文件 >saveamberparm comp 1t4j.prmtop 1t4j.inpcrd Checking Unit. Building topology. Building atom parameters. Building bond parameters. Building angle parameters. Building proper torsion parameters. Building improper torsion parameters. total 1 improper torsion applied Building H-Bond parameters. Not Marking per-residue atom chain types. Marking per-residue atom chain types. (Residues lacking connect0/connect1 - these don’t have chain types marked: res total affected CMET 1 ) (no restraints) >quit 现在准备好拓扑文件和坐标文件，接下来就可以开始能量优化和动力学模拟了。如果愿意的话还可以用ambpdb这个命令生成一 个pdb文件，直观地看一看生成了什么东西，命令如下： [snowyowls@localhost actualamber]$ ambpdb -p 1t4j.prmtop <1t4j.inpcrd> kankan.pdb | New format PARM file being parsed. | Version = 1.000 Date = 09/08/06 Time = 16:36:09 [snowyowls@localhost actualamber]$ 可以下载之后用看图软件欣赏，如果加了溶剂盒子的话，看的时候可要小心，会恨吓人的样子。 第四步：能量优化 用amber进行分子动力学模拟需要坐标和拓扑文件，这在上一步已经完成了，分别生成了1t4j.prmtop 和1t4j.inpcrd两个文件 ，下面一步就是动力学模拟之前的能量优化了。 由于我们进行的起始结构来自晶体结构或者同源模建，所以在分子内部存在着一定的张力，能量优化就是在真正的动力学之前， 释放这些张力，如果没有这个步骤，在动力学模拟开始之后，整个分子可能会因此散架。 能量优化由sander模块完成，运行sander至少需要三个输入文件，分别是分子的拓扑文件，坐标文件以及sander的控制文件。 现在分子的拓扑文件和坐标文件已经完成，需要建立输入文件，min_1.in Initial minimisation of our structures &cntrl imin=1, maxcyc=4000, ncyc=2000, cut=10, ntb=1,ntr=1, restraint_wt=0.5 restraintmask=’:1-283′ / 文件首行是说明，说明这项任务的基本情况； &cntrl与/之间的部分是模拟的参数 其中imin=1表示任务是能量优化，maxcyc=4000表示能量优化共进行4000步，ncyc=2000表示在整个能量优化的4000步中， 前 2000步采用最陡下降法，在2000步之后转换为共轭梯度法，如果模拟的时候不希望进行方法的转换，可以再加入另一 个关键词NTMIN，如果NTMIN =0则全程使用共轭梯度法，NTMIN=2则全程使用最陡下降法，此外还有=3和=4的选项，分别是xmin法和lmod法 ，具体情况可以看手册。 第二行的cut=10表示非键相互作用的截断值，单位是埃， ntb=1表示使用周期边界条件，这个选项要和前面生成的拓扑文件坐标文件相匹配， 如果前面加溶剂时候用的是盒子水，就设置ntb=1，如果加的是层水，那就应该选择ntb=0；ntr=1表示在能量优化的过程中要约束一些原子， 约束的是哪些原子呢？后面有。 第三行和第四行都是关于约束原子的信息，restraint_wt=0.5限定了约束的力常数，在这里约束原子就是把原子用一根弹簧拉在固定的 位置上，一旦原子偏离固定的位置，系统就会给他施加一个回复力，偏离的越远，回复力越大，回复力就是由这个力常数决定的，单位是 Kcal/(mol*A)。 restraintmask=’:1-283′表示约束的是从1到283号残基，在这个分子中，1-283号残基是蛋白质上的氨基酸残基，从284 号开始，就都是水了，所以这个命令的意思就是，约束整个蛋白质，自由地优化溶剂分子。因为溶剂分子是前面的tleap自动给加上的， 所以一定要额外多关注一些。 下面运行sander来执行这个优化： [snowyowls@localhost actualamber]$ sander -O -i min_1.in -p 1t4j.prmtop -c 1t4j.i npcrd -ref 1t4j.inpcrd -r 1t4j_min1.rst -o 1t4j_min1.out 命令中，-O表示覆盖所有同名文件，-i min_1.in表示sander的控制文件是min_1.in，-p 1t4j.prmtop表示分子的拓扑文件 ，-c 1t4j.inpcrd表示坐标文件，-ref 1t4j.inpcrd是参考坐标文件，只有在控制文件中出现关键词ntr=1的时候才需要给 sander指定-ref文件，这是约束原子的参考坐标，- ref 1t4j.inpcrd就是说以1t4j.inpcrd中的坐标为准进行约束原子的优化。 以上这四个是输入文件。-r 1t4j_min1.rst表示经过模拟之后新的原子坐标会输出到1t4j_min1.rst文件中，-o 1t4j_min1.out 则表示优化过程中的相关信息都会写入到1t4j_min1.out文件中。 运行起这个命令之后，等拿到 1t4j_min1.rst文件就意味着对溶剂的优化已经差不多了，显然下面还需要对蛋白质本身进行优化， 这个优化还要分两步进行，控制文件分别是min_2.in 和min_3.in min_2.in Initial minimisation of our structures &cntrl imin=1, maxcyc=5000, ncyc=2500, cut=10, ntb=1,ntr=1, restraint_wt=0.5 restraintmask=’:1-283@CA,N,C’ / 在这里发生变化的是约束原子的范围， ‘:1-283@CA,N,C’表示1-283号残基中名叫CA，N和C的原子，这些原子实际上是蛋白质主链上 的原子，也就是说这一次的优化是约束了蛋白质主链上的原子之后，对溶剂和侧链原子进行自由优化。 min_3.in Initial minimisation of our structures &cntrl imin=1, maxcyc=10000, ncyc=5000, cut=10, ntb=1, / 在这里不再进行约束原子的优化了，对整个分子进行全原子优化。 三步优化的命令如下： [snowyowls@localhost actualamber]$ sander -O -i min_1.in -p 1t4j.prmtop -c 1t4j.inpcrd -ref 1t4j.inpcrd -r 1t4j_min1.rst -o 1t4j_min1.out [snowyowls@localhost actualamber]$ sander -O -i min_2.in -p 1t4j.prmtop -c 1t4j_min1.rst -ref 1t4j_min1.rst -r 1t4j_min2.rst -o 1t4j_min2.out [snowyowls@localhost actualamber]$ sander -O -i min_3.in -p 1t4j.prmtop -c 1t4j_min2.rst -r 1t4j_heat0.rst -o 1t4j_min3.out 接下来就是真正的分子动力学模拟了…… |
回复此楼» 收录本帖的淘帖专辑推荐
 物理前沿 | xuexijisuan | 分子模拟要术 | amber |
| 分子模拟 | 精华帖转载 | 分子对接 | Amber |
| 计算化学 | 搜集 | 申请 | 精华 |
| 分子动力学模拟 |
» 本帖已获得的红花（最新10朵）
wl507586» 猜你喜欢
售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急
已经有8人回复
-
体制内长辈说体制内绝大部分一辈子在底层，如同你们一样大部分普通教师忙且收入低
已经有13人回复
-
售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急
已经有8人回复
-
售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急
已经有8人回复
-
售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急
已经有11人回复
-
为什么中国大学工科教授们水了那么多所谓的顶会顶刊，但还是做不出宇树机器人？
已经有8人回复
-
售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急
已经有5人回复
-
售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急
已经有6人回复
-
售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急
已经有9人回复
-
售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急
已经有4人回复
» 本主题相关商家推荐: (我也要在这里推广)
» 本主题相关价值贴推荐，对您同样有帮助:
- 怎么处理UNIT MISSING 的错误？
已经有4人回复
- 分子力场中非键能
已经有4人回复
- 关于pdb库使用问题
已经有8人回复
- 关于MS力场编辑
已经有11人回复
- Amber 分子动力学 求助！！！！
已经有10人回复
- 求计算一个简单小分子的尺寸-用Gaussian
已经有4人回复
- Gromacs中添加CHARMM力场
已经有41人回复
- 使用amber力场，一些参数无法获得
已经有5人回复
- 高斯产生小分子参数问题
已经有6人回复
- 关于用amber跑REMD计算的流程及输入文件的设置
已经有6人回复
- 如何将sybyl对接后的蛋白质和配体保存成一个PDB文件
已经有8人回复
- 怎么自己做力场-及气体力场的选择
已经有6人回复
- 一个桥连H的问题
已经有6人回复
- ATP分子力场参数怎么做
已经有4人回复
- GAUSSIAN计算语法错误
已经有7人回复
- 糖蛋白---分子动力学模拟
已经有4人回复
- pdb文件中加入氢原子
已经有7人回复
- 一个模拟体系中能否用两个不同的力场？
已经有7人回复
- Gaussian模拟SERS
已经有6人回复
- 【求助】Amber中水分子的问题
已经有4人回复
- 【讨论】gaussian优化过程中具体在执行什么操作
已经有3人回复
- 【求助】量子化学方法求amber力场参数
已经有10人回复
- 【讨论】Gaussian 09调用内存的问题
已经有3人回复
4楼2011-04-10 23:10:46
xd200620940
铜虫 (正式写手)
- 模拟EPI: 2
- 应助: 0 (幼儿园)
- 金币: 753.9
- 红花: 5
- 帖子: 508
- 在线: 137.7小时
- 虫号: 494350
- 注册: 2008-01-13
- 专业: 化学计量学与化学信息学
★ ★ ★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
御剑江湖(金币+3): 谢谢！ 2011-04-10 11:28:01
小木虫(金币+0.5):给个红包，谢谢回帖交流
御剑江湖(金币+3): 谢谢！ 2011-04-10 11:28:01
|
最后需要修改的是配体分子的原子名，这是工作量最大的事情，仔细观察可以发现，一般PDB文件中配体的各个原子名称，和我们上面通过antechamber 生成的49mod.mol2中原子名称并不一致，这会造成leap无法识别这些原子，无法为这些原子赋予力参数和静电参数， 因此需要手动更改蛋白质文件中配体分子的原子名称。 进行这一步可以同时用看图软件打开49mod.mol2和蛋白质文件，隐藏蛋白质文件中除配体分子以外的所有结构，旋转两个文件，让他们姿态相近，以方便观察，并且在各自均标识原子名称，然后用文字编辑软件打开蛋白质文件，核对看图软件中两个分子对应的原子名称，手动逐一修改。 这个是用antechamber做小分子力场最为困难的地方，相当的不人性化！ |
2楼2011-04-10 11:25:22
3楼2011-04-10 18:00:41
xd200620940
铜虫 (正式写手)
- 模拟EPI: 2
- 应助: 0 (幼儿园)
- 金币: 753.9
- 红花: 5
- 帖子: 508
- 在线: 137.7小时
- 虫号: 494350
- 注册: 2008-01-13
- 专业: 化学计量学与化学信息学
★ ★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
御剑江湖(金币+2): 期待专帖详细介绍 2011-04-12 12:15:25
小木虫(金币+0.5):给个红包，谢谢回帖交流
御剑江湖(金币+2): 期待专帖详细介绍 2011-04-12 12:15:25
|
因为之前只用GROMACS，对AMBER了解不多。最近GROMACS将AMBER力场包含了进来，对于配体小分子的转化在网站的教程中也给出了一个方案（但没有实例）。 根据该教程提到的方案，可以采用ACPYPE + ANTECHAMBER得到GROMACS格式的配体分子topology和GRO，因此上周我对该问题进行了系统地摸索。得到几点体会： 1.antechamber做小分子参数时转化成功与否与原子类型无关； 2.antechamber转化小分子参数时要求PDB文件是氢饱和的，即从PDB库中直接得到的小分子需加氢，才能进行参数计算； 3.若是采用绘图软件得到的PDB文件，即使氢已饱和，直接使用antechamber做参数也不会成功，卡在leap阶段。经过研读教程，再比较pdb文件，找到了问题的关键。也得出第1点的结论。 关于antechamber以及acpype具体的使用步骤，近两天我将开一个专帖详细介绍，敬请虫友们关注，谢谢！ [ Last edited by xd200620940 on 2011-4-11 at 19:32 ] |
5楼2011-04-11 19:30:47
