关于GAFF力场参数优化的一些问题

KiritsuguPapa

大家好，sobereva老师好。

我是8月份刚刚入学的直博生，现在在做GAFF的力场参数优化(Force field parameterization)，希望以后在MOF体系中使用。其方法是我刚入学不久就走了的一个博后留下的。也因为我刚来他就走了，组里其他人并不懂这个，我也多少有一点语言障碍，因此想来这里请教一下。

参数优化的基本思想是对QM计算进行拟合，我是使用Gaussian来计算QM和设计的MM模型。
测试用的体系有水分子，苯，甲醇，H2O2，和Paddle Wheel（如图，橙色是Cu）


分析结果的方法是对计算得到的力场参数进行opt/freq MM计算，并和QM计算得到的频率值比对。（最后将所有差值加起来）
我现在的主要困难是在二面角项参数的确定上，由得到的参数进行opt/freq计算得到的频率在低频处（扭转）和QM结果差距很大，相对误差可以达到300％以上（绝对值一般是小几百）
。不同的二面角项参数对计算键角、键长参数的结果影响并不大，所以在高频的地方相对误差可以到10％甚至0.1％以下，还算不错（大概）。

故关于二面角项有一些问题想问问大家和sobereva老师，拟合的基本步骤和sobereva老师的文章《1-4非键作用对拟合二面角势能项的重要影响》差不多，不过我们的方法不是拟合能量曲线而是选了一些别的性质拟合。

• 因为要考虑非键作用，输入文件需要确定各原子的电荷。我看了卢天老师的《原子电荷计算方法的对比》（doi: 10.3866/PKU.WHXB2012281）一文，我想对于我的体系和目的，使用RESP电荷或MK电荷会比较好，不知道对不对？Gaussian能直接得到这两种电荷吗？
• 我注意到卢天老师的上述文章中，是在HF/6-31G*级别下进行的计算（对阴离子为HF/6-31+G*）；同样，博后留下的例子中除了paddle wheel使用B3LYP/6-31G*，其余例子也是使用HF/6-31G*来做QM计算，向之拟合。我是刚刚本科毕业，还没有具体学习量化，只是记得HF不如DFT可靠。为什么对这些计算使用HF而不是用DFT？仅仅是因为计算量的原因吗，但对于这些小分子DFT计算也不太耗时啊？（况且对更大的paddle wheel反倒使用了DFT）
• sobereva老师的文章《Amber力场二面角参数的解释》中提到
  

像这样原子明确的二面角CT-CT-N-C，IDIVF都是1。而X-CT-CT-X这样两端可以是任意原子的，IDIVF是两边可以连接的原子数的乘积，比如这里CT每边都可以任意连3个，所以IDIVF是3*3=9。

对于苯分子，Gaussian中内建的AMBER参数是 * CA CA *，其IDIVF(Gaussian中称作Npaths)为2x2＝4。如果我想要将各二面角分开进行参数化，分为（CA CA CA CA) (CA CA CA HA) (HA CA CA HA)，那么此时IDIVF是否就应该为1？  另外我猜想IDIVF对我求参数并不重要，因为如果IDIVF应该为4，而我设置为1，那么求出的参数是原来的4倍。在使用时也使用IDIVF为1和4倍的参数值就可以了。


最后，因为我还处于“百废待兴”，需要恶补理论知识，冒昧请老师和大家能推荐一点有关这方面的教材、资料或文献，不胜感激！有其他的对新人的宝贵建议也很欢迎：）


一次提了这么多问题真是不好意思，请sobereva老师和大家见谅。若不方便可以仅提示我一个关键词，我再去搜索学习：）

sobereva

看不到图

KiritsuguPapa

sobereva 发表于 2015-11-11 19:59
看不到图

抱歉 之前图放在自己的linode日本服务器 忘了国内是被墙着的； 已经更新

KiritsuguPapa

问题1我差不多明白了。 GAFF参数化适用RESP电荷或AM1-BCC电荷，可以用Gaussian添加Pop=MK iop(6/33=2)得到的输出文件再用antechamber来计算。

另外两个问题还希望有人解答：）

sobereva

1 Gaussian能通过pop=MK直接得到MK电荷，RESP电荷需要Gaussian结合Ambertools里的antechamber。

2019-May-15补充：用ambertools计算RESP电荷的做法已经完全过时，目前算RESP电荷的最理想做法是用Multiwfn程序，极度快速、灵活和方便，见
RESP拟合静电势电荷的原理以及在Multiwfn中的计算
http://sobereva.com/441（http://bbs.keinsci.com/thread-10880-1-1.html）
计算RESP原子电荷的超级懒人脚本（一行命令就算出结果）
http://sobereva.com/476（http://bbs.keinsci.com/thread-12858-1-1.html）

2 HF和DFT计算耗时基本一样，之所以用HF是因为HF没有考虑电子相关作用，这会导致高估键的极性。而水溶剂效应本身也会极化溶质的电荷分布使键的极性增加，因此用HF在气相下拟合RESP电荷，可以等效地反映出水溶剂的这个效应。

3 是

低频误差大可能是你拟合的时候没有考虑二面角的耦合，而只是做了单独拟合。低频模式普遍是由大量二面角耦合而成的。像是paramfit等程序，拟合的时候是把所有要拟合的参数放在一起，取一大堆随机结构来拟合能量差，这样等效地把内坐标之间的耦合考虑进去了。

KiritsuguPapa

sobereva 发表于 2015-11-12 21:27
1 Gaussian能通过pop=MK直接得到MK电荷，RESP电荷需要Gaussian结合Ambertools里的antechamber。

2 HF和D ...

谢谢sob老师的回答。
关于低频误差大，我想是因为方法本身的问题。我们的方法只考虑了平衡构型，并选该构型的其他性质一起做拟合。对键长和键角因为是简单的简谐函数，而且考虑的范围小，还比较准确；在二面角项时因为有非键作用和二面角项的耦合，二面角项的“平衡位置”（cos的最低点）并不在耦合后的平衡构型上，只有一个点就不太好了。
以上，再次谢谢sob老师

smutao

ruixing wang?

KiritsuguPapa

smutao 发表于 2017-11-14 23:45
ruixing wang?

是我

smutao

可以谈谈你最近的jcc上的文章介绍的方法的实际用途吗？我读了 没太懂 谢谢

smutao

http://onlinelibrary.wiley.com/doi/10.1002/jcc.25100/full

KiritsuguPapa

smutao 发表于 2017-11-18 01:39
可以谈谈你最近的jcc上的文章介绍的方法的实际用途吗？我读了 没太懂 谢谢

做MM模拟有时缺少一些参数（含金属的经常缺），这个方法是用于确定键长/键角/二面角参数的。我们拟合的目标是QM计算的Hessian矩阵，也就是找一组参数来最佳的reproduce QM的Hessian。这对于键长和键角项应该都是不错的，但对于二面角应用有限（文章末尾有叙述）。

smutao

http://pubs.acs.org/doi/10.1021/acs.jctc.7b00785
和你的工作有点类似

KiritsuguPapa

smutao 发表于 2017-11-23 11:01
http://pubs.acs.org/doi/10.1021/acs.jctc.7b00785
和你的工作有点类似

嗯 也刚看到，他改进了一下很早的一个Seminario方法（我们文章里和Sem做过对比）。
不过我觉得键长键角项的参数再做努力意义不大。现在做参数的难点主要是二面角和非键项。

smutao

你们老板好像来我们这找工作了

beefly

KiritsuguPapa 发表于 2017-11-20 09:39
做MM模拟有时缺少一些参数（含金属的经常缺），这个方法是用于确定键长/键角/二面角参数的。我们拟合的目 ...

你的键角力常数是如何定义的？力场程序给的公式是d^2E/dA^2，但是一些很老很经典的书上，还要把上面的公式再除以两个相邻边的键长。优点是与键长单位一致，缺点是无法推广到其他的角度，比如二面角、各种内坐标相位角的力常数