基于Gaussian FCHT计算分解重组能和计算Huang-Rhys因子的Python小脚本

Kalinite

在研究光化学过程时经常需要考察哪些振动模式对重组能起到的明显的贡献，也会需要计算Huang-Rhys因子或Duschinsky矩阵。众所周知，Dushin是常用的分解重组能的辅助程序(http://sobereva.com/330)，也支持输出Huang-Rhys因子或Duschinsky矩阵，但是使用时颇不稳健，经常出现各种各样的报错。然而Gaussian实际上本身就支持给出Duschinsky rotation所定义的shift vector以及Duschinsky矩阵，同时也可以输出Huang-Rhys因子，但关键词深藏在FCHT模块中不为人知。一开始使用Dushin时遇到了很多麻烦，后来发现Gaussian这一功能后解放了双手，最近看到有人问Dushin的报错，故这里分享一个平时用的小脚本，可以结合Gaussian FCHT输出处理数据。
-----------------------------------------------------------------------------------------------------------------------------------------
Gaussian版本：Gaussian16 C.01（其他版本可能需要根据具体输出格式改，没调查过）；需要机子上安装Python。
用法：python gau_reorg.py output.log（假定运行Python的命令是python）
gau_reorg.py (4.75 KB, 下载次数 Times of downloads: 691)

2022-3-31 01:42 上传 Uploaded

点击下载
Click to download

程序会问三个问题：

1. Output decomposition of reorganization energy? Yes, 1; No, 0.
   
2. Output Huang-Rhys factor diagram? Yes, 1; No, 0.
   
3. Output Duschinsky matrix? Yes, 1; No, 0.

复制代码

一般第一个问题都是需要的，选1，会像Dushin一样将重组能分解结果输出在屏幕上；第二个问题控制输出Huang-Rhys因子图的数据（在当前目录下output.log.hr.txt），可以直接导入Origin中作图；第三个问题控制输出Duschinsky矩阵，会将Gaussian中输出Duschinsky矩阵的格式修改为n×n矩阵的格式（在当前目录下output.log.mat.txt），可以直接导入Origin中作图。事实上只要输入一个非0的数都被会认为是需要相关数据。由于Gaussian可以直接输出Huang-Rhys因子，而且在输出格式上比力常数好看，故脚本中实际上是用Huang-Rhys因子计算重组能的分解结果，与根据定义(S_i=\lambda_i/h\nu_i)算没有本质区别。
-----------------------------------------------------------------------------------------------------------------------------------------
如何在Gaussian16 C.01中输出Huang-Rhys因子？（原贴：http://bbs.keinsci.com/thread-28612-1-1.html 5楼）

1. %oldchk=final.chk
   
2. #p geom=allcheck freq(readfc,fcht,readfcht)
   
3. 
   
4. initial=source=chk final=source=chk spectroscopy=onephotonemission
   
5. print=(huangrhys,matrix=JK)
   
6. 
   
7. initial.chk
   
8. final.chk

复制代码

输入文件说明：
1. huangrhys/J/K分别控制输出Huang-Rhys因子、Duschinsky矩阵和shift vector；
2. spectroscopy默认是onephotonabsorption，根据实际来；
3. 没有跃迁偶极的情况，应该写freq(readfc,fc,readfcht)，不过即使写fcht也不会影响输出huangrhys和JK；
4. initial.chk和final.chk分别是初态和末态的chk，需要包含振动信息，用做freq的chk即可；
5. 一般是基于末态结构分解重组能，所以oldchk写final.chk。
多数情况这个任务会报错（比如low progression after class n之类），因为Gaussian会试图计算振动耦合的光谱，只想要这些数据不想要光谱的时候不管就行了。但是要注意：在计算振动耦合光谱时，因为会时常用到RedDim=ClearLowFreq=n来忽略振动波数小于n的低频振动模式的贡献，此时的Huang-Rhys因子与Duschinsky矩阵都是在忽略低频振动模式的基础上输出的，故使用RedDim时不能采用Gaussian给出的数据，另外算一个不加RedDim就完事了。如果看到下面这个报错，很大可能是原子顺序不匹配，需要修改原子顺序并保证两态计算中顺序一致。

1. Initial state structure is set in Eckart orientation.
   
2. Final state structure is superposed to it.
   
3. 
   
4. IAn wrong for default in FilMas.

复制代码

-----------------------------------------------------------------------------------------------------------------------------------------
相关资源：
1. 作图
振动模式：http://sobereva.com/567
Duschinsky矩阵：http://bbs.keinsci.com/thread-4624-1-1.html and http://bbs.keinsci.com/thread-2968-1-1.html
2. Gaussian的FCHT模块
https://gaussian.com/freq/（选择More->ReadFCHT Input）
http://sobereva.com/223
3. 最近的一些应用Huang-Rhys因子分析的文章，太多了，随便弄了三篇：
Lu Wang, Qi Ou, Qian Peng*, and Zhigang Shuai*, J. Phys. Chem. A 2021, 125, 7, 1468–1475
Rongji Zhu; Xi Chen; Na Shu; Yunlong Shang; Yichen Wang; Pu Yang; Yihan Tang; Fei Wang*; Jiawei Xu*, J. Phys. Chem. A 2021, 125, 47, 10144–10154
Teng-Fei He, Ai-Min Ren, Yuan-Nan Chen, Xue-Li Hao, Lu Shen, Bo-Hua Zhang, Tong-Shun Wu, Hong-Xing Zhang, and Lu-Yi Zou*, Inorg. Chem. 2020, 59, 17, 12039–12053

-----------------------------------------------------------------------------------------------------------------------------------------
另外歪个楼分享一下Gaussian16 C.01中计算振动耦合光谱的做法，主要是根据223号博文和与Dr. Lufeng Zou (cnhelp)的交流总结出来的。以荧光发射为例，最基本的输入如下：

1. %oldchk=initial.chk
   
2. #p geom=allcheck freq(readfc,fcht,readfcht)
   
3. 
   
4. initial=source=chk final=source=chk spectroscopy=onephotonemission
   
5. 
   
6. initial.chk
   
7. final.chk

复制代码

输入文件的说明参考前文。输出文件可以在GaussView6中看各个振动耦合峰，选择Results->Vibronic即可，光谱窗口选择左上角Plot->Transition Table可以一目了然地看到各个峰的信息。下面是一些报错和问题（不仅限于此，只遇到过这些，欢迎补充）：
1. checkpoint文件损坏或不存在

1. Data taken from checkpoint file "xxx.chk"
   
2. FileIO operation on non-existent file.

复制代码

2. “Low progression” 很常见也很麻烦

1. ERROR: Low progression after class 2. Total convergence =  0.1%.
   
2. The vibronic spectrum will likely be unreliable. Stopping.

复制代码

Advanced=(ForceFCCalc)
有时包括更多振动激发态时会有一定的progression。Advanced=(ForceFCCalc)关键词忽略该阈值并继续测试更多的振动激发态，并查看是否有明显的progression。
PreScreening=(MaxC1=XX,MaxC2=YY,MaxInt=ZZZ) [Default: XX=13, YY=20, ZZZ=100]
设置初筛的最大值，尤其是MaxInt需要设置一个大得多的值(>=1000)。如果使用Advanced=(ForceFCCalc)后仍然progression很小，那么可能Franck-Condon近似不适用，因为强度很小，即使强制进行FC计算，结果也很可能并不可靠。
RedDim=ClearLowFreq=N
如果不需要完整光谱，可以使用Reduced Dimensionality选项，只研究光谱中感兴趣的部分区域。有些低频振动虽然只是涉及很小的内坐标变化，但涉及的笛卡尔坐标却改变很大，造成难以匹配。例如对于100cm-1以下的频率不感兴趣，使用RedDim=ClearLowFreq=100.0。
3. G09和早期G16版本默认NoSaveNM，前面的freq计算需要写SaveNormalMode (SaveNM)来将振动信息存入.chk文件。官网上电子手册一直没有更新，至少在Gaussian16 C.01版本是默认SaveNM，不知道现在官网上更新了没。

DianaW

Kalinite老师，您好！我用ONIOM(QM/MM)模型优化并算了S0和T1态频率，后面用Gaussian16 A.03中输出Huang-Rhys因子的gjf内容具体如下:
%nprocshared=64
%mem=5000MW
%oldchk=t1-opt.chk
#p geom=allcheck Freq=(ReadFC,ReadFCHT,FC)

initial=source=chk final=source=chk spectroscopy=onephotonabsorption
print=(huangrhys,matrix=JK)

s0-opt.chk
t1-opt.chk

但是出现报错为：
Overlap in IMove.
Error termination via Lnk1e in /opt/soft/g16/l718.exe at Tue Jan 20 22:34:25 2026.
Job cpu time:       0 days  9 hours 54 minutes 47.3 seconds.
Elapsed time:       0 days  0 hours 10 minutes  8.5 seconds.
File lengths (MBytes):  RWF=    847 Int=      0 D2E=      0 Chk=   3350 Scr=      2
想向您请教一下，哪里出了问题？

杲杲出日

老师您好，我想分解S0和anion之间跃迁的重组能，先计算了S0→anion的重组能，输入文件如下：
%oldchk=Y6-anion-freq.chk
#p geom=allcheck freq(readfc,fc,readfcht,intmode) IOp(7/75=-1)

initial=source=chk final=source=chk spectroscopy=onephotonemission
print=(huangrhys,matrix=JK)

Y6-S0-freq.chk
Y6-anion-freq.chk

计算的结果为1072cm-1，该结果和我使用四点法或者MOMAP计算的结果有较大出入（979cm-1）
当我将输入文件的S0和anion互换位置，计算anion→S0的重组能，结果为988cm-1，该结果也和别的方法计算的差异很大（1080cm-1）。存在差异能理解，为什么两个重组能的大小都翻转了呢

杲杲出日

wal 发表于 2025-12-12 21:43
看见FileIO首先尝试把S0和S1的chk对调

神医妙手啊老师

wal

杲杲出日 发表于 2025-12-12 21:02
老师，我是用的是wb97xd泛函，所以将opt和freq分开来做，S0和S1的关键词分别为：
#p freq=saveNM wb97xd/6 ...

看见FileIO首先尝试把S0和S1的chk对调

杲杲出日

老师，我是用的是wb97xd泛函，所以将opt和freq分开来做，S0和S1的关键词分别为：
#p freq=saveNM wb97xd/6-31g(d,p) iop(3/107=0095000000,3/108=0095000000)
#p freq=saveNM wb97xd/6-31g(d,p) iop(3/107=0095000000,3/108=0095000000) td
两个任务顺利完成，没有虚频。当我接下来按照例子计算HR因子时，输入文件为：
%oldchk=Y6-S1-freq.chk
#p geom=allcheck freq(readfc,fcht,readfcht)

initial=source=chk final=source=chk spectroscopy=onephotonemission
print=(huangrhys,matrix=JK)

Y6-S0-freq.chk
Y6-S1-freq.chk

但是算到一般就会报错：Error termination in NtrErr:
NtrErr Called from FileIO.
请问老师这是哪里出错了

清明雨上

老师您好，我按照您的教程计算振动分辨的荧光谱的输入文件为：
%oldchk=s1.chk
#p geom=allcheck freq(readfc,fcht,readfcht)


initial=source=chk final=source=chk spectroscopy=onephotonemission Advanced=(ForceFCCalc) RedDim=ClearLowFreq=60 FORCEPRTSPECTRUM


s1.chk
s0.chk

最后计算得到谱只有一个振动贡献的峰（已上传），请问这正常吗？

陈AG

陈AG 发表于 2025-8-29 08:57
老师，我遇到了这个问题，log输出：ERROR:2 imaginary frequencies found in initial state.终端输出：No.  ...

好像是两个freq都有有虚频导致的，我解决这个问题先

陈AG

老师，我遇到了这个问题，log输出：ERROR:2 imaginary frequencies found in initial state.终端输出：No.        Frequencies/cm^-1        Huang-Rhys Factors        lambda_i/kJ*mol^-1        lambda_i/eV
1        13.5959        Traceback (most recent call last):
  File "/home/server023/Python/gau_reorg.py", line 58, in <module>
    print(huangrhys[i-1],end='\t')
IndexError: list index out of range
，我转成fch文件，查看了也是有振动信息的

LastXuan

wzkchem5 发表于 2025-6-30 21:43
如果要算S1频率，必须写td freq。这和计算目的是不是算黄昆因子无关

感谢老师回复，我按这个算一下

wzkchem5

LastXuan 发表于 2025-6-30 20:42
请问老师一个问题：如果我要计算S0→S1的黄昆因子，其中S1态在结构优化完之后的S1态频率计算关键词是写freq ...

如果要算S1频率，必须写td freq。这和计算目的是不是算黄昆因子无关

LastXuan

请问老师一个问题：如果我要计算S0→S1的黄昆因子，其中S1态在结构优化完之后的S1态频率计算关键词是写freq还是td freq呢

Kalinite

keweixu 发表于 2025-4-9 21:00
老师，您好，我想再请教一下，我用g16+python脚本得到吡啶T1-S0的重组能和黄昆因子的频率贡献（386 cm-1） ...

首先需要确认，你在使用g16 g09时是否保证了各种设置完全一致（如g16有没有写g09default等）。
不存在某一个态的重组能，只有某一个过程的重组能，如果你要计算S0->T1 T1->S0，那么做两次FC计算即可。

keweixu

老师，您好，我想再请教一下，我用g16+python脚本得到吡啶T1-S0的重组能和黄昆因子的频率贡献（386 cm-1），与之前按照sob老师基于g09+dushin分解得到S0的重组能，频率贡献（426 cm-1）差别很大，那个数据准一些？还是我计算的问题？还有个疑问就是dushin可以分解得到S0和T1的重组能，那么通过python脚本可以得到T1的重组能嘛？

keweixu

wal 发表于 2025-4-9 19:54
Low progress不影响huang-rhys因子的输出，这个是振动分辨光谱报错，楼里有提到过

好的，谢谢。