求助：双重态自由基td吸收波长计算和实验波长对比不上

mly03

体系是由C/H/N组成的90个原子左右的碳中心自由基，目前采用UB3LYP/6-31G*进行了基态结构优化，在此基础上进行了吸收谱计算，分析D0-D1的跃迁的轨道贡献，判断为电荷转移激发为主，然后采用wb97xd/def2-svp，cam-b3lyp/def2-svp重新计算吸收波长（给了隐式溶剂模型），得到的是D1态吸收波长分别为742nm和736nm，同实验570nm相差很大，并且文件中显示D1态对应<s**2>均大于0.82。
想请教一下是否是泛函和基组设置不当问题？或者是否应该换一种方法计算？目前看到有推荐BDF的x-tddft对于开壳层体系处理比较好。

wzkchem5

首先纠正一下：自由基的基态叫做D0不叫S0，第一二重态激发态叫做D1而非S1，以此类推。
这个自旋污染还是比较明显的，不妨用BDF算一下。即使没买BDF，也可以申请14天免费试用（并附送计算所需超算机时），对结果满意再购买授权即可。详见https://mp.weixin.qq.com/s/er7-XOt0MfLrN_-5itZEwQ
此外实验上看到的是第一个亮态，未必是D1态。如果D2的振子强度远大于D1，可能计算的D2才是和实验那个峰对应的。

zjxitcc

对于二重态基态的体系，它的基态是D0，最低几个激发态是D1, D2, Q1等等，不存在S0, S1，符号不能混用。你这种描述下，建议使用BDF的X-TDDFT。

我不知道你的具体体系，这里用二重态的氰基CN·作为示例，PySCF输入脚本为

1. from pyscf import gto, dft, lib
   
2. from mokit.lib import py2bdf
   
3. 
   
4. lib.num_threads(12) # 12 CPU cores
   
5. mol = gto.M()
   
6. mol.atom = '''
   
7. C   0.0   0.0   0.0
   
8. N   0.0   0.0   1.14660
   
9. '''
   
10. 
    
11. mol.basis = 'def2-TZVP'
    
12. mol.charge = 0
    
13. mol.spin = 1
    
14. mol.verbose = 4
    
15. mol.build(parse_arg=False)
    
16. 
    
17. mf = dft.ROKS(mol)
    
18. mf.xc = 'bhandhlyp'
    
19. mf.grids.atom_grid = (99,590)
    
20. mf.max_cycle = 128
    
21. mf.max_memory = 24000 #MB
    
22. mf.kernel()
    
23. 
    
24. if mf.converged is False:
    
25.     mf = mf.newton()
    
26.     mf.kernel()
    
27. 
    
28. mo, _, stable, _ = mf.stability(verbose=5, tol=1e-7, return_status=True)
    
29. if not stable:
    
30.     dm = mf.make_rdm1(mo, mf.mo_occ)
    
31.     mf.kernel(dm0=dm)
    
32.     mo, _, stable, _ = mf.stability(verbose=5, tol=1e-7, return_status=True)
    
33.     if not stable:
    
34.         raise OSError('PySCF ROHF/ROKS wfn is still not stable.')
    
35. 
    
36. py2bdf(mf, 'CN_bdf.inp', xtddft=True)
    

复制代码

提交PySCF任务，即运行

1. python CN.py >CN.out 2>&1

复制代码

正常结束后会产生适用于BDF的X-TDDFT输入文件、基组文件和轨道文件CN_bdf.inp, CN_BDF.BAS, CN_bdf.scforb。可以直接提交BDF任务，或者把计算文件发到有安装BDF的集群上进行X-TDDFT计算。例如运行

1. export BDFHOME=$HOME/software/bdf-pkg
   
2. export BDF_TMPDIR=$HOME/tmp/bdf
   
3. export OMP_NUM_THREADS=12
   
4. export OMP_STACKSIZE=8000M
   
5. $BDFHOME/sbin/bdfdrv.py -tmpdir $BDF_TMPDIR -keeptmpdir -r CN_bdf.inp >CN_bdf.out
   

复制代码

这个体系很小，几十秒结束。从BDF输出文件可以看到二重态基态和激发态信息

1.   No. Pair   ExSym   ExEnergies     Wavelengths      f     D<S^2>          Dominant Excitations             IPA   Ova     En-E1
   
2. 
   
3.     1   A    2   A    0.8754 eV       1416.38 nm   0.0029   0.0371  96.3% CO(bb):   A(   5 )->   A(   7 )   7.014 0.550    0.0000
   
4.     2   A    3   A    0.8754 eV       1416.38 nm   0.0029   0.0371  96.3% CO(bb):   A(   6 )->   A(   7 )   7.014 0.550    0.0000
   
5.     3   A    4   A    3.2825 eV        377.71 nm   0.0277   0.0916  86.7% CO(bb):   A(   4 )->   A(   7 )   8.506 0.706    2.4071
   
6.     4   A    5   A    7.5607 eV        163.99 nm   0.0020   1.9790  23.3% CV(aa):   A(   5 )->   A(   8 )  14.188 0.819    6.6853
   

复制代码

如果把上述PySCF输入脚本最后一行xtddft=True修改为socxtda=True，则产生的BDF输入文件可以同时计算二重态基态，二重态激发态，四重态激发态，以及D-Q旋轨耦合。

sobereva

“均大于0.82”描述不够确切。如果比如是0.821，那自旋污染还是在可以接受的程度

单纯想排除自旋污染带来的不良影响，可以用免费的ORCA的ROCIS/DFT，比较易用，并且北京科音高级量子化学培训班（http://www.keinsci.com/KAQC）非常详细讲ORCA的一节还专门给了具体讲解和例子。只不过这个没有解析梯度，因而对大体系优化激发态费劲，但如果如前所述当前自旋污染不大，那么UTDDFT做激发态优化还是可以用的。

当前用的基组可以改大点，诸如6-311G*（至少是电子激发关系密切的区域）

计算级别只是结果不好的可能性之一，注意构象、聚集状态、溶剂效应（不是光设个隐式溶剂模型就完事）等其它因素

mly03

wzkchem5 发表于 2026-3-20 17:02
首先纠正一下：自由基的基态叫做D0不叫S0，第一二重态激发态叫做D1而非S1，以此类推。
这个自旋污染还是比 ...

谢谢老师，受教，想再请教一下老师，我是想计算该分子的理论最大吸收波长，并同实验给出的值相比较，然后我认为做吸收谱计算得到的D1态对应的波长理论上应该是最大的吸收波长，所以我一直是拿D1态的波长数据和实验比对的，请问老师这个思路是正确的吗？如果不正确的话，要与实验作比较应该重点考虑哪些方面呢？

mly03

zjxitcc 发表于 2026-3-20 17:16
对于二重态基态的体系，它的基态是D0，最低几个激发态是D1, D2, Q1等等，不存在S0, S1，符号不能混用。你这 ...

谢谢老师，我后续用BDF尝试一下

mly03

sobereva 发表于 2026-3-20 19:31
“均大于0.82”描述不够确切。如果比如是0.821，那自旋污染还是在可以接受的程度

单纯想排除自旋污染带 ...

感谢sob老师，我再分析一下溶剂和构象的问题

wzkchem5

mly03 发表于 2026-3-21 14:26
谢谢老师，受教，想再请教一下老师，我是想计算该分子的理论最大吸收波长，并同实验给出的值相比较，然后 ...

不正确。我已经说了，如果D2（或D3, D4, ...）的振子强度明显高于D1，可能实验上看到的那个峰是D2而非D1。应当把更高的几个激发态也都算出来，做高斯展宽，再和实验谱图比较