计算化学公社

标题: 计算适用于OPLS-AA力场做模拟的1.2*CM5原子电荷的懒人脚本 [打印本页]

作者
Author: sobereva 时间: 2021-1-28 20:41
标题: 计算适用于OPLS-AA力场做模拟的1.2*CM5原子电荷的懒人脚本

计算适用于OPLS-AA力场做模拟的1.2*CM5原子电荷的懒人脚本

A lazy script to calculate 1.2*CM5 atomic charges suitable for OPLS-AA force field simulation

文/Sobereva@北京科音 2021-Jan-28

根据J. Phys. Chem. B., 121, 3864 (2017)中的测试，基于OPLS-AA力场的模拟很适合结合1.2*CM5原子电荷。虽然用1.2*CM5的时候水合自由能的计算精度稍逊于文中测试的另一种原子电荷1.14*CM1A-LBCC（因为这种原子电荷本来就是针对水合自由能训练的校正参数），但在计算其它属性上，如蒸发焓、密度，用1.2*CM5时的误差都更小，而且经验性更低、更有普适性。所以当大家用OPLS-AA力场研究新的小分子的时候，我比较推荐用1.2*CM5电荷。1.2*CM5电荷就是在Truhlar等人提出的原始的CM5电荷基础上乘上1.2，这相当于增大了原子电荷的数量级，等效地体现出了溶剂环境对溶质的极化作用。

笔者之前在《计算RESP原子电荷的超级懒人脚本（一行命令就算出结果）》（http://sobereva.com/476）和《RESP2原子电荷的思想以及在Multiwfn中的计算》（http://sobereva.com/531）中分别给出了超级懒人计算RESP和RESP2电荷的Linux脚本，完全不会用Gaussian的人都能轻松计算，已经有不少人在用。最近碰到完全不会用Gaussian的人试图计算1.2*CM5电荷向我求助，我遂又写了个类似的计算1.2*CM5电荷的懒人Linux脚本，在这里说一下。

这个脚本是Multiwfn文件包里的examples\scripts\1.2CM5.sh，在2021-Jan-28及之后更新的Multiwfn中才有，Multiwfn可以在其主页http://sobereva.com/multiwfn免费下载。

此脚本使用很简单：先确保Gaussian在当前机子里已经装了，见《Gaussian的安装方法及运行时的相关问题》（http://sobereva.com/439），也确保将Multiwfn按照手册2.1.2节的说明在机子里装了。之后假设1.2CM5.sh和一个结构文件phenol.xyz都在当前目录下，只需要运行./1.2CM5.sh phenol.xyz，之后在屏幕上就会看到脚本的运行过程：
Net charge was not defined. Default to 0
Spin multiplicity was not defined. Default to 1
Running optimization task via Gaussian...
Done!
Running formchk...
Running Multiwfn...
Finished! The optimized atomic coordinates with 1.2*CM5 charges (the last column) have been exported to phenol.chg in current folder

最终得到的phenol.chg用文本编辑器打开后可见，其中2、3、4列是优化后的XYZ坐标（埃），最后一列就是1.2*CM5电荷了。

脚本原理：这个脚本会将输入文件里的结构作为初猜结构用Gaussian在B3LYP-D3(BJ)/def2-SVP下做几何优化，然后调用formchk将得到的chk转化为fchk文件，然后调用Multiwfn计算CM5电荷，最后给出1.2*CM5电荷。

几点相关事项：
• 用的输入文件可以是任意含有结构信息而且Multiwfn支持的文件格式，比如pdb/mol/mol2/xyz/fch/gjf/wfn等等等等，见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）。
• 默认情况下电荷当成0，自旋多重度当成1。如果是比如不带电的二重态体系，应当写比如./1.2CM5.sh phenol.xyz 0 2。
• 如果运行脚本时提示没可执行权限，先运行chmod +x ./1.2CM5.sh
• 此脚本调用的是Gaussian 16，如果当前机子里装的是诸如Gaussian 09，需要将文件的第22行改为Gaussian=g09。
• 从脚本里的关键词可见，计算是在真空下进行的，这是刻意而为之，不要自行改成在溶剂模型下进行。
• 如果脚本运行失败，有这么几种可能：(1)Gaussian根本没恰当安装，应当确保能手动调用Gaussian运行一个最简单的计算任务 (2)Multiwfn没按照手册2.1.2节恰当安装。其中有些图形库不方便装的话，可以用Multiwfn的noGUI版，此时不需要装图形库 (3)Gaussian计算失败，可能是版本太老而不支持em=GD3BJ关键词，可将脚本中此关键词去掉。也可能几何优化或SCF不收敛，注意看计算中途产生的gau.out文件内容判断原因。
• 如果你手头已经有Gaussian或其它程序产生的Multiwfn支持的波函数文件了（fch/wfn/wfx/molden/mwfn等），就没必要再用这个脚本了。直接载入Multiwfn，依次输入7、16、1就能得到CM5电荷，再手动乘上1.2即可。
• 前述J. Phys. Chem. B.文章中只是考虑了中性体系，没有考虑带电体系。带电体系用什么原子电荷结合OPLS-AA没有统一说法，一般情况下不能用1.2*CM5电荷，因为此时总电荷都不是整数了。笔者建议此时用Multiwfn算RESP2电荷，虽然与OPLS-AA兼容性没有充分测试，但原理上问题不大。
• 如果使用此脚本计算1.2*CM5电荷，请按照Multiwfn程序包里How to cite Multiwfn.pdf文档的说明恰当引用Multiwfn。

如果你没买Gaussian，也可以用免费的ORCA量子化学程序结合Multiwfn算1.2*CM5原子电荷，笔者也提供了相应的傻瓜式脚本，见《ORCA结合Multiwfn计算RESP、RESP2和1.2*CM5原子电荷的懒人脚本》（http://sobereva.com/637）。

作者
Author: 牧生 时间: 2021-1-28 21:39
本帖最后由牧生于 2021-1-29 07:38 编辑

日期有笔误，应该是
2021-Jan-28及之后更新的Multiwfn中才有。

另有一个疑问，“带电体系用什么原子电荷结合OPLS-AA没有统一说法，一般情况下不能用1.2*CM5电荷，因为此时总电荷都不是整数了。笔者建议此时用Multiwfn算RESP2电荷，虽然与OPLS-AA兼容性没有充分测试，但原理上问题不大。”

我的理解，如果做表面活性剂自组装，比如十二烷基磺酸钠或别的表面活性剂，想用opls力场，此时应该是1.2*CM5电荷也可以（十二烷基磺酸根就用1.2*CM5脚本计算，得到总电荷就是-1.2了，钠离子就用ions自带的，但手动改为电荷+1.2，行吗？？），使用RESP2电荷也可以的？？

如果做不带电的体系油水分离，如甲苯和水，想要用opls力场，1.2*CM5和RESP2哪个会更好呢？

如果我要做聚合物的模拟，如聚丙烯酸钠，要用opls力场，那么，应该用1.2*CM5电荷计算三聚体，然后手动重复建立聚合物链，此时还能不能用RESP2电荷？

作者
Author: sobereva 时间: 2021-1-29 19:34

牧生发表于 2021-1-28 21:39
日期有笔误，应该是
2021-Jan-28及之后更新的Multiwfn中才有。

笔误改了

Na用+1.2明显说不通。用RESP2明显有意义得多

用1.2*CM5，毕竟有JPCB的文献支撑。水的原子电荷不要自己算，直接用水模型里定义的。

用1.2*CM5

作者
Author: Kangtor 时间: 2021-2-20 17:19

牧生发表于 2021-1-28 21:39
日期有笔误，应该是
2021-Jan-28及之后更新的Multiwfn中才有。

我觉得做聚合物用原力场电荷参数就挺好的。

作者
Author: exity 时间: 2021-3-2 08:23
请问楼主和社长研究小分子多聚现象的时候，GAFF小分子+amber99sb-ildn水好？还是OPLS-AA好？

作者
Author: nianbin 时间: 2021-3-2 19:49

exity 发表于 2021-3-2 08:23
请问楼主和社长研究小分子多聚现象的时候，GAFF小分子+amber99sb-ildn水好？还是OPLS-AA好？

GAFF2+amber14 SB ，RESP2电荷

作者
Author: exity 时间: 2021-3-3 08:30

nianbin 发表于 2021-3-2 19:49
GAFF2+amber14 SB ，RESP2电荷

谢谢

作者
Author: 牧生 时间: 2021-3-7 10:23
本帖最后由牧生于 2021-3-7 10:49 编辑

https://doi.org/10.1021/acs.jctc.9b00947

针对带电的SDS，RESP电荷比CM5好一些（文中没有使用1.2*CM5）

作者
Author: sobereva 时间: 2021-3-8 00:39

牧生发表于 2021-3-7 10:23
https://doi.org/10.1021/acs.jctc.9b00947

针对带电的SDS，RESP电荷比CM5好一些（文中没有使用1.2*CM5 ...

那篇文章对拟合静电势电荷的计算方式不好。MP2直接在气相下算没有考虑溶剂的等效极化效应，而在SMD模型下算时，又没有像RESP2(0.5)那样避免高估极化程度。而且体系是阴离子，也没带弥散函数。鉴于原子电荷对结果影响极大，应当测试得更严谨才对。

作者
Author: 牧生 时间: 2021-3-8 07:49
本帖最后由牧生于 2021-3-8 07:51 编辑

sobereva 发表于 2021-3-8 00:39
那篇文章对拟合静电势电荷的计算方式不好。MP2直接在气相下算没有考虑溶剂的等效极化效应，而在SMD模型下 ...

所以，考虑使用RESP2才是更好的选择吧

作者
Author: sobereva 时间: 2021-3-9 04:30

牧生发表于 2021-3-8 07:49
所以，考虑使用RESP2才是更好的选择吧

原理上是

作者
Author: Sasuke 时间: 2021-7-22 22:22
本帖最后由 Sasuke 于 2021-7-22 22:23 编辑

社长我用了这个脚本高斯版本是gaussian09 Multiwfn是目前最新版3.8 然后高斯软件是可以运行的因为成功试算了里面自带的算例.Multiwfn也是按照手册2.1.2装的。然后我用MS画了苯乙烯三聚体导出mol2文件作为这个脚本的输入文件，结果运行报错如图所示。gaussian09的out文件里面显示输入文件有误 A syntax error was detected in the input line. 是高斯版本的问题吗？

作者
Author: sobereva 时间: 2021-7-23 03:41

Sasuke 发表于 2021-7-22 22:22
社长我用了这个脚本高斯版本是gaussian09 Multiwfn是目前最新版3.8 然后高斯软件是可以运行的因为成功试 ...

Gaussian版本太老，不支持DFT-D3。去掉em=GD3BJ

作者
Author: mgqqlwq 时间: 2021-8-13 11:11
请问下社长，我的分子是个聚合物，原子数目有点多1050个，这个用56个核的双路服务器Gaussian还能算得动吗？
我提交了两次似乎都在运行一段时间后报错，错误信息是：

Error in internal coordinate system.
Error termination via Lnk1e in ***/g16/l103.exe at Thu Aug 12 10:50:23 2021.
Job cpu time: 5 days 16 hours 46 minutes 17.4 seconds.
Elapsed time: 0 days 2 hours 27 minutes 48.4 seconds.
File lengths (MBytes): RWF= 56574 Int= 0 D2E= 0 Chk= 2389 Scr= 2

复制代码

请问下这个错误我应该如何解决呢？非常感谢！

作者
Author: sobereva 时间: 2021-8-13 22:48

mgqqlwq 发表于 2021-8-13 11:11
请问下社长，我的分子是个聚合物，原子数目有点多1050个，这个用56个核的双路服务器Gaussian还能算得动吗？ ...

算不动
另外，这种报错可以在脚本里的opt关键词里加入cartesian关键词解决。

聚合物通常都是对重复单元来计算原子电荷，然后拼接成整体，而非直接对整个聚合物直接算。

作者
Author: mgqqlwq 时间: 2021-8-13 23:45

sobereva 发表于 2021-8-13 22:48
算不动
另外，这种报错可以在脚本里的opt关键词里加入cartesian关键词解决。

谢谢社长的回复！
还请问下社长，如果是对重复单元算的话，我可以认为重复单元是电中性的，所以自选多重度的设定也不需要改动对吧？
不过聚合物头部和尾部的单元，会多一个氢原子和羟基，这两个我需要单独算是吗？并且是有净电荷的对吧？请问下社长这个电荷和自选多重度应该怎么设定呢？
如果我前面的理解是对的话，最后聚合物的原子电荷就是头上的单元+N*重复单元+尾上的单元合并起来得到对吗？
谢谢社长！

作者
Author: sobereva 时间: 2021-8-14 04:31

mgqqlwq 发表于 2021-8-13 23:45
谢谢社长的回复！
还请问下社长，如果是对重复单元算的话，我可以认为重复单元是电中性的，所以自选多重 ...

整体净电荷为0

实际计算用聚合度较低的模型来算，比如五聚体，两头的单元的原子电荷作为模拟实际聚合物两头的单元的原子电荷，五聚体最中间那个单元的原子电荷通过简单复制作为实际聚合物中所有不处在两头的单元的原子电荷

作者
Author: mgqqlwq 时间: 2021-8-14 06:46

sobereva 发表于 2021-8-14 04:31
整体净电荷为0

实际计算用聚合度较低的模型来算，比如五聚体，两头的单元的原子电荷作为模拟实际聚合 ...

明白了，非常感谢社长的解答！

作者
Author: naoki 时间: 2021-12-10 20:04

sobereva 发表于 2021-1-29 19:34
笔误改了

Na用+1.2明显说不通。用RESP2明显有意义得多

Sob老师好，我想请教您一下，我想对聚合物结构单元计算1.2*CM5电荷，但是发现Multiwfn主功能7下的子功能-16计算1.2*CM5电荷不像RESP模块可以设置等价性约束，如果想约束的话该如何实现呢，谢谢您！

作者
Author: 喵星大佬 时间: 2021-12-10 22:41

naoki 发表于 2021-12-10 20:04
Sob老师好，我想请教您一下，我想对聚合物结构单元计算1.2*CM5电荷，但是发现Multiwfn主功能7下的子功能- ...

手动平均等价原子

作者
Author: naoki 时间: 2021-12-10 23:42

喵星大佬发表于 2021-12-10 22:41
手动平均等价原子

大佬好，您觉得对于聚合物OPLS-AA算1.2*CM5电荷手动平均等价原子更合理，还是等价约束算RESP2更合理呢，谢谢~

作者
Author: 喵星大佬 时间: 2021-12-11 00:23

naoki 发表于 2021-12-10 23:42
大佬好，您觉得对于聚合物OPLS-AA算1.2*CM5电荷手动平均等价原子更合理，还是等价约束算RESP2更合理呢， ...

我觉得这两个方法算出来的电荷看上去应该区别不大

不过聚合物的话你要怎么算这两种电荷

作者
Author: sobereva 时间: 2021-12-11 04:39

naoki 发表于 2021-12-10 20:04
Sob老师好，我想请教您一下，我想对聚合物结构单元计算1.2*CM5电荷，但是发现Multiwfn主功能7下的子功能- ...

只有拟合静电势电荷在原理上才能实现等价性约束，其它原子电荷都不能，只能之后求平均。
原理上人为平均不如等价性约束来得理想，但差距不大，可以接受

作者
Author: naoki 时间: 2021-12-11 13:48

喵星大佬发表于 2021-12-11 00:23
我觉得这两个方法算出来的电荷看上去应该区别不大

不过聚合物的话你要怎么算这两种电荷

我想按Sob老师在17楼的回复先试试，但是我的体系是2-3官能度教练聚合物，可能搞起来比较麻烦

作者
Author: naoki 时间: 2021-12-11 13:48

sobereva 发表于 2021-12-11 04:39
只有拟合静电势电荷在原理上才能实现等价性约束，其它原子电荷都不能，只能之后求平均。
原理上人为平均 ...

谢谢Sob老师解惑~

作者
Author: 七尺贱 时间: 2023-11-18 10:46

作者
Author: sobereva 时间: 2023-11-19 01:21

七尺贱发表于 2023-11-18 10:46
**** 本内容被作者隐藏 ****

有机体系基本都可以用OPLS-AA，但极少有非用OPLS-AA不可的时候

作者
Author: 七尺贱 时间: 2023-11-20 13:06

sobereva 发表于 2023-11-19 01:21
有机体系基本都可以用OPLS-AA，但极少有非用OPLS-AA不可的时候

谢谢社长

作者
Author: 一条君 时间: 2024-1-21 23:41

sobereva 发表于 2021-1-29 19:34
笔误改了

Na用+1.2明显说不通。用RESP2明显有意义得多

老师，（1）针对【Na用+1.2明显说不通】延展的一个问题，比如NaPF6，阴离子用了0.8*OPLS-2009IL 的电荷，也就是阴离子电荷为-0.8，那么对应的钠离子电荷用+0.8是没有问题的吧？
（2）针对【聚丙烯酸钠，要用opls力场——最后建议用1.2*CM5】延展的一个问题，此时聚丙烯酸阴离子是属于带电体系吧，仍然是建议1.2*CM5吗 {3# 3楼的回答1（Na用+1.2...）和回答3区别在哪里呢），谢谢老师

作者
Author: sobereva 时间: 2024-1-22 02:48

一条君发表于 2024-1-21 23:41
老师，（1）针对【Na用+1.2明显说不通】延展的一个问题，比如NaPF6，阴离子用了0.8*OPLS-2009IL 的电荷， ...

对一个体系，用诸如1.2、-0.8这样的净电荷，除非有文章证明过对此类体系这种情况的电荷能得到合理结果，否则不要随便用。

作者
Author: xishaofan 时间: 2024-11-14 17:34

sobereva 发表于 2021-1-29 19:34
笔误改了

Na用+1.2明显说不通。用RESP2明显有意义得多

sob老师，请问
1）针对带电的十二烷基磺酸根，RESP2电荷计算是更合理，有没有文献支撑呀
2）若十二烷基磺酸根用RESP2电荷计算，其抗衡离子钠离子，可以直接使用opls力场文件里面的ions.itp中的钠离子参数吗？如果不可以，也要用RESP2电荷去计算吗？

作者
Author: sobereva 时间: 2024-11-15 04:53

xishaofan 发表于 2024-11-14 17:34
sob老师，请问
1）针对带电的十二烷基磺酸根，RESP2电荷计算是更合理，有没有文献支撑呀
2）若十二烷基 ...

1 没有。只能告诉你从原理上这是合适的

2 Na+的原子电荷只能为1.0，这是显而易见的，没什么可“计算”的

作者
Author: xishaofan 时间: 2024-11-15 09:56

sobereva 发表于 2024-11-15 04:53
1 没有。只能告诉你从原理上这是合适的

2 Na+的原子电荷只能为1.0，这是显而易见的，没什么可“计算” ...

谢谢sob老师

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