计算化学公社

标题: 辅助生成自定义残基的单链高分子opls力场拓扑的python脚本 [打印本页]

作者
Author: yulaile 时间: 2024-9-15 03:14
标题: 辅助生成自定义残基的单链高分子opls力场拓扑的python脚本
本帖最后由 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对应。
仅限简单一维单链的分子拓扑

主要流程下图所示：

(, 下载次数 Times of downloads: 3)

安装库

pip install biopython biopandas
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为脚本文件所在目录

(, 下载次数 Times of downloads: 3)

生成苯乙烯的RTP步骤

./tools/1_depart_noorder.sh ./ps/MOL_23525C.pdb ./ps/frag.txt

复制代码

#参数1：ligpargen下载的pdb 参数2：残基分割原子序号从左到右；生成残基的信息保存在depart目录下

(, 下载次数 Times of downloads: 5)

./tools/2_crete_res.sh 'ps' PS

复制代码

#参数1：残基文件目录参数2：残基名（多个残基名的首个字母不能相同）；当前目录会输出带残基信息的test_res_ps.pdb

#执行vmd test_res_ps.pdb -e tools/resname_vmd 查看残基分割是否正确

conda activate ztop
./tools/3_ztop_depart.sh ps/ph.log A

复制代码

#参数1：计算寡聚物的高斯文件参数2：第一个残基在编号（字母A-Z）；生成ztop的残基信息

./tools/4_insert_equal_atom.sh test_res_ps.pdb test_res_ps.pdb ps/atom.txt

复制代码

#参数1：pdb文件参数2：pdb文件参数3：原子等价信息；第一列参数1文件的原子序号第二列为参数2文件的原子序号生成映射表等价原子的原子名映射

./tools/5_itp.sh ps/MOL_23525C.itp

复制代码

#参数1：指定残基所在的itp文件；ps/creat_itp目录下生成 RTP文件、atomtypes.atp、residuetypes.dat、ffnonbonded.itp 等文件

./tools/6_install.sh 'ps'

复制代码

#参数1：残基目录；检查生成的rtp是否正确进入在ps/pdb目录

#执行 gmx pdb2gmx -f order.pdb -ter 依次输入 15 8 3 3 ；看是否能正常生成top

./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文件的原子名映射表

第二个残基同理执行

./2_run_oh.sh

复制代码

(, 下载次数 Times of downloads: 4)

两个残基的RTP已添加到gromcas的Top目录，先使用ztop生成高分子链的pdb 再替换原子名生成与RTP匹配的PDB文件

./tools/8_creat_poly_pdb.sh A[BBFBB]10D MyChain #B为苯乙烯 F为乙烯醇 AD为端基

复制代码

./tools/9_writepdb.sh run/MyChain.pdb

复制代码

# gmx pdb2gmx -f ./run/output.pdb -ter 依次输入15 8 3 3 测试是否正确生成高分子链top

进行模拟

./4_run_gmx_change_nozero.sh

复制代码

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

由于高斯文件较大，相关文件放在github

git clone -b new https://github.com/Yulaile-01/MD.git

复制代码

(, 下载次数 Times of downloads: 4)

欢迎光临计算化学公社 (http://bbs.keinsci.com/)