关于CHELPG拟合静电势的输入问题

cahndengbin

背景：
我是MD模拟的，最近在利用OPLS_AA的力场参数来模拟碳酸二甲酯DMC，乙烯碳酸酯EC，丙烯碳酸酯PC等几种纯溶剂的物理性质，例如密度、扩散系数、粘度等。
调研文献，发现对于这几个体系研究者大都不直接采用OPLS_AA的原子电荷参数，而倾向选择通过量化的方法拟合气态分子的静电势来获取原子电荷。
如Joural of Molecular Structure:THEOCHEM 806 (2007) 23-34，Structural and thermodynamic properties of liquid ethylene carbonate and propylene carbonate by Monte Carlo Simulations。
我根据文献中的描述也利用Gasuuisan计算EC和PC的原子电荷，但是得到的结果相差较多，所以感到很困惑，烦请大家看看原因何在。

文献中的方法描述：
For consistency with the experimental data, the gas phase geometries of the EC and PC monomers were optimized constraining the ring to a planar conformation. The partial atomic charge distributions which were needed to describe the electrostatic potential in the simulations were computed using the CHELPG procedure. In order to evaluate the effect of electron correlation on charge distributions, the EPS partial charges were computed from both HF and MP2 densities.All the ab initio calculations were carried out with the GAUSSIAN-94 program。

我采用的是Gaussian 09, Revision A.02,PC和EC的输入文件分别如下
%mem=8GB
%nprocshared=24
%lindaworker=cu09
%chk=EC-esp.chk
#p opt=ModRedundant
   hf/6-31g(d)
   gfinput gfprint
   iop(6/42=6)
   pop=chelpg
   scf=(direct,tight)

EC

0 1
O                 -0.01810000    1.11950000   -0.11210000
O                 -0.01830000   -1.11960000    0.11150000
O                  1.98220000   -0.00010000    0.00040000
C                 -1.37108565    0.69073802   -0.07023094
C                 -1.37121561   -0.69094780    0.06774637
C                  0.76030000   -0.00010000   -0.00030000
H                 -2.00030509    1.31867081   -0.70647784
H                 -1.72438657    0.75714750    0.96412714
H                 -2.00110056   -1.31806968    0.70418316
H                 -1.72348115   -0.75860724   -0.96693418

D 3 6 1 4 F
D 3 6 2 5 F

最后两行是为了固定O=C-OS-C2和O=C-OS-C1在同一平面上，以保证环状在平面上。
算完的原子电荷如下
47160 points will be used for fitting atomic charges
Fitting point charges to electrostatic potential
Charges from ESP fit, RMS=   0.00278 RRMS=   0.09423:
Charge=   0.00000 Dipole=     6.0666     0.0000     0.0003 Tot=     6.0666
              1
     1  O   -0.496906
     2  O   -0.496848
     3  O   -0.641308
     4  C    0.240287
     5  C    0.240071
     6  C    1.061391
     7  H    0.023430
     8  H    0.023186
     9  H    0.023458
    10  H    0.023240
但在文献中的结果却是
     1  O   -0.490
     2  O   -0.490
     3  O   -0.610
     4  C    0.235
     5  C    0.235
     6  C    1.000
     7  H    0.030
     8  H    0.030
     9  H    0.030
    10  H   0.030
请教各位，可能的原因是什么？谢谢！

sobereva

差异只有一丁点，很正常
CHELPG结果对于构象很敏感，除非你确保你最终算CHELPG电荷时的结构和文献里绝对一模一样，否则这点差异完全可以接受
而且你的6/42的设定和作者的设定也未必一样。

cahndengbin

sobereva 发表于 2017-3-26 22:42
差异只有一丁点，很正常
CHELPG结果对于构象很敏感，除非你确保你最终算CHELPG电荷时的结构和文献里绝对一 ...

谢谢sobereva老师！
我是最近才开始学习高斯，所以一些想法大都是根据过去做VASP的一些经验。在VASP中，如果起始的晶格对称性一致，优化得到的结构都是一样的。因为文献中提到是固定环状在平面上优化，我通过二面角的固定也使得环状在平面上，然后才优化，所以认为优化的结构应该和文献一致。
IOP的设置(6/42=6)是参考卢天发表的文章《原子电荷计算方法的对比》中计算MK电荷的设置。
我刚入门量化，不太清楚多大的电荷差异是可以接受的？是小数点后一位，还是后两位？例如在我的情况下，第六号C原子和原文中计算的相差0.06，第三号O原子和原文中相差-0.03,其他的原子在小数点后三位。

sobereva

cahndengbin 发表于 2017-3-27 08:52
谢谢sobereva老师！
我是最近才开始学习高斯，所以一些想法大都是根据过去做VASP的一些经验。在VASP中， ...

这种差异完全可以接受，数值设定、程序版本，以及几何结构等方面上，不可能绝对做到和作者完全一样。很多因素都可能导致这种程度的数值差别。
至于你说结构和作者的完全一样，我认为是不可能的

cahndengbin

明白了，谢谢！