PSI4能量分解结果与高斯计算结果不符

xiaoshan

各位老师，我通过高斯计算的二体的结合能与能量分解结果里的总能量不符，这样的情况应该怎么办？
我优化的方法和基组是m06-l/aug-cc-pvdz(pp),能量分解的级别是Bronze (sSAPT0/jun-cc-pVDZ)，inp文件里基组使用的是def2-svp,再高级别的就算不动了，不知道是不是基组不合适导致的结果不符
图片左侧是计算的结合能，右侧是能量分解的总能量

zjxitcc

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

xiaoshan

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小节

zjxitcc

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

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

sobereva

“有过渡金属，aug-cc-pvdz是算不了的”明显是错误的
此基组又不是没定义3d金属

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

scf

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

413

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]

zjxitcc

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

体系有多参考特征么？

413

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

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

zjxitcc

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

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

413

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

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

喵星大佬

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

都可以，瓜丝也行

413

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

zjxitcc

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，哪个是你想要的？负号代表有利于结合？

413

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+')