光谱之家

如题：怎么采用MD方法 算某一体系的红外光谱 。
假如，我模拟的某一个溶剂体系。采用什么样的方法求的 这种溶剂分子的光谱？
计算的原理是什么？
我接触分子模拟不久，希望高手给指导下。不胜感激！！！
若金币不够，我可以再加！！！！

拿一束特定频率的红外光谱照射样品，那么红外光子有特定的能量E=hv, v是频率，如果此能量与样品中某分子振动态能级间的差相等，分子就吸收之，飞到振动激发态。绝大多数激发都是激发态向第一激发态跃迁。而第一激发能恰等于基态振动能，也就是吸收光子频率v等于分子基态振动频率v。所以，MD模拟的话，只需模拟基态即可。你的盒子里有哪些振动频率，实际就会对应哪些频率的吸收。
然而，除了频率，我们还需要知道吸收强度。吸收强度是什么决定的呢？首先当然是浓度，分子数目越多，吸收越强。其次是方向。光子由相互垂直的电场和磁场组成，当其中的电场方向与分子振动模式的跃迁偶极矩方向一致时，才有最大可能吸收。此外，还与溶液体系的折射率有关，而这个折射率会随频率变化。
综上分析，我们需要寻找体系的某个性质，这个性质的微观量的变化频率由振动频率变化决定，而且强度与跃迁偶极距成正比。这个性质就是体系的总偶极距。这里需要一个很好的想象力：你的盒子里的分子、原子在不停的随机蠕动，每个蠕动都引起相关原子基团偶极距的变化（想象一根C=O键的伸缩），所有基团蠕动形成的偶极距变化构成总偶极变化。其中，分子基团蠕动形成的偶极变化与蠕动本身的偶极矩有关。因此，如果我们把总偶极矩在每个时刻的值记录下来，我们就能把其中每个蠕动的贡献分解出来。怎么分解，有个数学工具叫傅立叶分析，专门把混合信号分解成单个的信号。
吸收强度I(w)随频率w变化的函数： I(w)*n(w) = 量子校正系数 * 积分{ exp(-iwt) *   }  
后面积分{}就是傅立叶变换；M是体系总偶极距， 是体系总偶极距的自相关函数。n(w)是折射率。
因此，我们要做的，就是建立一个盒子，跑MD，同时记录盒子总偶极变化，最后做傅立叶分解，得到红外光谱。
量子效应会改变吸收强度，即影响峰的形状，但是不会改变振动频率。有很多经验公式来校正这个峰的形状。相关论文有很多讨论。需要记住，我们做的MD是经典MD，即牛顿力学决定的演化，不能包含量子效应，因此需要经验校正。
此外，我们也不知道折射率随频率变化的取值。因此，一般MD得到的红外光谱，其吸收强度其实是吸收强度与折射率的乘积。我们假设折射率为常数，这样我们可以比较峰的形状，或者不同峰的相对高度，但不能比较绝对强度。
那么，我们最关心的可能是，要想准确模拟红外光谱，什么是最重要的？即MD模拟红外光谱的主要误差来源是什么呢？是化学模型(Model Chemistry)，简单说就是用来计算能量和偶极的力场，或者第一原理。只有能量的二级导计算（频率）和偶极矩计算（强度）准确，最后总的峰形才能大致与实验对应。Material Studio中的经典力场能不能准确给出频率和偶极？这取决于参数化过程的好坏，做模拟前要验证一下（即做单分子的Normal mode analysis的harmonic频率计算，比较力场和ab intio的结果）。即使参数很好，还有个问题，一般非极化力场不能描述分子间作用引起的极化，这些极化改变偶极距，影响峰形。大部分MD模拟红外光谱要用第一原理MD模拟。原因是DFT计算天然包括了极化，偶极距计算比较准确。其次，频率计算也远比力场要准确。
有了上述基本概念，希望对你进一步学习文献有帮助。

您好，大侠。看了您的帖子，我茅塞顿开，但是我有一点不大明白：
M是体系总偶极距， 是体系总偶极距的自相关函数。这个 什么意思，角括号代表什么含义？ 不胜感激！！

自相关函数是用来衡量以某一频率出现某信号的规律化程度的。比如说有根键的振动频率是1500cm-1 波数，那么应该约每50fs（没做计算，假设吧）振动一个周期，那应该在总信号中存在这样一个周期强度变化，那你怎么去总信号中找呢？你需要以50fs为间隔，测定轨迹上每个点与它之后50fs那个点之间是否相等对吧。如果存在这种以50fs为周期的信号相关，则记起相关系数为1，表完全相关；无此相关，则记为0。不过体系在波动中，一根键的伸缩周期会随着分子碰来碰去而变化，所以你会得到一个居于0和1之间的数字来衡量这种相关程度。
式中M(0)表50fs间隔前的偶极距信号强度，M(t=50)表50fs之后的信号强度，两个尖括号<>表示对所有这种间隔的信号对相关程度取平均。当然，为了得到整个谱，你需要依次对每个频率的相关做计算，也就是说，对不同间隔t求其相关系数得到全谱。
一列信号的自相关系数的计算公式见cuturl('http://zh.wikipedia.org/wiki/')自相关函数

收益匪浅。我有以下两个表述，不知道对不对，请教。
（1） BOMD和CPMD都可以得到偶极随时间的变化？
CPMD软件可以算CP动力学，也可以算BO动力学，结果应该是相近的吧？
（2）光是电磁波，它与物质振动之间能量传递有两个条件：振动频率一致，振动结构有偶极变化。如果这种振动的频率落在红外波数，就是红外光谱。不是所有声子都有偶极变化，所以不是所有振动都有光学响应。

另外：
如果通过CPMD软件的BO-/CP-动力学都能得到红外光谱，那是否可以通过VASP软件实现的BO-动力学得到红外光谱呢？
VASP-BOMD与CPMD-BOMD在计算结果上应该是等效的？

（1） Agree。BO与CP只是计算波函数（从而计算原子受力）的方式不同，但是可以认为两种方法得到的动力学结果是等价的。
（2） Agree。红外选律强调红外吸收的前提是存在偶极变化的振动。我们通过MD和总偶极的分解去模拟红外光谱，自然是符合这个选律的。

当然可以。cp2k, quantum-espresso都可以做MD红外光谱。只要这个第一原理软件支持MD和计算每一步的总偶极距即可。剩下的就可以用Fourier transform的小脚本处理一下即可得到红外光谱。
需要指出的，MD红外光谱的波形（位置，峰高和峰宽）与键之间的强度有关，与所计算的偶极距的准确程度相关。对于平面波凝聚态盒子体系，直到berry phase方法（modern polarization theory）出现以后才解决了这个问题。文献中给出的波形都比较合理。现在我在做的一项工作是对其准确程度给一个定性或定量的描述。

那么，VASP还有一种算Gamma点谐振近似得到频率的方法，比如算OH在金属氧化物表面吸附结构的谐振频率，这个数据对应红外谱上若干的峰位，但没有峰高、峰宽，也不能连成连续的光谱。
请问这种计算是否能通过某种方法变成连续的红外光谱？

你说的是正则分析得到的谐振频率。这是最经常用的一种计算红外光谱的方法，特别是用Gaussian计算分子和分子簇的红外光谱。VASP谐振频率计算原理和gaussian中完全一样。都是通过hessian矩阵变换得到频率。其峰高，即吸收强度，与某正则振动模式qi导致的偶极µ变化的平方成正比，
Ai = NA*pi/(3*c^2) * (∂µ/∂qi)^2
ref: JPC 92, 1739
大多数平面波code都不计算吸收强度，原因不清楚。可能是计算不准确，也可能是因为developer并不关心这个强度。不管怎样，需要自己按照上述公式去计算这个强度即可。
得到频率和强度，你就可以得到线光谱，没有峰宽，这是因为简谐频率计算计算的是0点下的振动，所以不存在broadening,即展宽。怎么展宽呢？下面是我的一篇blog：《高斯展宽，洛伦慈展宽，Voigt展宽的fortran code》(cuturl('http://hi.baidu.com/xijunw/item/2b581b1ad4f463fd746a841c'))
"展宽，是对量子化学计算出来的分子的线性光谱（红外，拉曼等）展开成带状峰。这是一种经验处理，只是为了视觉上更接近于实验得到的光谱。既然实验得到的是带状光谱，那为什么量子化学计算得到是线性光谱呢？原因很简单，量子化学计算的气相分子在0K下的振动光谱，是在势能面的最低点根据二次简谐近似得到的振动频率。这个频率-吸收强度数据对，体现到谱图上就是一根线。实际的光谱总是在一定的温度下的凝聚相（固体或液体）下的振动光谱，它是极大量分子（～Na, 10^23）的平均结果，每个分子都具有大小不同的能量，处于不断运动、碰撞中，分布在各种激发态能级中。这些都会或多或少地改变其指定模式的振动频率，使得它们或大或小，按照一定的概率分布在特征频率的两侧，形成峰带。展宽是多种效应的总和效果。"
"按照量子电动力学，这些展宽效应分为两类：非均展宽，和均匀展宽。这个均和非均，并非是指展宽的峰形是否对称，而是指效应本身是否对所有的分子造成的后果是否一样。比如，杂质效应(dopant effect)，是指分子受到其中混杂的少量异种分子的作用后能级的变化。这种效应必然随着与杂质分子的距离而变化，因而属于非均展宽。然而，如果是分子在体系压力下造成的能级变化，那么对于所有的分子来说，振动频率受到影响的机会是均等的，因而属于均匀展宽。另外一个均匀展宽的例子是多普勒效应，即高速向你运动的物体，测得的频率会升高，而远离你的物体，频率会降低。非均展宽的概率分布符合gaussian分布，均匀展宽的概率分布则是Lorentzian型的。如果多种不同效应之间互相独立，则综合效果符合Voigt峰形。实际上各效应之间的耦合是不可避免的。但是Voigt展宽仍然是最符合实际的。"