IRスペクトルを描きたい

 前回、Gaussianをネタにしたので、ついでにこれを書いておきたくなった。
 GaussianではIRとかVis-UVとかNMRとかいったスペクトルをシミュレートする機能がある。そこで得られたスペクトルというのは実測するのとは大分違った形で表記される。今回はIRを例にとって話を進める。
 Gausianプログラムで学ぶ情報化学・計算化学実験には、アセトアルデヒドのIRスペクトルの計算値と実測値との比較と称して次の表が載っている。

 折角なのでこれを使う。表のB3LYP/6-31G(D)にある振動子強度とscaledが必要な値なので、他の部分には目をつぶる。これを書き出すと以下のとおりとなる。

scaled 振動子強度
148.7 0.002806
487.5 0.128902
749.8 0.006655
859 0.082468
1093.1 0.246019
1100.4 0.014986
1348.9 0.185796
1390.5 0.123765
1432.3 0.169904
1441.5 0.093053
1771.6 1.56687
2784.6 1.3846
2926 0.025906
2980 0.093043
3042.1 0.104524

 "scaled"が横軸波数(cm-1)で、"振動子強度"が縦軸となる。振動子強度は縦軸をモル吸光係数としたときのピーク面積に比例する値[2]なのだが、ピークがシャープなのかブロードなのかで見た目がだいぶ変わるため、振動子強度にあまり多くの意味を込めることは出来ない。
 15個数値があるのは、アセトアルデヒドは原子の数が7個あり、分子の基準振動の数は(3n-6)で与えられるため。
 これをグラフにすると次のような棒グラフができる。

 このグラフの描き方はいつものようにアップしたエクセルのシートを見てもらえれば分かるのだけど、Aの列に横軸をテキトーに並べて、BにAの値が"scaled"の値を跨いだかを判別して、跨いでいたらCに振動子強度の値を入力する、跨いでいなければテキトーなマイナスの値を入力するとしている。
 この描き方はXRDで既知のチャートを一緒に表示するときなんかにも使える。
 論文とかの硬派なメディアで開陳するならこのきっちりと理想の数値をひとつだけ表す棒グラフがよいのだけど、ブログとかの軟派な場ではもっとキャッチーな表現をしたくなる。具体的には普通に測定したIRスペクトルのようなものを描きたい。
 Gauss Viewなんかを使えばそれらしいものを見せてくれるのだけど、マニュアルを見るとピークの幅なんかはテキトーに合わせているらしい。それなら自分でもエクセルで作れるんじゃないかと思ってやってみた。
 とはいっても、それぞれのピークについて正規分布のグラフを作って、最後に全部足し合わせるだけなので特に苦労はない。
 正規分布は、
{ \displaystyle f(x) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\lbrace- \frac{(x- \mu )^2}{2 \sigma^2} \rbrace }
で示される。
 x=μのときにexpの中は0となるので、e0=1であり、この関数は最大値、つまり(2πσ2-0.5となる。ここから、μとσ2を求めることができる。
 強度をI、ピークの波数をνとすると、
{ \displaystyle \frac{1}{\sqrt{2 \pi \sigma^2}} = I \\ \displaystyle \frac{1}{2 \pi \sigma^2} = I^2 \\ \displaystyle \sigma^2 = \frac{1}{2 \pi I^2} }
μ=ν
となる。
 これでそれっぽいグラフになるのだけど、普通IRは縦軸を透過率で表すので、そういう形にしたい。
 縦軸の値を100から引くだけでは強度が弱すぎてスペクトルが見えてこないので、最も吸収の大きい波数のところが15%くらいになるように強度を補正する。さらに、ベースラインが100%のところにあるのはあまりにも不自然だから98%がベースラインとなるように全体を縮めてみると次のようにIRのスペクトルを描くことができる。

 あとついでに、強度が変化しないようにσ2に手を加えるとピークをシャープにしたりダルくしたりできる。

 なお、アセトアルデヒドのIRスペクトルの実測データは下の通り[1]

 何をどう帰属して表4-1のようになったのかよくわからんけど、取り敢えず強度はそのままのスケールじゃなくて対数に取ったほうがいいのかなという気はする。Lambert-Beerの法則では対数が出てくるので、強度を対数表記したほうがそれっぽいのかもしれない。でも1以下の数値を対数に入れると何か良くないことが起こりそう。そこら辺はなんとかするとして、それと、振動強度が面積を示すんならルートを付けた値としたほうがいい気がする。

 そんなわけで、取り敢えず対数とルートを取り込んで吸収強度の格差の解消を試みてみた。

 ついでに、実測のグラフに合わせて縦線を追加してみた。
 少しは格差は減った感じはするのだが、見ての通り全然である。
 縦線を追加したせいで実測との数値の違いが明確になっている。
 まあ、現時点(とはいっても、参考資料は10年以上前の書籍だが)での計算性能はこんなもんということである。本エントリーはピークの数値からIRのチャートを描くことだから、計算性能はあんまり考慮に入れる必要はないので、そこは深く追求するところではない。

 いつものようにエクセルファイルを上げておく。
 IR 20171107.zip

参考文献
 [1]NIST Chemistry WebBook(魚拓)
 [2]富士通