使用Multiwfn绘制红外、拉曼、UV-Vis、ECD和VCD光谱图

使用Multiwfn绘制红外、拉曼、UV-Vis、ECD和VCD光谱图

文/Sobereva(3)
First release: 2014-Mar-19  Last update: 2016-Feb-5

1 前言

Multiwfn是功能最为强大的波函数分析程序,虽然光谱图的绘制并不是Multiwfn分内的事,但是由于光谱的绘制在计算化学中经常要涉及,所以早在2.0版的时候就加入了光谱绘制的功能,在后续版本中又逐步改进。Multiwfn的这个功能不仅比起GaussView绘制光谱的功能灵活得多得多,可以满足高端用户的需求,比起专业的产生光谱数据的Swizard也完全是有过之而无不及(Swizard根本都没法直接绘制图形来),更是完爆GaussSum。在此文将介绍在Multiwfn中绘制光谱的方法。由于很多人并不清楚理论计算产生光谱的基本原理,比如经常见到有人搞不懂振子强度和峰高、吸光度的关系,因此在第二节将首先介绍一些基本原理。如果急着马上就作出图来,可以直接看3.1节对输入文件的说明和第5节的例子。

本文所使用的Multiwfn为3.3.8(dev)版,可以在Multiwfn主页http://multiwfn.codeplex.com上免费下载。老版本Multiwfn不支持本文提到的一些特征,故务必使用最新版。不排除在未来的版本中在用法上出现细微变动的可能,届时以手册为准。另外,如果是第一次接触Multiwfn程序,建议先阅读《Multiwfn入门tips》(http://sobereva.com/167)了解一些基本信息。


2 理论光谱产生的原理

2.1 IR、VCD、ECD、UV/Vis光谱

物质的实验光谱是一条连续的吸收曲线,包含一大堆吸收峰,表现了对不同频率的光的吸收度的不同,一般用摩尔吸收系数(epsilon,ε)来衡量,含义是溶液浓度为1mol/L、液层厚度为1cm时的吸光度,单位是L/mol/cm。

理论光谱计算给出的是离散的跃迁数据。比如做电子激发计算,程序会给出基态到各个电子激发态的跃迁能以及振子强度。比如下图,黑色的一条条竖线的横坐标位置就是电子激发能,竖线高度就是振子强度(对应右侧坐标轴)。



显然理论计算给出的离散的跃迁数据和实验给出的连续的吸收曲线完全不同,若想和实验光谱对应上,以预测实际光谱或者对实验光谱进行解释,就需要对理论得到的跃迁数据进行“展宽”成为峰形。先对每个跃迁进行展宽,比如图中160nm处的竖线对应S0->S2的跃迁,展宽后就是红色曲线,S0->S11展宽后就粉色曲线(由于跃迁方式很多,为了避免太乱,图中只示意了把振子强度大于0.01的跃迁展宽后的曲线)。把所有跃迁都这么进行展宽成为X-Y曲线后,再把所有曲线的Y值相加,所得到的总的曲线,即图中的黑色曲线,这就是理论预测或者说理论模拟出的光谱图。

吸收曲线的每个峰往往对应于一个强度较大的跃迁。比如图中113nm左右有个峰,很显然就是对应于S0->S13的跃迁(振子强度为0.129)。但是峰的位置和具有较大强度的跃迁的能量却并不总是对应的,无论是实际光谱还是上述方式理论模拟出的光谱都是如此,因为所有跃迁都对附近范围的吸收曲线有贡献。例如上图中S0->S3的跃迁振子强度不算非常小,为0.036,但是黑色曲线在相应位置处(146.3nm)却没有峰。其原因从上图中容易理解,这是因为在S0->S3附近有个振子强度更大的跃迁S0->S5(振子强度为0.107),它对光谱的贡献(青色曲线所示)比起S0->S3的贡献(蓝色曲线所示)大得多,这导致S0->S3没有对应的峰,而被淹没进了最大值在138.2nm处的吸收峰了。这个例子也说明了考察每个跃迁对光谱的各自的贡献的用处,假设我们不把这些振子强度大的跃迁的贡献分别画出来,光谱图的内在结构是不容易搞清楚的,甚至导致对峰的本质的错误指认。这种绘制各个跃迁的独立的贡献的功能是GaussView这样的非专业程序所不具备的,在Multiwfn里则能实现。

上面这个简单的例子讨论的是紫外光谱,这个产生光谱的过程对于红外、ECD(电子圆二色谱)、VCD(振动圆二色谱)也是一致的。这些光谱的理论计算都会给出一个个跃迁模式的能量和强度值,都得把它们各自展宽成曲线并相加才能得到能够和实验结果相对比的光谱。不同类型光谱计算给出的强度有不同的叫法,红外叫红外强度,UV-Vis叫振子强度,VCD/ECD计算叫转子强度(有正有负,对应于左、右两种圆偏振光)。在下文中我们统称为跃迁的“强度”。


PS:这里顺带说一下单位。对于红外、拉曼和VCD谱,常用的单位是cm^-1;对于UV-Vis和ECD,常用的有1000cm^-1、eV和nm。1eV = 8.0655*1000cm^-1。eV为单位的能量的倒数乘以1240.7011就是nm为单位的能量,nm为单位的能量的倒数乘以1240.7011就是eV为单位的能量,因此nm和eV、cm^-1这样的能量单位并不是线性的转换关系。振子强度是没有量纲的,红外强度单位通常是km/mol(千米每摩尔),有时候也用1 esu^2*cm^2 = 2.5066 km/mol为单位。Gaussian输出的ECD转子强度单位是cgs (10^-40 erg-esu-cm/Gauss),VCD转子强度单位是10^-44 esu^2 cm^2。

给定一个跃迁的能量和强度,怎么把它像前面的图所示的那样展宽成峰形的曲线?这需要两个条件,一是展宽函数,它定义了曲线的函数形式;二是半高全宽FWHM(full width at half maximum),它决定了展宽出的峰在高度为一半的位置的峰的宽度。很多地方也用HWHM,是指峰高一半处宽度的一半,HWHM=FWHM/2。显然FWHM越大,峰看起来也就越宽,FWHM越小,则峰看起来就越窄,如果FWHM是一个无限小的值,那么这就不是峰了,而是一个竖线了,相当于没做展宽。下面的图展示了FWHM为不同值的时候的吸收曲线



可见,FWHM=4eV的时候,峰的特征根本看不见,整体就是一个大弧线,显然这样的光谱是没有意义的。FWHM减小为1eV的时候,峰的轮廓已经出来了,但还是比较含糊,或者说此时的光谱分辨率太低,连S0->S11和S0->S13跃迁对应的两个吸收峰都分辨不开。FWHM=0.5eV时,每个强度较大的峰都能看出来了,但是情况和之前一样,S0->S3和S0->S5跃迁区分不开,合并成了一个峰。FWHM进一步降低到0.35eV,这时候光谱分辨率足够高了,S0->S3和S0->S5跃迁都能看到各自的吸收峰了。实验光谱的分辨率是有限的,难以做到非常精细,所以理论模拟光谱的时候FWHM应当选择合适,以便于能够和实验比对。FWHM的选择带有一定任意性,如果不知道怎么设就用程序默认的就好了。顺带一提,对于UV-Vis光谱,如果分辨率特别高,不仅能看出电子态的跃迁,还同时可以看到不同振动态的跃迁,这称为振动分辨的电子光谱,有兴趣者可参考《振动分辨的电子光谱的计算》(http://sobereva.com/223)。

说完了FWHM,再来说展宽函数。常用的展宽函数有下面三个,高斯函数、洛伦兹函数和Pseudo-Voigt函数。



通常,电子光谱(UV-Vis、ECD)用高斯函数来展宽,振动光谱(红外、拉曼、VCD)用洛伦兹函数来展宽,洛伦兹函数衰减得比高斯函数要缓慢。Pseudo-Voigt函数则是高斯函数和洛伦兹函数的线性组合,组合系数是可调的。上面的公式中,ω就是光谱的横坐标,在给定了跃迁能量ω_i和FWHM之后,就立刻有了这个跃迁对应的吸收曲线的表达式。在ω=ω_i的位置处峰最高,而随着ω偏离ω_i函数值逐渐衰减。

上面给出的函数都是归一化的,也就是说函数的积分值为1。展宽时当然也要把跃迁的强度考虑进去。这里有个重点务必记住:跃迁的强度正比于它对应的峰的积分面积。因此,如果A跃迁的强度是0.5,B跃迁的强度是0.1,那么A展宽出的峰的面积也就应当为B的5倍,所以我们要把振子强度乘到展宽函数前面去。值得一提的是,当FWHM相同时峰高正比于峰面积,所以A的峰高也相应地为B的5倍。

我们有了展宽函数的数学形式,也知道了峰面积和跃迁强度的正比关系,展宽出来的吸收曲线怎么才能和实验光谱在定量上对应上?对于红外谱,通过比较红外强度单位和摩尔吸收系数的单位(见Multiwfn手册3.13.1节),可以推得如果以cm^-1为横坐标单位,L/mol/cm为纵坐标单位,则如果一个跃迁的红外强度为p,那么展宽出的曲线面积应当为100*p,也就是说把展宽函数再乘上100*p即可。对于其它类型的光谱则没有这样的形式上的关系,我们只能模拟出光谱的形状,其数值和实际光谱相差一个系数因子,只有通过将大量实验光谱和模拟光谱相对比才能得到这个因子。好在对于UV-Vis有人做过这件事,结论是:如果将1000cm^-1单位用于光谱的横轴,将L/mol/cm单位用于光谱的纵轴,那么1个单位的振子强度展宽出的曲线的积分面积应当为1/4.32*10^6。如果把eV作为横轴单位,则1单位振子强度应当展宽出面积为28700的曲线。有了这个关系,理论模拟的UV-Vis光谱就和实验光谱在定量上就有一定可比性了。不过,对于拉曼、VCD、ECD都缺乏这样的关系,所以只能简单地让p强度的跃迁展宽出积分面积为p的峰,得自己调刻度系数来使模拟的谱和实验谱吻合。

量子化学理论计算出的跃迁都是一个个离散的,为什么实测的光谱是连续的,而不是仅在入射频率恰当与激发能精确一致的位置才有才有吸收,从而观测到离散的谱线?这个问题在一些分子光谱书里都有介绍,导致谱线具有有限宽度的原因有很多,比如(1)不确定原理导致的加宽,即激发态寿命是有限的,故而激发态能量有不确定性(2)由于分子的运动产生的多普勒效应导致的加宽(3)分子间碰撞产生能级位移导致的加宽(4)高辐射强度导致低能态的布居数耗尽导致的饱和加宽(5)柔性分子具有大量可及构象。

理论模拟的光谱和实验光谱常有一定整体的偏差,为了能够尽量相符,我们往往需要一些调节。一是对光谱的高度进行scale,即乘上刻度系数,使模拟光谱的峰高能和实验光谱有较好的对应。通常想算准光谱的强度比起算准峰的位置更为困难,能定性符合就不错了,而且如上所述模拟光谱和实验光谱的本来就缺乏理论上的对应关系,所以做这样的高度的scale完全是合理且也是必要的。另外就是对模拟光谱的横坐标也进行scale或整体加减一个数值,以消除跃迁能量计算的系统性的偏差。比如CIS算的激发能通常偏高,一些研究中会被乘上0.72来修正,而振动光谱众所周知也需要乘上频率校正因子以解决计算方法的系统误差并同时等效地考虑非谐振效应,这点可参见《谈谈谐振频率校正因子》(http://sobereva.com/221)。另外,有时候还需要调节FWHM和展宽函数使结果更好地接近实验谱。这类调节并不算是弄虚作假,因为涉及到的问题是难以克服或者根本不可能克服的,这只是采取一些技巧以便于更好地分析和解释实验光谱。

2.2 拉曼光谱

拉曼光谱是散射光谱,横轴是散射光相对于入射光的频率,纵轴是散射光的强度。量子化学程序直接算出来的是每个振动模式的拉曼活性(Raman activity),单位一般是Å4/amu(amu是原子质量单位)。拉曼活性是每个振动模式自身的特征,和入射光频率及温度无关。原理上来说,获得模拟的拉曼光谱应当先把每个振动模式的拉曼活性转换为每个振动模式的拉曼强度(Raman intensity),拉曼强度是依赖于入射光频率和温度的。转换关系如下:

其中S_i、I_i、νi是第i个振动模式的拉曼活性、拉曼强度和振动频率(波数)。ν0是入射光频率(波数)。T是温度。h是普朗克常数,c是光速,k是玻尔兹曼常数。C是个常数,这无关紧要,我们感兴趣的并非吸收峰的绝对高度,故随便选取个值来让谱图纵轴数量级合适即可。

这样转换成拉曼强度后,再按照上一节所述方式通过洛伦兹函数进行展宽,就得到了可以和实验对照的拉曼光谱。注意,Gaussview等程序绘制拉曼光谱时并不先转换出拉曼强度,是直接基于拉曼活性进行展宽,这么做显然是不严格的。不过,基于拉曼强度和活性展宽出的图的峰的位置都一样,而且每个峰的拉曼强度正比于拉曼活性,因此直接基于拉曼活性展宽出的谱图倒也并非没意义,和实际拉曼光谱在基本特征上还是相似的。但显然,基于拉曼活性和拉曼强度展宽出的谱图的具体形状还是有所不同,峰高不同。
 

3 Multiwfn绘制光谱的输入文件

3.1 Gaussian03/09输出文件

Multiwfn能够直接从Gaussian03/09的输出文件中读取跃迁能和强度数据用于作图。绘制不同的光谱图所要求的关键词如下
(1)对于红外光谱,需要写freq关键词。从G09 D.01开始,freq=anharmonic非谐振计算不仅给出非谐振频率还给出对应的强度,对于这样的Gaussian输出文件,Multiwfn在读入时会提示你要载入非谐振的数据还是谐振的数据。
(2)对于拉曼光谱,写freq=raman关键词就可以了。(由于长期以来以讹传讹,无数人居然误以为Gaussian的freq任务默认就计算raman,还每次都特意写上noraman以为这样会节省时间,这是弥天大误,Gaussian只对HF的freq任务才默认计算raman!!!)
(3)对于VCD光谱,写freq=VCD关键词就可以了。
(4)对于UV-Vis和ECD光谱,直接用普通的TDDFT、TDHF、CIS、ZINDO的输出文件即可,不需要再加其它关键词。ECD的转子强度有两种表象,一种是长度表象,另一种是速度表象,前者依赖于原点而后者则不依赖,在完备基组下它们是一致的,但是在有限基组下则有些偏差,尽管定性一致,Gaussian会同时给出二者,在Multiwfn中可以选择读取哪种转子强度(GaussView作ECD图时用的是速度表象的转子强度)。

3.2 ORCA输出文件

Multiwfn能够直接从ORCA 3.0.3(其它版本未经测试)的输出文件中读取数据绘制IR、Raman、UV-Vis和ECD光谱。
(1)对于红外光谱,写freq关键词即可。
(2)对于Raman光谱,用诸如这样的关键词:
! b3lyp def2-SVP numfreq
%elprop Polar 1 end
(3)对于UV-Vis和ECD光谱,用诸如这样的关键词(CIS和ZINDO/S计算也支持):
! b3lyp def2-SVP
%tddft nroots 20 end
ECD转子强度在ORCA中用的是长度表象。

3.3 文本文件

考虑到很多人不是Gaussian或ORCA的用户,也为了灵活起见,以便于自定义,Multiwfn也支持从文本文件里直接读取跃迁能和强度,并且以这种方式输入数据还有个额外的好处,就是可以定义每个跃迁各自的FWHM,而不强制要求对所有跃迁在展宽时都用相同的FWHM。这样的文本文件的格式如下:
numdata inptype
energy strength [FWHM]                   //1号跃迁的能量、强度和FWHM。FWHM不是必需要写的
energy strength [FWHM]                   //2号跃迁
energy strength [FWHM]                   //3号跃迁
...
energy strength [FWHM]                   //numdata号跃迁

其中numdata是跃迁的数目,也就是条目数。intype为1时只读取跃迁的能量和强度,FWHM由程序自动设定;intype为2时还同时读取每个跃迁的FWHM。跃迁信息应当按照能量从小到大排列。

对于绘制红外、拉曼和VCD光谱,此文件里的能量和FWHM的单位必须都是cm^-1,而对于UV-Vis和ECD则必须以eV为单位。此文件里的强度单位对于红外、拉曼、ECD和VCD必须分别是km/mol, Å^4/amu, cgs (10^-40 erg-esu-cm/Gauss), 10^-44 esu^2 cm^2,这也正是Gaussian程序输出的默认的单位。而UV-Vis振子强度本身就是无量纲的。

这种文本文件的一个例子如下。将它作为Multiwfn启动时载入的文件即可。
6 2
  81.32920        0.72170    8.0
 417.97970        3.58980    8.0
 544.67320       21.06430    8.0
 583.12940       41.33960    8.0
 678.66900       91.47940    8.0
 867.37410        2.94480    8.0

4 使用方法和选项

启动Multiwfn并载入Gaussian/ORCA输出文件或记录着跃迁信息的文本文件后,进入主功能11,然后选择要绘制的光谱类型,然后就能看见一个菜单,选0就可以立刻绘制出光谱。图中既有模拟出的吸收曲线,对应左侧坐标轴;也显示出了各个跃迁,以一条条离散的竖线表示,数值对应于右侧坐标轴。另外还包含很多其它选项,下面依次介绍。

  -2 Export transition data to plain text file:把跃迁数据导出到当前目录下的transinfo.txt,格式和3.3节介绍的格式一样,因此用户可以自行修改里面的数值然后直接作为Multiwfn的输入文件。
  -1 Show transition data:这个功能把跃迁数据直接显示到屏幕上,例如
  Index  Freq.(cm^-1)  Intens.( km/mol   esu^2*cm^2)
     1     81.32870            0.72170     0.28792
     2    417.97970            3.58980     1.43214
     3    544.67310           21.06430     8.40353
     4    583.12940           41.33960    16.49230
     5    678.66890           91.47940    36.49541
很多人想把Gaussian输出文件里的跃迁数据给批量提取出来,但自己又不会写脚本去实现,那么就可以用Multiwfn的这个功能了。跃迁能量、强度整整齐齐地列出来,将它们从屏幕上直接拷贝下来,就可以很方便地进一步处理了。如果不知道怎么从屏幕上把数据拷贝下来,参见手册5.4节。
  1 Save picture:把光谱图输出到当前目录下,和选项0看到的图一样。图像格式可以通过settings.ini文件里的graphformat参数来修改,图像的宽和高由settings.ini里的graph1Dsize参数决定。(注:默认的是png格式。只有用wmf和pdf格式时,光谱图中的灰色网格才能显示出来)
  2 Export X-Y data set of lines and curves to plain text file:把模拟出的光谱曲线以X-Y数据点形式导出到当前目录下的spectrum_curve.txt中,而离散的竖线则导出到spectrum_line.txt里。把这两个文件导入到Origin之类的程序里作图,就可以得到和Multiwfn里产生的图同样的效果。由于Origin这样的程序可调选项更为丰富,因此如果你觉得Multiwfn直接给出的图的效果不满意,可以基于这两个文件在自己擅长的绘图软件中作出完全符合自己要求的图。
  3 Set lower and upper limit of X-axis:设定光谱图的X轴上、下限和刻度间隔,默认是自动确定。自行设定的上下限可以反过来,比如既可以输入0,1000,100也可以输入1000,0,-100,只不过图像左右反转了。
  4 Set left Y-axis、5 Set right Y-axis:设定左、右侧Y坐标轴。需要输入起始值、终止值和刻度间隔。
  6 Select broaden function:选择展宽函数,包括Gaussian、Lorentzian和Pseudo-Voigt。
  7 Set scale ratio for curve:光谱图的刻度因子,设为k的话,那么最后光谱图的数值就会被乘以k。
  8 Input full width at half maximum (FWHM):设定FWHM。
  9 Toggle showing discrete lines:在绘制的光谱图中是否把表示跃迁的离散的线显示出来。
  10 Switch the unit of infrared intensity / Set the unit of excitation energy:对于红外,在km/mol和esu^2*cm^2两个红外强度单位之间切换。对于UV-Vis和ECD,设定光谱图中的能量单位,可以为eV、nm或1000cm^-1。
  11 Set Gaussian-weighting coefficient:当展宽函数被选为Pseudo-Voigt时,就可以用这个选项来设定高斯函数的权重。
  12 Set shift value in X:设定模拟的光谱图的X坐标的位移值
  13 Set colors of curve and discrete lines:设定图中的曲线以及离散的竖线的颜色。
  14 Multiply the transition energies (or vibrational frequenices) by a factor:对跃迁能量乘以一个因子。可以指定只对某些跃迁乘上因子。
  15 Output contribution of individual transition to the spectrum:需要输入一个阈值,比如k。然后就像选项2一样会输出spectrum_curve.txt和spectrum_line.txt,只不过这次输出的spectrum_curve.txt里面还包含了所有强度大于k的跃迁产生的贡献。
  16 Find the positions of local minima and maxima:列出光谱图中的极大值和极小值的位置和数值,这对于确定峰的准确位置很有用。
  17 Toggle showing dashed grid lines:切换是否显示谱图上的虚线网格。
  19 Convert Raman activities to intensities:对于绘制拉曼光谱,默认情况下绘制的是基于拉曼活性展宽的光谱。如果选了这个选项,并且输入温度和入射光波数,就会按前述公式把拉曼活性转换为拉曼强度。之后再作图看到的就是基于拉曼强度展宽的结果了。若再次选这个选项,并且输入的频率和温度和之前一样,则会转换成原先的拉曼活性,相当于逆操作。
  20 Modify oscillator strengths/Raman activities/rotatory strengths:手动设定指定跃迁的强度数据。
 
Multiwfn绘制光谱的时候实际上是计算一系列点上的值,在当前的X轴范围内总共计算num1Dpoints个点。num1Dpoints是settings.ini文件里的参数,默认是3000。这个参数设定得越大,光谱越平滑,导出的X-Y数据点也就越多。通常默认数值足够了,再增加也看不出什么变化。

 

5 实例

5.1 红外光谱

这个实例绘制NH3BF3的红外光谱。启动Multiwfn,然后依次输入
examples\spectra\NH3BF3_freq.out   //程序自带的例子文件,是B3LYP/6-31G*下的振动分析
11
1         //红外光谱
0      //在屏幕上显示光谱
立刻可以看到下图



在图上点右键就可以关闭图像。

接下来我们随意地对图像进行一些调整,都不是必要的,纯粹只是示例一下而已
14  //对频率乘上校正因子
all   //选择所有模式
0.9614   //对这些模式的跃迁能乘上B3LYP/6-31G*下的频率校正因子0.9614(低频区域的校正因子和这个值是显著不同的,这里为了省事就没有单独对高频和低频部分单独乘上各自的校正因子)
9  //不显示离散的竖线
8  //修改FWHM
20  //改为20cm^-1
4 //修改左侧坐标轴
2400,-200,-200  //Y轴下端为2400,上端为-200,步长为-200。此时得到的图是上下反转的,比较符合实验图谱的习俗
0   //重新作图
结果如下所示



选择16,就可以把图中的极大点和极小点找出来,也即得到了峰的准确位置。极大点如下所示
 Local maximum X:      3439.8133      Value:       238.9983
 Local maximum X:      3321.1070      Value:        45.7052
 Local maximum X:      1629.8766      Value:       187.7325
 Local maximum X:      1307.1024      Value:       757.0280
 Local maximum X:      1252.4175      Value:      2203.5450
 Local maximum X:       849.6165      Value:       705.3174
 Local maximum X:       774.9250      Value:        40.1196
 Local maximum X:       637.5458      Value:       457.0291
 Local maximum X:       432.1440      Value:         5.0576
 Local maximum X:       405.4685      Value:         3.4018
 Local maximum X:       252.0840      Value:        74.3636
那些数值比较小的峰我们可以无视掉。从数据可见,最强的峰的位置在1252cm^-1处。

注意极值点的位置的定位精度和前面提到的num1Dpoints参数有关,此参数值越大,数据点的间距越小,对于特别是尖锐的峰的描述越精确,峰的最大值的位置的定位精度也就越高。
 

5.2 拉曼光谱

和绘制红外光谱过程没什么区别。用第3节方式得到的Gaussian/ORCA的输出文件或文本文件作为Multiwfn的输入文件,然后进主功能11,然后选Raman,再选0绘图即可。比如苯胺的图

上面的图是直接基于拉曼活性展宽产生的,若要更严格地绘制拉曼光谱,则应当先选择19,输入入射光波数和温度来转换成拉曼强度,然后再选0看到的就是基于拉曼强度展宽出的拉曼光谱了。

5.3 UV-Vis光谱

这一节我们我们用乙酸的TDDFT的输出文件为例子绘制UV-Vis谱,这也正是本文第二节讨论的那个光谱,用的是程序自带的例子文件,是B3LYP/cc-pVDZ TD(nstates=15)的输出。如果需要得到更完整的谱,nstates应该更大,但这里不追求得到准确的谱,只是示例罢了。启动Multiwfn后输入
examples\spectra\acetic_acid_TDDFT.out
11
3
0
然后就看到了图。

这里介绍下怎么把数据导出来用Origin程序来重新绘制,这样专业的绘图程序在绘图时拥有更多更灵活的可调选项。选择2,曲线数据就导出到了当前目录下spectrum_curve.txt下面,屏幕上也明确提示了其中每一列数据是什么含义。这里第一列是波长(nm),第二列是摩尔吸收系数。将spectrum_curve.txt直接拖到Origin窗口里(本文用的Origin 8),选择绘制曲线图,把A、B两列分别当成X和Y轴数据即可。Multiwfn还同时输出了spectrum_line.txt,把这里面的两列数据分别当成X和Y轴在Origin里绘制成曲线图,就得到了一条条离散的竖线表示出跃迁能和强度(实际上,每个竖线对应于一个往返的折线,比如对(5,0),(5,3.2),(5,0)这三个点绘制成曲线图,看起来就是在X=5的位置上出现了高度为3.2的竖线)。若在曲线图的那张图上点右键选New layer-(Linked) Right Y,在新增的这个Layer上绘制spectrum_line.txt的数据,经过一些细微的调整,就能得到下面的图


如果想绘制出第二节给出的那种显示每个跃迁各自贡献的图,需要选15,然后输入阈值,强度大于这个阈值的跃迁对每个点的贡献会被输出到spectrum_curve.txt的第三列及之后的列当中,从屏幕上可以看到列的编号和跃迁模式的编号的对应关系:
Column#   Transition#
     3           2
     4           3
     5           5
     6          11
     7          13
例如第6列的数据对应于基态到第11激发态的跃迁展宽出的曲线。而spectrum_curve.txt的前两列数据,以及同时输出的spectrum_line.txt,和选项2输出的完全一样。把这两个文件里的数据都放到Origin里作图,就能得到第二节的那种图了。examples\spectra目录下的acetic_acid_TDDFT.opj是那幅图的Origin 8的.opj文件,如果不知道怎么绘制可以直接参考这个文件。

我们还可以计算每个峰由每个跃迁的贡献量。首先我们用选项16来找出吸收曲线的每个峰的准确位置。对于此例,程序找出的极大值的位置,即峰的位置如下
Local maximum X:       113.1515      Value:      5680.8733
Local maximum X:       122.6050      Value:      8239.3520
Local maximum X:       138.1754      Value:      4667.8074
Local maximum X:       159.8627      Value:      2123.3503
Local maximum X:       213.9266      Value:        48.5312
假设这里我们要考察138.1754 nm的那个峰的成份。于是在spectrum_curve.txt里找到与这个波长对应的那一行,内容为
    138.17536   4667.80738      0.14608    310.33658   4272.69585      2.31729      0.00000
其中4667.80738显然就是这个峰的高度。第4列(S0->S3)贡献的值为310.33658,第5列(S0->S5)贡献的值为4272.69585,因此我们就可以说138.1754 nm的峰包含了310.33658/4667.80738=6.65%的S0->S3成份,以及4272.69585/4667.80738=91.53%的S0->S5成份。

顺带一提怎么用Multiwfn绘制荧光光谱。根据Kasha规则荧光是从S1发射的,因此和绘制吸收光谱的关键不同点是要把S0至S1以上激发态的振子强度都设为0。做法就是先优化出S1的结构,然后在这个结构下做一次电子激发计算(默认的nstates=3实际上已经够了),然后将输出文件载入Multiwfn,依次选11、3。然后进选项20,选择除了第一激发态以外的所有的态(比如算了15个态,就输入2-15),然后输入0使得它们的振子强度为0。然后再选0,看到的就是荧光光谱了。

5.4 振动圆二色谱(VCD)

这里来绘制(2S)-2-methyloxirane的VCD谱,用的是程序自带的例子文件,是B3LYP/6-31G*下由freq=VCD关键词计算得到的。启动Multiwfn后输入
examples\spectra\(2S)-2-methyloxirane_VCD.out
11
5    //绘制VCD
0    //显示光谱
假设我们感兴趣的是700~1700cm^-1部分,想只看这部分,于是选3,输入1700,700,再选0重新作图,得到下面的图



5.5 电子圆二色谱(ECD)

这里用天冬酰胺为例子绘制ECD谱。启动Multiwfn后输入
examples\spectra\Asn_TDDFT.out   //PBE1PBE/6-311G* TD(nstates=30)计算的输出
11
4   //绘制ECD
2   //读取速度表象的振子强度
0   //绘制图像
结果如下


两篇使用Multiwfn绘制天然产物ECD光谱的文章见Org. Lett.,15,3526(2013)和J. Nat. Prod.,77,346(2014)。



6 玻尔兹曼分布对光谱的影响


对于柔性分子,势能面有很多极小点。在有限温度下,分子并非只处于特定构象,而是在许多构象上都有一定的分布比例,当然能量越高的构象分布比例越低。分布比例可以通过玻尔兹曼分布公式来计算,具体方法见《根据Boltzmann分布计算分子不同构象所占比例》(http://sobereva.com/165)。不同构象下分子的光谱往往是存在差异的,因此对于柔性分子若想较为准确地模拟出它的实际光谱,就必须把构象的权重考虑进去。比如某有三种构象A、B、C,而且分布比例都不小,那么应当分别把A、B、C的光谱都用Multiwfn导出成文本文件,然后都导入到Origin里,把A、B、C的曲线数据乘上各自的玻尔兹曼分布比例再相加,然后再作图即可。

添加新评论