HiTRANで透過法のガススペクトルを計算しCSVで出力させる

背景

フーリエ変換型赤外分光光度計（FTIR）を使用して有機ガスの赤外スペクトルを測定した際、ガスの典型的な振動回転スペクトルが得られる。
塩化水素の振動回転スペクトルは学生実験でも採用しやすい典型的な振動回転スペクトルが得られる。(図１)

上記データは日本分光のFT/IR-6100、光路長が10cmのガラス製のガスセル、窓板がNaClで測定したデータ　...だったと思います。当時の実験ノートを見ても窓板は流石に記録してなかったです。
分解能は0.5cm-1、積算回数は32回で測定しています。データ見れば分かりますがDTGS検出器です。

同様に他の有機ガスで回転振動スペクトルを得ようとする。通常はガスセルを使うのだが、その窓板,光路長は種類が大量にあり、 予め選択しなければならない。

選択をするためにはモル吸光係数εと濃度が必要になる。
モル吸光係数εは量子化学計算を行えばある程度推定することができる。
しかし、量子化学計算は計算リソースが必要となり、ソフトによっては有料になる。

ガスの振動回転スペクトルの計算ぐらいであればそこまでリソースがいらないので、gaussianのライセンス料が払えるかが重要。　
GAMESSなどを使えば無料でも出来るが情報が少なめな点には注意。

量子化学計算を使わない方法として市販の定量用ガスライブラリを使う方法もある。
濃度と光路長の情報がライブラリに含まれるためモル吸光係数εを算出するやり方である。
これは有料ではあるが上記よりかは安価に行える。

しかしライブラリを使った方法ではアサイメントが難しくなる。そこでアサイメントとスペクトルの計算の両方ができるデータベース「HiTRAN」を使った方法を記述する。

HiTRANそのものは多くの論文によるデータベースではあるものの、アサイメントを含む論文をpublishする際はGaussian等の量子化学計算をリファレンスにしたものにすることを推奨します。

計算方法　

以降の作業は各種公式のリファレンスマニュアルを読むべきです。
コードはあくまで一例です。
Windowsでの作業を想定しています。Macでも原理的には動作するはずですが未検証です。

Python環境構築

Python公式ページからPythonのインストーラーをダウンロードし、インストール作業を行う。
インストールの際にパスを通すことを忘れずに。

Pythonは人気の言語なのでQiitaで検索すればより詳しいインストール方法が学べるはずです。HiTRANではPythonの実行環境が必須です。

インストールが嫌ならばGoogle Colaboratoryの使用をオススメいたします。

パッケージのインストール

HiTRAN公式ページから hapi.pyをダウンロードし、pythonのコードを置いている,置くつもりのディレクトリに設置する。

...もしくは
Google Colaboratoryなどローカル以外の環境でPythonを動かしている場合は以下のpipコマンドを使用してインストールする。 Google Colabの場合は!pip install hitran-apiと打つ

pip install hitran-api

importする

HAPIは以下のようにimportすれば使える。

from hapi import *

HiTRANのデータベースへの接続は以下の文で行う。

db_begin('.')

データベースへ接続後、fetch関数に任意の引数を入れて呼び出すことで必要なデータをダウンロードすることができる。

fetch('HCl', 14, 1, 400, 4000)

引数

第一引数　データ名　
上記の例では取得したデータはHCl.dataに格納される。

第二引数　moleculeID
HiTRANにおける分子のID　Line-by-Line SearchのページでIDを確認できる。

第三引数　isotopologueID
HiTRANにおける各分子の同位体のID　目的の分子にチェックを入れたあとSelect Isotopologuesをクリックしたあとに表示されるページから確認可能　HClの例

上記の例では原子量が35のHClのデータをダウンロードしている。（存在比75.76%）
実際のスペクトルと比べるときは、原子量35のHClと原子量37のHClのコンボリューションスペクトルを計算して比較するのが望ましい。

fetch関数を使用してデータを取得した場合はローカルにデータが残るので、同じ分子であれば２回目以降はfetch関数を呼び出さなくても良い。

データを処理する

得られたhogehoge.dataは多量のパラメーターを持ったデータ群である。このうち必要なものを抜き出してからデータ処理する必要がある。
今回は最終が波数対吸光度のグラフを得たいので、波数と吸光度を計算するためのパラメーターが必要になる。
ガスのスペクトルを計算するには、

• 圧力
• 温度
• 濃度

があれば吸収断面積（モル吸光係数）が求められる。

求めるにはabsorptionCoefficient_Lorentz()関数を呼び出す

absorptionCoefficient_Lorentz(SourceTables='HCl', 
                                         HITRAN_units=False, 
                                          WavenumberStep=0.1, 
                                          Environment={'p': 1.0, 'T': 298}, 
                                          Diluent={'air': 1.0-c,'HCl':c})

引数

引数SourceTables　利用するデータ名　
上記の例では取得したデータHCl.dataを利用する

引数HITRAN_units 論理型　TRUE or FALSE
TRUE： cm^2/molec
FALSE: cm-1
吸収スペクトルや透過率などの計算はFALSEに設定する

引数WavenumberStep 数値型
デフォルトは0.01cm-1
あとでSciPyライブラリで分解能のリサンプリングを行う方針であれば、0.001以下に設定する
（そもそもそのやり方の方がより厳密な計算が行える。）

引数Environment
第一引数　圧力[atm]　デフォルト＝1.0
第二引数　温度[K]　デフォルト＝296.0

引数Diluent 混合ガスの定義
A+B+C...=1になるように定義する

装置関数を使用（やりたければ）

装置関数を使用することで実際のスペクトルに近づけることができる。例えば装置のアポダイゼーションが三角関数である場合は以下のような関数を実行してデータを変換する。

nu,coef,i1,i2,Res = convolveSpectrum(nu, coef, SlitFunction=SLIT_TRIANGULAR,Resolution=1.0)

吸光度を計算する

吸光度はLambert-Beer則に基づいて計算する。下記は10.0cmガスセルを使用して濃度が100ppmとした時の計算

absorbance = coef * 10.0 * 1.0e-4

あとはnumpyライブラリ、matplotlibライブラリ、pandasライブラリうを使っていい感じにデータを表示して出力させる。

具体例 塩化水素のスペクトル

HiTRAN_HCl_csv.py

# 必要なパッケージをインポート
import numpy as np
from hapi import *
import matplotlib.pyplot as plt
import pandas as pd

# HITRANデータベースへの接続を開始
db_begin('.')

# 塩化水素（HCl）のデータをダウンロード
fetch('HCl', 15, 1, 400, 4000)

# パラメータを設定
P = 1.0 # 圧力 (atm)
T = 298.0 # 温度 (K)
L = 10.0 # 光路長 (cm)
c = 1.0e-4 # 濃度 (ppm)

# 吸収係数の計算
nu, coef = absorptionCoefficient_Lorentz(SourceTables='HCl', 
                                         HITRAN_units=False, 
                                          WavenumberStep=0.1, 
                                          Environment={'p': P, 'T': T}, 
                                          Diluent={'air': 1.0-c,'HCl':c})

# 装置関数の使用
nu,coef,i1,i2,Res = convolveSpectrum(nu, coef, SlitFunction=SLIT_TRIANGULAR,Resolution=1.0)

#吸光度の計算
absorbance = coef * L * c

# グラフを作成
plt.figure(figsize=(10, 6))
plt.plot(nu, absorbance)

# 軸のラベルを設定
plt.xlabel('wavelength (cm-1)')
plt.ylabel('Absorbance (Abs.)')

# グラフのタイトルを設定
plt.title('Absorbance spectrum of HCl')

# グラフを表示
plt.show()

# データフレームを作成
df = pd.DataFrame({
    'wavelength': nu,
    'Absorbance': absorbance
})

# CSVファイルに保存
df.to_csv('spectrum_HCl.csv', index=False)

得られた計算スペクトルは以下の通りです。

csvファイルも同じフォルダに出力されてていい感じです。

一方で最初の実験スペクトルとのAbs.解離に気が付いた方もいるかもしれません。これは実験でのHClガスの捕集方法として、ガスセルを真空引きして塩酸試薬瓶に近づけて捕集する方法を使用したので、おそらく100ppm以上あると思います。今回は適当に100ppmで計算をしましたが、他の分析方法例えばGC-MSやTOF-MSを使用して濃度をはっきりさせてから計算スペクトルを計算させるとかなり近い計算スペクトルが得られると考えられます