一种寻找团簇异构体的准动力学方法.doc

一种寻找团簇异构体的准动力学方法
——稳定团簇异构体能谱研究结题论文

摘 要：本文建立了一套普遍适用的寻找团簇异构体的准动力学方法，该方法能够迅速给出在一般气相生长条件下形成几率较大的异构体。我们用该方法得到了C21的异构体谱，并采用分子动力学方法模拟了21个自由的碳原子在氦气氛中形成稳定团簇的过程，表明动力学过程中形成几率较大的异构体都已包含在该方法所得到的异构体谱中。我们所得到的C21最稳定结构的势能远低于采用遗传算法所得到的结果[Chem. Phys. Lett. 364 213，2002]。
关键词：团簇，异构体，最低势能，全局优化，时间倒流

1．引言

团簇的物理和化学性质随团簇所含原子数目及团簇构型的不同而改变，而实验观察很难直接给出团簇的原子构型，于是理论研究便成了人们获知团簇异构体构型及最稳定构型的重要方法。从几何角度看，团簇异构体数量随着团簇所含原子数目呈指数增长[1]，理论上已不可能得到含有几十个原子的团簇的全部异构体。因此，如何找到那些在实际生长过程中出现几率较大的异构体，尤其是找到具有最低势能的异构体是团簇研究的重要问题之一。
一种有效的方法是先从几何构型上将异构体进行分类，再重点研究可能含有最低势能构型的整个异构体族。例如在以前关于富勒烯的研究中，人们只关注仅含五元环和六元环的经典异构体（经典富勒烯）[2]。但是，我们在C60的生长动力学研究中发现[3]，一些含有七元环缺陷的异构体在动力学生长过程中占据十分重要的地位，其出现几率甚至超过许多经典异构体。值得注意的是，含有七元环缺陷的异构体已超过数百万个[4]，对其逐个进行研究将是十分困难的；如果再考虑其它缺陷（如四元环缺陷或不成笼等）的异构体，将面临更大的困难。而事实上在一般实验条件下，只有为数不多的异构体以可观的几率形成。
为了避免涉及实际生长过程中很少出现的异构体，人们发展了很多全局优化的方法[1, 5-11]以期寻找低势能段的异构体及其最稳定构型。诸多算法中，模拟退火算法[1]和遗传算法[6]得以广泛应用。前者使团簇原子沿着类似“固体退火”的轨迹迅速到达势能面的一个势谷，产生一种异构体，其局部搜索能力很强，而全局搜索能力却较差；后者则按照“生物自然选择”的方式进行全局搜索，但是局部搜索能力却稍显不足。
本文建立了一套用于寻找团簇异构体及最稳定构型的时间倒流准动力学方法。在该方法中，原子按照接近实际的“团簇动力学生长轨迹”快速降温演化，并形成动力学中可几的异构体；同时，引入时间倒流的思想[12]，使到达势能极小值的团簇可以从势谷中跳出，向更低势能演化，兼顾了搜索的全局性和局域性。我们使用时间倒流准动力学方法给出了C21团簇的异构体势能谱（E-chart），同时获得了C21团簇的最低势能构型，其势能远低于Zhang等[13]采用遗传算法找到的最低势能。

2．理论模型

考虑到气相条件下团簇的生长过程可分为三个阶段，即高温环境中自由原子的产生，缓冲气体中原子的冷却以及最后稳定团簇的形成，我们建立如下的准动力学模型：在一个施加了周期性边界条件的立方体盒子中随机的放置一定数量的孤立原子，并初始赋予它们高温（显著高于熔点）条件下的麦克斯韦速度；在随后的演化过程中，我们每隔一定时间Δt随机的选择一个原子将其速度用一个新的速度替代[14,15]。这里
i =x，y，z (1)
其中是随机选取的室温（300 K）麦克斯韦速度，是介于0,1之间的参数。整个体系在演化过程中逐渐冷却，当其温度到达300 K后我们采用阻尼轨线方法[16]使体系迅速弛豫到0 K，从而得到团簇一个异构体的记录（原子构型和势能）。需要说明的是，为了提高计算效率，立方体盒子的线径应尽可能小，但应能够保证原子之间的初始距离显著大于原子间的作用力程。
重复运行上述准动力学模型，可以得到若干个团簇异构体。但一般说来它们只对应于势能面上的一些极小点，我们仍无法确定是否已经得到了具有最低势能的异构体。如果采用Γ相空间中的一个代表点表示上述体系的一个状态，则体系在温度降低的过程中向着动量p减小的方向沿某一轨迹运行到某一势能极小值。由于粒子在运动过程中受到随机作用力的影响，系统由相空间同一点P1出发可能到达不同的终点C1, C2，而同一终点C2可能对应许多不同的出发点P1，P2，从而在相空间中形成复杂交错的轨迹，如图1所示。我们引入的时间倒流方法的基本思想是[12]，选取相空间中的某一终点Ci（对应于势能面上一极小值点），使‘时间倒流’，系统将从Ci退回到‘过去’某一点Pi；然后再使‘时间顺流’，这时系统将重新向势能极小值演化。由于代表点在相空间中的轨迹是交错的，系统经过一次循环后不一定再回到终点Ci，而可能到达另一终点Cj。事实上，如果我们将势能面按高势垒进行划分，每两个高势垒间的区域回流到相空间中的等温面可能经历的路径形成一个分区，这些分区有一定的交集，而这些交集连通了整个势能面。在时间倒流顺流的循环过程中，体系可以以一定的几率回流到原分区的极小值，加强了该分区的局部搜索能力；与此同时，体系又有一定的几率通过“交集”到达其他分区的极小值（C1-P1-C2-P2-C3），从而加强了全局搜索能力。

图1 体系在Γ相空间中的轨迹示意图（贴近等温线的相空间“交集”连通整个势能面）
搜寻团簇异构体的程序分为两步。首先，应用准动力学模型得到若干个异构体，从中选取一个具有最低势能的构型；然后，赋予该构型300 K下的麦克斯韦速度，并使该体系多次经历时间倒流顺流循环，从而产生许多新的异构体。在此过程中，如果产生了具有更低势能的异构体，它将作为下一次倒流顺流循环的初始构型。
需要说明的是，随着时间倒流的进行，体系的温度将不断上升。我们把回流的最高温度设置为T（接近于相应体材料的熔点）以控制体系倒流的‘距离’。如果在进行了N次循环后仍没有产生势能更低的构型，便认为在（T，N）水平上获得了最稳定构型。我们的计算表明，在（3000，300）水平上能够迅速找到C18~C36的最低势能构型[17]。
值得进一步解释的是，采用时间倒流方法使体系回流到温度T与直接给体系赋同样温度的麦克斯韦速度是显著不同的。首先，直接给体系赋予较高的温度往往会造成体系的严重混乱，甚至可能将大团簇“蒸发”成散开的小团簇，从而大大影响了局部优化的效果；其次，直接给体系赋T K的温度，使得体系失去了在回流过程中选择其他路径的可能，而其降温演化过程无异于一次模拟退火，因此体系很难跳出较高的势垒，从而影响到全局搜索的效果。

3．C21异构体谱

在应用上述理论模型研究C21异构体的过程中，立方体盒子的边长取为2.5 nm，该线径能够保证21个碳原子之间的初始距离（0.8 nm）显著大于原子间的力程；碳原子之间的相互作用势由Brenner函数[18]描述，其经典轨迹采用标准的Verlet算法计算，积分步长取为0.2 fs；与(1)式相关的参数和Δt分别取为0.1和1.3 ps.
应用准动力学模型时，体系初始温度设置为5000 K，经过50次计算后产生了50个C21团簇，其中34个是不同的异构体。我们选取其中的最低势能构型作为出发点，进行了(3000, 300)水平下的时间倒流顺流循环计算。剔除计算中重复出现的异构体，其余不同异构体平均每原子的势能（PEPA）分布如图2所示，我们称其为E-chart。图3中的a-g显示E-chart中势能最低的七个构型。

图2 准动力学时间倒流方法得到的C21的E-chart

图3 C21团簇的部分异构体（a-g.准动力学方法得到的构型 h.构造的符合IPR规则的构型 i.构造的笼状构型 j. 遗传算法得到的最低势能构型[13]）

为了确认该方法是否已找到最低势能构型和动力学中可几的异构体，我们进行了如下的分子动力学模拟[3, 19-22]：在一个边长为7.2 nm并施加了周期性边界条件的立方体盒子里，随机放置21个碳原子和100个氦原子（相应的室温气压为1.01×106 Pa），碳原子之间的相互作用仍采用Brenner函数描述，而碳-氦作用和氦-氦作用则采用Lennard-Jones势函数描述，其中的参数根据Lorentz-Berthelot混合规则[23]确定。初始时，碳原子和氦原子被分别赋予5000和2500 K下的麦克斯韦速度。随着体系按照经典力学轨迹的演化，氦原子平均温度将不断升高，当其温度高于3000 K时给氦原子重置2500 K下的麦克斯韦速度，以使体系温度大致保持在2500 K。当碳原子的平均温度达到环境温度2500 K后，体系继续驰豫200 ns，期间每隔2 ps采集一个团簇‘样品’。然后，采用阻尼轨线的方法将每一‘样品’迅速冷却到0 K，由此给出一个异构体的记录（原子构型及相应的势能）。关于异构体的全部记录表明，上述动力学模拟中没有形成比E-chart(图2)中势能更低的构型，而出现几率大于1%（如图4所示）的12个异构体都已包含于C21的E-chart中。

图4 C21动力学演化过程中出现概率较大的异构体几率

上述关于各种异构体相对几率的统计，仅是体系达到平衡后200ns内的统计数据，我们尚不能确定各种异构体的相对几率是否会随体系的进一步演化而显著变化。为此，我们将图3.a所示的最稳定构型作为初始构型，使含有氦的体系在2500 K重新驰豫200 ns，并对0 K时的异构体记录进行统计，结果给出了与图4类似的几率分布。因此，可以将初始处于高温状态的碳原子体系与环境达到平衡后驰豫200 ns的统计结果近似视为无限长时间的驰豫结果。
进一步比较动力学模拟中记录的构型和E-chart（图2）中的构型发现，E-chart中的一些异构体在动力学模拟过程中没有出现；在我们更大碳团簇的研究中这一现象更为明显，甚至E-chart中最低的十几个异构体都没有在其2500 K下的动力学模拟中出现[17]。我们知道，团簇的生长动力学受到很多热力学因素的影响，尤其是势能最低的异构体在不同环境条件下不一定同样具有最低的自由能[24]，因此某种实验条件的动力学模拟只能给出相应条件下自由能较低的异构体而常遗漏最稳定构型。另一方面，动力学模拟会产生许多在E-chart中不存在的异构体，他们对应于动力学演化过程中的过渡态构型，因此是不稳定的。
我们注意到，Zhang等人在2002年用Brenner势和遗传算法找到的最低势能构型如图3.j所示[13]，其势能（-6.39357 eV）位于我们所得到的E-chart的第14与15位之间，而在动力学模拟中该构型出现的几率仅为0.0027%。此外，我们还构造了E-chart(图2)之外的两种构型，分别示于图3.h和图3.i。前者可被视为C60（Ih）的帽子，且符合IPR规则[25]，似乎应为最稳定构型；而后者是由五元环和六元环构成的笼，也应具有较低势能。然而，图3.h构型势能（-6.42004 eV）高于E-chart中最低的两个异构体能量，而图3.i构型势能(-6.44421 eV)略高于E-chart中最稳定构型势能(-6.44511 eV)。有趣的是，在动力学模拟中C60帽子出现的几率仅为万分之六，而笼状结构没有出现。我们还采用广义Stone-Wales变换(GSW)[4]得到了C21的若干个异构体。将这些异构体的势能与E-chart比较（图5），易见GSW变换异构体的势能普遍偏高。而我们模拟C21团簇在氦气氛中生长动力学的数据表明，这些异构体在动力学演化过程中很少出现。由此可见，我们的时间倒流准动力学方法找到了最低势能构型，但一般不产生在通常气相生长条件下不可几的异构体。

图5 GSW得到的势能谱和E-chart的比较
4．结论
本文建立的寻找团簇异构体的方法，能够快速有效的获得气相生长条件下形成几率较大的异构体势能谱及最低势能构型。由于目前大量团簇粒子是在气相条件下生成的，因此我们建立的方法可广泛应用于相关的研究。文中所给出的异构体谱对于进一步应用第一性原理研究C21异构体构型及其物理化学特性有指导意义。

参考文献：
1. S.Kirkpatrick et al 1983 Science 220,671
2. D.E. Manolopoulos et al 1992 J. Chem. Phys. 96 7603
3. P. Li et al 2004 J. Chem. Phys. 121 7701
4. 李鹏 2006 复旦大学博士学位论文 p80-p86
5. N. Metropolis et al 1953 J. Chem. Phys. 21 1087
6. D.M.Deaven et al 1995 Phys. Rev. Lett. 75 288
7. D.J.Wales et al 1999 Science 285 1368
8. D.L.Freeman et al 1985 J. Chem. Phys. 82 462
9. J. Lee et al 1997 J. Comp. Chem. 18 1122
10. T. Zhou et al 2005 Phys. Rev. E 72 016702
11. W.S. Cai et al 2002 J. Comp. Chem. 23 427
12. X.J. Ning et al 1999 J. Chem. Phys. 110 4920
13. C. Zhang et al 2002 Chem. Phys. Lett. 364 213
14. M.E. Riley et al 1988 J. Chem. Phys. 88 5934
15. Y. Sheng et al 2004 Chin. Phys. Lett. 21 2012
16. L.M. Raff et al 1991 J. Chem. Phys. 95 8901
17. W. Wei, D.J. Han, X.J. Ning et al , unpublished data
18. D.W. Brenner 1990 Phys. Rev. B 42 9458
19. Y. Sheng et al 2004 J. Chem. Phys. 121 2013
20. Sheng Y. et al 2004 Acta Phys. Sin. 53 1039 (in Chinese) [盛阳 等2004 物理学报 53 1039]
21. Wang Y. et al 2005 Acta Phys. Sin. 54 2847 (in Chinese) [王音 等2005 物理学报 54 2847]
22. Y. Wang et al 2007 J. Mol. Struc.:THEOCHEM 807 201
23. Allen M. P. et al 1987 Computer Simulation of Liquids (Oxford: Clarendon Press) p20
24. J. Gao et al, accepted by J. Chem. Phys.
25.S.J. Austin et al 1995 J. Chem. Phys. 99 8076

- 193 -