计算化学公社

 找回密码 Forget password
 注册 Register
Views: 161|回复 Reply: 0
打印 Print 上一主题 Last thread 下一主题 Next thread

[GROMACS] 辅助生成自定义残基的单链高分子opls力场拓扑的python脚本

[复制链接 Copy URL]

2

帖子

1

威望

1060

eV
积分
1082

Level 4 (黑子)

本帖最后由 yulaile 于 2024-11-14 17:35 编辑

主要功能:pdb2gmx+Ztop+LigParGen 生成opls力场的单链聚合物链拓扑
  • RTP文件制作:     LigParGen服务器下载的拓扑文件,输入切分残基和原子等价的原子序号,生成opls-aa力场的RTP文件
  • PDB文件制作:    这里借助钟成老师的ztop软件的大分子构建的功能生成pdb文件。因为原子命名与LigParGen的原子命名没有对应。需要通过ztop生成寡聚物的pdb(上传给LigParGen的寡聚物结构)与Ligpargen的pdb进行坐标重叠生成映射表。然后用ztop拼接大长链pdb文件,替换其原子命名与RTP对应。
  • 仅限简单一维单链的分子拓扑

主要流程下图所示:



安装库
  1. pip install biopython biopandas
  2. pip install pandas numpy
复制代码
     安装 Multiwfn vmd
     设置修改tools目录下的ztop_path文件的ztop根目录路径
     使用之前请备份gromacs的opls拓扑目录  路径:gromacs目录/share/gromacs/top/oplsaa.ff
使用步骤
          因为Ligpargen的上传文件要求原子数不能超过200 ,所有这里以两个分子苯乙烯和乙烯醇寡聚物为例,ps目录是苯乙烯分子的LigParGen下载的文件和原子分割的原子序号信号,同理oh为乙烯醇的相关文件 tools为脚本文件所在目录

生成苯乙烯的RTP步骤
  1. ./tools/1_depart_noorder.sh ./ps/MOL_23525C.pdb  ./ps/frag.txt
复制代码
#参数1:ligpargen下载的pdb   参数2:残基分割原子序号从左到右;生成残基的信息保存在depart目录下

  1. ./tools/2_crete_res.sh 'ps' PS
复制代码
#参数1:残基文件目录  参数2:残基名(多个残基名的首个字母不能相同);当前目录会输出带残基信息的test_res_ps.pdb  
#执行vmd test_res_ps.pdb -e tools/resname_vmd  查看残基分割是否正确

  1. conda activate ztop
  2. ./tools/3_ztop_depart.sh ps/ph.log A   
复制代码
#参数1:计算寡聚物的高斯文件 参数2:第一个残基在编号(字母A-Z); 生成ztop的残基信息

  1. ./tools/4_insert_equal_atom.sh  test_res_ps.pdb test_res_ps.pdb ps/atom.txt
复制代码
#参数1:pdb文件 参数2:pdb文件 参数3:原子等价信息;第一列参数1文件的原子序号 第二列为参数2文件的原子序号 成映射表等价原子的原子名映射

  1. ./tools/5_itp.sh ps/MOL_23525C.itp
复制代码
#参数1:指定残基所在的itp文件;ps/creat_itp目录下生成 RTP文件、atomtypes.atp、residuetypes.dat、ffnonbonded.itp 等文件

  1. ./tools/6_install.sh 'ps'
复制代码
#参数1:残基目录;检查生成的rtp是否正确 进入在ps/pdb目录
#执行 gmx pdb2gmx -f order.pdb -ter 依次输入 15 8 3 3 ;看是否能正常生成top

  1. ./tools/7_pdb_atom_map.sh 'ps'
复制代码
#参数1:残基目录;映射ztop生成的pdb原子名与Ligpargen拓扑的原子名
#由于ztop生成的pdb原子的坐标和原子名都与Ligpargen的不同,但是可以通过Bio库选取三个依次对应原子进行坐标变化,这样每个原子名就可一一对应。tools/creat_tran_by_two.py 会自动找到两个分子各自对应的三个原子。对一个分子的所有原子与各自中心坐标的距离列表进行分箱找到刚好有三个为1,然后通过3个距离找到各自三个对应原子; (运行后有报错执 crete_trans_file.py FRAG.pdb atom1,2,3 order.pdb atom1,2,3 手动标注原子序号)
#可以执行pdb2gmx检查pdb_atom_map/test_pdb2gmx.pdb,最终生成两个pdb文件的原子名映射表

第二个残基同理执行
  1. ./2_run_oh.sh
复制代码

两个残基的RTP已添加到gromcas的Top目录,先使用ztop生成高分子链的pdb  再替换原子名生成与RTP匹配的PDB文件   
  1. ./tools/8_creat_poly_pdb.sh  A[BBFBB]10D MyChain  #B为苯乙烯 F为乙烯醇 AD为端基
复制代码

  1. ./tools/9_writepdb.sh run/MyChain.pdb
复制代码
# gmx pdb2gmx -f ./run/output.pdb -ter 依次输入15 8 3 3 测试是否正确生成高分子链top     

进行模拟
  1. ./4_run_gmx_change_nozero.sh
复制代码

     至此,编写这些脚本的目的是想学习理解测试gromacs的pdb2gmx运行流程,利用高效灵活的ztop高分子残基组装功能,快速建立opls力场的高分子top。此外可以通过修改代码里面的不同力场的参数设置,不仅限于opls力场,pdb2gmx支持的立场都能以相同方法建立高分子链top。希望这个程序帮助有需要理解gromacs的rtp规则的同学,此外程序可能存在居多问题还请各位老师多多指教。

由于高斯文件较大,相关文件放在github
  1. git clone -b new https://github.com/Yulaile-01/MD.git
复制代码




评分 Rate

参与人数
Participants 2
威望 +1 eV +5 收起 理由
Reason
danyj + 5
sobereva + 1

查看全部评分 View all ratings

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2024-11-27 15:30 , Processed in 0.572350 second(s), 25 queries , Gzip On.

快速回复 返回顶部 返回列表 Return to list