Gaussianで染料の「色」を当てる ― TD-DFTでインディゴ／アリザリンのUV-VISスペクトルを計算する
「ジーンズの青はインディゴ、赤い染料といえばアリザリン」。教科書で耳にする話ですが、これらの色は、量子化学計算でどこまで再現できるのでしょうか？
本記事では、Gaussian の TD-DFT 計算でインディゴとアリザリンの UV-VIS 吸収スペクトルを求め、計算スペクトルから予測される色（補色）が、実際の染料の見た目と一致することを確かめます。
入力ファイル・出力の読み方・スペクトル描画用 Python スクリプトまで一通り示すので、Gaussian を触ったことがある方ならそのまま再現できます。
1. なぜ「吸収スペクトル」から「色」が分かるのか
可視光は約 380～780 nm。物質に当たる白色光のうち、ある波長帯が吸収されると、私たちは吸収された光の補色を物質の色として認識します。
| 吸収帯（およそ） | 吸収される色 | 見える色（補色） |
|---|---|---|
| 380–450 nm | 紫 | 黄緑〜黄 |
| 450–495 nm | 青 | 黄〜橙 |
| 495–570 nm | 緑 | 赤紫 |
| 570–590 nm | 黄 | 青 |
| 590–620 nm | 橙 | 青緑 |
| 620–780 nm | 赤 | 緑〜青緑 |
つまり、UV-VIS スペクトル上で 「どの波長が吸収されていないか」 を見れば、その物質が呈する色を理屈で予測できます。
2. TD-DFT で何が出てくるのか
UV-VIS の吸収は、分子内の電子が基底状態から励起状態へ遷移するときに起こります。これを計算で扱う代表的手法が 時間依存密度汎関数法（TD-DFT） です。
TD-DFT 計算からは、励起状態ごとに以下が得られます：
- 励起エネルギー $\Delta E_i$（eV または nm）
- 振動子強度 $f_i$（無次元、遷移の起きやすさ）
ただし TD-DFT が直接返すのは「線スペクトル」（離散的なピーク）なので、現実の幅広いスペクトルにするには、各遷移をガウシアン関数で 畳み込み（broadening） して足し合わせます。モル吸光係数 $\varepsilon$ は次式で表せます1：
$$
\varepsilon_i(\tilde{\nu}) = \frac{\sqrt{\pi}\cdot e^{2}\cdot N_\mathrm{A}}{1000 \cdot \ln(10) \cdot c^{2}\cdot m_\mathrm{e}}\ \frac{f_i}{\sigma}\ \exp\left[-\left(\frac{\tilde{\nu} - \tilde{\nu}_i}{\sigma}\right)^{2}\right]
$$
ここで $\tilde\nu$ は波数（cm⁻¹）、$\sigma$ は半値幅、$f_i$ は振動子強度、$N_\mathrm{A}$ はアボガドロ定数。本記事では $\sigma = 0.4$ eV を採用します（実測スペクトルとよく合う経験値）。
3. 計算の流れ
実施するステップは以下のとおりです：
- SMILES から3D初期構造を生成（OpenBabel など）
- Gaussian で B3LYP/6-31G(d) による構造最適化
- 最適化構造に対して TD-B3LYP/6-31G(d) の励起状態計算（励起状態数 20）
- 出力ファイルから $(\Delta E_i, f_i)$ を抽出
- Python でガウシアンを畳み込みスペクトル描画
- 補色と実際の色を比較
計算リソースの目安：本記事の系（30〜40原子）なら、4コア・8GB メモリの計算機で最適化〜TD-DFT 合わせて 1分子あたり10分程度で終わります。
4. 入力ファイル
4-1. 初期構造の準備
SMILES から3D構造を生成し XYZ で保存します（OpenBabel を使用する例）：
# インディゴ
obabel -:"O=C1C(\Nc2ccccc12)=C3/Nc4ccccc4C3=O" -O indigo_init.xyz --gen3D
# アリザリン (1,2-dihydroxyanthraquinone)
obabel -:"Oc1ccc2C(=O)c3ccccc3C(=O)c2c1O" -O alizarin_init.xyz --gen3D
GaussView や Winmostar でも同様に作れます。生成した座標を、後述の .gjf の座標部分に貼り付けてください。
4-2. 構造最適化（B3LYP/6-31G(d)）
%mem=4GB
%nprocshared=4
%chk=indigo_opt.chk
# B3LYP/6-31G(d) Opt Freq
Indigo geometry optimization
0 1
C ... ... ...
N ... ... ...
（初期構造として作成した3D座標を貼り付け）
Freq を付けて虚振動が無い（極小点である）ことを必ず確認します。アリザリンも同様に alizarin_opt.gjf を用意します。
ファイル名の使い分けに注意: アリザリン用の入力では、%chk=alizarin_opt.chk のように .chk の名前も分子ごとに変えてください。インディゴと同じ indigo_opt.chk のままにすると、ファイルが上書きされたり、次の TD-DFT 計算（%oldchk で参照）で別の分子の構造を読み込んでしまうため、結果を取り違えます。後段の *_td.gjf でも %oldchk / %chk を対応する分子名に合わせます。
4-3. TD-DFT 計算（励起状態 20 個まで取得）
最適化後の .chk を読み込んで TD-DFT のシングルポイントを実行します：
%mem=4GB
%nprocshared=4
%oldchk=indigo_opt.chk
%chk=indigo_td.chk
# TD(NStates=20)/B3LYP/6-31G(d) geom=check guess=read
Indigo TDDFT
0 1
NStates=20 は、可視〜近紫外領域までの主要遷移を漏らさないための余裕を持った設定です。系が大きい場合は 30〜100 まで増やします。
座標について: geom=check を指定しているため、入力ファイルに座標は記述しません。最適化後の構造が indigo_opt.chk から自動で読み込まれます。UV-VIS（垂直励起）は 基底状態の安定構造 で計算するのが正しい手順なので、必ず opt 後の構造を使います（Franck-Condon 原理：電子遷移は核が動かない瞬間に起こるとみなす）。
.chk を使わずに座標を直接書く場合: %oldchk と geom=check guess=read を外し、電荷・スピン多重度の行（0 1）の下に座標を貼り付けます。その際は opt 計算後の最適化座標（opt の .log 末尾の "Standard orientation" などの座標）を使ってください。初期構造（最適化前）ではありません。
Tip: geom=check guess=read で最適化済み構造と波動関数を引き継ぐと、TD 計算が高速かつ安定します。
5. 出力からの励起情報の取り出し
Gaussian の出力（.log）には、励起状態ごとに次のような行が出ます：
Excited State 1: Singlet-A 2.0541 eV 603.61 nm f=0.2913 <S**2>=0.000
これを正規表現で抜き出すスクリプト：
import re
def parse_td(logfile):
"""Gaussian TDDFT log から (励起エネルギー[eV], 振動子強度 f) のリストを返す"""
results = []
pattern = re.compile(
r'Excited State\s+\d+:\s+\S+\s+([\d.]+)\s+eV\s+[\d.]+\s+nm\s+f=([\d.]+)'
)
with open(logfile) as fp:
for line in fp:
m = pattern.search(line)
if m:
results.append((float(m.group(1)), float(m.group(2))))
return results
if __name__ == '__main__':
import sys
for E, f in parse_td(sys.argv[1]):
print(f'{E:8.4f} eV f = {f:7.4f}')
実行例：
$ python parse_td.py indigo_td.log
2.0541 eV f = 0.2913
3.4178 eV f = 0.0021
3.7625 eV f = 0.0044
...
6. ガウシアン畳み込みでスペクトル描画
各励起をガウシアンで畳み込んで波長軸（nm）の $\varepsilon(\lambda)$ にし、第一軸（左）にモル吸光係数 $\varepsilon$、第二軸（右）に各遷移の振動子強度 $f$（スティックスペクトル） を重ねて描画します：
import numpy as np
import matplotlib.pyplot as plt
from parse_td import parse_td
# パラメータ
SIGMA_EV = 0.4 # ガウシアンの半値幅 (eV)
EV_TO_CM = 8065.544 # 1 eV = 8065.544 cm-1
EV_TO_NM = 1239.84193 # λ[nm] = 1239.84/E[eV]
# f↔積分吸収の関係 f = 4.319e-9 ∫ε dν~ に由来する標準係数。
PREFACTOR = 1.3062974e8 # = (1/4.319e-9)/√π [L mol-1 cm-1 ・ cm-1]
def spectrum(excitations, x_nm):
"""波長配列 x_nm に対するモル吸光係数 ε(λ) を返す"""
sigma_cm = SIGMA_EV * EV_TO_CM
x_cm = 1e7 / x_nm # 波数 (cm-1)
eps = np.zeros_like(x_nm)
for E_eV, f in excitations:
nu_i = E_eV * EV_TO_CM # 遷移の波数 (cm-1)
eps += f * np.exp(-((x_cm - nu_i) / sigma_cm) ** 2)
return eps * PREFACTOR / sigma_cm
# 励起情報の読み込み
indigo = parse_td('indigo_td.log')
alizarin = parse_td('alizarin_td.log')
# 300–800 nm
x = np.linspace(300, 800, 1201)
e_indigo = spectrum(indigo, x)
e_alizarin = spectrum(alizarin, x)
# 描画
fig, ax = plt.subplots(figsize=(8, 5))
# 可視光帯の色帯（補色を直感的に見るため）
visible_bands = [
(380, 450, 'violet'),
(450, 495, 'royalblue'),
(495, 570, 'mediumseagreen'),
(570, 590, 'gold'),
(590, 620, 'orange'),
(620, 750, 'red'),
]
for lo, hi, c in visible_bands:
ax.axvspan(lo, hi, color=c, alpha=0.12)
# 第一軸（左）: モル吸光係数 ε（ブロードニング済みの連続スペクトル）
ax.plot(x, e_indigo, label='Indigo', color='navy', lw=2)
ax.plot(x, e_alizarin, label='Alizarin', color='crimson', lw=2)
ax.set_xlabel('Wavelength (nm)')
ax.set_ylabel(r'$\varepsilon$ (L mol$^{-1}$ cm$^{-1}$)')
ax.set_xlim(300, 800)
ax.set_ylim(bottom=0)
ax.legend(loc='upper right')
ax.grid(alpha=0.3)
# 第二軸（右）: 振動子強度 f（各遷移をスティックで表示）
ax2 = ax.twinx()
for excitations, color in [(indigo, 'navy'), (alizarin, 'crimson')]:
lam = np.array([EV_TO_NM / E for E, f in excitations]) # 各遷移の波長 (nm)
osc = np.array([f for E, f in excitations]) # 振動子強度
ax2.vlines(lam, 0, osc, color=color, alpha=0.45, lw=1)
ax2.set_ylabel('Oscillator strength $f$')
ax2.set_ylim(bottom=0)
fig.tight_layout()
fig.savefig('uvvis_indigo_alizarin.png', dpi=150)
実行方法: parse_td.py・plot_uvvis.py・indigo_td.log・alizarin_td.log を同じフォルダに置き、必要なライブラリを入れてから実行します。
pip install numpy matplotlib # 初回のみ
python plot_uvvis.py
引数は不要です（ファイル名はスクリプト内に記述済み）。成功すると同じフォルダに uvvis_indigo_alizarin.png が出力されます。
plot_uvvis.py は冒頭で from parse_td import parse_td としているため、parse_td.py を同じフォルダに置いてください。ログファイルも相対パスで読み込むので、同フォルダに置くのが簡単です。
スクリプト中のスケーリングファクター 1.3062974e8 は、振動子強度と積分吸収を結ぶ実用関係式 $f = 4.319\times10^{-9}\int\varepsilon(\tilde\nu),d\tilde\nu$（$\varepsilon$ は L mol⁻¹ cm⁻¹、$\tilde\nu$ は cm⁻¹）に由来します。ガウシアン1本の面積が $\varepsilon_\mathrm{max},\sigma\sqrt{\pi}$ であることから $\varepsilon_\mathrm{max} = \dfrac{1}{4.319\times10^{-9}\sqrt{\pi}}\cdot\dfrac{f}{\sigma} = 1.3062974\times10^{8}\cdot\dfrac{f}{\sigma}$ となり、$\sqrt{\pi}$ はこの係数にすでに含まれています（GaussSum / Multiwfn と同じ定義）。
7. 結果
得られたスペクトル（イメージ：実際の計算結果は元記事を参照）：
図の出典：色素分子のUV-VISスペクトル ― TSテクノロジー テクニカルノート
上図および本記事の計算条件（TD-B3LYP/6-31G(d)、$\sigma = 0.4$ eV）は、TSテクノロジーのテクニカルノートを元にしています。各分子の構造・3Dモデル・実サンプル写真を含む詳細はリンク先をご覧ください。
TD-DFT の理論は M. E. Casida, J. Mol. Struct. THEOCHEM 914, 3 (2009)、broadening の式は GaussSum / Multiwfn の各マニュアルに準拠しています。
7-1. インディゴ（青色染料）
- 長波長側のピーク：約 600 nm 付近（橙〜赤の吸収）
- 380〜430 nm の紫の領域は吸収が小さい
→ 吸収の「窓」が紫〜青に開いているため、補色である紺〜青紫色を呈する。実際のインディゴの「藍色」と一致。
インディゴのサンプル 写真出典：Wikimedia Commons
7-2. アリザリン（赤色染料）
- 約 430 nm 付近にピーク（青〜紫の吸収）
- 500 nm 以降の緑〜橙〜赤の領域は吸収が小さい
→ 補色である橙〜赤色を呈する。実際のアリザリンの赤と一致。
アリザリンのサンプル 写真出典：Wikimedia Commons
7-3. まとめると
| 化合物 | 主吸収帯（計算） | 吸収の「窓」 | 予測色（補色） | 実際の色 | 一致 |
|---|---|---|---|---|---|
| インディゴ | ~600 nm | 紫〜青 | 紺〜青紫 | 藍色 | ✓ |
| アリザリン | ~430 nm | 緑〜赤 | 橙〜赤 | 赤 | ✓ |
TD-B3LYP/6-31G(d) という比較的軽い計算レベルでも、染料の色の傾向はきちんと再現できることが分かります。
8. 注意点と精度向上のヒント
今回の計算は気相・単一構造の計算ですが、実用上は以下を意識すると、より定量的に実験と合います。
8-1. 溶媒効果（PCM）
実測 UV-VIS は溶液中。Gaussian なら 1 行追加するだけで PCM（連続誘電体モデル）が使えます：
# TD(NStates=20)/B3LYP/6-31G(d) SCRF=(PCM,solvent=ethanol)
極性溶媒では n→π* 遷移が短波長シフトしたり、π→π* が長波長シフトするので、観測されている溶媒系に合わせるのが望ましいです。
8-2. 汎関数の選択
B3LYP は安価で扱いやすいですが、電荷移動遷移（CT）やライドベルク状態で励起エネルギーを過小評価しがちです。問題があるときは：
- CAM-B3LYP や ωB97X-D（長距離補正型）
- M06-2X（メタ GGA、励起状態でも比較的良好）
を試すと良いでしょう。インディゴの長波長ピークもメチン鎖を含む CT 性があるため、CAM-B3LYP との比較は面白いと思います。
9. まとめ
- Gaussian の TD-B3LYP/6-31G(d) という比較的軽い計算で、インディゴ・アリザリンの UV-VIS スペクトルが得られた
- 計算スペクトルから予測される補色が、実際の染料の色と一致することを確認
- 染料・色材設計、光反応の波長設計、新規発色団のスクリーニングなど、「色」を設計する研究に量子化学計算は十分使えるレベルにある
次のステップとして、
- 溶媒効果（PCM）
- 長距離補正型汎関数（CAM-B3LYP/ωB97X-D）
を取り入れていくと、より定量的な議論ができるようになります。
-
$\sigma$ を eV で扱う場合は $\sigma_\mathrm{cm^{-1}} = \sigma_\mathrm{eV} \times 8065.544$ で波数に換算してから用います。 ↩