• 2017-08-16 03:54:21 整理 王喜军 的几篇博客文章, 作为分子动力学模拟红外光谱的基础知识.

化学中经常用红外光谱来分析溶液的组成和变化，因为某些分子基团有红外特征指纹。问题是，溶剂和溶质的峰常常叠在一起，分析起来甚是棘手。所以，我们可以借助于分子动力学模拟来模拟溶剂的红外光谱，以便帮助分析整个溶液的红外光谱。

要想计算一种物质的红外光谱，最简单的方法是用量子化学计算气相中的一个单分子。即通常用高斯跑一个优化和频率(opt + freq) 计算。这种计算是在简谐近似（harmonic approximaiton，即弹簧近似）下求能量对原子位移的二阶导数得到频率，并按照吸收强度与对应跃迁偶极距的平方成正比的关系计算每个振动模式的相对强度。这种计算有两个弊端，一是简谐近似对于非简谐性强的体系，比如水分子簇、含有可旋转基团的柔性分子，误差较大；其次，气相单分子中没有分子间作用，而溶液中有分子间作用，也可引起较大误差。

怎么办呢？用凝聚态分子动力学模拟溶液体系，然后对体系极性演化的自相关函数做傅立叶分解，可得到既包含分子间作用效应，又不含简谐近似的红外光谱。那么它的原理是什么呢？这得从物质对红外光子的吸收谈起。

拿一束特定频率的红外光照射样品，那么红外光子有特定的能量 E=hνE=hν, <span class="MathJax" id="MathJax-Element-2-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">&#x03BD;" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;text-align:left;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">νν 是频率，如果此能量与样品中某分子振动态能级间的差值相等，分子就吸收之，跃迁到振动激发态。绝大多数激发都是基态向第一激发态跃迁，而第一激发能恰等于基态振动能，也就是说，基态频率=吸收频率。所以，MD模拟的话，只需模拟溶液基态，不需要考虑激发态。

除了频率，我们还需要知道吸收强度。吸收强度是由什么决定的呢？首先当然是浓度，分子数目越多，吸收越强。其次是跃迁偶极矩，即伴随振动引起的偶极矩改变。只有引起偶极矩变化的振动模式才可能吸收红外光子。为什么呢？光子由相互垂直的电场和磁场组成，当其中的电场方向与分子某一振动模式的跃迁偶极矩方向一致时，才可能形成有效吸收，因此，不改变分子偶极矩的振动是不会产生吸收的。此外，还与溶液体系的折射率有关，并且这个折射率会随频率变化。

因此，如果我们想探测溶液中的分子有哪些引起红外吸收的振动频率，需要寻找一个可观测或者可计算的性质，这个性质的微观量能与红外吸收相关，这个性质就是盒子(被模拟的理想体系)的总偶极矩。这里需要一点想象力：你的盒子里的分子、原子在不停的随机蠕动，每个蠕动都引起相关原子基团偶极矩的变化（想象一根C=O键的伸缩），所有基团蠕动形成的偶极矩变化构成总偶极变化，我们把总偶极矩在每个时刻的值记录下来，然后再把偶极矩信号转变为频率-强度信号，就能得到红外光谱。怎么转换？有个数学工具叫傅立叶变换，专门把时域信号转变为频域信号。

因此，总的流程就是，建立一个盒子，跑MD，同时记录总偶极-时间序列信号，做自相关函数，然后做傅立叶变换，即可得到吸收强度随频率变化的红外光谱。 从“总偶极-时间序列信号”得到红外光谱“频率-吸收强度”公式如下：

吸收强度 <span class="MathJax" id="MathJax-Element-3-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">I<mo stretchy="false">(&#x03C9;<mo stretchy="false">)" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;text-align:left;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">I(ω)I(ω) 随频率 <span class="MathJax" id="MathJax-Element-4-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">&#x03C9;" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;text-align:left;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">ωω 变化的函数:

后面积分过程就是傅立叶变换; <span class="MathJax" id="MathJax-Element-6-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML"><mi mathvariant="bold-italic">M" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;text-align:left;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">MM 是体系的总偶极矩, <span class="MathJax" id="MathJax-Element-7-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">&#x003C;<mi mathvariant="bold-italic">M<mo stretchy="false">(t<mo stretchy="false">)&#x22C5;<mi mathvariant="bold-italic">M<mo stretchy="false">(0<mo stretchy="false">)&#x003E;" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;text-align:left;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;"><M(t)⋅M(0)><M(t)⋅M(0)> 是体系总偶极矩的自相关函数。 <span class="MathJax" id="MathJax-Element-8-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">n<mo stretchy="false">(&#x03C9;<mo stretchy="false">)" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;text-align:left;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">n(ω)n(ω)是折射率。

量子效应会改变吸收强度，即影响峰的形状，但是不会改变振动频率。有很多经验公式来校正这个峰的形状。相关论文有很多讨论。需要记住，我们做的MD是经典MD，即牛顿力学决定的演化，不能包含核运动的量子效应，因此需要经验校正。

此外，我们也不知道折射率随频率变化的取值。因此，一般MD得到的红外光谱，其吸收强度其实是吸收强度与折射率的乘积。我们假设折射率为常数，这样我们可以比较峰的形状，或者不同峰的相对高度，但不能比较绝对强度。

那么，我们最关心的可能是，要想准确模拟红外光谱，什么是最重要的？即MD模拟红外光谱的主要误差来源是什么呢？是化学模型(Model Chemistry)，简单说就是用来计算能量和偶极的力场，或者第一原理方法。只有能量的二阶导数计算（频率）和偶极矩计算（强度）准确，最后总的峰形才能大致与实验对应。Material Studio中的经典力场能不能准确给出频率和偶极？这取决于参数化过程的好坏，做模拟前要验证一下（即做单分子的harmonic频率计算，比较力场和ab intio的结果）。即使参数很好，还有个问题，一般非极化力场不能描述分子间作用引起的极化，这些极化改变偶极矩，影响峰形。所以大部分MD模拟红外光谱的论文是用第一原理MD模拟，DFT计算天然包括了极化，偶极矩计算比较准确。其次，频率计算也远比力场要准确。

先用MD模拟体系，并对总偶极距的时间自相关函数做傅里叶变换，并经适当量子校正和平滑化，从而得到红外光谱的方法，可以称为MD-VDOS方法。这个名称比称为MD-IR更准确，因为实际的强度并非红外强度，而是红外强度和折射率的乘积。由于折射率随频率变化，且与体系的具体状态有关，所以无法确定。因此，可以退一步，称为振动态分布(vibrational density of states, VDOS)。这个方法的好处主要是能够包含非谐效应（模式耦合）、温度效应（展宽）、可适用于任意大小的体系。

上次CSTCC会议碰到一个做拟合势能面，然后用VSCF计算振动激发态的。他听了我的介绍后，认为我的方法不包含非谐效应。我稍微迟疑和谦虚了一下，他便一口咬定我的方法是错误的。这是他的误解，也说明像他这样研究红外多年的老手也不深刻理解什么是非谐效应。一开始我也曾犯过迷糊，后来有一次我跟同事一起深入讨论过这个问题，就弄明白了。所谓非谐效应，就是你用简谐近似方法（即Gaussian中的freq计算，又称为简正模式分析, normal mode analysis, NMA）得到的红外光谱与真实光谱之间的所有差异，都叫非谐效应。这不是一种真实的物理效应，而是对简谐近似这种方法造成的误差的统一描述。这就是为什么搞实验的人有时候也不明白“非谐效应”的原因，如果它们不做计算，根本谈不上简谐近似，非谐效应有什么意义？

简谐近似的根本优点在于，它能把体系振动简化成3N-6个基本振动模式(normal modes)，便于分析。它的根本缺陷在于，第一，它不能描述这些基本模式的相互影响或作用，这种影响称为耦合(modes coupling)。如果你的讨论不从3N-6个简正模式出发，你就谈不上任何简谐近似。这也就是为什么MD-VDOS方法得到的光谱与实验更接近，因为这个方法从头到尾这里面都没有涉及到任何简谐近似。简谐近似为什么能够简化问题？还是因为它把所在的势能面展开到了二阶；如果展开到三阶，那你需要考虑模式两两之间的作用；如果展开到四阶，就要考虑模式之间的三体作用。。。这跟HF/MP2/MP3完全是一个道理。所以我说，凡是用Gaussian做红外的，不管是展开到二阶还是多阶，他们都是基于简谐近似。简谐近似是他们难以摆脱的阴影。对于势能面平缓的体系，简直就是噩梦。

除了模式耦合，振动从基态到第二第三激发态造成的倍频、合频也是非谐效应的来源。VSCF能够很容易地对简谐振子进行基本量子化处理，然后讨论振动激发态，这是MD-IR不能比的。MD中原子核的运动包括振动、转动、平动是彻底经典的，不能给出任何形式的倍频、合频吸收。不过问题不大，因为常温下的中远红外光谱中，其中的倍频效应、合频效应这些涉及第二第三振动激发态的效应有多大呢？可以忽略。而且，理论计算作为一种分析方法，关注的主要是基本模式，用以帮助指认光谱和分析结构变化。就这个目的而言，问题不大。

除了非谐效应，两种方法还都受制于势能面的准确性。振动频率准不准，要看电子结构计算的绝热势能面。吸收强度计算则与偶极距的准确性有关。最后，计算流程的便利性也大大影响两种方法的受欢迎程度。这一方面，NMA稍胜一筹，像B3LYP+中等基组的势能面已经能准确计算小分子的振动频率和红外强度了，计算量也很小。与实验比，其误差与体系有关。比如，对非谐效应不显著的有机小分子而言，平均误差大概是系统性高估约5%。这样，稍微加个scaling factor，就能大大降低误差。而MD-IR目前普遍用平面波DFT或者半经验来跑MD，计算量比较大，计算过程稍微复杂，并未标准化。其中，MD采样有着很大技巧性，比如，初始构型的搭建，预平衡，总能守恒，偶极矩计算方法，量子校正方法的选择等等，都影响最后的结果。

很多第一性原理分子动力学软件可以计算固体的红外光谱。因为红外吸收是电磁波与周期性振动的偶极矩产生的电磁场相互作用产生的，因此，只要获得你这个固体体系的偶极矩随时间变化的曲线，然后做自相关函数的傅立叶变化，即可将时间-偶极信号，转变为频率-吸收强度信号。

要获得偶极矩随时间变化的曲线，可以使用CPMD, cp2k等基于DFT密度泛函的分子动力学模拟软件，生成一个大约20 ps-40 ps的动力学轨迹，每个运动步计算体系的偶极矩即可得到偶极矩曲线。

注意

1. 这里的体系偶极矩并非每个分子的偶极矩，而是整个模拟体系的偶极矩，它的周期性变化包含了各个基团的振动对总偶极矩的贡献；

2. 准确计算体系的偶极矩是关键，要选取合适的电子相关处理方法（广域梯度近似GGA比局域密度近似LDA要好，因为偶极计算需要准确的电子密度的梯度信息）和合适的交换相关泛函(PBE, BLYP, BP)，并测试所用的赝势。

得到体系偶极矩的演化轨迹以后，即可计算其自相关函数的傅立叶变化，将时域信号转换成频域信号，而频域信号的傅立叶变化在整个时间上的积分则可以得到在各频率上的吸收强度。具体理论可以参考文献 CPMD计算红外光谱, M.P. Gaigeot, JCTC, 1, 772-789, 2005

推荐使用CPMD，关于偶极轨迹的自相关-傅立叶变化的处理见：http://groups.google.ca/group/cpmd-list-cn/browse_thread/thread/a6a15953029189f9

自相关函数是用来衡量以某一频率出现某信号的规律化程度的。比如说有根键的振动频率是1500 cm^-1^波数，那么应该约每22 fs振动一个周期，那应该在总信号中存在这样一个周期强度变化，那你怎么去总信号中找呢？你需要以22 fs为间隔，测定轨迹上每个点与它之后22 fs那个点之间是否相关。如果存在这种22 fs的周期变化信号，则记相关系数为1，表示完全相关；无此相关，则记为0。不过体系在波动中，一根键的伸缩周期会随着分子碰来碰去而变化，所以你会得到一个居于0和1之间的数字来衡量这种相关程度。各种不同频率的信号交杂在一起，相关程度也各不相同。

式中 <span class="MathJax" id="MathJax-Element-10-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML"><mi mathvariant="bold-italic">M<mo stretchy="false">(0<mo stretchy="false">)" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;text-align:left;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">M(0)M(0) 表示22 fs间隔前的偶极矩信号强度, <span class="MathJax" id="MathJax-Element-11-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML"><mi mathvariant="bold-italic">M<mo stretchy="false">(t=22<mo stretchy="false">)" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;text-align:left;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">M(t=22)M(t=22) 表示22 fs之后的信号强度，两个尖括号表示对所有这种间隔的信号对相关程度取平均。当然，为了得到整个谱，你需要依次对每个频率的相关做计算，也就是说，对不同间隔t求其相关系数得到全谱。

自相关系数的计算公式见：维基百科: 自相关函数

通过体系总偶极的时间相关函数的傅里叶变换，可以得到红外光谱的谱密度I(ω)，并进而得到吸收系数α(ω)。其中ω是频率。问题是，这个时间相关函数是通过经典分子动力学模拟得到的。而推导谱密度和吸收系数的公式，则是基于量子跃迁速率公式得到。在量子态跃迁速率公式中，使用的是量子的时间相关函数。二者有显著的差异，并直接导致未经校正的经典时间相关函数得到谱形(shape function)与真实的谱存在明显的差异。

那么，经典的和量子的时间相关函数二者究竟有什么不同呢？首先，量子的时间相关函数是一个复函数(complex function)，而经典的，则是一个实函数。（这意味着什么？）其次，量子时间相关函数存在着细致平衡假定。所谓细致平衡，即量子体系的态分布处于平衡态，或者说，在任意时刻内，从某态跃迁而离开的粒子的数目，等于从其它态跃迁至该态的粒子数目，从而处于各种状态的粒子数目是固定不变的。这是一个热平衡态的量子描述。在经典热平衡理论中，没有态这个概念。所有量的变化都是连续的，即体系处在连续的变化中。经典平衡态是指体系的各种性质在平均值附近波动。这两者显然有着很大的不同。

因此，需要引入一些假定，进行数学推导，得到对经典时间相关函数的量子校正公式，来校正得到的谱密度，或者说shape function。这个容易理解。常见的量子校正因子包括：(1) Oxtoby, Berne, Litoviz等人给出的标准修正；(2) Wilson, Frommhold, Berne等人给出的简谐修正；(3) Schofield修正；(4) Skinner的简谐+Scholdfieild修正 (5) Egelstaff修正。

那么，为什么存在多种不同的量子校正呢？为什么这个校正不是唯一的？文献中，对这个问题的答案是，经典极限不是唯一的。

The well-known problem is that there is no unique classical limit available in view of the fact that the correspondence of the different quantum correlation functions to a single classical correlation function leads to different quantum correction schemes.

什么是经典极限classical limit？各个修正之间的根本不同在哪里？

经典极限classical limit是量子力学中的一个概念，是说从量子力学的理论公式出发，推导出经典力学的理论公式。怎么推导呢？令量子力学公式中的普朗克常数ℏ趋于0。普朗克常数用于描述量子的大小，比如一个振动频率为ω的振子的量子能级之差为ℏω。当ℏ=0时，量子态之间的差为0, 即成为连续。这是理论上的理解。实际上量子体系的时间相关函数是很难得到的，基本上还没人得到。怎么办呢？一般是把一些简化的理想体系，比如一个简谐振子，进行量子化处理，得到近似的量子时间相关函数，然后对比其经典的时间相关函数，可以得到其间的关系，即可得到校正因子。这样，通过对不同的理想模型，以及其间作用的不同假定，可以得到不同的校正因子。比如，简谐振子+线性耦合可以得到Wilson简谐修正因子，而简谐振子+指数耦合则得到Skinner校正因子。

在CPMD附带的由Axel等编写的fourier变换程序中，第一列是<frequencies>，后面2到7列<value1>~<value6>是各种结果。其中

<value1>和<value6>没有意义；

<value2>(即第3列)是不含任何修正

<value3>(即第4列)是标准修正，又称Kubo修正；

<value4>(即第5列)是简谐修正harmonic；

<value5>(即第6列)是Scholfield修正；

参考

1. 田国才，邵久书. 凝聚相中振动能量弛豫的理论研究

2. Rafael Ramırez et al, Quantum correction to classical time-correlation functions

3. J. L. Skinner, Quantum dynamics in simle fluid

◆本文地址: https://jerkwin.github.io/2017/08/16/整理-使用分子动力学模拟红外光谱/, 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆