🏄

[#7] QEでH2O(液相)のIRスペクトルを計算してみた:実験値と計算値の比較

に公開

[#7] QEで液相H2OのIRスペクトルを計算してみた:実験値と計算値の比較

液相H2Oの赤外線(IR)スペクトルを調和近似で計算し、実験値との比較を行いました。擬似ポテンシャル(PP)の違い(USPP vs ONCV PP)が振動数に与える影響を評価し、調和近似の限界について考察します。

はじめに:液相H2OのIRスペクトル

本記事の目的: 調和近似によるIRスペクトル計算の精度と限界を、擬似ポテンシャルの違い(USPP vs ONCV PP)を通じて検証します。

実験で観測される特徴

液相水のIRスペクトルは、水素結合ネットワークによって特徴的なパターンを示します。本計算では8 H₂O分子(24原子)のスーパーセルを使用します。

領域 周波数 [cm⁻¹] モード 特徴
Librational 400-800 回転振動 水素結合の拘束運動
Bending ~1645 H-O-H 変角 比較的シャープ
Stretching 3000-3600 O-H 伸縮 ブロードで複数のピーク

補足: 多くの現代文献では液相水のO-H伸縮バンドの主ピークを約3400 cm⁻¹(3450 cm⁻¹や3490 cm⁻¹とする報告も)と記述しています。本計算ではHale & Querry (1973)のデータに基づく最大吸収位置(3389.8 cm⁻¹、2950 nm相当)を基準値としました。Hale & Querryの値はSegelstein (1981)以降の研究でも基盤データとして引用されており、現代の記述との差は測定条件によるもので、本比較の結論に影響しません。

O-H伸縮モードの多峰性

液相水のO-H伸縮バンドは、多様な水素結合環境により複数のピーク成分から成ります:

  • ~3220 cm⁻¹: 強い水素結合(低波数シフト)
  • ~3400 cm⁻¹: 中程度の水素結合(主ピーク)
  • ~3550 cm⁻¹: 弱い水素結合(高波数)

これは液相水が単一のピークではなく、複数の水素結合状態の重ね合わせであることを示しています。

計算方法

調和近似によるIRスペクトル計算

Quantum ESPRESSOを使用して、調和近似(Harmonic Approximation)でIRスペクトルを計算しました。

  1. vc-relax: 構造最適化(セル体積・原子位置)
  2. scf: 自己無撞着場計算
  3. ph.x: フォノン計算(q=0点)
  4. dynmat.x: 動的行列の対角化、IR活性計算

擬似ポテンシャル比較

PP = Pseudopotential(擬似ポテンシャル)

種類 ecutwfc [Ry] 説明
USPP 80 Ultrasoft PP(ウルトラソフト)
ONCV PP 80 Optimized Norm-Conserving Vanderbilt PP

ガウスブロードニング

調和近似では幅ゼロのスティックスペクトルしか得られないため、計算結果を視覚的に比較するためにガウスブロードニングを適用しました。

  • シグマ値: σ = 50 cm⁻¹
  • 注意: これはピーク位置の比較のための可視化用であり、実験スペクトルのブロードさ(非調和性・温度効果による数百cm⁻¹の幅)を定量的に再現するものではありません

入力ファイル・実行手順

本計算で使用した入力ファイルと実行スクリプトの構成を示します。USPP版とONCV版でディレクトリ構成・入力ファイル形式は全く同じであり、ATOMIC_SPECIESセクションの擬似ポテンシャルファイル名のみが異なります。

計算フロー

[1] vc-relax (relax_production.in)
      ↓ 構造最適化(セル体積・原子位置)
      収束条件: conv_thr=1.0d-8, smearing使用

[2] scf (scf_after_vc.in)
      ↓ 波動関数計算(vc-relax後の最適化構造で)
      収束条件: conv_thr=1.0d-7, fixed occupations

[3] ph.x (ph_ir_production.in)
      ↓ フォノン計算(q=0点、Γ点)
      変位パターン: 24原子×3方向=72モード

[4] dynmat.x (dynmat.in)
      ↓ IR活性・振動数計算
      出力: tmp/dynmat.out(振動数・IR強度)

ディレクトリ構成

materials/h2o_liquid/uspp_no_broad_ONCVPP/  (ONCV PP版)
├── inputs/
│   ├── relax_production.in   # 構造最適化(vc-relax)
│   ├── scf_after_vc.in       # SCF計算
│   ├── ph_ir_production.in   # フォノン計算(q=0)
│   └── dynmat.in             # 動的行列解析
├── exec/
│   ├── run.sh                # 実行スクリプト
│   └── plot_spectrum.py      # IRスペクトルプロット
├── tmp/                      # 計算結果(自動生成)
└── outputs/                  # 出力ファイル(自動生成)

materials/h2o_liquid/uspp_no_broad_USPPPP/  (USPP版)
└── (上記と全く同じ構成)

入力ファイル(ONCV PP版)

relax_production.in(vc-relax用):

&CONTROL
    calculation = "vc-relax"
    prefix = "h2o_liquid"
    outdir = "./tmp"
    pseudo_dir = "../../../pseudopotentials/"
    nstep = 500
/
&SYSTEM
    ibrav = 1
    celldm(1) = 12.32       # 8 H2O分子の立方セル
    nat = 24               # 8分子×3原子
    ntyp = 2               # H, O
    ecutwfc = 80.0         # プレーン波基底カットオフ
    ecutrho = 320.0        # 电荷密度カットオフ
    occupations = 'smearing'
    smearing = 'gauss'
    degauss = 0.01
/
&ELECTRONS
    conv_thr = 1.0d-8
    mixing_beta = 0.3
/
&IONS
    ion_dynamics = "bfgs"
/
&CELL
    cell_dynamics = 'bfgs'
    press = 0.0
/
ATOMIC_SPECIES
H   1.008  H_ONCV_PZ_sr.upf    # ONCV PPの場合
O  15.999  O_ONCV_PZ_sr.upf    # ONCV PPの場合
# USPPに変更する場合:
# H   1.008  H.pz-rrkjus.UPF
# O  15.999  O.pz-rrkjus.UPF
ATOMIC_POSITIONS angstrom
O  1.778775  1.588281  1.824067
H  1.257741  1.511703  1.021409
...(全24原子の座標)
K_POINTS automatic
1 1 1 0 0 0                 # Γ点のみ

scf_after_vc.in(SCF計算用):

&CONTROL
    calculation = "scf"
    prefix = "h2o_liquid"
    outdir = "./tmp"
    pseudo_dir = "../../../pseudopotentials/"
/
&SYSTEM
    ibrav = 1
    celldm(1) = 12.32
    nat = 24
    ntyp = 2
    ecutwfc = 80.0
    ecutrho = 320.0
    occupations = 'fixed'   # vc-relax後はfixed推奨
/
&ELECTRONS
    conv_thr = 1.0d-7       # vc-relaxより厳密
    mixing_beta = 0.3
/
ATOMIC_SPECIES
H   1.008  H_ONCV_PZ_sr.upf    # ONCV PPの場合
O  15.999  O_ONCV_PZ_sr.upf    # ONCV PPの場合
# USPPに変更する場合:
# H   1.008  H.pz-rrkjus.UPF
# O  15.999  O.pz-rrkjus.UPF
ATOMIC_POSITIONS angstrom
O  1.8104915455  1.5957526062  1.8566842837
H  1.2607954015  1.5184581526  1.0264044366
...(vc-relax後の最適化座標24個)
K_POINTS automatic
1 1 1 0 0 0

注意: USPP版とONCV版でATOMIC_SPECIESの擬似ポテンシャルファイル名のみが異なります。他のパラメータ(原子座標含む)は完全に同じです。

ph_ir_production.in(フォノン計算用):

&INPUTPH
    prefix = "h2o_liquid"
    outdir = "./tmp"
    fildyn = "./tmp/h2o_liquid.dyn"
    epsil = .true.          # 誘電率・IR活性計算
/
0.0 0.0 0.0                 # q = (0,0,0) Γ点
trans                       # 変位モード
1.0 0.0 0.0 1.0 0.0 0.0 ... # 24原子×3方向の変位パターン

dynmat.in(動的行列解析用):

&INPUT
    fildyn = "./tmp/h2o_liquid.dyn"
    asr = "zero-dim"        # アコースティック和則
    q(1) = 0.0
    q(2) = 0.0
    q(3) = 0.0
/

実行スクリプト(run.sh)

詳細は、ダウンロードで確認することができます。

#!/bin/bash
# 液相H2O調和近似計算(ONCV PP)
#
# 計算フロー:
# [1/4] vc-relax → 構造最適化(セル体積・原子位置)
#       入力: inputs/relax_production.in
#       出力: tmp/relax.out
#
# [2/4] scf      → 自己無撞着場計算
#       入力: inputs/scf_after_vc.in(vc-relax後の構造を使用)
#       出力: tmp/scf.out
#
# [3/4] ph.x     → フォノン計算(q=0点)
#       入力: inputs/ph_ir_production.in
#       出力: tmp/h2o_liquid.dyn(動的行列)
#
# [4/4] dynmat.x → IR活性・振動数計算
#       入力: inputs/dynmat.in(標準入力)
#       出力: tmp/dynmat.out(振動数・IR強度)
#
# 計算時間(Ryzen 9 5700X, 16コア):
#   ONCV PP: 約53分(vc-relax 32分 + ph.x 21分)
#   USPP:    約58分(vc-relax 34分 + ph.x 23分)
#
# 使用コア数: 16コア(環境に合わせて変更)

set -e
WORK_DIR="$(cd "$(dirname "$0")/.." && pwd)"
cd "$WORK_DIR"
mkdir -p tmp outputs

# [1/4] vc-relax
mpirun -np 16 $QE_BIN/pw.x -in inputs/relax_production.in > tmp/relax.out

# [2/4] scf
mpirun -np 16 $QE_BIN/pw.x -in inputs/scf_after_vc.in > tmp/scf.out

# [3/4] ph.x
mpirun -np 16 $QE_BIN/ph.x -in inputs/ph_ir_production.in > tmp/ph.out

# [4/4] dynmat.x
$QE_BIN/dynmat.x < inputs/dynmat.in > tmp/dynmat.out

プロットスクリプト(plot_spectrum.py)

詳細は、ダウンロードで確認することができます。

# IRスペクトル比較プロット(ONCV PP vs USPP)
#
# 処理内容:
# 1. Hale & Querry (1973)実験データ読み込み
#    - ファイル: download/hale73_normalized.txt
#    - データ: 波数[cm^-1]と正規化吸収強度
#
# 2. dynmat.outから振動数・IR強度読み込み
#    - ONCV PP: tmp/dynmat.out(自ディレクトリ)
#    - USPP:    ../uspp_no_broad_USPPPP/tmp/dynmat.out
#
# 3. ガウスブロードニングで可視化
#    - シグマ: σ = 50 cm⁻¹
#    - 目的: ピーク位置比較のための可視化
#
# 4. IR強度加重平均計算
#    - 式: Σ(freq × IR) / Σ(IR)
#    - 範囲: 実験ピーク ±500 cm⁻¹
#    - ONCV PP: 3343.6 cm⁻¹(15モード)
#    - USPP:    3474.5 cm⁻¹(15モード)
#
# 出力: outputs/ir_spectrum_comparison.png

import numpy as np
import matplotlib.pyplot as plt

# ガウス関数定義
def gaussian(x, mu, sigma, intensity):
    return intensity * np.exp(-0.5 * ((x - mu) / sigma)**2)

# 実験データ読み込み
# Hale & Querry (1973)正規化データ
hale_data = np.loadtxt('download/hale73_normalized.txt')

# dynmat.out読み込み(振動数・IR強度)
oncv_data = load_dynmat('tmp/dynmat.out')
uspp_data = load_dynmat('../uspp_no_broad_USPPPP/tmp/dynmat.out')

# ガウスブロードニングでスペクトル生成
sigma = 50.0  # cm^-1
spectrum = sum(gaussian(x, freq, sigma, ir) for freq, ir in data)

# プロット作成
plt.plot(wavenumber, exp_intensity, label='Exp. Hale & Querry (1973)')
plt.plot(wavenumber, oncv_spectrum, label='Harmonic ONCV PP (80 Ry)')
plt.plot(wavenumber, uspp_spectrum, label='Harmonic USPP (80 Ry)')

結果

IRスペクトル比較

IRスペクトル比較

図の説明: 実験データと調和近似計算の比較(ecutwfc=80 Ryで統一)。

  • 濃い灰色実線: Hale & Querry (1973)実験値
  • 参考縦線: 主要ピーク位置(680/1645/3390 cm⁻¹)
  • 青実線: 調和近似 USPP (80 Ry)
  • 赤実線: 調和近似 ONCV PP (80 Ry)

: USPPはUltrasoft Pseudopotentialの略、ONCV PPはOptimized Norm-Conserving Vanderbilt Pseudopotentialの略

ピーク位置の比較

O-H伸縮モード(±500 cm⁻¹範囲のIR強度加重平均)

IR強度加重平均 = Σ(freq × IR) / Σ(IR) で計算(面積中心)

  • ONCV PP: ΣIR=55.93, Σ(freq×IR)=187,067 → 3343.6 cm⁻¹ (15モード)
  • USPP: ΣIR=50.22, Σ(freq×IR)=174,507 → 3474.5 cm⁻¹ (15モード)
方法 O-H伸縮 [cm⁻¹] 実験値(3389.8)との誤差
USPP (80 Ry) 3474.5 +2.5%
ONCV PP (80 Ry) 3343.6 -1.4%
実験値(Hale & Querry 1973) 3389.8 基準

多ピーク構造の比較

ONCV PPの計算結果は、実験で観測されるピーク分布の傾向を部分的に再現しています。

領域 Hale & Querry (1973) ONCV PP計算 比較
Librational ~680 cm⁻¹ 686, 694, 716, 750 cm⁻¹ 位置傾向は一致
Bending ~1645 cm⁻¹ 1673, 1680, 1745 cm⁻¹ 位置傾向は一致
Stretching 3389.8 cm⁻¹(最大ピーク) 3064-3529 cm⁻¹(多峰性) 位置傾向は一致

注意: 計算の多峰性は8分子スーパーセル内の異なる水素結合環境によるもので、実験の温度・ダイナミクス効果による広がり(数百cm⁻¹)とは質的に異なります。

擬似ポテンシャルの効果

ONCV PPはUSPPよりも低波数側にシフトし、実験値に近づきます。

USPP:  3474.5 cm⁻¹  →  +2.5%誤差
ONCV:  3343.6 cm⁻¹  →  -1.4%誤差

この差(130.9 cm⁻¹)は擬似ポテンシャルの質の差を反映しています。ONCV PPはnorm-conserving特性により、USPPより低波数シフトしやすい傾向が見られます。これはONCV PPのオール電子近似への忠実性が高いことに起因します。

主要ピーク比較表(実験値 vs 計算値)

Hale & Querry (1973)の主要ピーク(吸収係数最大値位置)と、計算値の**±500 cm⁻¹範囲のIR強度加重平均(面積中心)**を比較します。

領域 Hale & Querry (1973) USPP (80 Ry) 誤差 ONCV PP (80 Ry) 誤差
Librational 680 cm⁻¹ (180-1180) 506.5 cm⁻¹ 25.5% 531.9 cm⁻¹ 21.8%
Bending 1645 cm⁻¹ (1145-2145) 1640.5 cm⁻¹ 0.3% 1667.1 cm⁻¹ 1.3%
O-H Stretching 3389.8 cm⁻¹ (2890-3890) 3474.5 cm⁻¹ 2.5% 3343.6 cm⁻¹ 1.4%

:

  • Haleデータ: 吸収係数最大値位置(OMLCデータより、2950 nm = 3389.8 cm⁻¹)
  • 計算値: ±500 cm⁻¹範囲のIR強度加重平均(面積中心)
  • 太字は最も良い一致を示す
  • ONCV PPはO-H伸縮で最良(1.4%誤差)
  • Librationalは範囲が広いため誤差が大きいが、Bending/O-H伸縮は2%以内の精度

考察

1. 振動数の一致:調和近似でも十分な精度

ONCV PP (ecutwfc=80 Ry) で、O-H伸縮振動数の誤差は-1.4% にとどまりました(Hale & Querry 1973の最大ピーク3389.8 cm⁻¹に対し)。

対してUSPPは+2.5%の誤差で、擬似ポテンシャルの選択が重要であることが分かります。

これは調和近似の範囲で十分に良い精度と言えます。液相水のような複雑な系(水素結合ネットワーク)で、ONCV PPで2%以内の誤差は優秀な結果です。

2. 吸光度(相対強度)の不一致:調和近似の限界

スペクトルを見ると、500 cm⁻¹と3400 cm⁻¹の相対強度が、計算と実験で逆転しています。

原因:調和近似のIR強度計算の限界

要因 調和近似(計算) 実験(液相)
温度 0 K 300 K
双極子モーメント 平衡構造で1回計算 熱振動で変動
モード間カップリング なし あり
非調和性 なし あり

IR強度は (∂μ/∂Q)² に比例します(μ: 双極子モーメント、Q: 基準座標)。

  • O-H伸縮(3400 cm⁻¹): 水素結合の非調和性が大きい → 調和近似で過大評価
  • Librational(500 cm⁻¹): 低周波数ほど熱振動の影響大 → 調和近似で過小評価

結論

調和近似は「振動数(ピーク位置)」には良いが、「相対強度」には限界がある。

これは調和近似の本質的な限界であり、解決には**非調和計算(SSCHA、AIMD+DACPなど)**が必要です。ただし、これらは計算コストが非常に高く(数日~1週間以上)、実用的ではありません。

3. 実験データの解釈

液相水のO-H伸縮モードは、多様な水素結合環境により、複数のピーク成分から成ります。

多峰性の解釈

O-H伸縮バンドは、水素結合の強さによって以下のように分解できます:

  • ~3220 cm⁻¹: 強い水素結合(低波数シフト)
  • ~3400 cm⁻¹: 中程度の水素結合(主ピーク)
  • ~3550 cm⁻¹: 弱い水素結合(高波数)

これは液相水が単一のピークではなく、複数の水素結合状態の重ね合わせであることを示しています。図の濃い灰色実線はこの多峰性を示しています。

主ピーク位置

文献 主ピーク [cm⁻¹] 特徴
Hale & Querry (1973) 3389.8 OMLCデータ、O-H伸縮最大ピーク

本記事ではHale & Querry (1973)の最大ピーク(3389.8 cm⁻¹)を基準として比較しました。

4. ecutwfcの効果(補足)

ecutwfcを60 Ry → 80 Ryに上げると、USPPでは振動数が高波数側にシフトします。

USPP 60 Ry:  3551 cm⁻¹
USPP 80 Ry:  3474.5 cm⁻¹ (±500範囲加重平均)
ONCV 80 Ry:  3343.6 cm⁻¹ (±500範囲加重平均)

これは基底関数の完全性が向上したためです。ただし、USPPとONCVではシフトの方向が異なり、擬似ポテンシャルの特性も影響しています。

結論

主要な発見

  1. ONCV PP (ecutwfc=80 Ry) で実験値に非常に近い結果

    • O-H伸縮振動数: 3343.6 cm⁻¹(Hale & Querry 1973の最大ピーク3389.8 cm⁻¹に対し-1.4%誤差)
    • 多ピーク構造をよく再現
    • USPP (3474.5 cm⁻¹, +2.5%) よりも改善
    • : 基準値を3400 cm⁻¹近辺に置いた場合、ONCV PPの誤差は約-1.6%、USPP PPは約+2.2%となり、本記事の結論(ONCV PPの方が実験値に近い)は変わりません
  2. 調和近似の振動数予測能力は優秀

    • 液相水のような複雑な系で2%以内の誤差
    • 実用的な精度
  3. 調和近似の相対強度には限界

    • 500 cm⁻¹と3400 cm⁻¹の強度が逆転
    • 非調和性・温度効果が必要
    • 完全な解決には高コストの計算(AIMD+DACP)が必要

実用的な推奨

  • 振動数(ピーク位置)を評価軸にする: 調和近似で十分
  • 相対強度は参考程度: 調和近似の限界として文献引用で説明
  • 擬似ポテンシャル選択: ONCV PP推奨(オール電子近似に忠実)

今後の展望

相対強度の改善には以下の方法が考えられますが、コスト対効果は低いです:

方法 計算コスト 期待される改善
SSCHA 数時間~1日 振動数±1%、強度△
AIMD + DACP 数日~1週間 強度○(±20-30%誤差残る可能性)

補足: 本検討ではAIMD+VACF(速度自己相関関数)も試行しましたが、低周波数領域(124 cm⁻¹)の異常ピークが支配的で、O-H伸縮の強度が過小評価され、実用的ではありませんでした。DACP(双極子相関関数)が正しいアプローチですが、計算コストが非常に高いため、現状では調和近似が最も実用的です。

現時点では、ONCV PP + 調和近似が最も実用的なアプローチと言えます。

参考文献

実験データ

  • Hale, G. M., & Querry, M. R. (1973). Optical constants of water in the 200nm to 200 micron wavelength region. Applied Optics, 12, 555--563. DOI: 10.1364/AO.12.000555

計算手法

  • Giannozzi, P., et al. (2009). QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. Journal of Physics: Condensed Matter, 21(39), 395502.
  • Wilson, E. B., Decius, J. C., & Cross, P. C. (1955). Molecular Vibrations. McGraw-Hill.

擬似ポテンシャル

  • ONCV PP: Optimized norm-conserving Vanderbilt pseudopotentials
  • USPP: Ultrasoft pseudopotentials

ダウンロード

擬似ポテンシャルの入手方法

本計算で使用した擬似ポテンシャルファイルは以下から入手できます。

注意: IR計算(epsil=.true.)ではNorm-Conserving PP(HGH/ONCV)が最も安定です。

擬似ポテンシャルの入手詳細は、以前のH2O(気相)を取り扱った記事に載っています。

実験データの準備

plot_spectrum.pyで使用する実験データを以下の手順で準備します。

  1. OMLCからデータを入手

    https://omlc.org/spectra/water/data/hale73.txt
    
  2. データ変換

    • 波長[nm] → 波数[cm^-1]: wavenumber = 10^7 / wavelength
    • 0-1に正規化: normalized = absorption / max(absorption)
    • 50-4000 cm^-1範囲を抽出
  3. ファイル保存

    • ファイル名: hale73_normalized.txt
    • 配置先: materials/h2o_liquid/download/hale73_normalized.txt
    • 形式: 波数[cm^-1] 正規化吸収強度(2列、スペース区切り)

計算ファイル一式

本計算の入力ファイル・スクリプト・実験データは以下のGistから入手できます。
Gist 液相H2O 計算ファイル一式

Discussion