求助：过渡金属配合物TDDFT计算出现异常，stable=opt后guess=read计算多重态难以判断

SkyHou

各位老师好！我最近在计算一例Pd(II)配合物的激发态时遇到了问题，想请各位老师帮忙看一看问题在哪，万分感谢！
分子结构是这样的：

由于这个分子在实验中表现出了光催化的活性，所以我想计算其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)
   
5. 
   
6. 240110a_0m
   
7. 
   
8. 0 1
   
9. ———基态原子坐标———
   
10. 
    
11. Pd  I 0
    
12. SDD
    
13. ****
    
14. C H N F 0
    
15. 6-311G**
    
16. ****
    
17. 
    
18. Pd I 0
    
19. SDD
    

复制代码

但是输出文件中出现的S1能量非常低，T1甚至为负，让我觉得很困惑（节选了前10个方便各位老师检查）：

1. Excitation energies and oscillator strengths:
   
2. 
   
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.
   
9. 
   
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
    
13. 
    
14. Excited State   3:      Triplet-AU     1.1264 eV 1100.69 nm  f=0.0000  <S**2>=2.000
    
15.      252 -> 254        0.70688
    
16. 
    
17. Excited State   4:      Singlet-AU     1.1283 eV 1098.87 nm  f=0.0012  <S**2>=0.000
    
18.      252 -> 254        0.70702
    
19. 
    
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
    
24. 
    
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
    
29. 
    
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
    
33. 
    
34. Excited State   8:      Singlet-AG     1.6918 eV  732.86 nm  f=0.0000  <S**2>=0.000
    
35.      252 -> 255        0.70612
    
36. 
    
37. Excited State   9:      Singlet-AU     1.6918 eV  732.84 nm  f=0.0133  <S**2>=0.000
    
38.      252 -> 256        0.70609
    
39. 
    
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
   
6. 
   
7. 240110a_0m
   
8. 
   
9. 0 1
   
10. ————基态原子坐标————
    
11. 
    
12. Pd  I 0
    
13. SDD
    
14. ****
    
15. C H N F 0
    
16. 6-311G**
    
17. ****
    
18. 
    
19. Pd I 0
    
20. SDD
    

复制代码

这次计算也是没有报错，正常结束的，但是输出的激发态的多重态很出乎意料，完全没有找到单重态，反而出现了很多近似三重态的激发态：

1. Excitation energies and oscillator strengths:
   
2. 
   
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.
    
14. 
    
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
    
19. 
    
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
    
24. 
    
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
    
41. 
    
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能量很高，所以前五个没有覆盖到吗？想了很久也没想到该怎么解决，谢谢各位老师了！

sobereva

你先看看stable=opt之后波函数变成了什么状态，自旋密度和<S^2>如何。当前肯定已经不是闭壳层单重态了

hebrewsnabla

<S**2> 小于 2的是有自旋污染的单重态。

SkyHou

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原子的自旋密度有点反常，是说明这个所谓稳定的波函数不是闭壳层了吗？所以是不合理的吗？

zjxitcc

无法解决，这是U-TDDFT方法的理论局限性，无法通过“我加一两个fancy的关键词”解决。基态是开壳层单重态时，<S^2>大于0，不是单重态纯态；基于此UHF/UKS波函数做TDDFT计算，每个激发态都有自旋污染，偏离理想值，往往无法分辨激发态所属自旋多重度。同时TD()括号中写singlet/triplet均没用；这也不是某个软件的局限，任何有U-TDDFT方法的软件都有这个问题。这在本论坛上过往“负激发能”的帖子中一般都会顺带提及。解决办法一句话总结：使用更严谨的理论方法，例如MRSF-TDDFT方法（注意，不要使用SF-TDDFT方法），SA-SF-DFT方法，CASSCF-NEVPT2方法等等。

zjxitcc

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>足矣，不需要看其他量。要分清每个量的用途。

SkyHou

hebrewsnabla 发表于 2024-12-28 14:34
小于 2的是有自旋污染的单重态。

谢谢老师的建议！老师您的意思是：对于<S^2>小于2的态，我可以认为它是单重态激发态的垂直激发能吗？那这个数据有可信度吗？

SkyHou

zjxitcc 发表于 2024-12-28 16:06
无法解决，这是U-TDDFT方法的理论局限性，无法通过“我加一两个fancy的关键词”解决。基态是开壳层单重态时 ...

好的好的！我想请教一下老师，在计算前通过分子结构可以判断出它是闭壳层还是开壳层单重态吗？

SkyHou

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信号，说明其基态应当不是多重态。那这样的计算结果是说明物质很可能是开壳层单重态了吗？谢谢老师！

zjxitcc

SkyHou 发表于 2024-12-28 16:14
好的好的！我想请教一下老师，在计算前通过分子结构可以判断出它是闭壳层还是开壳层单重态吗？

除少数结构奇特的体系（比如一个分子里明显少了2个氢原子，烷烃砍掉2个氢原子；历史上被详细研究过的体系，例如O2的最低单重态），一般无法仅凭肉眼预测体系的双自由基特征强弱。

SkyHou

zjxitcc 发表于 2024-12-28 16:19
除少数结构奇特的体系（比如一个分子里明显少了2个氢原子，烷烃砍掉2个氢原子；历史上被详细研究过的体系 ...

我明白了！太谢谢老师了！耽误您的时间解答我的疑惑真是辛苦您了！

hebrewsnabla

SkyHou 发表于 2024-12-28 16:10
谢谢老师的建议！老师您的意思是：对于小于2的态，我可以认为它是单重态激发态的垂直激发能吗？那这个数 ...

我觉得是可以的。如果用(MR)SF-TDDFT的话，<S^2>的偏离通常会更少，激发能也更准。

SkyHou

hebrewsnabla 发表于 2024-12-28 16:26
我觉得是可以的。如果用(MR)SF-TDDFT的话，的偏离通常会更少，激发能也更准。

好的明白了！谢谢老师！

hebrewsnabla

SkyHou 发表于 2024-12-28 16:10
谢谢老师的建议！老师您的意思是：对于小于2的态，我可以认为它是单重态激发态的垂直激发能吗？那这个数 ...

抱歉，我仔细推了一下，基态为双自由基的三重态激发态S^2也是可以小于2的，并且与对应的单重态激发态差不多。所以并不能用是否小于2来判断。

SkyHou

hebrewsnabla 发表于 2024-12-28 18:09
抱歉，我仔细推了一下，基态为双自由基的三重态激发态S^2也是可以小于2的，并且与对应的单重态激发态差不 ...

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