使用ztop拟合二面角参数

ggdh

ztop的安装以及使用ztop产生以及优化小分子力场见：
构建及优化[小大]分子力场参数脚本ztop.py
用ztop制作大分子top文件可以参考：
用ztop构建大分子/聚合物的拓扑文件

背景知识
下面很多东西是基于我自己的理解写的，大家有不同看法欢迎批评讨论。
a)为什么要拟合二面角参数
在力场的各种参数中，VDW，电荷，以及二面角这个几个参数的准确度通常比键长键角的参数重要，因为这些作用比较弱，对应的力常数小，反映在动力学上就是对轨迹的影响更大。其中二面角参数对单分子的构象影响很大，而力场通过原子类型赋予的二面角可能不准，因此在模拟柔性分子的时候，通常需要优化一下二面角以获得更准确的分子构象。
常用的做法是用量化方法去扫描二面角，然后根据量化的结果去拟合二面角的参数，具体逻辑先参考sob在古代写的博文：
1-4非键作用对拟合二面角势能项的重要影响
b)为什么不能从Hessian矩阵直接获得二面角参数
实际上二面角发生变化时，通常对应于这个二面角对应于两个片段的相对旋转，而用Seminario方法获得的二面角参数对应于1,4两个原子的相对旋转，这里的区别是，如果原子旋转，而片段不旋转，就意味着这两个原子会脱离他们本应属于的片段，这也意味着包含这两个原子的一系列键长键角的改变，因此Seminario方法获得力常数会非常的大，通常比正常的二面角力常数大几个数量级。
另外，在动力学中，实际决定一个二面角扭转情况的参数包括1,4原子的电荷，1,4原子的VDW参数，1,4作用的缩放系数(比如Gromacs中的fudgeQQ和fudgeLJ)，以及二面角参数。
因此，在拟合二面角参数之前，一定要先确定好VDW参数，电荷，以及1,4缩放系数。
这也是为什么Hessian矩阵投影方法（Seminario method）原则上不准的原因之一，因为该方法算二面角参数时把电荷，范德华力的作用的统统计算在了二面角参数之内而没有分离出来。但是用这个方法做键长键角参数则没有此问题，因为1,2和1,3的非键作用在动力学中通常被排除掉了(比如gromacs中的nrexcl参数通常=3)。
另外一个原因是二面角的运动范围可能很大，而Seminario所需的Hessian矩阵是在平衡位置获得的，随着二面角的扭转，Hessian矩阵势必也会改变，因此当二面角偏离平衡位置较远时，根据平衡位置Hessian矩阵获得的二面角参数就可能不准了。
c)什么时候二面角参数没那么重要
这里说不重要的意思是可以设为0，或者设一个很小的值，而不是可以乱设一个巨大的值。。。
首先对于环状体系（尤其是刚性共轭环，比如苯环），环内的二面角的的扭转通常会伴随着键角的改变，而键角的力常数通常远大于二面角的力常数，并且环内二面角的扭转通常不影响分子构象。此时二面角参数不是很重要。其次improper二面角也是这种情况，improper二面角的扭转也会带来键角的改变，因此也不怎么重要。
d)二面角参数的含义
如果电荷和范德华作用足够准的话，二面角参数就可以认为是对应于共轭/超共轭作用。而实际上电荷和范德华参数不可能准，部分原因是没有考虑电荷极化效应，vdw力各向异性，原子类型不足等。因此二面角参数还承担了补偿电荷和范德华作用不准的责任。这大概也是1,4缩放系数的意义——减少1,4的电荷和VDW作用，给二面角参数足够的操作空间。
e)影响二面角参数的其他因素
1,n非键作用，其中n>4，有时候因为分子结构的原因，1,n非键作用（比如氢键，位阻）会强烈影响二面角，而1,n作用中由于n>4，没有被缩放。使得二面角的代偿空间变小。这会导致二面角参数拟合的困难，后面会增加Tutorial的去详细讨论这种情况。

类似软件
除了sob介绍的手动方法以外，还有一些软件也能做二面角拟合，我知道的有：
AmberTools里面的mdgx
AmberTools里面的Paramfit
这两个我都曾尝试过，这两个用的方法可能更先进吧？好像是可以做到一次可以拟合多个二面角。
但我个人不是很喜欢用，首先不够懒，步骤n长，还要写各种输入文件，令人头大。
其次方法不够透明，随机产生结构(需要大量的初始结构，否则需要检查是否有随机bias)，一次拟合多个二面角，甚至mdgx还把二面角和键角搞到一起拟合，让人有点无法放心结果。
机器学习里面也有类似的问题，样本量不够大，还弄一个不透明不可解释的模型，搞出来的结果就不可信.
FFTK
这个是VMD的插件，有GUI，好像主要针对charmm力场的，我没用过。
本脚本拟合二面角的逻辑和和sob那篇博文中写的一样，朴素直白，用最小二乘拟合量化和动力学能量的差值，不过加了一点点细节。具体的细节在下面的Tutorial会讲到.

Tutorial
准备工作
ztop的基本介绍和安装方法见：
构建及优化[小大]分子力场参数脚本ztop.py
这里推荐用git下载最新版本的ztop，本文中的例子也会一并下载，一些bug的修复会及时的加入最新版。

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

复制代码

下载后进入coordmagic目录下的ztop_test目录

1. cd coordmagic/ztop_test

复制代码

体系介绍——简单的PI共轭体系
体系如下图，我们需要拟合其中的8-15可旋转键对应的二面角参数，这里我们先过一遍基本流程，后面在介绍一些可变参数。


1）产生top文件
把Gaussian优化得到的log和fchk拷贝到当前目录下，并用它们制作top和gro文件。

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

复制代码

这里用的是默认resp电荷，Gaff2力场，并用log中的结构优化平衡键长键角。

2）产生Gaussian输入文件
这一步产生扫描的Gaussian文件，为了方便并行，每个角度是一个单独的输入文件

1. ztop.py -p PPI.top -x PPI.gro -d 8-15

复制代码

运行这一步后，会在当前目录下生成一个名为PPI_DFit_8-15的目录，里面会有0-360度，每隔5度旋转一次，所以共有73个的单点计算Gaussian输入文件。
这里默认的是刚性扫描（因为速度快），默认的的级别是B3LYP def2SVP em(gd3bj) SCRF(solvent=water)， nprocshared=8， mem=4GB，这些参数都可以在命令行中修改，稍后介绍。

3）提交任务计算
将这些Gaussian任务批量提交到服务器上运算，可以参考：
使用Gaussian时的几个实用脚本和命令
批量提交Gaussian任务至pbs队列
注意最后得到的log文件必须在PPI_DFit_8-15目录中。

4）拟合结果
如果没有运行第三步的，可以先把预先准备好的结果文件拷贝出来用（目前git下载的才有此文件）

1. cp PI_poly_init/PPI_DFit_8-15.xyz ./

复制代码

然后输入和第二步同样的命令进行拟合，

1. ztop.py -p PPI.top -x PPI.gro -d 8-15

复制代码

脚本的逻辑是，先在当前目录下找一个名为PPI_DFit_8-15.xyz的文件，该文件包含了所有的坐标和对应的量化能量，如果找到了就利用该文件中的数据进行拟合。如果该文件找不到，就会寻找一个名为PPI_DFit_8-15的目录，如果没有，就创建该目录和里面的gjf文件，如果有，就在该目录中找正常结束的log文件，如果log文件数量足够(并不需要全部的gjf都正常结束)，就会通过这些log产生PPI_DFit_8-15.xyz，并用其进行拟合。屏幕上会输出：

这里K2的意思是周期为2的力常数，同时会有弹窗显示拟合的曲线图(如果是远程，需要用带xserver的比如xmanager或者xmobaterm)：

可以看到原始的力常数K比较大，导致旋转势垒非常高，这是比较容易出现的状况，因为这里的力常数是根据原子类型来设定的，原子类型肯定都是共轭原子，但是仅通过原子类型是共轭这个信息无法确定两个原子之间到底是单键还是双键，因此这里应该是AmberTools(甩锅动作)把这个键判定为双键了，导致力常数K很大。拟合后的力常数明显变小，势能面也基本符合了。
注意这里其实有4个二面角在同时旋转，这里默认选取第一个二面角来作图。这里ztop会有个内置的打分函数，按照二面角位阻从大到小的顺序排序(不准)。

5）保存结果
确定拟合结果满意后，把拟合结果保存到top文件中

1. ztop.py -p PPI.top -x PPI.gro -d 8-15 -o PPI.top

复制代码

这里用了跟输入文件一样名字的top，其实用不同名字的也行

参数说明
上面例子的命令很简单，用了默认的参数，这里其实有不少参数可以调节，参数的格式为
原子1-原子2;参数1=值;参数2=值......
下面结合上面的这个例子将参数分为两类进行介绍：

产生输入文件用到的参数

1. ztop.py -p PPI.top -x PPI.gro -d "8-15;angle=e5,360,0;opt=no;regen=no"  --gjf "nproc:8;mem:4GB;charge:0;spin:1;method:b3lyp;basis:def2svp;solvent:water"

复制代码

上面是产生Gaussian输入文件用到的参数和对应的默认值
angle=e5,360,0 转8-15这个键0 5 10 15 ...... 360度
opt=no 不做优化算单点，opt=fixone 固定一个二面角做柔性扫描，opt=fixall 固定全部二面角做柔性扫描
regen=no 不重新产生输入文件， regen=gjf 重新产生gaussian输入文件，rengen=xyz 重新读取log文件产生xyz文件
--gjf 是对Gaussian输入文件的参数做一些设定，注意这里用的是冒号":"而不是等号去分隔参数名和值

拟合结果用到的参数

1. ztop.py -p PPI.top -x PPI.gro -d "8-15;period=1,2,3,4;nd=4,3,2,1;eth=100;ar=-180~180"

复制代码

period=1,2,3,4
设定拟合的时候用几个周期，这个是最重要的参数，周期的含义在下面二面角函数形式中有详细介绍，1,2,3,4就是同时使用这4个周期，默认只会用一个自动选择的周期从而避免过拟合，如果效果不好可以考虑增加周期数，但要注意过拟合，下面会详细介绍过拟合的情况。
nd=4,3,2,1
设定怎么把参数分配到二面角上，比如这里8-15旋转键对应了4个二面角，如果nd=4，那么这4个二面角的参数就会一模一样；
如果是nd=1，那么只有第1个二面角会有参数，其他的二面角参数都设为0；
如果是nd=2，那么第1，第2个二面角会有相同参数，其他的二面角参数都设为0；
这里说的第n个二面角是按照拟合后屏幕上输出的二面角的顺序，通常是按位阻从大到小的顺序来。
如果是nd=4,3,2,1则会进行4步拟合，第一步拟合让4个二面角参数一样，第二步用前3个二面角拟合去弥补第一步残留误差，第二步用前2个二面角拟合上一步剩下的误差。。。这样的结果就是这4个二面角的参数大体一样，会有一点小差别，但是会比nd=4更准，因为拟合了4次。
默认值是nd=max...1就是从所有的二面角开始，逐步拟合，每步减少一个二面角，直到第1个。
假如这里我希望2，3号二面角的参数一样（因为它们都是S-C-C-C），那么这里我可以输入nd=4,2.3,1  (注意这里中间是个句号)意思是先4个一起拟合，然后2，3号一起拟合，最后拟合1号。
eth=100
如果某个结构能量的相对量化能量高于100kj/mol，那么放弃这个结构。有时候旋转二面角会产生一些位阻特别大，能量特别高的结构，拟合这些能量没有意义，因此这里做下限制。
需要注意的是限制之后，扭角的分布变窄，样本量变少，更加容易出现过拟合的情况。
ar=-180~180
和eth类似，也是减少扭角的分布，比如某个结构0度和180度附近位阻特别大，那么可以输入ar=-150~-30,30~150

过拟合的问题
拟合二面角非常容易出现过拟合的问题，而判断过拟合的严重程度是看k值的大小。
还是以上面的分子为例，这次我们直接开4个周期进行拟合(上面不加period关键词默认用的是周期2进行拟合，因为旋转键的两边各有2个原子)

1. ztop.py -p PPI.top -x PPI.gro -d "8-15;period=1,2,3,4"

复制代码

结果如下

可以看到比之前拟合的结果更好，但是注意这里的k1和k3都比较大，比原本的k2大，并且符号相反，这就说明k1和k3的效果相互抵消，出现了不合理的大值。虽然在当前曲线下表现完美。但是动力学下的构象空间要大的多。无法确定某个构象下k1和k3不再对抗而是突然联合起来产生一个不合理的结果。这时候我们减少一个周期3再试试：

1. ztop.py -p PPI.top -x PPI.gro -d "8-15;period=1,2,4"

复制代码

结果如下

可以看到此时没有出现很大的k并且拟合结果甚至还有一点点提高，那这种结果就能够接受了。

无法代偿的情况
上面的例子二面角曲线拟合的很完美。但是很多情况下，特别是有1,n位阻作用时，由于位阻容易特别大，和二面角不是一个量级的，用二面角去代偿，就好像以卵击石。几乎不可能拟合的很好，如果硬要强行拟合，会产生很大的K去对抗位阻，而这么大的K在我们扫描的势能面之外会不会干坏事根本无法保证，这时候我们要保持冷静，认清现实，不能一味追求拟合结果，此时我们首先可以通过eth和ar参数去掉哪些能量特别高的点，然后再进行拟合时应该优先确保k值不能太大（我个人认为小于5比较理想），然后确保平衡位置(能量最低点)尽量靠近，而势垒的陡峭程度的优先级放到最低。（后面应该会出一个例子来讲这个情况）

二面角函数形式
之前介绍了二面角参数的意义，结合具体的函数形式去理解二面角参数的含义有利于帮助我们拟合的时候选择合适的参数。
二面角参数有多种函数形式，可以参考Gromacs手册中的说明，而ztop目前只支持周期形式(periodic type)，如下图所示。
其中的α是力常数，也就是上面拟合出来的k，δ是相位，n是周期，通常相位由周期n决定，不需要拟合，因此我们需要获得的是α和n，而二面角参数的特殊之处是n可以同时取多个，类似于傅里叶级数。
这里的n也不能乱取，容易出现过拟合的问题。因此理解不同n的含义，对于帮助我们选n具有一定的意义。
下图中就概括了不同周期下相位的默认值，能量最小和能量最大时二面角的样子，以及对应的化学意义。

n=1时相位为0，这时二面角处于反式共平面能量最低，而顺式共平面能量最高。由于这里假设排除了1,4非键作用，这种势能面可以看出是σ-σ超共轭效应的结果。
n=2时相位为180，这时二面角共面能量最低，而垂直时能量最高，对应于π-π共轭效应。
n=3时相位为0，这时二面角在60 180 300 能量最低，而在0 120 240能量最高，可以对应于经典的乙烷的构象，前者是交叉式，后者是重叠式，这里可以看成是3个超共轭共同作用的结果。
n=4时相位为180，这时二面角在0 90 180能量最低，而在45 135能量最高，化学中没有这样的作用。
上述能量最低角度和能量最高角度可以反过来，只需要把对应的力常数α设为负值就行了。


体系二，带有位阻的共轭体系：
后面的例子慢慢补充

xiaowu759

谢谢分享！

小范范1989

钟老师威武

HZW

钟老师厉害

Entropy.S.I

请问钟老师，需要拟合多个二面角的情况怎么做呢？

ggdh

Entropy.S.I 发表于 2021-8-14 00:35
请问钟老师，需要拟合多个二面角的情况怎么做呢？

一个一个做即可
如果耦合特别强，还没有什么好办法，目前还不支持cmap

函数与激情

十分实用的教程想请教钟老师下得到了小分子拟合的top，我想生成小分子和水盒子的top和inpcrd文件，这是在tleap中实现吗或者有什么其他办法吗，用的是AMBER软件

ggdh

好像不方便，因为tleap好像无法载入prmtop文件

你可以用gromacs 加完水以后 再转成amber（用ztop就能直接转）

函数与激情

ggdh 发表于 2021-11-6 09:24
好像不方便，因为tleap好像无法载入prmtop文件你可以用gromacs 加完水以后 再转成amber（用ztop就能直接转 ...

谢谢钟老师，我来试试哈~

小松鼠嘟嘟

感谢大神分享！非常有用。
如果我想要调用的是OPLS-AA力场，是不是就不可以用ztop来模拟二面角参数了。。

snljty

钟老师好，我想请教一个问题，如果我想模拟某个分子在溶液中的动力学行为，那么我计算模拟使用的原子电荷时，可以考虑RESP2(0.5)。那我在去掉取代基，拟合二面角的时候，应该怎么考虑能量和原子电荷呢？之前和导师讨论的结果，是对二面角做柔性扫描，和对扫描的每个点做高精度单点计算的时候用真空下的结果即可。那就有个问题，拟合二面角我用的原子电荷应该是计算的真空下的原子电荷，还是诸如RESP2(0.5)等方法描述的溶液中的原子电荷呢？谢谢老师~

ggdh

snljty 发表于 2022-4-20 17:42
钟老师好，我想请教一个问题，如果我想模拟某个分子在溶液中的动力学行为，那么我计算模拟使用的原子电荷时 ...

因为最后二面角参数需要和RESP2(0.5)去匹配(共同决定二面角势能面)
所以拟合二面角参数需要用RESP2(0.5)

snljty

ggdh 发表于 2022-4-20 19:15
因为最后二面角参数需要和RESP2(0.5)去匹配(共同决定二面角势能面)
所以拟合二面角参数需要用RESP2(0.5)

明白了，谢谢老师！

Medivan

请问老师，这种情况应该如何选择周期参数呢，period的1234排列组合都试过了，都不能把中间的峰压下去

enthalpy

很好的工具。请教以下，如果我已经得到QM以及MM(待拟合二面角的K设置为0)的能量，是不是也可用手动去调用脚本中函数fit_k_only和 fit_least_square 来拟合得到相关的参数？