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

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

TANGYIHAN

多谢Kalinite老师！

qidong8

Kalinit老师你好，我参照您的的输入文件计算重组能，但总出现内存错误

1. Error: segmentation violation
   

复制代码

末态为基态，优化且做了频率计算

1. #p opt freq pbe1pbe/6-31g** scrf=(pcm, solvent=toluene)

复制代码

初态为S1激发态，同样做了优化和频率计算

1. #p opt td(nstates=6) freq=saveNM pbe1pbe/6-31g** scrf=(pcm,solvent=toluene)

复制代码

gaussian版本为g16 C01，想向您请教一下，哪里出了问题？

Kalinite

qidong8 发表于 2022-4-20 21:28
Kalinit老师你好，我参照您的的输入文件计算重组能，但总出现内存错误

末态为基态，优化且做了频率计算
...

Error: segmentation violation 不能看出是什么错，可以提供一下输出文件末尾的报错信息看看是什么情况。

qidong8

Kalinite 发表于 2022-4-20 21:57
Error: segmentation violation 不能看出是什么错，可以提供一下输出文件末尾的报错信息看看是什么情况。

log末尾是

1. Using the following non-standard input for vibronic spectroscopy:
   
2. initial=source=chk final=source=chk spectroscopy=onephotonemission
   
3. print=matrix=(huangrhys,JK)
   
4. Unrecognized value for Spectra: (
   
5. Error termination via Lnk1e in /home/liujunyuan/Software/g16/l718.exe at Wed Apr 20 21:19:57 2022.
   
6. Job cpu time:       0 days  0 hours  0 minutes  1.0 seconds.
   
7. Elapsed time:       0 days  0 hours  0 minutes  0.9 seconds.
   
8. File lengths (MBytes):  RWF=     15 Int=      0 D2E=      0 Chk=     89 Scr=      1

复制代码

Kalinite

qidong8 发表于 2022-4-21 16:19
log末尾是

没见过这个问题，可能和chk有关，可以发邮件给cnhelp问问

Kalinite

qidong8 发表于 2022-4-21 16:19
log末尾是

应该是print=(huangrhys,matrix=JK)，之前打错了，抱歉

snljty

老师您好，我这边用您的方法做测试的时候，还是提示不存在文件。测试用的卢老师那篇博文里吡啶的例子。我把文件传上来了，可以麻烦您抽空看一下么？十分感谢！ test_dushin.zip (1.04 MB, 下载次数 Times of downloads: 19)

2022-4-26 14:28 上传 Uploaded

点击下载
Click to download
chk对应的Gaussian版本可能和您的不一样，那样的话麻烦您重跑一个S0.gjf和T1.gjf了，应该是几秒钟以内。谢谢！

Kalinite

snljty 发表于 2022-4-26 14:28
老师您好，我这边用您的方法做测试的时候，还是提示不存在文件。测试用的卢老师那篇博文里吡啶的例子。我把 ...

你算的是S0-T1之间的，这是自旋禁阻的过程，Gaussian不能给出跃迁偶极，而FCHT中的“HT”需要跃迁偶极导数，因此在读取.chk中相关内容时报错“File operation on non-existent file”。虽然计算Huang-Rhys因子不需要这部分信息，但是Gaussian的这个模块是做光谱的，仍然会试图读取，并不影响Huang-Rhys因子和Duschinsky矩阵的输出。
改成“Freq=(ReadFC,ReadFCHT,FC)”可以正常输出你需要的信息，报的是光谱计算的错误。
FCHT.log (60.84 KB, 下载次数 Times of downloads: 18)

2022-4-26 18:23 上传 Uploaded

点击下载
Click to download

snljty

Kalinite 发表于 2022-4-26 18:24
你算的是S0-T1之间的，这是自旋禁阻的过程，Gaussian不能给出跃迁偶极，而FCHT中的“HT”需要跃迁偶极导 ...

收到，谢谢老师！

xyran

Kalinite老师，您好！
    我最近在做重组能分解。想用您给的python方法进行分解。但是在分解过程中不知道如何运用python读取我已优化好的文件。老师能告知一读取的相应代码吗？谢谢！

Kalinite

xyran 发表于 2022-5-2 14:03
Kalinite老师，您好！
    我最近在做重组能分解。想用您给的python方法进行分解。但是在分解过程中不知道 ...

假定你的电脑上运行python的命令是python，输出文件是output.log，那么运行python gau_reorg.py output.log即可。

watfros

很棒的内容，之前接触过一点类似的，现在感觉gaussian还是很多隐藏功能的

anlancx

老师，对于黄昆因子，输出的结果是这样：
==================================================
                     Huang-Rhys Factors
     ==================================================

     Mode num.      1 - Factor:   0.684549D+01
     Mode num.      2 - Factor:   0.103728D-12
     Mode num.      3 - Factor:   0.168417D-15
     Mode num.      4 - Factor:   0.234325D+00
     Mode num.      5 - Factor:   0.137321D+01
     Mode num.      6 - Factor:   0.610608D-01
     Mode num.      7 - Factor:   0.155324D-13
     Mode num.      8 - Factor:   0.135165D+00
     Mode num.      9 - Factor:   0.163153D+00
     Mode num.     10 - Factor:   0.480870D-15
     Mode num.     11 - Factor:   0.235686D-02
     Mode num.     12 - Factor:   0.180464D+00
     Mode num.     13 - Factor:   0.316307D+01
     Mode num.     14 - Factor:   0.474745D-14
     Mode num.     15 - Factor:   0.993685D-15
     Mode num.     16 - Factor:   0.596251D+00
     Mode num.     17 - Factor:   0.203202D-15
     Mode num.     18 - Factor:   0.935426D-18
     Mode num.     19 - Factor:   0.560257D-15
     Mode num.     20 - Factor:   0.165495D+00
     Mode num.     21 - Factor:   0.937353D-16
     Mode num.     22 - Factor:   0.275564D+00
     Mode num.     23 - Factor:   0.177246D-03
     Mode num.     24 - Factor:   0.325253D-15
     Mode num.     25 - Factor:   0.226731D-13
     Mode num.     26 - Factor:   0.507573D-01
     Mode num.     27 - Factor:   0.346544D-02

对于这么多factor，我们该看哪个呢？

谢谢老师。

Kalinite

anlancx 发表于 2022-5-6 21:20
老师，对于黄昆因子，输出的结果是这样：
==================================================
         ...

Huang-Rhys因子定义为：S_i = \lambda_i / (h \nu_i)。
一般看振动模式所贡献重组能的大小，即\lambda_i的大小。贡献较大的振动模式往往是与始末态结构变化有一定对应关系的，当然也可能会有多个振动模式都具有比较明显的贡献。虽然由谐振子近似给出的重组能分解结果可能并不合理，因为谐振子模型对势能面近似地不一定好，但相对大小一般还是具有参考意义的。可以参考一下前面几篇文献里的做法。