乱弹辛积分动力学，辛流形，对称，与守恒律
2011-01-09 09:23
|
看到一篇用SPC/E water model计算水在-4‘C和80’C下光谱红移的文章[1]，(J. Phys. Chem. A 2004, 108, 11056-11062)，其中为了获得能量和角动量守恒，他们的动力学积分采用了一种很特别的算法：simplectic integration (split integration symplectic method (SISM))，中文干活地话叫辛积分，纠结到不行。于是花了一天时间看到底什么是辛，什么是辛积分。虽然没有完全懂，但是对动力学产生了一些更多的理解。
首先，辛这个词产生于数学上的"辛流形"。那什么是“流形”？这个词如雷贯耳，但是总是不甚了了。令人敬仰的丘成桐就是搞流形的。据说这个东西能化腐朽为神奇。看wiki半天，得出一个结论，流形就是图集，由n维图构成的一个n+1维的图集，就是流形。比如，把一个地球仪划分成很多小块块，每个小块块可以看成一个个简单的二维平面（图），然后把它们拼在一起可以构成一个三维的图集，就构成一个流形。为毛要这样搞呢? 因为计算可以在低维度进行，然后把结果传回流形，因为它们是“拓扑同胚”的（化学生到此止步）。 那啥是辛流形呢？Wiki: “数学上，一个辛流形是一个装备了一个闭、非退化2-形式ω的光滑流形，ω称为辛形式。” ( 辛形式看的人很辛苦啊) 想象上面说的地球仪那个流形(图集)，辛流形是说这个拼起来的图啊，不能缺一块，或者有棱角，就叫封闭光滑(可微可积)。啥是非退化，啥是2形式呢？k形式是在一个光滑的流形的子集（图）元素之间定义计算规则，定义了正负叫1-形式，定义点积叫2-形式，定义乘积叫3-形式；非退化，是指矩阵是满秩的。满秩是可逆的。所以综合来讲，非退化就是有确定维数，能进行互逆操作的空间（满足对称）。 这样辛苦的折腾出来一个辛流形有啥用呢？描述相空间。这下搞过分子动力学的都明白，相空间一堆粒子N个，每个粒子每个时刻都有确定的坐标(3N)和动量(3N)，那么整个体系由所有粒子的坐标变量和动量变量构成的6N为空间称为相空间。体系在一些特定的约束下运动（不同的约束构成不同的系综），由于不同种粒子之间的作用不同，所以体系在相空间每个点的能量也不同，如果计算出来在相空间每个点的能量大小（遍历相空间），那么根据统计热力学的配分函数，就可以计算出来体系的所有性质，简单说就是某个状态能量高，则它对体系平均性质贡献就小，如果能量低，则它对体系平均性质贡献就大，成负指数关系，所以平均性质按能量的负指数函数做加权平均。(如果只考虑N粒子坐标变量，构成的3N维空间叫位形空间，又叫组态空间) 辛流形跟相空间有啥关系呢？wiki: "一个系统的所有组态的空间(位形空间)可以用一个流形建模，而该流形的余切丛描述了该系统的相空间。"原来位形空间是一个流形，是一个图集。啥是余切丛呢？wiki: " 微分几何中，流形的余切丛是流形每点的切空间组成的向量丛。可以在余切丛上定义一组特殊的坐标系；这些被称为正则坐标。因为余切丛可以视为辛流形，任何余切丛上的实函数总是可以解释为一个哈密顿函数；这样余切丛可以理解为哈密顿力学讨论的相空间。" ( 相空间是6N维，这里又说只是3N维的余切丛，wiki不是坑爹吗) 姑且这样理解，整个体系在3N维位形空间流动，每个点都有一条切线（光滑的嘛），这个切线的斜率就是这点的一阶导（受力），但是这个一阶导实际上是由所有的粒子的一阶导（受力）线性加和构成的，所有粒子的一阶导(3N个分量fix, fiy, fiz)构成了向量丛。每个粒子都有坐标向量和动量向量，构成正则坐标，正则坐标构成相空间，哈密顿为此相空间的实函数。而这个相空间就是一个辛流形，是位形空间上的满足封闭，光滑，可点乘，不退化，满足空间对称和时间对称的哈密顿函数。所以辛流形描述了相空间与哈密顿之间的数学特性。 来看看时间反演对称和空间反演对称。这个容易理解，因为在经典动力学范畴里，往天上抛一个球，做成录像，正着放和倒着放一样，因为每个时刻遵循相同的力学原理。不过，最奇异的是Noether(内特)原理[3]，德国数学家Emmy Noether说：“：一个力学系统动力学行为的每一个对称性都意味着该系统的一个守恒律”。惯性封闭质点系中，时间均匀性(或时间平移对称性)意味着能量守恒。空间平移对称意味着动量守恒，空间转动对称意味着角动量守恒。“倘若我们再依据因果律，把时空均匀性和各向同性，即时空平移对称性和转动对称性，看作原因的对称性，而系统的能量、动量和角动量守恒律看作结果的对称性，则可引出结论：原因中的对称性必反映在结果中，这就是对称性原理，首先由P.居里于1894年提出.”。神奇。这就是个别文献上提到的一个词causally (因果地) 的含义。 再回到本文开头的辛积分动力学。[2]的介绍很好，转录如下：“长期以来，对经典轨迹法(QCT)的结果的改进主要来自于拟台的更精确的势能面而经典轨迹法的计算结果本身也被当作势能面拟合得好不好的一个判断方法。但是，在我们的研究中我们发现，对于很多体系来说，采用耗散的积分方法(Runge-Kutta是其中的一个代表，自然还包括其他一些方法)，当轨迹时间较长时，发生能量损失现象。换句话说，在长时间的模拟时，体系表现出能量耗散现象。这一现象在研究体系非线性行为，体系IVR过程，以及大分子动力学过程等需要长时问模拟计算的系统时，将会变得十分严重，值得关注。从数学上来看，准经典轨迹计算的质量，显然不仅与势能面的选择有关，而且与所用的积分算法有关。也就是说，并不是采用的势能面更精确，经典轨迹法计算的结果就会更好，积分方法在其中也起到了一定的作用。所以说，改进经典轨迹法的计算结果包括两个方面：势能面的改进，以及积分方法的改进。近年来，以辛几何为理论框架发展了一种新的 Hamilton系统数值积分方法，其插分格式可以精确地反映体系Hamilton相流最重要的几何特性，即辛结构的保持，而辛结构是Hamilton相流的一种整体结构。这种结构能否保持，往往决定了长期演化计算的有效与否。具体到准经典轨迹计算所研究的Hamilton体系来说，即决定了体系长时间演化时能量和角动量是否守恒。在辛算法出现以前，Hamilton系统的一些长时间演化的研究所采用的各种数值方法都因为其截断误差的累积，导致人工耗散而不能保持其辛结构，从而引进一些非Hamilton相流的特征，甚至完全歪曲了运动性质，导致QCT研究中体系长时间演化时能量和角动量不守恒．因此，在动力学研究中，辛算法的应用有着很大的价值。“ 他们计算H+H2-->H2+H的化学反应动力学表明，采用四级Runge—Kutta算法，60000(fs)步以后能量耗散(降低)约为0.001eV (~0.1kJ/mol)，而采用4级或6级辛算法，则观察不到能量耗散。看上去4阶Runge-Kutta算法产生的能量耗散0.1kJ/mol/60ps (每ns则是1.7kJ/mol) 也不严重啊，对体系总体性质的影响作者说不好估计，因为结果产物分布还受到初始构型的影响。 回到文献[1]，他们打算模拟液态水振动光谱的bending band随温度升高而变窄的实验现象(J. Mol. Struct. 1994, 322, 105)，为了能用经典模拟捕捉这个与直觉相悖的却与核运动的量子效应无关(?)的现象，他们采用了SPC和SPC/E两种water model。其中为了避免能量耗散，他们采用了分裂积分辛方法(split integration symplectic method, SISM). 首先哈密顿分裂成2项H=H0+Hr, H0是可解析的部分，就是力场中为坐标q的函数的能量项比如bond, angle, cross terms, vdW(LJ12-6), coulomb，Hr是剩余的部分，即力场中并非坐标显函数的部分，包括动能，静电作用的长程校正(反应场模型)，分别对这2个部分进行解析积分和数值积分。然后采用经典传播子(classical propagator, exp(Δt*LH)) 的二级近似即广义蛙跳算法(generalized leapfrog scheme)进行传播(含时演进)[4-9]。由于作者这一部分一带而过，还要仔细看他的引用文献，看看公式是怎样的，如何传播，跟常用velocity verlet有何不同。
他们的结论很有趣，因为折腾了半天，也没有模拟出来那个传说中液态水振动光谱bending band随温度升高变窄的灵异现象，倒是把其它一堆性质预测的很好，比如stretching 蓝移，分子内bending红移，以及介电常数等。作者分析说，这个现象跟量子效应无关，所以跑quantum dynamics也未必有用；说“高温下水结构更有序”的说法也可以排除；作者不太确定地说，大概势能面有问题，说不定用个比SPC/E更好的model能把bending带变窄模拟出来，比如可极化力场，并且参数化要结合反应场模型来搞(指长程校正)。或者用CPMD也许能行。 1. Matej Praprotnik, Dusˇanka Janezˇicˇ, and Janez Mavri*, Temperature Dependence of Water Vibrational Spectrum: A Molecular Dynamics Simulation Study. J. Phys. Chem. A 2004, 108, 11056-11062 2. 吴韬等，辛算法在准经典轨迹计算中的应用，中国科学B, V32, P46. 2002 3. http://kejian.tzc.edu.cn/ckwx/pdf/268.pdf 4. Strang, G. SIAM J. Numer. Anal. 1968, 5, 506. 5. Grubmu¨ller, H.; Heller, H.; Windemuth, A.; Schulten, K. Mol.Simul. 1991, 6, 121. 6. Wilson, E. B.; Decius, J. C.; Cross, P. C. Molecular Vibrations; McGraw-Hill Book Company, Inc.: New York, 1955. 7. Brooks, B. R.; Janezˇicˇ, D.; Karplus, M. J. Comput. Chem. 1995,16, 1522. 8. Janezˇicˇ, D.; Brooks, B. R. J. Comput. Chem. 1995, 16, 1543. 9. Janezˇicˇ, D.; Venable, R. M.; Brooks, B. R. J. Comput. Chem. 1995,16, 1554. |
