镁离子导致了SAPT报错

ZetaFunction

给位老师好，我最近在尝试对一个含镁离子结合二磷酸根的体系做SAPT，由于体系包含原子数量较多而且对精度没有特别高的需求，因此只需要sapt0就足够了。我按照使用PSI4做对称匹配微扰理论(SAPT)能量分解计算 - 思想家公社的门口：量子化学·分子模拟·二次元 (sobereva.com)中的方式选取了“铜”级别精度产生输入文件，但是运行后却遇到报错：

1. Traceback (most recent call last):
   
2.   File "/lustre/home/acct-clsst/clsst/Software/psi4conda/bin/psi4", line 247, in <module>
   
3.     exec(content)
   
4.   File "<string>", line 82, in <module>
   
5.   File "/lustre/home/acct-clsst/clsst/Software/psi4conda/lib//python2.7/site-packages/psi4/driver/driver.py", line 454, in energy
   
6.     wfn = procedures['energy'][lowername](lowername, molecule=molecule, **kwargs)
   
7.   File "/lustre/home/acct-clsst/clsst/Software/psi4conda/lib//python2.7/site-packages/psi4/driver/procrouting/proc.py", line 3263, in run_sapt
   
8.     "RIFIT", core.get_global_option("BASIS"))
   
9.   File "/lustre/home/acct-clsst/clsst/Software/psi4conda/lib//python2.7/site-packages/psi4/driver/p4util/python_helpers.py", line 64, in pybuild_basis
   
10.     key, target, fitrole, other, return_atomlist=return_atomlist)
    
11.   File "/lustre/home/acct-clsst/clsst/Software/psi4conda/lib//python2.7/site-packages/psi4/driver/qcdb/libmintsbasisset.py", line 641, in pyconstruct
    
12.     return_atomlist=return_atomlist)
    
13.   File "/lustre/home/acct-clsst/clsst/Software/psi4conda/lib//python2.7/site-packages/psi4/driver/qcdb/libmintsbasisset.py", line 793, in construct
    
14.     shells, msg = parser.parse(entry, lines)
    
15.   File "/lustre/home/acct-clsst/clsst/Software/psi4conda/lib//python2.7/site-packages/psi4/driver/qcdb/libmintsbasissetparser.py", line 216, in parse
    
16.     raise ValidationError("""Gaussian94BasisSetParser::parse: Unable to match an exponent with one contraction: line %d: %s""" % (lineno, line))
    
17. 
    
18. ValidationError: Gaussian94BasisSetParser::parse: Unable to match an exponent with one contraction: line 721: 217       1.0
    
19. 
    
20. *** Psi4 encountered an error. Buy a developer more coffee!
    
21. *** Resources and help at github.com/psi4/psi4.

复制代码

好像是因为使用的基组不支持镁离子？同一体系用相同方式生成的两个不含镁离子的片段可以正常计算。请问遇到这样的问题要如何解决？由于原子数量比较多，大基组的耗时太大，现考虑将jun-cc-pVDZ基组更换为差不多级别的6-31+G**，请问这样做是否合理？

其实我的体系可以将结合着二磷酸根的镁离子删去，但是遇到了SAPT不收敛的问题，在尝试了更换初猜，添加“soscf true”和“soscf_max_iter 30”关键词等方法后依然没能解决不收敛的问题，考虑可能因为带三个单位负电荷的二磷酸根导致了难收敛，所以才想加上镁离子平衡一下电荷。

ZetaFunction

更新：尝试了6-31+G**，还是报相同的错，可以确定就是镁离子的问题，含有镁离子的体系要怎么设置SAPT输入文件？

sobereva

用sobEDA能量分解就完了，又好又快又方便，完全没必要跟SAPT0较劲
sobEDA原文：
Tian Lu, Qinxue Chen, Simple, Efficient, and Universal Energy Decomposition Analysis Method Based on Dispersion-Corrected Density Functional Theory, J. Phys. Chem. A, 127, 7023 (2023) https://doi.org/10.1021/acs.jpca.3c04374

教程：http://sobereva.com/soft/sobEDA_tutorial.zip

ZetaFunction

sobereva 发表于 2023-8-27 21:37
用sobEDA能量分解就完了，又好又快又方便，完全没必要跟SAPT0较劲
sobEDA原文：
Tian Lu, Qinxue Chen, S ...

感谢sobereva老师回复，我简单看了一下你新发的文章，就测试数据而言和大基组下的SAPT2+符合度很好，看上去确实适合代替SAPT，但是鉴于距离修回期限已经不到一周了，之前使用SAPT计算了不少结构现在就剩下最后几个含镁离子和二磷酸根的结构不收敛得不到结果，所以如果可能我还是希望能够使用SAPT算完全程。要如何设置才能使得含有镁离子的体系可以做SAPT？

如果要使用sobEDA的话我有以下疑问：首先是基组和泛函选择是否需要和几何优化或者单点能计算保持一致？其次是sobEDA是一种新的方法，审稿人可能会不认识，这方面会不会有所影响？

如果不考虑分解直接计算两个分子之间的相互作用，可否简单的将两个片段的单点能相加后减去总的单点能？Houk提出的的DA分解就是这么做的，请问这样的做法有道理吗？

另外我还有一个疑问，如果我像计算的是一个大环分子首尾两端基团之间的弱相互作用，而这个大环分子还没有大到可以忽略把首尾分别截断带来的误差，请问用什么方法比较合适？

ZetaFunction

sobereva 发表于 2023-8-27 21:37
用sobEDA能量分解就完了，又好又快又方便，完全没必要跟SAPT0较劲
sobEDA原文：
Tian Lu, Qinxue Chen, S ...

尝试了一下sobEDA，发现计算速度比预想的要慢不少，请问脚本会自行调用当前可用的所有CPU核心吗？还是说需要我在template.gjf里写上可用核数和内存？

更正：在template.gjf前加上了%nprocshared=40和%mem=160GB，计算速度正常了。

zjxitcc

PSI4 SAPT0可以算，而且不需要换基组、删镁离子，也不需要“尝试了更换初猜，添加'soscf true'和'soscf_max_iter 30'关键词等方法”，也不会碰到SCF不收敛，这些都太耗人力、物力了。写一个最常见的gjf文件，划分好片段，例如

1. %mem=150GB
   
2. %nprocshared=32
   
3. #p HF/jun-cc-pVDZ guess(fragment=2)
   
4. 
   
5. {sapt,bronze}
   
6. 
   
7. 0 1 0 1 0 1
   
8. Mg(fragment=1)   -2.317434   -1.047266    0.166227
   
9. P(fragment=1)   -1.056699    1.019534   -0.011045
   
10. O(fragment=1)    0.398543    1.308235   -0.012175
    
11. O(fragment=1)   -1.924357    2.346760   -0.150724
    
12. H(fragment=1)   -1.331710    3.132213   -0.226824
    
13. O(fragment=1)   -1.575132    0.023841   -1.147757
    
14. O(fragment=1)   -1.617042    0.257333    1.277789
    
15. O(fragment=2)    0.255200    3.986246   -0.309545
    
16. H(fragment=2)    0.536825    4.438329    0.491828
    
17. H(fragment=2)    0.606860    3.067152   -0.234731
    
18. 
    

复制代码

可以用GaussView可视化，检查划分方式。提交

1. frag_guess_wfn MgHPO4_h2o.gjf >MgHPO4_h2o.log 2>&1

复制代码

如果是在集群上，就把上面那句话写到提交任务的脚本里。frag_guess_wfn小程序会调用Gaussian做HF计算，生成PSI4的SAPT0输入文件，传轨道给PSI4，使PSI4各个任务的SCF极速收敛。现在我们有了MgHPO4_h2o.inp文件，不用修改，全都写得好好的，提交

1. psi4 MgHPO4_h2o.inp MgHPO4_h2o.out -n 32

复制代码

搞定。我这个例子结果是

1.   Total HF                        -16.42811437 [mEh]     -10.30879740 [kcal/mol]     -43.13200833 [kJ/mol]
   
2.   Total SAPT0                     -23.83550975 [mEh]     -14.95700818 [kcal/mol]     -62.58012222 [kJ/mol]
   
3. 
   
4.   Special recipe for scaled SAPT0 (see Manual):
   
5.     Electrostatics sSAPT0         -44.88379938 [mEh]     -28.16500933 [kcal/mol]    -117.84239903 [kJ/mol]
   
6.     Exchange sSAPT0                44.45351567 [mEh]      27.89500222 [kcal/mol]     116.71268931 [kJ/mol]
   
7.     Induction sSAPT0              -15.48566734 [mEh]      -9.71740297 [kcal/mol]     -40.65761401 [kJ/mol]
   
8.     Dispersion sSAPT0              -7.31142073 [mEh]      -4.58798577 [kcal/mol]     -19.19613248 [kJ/mol]
   
9.   Total sSAPT0                    -23.22737178 [mEh]     -14.57539584 [kcal/mol]     -60.98345621 [kJ/mol]

复制代码

全程半自动化，不用各种尝试，一次成功。frag_guess_wfn是MOKIT里的一个小程序，不仅支持SAPT0，还支持GKS-EDA等多种能量分解计算自动化、加速，见更多示例。如果你只需要frag_guess_wfn，不需要MOKIT其他功能，直接下载压缩包，解压写好环境变量就可以用了。

sobereva

ZetaFunction 发表于 2023-8-28 14:38
尝试了一下sobEDA，发现计算速度比预想的要慢不少，请问脚本会自行调用当前可用的所有CPU核心吗？还是说 ...

显然Default.Rou里设默认用的核数就完了
sobEDA产生的中间的gjf文件都直接给出了，稍有Gaussian使用常识就知道怎么并行计算

sobereva

ZetaFunction 发表于 2023-8-28 00:41
感谢sobereva老师回复，我简单看了一下你新发的文章，就测试数据而言和大基组下的SAPT2+符合度很好，看上 ...

干嘛跟优化一致，毫无一致的意义。显然能用多好用多好

sobEDA的文章都正式发表了，有什么认识不认识的，倘若不认识就让审稿人认识不就完了

看不懂

ZetaFunction

sobereva 发表于 2023-8-28 15:51
干嘛跟优化一致，毫无一致的意义。显然能用多好用多好

sobEDA的文章都正式发表了，有什么认识不认识的 ...

在template.gjf里设置了核数和内存，速度正常了，得到的结果对比我之前用“银级”SAPT2+得到的结果，分解后的静电交换色散能量的相一致程度很高，差距基本在1kcal/mol左右，但是计算效率提高的实在是太多了，SAPT2+需要算二十个小时才能出的结果居然五分钟就算完了，快到我都以为程序是不是出错了，由此看来SAPT要走入历史了。恭喜sobervera老师，这篇文章未来被引数至少要以千计！

然而我却还有一些疑问：
首先是总相互作用能sobEDA的结果和SAPT2+出入甚大，超过5kcal/mol，分解的各项能量差距却可以忽略，是否是对于总相互作用能的定义有差异？两个片段之间的总相互作用能=片段一单点能+片段二单点能-总单点能，这样的定义是被普遍接受的吗？
其次是SAPT毕竟是用了许多年非常成熟的方法，也是计算片段间弱相互作用的首选，而sobEDA虽然初次体验下来效果很好但是距离被普遍接受还需一段时日，我主要担心审稿人会觉得SAPT用了这么多年可靠性经得起检验为什么要另辟蹊径。不过我倒是很愿意身列第一批引用sobEDA的文章。

sobereva

ZetaFunction 发表于 2023-8-28 16:05
在template.gjf里设置了核数和内存，速度正常了，得到的结果对比我之前用“银级”SAPT2+得到的结果，分解 ...

定义没有差异，只不过计算方法有差异（一个是超分子方式，一个是分子间微扰方式）。sobEDA算的总相互作用能和自己手动做DFT算的总相互作用能是一致的，DFT是什么精度sobEDA就是什么精度。不应认为SAPT2+的总能量精度比DFT更高，因为SAPT2+的精度本来就那么回事，不算很高阶微扰（特别是牵扯电荷转移、共价作用显著时精度更堪忧），而且如sobEDA原文所示，B3LYP-D3(BJ)结合适当基组算的结果和S66测试集的高精度参照数据已经吻合相当好了。

ZetaFunction

谢谢老师解答，“手动做DFT算的总相互作用能”是否是指两个片段分别的单点能之和减去总能量再减去基组重叠误差？这是否是被普遍接受的相互作用能的基本定义？只计算总相互作用能的话DFT计算是否比SAPT在精度和效率上更有优势？
如果是这样那么鉴于实际上相比于分解后的各项能量总相互作用能对于我的文章而言才是更值得关心的，而sobEDA分解后的各项能量和SAPT2+相差很小，再考虑到超高的计算效率，我打算在正在修回即将发表的文章中使用你的方法，如果审稿人问起来我就用上述内容回复。

sobereva

ZetaFunction 发表于 2023-8-28 18:55
谢谢老师解答，“手动做DFT算的总相互作用能”是否是指两个片段分别的单点能之和减去总能量再减去基组重叠 ...

>99%的文章都是这么算的相互作用能
具体看什么泛函和什么档次的SAPT比，没前提没法说。如果是比如SAPT0那种级别最low的级别，一大批色散校正的泛函都能吊打之

之后专门回复我的帖子最好点击下方的回复按钮，否则>95%的概率我注意不到。

2532251334

sobereva 发表于 2023-8-27 21:37
用sobEDA能量分解就完了，又好又快又方便，完全没必要跟SAPT0较劲
sobEDA原文：
Tian Lu, Qinxue Chen, S ...

老师我想问下sobEDA能量分解求电荷转移项有没有具体的教程？

sobereva

2532251334 发表于 2023-12-27 14:19
老师我想问下sobEDA能量分解求电荷转移项有没有具体的教程？

通常没必要单独讨论