基于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: 596)

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，不知道现在官网上更新了没。

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因子的输出，这个是振动分辨光谱报错，楼里有提到过

好的，谢谢。

wal

keweixu 发表于 2025-4-9 19:42
老师您好，按照您的脚本方法，我重复了一下，sob老师原帖子计算吡啶的重组能和黄昆因子。在用g16计算黄昆因 ...

Low progress不影响huang-rhys因子的输出，这个是振动分辨光谱报错，楼里有提到过

keweixu

老师您好，按照您的脚本方法，我重复了一下，sob老师原帖子计算吡啶的重组能和黄昆因子。在用g16计算黄昆因子提交任务一直报错如下图，想请教一下原因？计算文件在压缩包里了麻烦老师帮忙看下。

Kalinite

tzr 发表于 2024-12-13 18:22
老师，我按照您的输入文件写了之后不能读取initial.chk，请问我是哪里出了错误呀？可以在哪里进行改进嘛？
...

在Windows下用Gaussian，不写文件路径默认在scratch目录，你在计算时initial.chk final.chk是放在了计算目录？

tzr

老师，我按照您的输入文件写了之后不能读取initial.chk，请问我是哪里出了错误呀？可以在哪里进行改进嘛？
报错的最后显示为
Initial State
-------------
Data taken from checkpoint file "initial.chk"
Error termination in NtrErr:
ntran open failure returned to fopen.

附上我的chk和输入文件，请您帮我看一看，谢谢！

任嘉人amber

Kalinite 发表于 2024-10-25 21:34
荧光发射和重组能这两个计算没有关系。
发射计算的是垂直过程，始态和末态没有结构上的区别。

好的，谢谢老师

Kalinite

任嘉人amber 发表于 2024-10-20 14:38
老师，您好，我看您上文提到一般用末态结构分解重组能，但是这个荧光发射例子中oldchk用的是始态的chk，这 ...

荧光发射和重组能这两个计算没有关系。
发射计算的是垂直过程，始态和末态没有结构上的区别。

任嘉人amber

老师，您好，我看您上文提到一般用末态结构分解重组能，但是这个荧光发射例子中oldchk用的是始态的chk，这是什么原因呢，是不是对结构不太敏感，谢谢老师

身在盘丝洞

Kalinite 发表于 2024-9-20 00:01
我觉得可以用四点法计算各个分子的重组能，按照Gaussian或FCClasses分解结果的比例，把四点法得到的重组 ...

老师可以分享一下这篇文献吗，我的也是分解重组能和四点法相差太大