背景导读
自20世纪80年代以来，冷冻电镜（cryo-EM）已成为生物学领域最有用的工具之一。最近，冷冻电镜也被用于材料研究。与传统的透射电镜（TEM）相比，冷冻电镜通过液氮冷冻样品，可以有效减轻电子束对样品的辐射损伤。然而，关于这种冷冻保护机制一直存在不同观点。
文章简介
关于cryo-EM中的冷冻保护机制存在多种假设，比如较低温情况下分子电离后解离的可能性更小，或者冷冻的水可以通过防止解离的碎片飞走而具有“笼”效应。近日，国防科技大学戴佳钰教授团队和中科院半导体所汪林望教授团队针对上述问题展开合作，通过实时的含时密度泛函理论分子动力学（rt-TDDFT-MD）及其新开发的自然轨道分支（NOB）算法来研究电离后的分子解离过程，以及从头算分子动力学（AIMD）研究碎片的动力学过程。研究表明，不同温度下，分子的初始解离过程是相似的；并且解离的碎片在室温下会扩散，低温时会留在同一个“笼子”里。该研究结果为“笼”效应机制提供了直接支持。相关研究成果发表在了npj Computational Materials上，论文第一作者为国防科技大学和中科院半导体所的联合培养研究生李毅。
计算方法
本文的计算采用PWmat完成。要模拟水分子环境中分子的辐射解离过程，需要解决两个难点：一是要提高非绝热的分子动力学模拟的空间和时间尺度，二是要能合理地再现分子分解过程中轨迹的随机性。PWmat的rt-TDDFT模块采用线性哈密顿量近似方法，可以有效的将模拟步长从0.1 as 量级提高到 0.1 fs 量级附近，从而允许模拟百原子体系在皮秒量级的动力学过程。而针对平均势能面方法中欠缺的细致平衡问题，退相干问题以及轨迹的随机性问题，PWmat中提供了NOB方法。NOB方法是基于密度矩阵的非绝热分子动力学的延伸，该方法的关注点从基于波函数的含时的薛定谔方程转移到了密度矩阵的含时演化方程，而自然轨道正是密度矩阵的本征矢量。相比于Boltzmann方法的rt-TDDFT，NOB方法能摆脱平均势能面的局限性，对关注重点为断键或者化学反应的多轨迹随机动力学过程有非常出色的模拟效果。想更详细了解相关算法的读者可以阅读10.1021/acs.jpca.0c06367。
主要内容
在模拟之前，首先需要构建不同温度下的平衡结构。为了合理地体现水分子环境对有机分子的影响，本研究构建了一个大小为14*14*14 A3的模型，其中包含1个碳酸乙烯酯（EC）分子和88个水分子。
该结构如图1（a）所示，水分子聚集在EC分子周围，密度约0.94g/cm3，是实验上冷冻样品中低密度无定形冰（LDA）的密度。为了获得合理的无序结构，先将该体系加热到1000 K（在这个过程中没有水分子断键），随后将无序结构分别退火到100 K（冷冻条件，接近液氮温度）和300 K（室温条件）。这样获得的结构的RDFs与Alexander等人的文章结果是一致的（图1（b））。
图 1（a）冷冻条件下模拟样品的结构；红色和白色部分表示水分子，棕色和黄色部分表示EC分子。（b）O原子之间的RDFs的比较。浅蓝色和红色实线表示本研究的结果，而星星和圆点是Alexander等人的结果。（c）电中性EC分子和整个模拟体系的DFT本征能量谱。（d）EC分子的最低能量本征态中的一些的波函数（编号为1、2、3、92、93和94），而那些较高能量本征态（C和O原子的p轨道电子）倾向于与水混合。
01.排除温度对初始分子解离过程的影响
建立平衡结构之后，通过rt-TDDFT-MD计算来模拟分子解离过程。由于电子束辐照使一个小分子内产生多次激发的可能性较小，所以该研究只考虑单电子电离情况。
本研究在每个温度下进行了21次NOB模拟，每次模拟的初始条件和时间相同。图2是4次NOB模拟的结果，其中（a）-（d）处于冷冻条件下，（e）-（h）处于室温下。模拟开始时，位于价带的state-1被一个空穴占据。随着空穴跃迁到不同的能级，键长开始变化，4次模拟中均有三种键长显著增加（表明这种键的断裂）：C1-O2，C1-O3 和 C2-C3（图 2（a）、（c）、（e）和（g））。其中，4次模拟中的断键时间和断键模式不同，证明了模拟的随机性。并且，当空穴从state-1跃迁到state-2 或state-3时，键长没有显著变化。直到跃迁到state-92，键长才开始急剧变化。上述结果表明，断键总是由WFC引起的，且通常发生在两个能量相差比较大的深能级之间。
图 2 左图（a，c，e，g）显示了 EC 分子中键长的变化。在100 K 下（a）和（c）随机数不同，在 300 K 下（e）和（g）随机数不同，导致不同的坍缩。EC 分子模型及其原子编号如（c）所示。右图（b，d，f，h）表示整个系统的绝热态电子占据和能级随时间的变化。线表示能级的变化，红点表示空穴的占据，点的大小与占据的空穴数成正比。
总结所有42次模拟，结果如图3（a）所示。体系的能级可以分为四个区域。其中，空穴从第一个区域跃迁到第二个区域几乎没有键断裂；跃迁到第四个区域，由于大的能量差导致结构发生很大变化，不确定哪些键会发生变化；跃迁到第三个区域的不同能级上，产生不同类型的键断裂。具体来说，两个温度下，空穴跃迁到其他三个区域的概率相近（图3（a）），产生断键类型相同，而且分解产生的不同类型碎片概率也相近（图3（b））。这说明，两个温度下体系具有非常相似的分子解离模式，即温度对于电子束辐照使分子电离导致解离过程的影响是差不多的
图 3（a）所有模拟中空穴跃迁的路径和分子解离的模式; 红点代表空穴的占据。左图列出了在两个温度下空穴从第一个区域跃迁到另一个区域的概率。右图显示了由于空穴跃迁到第三个区域而产生的分子碎片图像。（b）碎片的类型和发生键断裂的百分比。“Other”是指不太常见的片段。“Break”表示在200-fs模拟期间至少有一个键断裂的概率。
02.确定冷冻保护机制是“笼”效应
NOB模拟期间，空穴跃迁产生的能量会转移到分子碎片中，因此，碎片具有很大的动能。为了探究碎片的动力学过程，将rt-TDDFT-MD的最终原子位置和速度作为初始条件，进行了2 ps的AIMD模拟。
首先以单个“O”原子碎片的演变来看，在300 K下，“O”原子碎片向远离原始EC位置的方向移动，容易与水分子结合形成水合物，甚至将水分子的H原子吸附到自己身上；在100 K， “O”原子碎片与EC分子和其它碎片保持很近的距离，可以与其它碎片反应形成新的碎片，如图4（a）和（b）所示。
“笼”效应也发生在其它碎片上。例如，在100 K时，碎片"C2-O2"和"O1- C1-O3-C1"将结合形成长链，且具有类似于原始EC的弯曲形状，如图4（c）、（d）和（e）所示。但是，在2 ps AIMD中，长链并没有完全重组成EC。为了估计长链完全重组成EC的可能性，使用NEB方法计算了形成EC所需的反应势垒约1.5 eV。在2 ps AIMD结束时，分子具有800 K的温度。如果该温度保持恒定，则重组为EC的反应可以在大约0.003秒发生。
图 4 碎片的动力学过程。（a）（b）AIMD模拟中单个“O”原子碎片的路径，（a）发生在100 K，（b）发生在300 K；（c）（d）（e）EC分子在冷冻和室温条件下电离2 ps后的动力学。（c）电子束辐射导致键断裂，断键后如（d）；（e）显示了两种温度条件的差异，其中碎片在冷冻过程中由于“笼”效应而重组，在室温下没有“笼”效应而飞散。黄色原子表示EC分子中的O原子，棕色原子表示EC中的C原子，红色原子表示水分子中的O原子，白色原子表示水分子中的氢原子。EB是电子束（electron beam）的缩写。
综上所述，分子在初始电离后的解离过程几乎是相同的。在初始解离过程之后，低温下冻结的水分子可以防止碎片飞散和与水分子反应（如在高温情况下）。这种“笼”效应将碎片保持在一起，从而促进它们的修复反应。因此，本研究为冷冻电镜中的“笼”效应提供了理论支持。
