チュートリアル: DFT 構造最適化から物性まで

チュートリアル: DFT 構造最適化から物性まで#

このチュートリアルは Part I を 1 つの現実的な課題でまとめます:分子の大まかな推定から出発し、DFT に平衡構造を見つけさせ、その構造を解析します。これは量子化学研究の日常の形です — 最適化してから解釈する。数値を確かめやすいよう、ここでも水を使います。

最後まで進めば、中核概念のチェックポイント/ステップ/実行モデルを、複数ステップの計算で使ったことになります。

ステップ 1 — （わざと不完全な）分子を作る#

悪い推定から始めます:明らかにずれた O–H 距離と角度。最適化の目的は、それを直すことです。

import qc

mol = qc.chk.new(
    atom="O 0.00 0.00 0.00; H 0.00 0.80 0.55; H 0.00 -0.80 0.55",   # 粗い、平衡ではない
    ao="def2-svp",
    unit="angstrom",
)

ステップ 2 — DFT で構造最適化#

SCF ステップ（B3LYP）を追加し、それに .opt() を付けて実行します。.opt() はエネルギーとその勾配（各原子にかかる力）を繰り返し計算し、力が消えるまで核を坂下へ動かします — それが平衡構造です。

mol = mol.scf(ref="r", xc="b3lyp").opt().run()

print("converged:", mol.opt.converged)          # True
print("energy    :", round(mol.opt.energy, 6))   # -76.358316  (hartree)

.run() は既定では沈黙します;run(log="stdout") を加えると最適化の進捗を流します — サイクルごとのエネルギーとその勾配の表で、数ステップ後に Converged! で終わります。（同じ log= スイッチがあらゆる実行を制御します;log_style="modern"/"orca" で見た目を変え、plot=True で軌跡をインライン描画 — クイックスタートの tip とログと出力参照。）最適化された水は

O–H 結合長 ≈ 0.967 Å（粗い推定の 0.97 前後から、今度は両側で整合）、
H–O–H 角 ≈ 103.1° — 実験値 104.5° に近く、B3LYP/def2-SVP の典型的な結果。

注釈

ワークフロー的には何が起きたか .scf(...).opt() は連結した 2 つの pending ステップを追加し、.run() がそれらを解決 — 各最適化サイクルは、少し動いた構造での SCF です。最後のチェックポイント mol は、最適化された構造とその電子状態を保持し、解析の準備が整っています。

Tip

コアを増やす:nthread= 最適化は SCF を何回も（サイクルごとに 1 回）走らせるので、マシンをより多く使う自然な場面です。qc-rs は既定で CPU コア 1 つを使いますが、run(nthread=8) を加えると各 SCF を 8 コアに分散できます — クイックスタートと同じ共有メモリのスレッド並列で、ノート PC や単一ノードでの手軽な高速化です。（ジョブを複数のマシンに分散する MPI プロセス並列 nmpi= はより高度な話題で、後の並列計算と HPC で扱います; ここでは不要です。）

ステップ 3 — 最適化構造を解析する#

最適化されたチェックポイントに物性を尋ねます。まず Mulliken 原子電荷:

q = qc.prop.chrg.mulliken(mol)
print(q["charges"])       # [-0.2945, 0.1472, 0.1472]

酸素が負、水素が正 — 水のおなじみの極性です。この極性は双極子モーメントで定量化できます:

mp = qc.prop.mpol.molecular(mol)
print(round(mp["dipole_magnitude_debye"], 2))    # 2.00  (debye)

約 2.0 D の双極子は水として妥当です（実験値 ≈ 1.85 D;差は基底関数と手法によるもので、基礎のテーマです）。

全体像#

まとめると、この研究は 7 行です:

import qc

mol = qc.chk.new(
    atom="O 0.00 0.00 0.00; H 0.00 0.80 0.55; H 0.00 -0.80 0.55",
    ao="def2-svp", unit="angstrom",
).scf(ref="r", xc="b3lyp").opt().run()

print("E =", round(mol.opt.energy, 6), "hartree, converged:", mol.opt.converged)
print("charges =", qc.prop.chrg.mulliken(mol)["charges"])
print("dipole  =", round(qc.prop.mpol.molecular(mol)["dipole_magnitude_debye"], 2), "D")

これで「最適化してから解析する」ワークフローが完結です:分子を入れると、平衡構造とその物性が出てきます。

学んだこと、そして次へ#

あなたは今:

分子を作り、手法と基底関数を選び;
ステップを連結し（scf → opt）実行し;
結果を読み（mol.opt.energy）物性を計算しました（mulliken、双極子）。

これが qc-rs の中心的なループです。Hartree–Fock と DFT が実際に何を計算しているのか — そしてなぜ双極子が少し大きく出たのか — を理解するには、Part II — 基礎を読んでください。個々のステップを深掘りするには、ユーザーガイドに各章があり、最後は完全な分子物性スイートです。