计算化学公社

标题: PSI4能量分解结果与高斯计算结果不符 [打印本页]

作者
Author: xiaoshan 时间: 2021-11-23 09:20
标题: PSI4能量分解结果与高斯计算结果不符
各位老师，我通过高斯计算的二体的结合能与能量分解结果里的总能量不符，这样的情况应该怎么办？
我优化的方法和基组是m06-l/aug-cc-pvdz(pp),能量分解的级别是Bronze (sSAPT0/jun-cc-pVDZ)，inp文件里基组使用的是def2-svp,再高级别的就算不动了，不知道是不是基组不合适导致的结果不符
图片左侧是计算的结合能，右侧是能量分解的总能量 (, 下载次数 Times of downloads: 41)

作者
Author: zjxitcc 时间: 2021-11-23 09:29
本帖最后由 zjxitcc 于 2021-11-23 09:31 编辑

“能量分解的级别是Bronze (sSAPT0/jun-cc-pVDZ)，inp文件里基组使用的是def2-svp”
那你到底用啥基组做的能量分解？Bronz指定了jun-cc-pVDZ基组，你却用def2SVP去做？

作者
Author: xiaoshan 时间: 2021-11-23 09:38

zjxitcc 发表于 2021-11-23 09:29
“能量分解的级别是Bronze (sSAPT0/jun-cc-pVDZ)，inp文件里基组使用的是def2-svp”
那你到底用啥基组做的 ...

老师，能量分解生成inp文件的时候就三种级别，我只能选Bronze (sSAPT0/jun-cc-pVDZ)这个级别较低的，但是我的体系里有过渡金属，aug-cc-pvdz是算不了的，所以，生成的inp文件里的基组改成了def2-svp去算。卢天老师的帖子里也有一个类似的算例http://sobereva.com/526里4.3小节

作者
Author: zjxitcc 时间: 2021-11-23 10:06
本帖最后由 zjxitcc 于 2021-11-23 10:08 编辑

xiaoshan 发表于 2021-11-23 09:38
老师，能量分解生成inp文件的时候就三种级别，我只能选Bronze (sSAPT0/jun-cc-pVDZ)这个级别较低的，但是 ...

def2-svp太小了，有更好更合理的选择。建议使用混合基组，即非金属原子用jun-cc-pVDZ基组，过渡金属原子用def2-TZVP，算算看结果如何。

作者
Author: sobereva 时间: 2021-11-23 10:09
“有过渡金属，aug-cc-pvdz是算不了的”明显是错误的
此基组又不是没定义3d金属

要想得到准确答案，就把体系截图、Gaussian输入文件和SAPT的输入文件全都发出来。仔细看
在网上求助计算化学问题的时候必须把问题描述得详细、具体、准确、清楚
http://sobereva.com/620（http://bbs.keinsci.com/thread-25787-1-1.html）

作者
Author: scf 时间: 2021-11-23 12:28
本帖最后由 scf 于 2021-11-23 12:29 编辑

SAPT并不是DFT的超分子计算（算二聚体能量减去单体能量），比如用B3LYP作为单体算SAPT能有色散相互作用，但用B3LYP算二聚体基本是算不出来的，而且SAPT没BSSE。所以并不能保证算出的结合能严格一致（即使积分格点收敛标准基函数数量泛函定义等等都一致）

作者
Author: 413 时间: 2022-2-18 14:57

zjxitcc 发表于 2021-11-23 10:06
def2-svp太小了，有更好更合理的选择。建议使用混合基组，即非金属原子用jun-cc-pVDZ基组，过渡金属原子 ...

请教个做sapt的问题哈。我用高斯pbe0/aug-cc-pvtz级别的算BSSE得到的相互作用能是  86 kcal/mol。
然后用psi4做了sapt为啥跟高斯的结果差这么多呢？
set {
scf_type DF
MAXITER 256
soscf true
guess GWH
freeze_core True
basis aug-cc-pVDZ
}

energy('sapt2+')

结果出来是
  Total HF                      257.95364206 [mEh]    161.86835419 [kcal/mol]    677.25719393 [kJ/mol]
  Total SAPT0                   -314.87458790 [mEh] -197.58678696 [kcal/mol] -826.70311664 [kJ/mol]
  Total SAPT2                   -338.87548447 [mEh] -212.64757693 [kcal/mol] -889.71746189 [kJ/mol]
  Total SAPT2+                -318.83906288 [mEh] -200.07453257 [kcal/mol] -837.11184426 [kJ/mol]

  Special recipe for scaled SAPT0 (see Manual):
Electrostatics sSAPT0    -1222.40785860 [mEh] -767.07251209 [kcal/mol] -3209.43139060 [kJ/mol]
Exchange sSAPT0             2545.77511279 [mEh] 1597.49800138 [kcal/mol] 6683.93163779 [kJ/mol]
Induction sSAPT0          -1549.58832402 [mEh] -972.38135377 [kcal/mol] -4068.44358419 [kJ/mol]
Dispersion sSAPT0          298.16850073 [mEh]    187.10355899 [kcal/mol]    782.84129080 [kJ/mol]
  Total sSAPT0                   71.94743090 [mEh]    45.14769450 [kcal/mol]    188.89795380 [kJ/mol]

作者
Author: zjxitcc 时间: 2022-2-18 15:10

413 发表于 2022-2-18 14:57
请教个做sapt的问题哈。我用高斯pbe0/aug-cc-pvtz级别的算BSSE得到的相互作用能是 86 kcal/mol。
然后 ...

体系有多参考特征么？

作者
Author: 413 时间: 2022-2-18 17:54

zjxitcc 发表于 2022-2-18 15:10
体系有多参考特征么？

苯环的共轭体系和水的吸附

作者
Author: zjxitcc 时间: 2022-2-18 17:56

413 发表于 2022-2-18 17:54
苯环的共轭体系和水的吸附

那你算下看看复合物的RHF/aug-cc-pVDZ波函数是不是稳定的，SAPT2基于RHF来的。

作者
Author: 413 时间: 2022-2-18 18:47
本帖最后由 413 于 2022-2-18 18:49 编辑

zjxitcc 发表于 2022-2-18 17:56
那你算下看看复合物的RHF/aug-cc-pVDZ波函数是不是稳定的，SAPT2基于RHF来的。

请问是用psi4算整个complex的单点吗？

作者
Author: 喵星大佬 时间: 2022-2-18 19:35

413 发表于 2022-2-18 18:47
请问是用psi4算整个complex的单点吗？

都可以，瓜丝也行

作者
Author: 413 时间: 2022-2-18 20:16
本帖最后由 413 于 2022-2-18 20:43 编辑

zjxitcc 发表于 2022-2-18 17:56
那你算下看看复合物的RHF/aug-cc-pVDZ波函数是不是稳定的，SAPT2基于RHF来的。

高斯算出来不稳定？请问在psi4利用sapt做能量分解的时候，哪个命令是优化波函数到stable呢？这样可以吗？
set {
scf_type DF
MAXITER 256
soscf true
guess GWH
freeze_core True
basis aug-cc-pVDZ
reference uhf
STABILITY_ANALYSIS follow
}
energy('sapt2+')
补充：用了上述命令后，报错“!  Only SAPT0 supports a reference different from "reference rhf". !”

求解

Stability analysis using <AA,BB:AA,BB> singles matrix:
***********************************************************************

Eigenvectors of the stability matrix:

Eigenvector 1:    Triplet-A Eigenvalue=-0.0873625  <S**2>=2.000
   69 -> 96       -0.11132
   72 -> 80       -0.13720
   72 -> 84       0.12993
   74 -> 79       0.12584
   74 -> 85       -0.16502
   75 -> 79       -0.27431
   75 -> 85       -0.13659
   75 -> 86       -0.10448
   76 -> 80       0.27812
   76 -> 84       0.16752

Eigenvector 2:    Triplet-A Eigenvalue=-0.0869265  <S**2>=2.000
   69 -> 95       0.11523
   70 -> 79       0.10181
   72 -> 79       0.13602
   72 -> 85       0.11130
   74 -> 80       0.14582
   74 -> 84       0.12963
   75 -> 80       0.26004
   75 -> 84       -0.16997
   76 -> 79       0.25583
   76 -> 85       -0.19565
   77 -> 80       0.11286

Eigenvector 3:    Triplet-A Eigenvalue=-0.0741611  <S**2>=2.000
   73 -> 99       -0.11826
   76 -> 79       -0.15264
   77 -> 80       0.24263
   77 -> 87       -0.11154
   77 -> 89       -0.19715
   78 -> 79       0.32388
   78 -> 86       -0.24102

Eigenvector 4:    Singlet-A Eigenvalue= 0.0362960  <S**2>=0.000
   77 -> 80       0.41242
   77 -> 89       -0.10215
   78 -> 79       0.51662
   78 -> 86       -0.12253

Eigenvector 5:    Singlet-A Eigenvalue= 0.0571764  <S**2>=0.000
   75 -> 79       -0.10487
   76 -> 80       0.23225
   77 -> 79       -0.11852
   78 -> 80       0.61218

Eigenvector 6:    Singlet-A Eigenvalue= 0.0607549  <S**2>=0.000
   76 -> 79       0.18277
   77 -> 80       0.49703
   78 -> 79       -0.39335
SavETr:  write IOETrn= 770 NScale= 10 NData=  16 NLR=1 NState= 6 LETran=    118.
The wavefunction has an RHF -> UHF instability.
Leave Link  914 at Fri Feb 18 20:08:18 2022, MaxMem=  2952790016 cpu:    14797.7

作者
Author: zjxitcc 时间: 2022-2-18 21:10

413 发表于 2022-2-18 20:16
高斯算出来不稳定？请问在psi4利用sapt做能量分解的时候，哪个命令是优化波函数到stable呢？这样可以吗？ ...

这是RHF波函数不稳定的意思。但仅凭这也无法断定该体系的SAPT2计算一定有问题。建议你另开一贴，晒出SAPT输入文件、你的分子结构等信息，供大家详细分析。

PSI4目前的SAPT2（或更高级别）只支持闭壳层计算，所以你没法用UHF参考态，只能用RHF。最低级的SAPT0或sSAPT0虽说支持开壳层计算，但精度往往不够。

我不知道你怎么定义相互作用能的，你说吸附，那SAPT2+的 -200 kcal/mol，与PBE0/aug-cc-pVTZ的86 kcal/mol，哪个是你想要的？负号代表有利于结合？

作者
Author: 413 时间: 2022-2-18 21:16

zjxitcc 发表于 2022-2-18 21:10
这是RHF波函数不稳定的意思。但仅凭这也无法断定该体系的SAPT2计算一定有问题。建议你另开一贴，晒出SAPT ...

PBE0/aug-cc-pVTZ的86 kcal/mol这个结果是想要的，可以确定这个是排斥>吸引，所以这个complex是正的相互作用能
结构没问题，mp2以及不同泛函都能优化出来

输入文件为：
memory 160 GB

molecule dimer {
0 1
C    1.51260000 -1.69570000 -0.51300000
...
C    2.25550000 -0.69800000 0.19480000
--
0 1
O    0.00230000 0.04350000 0.00050000
H    -0.09760000 -0.51740000 -0.81680000
H    0.04370000 -0.53200000 0.81470000
...
}

set {
scf_type DF
MAXITER 256
soscf true
guess GWH
freeze_core True
basis aug-cc-pVDZ
}
energy('sapt2+')

作者
Author: zjxitcc 时间: 2022-2-20 00:08

413 发表于 2022-2-18 21:16
PBE0/aug-cc-pVTZ的86 kcal/mol这个结果是想要的，可以确定这个是排斥>吸引，所以这个complex是正的相互 ...

你可以试试DFT-SAPT或GKS-EDA方法，前者有PSI4, Molpro, Qchem程序支持，后者是XEDA（发邮件问苏培峰老师索取）程序支持。如果结果合理，那你就用DFT-SAPT或GKS-EDA好了。

作者
Author: 413 时间: 2022-2-20 00:26

zjxitcc 发表于 2022-2-20 00:08
你可以试试DFT-SAPT或GKS-EDA方法，前者有PSI4, Molpro, Qchem程序支持，后者是XEDA（发邮件问苏培峰老师 ...

谢谢哈。我在试gamess里的EDA方法，还在反复试设置多少MWORD。有问题请教您哈

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