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

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

任嘉人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分解结果的比例，把四点法得到的重组 ...

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

wal

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

好的好的，谢谢老师，我再自己琢磨琢磨
老师您有空的话可以告诉我按照总重组能分配结果的文献是哪篇嘛，如果到时候真这么做了可能还要跟导师或是审稿人解释解释

Kalinite

wal 发表于 2024-9-19 23:24
老师，我是需要分解重组能的，而且需要对比几个分子的huang-rhys因子，若只讨论相对大小可能不太方便横向 ...

我觉得可以用四点法计算各个分子的重组能，按照Gaussian或FCClasses分解结果的比例，把四点法得到的重组能分配到各个振动模式。虽然这么做可能缺乏理论依据，但确实曾见过有文献这么做，不失为一种办法。

如果基团旋转造成的结构变化很大的话，我估计考虑非谐振的改进也很有限，不过值得一试吧。
Gaussian的FCHT计算有很多关于非谐振的选项，但是手册上写的模棱两可，我没法提供更多帮助了。

wal

Kalinite 发表于 2024-9-19 20:39
只要重组能的话，用四点法算就行了，没必要用Gaussian或FCClasses。
如果要分解重组能，此时用Gaussian ...

老师，我是需要分解重组能的，而且需要对比几个分子的huang-rhys因子，若只讨论相对大小可能不太方便横向对比，很苦恼
您觉得做非谐振频率分析对改进结果可能会有帮助吗？

Kalinite

wal 发表于 2024-9-19 10:54
老师，我也遇到了类似的问题，用gaussian和FCclasses内坐标模式得到的重组能比较相似，而且都是大得惊人 ...

只要重组能的话，用四点法算就行了，没必要用Gaussian或FCClasses。
如果要分解重组能，此时用Gaussian或FCClasses得到的结果，从相对大小上还是可以讨论的。

有基团发生旋转，势能面偏离谐振近似较多，大概率就是造成结果不合理的原因，即使换了MOMAP，也未必会得到满意的结果。

wal

Kalinite 发表于 2024-7-16 02:50
如果是这样的话，那可以把两个优化过结构在VMD里align一下，重新做一个freq。如果此时笛卡尔坐标系下能得 ...

老师，我也遇到了类似的问题，用gaussian和FCclasses内坐标模式得到的重组能比较相似，而且都是大得惊人。我尝试在VMD里把两个结构叠起来，发现两个结构基本是重合的（除了在激发过程中发生旋转的那个基团以外）。我还不太会用VMD做align，想问一下老师这种情况的话去align得到正确结果的希望是不是比较渺茫
我向momap提出过申请试用，人家没搭理我，现在这个情况我还有没有别的软件可以试，或者还有什么计算方法可以拯救一下

Kalinite

身在盘丝洞 发表于 2024-9-12 16:20
老师我按照你这篇文章分解了重组能末态为S0 初态为S1是成功的但是结果不太符合，我想看一下S2为初态的情况 ...

大概率是原子顺序不匹配造成的。

身在盘丝洞

老师我按照你这篇文章分解了重组能末态为S0 初态为S1是成功的但是结果不太符合，我想看一下S2为初态的情况在优化时相比于S1加了root=2其余都没变，但是用FCHT却失败了
%oldchk=optfreqs0.chk
%nprocshared=64
%mem=400GB
#P Geom=AllCheck Freq=(ReadFC,FCHT,ReadFCHT)

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

optfreqs2.chk
optfreqs0.chk
下面是报错信息老师是哪里有问题呢
Initial state structure is set in Eckart orientation.
Final state structure is superposed to it.

IAn wrong for default in FilMas.
Error termination via Lnk1e in /public/software/g16/l718.exe at Fri Mar 24 12:20:52 2023.
Job cpu time:       0 days  0 hours  0 minutes 21.5 seconds.
Elapsed time:       0 days  0 hours  0 minutes  1.2 seconds.
File lengths (MBytes):  RWF=     31 Int=      0 D2E=      0 Chk=    324 Scr=      1

Kalinite

Joel 发表于 2024-7-11 12:39
可能是坐标系的问题，我先用您的脚本算了一下，得到的重组能（1.8 eV）明显高于四点法（0.35 eV）算出来 ...

如果是这样的话，那可以把两个优化过结构在VMD里align一下，重新做一个freq。如果此时笛卡尔坐标系下能得到较好结果，那大概率能说明是坐标系的问题了。
不过按理说FCClasses用内坐标就能解决这种问题，97楼说他的例子FCClasses不行，不知道他的例子MOMAP行不行。根据笛卡尔坐标产生内坐标的方式，不同程序之间的差异可能太大了。

Joel

Kalinite 发表于 2024-6-19 20:10
FCClasses EMI/AH modal/FC/INTERNAL COORDS也不能给出合理结果的话，那换用MOMAP大概率也没法解决问题， ...

可能是坐标系的问题，我先用您的脚本算了一下，得到的重组能（1.8 eV）明显高于四点法（0.35 eV）算出来的重组能，然后用momap算了一下，momap给了两个不同坐标下的重组能，基于笛卡尔坐标得到的重组能和用您的脚本也就是用高斯算出来的结果接近，但是基于内坐标算出来的结果就和四点法算出来的结果相近，个人猜测可能是结构优化过程中坐标数值出现了较大的变动？

Kalinite

cokie 发表于 2024-6-17 22:17
老师您好！关于【输出Huang-Rhys因子的方法】有几个问题想向您讨教。

在使用Gaussian 16 A.0.3按您所述 ...

FCClasses EMI/AH modal/FC/INTERNAL COORDS也不能给出合理结果的话，那换用MOMAP大概率也没法解决问题，可能是体系比较特殊，也并不罕见。
个人认为此时讨论相对大小即可，如考察振动模式对总重组能贡献的比例；总重组能仍使用四点法计算的结果，这个是没有问题的。