请选择 进入手机版 | 继续访问电脑版

计算化学公社

 找回密码
 现在注册!
查看: 4988|回复: 17

[综合交流] 使用Dushin分解重组能和计算Huang-Rhys因子

[复制链接]

1万

帖子

25

威望

1万

eV
积分
35687

管理员

公社社长

发表于 2016-5-31 07:35:52 | 显示全部楼层 |阅读模式
使用Dushin分解重组能和计算Huang-Rhys因子

文/Sobereva @北京科音   2016-May-31


经常有做电荷转移的人问怎么用Dushin程序、怎么算重组能、怎么把重组能分解为各个振动模式的贡献、怎么计算Huang-Rhys因子,本文就结合实例专门说说这些问题。


1 重组能(reorganization energy)的计算

重组能是基于Marcus理论计算电子转移速率的关键的量,具体分为内重组能和外重组能,前者衡量电子得/失后(或更广义来说,电子态改变后)因几何结构的弛豫导致的体系能量变化,外重组能则对应将周围环境分子重新极化所花费的能量。外重组能不好算,溶液下往往经验性地取0.2~0.6eV,晶体下则会小一个多数量级,也有些文章专门讨论怎么算,如JPCL,1,941(2010)、JACS,130,12377(2008),本文我们只说内部重组能。

看这张示意图,有中性和带电状态两个势能面



1.png

图中重组能有两部分,λ(I)和λ(II)。做Marcus理论计算时说的重组能是指二者绝对值之和。

具体来说,如果charged态是带+1电荷的情况,算出来的重组能用λ+表示,也称空穴重组能λh。如果charged态是带-1电荷,算出来的重组能用λ-表示,也称电子重组能λe。

根据上面的示意图,很容易知道怎么计算λh和λe。例如计算λh,对应的示意图:

2.png


需要以下步骤
(a) 优化中性状态结构
(b) 在(a)的结构下算中性状态的能量E1
(c) 在(a)的结构下算阳离子状态的能量E2
(d) 优化阳离子状态结构
(e) 在(d)的结构下算中性状态的能量E3
(f) 在(d)的结构下算阳离子状态的能量E4
重组能为
λ(I)=|E3-E1|
λ(II)=|E2-E4|
λh= λ(I)+λ(II)

类似地,我们可以算电子态改变过程的重组能。本文我们实际算一个例子,考察吡啶的基态单重态S0和第一三重态激发态T1之间的两种重组能λ(I)和λ(II),示意图如下:

3.png
可见,S0和T1态优化的结构有明显不同,前者原子是在同一个平面上的,是C2v对称性,后者就弯折了,是Cs对称性。

通常做单点计算比几何优化用的级别要高,但为了省事,此例优化和单点计算都用B3LYP/def2SVP在G09 D.01下做。实际上,也只有几何优化和单点都在同一级别时算出的重组能和后文提到的用dushin做重组能分解的结果才有可比性。算出来的图上四个点的能量:
E1=-248.1053119 Hartree
E2=-247.9505667 Hartree
E3=-248.0570726 Hartree
E4=-247.9706883 Hartree
重组能为
λ(I)=E3-E1=30.27 kcal/mol
λ(II)=E2-E4=12.63 kcal/mol
可见,从S0平衡结构垂直跃迁到T1后,结构弛豫对应的重组能为12.63 kcal/mol,这个值比从T1垂直返回S0之后结构弛豫对应的重组能30.27 kcal/mol小很多。这也暗示了,在S0-T1极小点结构之间的方向上,S0极小点处的势能面曲率大于T1极小点处的,因为以相同方式位移,S0态能量变化λ(I)比T1态能量变化λ(II)大得多。


2 重组能的分解

对于上面吡啶的例子,在S0极小点结构下做freq任务,可以得到S0下的3N-6个正则坐标。S0和T1两个极小点之间的位移ΔQ可以分解为这些正则坐标的贡献,第i个模式的贡献记为ΔQ_i。谐振势公式是E=(1/2)*k*x^2,基于这个公式,我们也可以将λ(I)分解为各个正则模式的贡献,即λ_i=(1/2)*k_i*(ΔQ_i)^2,这里k_i代表相应振动模式的力常数。当然了,谐振模型终究只是对实际势能面的近似,所以∑λ_i并不精确等于λ(I),但能定性相符。类似地,在T1极小点结构下做freq得到3N-6个正则坐标后,也可以得到它们对ΔQ和λ(II)的贡献。

S0和T1极小点下分别得到的3N-6个正则坐标是不同的,彼此间是线性变换关系,这个关系也叫Duschinsky旋转或者称Duschinsky混合效应,可表达为Q''=J*Q'+ΔQ,这里Q'和Q''分别代表两个电子态极小点下的正则模式,J称为Duschinsky矩阵,记录了Q'与Q''间线性变换的组合系数。如果两个态极小点下振动模式完全一致,则J是单位矩阵,说明没有混合;J偏离单位矩阵越大说明混合越强烈。

提醒一下,对重组能分解时候都是基于末态极小点结构的正则模式来做分解,所以分解λ(I)要用S0极小点处的正则模式,分解λ(II)要用T1极小点处的正则模式。而不能用比如T1极小点处的正则模式分解λ(I),因为原理上说不通。

做重组能的分解和计算Duschinsky矩阵并不复杂,也有现成的程序可以做,下面介绍的Dushin就是其一。


3 Dushin程序的使用

Dushin程序由Reimers开发的,关于程序的一些原理细节可以看JCP,115,9103(2001),引用Dushin程序也应当引这篇文章。Dushin程序源代码包可以在这里下载:http://bbs.keinsci.com/forum.php?mod=viewthread&tid=319。Dushin程序有自带的说明文档README,但写得比较抽象,这里用人话说一下:

3.1 编译方法

Dushin是Fortran写的,机子里得有Fortran编译器才能编译,Dushin默认的是ifort编译器(gfortran应该也行,我没试)。

在Linux下编译的过程是:下载后解压,把compile里的ifort -g plot-modes.for subs.o recalc-freq.o proj0freq.o bmatred.o ddiag.o dmpower.o dmatinv.o -o ~/bin/plot-modes这一行前头加上#给注释掉,因为压缩包里并没有带plot-modes.for文件,这个本身也用不上。然后执行./compile,就会调用ifort进行编译,编译好的可执行文件会产生在用户主目录的bin目录下面,包括dushin、displace和compare-geom三个文件。然后可以在控制台直接输入dushin看是否能启动,如果不能则需要把这个bin目录添加到$PATH环境变量里。

为了便于那些不会编译或不会Linux的人使用,笔者编译了一份Windows版,下载链接: Dushin_win32_bin.rar (1.1 MB, 下载次数: 89)

评分

参与人数 18eV +88 收起 理由
Puying + 4 好物!
函数与激情 + 5 膜拜~
哇哇吐 + 4
让你变成回忆 + 4 经常读这样的帖子学到了太多东西。感谢!
aqhuangry + 10 谢谢
clvn + 5 好物!
夏之大三角 + 3 赞!
wyf + 4 好物!
小苹果 + 4 谢谢分享
yjmaxpayne + 5 好物!
bomsaude + 5 谢谢
stecue + 5 赞!
liyuanhe211 + 5
卡开发发 + 5 涨姿势
三石草祭 + 5 牛!
978142355 + 5 拜读之
hzfish + 5 赞!
小范范1989 + 5 好物!

查看全部评分

北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)
计算化学公社论坛:http://bbs.keinsci.com(高水平、高人气、综合性计算化学交流论坛)
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加,以免浪费宝贵的空位

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

8

帖子

0

威望

499

eV
积分
507

Level 4 (黑子)

发表于 2016-5-31 07:55:33 | 显示全部楼层
本帖最后由 初初 于 2016-5-31 07:57 编辑

感谢sob老师!

2363

帖子

9

威望

4027

eV
积分
6570

Level 6 (一方通行)

首席卖萌官

发表于 2016-5-31 08:19:56 | 显示全部楼层
板凳!虽然现阶段用不上,但是抢占一个好席位强势围观
She doesn't love me.
Even so,
my heart has been taken away by her.

37

帖子

0

威望

265

eV
积分
302

Level 3 能力者

发表于 2016-5-31 09:56:05 | 显示全部楼层
sob老师太赞了

230

帖子

0

威望

843

eV
积分
1073

Level 4 (黑子)

发表于 2016-5-31 11:03:41 | 显示全部楼层
赞!最近用不到,收藏之先。

83

帖子

0

威望

975

eV
积分
1058

Level 4 (黑子)

发表于 2016-6-2 01:29:06 | 显示全部楼层
Dushin最近我玩了不少,确实是很好用的一个程序。 Sob讲的很仔细到位,赞!

29

帖子

0

威望

241

eV
积分
270

Level 3 能力者

发表于 2016-6-14 06:47:45 | 显示全部楼层
再学习,谢谢群主分享!

61

帖子

0

威望

1331

eV
积分
1392

Level 4 (黑子)

发表于 2016-6-19 17:26:11 | 显示全部楼层
谢谢!

184

帖子

0

威望

1915

eV
积分
2099

Level 5 (御坂)

发表于 2016-6-22 17:57:44 | 显示全部楼层
sob老师,帖子中Huang-Rhys因子计算公式为S_i=λ_i/(h*ν_i),但我看到林圣贤的书籍说的公式,这两个公式怎样转化?
QQ图片20160622174914.png
而且,以前实验室同学说也存在这样的关系:λ_i=(1/2)*v_i*(Q_i)^2;S_i=(1/2)*(Q_i)^2
使用这种简单公式计算的Huang-Rhys因子=0.5*4.112^2=8.45
这样做是否合理?

1万

帖子

25

威望

1万

eV
积分
35687

管理员

公社社长

 楼主| 发表于 2016-6-22 18:53:56 | 显示全部楼层
赵云跳槽 发表于 2016-6-22 17:57
sob老师,帖子中Huang-Rhys因子计算公式为S_i=λ_i/(h*ν_i),但我看到林圣贤的书籍说的公式,这两个公式怎 ...

不要看那个书
看这些,都是高引用的:
DOI: 10.1063/1.1687675
Organic Electronics 11 (2010) 1701–1712
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)
计算化学公社论坛:http://bbs.keinsci.com(高水平、高人气、综合性计算化学交流论坛)
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加,以免浪费宝贵的空位

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

184

帖子

0

威望

1915

eV
积分
2099

Level 5 (御坂)

发表于 2016-6-22 21:24:08 | 显示全部楼层
谢谢sob老师

44

帖子

0

威望

584

eV
积分
628

Level 4 (黑子)

发表于 2016-6-25 18:22:03 | 显示全部楼层
非常感谢社长,momap不算好用的软件

98

帖子

0

威望

867

eV
积分
965

Level 4 (黑子)

发表于 2017-6-9 20:26:51 | 显示全部楼层
sob 老师,我用dushin的输出结果,按照λ=∑ћSiωi计算得到的值 和 30.27和12.63 kcal/mol 相差太多,
能请教一下怎么计算出30.27和12.63 kcal/mol的吗?(程序可以直接输出,但是用λ=∑ћSiωi手动算不出这个结果来)

1万

帖子

25

威望

1万

eV
积分
35687

管理员

公社社长

 楼主| 发表于 2017-6-9 21:56:25 | 显示全部楼层
ZHANGZY 发表于 2017-6-9 20:26
sob 老师,我用dushin的输出结果,按照λ=∑ћSiωi计算得到的值 和 30.27和12.63 kcal/mol 相差太多 ...

30.27和12.63 kcal/mol不是程序输出的,是第一节里直接根据单点能算的
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)
计算化学公社论坛:http://bbs.keinsci.com(高水平、高人气、综合性计算化学交流论坛)
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加,以免浪费宝贵的空位

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

98

帖子

0

威望

867

eV
积分
965

Level 4 (黑子)

发表于 2017-8-4 19:17:34 | 显示全部楼层
sobereva老师,我注意到,您在这个例子中算T1的frequency,不是按照激发态频率分析计算的,而是利用T1态的优化结构,按照T1态的优化结构对应的基态频率计算。这样做是忽略了电子结构的影响,只考虑原子核结构对频率的影响,这样做的依据是什么?这样做的误差是可以接受的吗?
您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

手机版|北京科音自然科学研究中心|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949-1号 )

GMT+8, 2018-11-15 08:59 , Processed in 0.191501 second(s), 28 queries .

快速回复 返回顶部 返回列表