计算化学公社
标题: 使用PSI4做对称匹配微扰理论(SAPT)能量分解计算 [打印本页]
作者Author: sobereva 时间: 2019-12-24 05:30
标题: 使用PSI4做对称匹配微扰理论(SAPT)能量分解计算
后记：笔者在2023年提出了sobEDAw能量分解方法，用起来非常方便，计算耗时以及对内存和硬盘的要求都远远低于SAPT0，而精度可以达到较高阶SAPT（明显好于SAPT0）的水平，因此强烈建议用sobEDAw代替本文介绍的SAPT！介绍和用法见《使用sobEDA和sobEDAw方法做非常准确、快速、方便、普适的能量分解分析》（http://sobereva.com/685）
使用PSI4做对称匹配微扰理论(SAPT)能量分解计算
Using PSI4 to perform symmetric-adapted perturbation theory (SAPT) energy decomposition calculations
文/Sobereva@北京科音 First release: 2019-Dec-24 Last update: 2020-Oct-15
1 前言
能量分解方法之前我在《Multiwfn支持的弱相互作用的分析方法概览》（http://sobereva.com/252）中简单介绍过。能量分解就是将片段间相互作用能分解成不同物理成分，从而能从能量角度更深入地了解相互作用的本质。其中最流行、被接受程度最高的一种是对称匹配微扰理论（Symmetry-Adapted Perturbation Theory, SAPT），其物理意义较严格，也已经有很多程序支持。经常有人问SAPT怎么做，本文就通过例子简单说一下如何通过免费的PSI4程序来实现。
本文不打算涉及太多SAPT的原理，这部分笔者会在量子化学波函数分析与Multiwfn程序培训班（http://www.keinsci.com/workshop/WFN_content.html）里面专门去详细介绍。如果想自己多了解一些原理的话，可以看以下综述：
Computational Molecular Spectroscopy (2000)第3章
Chem. Rev., 94, 1887 (1994)
WIREs Comput. Mol. Sci., 2, 254 (2012) DOI: 10.1002/wcms.86
WIREs Comput. Mol. Sci., 2, 304 (2019) DOI: 10.1002/wcms.84
WIREs Comput. Mol. Sci., 10, e1452 (2019) DOI: 10.1002/wcms.1452
笔者之前的一些论文使用了SAPT考察了不同问题，很建议一看，可以作为SAPT的应用例子进行参考，了解怎么对SAPT的结果进行分析讨论，也非常欢迎大家引用：
• J. Comput. Chem., 40, 2868 (2019) DOI: 10.1002/jcc.26068。主要内容介绍：《透彻认识氢键本质、简单可靠地估计氢键强度：一篇2019年JCC上的重要研究文章介绍》（http://sobereva.com/513）
• Carbon, 171, 514-523 (2021) DOI: 10.1016/j.carbon.2020.09.048。文章内容详细介绍+大量评注：《全面探究18碳环独特的分子间相互作用与pi-pi堆积特征》（http://sobereva.com/572）
• J. Mol. Model., 19, 5387 (2013) DOI: 10.1007/s00894-013-2034-2。主要内容介绍：《静电效应主导了氢气、氮气二聚体的构型》（http://sobereva.com/209）
• Comput. Theor. Chem., 1195, 113090 (2021) DOI: 10.1016/j.comptc.2020.113090。此文里笔者用SAPT结合def2-TZVPP基组计算了一大批AtXn(X=Cl, Br, I; n=1, 3, 5)∙∙∙H2O/H2S类型的体系，是SAPT研究含重元素体系的很好的实例。
本文只涉及对单个结构下做SAPT计算，如果你要考察SAPT能量项随分子间距离的变化，参看《考察SAPT能量分解的能量项随分子二聚体间距变化的简单方法》（http://sobereva.com/469）。
2 一些SAPT相关的关键性的知识
SAPT在上世纪七十年代最早提出。它既是将片段间作用分解为不同物理成份的方法，本身也是一种计算片段间相互作用能的方法，而且没有BSSE问题。SAPT只能用于研究弱相互作用，不能用来考察化学键这种强相互作用的物理成份，因为SAPT的基本思想是片段间微扰，而强相互作用已经违背了微扰理论的前提。
SAPT可以将相互作用能分解为四部分：
• 静电(electrostatics)：描述片段间经典库仑作用，数值可正可负
• 交换(exchange)：描述片段间近距离出现的交换互斥作用，数值为正（即不利于结合）
• 色散(dispersion)：数值为负，起到吸引作用
• 诱导(induction)：体现片段间电荷相互极化和相互转移的作用，数值为负
如果你对分子间相互作用了解不多的话，建议看此文中的相关介绍：《谈谈“计算时是否需要加DFT-D3色散校正？”》（http://sobereva.com/413）。
SAPT理论涉及片段内的微扰和片段间的微扰，随着微扰阶数增加结果也越来越好。高阶SAPT可以达到接近金标准CCSD(T)的弱相互作用计算精度。原理上来说（按照考虑的微扰阶数逐级增高），精度关系是SAPT0<SAPT2<SAPT2+<SAPT2+(3)<SAPT2+3。SAPT0对中等、中等偏大体系还能算得动，而SAPT2+(3)这样的就只算得动小体系了。
为了让SAPT算的相互作用能精度更好，PSI4程序中的SAPT0还包含了δHF项，它体现了高阶项诱导作用。对于SAPT2+、SAPT2+(3)、SAPT2+3这样的高阶SAPT，还可以加上δMP2项来考虑诱导和色散间的高阶耦合作用，诸如SAPT2+(3)结合δMP2被称为SAPT2+(3)δMP2。但δ项的物理意义不是很清楚，也无法进一步划分。δ项的数值通常较小，一般把它归入诱导项。
在实际中，并非耗时越高的越高阶SAPT结果就越好。在J. Chem. Phys., 140, 094106 (2014)中，通过不同基组结合不同级别SAPT方法的组合测试，发现主要有三个组合最值得推荐，并且用金属贵重程度进行命名：sSAPT0/jun-cc-pVDZ（铜）、SAPT2+/aug-cc-pVDZ（银）、SAPT2+(3)δMP2/aug-cc-pVTZ（金），耗时依次升高，精度也依次升高。除了这些以外的组合就完全没必要考虑了。其中sSAPT0代表scaled SAPT0，是对SAPT0的交换项的经验校正版本，耗时和SAPT0相同，但结合jun-cc-pVDZ时由于误差抵消较好的原因，有明显更好的精度（不要试图将sSAPT0和其它基组结合，比如用更大的aug-cc-pVTZ。因为这样会导致误差抵消得没那么好，结果反倒会整体变差）。在笔者的J. Comput. Chem., 40, 2868 (2019)中的氢键测试中，SAPT2+(3)δMP2/aug-cc-pVTZ计算的弱相互作用能的精度与公认的高精度CCSD(T)/jun-cc-pVTZ结合一半counterpoise校正的结果相比，相对误差只有3 %，对多数体系绝对误差小于0.1 kcal/mol，可见已经非常理想了。
能做SAPT计算的程序不少。PSI4是其中做得最好的，开源免费，速度快，输入文件不复杂，还可以做到很高阶的SAPT2+3。Molpro也能做SAPT，但只支持最低阶的SAPT0，而且价钱很贵，如今还按年收费。SAPT最初提出者Szalewicz等人的SAPT代码虽然功能也挺强，但是安装和使用麻烦。阿Q只能做最低阶的SAPT0，还要钱，没有使用价值。
SAPT还有基于DFT描述单体内作用的变体，包括Molpro能做的DFT-SAPT，以及Szalewicz的SAPT程序和CamCASP程序能做的SAPT(DFT)，它们的思想和结果相仿佛。它们比起原本基于微扰方式描述单体内作用的SAPT在原理上能以更低的耗时达到更高的精度，但麻烦的一点是需要额外提供分子的电离能来校正DFT交换势的渐进行为。前述的J. Mol. Model., 19, 5387 (2013)文中就用DFT-SAPT考察了氢气和氮气二聚体，结果很不错。目前PSI4的SAPT(DFT)搞得还不成熟，预计以后会完善。
PSI4的SAPT0速度很快（远远快于高阶SAPT），但sSAPT0/jun-cc-pVDZ对很大体系照样算不动。这种情况就别指望用基于量子化学的方式做能量分解了，此时最适合的是用Multiwfn做基于力场的能量分解，看《使用Multiwfn做基于分子力场的能量分解分析》（http://sobereva.com/442）。
PSI4的各种SAPT方法中只有SAPT0支持开壳层分子间以及开壳层与闭壳层分子间的计算，但研究开壳层弱相互作用的场合不多，本文不举例子。
有文章专门把SAPT里的极化作用和电荷转移作用能拆分开，以讨论更透彻。在PSI4里也支持，但此时没法用δMP2改进结果，而且需要额外的计算耗时，因此一般不特意考虑这个。
PSI4里SAPT2+及以上的高阶SAPT还可以用CCD方法算色散作用的多体项，需要更多的耗时，据说对于色散作用很强的时候会有改进，本文不做涉及。诸如PSI4里sapt2+(3)(ccd)dmp2关键词就是指SAPT2+(3)δMP2结合这种处理。
PSI4还支持其开发者搞的原子SAPT（A-SAPT）、官能团SAPT（F-SAPT）以及I-SAPT，它们底层都是基于SAPT0的。这些变体不属于本文内容范畴。
其它量子化学程序也支持一些五花八门的能量分解方法，但要么原理上不如SAPT严格（比如LMO-EDA里色散部分只是靠DFT-D色散校正能估算的），要么有这样或那样的缺点（诸如局限性大、昂贵，分解出的物理成分意义不清楚、缺乏主流程序支持等等），对于研究弱相互作用来说都不如SAPT理想。
值得一提的是在J. Comput. Chem., 40, 2868 (2019)中笔者发现对氢键作用体系，基于总的氢键作用能直接就能估计出SAPT各个物理成份的大小，并给出了拟合出来的公式。靠这个都能免得专门去做SAPT分析了，对氢键研究很有实际用处。详见《透彻认识氢键本质、简单可靠地估计氢键强度：一篇2019年JCC上的重要研究文章介绍》（http://sobereva.com/513）。
3 PSI4的安装
PSI4可以在官网http://psicode.org下载。可以直接下载installer包。如果你机子里已经有conda了，也可以按网站上的说明用conda装。也可以下载源代码包自行编译，但没必要这么折腾。
以下是笔者使用PSI4 1.3.2，用installer包的安装过程：
下载Psi4conda-1.3.2-py36-Linux-x86_64.sh后，用chmod +x加上可执行权限，然后运行./Psi4conda-1.3.2-py36-Linux-x86_64.sh启动之。然后输入安装路径，比如/sob/psi4_132，之后开始安装。等安装好后输入yes即可。安装器会在~/.bashrc文件末尾写入一大堆乱七八糟的东西，将之都删掉，然后加上以下两句
export PATH=$PATH:/sob/psi4_132/bin
export PSI_SCRATCH=/sob2/psi4scr
这里/sob2/psi4scr是笔者专门建立的存放PSI4运行过程临时文件用的目录。
之后重新进入终端，PSI4就可以通过输入psi4命令用了。
4 SAPT计算实例
为了令初学者也能很容易地上手SAPT计算，笔者在Multiwfn中加入了相应的选项产生PSI4的SAPT任务的输入文件，只需输入两个片段里原子序号、选择计算级别，就可以马上得到现成的SAPT任务的输入文件。此功能从2019-Dec-23及之后更新的Multiwfn中才有，最新版本可以在官网http://sobereva.com/multiwfn免费下载。对于这个功能，给Multiwfn用的输入文件可以是Multiwfn支持的任何含有结构信息的文件，诸如.xyz、.mol2、.pdb、.fch等等，详见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）。另外，在Multiwfn的主功能0里还提供了方便的获得整个片段里原子序号功能，使得产生SAPT输入文件时输入原子序号非常容易。
本文用到的文件都可以在此下载：http://sobereva.com/attach/526/file.rar。
4.1 例1：水-氨气二聚体
这里首先我们用一个很简单的体系，水-氨气二聚体为例介绍演示SAPT的计算过程。这个体系是笔者的J. Comput. Chem., 40, 2868 (2019)文中计算的第24号体系，当时用B3LYP-D3(BJ)/ma-TZVP级别优化出来的xyz文件是本文文件包里的H2O_NH3.xyz。
启动Multiwfn，然后输入
H2O_NH3.xyz
100 //其它功能 (Part 1)
2 //导出文件或产生量子化学程序输入文件
15 //产生PSI4的输入文件
sapt.inp //产生的文件的名字
1 //SAPT任务
3 //当前体系很小，因此我们选最昂贵的，即“金”级别的SAPT2+(3)δMP2/aug-cc-pVTZ
1-3 //第一个片段（水分子）里原子的序号
4-7 //第二个片段（氨气）里原子的序号
当前目录下此时已经有了sapt.inp文件了，可以关闭Multiwfn窗口了。文件内容如下
memory 6000 MB
molecule dimer {
0 1
O -0.03836200 1.54613400 0.00000000
H 0.06356500 0.57566500 0.00000000
H 0.85354600 1.90026300 0.00000000
--
0 1
N -0.03836200 -1.37627100 0.00000000
H -1.03728700 -1.54652400 0.00000000
H 0.34780000 -1.83229000 0.81779500
H 0.34780000 -1.83229000 -0.81779500
}
set {
scf_type DF
freeze_core True
basis aug-cc-pVTZ
}
energy('sapt2+(3)dmp2')
这个文件里memory应根据自己机子的实际情况进行修改，有更大内存的话可以设得更大。Multiwfn产生的SAPT输入文件里每个片段都假定电荷为0、自旋多重度为1，如果与实际情况不符需自行更改。PSI4默认不冻核，这里用freeze_core True启动冻核近似可以令计算更快。scf_type DF代表利用密度拟合加速SCF过程（其实这行可以不写，因为是默认的）。在官网的某些SAPT例子里还通过df_basis_sapt关键词指定SAPT过程用RI辅助基组，其实是多余的，因为这本来就是默认的。
然后运行psi4 sapt.inp sapt.out -n 36执行这个输入文件，得到sapt.out。-n后面接的是并行核数。在笔者的Intel 2*2696v3双路36核机子下此任务用了24秒就算完了。
从末尾可以看到SAPT结果汇总：
SAPT Results
--------------------------------------------------------------------------------------------------------
Electrostatics -18.65467700 [mEh] -11.70598655 [kcal/mol] -48.97784771 [kJ/mol]
Elst10,r -19.13938171 [mEh] -12.01014335 [kcal/mol] -50.25043976 [kJ/mol]
Elst12,r 0.07639831 [mEh] 0.04794066 [kcal/mol] 0.20058373 [kJ/mol]
Elst13,r 0.40830641 [mEh] 0.25621614 [kcal/mol] 1.07200833 [kJ/mol]
Exchange 20.14110989 [mEh] 12.63873727 [kcal/mol] 52.88047673 [kJ/mol]
Exch10 18.33496859 [mEh] 11.50536649 [kcal/mol] 48.13845340 [kJ/mol]
Exch10(S^2) 18.10056896 [mEh] 11.35827850 [kcal/mol] 47.52303725 [kJ/mol]
Exch11(S^2) 0.43063271 [mEh] 0.27022611 [kcal/mol] 1.13062603 [kJ/mol]
Exch12(S^2) 1.37550859 [mEh] 0.86314467 [kcal/mol] 3.61139731 [kJ/mol]
Induction -6.63886139 [mEh] -4.16594842 [kcal/mol] -17.43032818 [kJ/mol]
Ind20,r -8.79354765 [mEh] -5.51803446 [kcal/mol] -23.08745617 [kJ/mol]
Ind22 -0.83815403 [mEh] -0.52594959 [kcal/mol] -2.20057310 [kJ/mol]
Exch-Ind20,r 5.13721139 [mEh] 3.22364881 [kcal/mol] 13.48774663 [kJ/mol]
Exch-Ind22 0.48965157 [mEh] 0.30726100 [kcal/mol] 1.28558002 [kJ/mol]
delta HF,r (2) -2.69965724 [mEh] -1.69406050 [kcal/mol] -7.08794911 [kJ/mol]
delta MP2,r (2) 0.06563457 [mEh] 0.04118632 [kcal/mol] 0.17232354 [kJ/mol]
Dispersion -5.16221417 [mEh] -3.23933830 [kcal/mol] -13.55339144 [kJ/mol]
Disp20 -5.91459173 [mEh] -3.71146234 [kcal/mol] -15.52875845 [kJ/mol]
Disp30 0.25910602 [mEh] 0.16259148 [kcal/mol] 0.68028275 [kJ/mol]
Disp21 0.10497528 [mEh] 0.06587299 [kcal/mol] 0.27561257 [kJ/mol]
Disp22 (SDQ) 0.06225733 [mEh] 0.03906707 [kcal/mol] 0.16345661 [kJ/mol]
Disp22 (T) -0.79231330 [mEh] -0.49718410 [kcal/mol] -2.08021828 [kJ/mol]
Est. Disp22 (T) -0.94442417 [mEh] -0.59263511 [kcal/mol] -2.47958531 [kJ/mol]
Exch-Disp20 1.27046309 [mEh] 0.79722763 [kcal/mol] 3.33560039 [kJ/mol]
Total HF -7.16040663 [mEh] -4.49322299 [kcal/mol] -18.79964501 [kJ/mol]
Total SAPT0 -11.80453526 [mEh] -7.40745771 [kcal/mol] -30.99280306 [kJ/mol]
Total SAPT2 -10.27049811 [mEh] -6.44483487 [kcal/mol] -26.96518908 [kJ/mol]
Total SAPT2+ -11.04768966 [mEh] -6.93252993 [kcal/mol] -29.00570521 [kJ/mol]
Total SAPT2+(3) -10.38027724 [mEh] -6.51372231 [kcal/mol] -27.25341413 [kJ/mol]
Total SAPT2+dMP2 -10.98205509 [mEh] -6.89134361 [kcal/mol] -28.83338166 [kJ/mol]
Total SAPT2+(3)dMP2 -10.31464267 [mEh] -6.47253599 [kcal/mol] -27.08109059 [kJ/mol]
Special recipe for scaled SAPT0 (see Manual):
Electrostatics sSAPT0 -19.13938171 [mEh] -12.01014335 [kcal/mol] -50.25043976 [kJ/mol]
Exchange sSAPT0 18.33496859 [mEh] 11.50536649 [kcal/mol] 48.13845340 [kJ/mol]
Induction sSAPT0 -6.15381951 [mEh] -3.86158004 [kcal/mol] -16.15685089 [kJ/mol]
Dispersion sSAPT0 -4.59412980 [mEh] -2.88285997 [kcal/mol] -12.06188612 [kJ/mol]
Total sSAPT0 -11.55236243 [mEh] -7.24921687 [kcal/mol] -30.33072337 [kJ/mol]
程序把SAPT2+(3)δMP2定义的所有项都输出了，并且为了考察方便，把不同的项进行加和成为Electrostatics、Exchange、Induction、Dispersion四部分便于考察。其实有些项的归属是有任意性的，诸如Exch-Ind22这种交换和诱导的耦合项归属到交换也可以，归属到诱导也可以，PSI4的开发者将之归到了诱导项里。
从上面的结果我们可知当前级别总相互作用能是-27.08 kJ/mol，即-6.47 kcal/mol，和J. Comput. Chem., 40, 2868 (2019)的表1中相应的数据（BE-2）完全一致，而文中用高精度的CCSD(T)/jun-cc-pVTZ结合一半counterpoise校正能的结果为-6.41 kcal/mol，可见我们当前用的SAPT级别相当精确。总相互作用中静电作用贡献了-48.98 kJ/mol，交换作用贡献了52.88 kJ/mol，诱导作用贡献了-17.43 kJ/mol，色散作用贡献了-13.55 kJ/mol。可见此体系结合的主要贡献者是静电作用，而色散和诱导作用相对次要却也不可完全忽略。
由于当前我们算的是高阶SAPT，因此更低阶的SAPT相互作用能也都顺带给出了。另外还给出了sSAPT0能量，由于这个结果肯定比我们当前的SAPT2+(3)δMP2糙得多，因此没必要管它。
输出文件末尾显示Buy a developer a beer!，说明运行正常结束了。如果开发者让你买咖啡，说明任务失败了。
4.2 例2：尿素晶体中的相邻尿素间的相互作用
此例我们考察尿素晶体中的相邻尿素间的相互作用。之前笔者录过一个视频《基于分子晶体cif文件抠出分子团簇结构》（https://www.bilibili.com/video/av35864488/），按照视频里的做法，基于尿素的cif文件我们可以得到下图的团簇，其结构文件是本文文件包里的Urea_cluster.pdb。
(, 下载次数 Times of downloads: 143)
此例我们想获得中间那个尿素与它上方那个尿素之间的相互作用能和各种能量成份。注意由于X光衍射一般没法得到确切的氢的位置，做SAPT计算之前理应先固定此体系的重原子而优化所有氢原子，做法可参看《在Gaussian中做限制性优化的方法》（http://sobereva.com/404）。但本例仅是示例目的，所有忽略掉了这一步。
由于此例的团簇中尿素很多，分子内原子序号又不连着，所以我们首先应当用Multiwfn把要考察的分子中的原子序号得到。启动Multiwfn，输入Urea_cluster.pdb的路径，然后进入主功能0。在界面上方菜单栏中选Tools - Select fragment，然后输入上图的中间的尿素上随便一个原子序号，比如2，之后会看到整个分子都被高亮显示了，并且在文本框里显示了其中原子的序号，如下所示（如果是Linux版，序号在文本窗口里显示）。将显示的原子序号2,13,16,34,36,55,58,77拷出来备用。
(, 下载次数 Times of downloads: 120)
类似地，再次用上述功能，我们得到团簇最上方的尿素的序号，为3,14,17,35,37,56,59,78。
点击Return按钮返回Multiwfn主菜单，然后输入
pi //进入产生PSI4输入文件的快捷命令，意为PSI4 input
sapt.inp
1 //SAPT任务
1 //作为示例，这里我们用很便宜的sSAPT0/jun-cc-pVDZ级别
2,13,16,34,36,55,58,77 //中间尿素的原子序号
3,14,17,35,37,56,59,78 //上方尿素的原子序号
运行此命令开始计算：psi4 sapt.inp -n 36。用了15秒就算完了。像这样不指定输出文件名的话，程序会将输入文件名自动加上.dat后缀，因此我们得到了sapt.inp.dat。其中SAPT部分的结果如下
Special recipe for scaled SAPT0 (see Manual):
Electrostatics sSAPT0 -24.34792103 [mEh] -15.27855111 [kcal/mol] -63.92545785 [kJ/mol]
Exchange sSAPT0 15.70145759 [mEh] 9.85281339 [kcal/mol] 41.22417122 [kJ/mol]
Induction sSAPT0 -5.93603581 [mEh] -3.72491871 [kcal/mol] -15.58505986 [kJ/mol]
Dispersion sSAPT0 -4.77805738 [mEh] -2.99827627 [kcal/mol] -12.54478793 [kJ/mol]
Total sSAPT0 -19.36055663 [mEh] -12.14893270 [kcal/mol] -50.83113444 [kJ/mol]
即相互作用能为-50.83 kJ/mol。这个结果是比较靠谱的，差不多是两个中等强度氢键的能量。从前面的图来看，中间和上头的尿素间就是形成两个氢键。
注意，上面我们得到的并不能认为是绝对严格的在晶体环境中的中间的尿素和上方的尿素的相互作用能，因为我们忽略了周围尿素对它们电子结构的极化、电荷转移等效应。实际上也没有办法绝对严格地去计算凝聚相环境中的两个分子间的相互作用能，因为体系中存在不可分割的多体耦合作用。这种效应导致比如对于三聚体而言，三聚作用能E(ABC)-E(A)-E(B)-E(C)并不等于每一对二聚体作用能的加和，即[E(AB)-E(A)-E(B)] + [E(AC)-E(A)-E(C)] + [E(BC) - E(B) - E(C)]。尿素分子是显著极性的，但是是电中性的，所以多体耦合作用对相邻尿素之间的作用能的影响虽然不能完全忽略，但也不至于有定性程度的影响。但如果某一个二聚体旁边有一个离子，那么耦合作用就强到不可忽略了，因为它与这两个分子间会有显著的电荷转移，并且产生严重的电子密度的极化，进而明显影响两个单体间的相互作用能。
4.3 例3：碘化氢与水的相互作用
此例我们用SAPT2+(3)δMP2研究下图的碘化氢与水的二聚体。之前笔者用Gaussian通过B3LYP-D3(BJ)结合def-TZVP（氧和氢）和Lanl08(d)（碘）对此体系已做了几何优化，Gaussian输入输出文件已提供在了本文的文件包里。
(, 下载次数 Times of downloads: 123)
用GaussView打开H2O_HI_opt.out，直接保存成optimized.gjf文件。此文件就包含了优化后的结构信息，并且可以直接被Multiwfn读入。
启动Multiwfn，输入以下命令
optimized.gjf
pi
sapt.inp
1 //SAPT任务
3 //SAPT2+(3)δMP2/aug-cc-pVTZ级别
1-3 //水中的原子序号
4,5 //碘化氢中的原子序号
当前的sapt.inp不能像之前那样直接跑，因为aug-cc-pVTZ对第四周期之后的元素都没有定义！对于这种情况，最省事的做法就是把基组改成def2系列。如《谈谈赝势基组的选用》（http://sobereva.com/373）所述，def2对第四周期之后是赝势基组。对于当前体系用def2系列的话，PSI4自动就会对碘用对应的赝势。我们把当前输入文件里的basis aug-cc-pVTZ改成basis def2-QZVP。之所以用def2-QZVP，是因为比它更低的def2-TZVP不带弥散函数，质量比aug-cc-pVTZ差，因此我们得提高zeta数来弥补其不足。QZ档次基组计算弱相互作用精度就已经相当高了，除非是超高精度计算或者是阴离子体系，否则完全可以不带弥散函数（虽说aug-cc-pVTZ-PP对碘有定义，和aug-cc-pVTZ最搭，但在PSI4里没有自带，自定义用着麻烦，所以这里不用）。如果def2-QZVP算不动，改为def2-TZVPP也可以，耗时能低一个数量级甚至更多。
然后我们用PSI4对此体系进行计算，SAPT2+(3)δMP2/def2-QZVP计算的总作用能为-14.08 kJ/mol。静电、交换、诱导、色散分别贡献-31.67、45.99、-15.39、-13.02 kJ/mol。之前笔者用MP2/jun-cc-pVTZ（对碘用jun-cc-pVTZ-PP)结合counterpoise校正，结果是-13.6 kJ/mol，当前SAPT2+(3)δMP2/def2-QZVP的结果与其很接近（注：MP2算弱相互作用精度很平庸，但唯独算氢键很好。当前数据从侧面体现了SAPT2+(3)δMP2/def2-QZVP结果的合理性）。
5 计算CT对结合贡献的方法
如果想把电荷转移(CT)对结合能的贡献单独得到，把SAPT方法名后面加上-ct即可，在Multiwfn产生SAPT输入文件的界面里也可以看到有两个预置的：
10 SAPT0/aug-cc-pVDZ with explicit charge-transfer energy
11 SAPT2+(3)/aug-cc-pVTZ with explicit charge-transfer energy
如果对4.1节的例子选择上面第11号选项的话，就得到了本文文件包里的sapt_CT.inp。运行后输出的CT部分信息如下：
SAPT Charge Transfer Analysis
------------------------------------------------------------------------------------------------
SAPT Induction (Dimer Basis) -4.0048 [mEh] -2.5131 [kcal/mol] -10.5147 [kJ/mol]
SAPT Induction (Monomer Basis) -2.5741 [mEh] -1.6153 [kcal/mol] -6.7584 [kJ/mol]
SAPT Charge Transfer -1.4307 [mEh] -0.8978 [kcal/mol] -3.7563 [kJ/mol]
可见CT部分贡献只有-3.76 kJ/mol，这对应于上面6.7584和10.5147的差值，也即在Monomer basis和Dimer basis下不考虑δHF贡献的Induction项的差值。
当前任务实际上是在Dimer basis和Monomer basis下分别做了一次SAPT2+(3)计算，前者等于单独做SAPT2+(3)给出的结果，其中induction作用能（含δHF在内）为-17.6 kJ/mol。CT贡献的-3.76 kJ/mol只占其很小部分，可见对于弱相互作用二聚体来说，电荷转移作用对结合的贡献通常是非常小的。
6 其它
PSI4里一般自动用SAD方式产生初猜波函数，时候难收敛，此时可以在set { ... }里写上guess GWH换初猜方式，可能就收敛了。不过，有时候GWH会收敛到非基态波函数导致算出的SAPT作用能明显不对。为判断是否有这种可能，可以和Gaussian等程序对比一下HF的计算结果进行判断。另外，在难收敛时也可以尝试加上soscf true和soscf_max_iter 30。
笔者发现在个别时候，用SAPT2+(3)δMP2结合冻核近似时，诱导项里的delta MP2,r (2)项的数量级异常之大，导致SAPT2+(3)δMP2的结合能明显不合理，偏离SAPT2+(3)的很多。这种情况应当去掉freeze_core True避免用冻核（当然，代价是对于含有较多周期较靠后的原子的体系耗时会猛增N倍）。我认为这是PSI4的冻核考虑的bug所致。什么情况下会有这种问题笔者还没摸索出规律，至少之前用SAPT2+(3)δMP2/def2-QZVP计算AtBr5与H2S形成的卤键复合物的时候遇到过这个问题。
有人问做SAPT的时候怎么考虑溶剂效应，实际上SAPT不能直接结合隐式溶剂模型用。溶剂效应对结合能的贡献应当估计为：使用常规的超分子方法（即二聚体减两个单体能量之和）在气相下和溶剂模型下分别计算结合能，二者求差，由此作为溶剂效应对结合能的影响。注意这是一个独立的项，不要纳入到SAPT的任何一个物理成分中。
如果Multiwfn产生SAPT输入文件的功能对你的研究起到了帮助，希望在文章中提及，如The SAPT input files of PSI4 were generated with help of Multiwfn program，并引用Multiwfn程序启动时显示的原文。
作者Author: captain 时间: 2019-12-24 09:39
本帖最后由 captain 于 2019-12-25 10:08 编辑
请问大神两个问题
1 研究两个离子，比如Na+和Cl-，当它们距离比较近的时候，二者间的作用属于强相互作用，此时是否就不能用SAPT计算研究二者间的作用了？
2 第三个例子中用 def2-TZVPP 基组如何？
谢谢指导！
作者Author: mars936 时间: 2019-12-24 19:23
“进入主功能0。在界面上方菜单栏中选Tools - Select fragment”，老师，怎么没看到这个工具？只有Measure geometry和batch plotting orbits
作者Author: wuzhiyi 时间: 2019-12-24 22:37
本帖最后由 wuzhiyi 于 2019-12-24 22:45 编辑
想问一下sapt有没有办法计算两个片段之间的相互作用在第三个片段的影响下会有何变化？
我试着在molecule dimer 里面插入三个片段
得到错误：PsiException: SAPT requires active molecule to have 2 fragments, not 3.
查了查文档说three-body SAPT which is not currently supported by Psi4.
想问一下sob老师有什么推荐的解决方式
作者Author: captain 时间: 2019-12-25 10:18
本帖最后由 captain 于 2019-12-25 10:19 编辑
1 SAPT2008 以后的程序支持Three-body SAPT(DFT) 。
http://www.physics.udel.edu/~sza ... se2.html#x5-40002.1
2 PSI4里的F-SAPT功能，可以将分子划分为若干官能团片段，从而计算任意片段之间的相互作用。
https://pubs.acs.org/doi/abs/10.1021/ct500724p
作者Author: sobereva 时间: 2019-12-25 12:06
出行之前文件把windows版文件传错了，等1月2号回北京后我再上传。可以先用linux版
作者Author: sobereva 时间: 2019-12-25 12:08
1 如果都是离子键范畴的作用就不适合了。而且这时候会伴有强烈的电荷转移，不属于原本SAPT框架下能充分考虑的范围
2 不如def2-QZVP。但如果嫌这个耗时太高，对精度要求又不是特别高的情况，用def2-TZVPP将就一下也说得过去
作者Author: captain 时间: 2019-12-25 12:10
明白了，谢大神指导！
作者Author: wuzhiyi 时间: 2019-12-25 18:43
谢谢
作者Author: hzfish 时间: 2019-12-25 19:19
可以参考如下文献：
Application of Symmetry-Adapted Perturbation Theory to Small Ionic Systems
作者Author: captain 时间: 2019-12-25 19:43
本帖最后由 captain 于 2019-12-25 20:51 编辑
十分感谢！
作者Author: captain 时间: 2019-12-25 20:51
请问大神
看了hzfish推荐的这篇文献 https://pubs.acs.org/doi/10.1021/jp302548u
用SAPT计算了M+X− (M = Li, Na, K, Rb, and X = F, Cl, Br, I)的相互作用能，看起来结果还不错。
是不是意味着SAPT也可以用于这种离子间的强相互作用呢？
作者Author: hzfish 时间: 2019-12-25 23:06
本帖最后由 hzfish 于 2019-12-25 23:29 编辑
F-SAPT功能在psi4中使用fisapt0关键字设置，只能分为2个片段（针对分子间）或3个片段（针对分子内）。
设置4个片段报错如下：
RuntimeError:
Fatal Error: FISAPT: Molecular system must have 2 (A+B) or 3 (A+B+C) fragments
作者Author: sobereva 时间: 2019-12-26 22:07
我暂时不方便看此文章。
如果是以总相互作用能来评判的话，考虑了delta-HF项之后原理上也能用于典型的离子键作用，不过这个项不属于原本意义上的SAPT一部分
作者Author: captain 时间: 2019-12-26 23:30
谢大神指点！
作者Author: captain 时间: 2019-12-29 21:06
还得请问大神
进行Charge-Transfer SAPT扫描计算
最后输出各项能量仿照您的例子
energy('sapt2+(3)-ct')
E_disp = get_variable('SAPT DISP ENERGY') * psi_hartree2kcalmol
E_elst = get_variable('SAPT ELST ENERGY') * psi_hartree2kcalmol
E_exch = get_variable('SAPT EXCH ENERGY') * psi_hartree2kcalmol
E_ind = get_variable('SAPT IND ENERGY') * psi_hartree2kcalmol
E_tot = get_variable('SAPT TOTAL ENERGY') * psi_hartree2kcalmol
E_ct = get_variable("SAPT CT ENERGY") * psi_hartree2kcalmol
psi4.print_out("\n")
psi4.print_out(" Summary of SAPT result (kcal/mol)\n")
psi4.print_out(" Distance E_tot E_elst E_exch E_disp E_ind E_ct\n")
psi4.print_out("%s %6.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f\n" % ("R=",dimerdist,E_tot,E_elst,E_exch,E_disp,E_ind,E_ct))
Summary of SAPT result (kcal/mol)
Distance E_tot E_elst E_exch E_disp E_ind E_ct
R= 2.000 0.000 0.000 0.000 0.000 0.000 239.376
最后输出结果为中，除了E_ct正常，其它能量项都没有正常输出。
是否在进行Charge-Transfer SAPT计算时，有两种能量，一种是SAPT (Dimer Basis)，一种是 SAPT (Monomer Basis)，所以按上面的方式设置输出变量就不合适了？
感谢指点！
作者Author: sobereva 时间: 2019-12-29 23:25
我那个扫描随距离变化的文章的输入文件的写法不适合Charge-Transfer SAPT
作者Author: captain 时间: 2019-12-30 08:27
明白了 谢谢大神！
作者Author: mars936 时间: 2019-12-31 17:50
好的，多谢老师，祝旅途愉快，新年快乐！
作者Author: wuzhiyi 时间: 2020-1-2 03:32
挺期待社长什么时候讲一下cp2k的，特别是关于blyp/Aldrich加周期性边界的
作者Author: 任豹 时间: 2020-1-2 10:53
社长您好，我现在在VASP里中做小的药剂分子在高岭石表面吸附研究，假如两者之间是弱相互作用
因为看得到文献https://doi.org/10.1021/jp412765t研究的是稀有气体和TiO2（周期性晶胞）的相互作用，这个作者把TiO2构建成团簇模型，并使用SAPT（DFT）的方法分解了稀有气体和TiO2的相互作用。作者是在MOLPRO中进行计算的。
想问一下，PSI4是否能实现上述文献里SAPT（DFT）的功能？
假如可以，我是否也可以把高岭石（周期性晶胞）简化成团簇模型，从而研究高岭石团簇和小分子药剂的相互作用分解？
因为现在看到的PSI4的相关文献都是研究两个有机体之间的相互作用，而且您在上面也说了PSI4的SAPT(DFT)搞得还不成熟。
作者Author: sobereva 时间: 2020-1-3 12:03
帖子里说了，PSI4这块还没整好。用常规的SAPT就完了，可以用于这类体系
作者Author: 任豹 时间: 2020-1-3 21:08
好嘞，谢谢社长。刚开始研究SAPT，您给指明了方向。
作者Author: alonewolfyang 时间: 2020-1-9 22:02
请教老师，这程序在M062X/6-311+G**水平能做多大体系的能量分解，我的体系有60个原子，不知耗时多久
作者Author: sobereva 时间: 2020-1-10 01:59
文中说了，PSI4的SAPT(DFT)还没整好
别急着算，先把文章仔细看了，搞清楚SAPT的原理
作者Author: Peter_zhong 时间: 2020-1-10 09:42
老师 新手求助，看了他们的评论 我有点蒙。按照您文中的“3.2 例2：尿素晶体中的相邻尿素间的相互作用”例子，我先对某物质7分子团簇进行M062X/6-311G**水平优化，再进行sSAPT0/jun-cc-pVDZ级别下做SAPT能量分解计算，考察相邻分子间的相互作用，1、这样做是对 还是不对？ 还是说优化和能量分解要用相同的级别？ 2、上面评论的那个同学的意思是 想把sSAPT0/jun-cc-pVDZ级别换成M062X/6-311+G**水平？
作者Author: alonewolfyang 时间: 2020-1-11 08:29
ok，谢谢老师，时刻关注本帖子，我先去把博文和涉及的基本理论学习一下
作者Author: alonewolfyang 时间: 2020-1-11 08:33
我本意是想这做，先用M062X/6-311G** em=gd3下优化出二聚体构型，然后在M062X/6-311+G** em=gd3下做能量分解
作者Author: sobereva 时间: 2020-1-11 20:40
1 没有问题。显然不需要和优化相同级别
作者Author: hunterpyj 时间: 2020-1-15 09:17
PSI4还支持其开发者搞的原子SAPT（A-SAPT）、官能团SAPT（F-SAPT）以及I-SAPT，它们底层都是基于SAPT0的。这些变体不属于本文内容范畴。
sob老师有时间能否简单介绍一下以上几种方法（主要是操作步骤和结果分析），能力有限，自己搞不定
作者Author: sobereva 时间: 2020-1-15 23:34
官网有例子。
短期内没有写的打算。
作者Author: hunterpyj 时间: 2020-1-16 11:20
看到官网的例子了，不知那个.py 咋用
作者Author: Andie 时间: 2020-2-13 16:28
感谢大神的分享，受益匪浅。
参照老师帖子里的介绍，利用PSI4在sapt2+/lanl2dz水平下对惰性气体Xe与MOF中的金属簇进行能量分解分析，得出的相互作用能达到-143KJ/mol.
存在一些疑问
1：sapt2+/aug-cc-PVDZ中由于存在Xe，将基组改为lan2dz,这种做法合理吗？（尝试过def2-QZVP,一直不收敛）
2：弱相互作用能量达到-143KJ/mol，合理吗？请教弱相互作用的结合能通常在什么范围内。
感谢指导！
作者Author: sobereva 时间: 2020-2-13 17:33
用lan2dz这种又烂又不带弥散函数的基组结果属于搞笑
明显不合理，高了一个数量级。建议你先用一般的量子化学方法算算看大概是多少。对这个做能量分解也完全没什么意思，凭直觉就知道主导因素肯定是色散作用。
作者Author: Andie 时间: 2020-2-13 22:00
感谢大神指导，我换成自定义基组，其中对Xe采用aug-cc-PVDZ-PP,其他元素采用aug-cc-PVDZ，计算结果只有-5.2KJ/mol，这个结果又太小了，是不是表明几乎没有作用？
作者Author: sobereva 时间: 2020-2-14 11:32
Xe用DZ级别太小，还不足以得到准确作用能
这种大小再正常不过。区区一个原子，还是色散吸引作用占主导，能有多大？
作者Author: Andie 时间: 2020-2-14 22:27
好的，谢谢大神指导！
作者Author: lich502 时间: 2020-3-12 23:14
本帖最后由 lich502 于 2020-3-16 08:08 编辑
参照老师帖子里的介绍，用老师的例子能算成功，换成自己的体系。出现错误，请问出现这种情况是什么原因？上传了输入和输出文件。麻烦帮忙看看
作者Author: snljty 时间: 2020-3-12 23:25
不贴输入文件谁能猜到...
作者Author: snljty 时间: 2020-3-12 23:27
想请问老师一下，金银铜三个级别在比如四核机器上大概分别能算动多少个原子，含有MP2的方法内存消耗和MP2大概是一个级别么？谢谢
作者Author: sobereva 时间: 2020-3-13 04:51
凡是遇到报错，光贴终端显示的信息毫无用处，必须看输出文件判断
作者Author: sobereva 时间: 2020-3-13 04:52
SAPT的代码和MP2的代码消耗内存没有明显可比性
铜的话算几十个应该也可以，金的话两位数都悬，我没怎么在四核机子上专门测过。
作者Author: snljty 时间: 2020-3-13 08:47
谢谢老师！
作者Author: xiaos 时间: 2020-4-21 11:47
老师您好，对于主客体复合物体系，使用PSI4分解结合自由能时，按以下步骤进行是否可以：
1. 使用G09在M06-2X-D3/6-311G**水平下优化基态分子结构
2. 使用PSI4在SAPT0水平下使用jun-cc-pVDZ基组分解结合自由能
该过程是否合理？
作者Author: sobereva 时间: 2020-4-22 04:22
基本合理
作者Author: xiaos 时间: 2020-4-22 08:13
谢谢老师
作者Author: captain 时间: 2020-6-19 15:28
请问大神，关于文末SAPT Charge Transfer的例子，
文中提到“此任务输出的诱导作用能是-16.04 kJ/mol”，这应该是Monomer Basis SAPT计算里的结果，
我认为严格来说此任务的诱导能应该读取的是Dimer Basis SAPT计算的结果，即-17.60 kJ/mol，虽然二者相差也不大。
作者Author: sobereva 时间: 2020-6-21 03:05
是，笔误已改
作者Author: captain 时间: 2020-6-21 08:18
谢大神指导！
作者Author: captain 时间: 2020-9-18 15:35
请问大神，
使用SAPT-CT方法计算的两个片段的SAPT Charge Transfer 为很小的正值（0.1 kcal/mol），
这是否说明电荷转移作用对于二者的结合起不利作用？
作者Author: sobereva 时间: 2020-9-19 09:12
原理上是如此。
不过数值这么小，算法问题、数值误差之类的因素也可能产生这样数量级的大小，我建议不要过度解释
作者Author: captain 时间: 2020-9-19 09:47
好的，大神，谢谢，
我也怀疑可能是方法的精度问题。
作者Author: xxzj 时间: 2021-6-26 15:29
本帖最后由 xxzj 于 2021-6-26 15:34 编辑
老师，使用PSI4做对称匹配微扰理论(SAPT)能量分解计算有两个问题想请教老师：1. 进行SAPT分解前是需要对分子先优化吗，从晶体中直接选取的结构可以不进行优化直接使用吗？
2. 老师，我想探究的是两个分子间的相互作用，但是每个分子大约80个原子，这个体系算大体系吗，使用该方法是否适合？
作者Author: zjxitcc 时间: 2021-6-26 15:32
如果研究问题就是尽可能与晶体结构接近，模拟晶体环境，那么可以不优化总体的结构（SAPT不要求力为零），但是H原子的位置必须优化。
作者Author: xxzj 时间: 2021-6-26 15:36
知道了，老师，还有一个问题是我探究的是两个分子间相互作用，每个分子大约有80个原子，用这种方法可以计算吗，是否算大体系啦
作者Author: zjxitcc 时间: 2021-6-26 15:50
是大体系，你可以用文中的SAPT0试试，搞台大内存机器。实在算不动的话，你只能换用另一种能量分解方法GKS-EDA了。
作者Author: 喵星大佬 时间: 2021-6-26 18:12
Morokuma类型的能量分解及其衍生品种都是算的动单点就算的动的。
包括不限于LMO-EDA，GKS-EDA，ALMO-EDA
作者Author: xxzj 时间: 2021-6-26 19:15
好的，谢谢老师
作者Author: sobereva 时间: 2021-6-26 22:45
通常需要优化氢原子，看
实验测定分子结构的方法以及将实验结构用于量子化学计算需要注意的问题
http://sobereva.com/569
对于SAPT来说已经很大了，一般也就用scaled SAPT0
作者Author: xxzj 时间: 2021-6-28 15:24
收到，谢谢老师
作者Author: Scarlett-ww 时间: 2021-7-7 12:40
Sob老师，含贵金属的体系，sSAPT0/def2-TZVPP组合合适吗？还是用SAPT0/def2-TZVPP？
作者Author: sobereva 时间: 2021-7-8 12:27
应当说清楚体系具体特征，什么元素、多少原子、两个片段间是什么类型的作用，都得说清楚。描述越含糊，越难以得到准确回复
作者Author: Scarlett-ww 时间: 2021-7-9 14:01
本帖最后由 Scarlett-ww 于 2021-7-9 14:37 编辑
体系是炔基银聚集体，一个炔基，Ag=2-4，银的辅助配体为胺基配体。元素有C,H,N,Ag，主要考察C-Ag和Ag-Ag之间的相互作用，前者属于缺电子的sigma配位或pi配位，后者作用比较弱。因为文章中两种作用不方便用多种方法分析，所以用的PSI4，关于Ag-C键再用CDA进一步分析。也会尝试下用EDA-NOCV方法分析一下。
作者Author: sobereva 时间: 2021-7-10 09:39
SAPT不适合分解化学键程度的作用
作者Author: Scarlett-ww 时间: 2021-7-10 18:48
Sob老师，如果想认真的评估Ag-C的作用类型和Ag-Ag之间的作用，是不是分别用EDA和PSI4理论更严谨一些？还是说，有其它更好的能同时分析这两种键型的方法？
作者Author: sobereva 时间: 2021-7-11 12:53
本来PSI4支持的SAPT就不适合这种化学键程度的作用
EDA没有明确指代，把方法名说清楚
作者Author: zhaoyangwx 时间: 2021-7-17 17:53
本帖最后由 zhaoyangwx 于 2021-7-17 17:56 编辑
Sob老师，我在使用psi4时发现大部分时间只有一个核心满载，一核有难多核围观，只有运行到DF-RHF iter时能偶尔跑满。输出文件看起来线程数设置没问题，请问这是正常现象吗？
作者Author: sobereva 时间: 2021-7-17 18:30
正常
PSI4好多模块并行效率很差。
作者Author: dgwfly 时间: 2021-10-15 21:03
想问下，输入文件里的自旋多重度一般怎么确定呢？我一般设置成1来算，大部分和文献吻合的还可以。看手册SAPT对open-shell的原子的计算收敛度也不是很好。。但一直不太清楚这个量在SAPT里的影响有多大
作者Author: sobereva 时间: 2021-10-16 03:17
显然看你的体系
用什么量子化学程序不得设自旋多重度
作者Author: wuzhiyi 时间: 2021-11-25 00:53
本帖最后由 wuzhiyi 于 2021-11-25 05:19 编辑
感觉现在距离SAPT（DFT）的提出已经有一定时间了。
想问问sob老师觉得算不动sapt2+/aug-cc-pVDZ时，SAPT(DFT：PBE0)/aug-cc-pVDZ或者SAPT(DFT：PBE0)/jun-cc-pVDZ是不是一个好的替代呢？https://psicode.org/psi4manual/master/sapt.html#sapt-dft
或者因为是基于DFT，所以可以砍掉大量的弥散和高阶角动量，所以用SAPT(DFT：PBE0)/def2-DZVP或者SAPT(DFT：PBE0)/def2-TZVP也不错？我看到这里金银铜都带了弥散，是不是算SAPT就需要加弥散？比如用def2-TZVPD？
系统是CNOH带上一两个F，Cl，Br。
psi4这里默认用了PBE0，我别的部分的计算都是用B3LYP的，想问一下，这里PBE0有啥优势么？还是我之间全用B3LYP就行。
作者Author: sobereva 时间: 2021-11-25 10:16
因为目前还没有和金银铜级别的SAPT的对比测试，所以我暂时不纳入推荐。虽然原理上来说，大概率能找到一个SAPT(DFT)级别比SAPT金银铜的性价比更高的。
DFT只是起到描述片段自身的作用，片段间相互作用还是依赖于分子间微扰理论，所以高角动量基函数还是不能少。弥散函数依然得有。
作者Author: wuzhiyi 时间: 2021-11-25 20:49
本帖最后由 wuzhiyi 于 2021-11-25 20:53 编辑
谢谢sob老师的指点，发现好像哪怕是DFT，还是算不懂aug-cc-pVTZ，想问一下，aug-cc-pVDZ和def2-TZVPD之间，如何选择呢？对sapt(DFT)是加更多的弥散比较重要还是加到3 zeta比较重要？还是说这个答案只有做对比测试才可以知道？
系统是CNOH带上一两个F，Cl，Br。
作者Author: sobereva 时间: 2021-11-26 10:58
这个不一定，有理论方法和基组的误差巧合抵消的问题，不是基组越大一定越好，只能测试，比如和不错的级别以超分子方式算的结合能对照来判断。
作者Author: wuzhiyi 时间: 2021-11-26 16:16
我也觉得可能要做对比测试才知道，谢谢sob老师。
作者Author: dgwfly 时间: 2022-1-23 04:33
在计算离子SCN-和分子C54H18间SAPT能量分解出现了一些问题，在短距离（小于0.3nm）时，用bronze standard 计算出现的induction energy很奇怪，甚至出现了正值，可能是 scale的问题，对E_exch-ind 的计算不是很好。有的文献也指出SAPT对含离子的体系，又是强结合的情况计算的不是很好，更高级的SAPT2+(3)可能能得更好的效果。但我的体系太大，不能用更高级别的SAPT去做测试。大家对计算含离子短距离的体系有什么建议吗？
附件是我的输入输出文件
作者Author: sobereva 时间: 2022-1-23 13:25
尝试用SAPT(DFT)，这种大小做得动
作者Author: Lym 时间: 2023-1-12 16:27
您好，请问我想利用SAPT讨论溶剂对探针对碱土金属离子结合的贡献，因此计算探针离子(-1;1）和剩余部分(碱土金属+溶剂: +2;1）间SAPT能量分解。在自己笔记本的虚拟机上进行计算的，计算了一半屏幕出现提示出错，想知道如何解决这个问题，麻烦老师指点一下
(, 下载次数 Times of downloads: 2)
(, 下载次数 Times of downloads: 1)
屏幕提示:
(, 下载次数 Times of downloads: 51)
作者Author: sobereva 时间: 2023-1-13 06:24
可能临时文件把硬盘满了
另外，SAPT0结合def2-TZVPP是很不好的选择，更何况单体还有阴离子的。仔细看下文了解该用什么
使用PSI4做对称匹配微扰理论(SAPT)能量分解计算
http://sobereva.com/526（http://bbs.keinsci.com/thread-15902-1-1.html）
作者Author: 七尺贱 时间: 2023-4-26 13:58
想问一下卢老师，SAPT0做能量分解，SAPT0算不算DFT范畴
作者Author: sobereva 时间: 2023-4-27 02:24
明显不算，毫无关系
SAPT(DFT)或DFT-SAPT才跟DFT有关系
作者Author: 七尺贱 时间: 2023-4-27 12:28
谢谢老师
作者Author: 自由风 时间: 2023-6-29 20:00
我刚安装的Psi4计算时报错（见下图），是我未安装好的问题吗？我用sob老师给的算例是可以完成计算的。
作者Author: sobereva 时间: 2023-6-29 23:32
是
写明了没装libecpint模块
作者Author: 自由风 时间: 2023-6-30 09:13
谢谢回复！我是按照您帖子（http://sobereva.com/526）用installer包安装的，没有单独安装libecpint模块，像这种问题还有补救的措施吗？
作者Author: 自由风 时间: 2023-6-30 13:41
我用您网盘提供的1.5版本，安装之后就可以了
作者Author: 自由风 时间: 2023-7-1 09:37
我用Psi4做能量分解分析报错，找了很久也没找到类似的情况
作者Author: sobereva 时间: 2023-7-2 11:16
甭用M06-2X这种meta系列的泛函
作者Author: DLJS 时间: 2023-8-3 11:38
sob老师，您好，我按照老师的方法做能量分解总是中断，没有明显报错，体系七百多个原子，给了180个G是不是太大了算不了导致的呢？
作者Author: sobereva 时间: 2023-8-3 17:41
强烈建议用sobEDA代替SAPT做大体系能量分解，速度又快（比DFT算单点耗时只高一倍多）又不耗内存和硬盘，结果还非常合理，看原文以及里面的教程：
Tian Lu, Qinxue Chen, A simple, efficient and universal energy decomposition analysis method based on dispersion-corrected density functional theory
ChemRxiv (2023) DOI: 10.26434/chemrxiv-2023-n79rz
作者Author: DLJS 时间: 2023-8-4 15:46
好的，非常感谢老师~
作者Author: 牧生 时间: 2023-9-18 18:25
本帖最后由 牧生 于 2023-9-18 18:44 编辑
我们没有买Gaussian，所以只能用免费的软件。
近日遇到了一个疑问，氯化钠水溶液中，水分子之间可以形成氢键，
(, 下载次数 Times of downloads: 36)
三个氢键的能量分解后的结果分别是
Total SAPT2+(3)dMP2 -6.04419615 [mEh] -3.79279035 [kcal/mol] -15.86903482 [kJ/mol]
Total SAPT2+(3)dMP2 -3.88938541 [mEh] -2.44062619 [kcal/mol] -10.21157999 [kJ/mol]
Total SAPT2+(3)dMP2 -5.04728701 [mEh] -3.16722042 [kcal/mol] -13.25165023 [kJ/mol]
(, 下载次数 Times of downloads: 37)
另一方面，从等值面看起来，带正电的钠离子和水分子中的带负电氧原子之间电荷相反，有相互吸引的趋势，是范德华相互作用。四个区域做能量分解后得到
Total SAPT2+(3)dMP2 -23.36736356 [mEh] -14.66324201 [kcal/mol] -61.35100458 [kJ/mol]
Total SAPT2+(3)dMP2 -27.33041596 [mEh] -17.15009493 [kcal/mol] -71.75599721 [kJ/mol]
Total SAPT2+(3)dMP2 -32.38181808 [mEh] -20.31989762 [kcal/mol] -85.01845164 [kJ/mol]
Total SAPT2+(3)dMP2 -31.68921031 [mEh] -19.88527969 [kcal/mol] -83.20001022 [kJ/mol]
从数值来看，氢键作用弱于范德华相互作用。
疑问：潜意识中，我认为氯化钠水溶液中，水分子间的氢键，应该会高于钠离子和水分子的结合能啊？或者我的出发点就错了，不能用这个方法来判断钠离子与水分子之间的结合能？？
作者Author: sobereva 时间: 2023-9-18 18:52
钠离子的半径很大，因此价层区域电子密度较小，故IGMH等值面中间的颜色不会特别蓝
要搞清楚着色的原理，不要教条地指认
作者Author: Oliviaw 时间: 2024-4-8 14:55
本帖最后由 Oliviaw 于 2024-4-8 17:19 编辑
SAPT理论在70年代没有考虑BSSE，那PSI4有没有把BSSE算进去？我的体系BSSE不可忽略，所以energy decomposition analysis的时候最好是考虑BSSE.
是不是CT项属于energy decomposition analysis的 delta E_int 的诱导项的一部分？
作者Author: sobereva 时间: 2024-4-9 10:20
SAPT原理上不存在BSSE的问题
EDA形式的，比如sobEDAw的能量分解，才可能有BSSE问题，sobEDA/sobEDAw可以在计算时以类似counterpoise的方式考虑BSSE问题，下文说了
使用sobEDA和sobEDAw方法做非常准确、快速、方便、普适的能量分解分析
http://sobereva.com/685（http://bbs.keinsci.com/thread-39446-1-1.html）
不同EDA有不同的定义，诸如CT在sobEDA/sobEDAw里体现在极化项里
作者Author: Oliviaw 时间: 2024-4-9 15:26
本帖最后由 Oliviaw 于 2024-4-9 16:38 编辑
sobEDA/sobEDAw的贴子里只说了可以与g16联用。但是我的体系全部用orca算的，sobEDA/sobEDAw可以和orca联用吗？请问是不是改脚本就可以？脚本改哪几个地方？
贴子里说“在sobEDA或sobEDAw能量分解计算中，如果计算每个单体都使用整个体系里所有原子的基函数，则给出的ΔE_int及其各个能量成份就都相当于是counterpoise方法校正后的。” 请问是怎样操作？
作者Author: sobereva 时间: 2024-4-10 03:13
不支持ORCA
帖子里写得非常清楚
• iCP=0代表不用counterpoise校正，=1代表用（对应基组带w.CB后缀的情况）。分解弱相互作用的时候强烈建议用iCP=1以改进结果，分解化学键作用时则不建议iCP=1，不仅浪费很多时间，还可能令结果更差，见《计算化学键键能时考虑BSSE不仅是多余的甚至是有害的》（http://sobereva.com/381）。
作者Author: yufeng 时间: 2024-10-5 10:43
请问sob老师，PSI4怎么多节点并行呢？
作者Author: sobereva 时间: 2024-10-6 01:29
没法多节点
作者Author: yufeng 时间: 2024-10-7 17:01
请问sob老师，没有高斯能用sobEDA吗？
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |