sobEDAw计算时非常大的w值导致的色散项结果异常原因分析及解决办法

CrysW555

最近我使用sobEDAw对一种平面共轭的有机小分子二聚体中π–π堆积构型的弱相互作用进行能量分解分析。在对分子间面对面距离进行扫描时，发现sobEDAw在处理面对面距离比较远的构象（约6 Å）时出现异常。具体表现为result.txt中记录的w远超正常范围。例如，输出中出现如图1异常信息：

Note: 3598609999999999712910038715219312640.00% DFT correlation is combined with dispersion correction to yield a SAPT-like dispersion term

根据sobEDAw方法的原理及社长博文中的阐述（http://sobereva.com/685) ,  w应处于[0,1] 区间内。对于我的分子来说，面对面距离比较远时，静电项 ΔE_els为 较小的正值（例如为 0.47 kcal/mol），而色散项 ΔE_disp 为比较大的负值（例如为 –14.63 kcal/mol）时，其比值 ΔE_disp / ΔE_els 为负，远小于参数 a = 0.071。根据 w的计算公式，将导致w是一个非常非常大的数。
在sobEDA脚本中第400到403行给出了解决方法，也就是当w远大于1时，强制设定到1.0。

1.   if [[ $(echo "$scale > 1" | bc) -eq 1 ]] ; then
   
2.     echo | awk '{printf ("%s%6.2f%s\n","Warning: Current w (",v,") is larger than 1.0, scale it to 1.0")}' v=$scale
   
3.     scale=1
   
4.   fi

复制代码

然而实际使用中发现代码在运行这部分时出现语法错误。我使用的计算水平是B3LYP-D3(BJ)/6-31+G(d,p)   ，sobEDA中 iCP=1   从sobEDA脚本中可以读取参数应当为parm_c=0.575;parm_a=0.071;parm_r=2.571。检查Gaussian计算的文件并没有发现任何异常。

我猜测原有脚本中判断 w＞1 时，没有正确解析这种科学计数法（例如 3.5986E+36），导致本应执行的强制缩放逻辑被跳过，w 未被限制在 1.0 ，最终使得色散项计算出现异常。

为了解决这个问题我将这部分修改为

1.   if awk -v val="$scale" 'BEGIN {exit !(val > 1)}'; then
   
2.     echo "Warning: w is too large ($scale), force set to 1.0"
   
3.     scale=1
   
4.   fi

复制代码

重新运行后发现能够正常缩放，结果如图2。

但仔细想了一下，好像改不改都没差别，对于w特别大的情况，sobEDAw的色散项直接就使用上一步得到的E_c = E_DFTc + E_dc就好了....

sobereva

把具体结构文件给出来
难以置信相距6埃时能有“色散项 ΔE_disp 为比较大的负值（例如为 –14.63 kcal/mol）”