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

yulaile

主要功能：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

复制代码

     安装ztop http://bbs.keinsci.com/thread-22171-1-1.html

     安装 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
   

复制代码