计算化学公社

标题: 求助：过渡金属配合物TDDFT计算出现异常，stable=opt后guess=read计算多重态难以判断 [打印本页]
作者
Author: SkyHou    时间: 2024-12-28 14:12
标题: 求助：过渡金属配合物TDDFT计算出现异常，stable=opt后guess=read计算多重态难以判断
各位老师好！我最近在计算一例Pd(II)配合物的激发态时遇到了问题，想请各位老师帮忙看一看问题在哪，万分感谢！
分子结构是这样的：
(, 下载次数 Times of downloads: 8)
由于这个分子在实验中表现出了光催化的活性，所以我想计算其S1和T1的垂直激发能，从而研究其光催化的性能。
我的操作步骤是这样的：
1：用PBE0-D3(BJ)泛函，Pd和I使用SDD赝势，其它原子全部使用6-31G**做opt结构优化，得到基态的结构后在6-311G**下进行TDDFT计算。输入代码如下：
  1. %chk=C:\Users\houti\Desktop\modify-S0-TDDFT.chk
  2. %mem=32GB
  3. %nprocshared=16
  4. #p genecp pbe1pbe scf=vshift=300 symm=loose TD(nstates=8,50-50)

  6. 240110a_0m

  8. 0 1
  9. ———基态原子坐标———

  11. Pd  I 0
  12. SDD
  13. ****
  14. C H N F 0
  15. 6-311G**
  16. ****

  18. Pd I 0
  19. SDD
复制代码
但是输出文件中出现的S1能量非常低，T1甚至为负，让我觉得很困惑（节选了前10个方便各位老师检查）：
  1. Excitation energies and oscillator strengths:

  3. Excited State   1:      Triplet-AU     -0.0312 eV    -39742.26 nm  f=-0.0000  <S**2>=2.000
  4.      252 -> 253       -3.96485
  5.      252 <- 253        3.90124
  6. This state for optimization and/or second-order correction.
  7. Total Energy, E(TD-HF/TD-DFT) =  -3092.98193931   
  8. Copying the excited state density for this state as the 1-particle RhoCI density.

  10. Excited State   2:      Singlet-AU     0.0400 eV 31008.92 nm  f=0.0021  <S**2>=0.000
  11.      252 -> 253       -4.48062
  12.      252 <- 253        4.42452

  14. Excited State   3:      Triplet-AU     1.1264 eV 1100.69 nm  f=0.0000  <S**2>=2.000
  15.      252 -> 254        0.70688

  17. Excited State   4:      Singlet-AU     1.1283 eV 1098.87 nm  f=0.0012  <S**2>=0.000
  18.      252 -> 254        0.70702

  20. Excited State   5:      Triplet-AU     1.5943 eV  777.68 nm  f=0.0000  <S**2>=2.000
  21.      249 -> 254        0.14934
  22.      251 -> 254        0.68005
  23.      251 <- 254        0.16273

  25. Excited State   6:      Triplet-AU     1.6222 eV  764.29 nm  f=0.0000  <S**2>=2.000
  26.      252 -> 256        0.66542
  27.      252 -> 257        0.20502
  28.      252 -> 259       -0.11013

  30. Excited State   7:      Triplet-AG     1.6223 eV  764.27 nm  f=0.0000  <S**2>=2.000
  31.      252 -> 255        0.66580
  32.      252 -> 258        0.23162

  34. Excited State   8:      Singlet-AG     1.6918 eV  732.86 nm  f=0.0000  <S**2>=0.000
  35.      252 -> 255        0.70612

  37. Excited State   9:      Singlet-AU     1.6918 eV  732.84 nm  f=0.0133  <S**2>=0.000
  38.      252 -> 256        0.70609

  40. Excited State  10:      Triplet-AU     1.8505 eV  670.01 nm  f=0.0000  <S**2>=2.000
  41.      252 -> 253       -0.12344
  42.      252 -> 256       -0.23432
  43.      252 -> 257        0.58959
  44.      252 -> 259       -0.28913
  45.      252 <- 253        0.12425
复制代码
2：搜索之后我听从了sob老师的建议，尝试先对基态单点stable=opt。计算后我在log文件中确定波函数稳定了，便读取.chk，并且再进行单点TDDFT计算：
  1. %oldchk=C:\Users\houti\Desktop\modify-S0-sp-stable.chk
  2. %chk=C:\Users\houti\Desktop\S1.chk
  3. %mem=24GB
  4. %nprocshared=16
  5. #p genecp pbe1pbe scf=vshift=300 TD(nstates=5) guess=read

  7. 240110a_0m

  9. 0 1
  10. ————基态原子坐标————

  12. Pd  I 0
  13. SDD
  14. ****
  15. C H N F 0
  16. 6-311G**
  17. ****

  19. Pd I 0
  20. SDD
复制代码
这次计算也是没有报错，正常结束的，但是输出的激发态的多重态很出乎意料，完全没有找到单重态，反而出现了很多近似三重态的激发态：
  1. Excitation energies and oscillator strengths:

  3. Excited state symmetry could not be determined.
  4. Excited State   1:  3.612-?Sym    1.6372 eV  757.29 nm  f=0.0000  <S**2>=3.011
  5.     249A -> 253A       0.13724
  6.     251A -> 253A       0.67726
  7.     249B -> 253B       0.13716
  8.     251B -> 253B       0.67707
  9.     251A <- 253A       0.15755
  10.     251B <- 253B       0.15764
  11. This state for optimization and/or second-order correction.
  12. Total Energy, E(TD-HF/TD-DFT) =  -3092.98717371   
  13. Copying the excited state density for this state as the 1-particle RhoCI density.

  15. Excited state symmetry could not be determined.
  16. Excited State   2:  2.288-?Sym    2.7595 eV  449.31 nm  f=0.0010  <S**2>=1.058
  17.     252A -> 253A       0.71506
  18.     252B -> 253B      -0.69668

  20. Excited state symmetry could not be determined.
  21. Excited State   3:  2.288-?Sym    2.7598 eV  449.25 nm  f=0.0000  <S**2>=1.058
  22.     252A -> 253A       0.69668
  23.     252B -> 253B       0.71508

  25. Excited state symmetry could not be determined.
  26. Excited State   4:  3.098-?Sym    2.9709 eV  417.33 nm  f=0.0003  <S**2>=2.149
  27.     230A -> 254A       0.11884
  28.     244A -> 254A       0.12286
  29.     249A -> 254A       0.27409
  30.     249A -> 257A       0.11650
  31.     250A -> 254A      -0.55010
  32.     250A -> 256A      -0.20047
  33.     250A -> 262A       0.13064
  34.     251A -> 254A      -0.20505
  35.     249B -> 254B       0.16621
  36.     249B -> 257B      -0.19244
  37.     250B -> 254B      -0.33399
  38.     250B -> 256B      -0.12168
  39.     250B -> 257B      -0.11849
  40.     251B -> 254B      -0.12458

  42. Excited state symmetry could not be determined.
  43. Excited State   5:  3.095-?Sym    2.9710 eV  417.31 nm  f=0.0043  <S**2>=2.145
  44.     249A -> 254A       0.16635
  45.     249A -> 257A      -0.19268
  46.     250A -> 254A      -0.33227
  47.     250A -> 256A      -0.12149
  48.     250A -> 257A      -0.11486
  49.     251A -> 254A      -0.12388
  50.     230B -> 254B       0.11892
  51.     244B -> 254B       0.12295
  52.     249B -> 254B      -0.27404
  53.     249B -> 257B      -0.11739
  54.     250B -> 254B       0.54971
  55.     250B -> 256B       0.19991
  56.     250B -> 262B       0.13067
  57.     251B -> 254B       0.20470
复制代码
想请教一下各位老师，在这种情况下我应该如何操作，来计算得到S1和T1的垂直激发能呢？是因为S1能量很高，所以前五个没有覆盖到吗？想了很久也没想到该怎么解决，谢谢各位老师了！



作者
Author: sobereva    时间: 2024-12-28 14:34
你先看看stable=opt之后波函数变成了什么状态，自旋密度和<S^2>如何。当前肯定已经不是闭壳层单重态了
作者
Author: hebrewsnabla    时间: 2024-12-28 14:34
<S**2> 小于 2的是有自旋污染的单重态。
作者
Author: SkyHou    时间: 2024-12-28 16:06
本帖最后由 SkyHou 于 2024-12-28 16:08 编辑
sobereva 发表于 2024-12-28 14:34
你先看看stable=opt之后波函数变成了什么状态，自旋密度和如何。当前肯定已经不是闭壳层单重态了

感谢sob老师的建议！我去stable=opt的输出文件中找了一下，发现原文是这样的：
  1. SCF Done:  E(UPBE1PBE) =  -3093.04733978     a.u. after    7 cycles
  2.             Convg  =    0.3154D-07                    73 Fock formations.
  3.               S**2 =  1.0110                  -V/T =  2.0643
  4. <Sx>= 0.0000 <Sy>= 0.0000 <Sz>= 0.0000 <S**2>= 1.0110 S= 0.6229
  5. <L.S>= 0.000000000000E+00
  6. Annihilation of the first spin contaminant:
  7. S**2 before annihilation     1.0110,   after     0.0884
  8. Leave Link  508 at Sat Dec 14 21:46:49 2024, MaxMem=  4294967296 cpu:           11368.0 elap:           11368.6
复制代码
<S^2>的值确实很接近0，应该是单重态？
各个原子的自旋密度是这样的：
  1. Mulliken charges and spin densities with hydrogens summed into heavy atoms:
  2.                1          2
  3.      1  Pd   0.195619  -0.072914
  4.      2  I   -0.378884  -0.025316
  5.      3  N   -0.452642   0.001374
  6.      4  N   -0.458608   0.001417
  7.      5  C   -0.088439  -0.001084
  8.      6  C   -0.269191   0.005537
  9.      7  C    0.175458  -0.031568
  10.      8  C    0.085846   0.000554
  11.     10  C   -0.288767  -0.000615
  12.     11  C    0.069140   0.000090
  13.     12  C   -0.035275  -0.000039
  14.     13  C    0.230557  -0.008452
  15.     16  C   -0.057761  -0.000421
  16.     18  C    0.065342   0.000492
  17.     20  C   -0.037184   0.000016
  18.     21  C    0.173507  -0.009933
  19.     24  C    0.184343  -0.001359
  20.     25  N   -0.473817   0.019980
  21.     26  C    0.305820   0.000491
  22.     30  C    0.046250  -0.000031
  23.     32  C    0.204714   0.000138
  24.     33  C    0.055451  -0.000017
  25.     35  C    0.173841  -0.001046
  26.     36  C    0.209608   0.000098
  27.     37  C    0.087455  -0.005048
  28.     41  C    0.308603   0.000163
  29.     45  C    0.136429  -0.875513
  30.     47  C    0.117389   0.003040
  31.     51  C   -0.088439   0.001084
  32.     52  C   -0.269191  -0.005537
  33.     53  C    0.085846  -0.000554
  34.     55  C   -0.288767   0.000615
  35.     56  C    0.069140  -0.000090
  36.     57  C   -0.035275   0.000039
  37.     58  C    0.230557   0.008452
  38.     61  C   -0.057761   0.000421
  39.     63  C    0.065342  -0.000492
  40.     65  C   -0.037184  -0.000016
  41.     66  C    0.173507   0.009933
  42.     69  N   -0.473817  -0.019980
  43.     70  C    0.087455   0.005048
  44.     74  C    0.136429   0.875513
  45.     76  C    0.117389  -0.003040
  46.     80  Pd   0.195619   0.072914
  47.     81  I   -0.378884   0.025316
  48.     82  N   -0.452642  -0.001374
  49.     83  N   -0.458608  -0.001417
  50.     84  C    0.175458   0.031568
  51.     85  C    0.184343   0.001359
  52.     86  C    0.305820  -0.000491
  53.     90  C    0.046250   0.000031
  54.     92  C    0.204714  -0.000138
  55.     93  C    0.055451   0.000017
  56.     95  C    0.173841   0.001046
  57.     96  C    0.209608  -0.000098
  58.     97  C    0.308603  -0.000163
  59.    101  F   -0.207446  -0.000009
  60.    102  F   -0.207485  -0.000002
  61.    103  F   -0.207446   0.000009
  62.    104  F   -0.207485   0.000002
  63.    113  C    0.065804  -0.000005
  64.    114  C   -0.001196   0.000018
  65.    115  C   -0.001432  -0.000013
  66.    116  C    0.066952   0.000042
  67.    117  C    0.066952  -0.000042
  68.    118  C   -0.001432   0.000013
  69.    119  C   -0.001196  -0.000018
  70.    120  C    0.065804   0.000005
复制代码

好像45、74号C原子的自旋密度有点反常，是说明这个所谓稳定的波函数不是闭壳层了吗？所以是不合理的吗？
作者
Author: zjxitcc    时间: 2024-12-28 16:06
无法解决，这是U-TDDFT方法的理论局限性，无法通过“我加一两个fancy的关键词”解决。基态是开壳层单重态时，<S^2>大于0，不是单重态纯态；基于此UHF/UKS波函数做TDDFT计算，每个激发态都有自旋污染，偏离理想值，往往无法分辨激发态所属自旋多重度。同时TD()括号中写singlet/triplet均没用；这也不是某个软件的局限，任何有U-TDDFT方法的软件都有这个问题。这在本论坛上过往“负激发能”的帖子中一般都会顺带提及。解决办法一句话总结：使用更严谨的理论方法，例如MRSF-TDDFT方法（注意，不要使用SF-TDDFT方法），SA-SF-DFT方法，CASSCF-NEVPT2方法等等。
作者
Author: zjxitcc    时间: 2024-12-28 16:10
本帖最后由 zjxitcc 于 2024-12-28 16:16 编辑
SkyHou 发表于 2024-12-28 16:06
感谢sob老师的建议！我去stable=opt的输出文件中找了一下，发现原文是这样的：
的值确实很接近0，应该是 ...

您对输出文件存在错误解读，单重态<S^2>理想值为0，而你这个任务算出来是<S**2>= 1.0110，说明偏离程度很大，体系基态可能存在中等或较强双自由基特征（具体还需结合UNO轨道占据数下定论）。注意，一个out/log文件中可能出现多次<S^2>值，要看最后一次出现的<S^2>，那才是最后一帧结构的 经过稳定性检测的 稳定波函数。

Mulliken charges and spin densities with hydrogens summed into heavy atoms其下第二列数据叫自旋布居，不叫自旋密度，请勿被很多不严谨的文章和人误导。当然，也存在自旋密度一词，那是三维空间中每个位置的未成对电子 的分布。二者均可用于判断体系未成对电子情况。但如果仅仅是判断 是否为开壳层单重态，看<S^2>足矣，不需要看其他量。要分清每个量的用途。
作者
Author: SkyHou    时间: 2024-12-28 16:10
hebrewsnabla 发表于 2024-12-28 14:34
小于 2的是有自旋污染的单重态。

谢谢老师的建议！老师您的意思是：对于<S^2>小于2的态，我可以认为它是单重态激发态的垂直激发能吗？那这个数据有可信度吗？
作者
Author: SkyHou    时间: 2024-12-28 16:14
zjxitcc 发表于 2024-12-28 16:06
无法解决，这是U-TDDFT方法的理论局限性，无法通过“我加一两个fancy的关键词”解决。基态是开壳层单重态时 ...

好的好的！我想请教一下老师，在计算前通过分子结构可以判断出它是闭壳层还是开壳层单重态吗？
作者
Author: SkyHou    时间: 2024-12-28 16:17
zjxitcc 发表于 2024-12-28 16:10
您对输出文件存在错误解读，单重态理想值为0，而你这个任务算出来是= 1.0110，说明偏离程度很大，体系基 ...

对不起老师！是我学艺不精，我看到输出文件中写着：
  1. S**2 before annihilation     1.0110,   after     0.0884
复制代码

我误以为0.0884才是它处理出的<S^2>值，是我的错误！但我对这个物质进行了EPR实验，发现它基态并没有EPR信号，说明其基态应当不是多重态。那这样的计算结果是说明物质很可能是开壳层单重态了吗？谢谢老师！
作者
Author: zjxitcc    时间: 2024-12-28 16:19
SkyHou 发表于 2024-12-28 16:14
好的好的！我想请教一下老师，在计算前通过分子结构可以判断出它是闭壳层还是开壳层单重态吗？

除少数结构奇特的体系（比如一个分子里明显少了2个氢原子，烷烃砍掉2个氢原子；历史上被详细研究过的体系，例如O2的最低单重态），一般无法仅凭肉眼预测体系的双自由基特征强弱。
作者
Author: SkyHou    时间: 2024-12-28 16:26
zjxitcc 发表于 2024-12-28 16:19
除少数结构奇特的体系（比如一个分子里明显少了2个氢原子，烷烃砍掉2个氢原子；历史上被详细研究过的体系 ...

我明白了！太谢谢老师了！耽误您的时间解答我的疑惑真是辛苦您了！
作者
Author: hebrewsnabla    时间: 2024-12-28 16:26
本帖最后由 hebrewsnabla 于 2024-12-28 16:27 编辑
SkyHou 发表于 2024-12-28 16:10
谢谢老师的建议！老师您的意思是：对于小于2的态，我可以认为它是单重态激发态的垂直激发能吗？那这个数 ...

我觉得是可以的。如果用(MR)SF-TDDFT的话，<S^2>的偏离通常会更少，激发能也更准。
作者
Author: SkyHou    时间: 2024-12-28 16:51
hebrewsnabla 发表于 2024-12-28 16:26
我觉得是可以的。如果用(MR)SF-TDDFT的话，的偏离通常会更少，激发能也更准。

好的明白了！谢谢老师！
作者
Author: hebrewsnabla    时间: 2024-12-28 18:09
SkyHou 发表于 2024-12-28 16:10
谢谢老师的建议！老师您的意思是：对于小于2的态，我可以认为它是单重态激发态的垂直激发能吗？那这个数 ...

抱歉，我仔细推了一下，基态为双自由基的三重态激发态S^2也是可以小于2的，并且与对应的单重态激发态差不多。所以并不能用是否小于2来判断。
作者
Author: SkyHou    时间: 2024-12-28 20:26
hebrewsnabla 发表于 2024-12-28 18:09
抱歉，我仔细推了一下，基态为双自由基的三重态激发态S^2也是可以小于2的，并且与对应的单重态激发态差不 ...

那这样的话<S^2>小于2也不太能认为是单重态激发态的吗？



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3