【文献复刻】计算vdW异质结中热载流子冷却
Phys. Chem. Chem. Phys
01
NAMD 背景介绍
01
NAMD 背景介绍
在各类光物理与光化学过程当中，均会牵涉到激发态载流子动力学过程，诸如电荷弛豫、复合以及输运等等。光激发或者电子注入将初始的平衡状态打破，所产生的热载流子在其演化进程中，会与原子核产生强烈耦合。此时，将电子自由度与原子核自由度分开考虑的绝热近似不再适用，需要引入非绝热过程。
根据计算精度和体系大小，人们可以从不同的理论框架下发展非绝热方法。对小分子体系，人们可以使用将电子和原子核都放在量子力学框架下的全量子方法。对于大小适中的体系，人们往往让电子的波函数继续遵循量子力学的演化规律，而让原子核遵循经典力学的量子-经典混合方法（mixed quantum-classic, MQC）。在这个框架下原子核势能面可由平均势能面予以替代或者在不同电子态势能面之间进行跳跃。对于更大体系的激发态载流子动力学，有时可以在MQC的基础上更近一步，原子核一直在基态的势能面上运动，忽略激发态载流子对原子核运动的影响，即所谓的经典路径近似（classical path approximation, CPA）。这种忽略电子运动给核运动带来的影响的近似方法，在大体系计算中具备一定的合理性。一方面，在大体系中存在大量电子，单一的热电子波函数并不会对核运动产生显著的影响；另一方面，由于研究重点并非核运动而是电子动力学，这种近似能够节省相当多的计算资源，从而使得在第一性原理密度泛函理论（DFT）水平上对具有数百个原子的系统的电子动力学进行计算成为可能。由于一些历史原因，我们在之后的内容中所说的非绝热分子动力学（nonadiabatic molecular dynamics, NAMD）都是指的第三种框架。
绝大多数的NAMD方法都没有直接包含细致平衡（detail balance）和退相干（decoherence）过程。为了解决这两个问题，我们使用的是中国科学院半导体研究所的汪林望教授和北京计算科学研究中心康俊特聘研究员合作发展的一种基于密度矩阵形式的NAMD方法（P-matrix），通过密度矩阵的形式来引入退相干和细致平衡。具体来说，首先利用基态Born−Oppenheimer 分子动力学（BOMD）模拟得到核轨迹，然后将核轨迹用于含时薛定谔方程中得到波函数的含时演化，其中含时薛定谔方程采用密度矩阵形式进行重新表述。在此方法中，退相干通过密度矩阵非对角项衰减来描述。此外，这些非对角项被分为两部分，分别对应能量增加和能量降低的电子跃迁。通过将玻尔兹曼因子应用于能量增加的跃迁来实现细致平衡。
参考文献：
[1] Jun Kang,Lin-Wang Wang，Phys. Rev. B 99, 224303 (2019)
[2] Fan Zheng, Lin-Wang Wang, J. Phys. Chem. Lett. 2019, 10, 6174-6183.
02
文献案例复现
02
文献案例复现
在PWmat中NAMD部分作为BOMD的后处理来完成，教程及涉及的后处理程序安装包见官网的module 21。下面以复现文献为例展示计算流程。Module（包含后处理程序）和文献链接置于文末，点击阅读全文获取案例文件。
1. run PWmat——进行BOMD计算
(1) 准备需要的输入文件：
atom.config etot.input *.UPF slurm.sh
下面简要介绍下输入文件的内容：
atom.config（结构文件）：PtS2/GaSe 异质结
etot.input（参数文件）：
8 1
IN.ATOM = atom.config
JOB = NAMD
MD_DETAIL = 1, 1000, 1, 500, 300
NAMD_DETAIL = 100, 330, 5
in.atom = atom.config
in.psp1 = Ga.SG15.PBE.UPF
in.psp2 = Se.SG15.PBE.UPF
in.psp3 = S.SG15.PBE.UPF
in.psp4 = Pt.SG15.PBE.UPF
mp_n123 = 1 1 1 0 0 0 2
NUM_BAND = 331
ECUT = 60
VDW = DFT-D3
OUT.WG = T
#关键参数设置及含义：
设置job=NAMD，添加NAMD_DETAIL控制波函数的选取和输出情况
NAMD_DETAIL = m1 m2 nstep_out
m1 m2: 在NAMD计算中，执行传统的BOMD, 但输出绝热本征态在绝热窗口 [m1 m2] 内连续时间步长的重叠波函数。
nstep_out: 波函数 (在窗口 [m1 m2] 内) 输出到ugio.allxxxxxx 文件的步数间隔。
NAMD_DETAIL = 100, 330, 5
含义：输出100~330之间的231（m2-m1+1）个绝热态作为后续展开含时波函数的基组，并且每隔5步输出一次。
*.UPF：各元素的赝势文件，有多少种元素就包含多少个赝势文件
slurm.sh: 提交任务的脚本文件
(2) 提交计算：sbatch slurm.sh
(3) 关键的输出文件：OUT.NAMD ugio.allxxxxxx
2. 结果后处理
2.1 -update OUT.NAMD——在OUT.NAMD基础上重新选取一个绝热态范围
(1) 准备需要的输入文件：
OUT.NAMD read_outNAMD.x（后处理程序） job.sh (脚本文件,名字可以自定义)
下面简要介绍下输入文件的内容：
OUT.NAMD：job=NAMD的生成文件
read_outNAMD.x：后处理程序，从module 21的文件下载安装包中获取
job.sh脚本文件：
#!/bin/sh
#SBATCH --partition=cpu
#SBATCH --job-name=namd
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=64
module load intel/2016
export OMP_NUM_THREADS=64
export OMP_STACKSIZE=512M
./read_outNAMD.x 0 150 231 1000 231
#后处理命令相关参数的含义：
read_outNAMD.x istart iband1 iband2 nstep nstate
istart：将BOMD的第istart步当做NAMD的第1步，OUT.NAMD只包含只从MD的istart步之后输出的信息；
iband1：OUT.NAMD_update中的绝热态范围初始位置，iband1 ∈ （1 ~ nstate）
iband2：OUT.NAMD_update中的绝热态范围截止位置，iband2 ∈ （iband1 ~ nstate）
nstep：BOMD的步数
nstate：OUT.NAMD中的绝热窗口[m1,m2]中的本征态数量，330-150+1=231
read_outNAMD.x 0 150 231 1000 231
含义：从BOMD的第0步开始输出OUT.NAMD的信息，共输出1000步BOMD的信息。将其中含有231个绝热态信息的OUT.NAMD文件重新裁剪，选取其第150至第231态之间的82个绝热态信息并将其保存在OUT.NAMD_update文件中。
#注：这里是把后处理命令放在脚本文件中，是为了将后处理提交至cpu计算节点上计算（后续的后处理操作同理），最终结果等同于module 21教程中直接执行后处理命令 read_outNAMD.x
(2) 提交后处理计算：sbatch job.sh
(3) 输出文件：OUT.NAMD_update
2.2 -运行NAMD后处理程序
2.2.1-electron cooling
(1) 准备需要的输入文件：
OUT.NAMD_update NAMD.input deco_time NAMD.cc_init namd_dm.x（后处理程序） job.sh（脚本文件）
下面简要介绍下输入文件的内容：
OUT.NAMD_update: 从上一步2.1中获得
NAMD.input：后处理相关的参数文件
1000
82
1
1, 300, 0.1
1
250, 331
0d0
40
4000
100
#参数含义：
1000 ：BOMD的步数
82：OUT.NAMD_update中绝热窗口中的本征态数量
1 ：MD_DETAIL中设置的BOMD时间步长
1, 300, 0.1：空穴(-1)/电子(1) 冷却；温度，可与MD_DETAIL中的温度不同；win_Boltz (eV)，当|Ei1-Ei2|＜win_Boltz时，会使用玻尔兹曼因子来保持细致平衡
1 ：输出NAMD.cc_out内容的频率 (better to set to 1)
250, 331 ：OUT.NAMD_update中绝热窗口的范围
0d0 ：(not used)，但需要写上
40 ：(not used)，但需要写上
4000：ndt-两步BOMD间被划分为多少个dt (needs to be converged)
100：nit-两步BOMD间进行多少次对角化(needs to be converged)
deco_time：退相干时间
1 2 40
1 2 40
...
1 2 40 #82*82行，40为退相干时间，单位为fs
NAMD.cc_init：初始时刻的电子占据情况
82 #nstate
0 0
0 0
...
10 #要研究的热电子的位置
0 0
0 0
0 0
0 0
namd_dm.x：后处理程序，从module 21的文件下载安装包中获取
job.sh脚本文件：
#!/bin/sh
#SBATCH --partition=cpu
#SBATCH --job-name=namd
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=64
module load intel/2016
export OMP_NUM_THREADS=64
export OMP_STACKSIZE=512M
./namd_dm.x 0 > log.namd
(2) 提交后处理计算：sbatch job.sh
(3) 输出文件：
NAMD.graph.aveE NAMD.graph.eigen log.NAMD NAMD.cc_out NAMD.DM_out
注：
NAMDgraph.aveE：被激发的电子的能量平均值随时间的演化
NAMD.graph.eigen：所有的本征态能量随时间的演化
2.2.2-hole cooling
(1) 准备需要的输入文件：
OUT.NAMD_update NAMD.input deco_time NAMD.cc_init namd_dm.x（后处理程序） job.sh（脚本文件）
与electron cooling类似，所以下面简要介绍下有差异的输入文件内容：
NAMD.input:后处理相关的参数文件
1000
82
1
-1, 300, 0.1
1
250, 331
0d0
40
4000
100
NAMD.cc_init：初始时刻的电子占据情况
82 #nstate
1 0 #要研究的热空穴的位置
0 0
...
0 0
0 0
0 0
0 0
0 0
(2) 提交后处理计算：sbatch job.sh
(3) 输出文件:
NAMD.graph.aveE NAMD.graph.eigen log.NAMD NAMD.cc_out NAMD.DM_out
注：
NAMDgraph.aveE：激发产生的空穴的能量平均值随时间的演化
NAMD.graph.eigen：所有的本征态能量随时间的演化
3. 绘图
将NAMD.graph.eigen，电子的NAMD.graph.aveE和空穴的NAMD.graph.aveE绘制在同一张图上。
计算结果如图（1000步）:
文献中的结果（2000步）:
文献 DOI：
10.1039/d1cp02436a
Module 21（含后处理程序安装包）：
