【龙讯module小课堂】“光”怪陆离：PWmat计算光学性质（三）

最新推荐文章于 2024-03-22 10:35:54 发布

龙讯旷腾

最新推荐文章于 2024-03-22 10:35:54 发布

阅读量444

点赞数

分类专栏：材料计算module小课堂文章标签：量子力学物理学

本文链接：https://blog.csdn.net/PWmat/article/details/124456725

版权

材料计算module小课堂专栏收录该内容

5 篇文章 2 订阅

订阅专栏

往期 PWmat计算光学性质前两期回顾

点击查看第一期

点击查看第二期

前言

在上一期中（链接🔗），我们已经介绍了TDDFT方法和它在线性响应下的形式LR-TDDFT。这种方法直接求解了体系的响应函数在频域的信息，可以较好的描述小分子或团簇在弱场激发下的吸收峰的位置和极化强度等特性。TDDFT还可以直接求解时域内电子-原子核的含时演化，即real time TDDFT(rt-TDDFT)。这种做法不仅可以用来求解光吸收，还能处理许多强场下电子的动力学问题和离子的运动过程。本文将从rt-TDDFT的基本形式开始，介绍PWmat中的rt-TDDFT，再着重介绍rt-TDDFT在光学计算中的应用。在最后给出PWmat中利用rt-TDDFT计算光学性质的module。

PWmat中的rt-TDDFT

相比于LR-TDDFT，rt-TDDFT更关注TDDFT在动力学的那一面。因此，我们也可以将rt-TDDFT归为一种广义的非绝热分子动力学（nonadiabatic molecular dynamic, NAMD）。目前最常用的NAMD是量子-经典相结合的NAMD，即电子的波函数遵循量子力学的含时演化，而原子核依然遵循经典力学。与Born-Oppenheimer (BO)近似下的分子动力学 (MD)不同，在NAMD中，原子核不再只属于一个确定的势能面，它的势能面可以由一个平均势能面代替（Ehrenfest dynamics）1或者在不同电子态的势能面之间跳跃（surface hopping）2-4。

相比于NAMD，rt-TDDFT具有自己身为TDDFT的传统艺能，即可以描述电子-电子相互作用以及体系对外场的响应。在rt-TDDFT的框架下，我们不仅可以自己施加外场并模拟体系的电子结构对外场的响应（光激发），还能模拟被激发的载流子的动力学过程甚至是载流子的动力学影响下的原子核的动力学（这一点在NAMD中是被忽略的）。不过，作为代价，rt-TDDFT可以模拟的体系比NAMD小（一般少于100原子）。值得一提的是，虽然NAMD中原子核的运动轨迹都是用基态的分子动力学模拟的，但是这种做法在大体系中具备一定的合理性。其一，当电子数目很多时，少量电子的激发并不能显著影响体系的性质；其二，载流子在能带上的迁移会导致单一能级上的载流子寿命很短，原子核还来不及响应。

在TDDFT的理论框架下，系统的哈密顿量仍然可以由体系的电荷密度确定

电荷密度为

在传统的Ehrenfest dynamics下有

当我们求得了t 时刻的波函数之后，该时刻原子核的受力就确定了，即t+dt 时刻的原子核的位置也确定了，这一点与BOMD类似。而t+dt 时刻的波函数则可以有如下的演化方式

为了能进行数值计算，其中的指数项需要用Taylor展开，而保持展开精度的前提是

在平面波的基组下dt<10-3 fs，也就是说，通常的rt-TDDFT只能模拟阿秒（as）量级的过程，如果想与实验对比（时间尺度在fs或者ps量级），则需要支付千倍以上的计算时间，这严重制约了rt-TDDFT的应用。针对这个问题，PWmat中的rt-TDDFT做了如下的改进。

首先，用t 时刻的绝热本征态来展开含时的波函数

其中绝热的本征态是由该时刻的Kohn-Sham(KS)方程决定的，即

带入含时的哈密顿量

两端展开则有

在等式两边

由于KS的本征态是正交归一的

于是

通过调整作为基组的绝热波函数的数量，我们可以有效地减少求解 $C_{ji}\left ( t \right )$ 的矩阵的维度。由于能量的本征值只出现在对角项，而对角项中不含dt。因此在这种形式下，dt 的大小不再依赖于哈密顿量。但是，这个等式中包含了对时间的导数项，为了精确求解对时间的导数项，步长dt 依然不能取大。首先是等式左边

在t 时刻的展开系数与t-dt 的展开系数有关。非对角项中的 $V_{il}\left ( t \right )$ 也与dt 有关

这一项让t 时刻与t-dt 的非绝热波函数联系起来，要求解t 时刻的绝热展开系数，我们需要上一个时刻的展开系数和上一个时刻用于展开的所有绝热波函数。此外，还有一个根本性的问题没有解决，那就是H(t)怎么求，如果还是要依靠公式(1)，则dt 的量级又会被H(t)绑架，到头来我们只是缩小了问题的维数。

为了能使用更长的步长，我们采用哈密顿量随时间演化的线性近似5, 6。考虑体系从 $t_{1}$ 到 $t_{2}$ 的含时演化过程，假设哈密顿量在 $t_{1}$ 时刻由该时刻的基组 $\left \{\phi _{i} \left ( t_{1} \right )\right \}$ 展开后的矩阵是对角化的，且哈密顿量在 $t_{2}$ 时刻以相同的基组 $\left \{ \phi _{i}\left ( t_{2} \right ) \right \}$ 展开一样是对角化的，则称哈密顿量在 $t_{1}$ 到 $t_{2}$ 的时间段的演化是线性的。此时，对于 $t_{1}$ 到 $t_{2}$ 之间的任意时刻t 有

它相当于一种平均化处理，把两个时刻之间复杂的演化过程用一个平均的线性过程来代替，既然是要平均化，那么这两个时刻之间的差Δt 就可以适度的拉长（接近0.1 fs量级）。需要注意的是，为了保证(4)中求导的准确性，我们仍然需要在 $t_{1}$ 到 $t_{2}$ 之间不断使用小dt（~0.001 fs）来实现波函数随时间的演化。但是与之前不同的是，我们只需要知道 $\left \{\phi _{i} \left ( t_{1} \right )\right \}$ $\left \{ \phi _{i}\left ( t_{2} \right ) \right \}$ 以及对应的本征值，就可以直接靠对角化（5）中的H(t)来得到任意时刻的 $\left \{\phi _{i}\left ( t \right ) \right \}$ ，再求解(4)代入(2)即可。不必再对每一个H(t)都做一次(3)的对角化。这就相当于节省了N个DFT计算（N的数量级~100）。如果真的能实现上述算法，则计算速度会暴增。

现在的问题是，我们如何在 $t_{1}$ 的时候就得到H( $t_{2}$ )呢?由于时间间隔较长，显然不能从 $t_{1}$ 的波函数利用(1)得到 $t_{2}$ 时刻的波函数了。为此，我们使用 $t_{1}$ 时刻的电荷密度n( $t_{1}$ )通过数值方法线性外推得到n( $t_{2}$ )，然后利用n( $t_{2}$ )得到H( $t_{2}$ )，用这个粗糙的H( $t_{2}$ )在 $t_{1}$ 到 $t_{2}$ 之间以小dt 求解(4)，这样逐步演化后，在t= $t_{2}$ 时会得到一个新的 $n^{''}\left ( t_{2} \right )$ 。我们混合n( $t_{2}$ )和 $n^{''}\left ( t_{2} \right )$ ，得到 $n^{'}\left ( t_{2} \right )$ ，再重复之前的流程，这样做自洽直到两次循环的 $n^{''}\left ( t_{2} \right )$ 之差足够小。虽然进行了自洽循环，但是这种计算方法还是能远快于对每个小dt 都求解一次KS方程(3)。这种看起来简单粗暴的近似却异常的有效，它可以在不牺牲准确性的前提下，把步长的数量级提高到0.1 fs（通常在0.2 fs内比较准确），有时甚至可以达到0.5 fs6。

以上的方法保证了PWmat具备非常高效的rt-TDDFT计算能力。在此基础之上，我们还提供了自旋轨道耦合下的含时哈密顿量7, 8；我们引入了退相干时间9，Boltzmann因子10, 11和自然轨道分支（natural orbital branching, NOB）10, 11方法以弥补Ehrenfest dynamics的缺陷。基于以上高效、完善且强力的rt-TDDFT，PWmat在超快动力学模拟的领域拥有了丰富的计算成果，如图1所示，包括模拟：团簇的等离激元12，光致超快相变13，光致退磁过程7, 8，高能粒子撞击14和辐照分子损伤过程11。

图1 PWmat在超快动力学领域的代表性应用成果

提到PWmat，人们总是联想到快得离谱的杂化泛函计算以及获得戈登贝尔奖的分子动力学。在此，笔者拼着跑题的危险，隆重的安利一波PWmat的rt-TDDFT。至此，笔者终于想起来了这一期的重点其实是rt-TDDFT计算光学性质，于是诸如Boltzmann因子, 自旋轨道耦合和NOB等操作将会放到以后相应的专题中详细介绍（坑好像越挖越多）。让我们开始进入rt-TDDFT与光学性质的纠葛。

用rt-TDDFT计算光吸收

上一部分在跑题和没跑题的奇怪叠加态中讲述了PWmat是怎么高效地描述波函数和哈密顿量的实时演化的。在光学性质的第一期中，我们已经说明了计算光学性质的核心在于计算介电函数；在上一期中，我们说明了计算介电函数的过程可以简化成求解电极化率的虚部。在rt-TDDFT的框架下，我们可以自己给体系施加一个随时间变化的外场E(t)，然后计算体系的极化随时间的演化P(t)，然后做一个Fourier变换

简单粗暴又直接。在计算分子和团簇等孤立体系时，这种做法是非常有效的,这是因为无论是给孤立体系施加的电场还是孤立体系在电场下的极化都能有正确的定义。在PWmat之类的平面波计算软件中，我们采用了带有真空层的大晶胞来模拟孤立体系，假设电场为 $t_{0}$ 时刻沿着x方向的脉冲，晶格沿着x方向的起始位置为 $x_{0}$ ，沿着x方向的长度为 $L_{x}$ ，则由于施加外电场带来的外加势能为

如果要保证跨过晶格的势能是连续的，则电场强度需要满足

显然，电场在晶格的边界处会突然发散一次。由于晶格的边界处是真空层，当真空层取的足够大时，这个突变对计算结果没有影响。

在将上述方法扩展到晶体中时，马上就出问题了。首先，晶体在晶格的边界是连续的，如果交接处发生势能跳变，那么不用算都知道肯定是错误的。你或许可以尝试加一个强度随r 变化并保证在晶格边界处强度为0的电场，但随后你还需要给这种势场找一个实际意义。针对这个问题，解决思路是引入矢势A(t)，将位移空间的规范转移到动量空间。刚才那句的人话版本：由于晶格的周期性，外场的影响不能再加到势能项上了，此时需要在动能项上做文章，即加一个矢势将动量变为正则动量，在原子单位制下有

类似的变换在电磁学中很常见，此时的电场仍然是一个δ函数，在数值上，我们取的是高斯波包

解决了电场的问题之后，我们再来看极化强度。被铁电材料虐过的读者应该瞬间明白了我要说什么，在周期性边界条件下，根本不能定义极化强度P，它会随着晶格的取法而改变。参考电场的解决方案，我们是否也能找到一个物理量来代替P 呢？解决之道隐藏在光学性质的第一期之中

由介质的连续方程，有

于是

只要能解除j 对r 的依赖，我们就能用电流密度来代替极化强度。这很容易做到，只要用Bloch波函数展开电流密度即可

最终

需要注意的是，电介质中是不会有净电流密度的，因此j 在整个演化过程中的期望应该是0，在最后求积分和傅里叶变换的过程中应该注意这一点。

PWmat中利用rt-TDDFT计算光吸收的module

PWmat针对孤立体系和周期性体系提供了两个module以使用rt-TDDFT计算光吸收：module 2（孤立体系）和module 18 (周期性体系)。

两个module的差别主要表现在外场的施加方式。正如上一节中讲到的，孤立体系可以直接施加实空间的外场

求得相应的极化强度 $P_{\mu }\left ( t \right )$ ，然后得到

在module 2中，只需要进行一个提供外场的高精度的job=TDDFT即可。为了精确的描述 $P_{\mu }\left ( t \right )$ ，步长控制在dt =0.01 fs较为合适，将MD的初始和终末温度都设为0 K，这样原子核就不会移动。外场的形式则由IN.TDDFTOPT文件提供，PWmat中为孤立体系提供了多种外场形式，具体内容可以阅读Manual上的B.16。

如图2所示，以 $CH_{4}$ 分子为例，计算得到的(a) $P_{x}\left ( t \right )$ 和(b)介电函数虚部的xx分量。

图2 利用rt-TDDFT计算CH4单分子的光吸收

在周期性体系中，则需要提供矢势A(t)和电流密度 j(t)，注意，这两个物理量都属于倒空间，因此IN.TDDFTOPT中需要OUT.MDDIPOLE.KSPACE=T。如图3所示，以Si的原胞为例，(a)为施加的矢势，由于它源自一个类δ函数（窄高斯波包）对时间的积分，因此其形式接近一个阶梯函数。为了精确的计算跳变点附近的响应性质，在前200步内，采用0.001 fs的步长；对此后矢势较为平坦的区域则使用0.01 fs的步长再续算1000步。图3(b)为距离分离的杂化（range separated hybrid, RSH）泛函计算得到的响应电流。最终求得的介电函数的虚部则在图3(c)中，我们将PBE的结果也放在其中作为对比。

图3 利用rt-TDDFT计算Si原胞的光吸收

1. Li, X.; Tully, J. C.; Schlegel, H. B.; Frisch, M. J. The Journal of Chemical Physics 2005, 123, (8), 084106.

2. Tully, J. C. The Journal of Chemical Physics 2012, 137, (22), 22A301.

3. Wang, L.; Akimov, A.; Prezhdo, O. V. The Journal of Physical Chemistry Letters 2016, 7, (11), 2100-2112.

4. Larsen, R. E.; Bedard-Hearn, M. J.; Schwartz, B. J. The Journal of Physical Chemistry B 2006, 110, (40), 20055-20066.

5. Wang, Z.; Li, S.-S.; Wang, L.-W. Physical Review Letters 2015, 114, (6), 063004.

6. Ren, J.; Vukmirović, N.; Wang, L.-W. Physical Review B 2013, 87, (20), 205117.

7. Chen, Z.; Wang, L.-W. Science Advances 5, (6), eaau8000.

8. Chen, Z.; Luo, J.-W.; Wang, L.-W. Proceedings of the National Academy of Sciences 2019, 116, (39), 19258-19263.

9. Bedard-Hearn, M. J.; Larsen, R. E.; Schwartz, B. J. The Journal of Chemical Physics 2005, 123, (23), 234106.

10. Wang, L.-W. The Journal of Physical Chemistry A 2020, 124, (43), 9075-9087.

11. Cai, Z.; Chen, S.; Wang, L.-W. Chemical Science 2019, 10, (46), 10706-10715.

12. Liu, W.-H.; Luo, J.-W.; Li, S.-S.; Wang, L.-W. Physical Review B 2020, 102, (18), 184308.

13. Ma, J.; Wang, Z.; Wang, L.-W. Nature Communications 2015, 6, (1), 10107.

14. Bi, G.; Kang, J.; Wang, L.-W. Physical Chemistry Chemical Physics 2017, 19, (13), 9053-9058.

【龙讯module小课堂往期回顾】

【龙讯module小课堂】浅谈对gap的认识：PWmat中修正gap的module

【龙讯module小课堂】“光”怪陆离：PWmat计算光学性质（一）

【龙讯module小课堂】“光”怪陆离：PWmat计算光学性质（二）

—E N D—

北京龙讯旷腾科技有限公司是成立于2015年的国家高新技术企业，是国内材料计算模拟工具软件研发创新的领导者，致力于开发满足“工业4.0”所需的原子精度材料研发Q-CAD（quantum-computer aided design）软件。公司自主开发的量子材料计算软件PWmat（平面波赝势方法并基于GPU加速）可以进行电子结构计算和从头算分子动力学模拟，适用于晶体、缺陷体系、半导体体系、金属体系、纳米体系、量子点、团簇和分子体系等。

官网为www.lonxun.com

BBS论坛：www.bbs.pwmat.com