用ztop构建大分子/聚合物的拓扑文件

ggdh

这个帖子介绍如何用ztop构建大分子的拓扑文件。ztop的安装以及基本用法参见：构建及优化[小大]分子力场参数脚本ztop.py
用ztop拟合二面角可以参考
使用ztop拟合二面角参数


程序基本逻辑和输入输出：一句话逻辑：用户给小分子的拓扑和结构文件，从中切割出片段，然后用这些片段组合成大分子
输入：
包含片段结构的小分子的拓扑文件和坐标文件，支持gromacs，amber，charmm格式, 本程序可以直接生成gaff2 和amber力场的拓扑文件，如果是别的力场，可以用其他手段产生小分子拓扑文件，然后用本程序合并成大分子。
程序运行：
用户通过指定的键去切割输入的分子，得到片段，同一个分子可以有不同的切法得到不同的片段，也可以从不同的分子得到不同的片段。这些片段可以立马使用，或者存在硬盘上日后再用。程序通过用户的指令将不同的片段组合起来产生大分子（同时产生拓扑文件和坐标文件）。
输出：
大分子的拓扑文件和坐标文件，支持gromacs，amber，charmm格式。

Tutorial
这里用一个实例来详细讲解操作过程，过程虽然很长，但是逻辑很简单，兄弟们跟着做一遍应该就会了。
1.所需示例文件
需要用到的文件在程序包里面的ztop_test文件夹下的PI_poly_init目录下，如果是从本论坛下载的压缩包，会缺乏两个fchk文件，需要用户自己提交里面对应的gjf文件获得。有条件的直接运行

1. git clone https://gitee.com/coordmagic/coordmagic.git

复制代码

可以得到完整的例子文件。


2. 系统介绍

目标分子是如图所示的共轭聚合物，假设我们这里要做一个10聚体。我们可以制作10个单体片段 (由上图中4个共轭片段组成DON-PPI-ACC-PPI) 然后拼接起来。但是这里为了演示程序功能。我们把每一个独立的共轭片段作为一个片段来进行组合，因此这个聚合物由这里的DON，PPI和ACC三种片段构成。也就是所谓的D-PI-A-PI型聚合物。这里使用共轭片段的另外一个好处是每个共轭片段可以有自己的resname，方便后期vmd可视化或者是轨迹分析时利用。此外除了这三个种重复片段以外，我们还需要一个头和一个尾，因为这里做的不是周期性分子，而是寡聚物。

3.制作片段DON

ztop_test\PI_poly_init目录下的DON.gjf, DON.log, DON.fchk对应的结构如图所示，我们在DON的片段两边加上噻吩模拟这个片段在聚合物中的情况，其实更严谨的做做法是从更长的一段寡聚物中间挖出这个DON片段，比如从DPAPDPAPD中把其中标红的DON基团挖出来，这样得到的电荷和力场参数可能更符合聚合物中的情况。但是这里为了演示，用了一个比较短的片段。
首先产生top和gro文件，顺便优化一下平衡键长，这里使用了默认的resp电荷和gaff2力场。

1. cp PI_poly_init/DON.log ./
   
2. cp PI_poly_init/DON.fchk ./
   
3. ztop.py -g "DON.log" -r e -o DON.top,DON.gro

复制代码

然后从top和gro文件出发，切割片段，命令如下

1. ztop.py -f "D;p=DON.top;x=DON.gro;site=D12:11-22,D11:8-21;resname=DON;comment=benzodithiophene" --savelib

复制代码

-f 选项是开启片段制作模式，后面一长串指定如何制作片段，第一个字母D是片段的标签，用一个字母代表，后面用分号分隔不同的选项：
p=DON.top    指定top文件
x=DON.gro    指定gro文件
site=D12:11-22,D11:8-21    这里产生了两个链接位点，切断11-12键产生D12位点，切断8-21键产生D11位点，这里需要注意键的原子顺序是片段内-片段外的顺序，另外位点的名称决定了后面不同基团之间接是如何链接的。对于这个分子而言它的两个位点是等价的，但是本程序不允许两个位点名称相同，因此分别命名为D11和D12，这样这两个位点是1级等价的，参见后面的位点匹配规则。
resname=DON    设定这个片段的resname，会写到后面输出的top文件中
comment=benzodithiophene     加上注释方便用户记住这个片段是个啥，尤其片段多，或者是保存片段下次用的情况，注释不影响输出文件。
--savelib 选项把-f制作好的片段保存在当面目录下的FRAG_LIB下，其实就是把DON.top, DON.gro拷贝到该目录下，然后创建一个名为DON.frg的文件，把-f后面那一长串写进去。

程序最后会告诉你，当前库中有哪些片段了。

4.连接位点命名和连接规则
从分子产生片段时，每切断一个键就会产生一个连接位点，这些连接位点用来和其他片段连接。
两个片段之间用什么位点连接，取决于连接位点的命名，规则如下。
a.第一个字符通常和片段的标签一致，不用来匹配连接位点
b.从第二个字符开始，采取最短匹配规则，比如D1和A1，A12，A123都可以匹配而D2不行。 D13可以和A1，A13，A131匹配，但是不能和A2，A12和A123匹配。因此，一个连接位点标签的长度越短，这个位点就越随便。如果一个位点的名称是A，那么它可以去和所有其他位点匹配。
c.一个片段上不能有名字一样的位点
d.比如一个片段上有位点A1，A2，另外一个片段上有位点A1，B2，这时候第一个字符不同的优先度更高，也就是A2-B2的优先度比A1-A1的优先度高。

5.制作片段ACC

如果之前没有产生top文件，那么首先还是做top

1. cp PI_poly_init/ACC.log ./
   
2. cp PI_poly_init/ACC.fchk ./
   
3. ztop.py -g "ACC.log" -r e -o ACC.top,ACC.gro

复制代码

然后产生片段

1. ztop.py -f "A;p=ACC.top;x=ACC.gro;site=A21:4-11,A22:1-10;resname=DON;comment=benzothiadiazole" --savelib

复制代码

注意这里连接位点的编号是A21和A22，根据前面讲的规则，这两个位点是不会直接和DON片段上的D11和D12位点连接的。

6.制作PPI片段和头尾片段

这里我们用一个分子来制作3个片段，这个是一个D-PI-A结构的分子，中间的PI我们切下来作为聚合物中间的重复PI片段。左边的D我们切下来作为聚合物的头，右边的A我们切下来作为聚合物的尾巴。
命令如下，如果之前没有产生PPI的top文件，那么第一步还是制作PPI的top文件：

1. cp PI_poly_init/PPI.log ./
   
2. cp PI_poly_init/PPI.fchk ./
   
3. ztop.py -g "PPI.log" -r e -o PPI.top,PPI.gro

复制代码

然后开始制作3个片段：

1. ztop.py -f "P;p=PPI.top;x=PPI.gro;site=P1:15-8,P2:18-20;resname=PPI;comment=thiophene" \
   
2. -f "H;p=PPI.top;x=PPI.gro;site=H1:8-15;resname=DON;comment=D as head" \
   
3. -f "T;p=PPI.top;x=PPI.gro;site=T2:20-18;resname=ACC;comment=A as tail" --savelib

复制代码

这里分别用3个-f来制作3个片段，输入这么长的命令容易出错，其实也可以一个片段一个片段的savelib，这是更稳妥的做法。
中间PPI片段的两个位点分别是P1和P2，根据之前讲的规则P1会和D11或D12连接，而P2会和A21和A22连接。
注意这里我们的头片段用的标签是H，但是他的resname还是DON。因为H是待会建聚合物时需要用的代号，需要和中间的D片段区别开，而从化学结构上来说头片段和聚合物中间的D片段一样，因此可以用相同的resname，当然这里完全可以用不同的resname，比如你特别关心头尾的之间的距离，那么可以把头尾设为单独的resname方便测量。

7.查看片段库
到这一步，我们的5个片段都做好了，分别是H，D，P，A，T， 我们可以把它们全部载入到一起看看，输入：

1. ztop.py --loadlib

复制代码

程序会执行FRAG_LIB目录中所有的frg文件中对片段的定义，最后输出

这里列出了所有片段的标签，分子式，位点，以及注释。
这里解释一下位点的意思，比如第一行的T2:C2-C2的意思是，这个位点名是T2，他是切断一个C-C键得到的，其中内部的C连了2个原子，外部碳也连了2个原子。这里部分的反应了内外接头的信息。而内外接头是连接不同片段时把新片段放到合适位置用的，下面会详细介绍。

8.见证奇迹
下面我们开始做聚合物，运行命令：

1. ztop.py --loadlib -b HPAPDPT

复制代码

尝试得到一个小短链，看能否顺利运行，运行命令

1. ztop.py --loadlib -b HP[APDP]10T -o poly.top,poly.pdb

复制代码

得到一个10聚体，用gaussview打开pdb可以看到


9.运行逻辑
这里来复述一下机器的运行逻辑，方便大家了解背后过程
1）从H片段开始，把H片段作为基础片段，H片段上有接口H1。
2）载入P片段，P片段上有接口P1，P2，此时发现P1和H1适配，进行连接，产生新的基础片段，此片段上有接口P2
3）载入A片段，A片段上有接口A21，A22，这两个接口都和基础片段上的P2适配，此时随意选择A21和P2连接，产生新的基础片段，片段上剩余的接口有A22
4）载入P片段，P片段上的P2和A22连接，还剩P1接口
5）载入D片段，D片段上的D11和P1连接，还剩D12接口
6）载入P片段。。。。。。
据此逻辑，大家可以尝试用已有的片段构建不同的大分子，比如：

1. ztop.py --loadlib -b PPPPPPP -o poly.top,poly.pdb

复制代码

得到聚噻吩，当然这里会遗留开放接口，开放接口如果不封闭默认会用片段的来源分子封闭。

1. ztop.py --loadlib -b HPAPA -o poly.top,poly.pdb

复制代码

这个会报错，大家根据刚才讲的组合规则想想为什么。
以上是举的一个线性的例子。其实也可以有树枝状的例子，只要新载入的片段和基础片段上有能够接上的接口，程序就能运行下去，否者会报错。
下面举一个带分叉的例子。
假如上面这个聚合物中DON上的烷基链很长，我想把这两个烷基链从DON上切下来作为独立的片段R，这样片段DON上会出现4个接口
分布是D11，D12和D31，D32，然后烷基链片段的定义为R，接口命名为R3，那么最后建造10聚物的命令是：

1. ztop.py --loadlib -b HRRP[APDRRP]10T -o poly.top,poly.pdb

复制代码

中间我省略了制作片段的步骤，大家可以练习一下。
另外也可以制作dendrimer，比如核心，1代，2代，3代的片段分别是C，L，M，N，接口分别是C11，C12，C13，L1，L21，L22，M2，M31，M32，N3
那么最后的制作命令应该是

1. ztop.py --loadlib -b CLLL[M]6[N]12 -o poly.top,poly.pdb

复制代码

10.接口参数处理，以及电荷处理
连接处的键参数取两个片段的平均。
键角各取各的。
二面角4个原子，如果只有一个原子在对方片段上，则取本片段参数。
如果有2个原子在对方片段上，则取先加入片段的二面角参数，比如-b DP  这个命令产生的片段中连接DP的二面角的参数就是来自D中的二面角参数。
另外可以用参数转移功能进一步修正接口处的参数，详见下一节
由于片段是从分子中切割下来的。因此各种片段的电荷不一定能中和彼此，导致最后的聚合物产生净电荷。
本脚本目前的做法是简单粗暴的将这个净电荷平均分配到每个原子之上。通常每个原子的电荷调整在0.001以下
这个电荷调整会输出在屏幕上，如果发现这个电荷很大不能接受，可以尝试
1）从更大的分子中切割得到片段。
2）利用参数转移功能将较小分子算得的电荷转移过去（比如算三聚体中部单体的电荷，然后转移到十聚体上）下面会介绍参数转移功能如何用。
3）把情况告诉我，开发更好的算法。

11. 参数转移

如图所示，我拟合了8-15这个可旋转键的二面角，现在我想把这个二面角的参数转移到刚才建立的那个10聚体上。
那个10聚体上有21个这样的二面角（每个DON边上有2个噻吩，加上开头的DON边上的一个噻吩），本程序可以通过结构匹配来实现自动转移功能，操作如下：
首先我们产生一个拟合好二面角的top文件

1. cp PI_poly_init/PPI.log ./
   
2. cp PI_poly_init/PPI.fchk ./
   
3. cp PI_poly_init/PPI_DFit_8-15.xyz .
   
4. ztop.py -g PPI.log -r e -d 8-15 -o PPI.top,PPI.gro

复制代码

这里PPI_DFit_8-15.xyz 这个文件中包含了扫描二面角获得的能量信息，如何获得这个文件在后面专门讲二面角拟合的帖子详细介绍。
如果之前10聚体的top搞没了，我们再重做一个，否者下面这步省略：

1. ztop.py --loadlib -b HP[APDP]10T -o poly.top,poly.pdb

复制代码

确定匹配部分，上图中绿色圈圈，圈住的部分是匹配的部分，我们打开Gaussview，选中这部分原子，用Tools --> Atom Selection...获得其编号为：
1-2,7-9,15-19,29,38-39
然后我们想转移参数的原子是图中蓝色有蓝色光晕的原子，注意这里转移参数的原子必须包含在匹配的原子之内，因为8-15这个键的旋转涉及到这6个原子，共4个二面角，同样通过Gaussview获得这部分的编号是：
7-9,15-16,19
准备好这些之后开始参数转移：

1. ztop.py -f "S;p=PPI.top;x=PPI.gro" -p poly.top -x poly.pdb  -t "S;m=1-2,7-9,15-19,29,38-39;a=7-9,15-16,19;p=d"

复制代码

其中-f定义片段，但是这里不用切割了，直接把PPI.top载入为一个名为S的片段。
-p和-x选项载入聚合物的拓扑和坐标文件
-t选项设定如何参数转移，第一个字符S表明需要从S片段来转移参数，
m=1-2,7-9,15-19,29,38-39意思是S片段中的这些原子用来进行结构匹配（m for match）
a=7-9,15-16,19意思是匹配结构中的这些S片段中的原子的参数将被转移到大分子上（a for atom）
p=d的意思是，转移包含这些原子的二面角参数（p for param， d for dihedral），这里我们可以用p=a,b,d,i,c来同时转移angle,bond,dihedral,improper,charge
最后程序会输出：

它这里告诉你一个找到了21个匹配，转移了21X4个二面角参数，比较了转移前后的优化能量，以及优化后结构。

12.对接规则
下面介绍一下当B片段接到A片段上以后，B片段的位置如何确定的。
大家先看CON片段的那个图，那个图中，我对连接位点D11标出了内外接口(inner adaptor 和out adaptor)。
这个接口就是由切断的那个键的原子和他的相邻的原子组成。

当把片段B和片段A拼接时，会首先比较B的外接口-A的内接口，以及A的外接口-B的内接口这两种情况的匹配度。
匹配度由接口原子的元素，数量，以及连接情况，是否相同等决定，优先选取匹配度高的情况。假如这里B的外接口和A的内接口匹配度更高，那么程序就会调整片段B的坐标，使之外接口的原子同A的内接口的原子重合。

13.库的管理
用-f选项制作的片段放在当前目录下的FRAG_LIB中，每个片段由3个文件组成，拓扑文件，坐标文件，以及frg文件。
frg文件中包含了当时制作片段的-f 选项后面的那一行字符串。
如果用户想重复使用某个片段，可以把片段对应的三个文件保存起来，需要使用的时候拷贝到新项目的FRAG_LIB的目录下。
修改frg文件中的片段label，以及site label来适配新的体系。
然后使用--loadlib ，此时程序会读取FRAG_LIB目录下所有的frg文件中的字符串。

14.同类软件介绍
polypargen
它的运行逻辑是你提交一个大分子给他，它帮你切片（切成<200)个原子的方便提交给LigParGen在线服务器，然后它再把这些片段组合起来交给你。
这个的自动化程度比我这个还要高。毕竟我这个要手动设计切片，手动组装。相应的我这个也更加灵活一些。另外这个只有在线服务器，某个聚合物能否成功比较看运气，参数质量如何不好控制。
pyPolyBuilder
逻辑和本程序类似，也是自己通过片段构建大分子，不过需要手动创建片段（它叫做BBs）的itp文件，有点像构建Gromacs的.rtp文件，比较麻烦。然后它让用户指定dendrimer的内核片段，中间片段，结尾片段，以及一共有几代，它会自己算每一代需要多少个什么片段，然后组装起来。对于其他复杂大分子，用户需要创建一个连接文件，里面要写什么片段和什么片段用什么原子连接。我感觉用起来比我这个麻烦一点。不过它可以做网状结构（就是包含环状的top），我这个目前还不行。
最重要的是，使用本软件构建大分子遇到问题可以享受作者在线反馈，另外如果你觉得目前的设计有什么sb的地方，或者有啥好的功能改进，方便化改进，都欢迎在本贴留言，或者加QQ：32589927

lonemen

钟老师终于肯发布啦。谢谢分享！

xiaowu759

谢谢分享！

naoki

好东西，mark一下

aw10279

钟老师，关于聚合体电荷的问题想请教下：
1. 如果单体本来就带电荷，聚合体中也希望其体现出电荷，是否就需要使用电荷转移功能来实现？
2. 如果只转移电荷这一项，是否也需要提前生成.xyz文件？
3. 选择待转移部分的原子时，是否只选择带电基团的原子就可以？比如单体中含有NH3+，那只选择NH3+这4个原子就可以吗？

ggdh

aw10279 发表于 2021-11-25 15:08
钟老师，关于聚合体电荷的问题想请教下：
1. 如果单体本来就带电荷，聚合体中也希望其体现出电荷，是否就 ...

1. 单体带电，如果是用帖子中的组装法，那么组装得到的聚合物的电荷会保留单体片段的电荷。如果是用参数转移法，才需要转移电荷。
2. 只要是参数转移，都需要转移的给体受体的坐标文件。
3. 电荷通常是离域的，如果有NH3+，那么这个正电荷可能分散到周围的很多原子上。所以只选NH3+是不够的。
但是注意，无论哪种方法，要注意片段电荷和整体电荷的正确，比如你的单体带一个正电。那么聚合物中，单体的所有原子的电荷加起来应该是+1，如果你的聚合物有10个单体，那么全部的电荷加起来就是+10

aw10279

ggdh 发表于 2021-11-25 20:41
1. 单体带电，如果是用帖子中的组装法，那么组装得到的聚合物的电荷会保留单体片段的电荷。如果是用参数 ...

谢谢老师的回复！因为看到您的说明讲到目前对单体间的电荷无法中和的处理规则是“将这个净电荷平均分配到每个原子之上。通常每个原子的电荷调整在0.001以下”，所以担心本来带电的单体聚合之后出现带电无法呈现的情况。如果按您讲的组装后聚合物中会保留单体片段的电荷，那只要单体的带电计算正确，直接组装就可以，无需再进行参数转移了是吧？
还有个问题想请教下就是，在得到了聚合体的amber文件prmtop和inpcrd后，该如何将其力场参数导入回tleap，以便在将其溶剂化，或者和其他分子整合后生成整体的拓扑和坐标文件呢？

最后想问下您，如果文章中使用了您的这个工具，该如何进行引用呢，谢谢老师！

ggdh

aw10279 发表于 2021-11-26 09:31
谢谢老师的回复！因为看到您的说明讲到目前对单体间的电荷无法中和的处理规则是“将这个净电荷平均分配到 ...

tleap不支持逆向读取prmtop，我建议用gmx的方式，把top转成itp，然后组合。组合好了以后，可以再转成prmtop。
目前文章还没写（主要是要有功能没完成），可以先引用gitee的网址。

aw10279

ggdh 发表于 2021-11-26 10:47
tleap不支持逆向读取prmtop，我建议用gmx的方式，把top转成itp，然后组合。组合好了以后，可以再转成prmt ...

明白了老师。其实我就是想用gromacs跑MD，但不知道在gromacs里该如何选择水模型对gaff力场的分子进行溶剂化。sob老师回复我说水模型就是针对水分子的力场，那应该就不需要考虑力场一致的问题了吧？另外想确认下：是不是只要单体的电荷计算正确，用组装法得到的聚合体也可以正确反映出单体的带电情况呢？烦请老师指教！

ggdh

aw10279 发表于 2021-11-26 10:58
明白了老师。其实我就是想用gromacs跑MD，但不知道在gromacs里该如何选择水模型对gaff力场的分子进行溶剂 ...

因为单体组装聚合物的时候，需要截掉部分原子，所以可以尝试做resp或者resp2的时候把重复单元的电荷固定为整数

xiaowu759

It seems good. Thank you for sharing.

aw10279

ggdh 发表于 2021-11-26 20:15
因为单体组装聚合物的时候，需要截掉部分原子，所以可以尝试做resp或者resp2的时候把重复单元的电荷固定 ...

感谢老师回复！想问下如果单独用Multiwfn来拟合resp，下一步该如何用您的ztop.py来生成gro和top文件呢？ztop.py调用Multiwfn计算resp的设置是否可以修改，或者自动运行Multiwfn的脚本呢？谢谢老师！

ggdh

aw10279 发表于 2021-11-29 09:08
感谢老师回复！想问下如果单独用Multiwfn来拟合resp，下一步该如何用您的ztop.py来生成gro和top文件呢？z ...

multiwfn 产生好电荷文件abc.chg后，用
ztop.py -p abc.top -x abc.log -r "c=abc.chg"  去载入该电荷文件中的电荷

aw10279

ggdh 发表于 2021-11-29 17:04
multiwfn 产生好电荷文件abc.chg后，用
ztop.py -p abc.top -x abc.log -r "c=abc.chg"  去载入该电荷文 ...

明白了老师，那是否可以用ztop.py -g "abc.log" -r "c=abc.chg" -o abc.top,abc.gro 来一次性生成呢？我是看了您的例4想到是不是有这种偷懒的办法
另外我运行了以下命令，出现报错如下，查看了下源码发现structure类里没有cell_vector属性而只有cell_vect，所以我修改了paramrefine.py里的对应位置，然后运行成功。您可以检查下源码

ggdh

aw10279 发表于 2021-11-30 09:18
明白了老师，那是否可以用ztop.py -g "abc.log" -r "c=abc.chg" -o abc.top,abc.gro 来一次性生成呢？我 ...

对 可以一次性生成。感谢报告bug！