基于经验势分子动力学模拟的N719染料/二氧化钛界面在室温离子液体溶剂中的振动光谱

摘要

采用室温离子液体（RTILs）对染料敏化太阳能电池进行溶剂化处理；通过优化的经验势函数，应用分子动力学模拟计算其振动特性。将所得振动光谱与实验结果及从头算分子动力学结果进行比较；多种经验势函数的光谱表明，离子液体的部分电荷参数化如何影响振动光谱的预测。

摘要

在物理化学领域，对与光吸收染料接触的光敏金属氧化物表面的振动光谱及其他结构、能量和光谱特性进行精确的分子模拟预测，仍是一项持续存在且棘手的挑战。基于此，本研究采用优化的经验势函数，对一种具有代表性且典型的染料敏化太阳能电池（DSC）体系进行了分子动力学（MD）模拟。该体系以广泛研究的室温离子液体（RTIL）作为溶剂，具体为[bmim]⁺[NTf2]^- RTIL，其中N719敏化染料吸附在锐钛矿相二氧化钛（101）表面。通过该模拟，深入揭示了以RTIL作为电解质空穴受体时，如何调控N719染料的动力学与振动特性，并通过MD得到的质量加权速度自相关函数的傅里叶变换，估算了DSC光敏界面的光谱。所获得的振动光谱与实验光谱以及从从头算分子动力学（AIMD）采样得到的光谱进行了比较；特别是，由MD生成的多种经验势光谱揭示了离子液体的部分电荷参数化方式对振动光谱预测的影响。总之，当经过AIMD和实验验证后，经验力场模型的精细拟合已被证明是研究DSC振动特性的有效工具。

引言

在染料敏化太阳能电池（DSCs）中，半导体的光学带隙由一种可吸收光的敏化染料桥接。DSCs需要持续再充电，因此必须使用氧化还原电解质以维持持续的电荷供应（通常为有机溶剂中的I^-/I^3-体系）。该过程促使空穴从敏化染料传递至电解质，同时注入的光生电子进入金属氧化物基底，并通过外电路传输，最终在阴极发生复合¹。DSCs之所以在多种实际应用中具有广阔前景，关键在于其制造工艺简单，无需高纯度原材料；这与硅基光伏器件所需的高资本投入和超高纯度形成鲜明对比。无论如何，通过用具有低挥发性的室温离子液体（RTILs）替代稳定性较差的电解质，有望显著延长DSCs的工作寿命，展现出巨大潜力。RTILs兼具类固态的物理特性与类液态的电学特性（如低毒性、低可燃性和低挥发性）¹，使其成为DSC应用中极为理想的候选电解质。

鉴于RTILs在DSC中的应用前景，近年来针对DSC原型N719染料/二氧化钛界面与RTILs体系的研究活动显著增加，这并不令人意外。特别是已有针对此类系统的重要研究工作^2,3,4,5，这些研究涵盖了一系列广泛的物理化学过程，包括染料中的电荷补充动力学^2,5、电子-空穴动力学与转移的机理步骤³，以及二氧化钛基底的纳米尺度特性对这些及其他过程的影响⁴。

现在，考虑到基于密度泛函理论（DFT）的分子模拟，特别是从头算分子动力学（AIMD）领域的显著进展⁶，作为材料科学中一种极为有用的原型设计工具，尤其适用于染料敏化太阳能电池（DSCs）^7,8,9,10,11，其中对最佳功能选择的严格评估至关重要^8,9，AIMD 方法此前在深入研究染料结构、吸附模式及在染料敏化太阳能电池半导体表面的振动特性所受的显著色散效应和显式离子液体溶剂化效应方面已证明非常有效。特别是，采用 AIMD 已在合理地半定量捕捉和预测重要电子性质（如带隙）以及结构结合方面取得一定成功¹³ 和振动光谱¹⁴ 在参考文献12-14中，对光敏N719染料分子吸附于锐钛矿相二氧化钛(101)表面的体系进行了广泛的AIMD模拟，评估了在含有[bmim]条件下体系的电子性质和结构性质⁺[NTf2]^{- 12,13} 和 [bmim]⁺[I]^{- 14} RTILs，以及[bmim]情况下的振动光谱⁺[I]^{- 14}特别是半导体表面的刚性¹⁵，除了其固有的光催化活性外，还在AIMD模拟中导致表面发生轻微变化，从而使(101)锐钛矿界面^12,13,14 一个合适的选择。如参考文献12所示，阳离子与表面之间的平均距离减少了约0.5 Å，阳离子与阴离子之间的平均间距减少了0.6 Å，染料周围第一层中的室温离子液体（RTILs）发生了显著改变，其中阳离子与染料中心的平均距离增加了1.5 Å，这些现象均由RTIL溶剂化体系中显式的色散相互作用直接导致。在真空中引入显式色散效应也会导致吸附态N719染料构型出现非物理性的弯曲。参考文献13分析了显式RTIL溶剂化及泛函选择对染料敏化太阳能电池（DSSCs）行为的影响，结论指出显式溶剂化和色散作用的处理均至关重要。参考文献14在掌握其他课题组高质量实验振动光谱数据的基础上，系统地对显式[bmim]⁺[I]^- 溶剂化作用及色散效应的精确处理，见参考文献12 & 13 对显著谱模特征的再现；这得出结论：显式溶剂化效应与色散相互作用的精确处理同样重要，此结论呼应了此前在显式溶剂中催化剂的 AIMD 建模中关于结构与动力学性质的发现¹⁶事实上，Mosconi 等人也对显式溶剂化效应在DFT处理DSC模拟中的影响进行了出色的评估。¹⁷. Bahers 等¹⁸ 研究了染料的实验吸收光谱及其在TD-DFT水平下的相关光谱；这些TD-DFT计算得到的光谱在跃迁特性方面与其实验光谱高度一致。此外，Preat等人在多种溶剂中研究了吡咯烷（PYR）衍生物的吸收光谱¹⁹，为染料的几何与电子结构提供了重要见解，并揭示了适当的结构修饰可优化基于PYR的染料敏化太阳能电池（DSSCs）的性能——这正是模拟驱动/理性化的“分子设计”理念的体现。

在明确确立了密度泛函理论（DFT）和从头算分子动力学（AIMD）在准确模拟染料敏化太阳能电池（DSC）的性质与功能方面的重要贡献之后——包括从结构、电子和振动角度对显式溶剂化以及色散相互作用的恰当处理等关键技术问题^{7,8,9,10,11,12,13,14}，本研究的关注点转向一个更具实践意义的问题：基于经验势的方法在多大程度上能够被优化，以合理而准确地预测此类典型DSC体系的结构与振动性质。本文以N719染料吸附在锐钛矿（101）表面并处于[bmim]⁺[NTf2]^-室温离子液体（RTIL）环境中的体系为例进行探讨。这一问题至关重要，不仅因为已有大量基于力场的分子模拟方法和理论工具可用于DSC的模拟研究⁷，以及更广泛地应用于金属氧化物表面体系，还因为相较于基于DFT的方法，经验势方法具有显著降低的计算成本，同时具备与偏置采样方法高效耦合的可能性，从而更有效地捕捉在常温下表现出类固态物理特性的高黏性RTIL溶剂中的相空间行为和结构演化。因此，受这一关于如何评估并优化力场方法的开放性问题所驱动，并结合DFT、AIMD以及实验振动光谱数据的指导¹⁴，我们着手开展一项紧迫任务：通过分子动力学（MD）模拟中对N719染料原子速度自相关函数（VACF）进行质量加权傅里叶变换，评估经验势在振动光谱预测方面的性能表现。其中一个关键问题是RTIL的不同部分电荷参数化方案如何影响振动光谱的预测结果，对此我们给予了特别关注，同时也系统考察了相对于实验和AIMD结果优化力场以实现最佳谱峰模式预测的更广泛目标²⁰。

方案

1. 使用 DL_POLY 进行分子动力学模拟

1. 构建N719染料吸附于溶剂化[ bmim ]+[ NTf2 ]环境中的锐钛矿型二氧化钛（101）表面的DSC系统的初始结构^- 来自先前的研究工作^12,13使用 VESTA 软件绘制所需结构。
2. 选择无抗衡离子的N719顺式-二(硫氰酸根)-双(2,2'-联吡啶-4-羧酸-4'-羧基)-钌(II)敏化染料，并确保其含有两个表面结合的质子，以维持整个系统的电荷中性。
   注意：事实上，De Angelis 等人的深入研究已证实，这构成了N719在锐钛矿型二氧化钛上吸附的真实表征²¹. 这是因为在参考文献21中，多个性质的理论结果与实验结果符合得最为令人信服；在实验体系中，这些 事实上的 表面质子化被认为源于离子液体的阳离子和阴离子，导致表面发生一定程度的电荷转移²¹.
3. 确保染料在化学上吸附到 TiO₂ 上₂ 通过两个羧酸根（即双齿配位）吸附在表面。这种初始的染料吸附构型与Schiffman等人所描述的I1构型相似。²²，该结果在考虑表面质子化的情况下被确定为最稳定的状态。参见参考文献12 & 13，了解具体操作方法，包括化学吸附坐标。
4. 确保敏化染料 N719（顺式-二（硫氰酸根）-双（2,2'-联吡啶-4-羧酸-4'-羧酸）-钌(II)）不含任何抗衡离子。为保持电荷中性，添加两个表面结合的质子，如参考文献12所示。 & 13.
5. 选择12对阳离子-阴离子组合的1-丁基-3-甲基咪唑双(三氟甲基磺酰)亚胺，共包含480个原子^12,13这些数据来自参考文献12。 & 13.
6. 使用 Lopes 等人经过充分验证的力场，通过经验势对 RTIL 构型进行弛豫。²³使用 Matsui-Akaogi (MA) 力场构建锐钛矿模型，并在弛豫过程中考虑二氧化钛的可移动性。根据下文步骤 2.1 中的 DL-POLY 参数，在 DL-POLY 中执行几何优化（而非分子动力学模拟），采用共轭梯度最小化方法，相对终止梯度设为 0.001。在此过程中，需在 FIELD 文件中指定优化计算，而非动力学模拟。
7. 对于锐钛矿表面，(TiO₂)₉₆由288个原子组成，确保其沿实验室坐标系的x轴和y轴方向具有周期性，从而为RTIL提供一对平行的(101)表面；x轴方向尺寸为23 Å，y轴方向为21 Å。该数据引自参考文献12。 & 13.
8. 确保包含显式溶剂的整个DSC系统由827个原子组成^12,13对于“真空”情况，即无RTIL溶剂化时，体系中应包含347个原子。

2. 使用 DL_POLY 进行基于力场的分子动力学模拟

1. 使用 DL-POLY 进行分子动力学模拟，采用多种不同的部分电荷集（见下文），在 NVT 系综下以 1 fs 时间步长运行 15 ps，温度保持在 300 K^24,25，使用 Lopes 等人的方法²³ RTIL 的力场参数和染料的通用 OPLS 模型²⁶，采用经过深入研究且可靠的Matsui-Akaogi势函数描述二氧化钛的相互作用²⁷ 力场。要在终端上运行 DL-POLY，请输入 DLPOLY.X & 输入文件所在的位置。
2. 通过 DL_POLY 中实现的上述指定经验力场进行经典分子动力学模拟²⁸此处无需使用软件中的图形用户界面（GUI），建议根据该软件详尽且易于遵循的手册输入相关参数²⁹在 CONTROL 文件中（输入文件详见补充信息），为 NVT 选择“Nose-Hoover”，并设置每 1 fs 输出一次位置-速度轨迹。
3. 在 FIELD 文件中，对于 Lennard-Jones 参数，采用 Lorentz-Berthelot 混合规则²⁵取Lennard-Jones（LJ）半径的算术平均值和LJ势阱深度的几何平均值（针对经验力场，详见参考文献25），并将其填入FIELD文件非键相互作用选项卡的底部区域。
4. 为处理长程静电相互作用，采用 Ewald 方法²⁵；使用非键截断长度 r_切割 = 10 Å。参见参考文献25。 & 30，以获取有关如何优化静电参数的精确细节。在 CONTROL 文件中，将 Ewald 方法的实空间衰减参数设置为约 3.14/r_切割，并选择 Ewald 波矢的数量以确保 Ewald 求和计算的相对容差为 1E-5；在 CONTROL 文件中指定该参数。
5. 在 CONTROL 文件中，确保声明 r_切割 = 10 Å；使用 REVCON 文件（重命名为 CONFIG）进行一系列势能计算，直至 OUTPUT 中的系统压力在百分之几范围内收敛，以确定 r_切割，但避免任何 r_切割 小于约2.5倍的最大LJ距离^25,30.
   注意：此15 ps的MD传播时间尺度的选择与参考文献9、10中具有相同初始构型的约8.5 ps的玻恩-奥本海默分子动力学（BOMD）模拟相近。 & 17，以便直接比较两种人工智能方法提供的振动光谱预测结果²⁰ 以及基于力场的分子动力学模拟（同时与实验结果进行比较和验证）。
6. 从 HISTORY 文件（根据 CONTROL 文件的指令，每个时间步长的速度和位置均已输出至该文件）中，使用 Python 提取 x、y、z 方向的速度 dye_atom_velocity_seperate.py （参见补充信息）在终端中执行，它将分离每一步的速度。
7. 使用以下方法计算速度自相关函数（VACF） vacf151005.py （参见补充信息）。在终端中输入 ./classical_dye_autocorr.sh；它将计算所有染料原子的VACF。从分子动力学（无论是AIMD）计算光谱^14,20 或基于力场）通过染料原子速度自相关函数（VACF）的质量加权傅里叶变换^31,32,33 通过使用 Python MWPS.py （参见补充信息）。在终端中输入 ./run_all_4.sh；它将计算质量加权的功率谱。
8. 使用常用软件对这些速度自相关函数（VACFs）进行傅里叶变换。
9. 请注意，为了与实验数据进行基准比较，采用最高质量的DFT处理方法（例如，使用显式溶剂化模型、精确处理色散作用，并谨慎选择泛函）在AIMD中至关重要^12,13,20 评估和调整经验势分子动力学（empirical-potential MD）的相对性能，特别是静电作用及部分电荷选取所造成的显著影响¹⁶参见参考文献12和13，深入研读以充分理解；若未来需进行AIMD（本研究中未涉及），应在必要时据此采取相应措施。

3. 比较每种力场的结果

注意：在步骤 2 中，对于基于经验势的分子动力学模拟，评估离子液体（RTIL）的部分电荷集非常重要，以便能够相互之间、与实验结果以及在显式 RTIL 溶剂中进行的从头算分子动力学（ab initio-MD）结果进行直接比较（采用 PBE 泛函结合 Grimme-D3 色散校正，因其在振动光谱预测方面表现出更优的性能）²⁰；这些部分电荷集如下：

1. 请注意，对于文献推导的RTIL电荷，阴离子电荷需通过扩展休克尔理论获得^34,35，基于AIMD轨迹²³，由于参考文献20中缺乏阴离子电荷参数化，阳离子电荷需取自Lopes等人。²³ 准备一份文献电荷表，并将其转换为 DL-POLY 的 FIELD 文件格式。
2. 请注意，Mulliken RTIL 电荷需通过 Mulliken 布居分析计算。通过对从头算分子动力学轨迹中的四个时间点取平均值来进行 Mulliken 分析。²⁰，重新归一化并制备文献电荷表，然后转换为 DL-POLY 的 FIELD 文件格式。
3. 请注意，扩展休克尔理论（EHT）电荷需根据 AIMD 轨迹的最终构型进行拟合²⁰，使用扩展休克尔理论（EHT），适用于RTIL的阴离子和阳离子。通过对从头算分子动力学轨迹的四个时间点取平均，进行EHT分析。²⁰在MOE软件包中，通过读取构型文件后选择“电荷分析”菜单来实现³⁵，重归一化并制备文献电荷表，然后转换为 DL-POLY 的 FIELD 文件格式。
4. 请注意，Hirshfeld RTIL 电荷应通过对从头算分子动力学轨迹中的四个时间点进行Hirshfeld电荷分析并取平均值得到²⁰，用于阴离子和阳离子²⁰在MOE软件包中，通过读取构型文件后选择“电荷分析”菜单来实现³⁵根据对这些获得的电荷进行重整化，以适当格式将其列入 DL-POLY FIELD 文件中。
5. 请注意，[bmim] 的不同电荷集⁺[NTf2]^- 原子以某种方式呈现 表1 & 表2，同时也显示了为考虑对称性和总电荷守恒而需要调整的部分原子电荷的平均变化量。
6. 请注意，最终的电荷设置应使阳离子的总电荷为+1，阴离子的总电荷为-1。阳离子和阴离子如图所示 图1a 和 图1b分别为。这些电荷集所依据的DFT采样光谱在很大程度上包含隐含的电荷转移，往往导致电荷值更接近±1。

结果

结合基序的结构特性
在分子动力学模拟进行15 ps后，四种不同部分电荷集的代表性结合基序如图2所示。在图2a中，对于上述文献推导的电荷，可以看出存在与表面质子显著的氢键相互作用。通过对轨迹的仔细分析发现，氢键主要与表面质子结合，而其余三种（基于AIMD的）²⁰电荷集则未表现出与表面质子如此强烈的库仑相互作用。参考文献20中的图1f（描绘了约8.5 ps后的AIMD弛豫染料结合构型），与表面质子形成强氢键的证据也较少，因此在定性上与本研究中基于力场的分子动力学模拟结果一致——即采用基于AIMD的电荷集可获得定性相似的底物结合构型（图2b-d）。事实上...

讨论

从头算模拟技术的计算成本较高，因此若要在更长的时间尺度上进行模拟，则至少需要对DSC系统的部分组件采用经验力场。为此，我们使用一种经验性的经典模拟力场，建立了溶剂化界面[bmim]⁺[NTf2]^-的等效原子级模型用于分子动力学（MD）模拟。其中，锐钛矿相采用Matsui-Akaogi（MA）力场建模，染料结构则使用OPLS参数处理。对于室温离子液体（RTIL），采用了四种不同的力场以评估RTIL电荷对染料振动动力学的影响。各力场的参数保持一致，仅改变电荷参数化方案，力场参数通过扩展Huckel理论（EHT）的自动拟合程序生成，该程序在MOE软件包中实现，并基于从头算轨迹进行拟合。

从头算结果在 300-400、625、750、825、1000、1300、1400、1500、1525、1575、1600、1700 和 2075 cm^-1 处显示出一组光谱峰。根据与从头算光谱的接近程度，经典电荷集按从近到远的顺序排列为：EHT、Hirshfeld、文献值和 Mulliken。在经典光谱中，基于文献导出电荷的主要峰位于 6...

材料

本文使用的材料清单

姓名	公司	目录编号	评论
这是一项分子模拟研究，因此未使用任何实验设备。
所用软件名为 DL-POLY（其“Classic”版本可通过 sourceforge 在 GnuPublic 许可下获取）