博文
如何用分子动力学来模拟液体的红外光谱？
|||
化学中经常用红外光谱来分析溶液的组成和变化，因为某些分子基团有红外特征指纹。问题是，溶剂和溶质的峰常常叠在一起，分析起来甚是棘手。所以，我们可以借助于分子动力学模拟来模拟溶剂的红外光谱，以便帮助分析整个溶液的红外光谱。
要想得到一种物质的红外光谱，最简单的方法是用量子化学分析气相中的一个单分子。即通常用高斯跑一个优化和频率(opt + freq) 计算。这种计算是在简谐近似（harmonic approximaiton，即弹簧近似）下求能量对原子位移的二阶导得到频率，并按照吸收强度与对应跃迁偶极距的平方成正比的关系计算每个模式的相对强度。这种计算有两个弊端，一是简谐近似对于非谐性强的体系，比如水分子簇、含有可旋转基团的柔性分子，误差较大；其次，气相单分子中没有分子间作用，而溶液中有分子间作用，也可引起较大误差。
怎么办呢？用凝聚态分子动力学模拟溶液体系，然后对体系极性演化的自相关函数做傅立叶分解，可得到既包含分子间作用效应，又不含简谐近似的红外光谱。那么它的原理是什么呢？这得从物质对红外光子的吸收谈起。
拿一束特定频率的红外光谱照射样品，那么红外光子有特定的能量E=hv, v是频率，如果此能量与样品中某分子振动态能级间的差相等，分子就吸收之，飞到振动激发态。绝大多数激发都是基态向第一激发态跃迁，而第一激发能恰等于基态振动能，也就是说，基态频率=吸收频率。所以，MD模拟的话，只需模拟溶液基态，不需要考虑激发态。
除了频率，我们还需要知道吸收强度。吸收强度是什么决定的呢？首先当然是浓度，分子数目越多，吸收越强。其次是跃迁偶极矩，即伴随振动引起的偶极矩改变。只有引起偶极矩变化的振动模式才可能吸收红外光子。为什么呢？光子由相互垂直的电场和磁场组成，当其中的电场方向与分子某一振动模式的跃迁偶极矩方向一致时，才可能形成有效吸收，因此，不改变分子偶极矩的振动是不会产生吸收的。此外，还与溶液体系的折射率有关，并且这个折射率会随频率变化。
因此，如果我们想探测溶液中的分子有哪些引起红外吸收的振动频率，需要寻找一个可观测或者可计算的性质，这个性质的微观量能与红外吸收相关，这个性质就是盒子(被模拟的理想体系)的总偶极距。这里需要一点想象力：你的盒子里的分子、原子在不停的随机蠕动，每个蠕动都引起相关原子基团偶极距的变化（想象一根C=O键的伸缩），所有基团蠕动形成的偶极距变化构成总偶极变化，我们把总偶极矩在每个时刻的值记录下来，然后再把偶极矩信号转变为频率-强度信号，就能得到红外光谱。怎么转换？有个数学工具叫傅立叶变换，专门把时域信号转变为频域信号。
吸收强度I(w)随频率w变化的函数： I(w)*n(w) = 量子校正系数 * 积分{ exp(-iwt) * <M(t)*M(0)> }
后面积分{ }过程就是傅立叶变换；M是体系总偶极距， <M(t)*M(0)>是体系总偶极距的自相关函数。n(w)是折射率。
量子效应会改变吸收强度，即影响峰的形状，但是不会改变振动频率。有很多经验公式来校正这个峰的形状。相关论文有很多讨论。需要记住，我们做的MD是经典MD，即牛顿力学决定的演化，不能包含量子效应，因此需要经验校正。
此外，我们也不知道折射率随频率变化的取值。因此，一般MD得到的红外光谱，其吸收强度其实是吸收强度与折射率的乘积。我们假设折射率为常数，这样我们可以比较峰的形状，或者不同峰的相对高度，但不能比较绝对强度。
那 么，我们最关心的可能是，要想准确模拟红外光谱，什么是最重要的？即MD模拟红外光谱的主要误差来源是什么呢？是化学模型(Model Chemistry)，简单说就是用来计算能量和偶极的力场，或者第一原理方法。只有能量的二级导计算（频率）和偶极矩计算（强度）准确，最后 总的峰形才能大致与实验对应。Material Studio中的经典力场能不能准确给出频率和偶极？这取决于参数化过程的好坏，做模拟前要验证一下（即 做单分子的harmonic频率计算，比较力场和ab intio的结果）。即使参数很好，还有个问题，一般非极化力场不能描述分子间作用引起的极化，这些极化改变偶极距，影响峰形。所以大部分MD模拟红外光谱的论文是用第一原理MD模拟，DFT计算天然包括了极化，偶极距计算比较准确。其次，频率计算也远比力场要准确。
【补充】自相关函数是怎么回事儿？
【原发小木虫，经修改。】
http://blog.sciencenet.cn/blog-38740-597173.html
要想得到一种物质的红外光谱，最简单的方法是用量子化学分析气相中的一个单分子。即通常用高斯跑一个优化和频率(opt + freq) 计算。这种计算是在简谐近似（harmonic approximaiton，即弹簧近似）下求能量对原子位移的二阶导得到频率，并按照吸收强度与对应跃迁偶极距的平方成正比的关系计算每个模式的相对强度。这种计算有两个弊端，一是简谐近似对于非谐性强的体系，比如水分子簇、含有可旋转基团的柔性分子，误差较大；其次，气相单分子中没有分子间作用，而溶液中有分子间作用，也可引起较大误差。
怎么办呢？用凝聚态分子动力学模拟溶液体系，然后对体系极性演化的自相关函数做傅立叶分解，可得到既包含分子间作用效应，又不含简谐近似的红外光谱。那么它的原理是什么呢？这得从物质对红外光子的吸收谈起。
拿一束特定频率的红外光谱照射样品，那么红外光子有特定的能量E=hv, v是频率，如果此能量与样品中某分子振动态能级间的差相等，分子就吸收之，飞到振动激发态。绝大多数激发都是基态向第一激发态跃迁，而第一激发能恰等于基态振动能，也就是说，基态频率=吸收频率。所以，MD模拟的话，只需模拟溶液基态，不需要考虑激发态。
除了频率，我们还需要知道吸收强度。吸收强度是什么决定的呢？首先当然是浓度，分子数目越多，吸收越强。其次是跃迁偶极矩，即伴随振动引起的偶极矩改变。只有引起偶极矩变化的振动模式才可能吸收红外光子。为什么呢？光子由相互垂直的电场和磁场组成，当其中的电场方向与分子某一振动模式的跃迁偶极矩方向一致时，才可能形成有效吸收，因此，不改变分子偶极矩的振动是不会产生吸收的。此外，还与溶液体系的折射率有关，并且这个折射率会随频率变化。
因此，如果我们想探测溶液中的分子有哪些引起红外吸收的振动频率，需要寻找一个可观测或者可计算的性质，这个性质的微观量能与红外吸收相关，这个性质就是盒子(被模拟的理想体系)的总偶极距。这里需要一点想象力：你的盒子里的分子、原子在不停的随机蠕动，每个蠕动都引起相关原子基团偶极距的变化（想象一根C=O键的伸缩），所有基团蠕动形成的偶极距变化构成总偶极变化，我们把总偶极矩在每个时刻的值记录下来，然后再把偶极矩信号转变为频率-强度信号，就能得到红外光谱。怎么转换？有个数学工具叫傅立叶变换，专门把时域信号转变为频域信号。
因此，总的流程就是，建立一个盒子，跑MD，同时记录总偶极-时间序列信号 ，做自相关函数，然后做傅立叶变换，即可得到吸收强度随频率变化的红外光谱。 从“总偶极-时间序列信号”得到红外光谱“频率-吸收强度”公式如下：
吸收强度I(w)随频率w变化的函数： I(w)*n(w) = 量子校正系数 * 积分{ exp(-iwt) * <M(t)*M(0)> }
后面积分{ }过程就是傅立叶变换；M是体系总偶极距， <M(t)*M(0)>是体系总偶极距的自相关函数。n(w)是折射率。
量子效应会改变吸收强度，即影响峰的形状，但是不会改变振动频率。有很多经验公式来校正这个峰的形状。相关论文有很多讨论。需要记住，我们做的MD是经典MD，即牛顿力学决定的演化，不能包含量子效应，因此需要经验校正。
此外，我们也不知道折射率随频率变化的取值。因此，一般MD得到的红外光谱，其吸收强度其实是吸收强度与折射率的乘积。我们假设折射率为常数，这样我们可以比较峰的形状，或者不同峰的相对高度，但不能比较绝对强度。
那 么，我们最关心的可能是，要想准确模拟红外光谱，什么是最重要的？即MD模拟红外光谱的主要误差来源是什么呢？是化学模型(Model Chemistry)，简单说就是用来计算能量和偶极的力场，或者第一原理方法。只有能量的二级导计算（频率）和偶极矩计算（强度）准确，最后 总的峰形才能大致与实验对应。Material Studio中的经典力场能不能准确给出频率和偶极？这取决于参数化过程的好坏，做模拟前要验证一下（即 做单分子的harmonic频率计算，比较力场和ab intio的结果）。即使参数很好，还有个问题，一般非极化力场不能描述分子间作用引起的极化，这些极化改变偶极距，影响峰形。所以大部分MD模拟红外光谱的论文是用第一原理MD模拟，DFT计算天然包括了极化，偶极距计算比较准确。其次，频率计算也远比力场要准确。
【补充】自相关函数是怎么回事儿？
自相关函数是用来衡量以某一频率出现某信号的规律化程度的。比如说有根键的振动频率是1500cm-1 波数，那么应该约每22fs振动一个周期，那应该在总信号中存在这样一个周期强度变化，那你怎么去总信号中找呢？你需要以22fs为间隔，测定轨迹上每个点与它之后22fs那个点之间是否相关。如果存在这种50fs的周期变化信号，则记起相关系数为1，表完全相关；无此相关，则记为0。不过体系在波动中，一根键的伸缩周期会随着分子碰来碰去而变化，所以你会得到一个居于0和1之间的数字来衡量这种相关程度。各种不同频率的信号交杂在一起，相关程度也各不相同。
式中M(0)表50fs间隔前的偶极距信号强度，M(t=50)表50fs之后的信号强度，两个尖括号<>表示对所有这种间隔的信号对相关程度取平均。当然，为了得到整个谱，你需要依次对每个频率的相关做计算，也就是说，对不同间隔t求其相关系数得到全谱。
自相关系数的计算公式见： http://zh.wikipedia.org/wiki/自相关函数
【原发小木虫，经修改。】
http://blog.sciencenet.cn/blog-38740-597173.html
No comments:
Post a Comment