分子の入力

分子の入力#

Part III — ユーザーガイドへようこそ。Part I・II でインストールと理論を済ませました;ここからは、実際に使うのと同じく機能ごとに進めます。あらゆる計算は同じ場所から始まります:分子を記述すること。本章はその最初の一歩 — 何を計算するのかを qc-rs にどう伝え、それが正しく理解されたかをどう確認するか — についてです。

理論:量子化学計算にとって分子を定めるものは何か#

Part II で Born–Oppenheimer 近似に出会いました:核は電子よりはるかに重いのでその場に固定し、核が作る場の中で電子問題を解く、というものです。この 1 つの考えが、計算の入力に必要なものを正確に決めます。3 つ、それだけです:

核配置（構造） — どの原子がどこにあるか。核は固定なので、その位置 \(\{\mathbf R_A\}\) はあなたが与える固定値です。これが電子の動く外部ポテンシャル \(v_{\text{ext}}(\mathbf r) = -\sum_A Z_A / |\mathbf r - \mathbf R_A|\) と、電子エネルギーに加わる定数の 核間反発 \(E_{\text{nuc}} = \sum_{A<B} Z_A Z_B / R_{AB}\) を決めます。
全電荷 — 電子が何個あるか。中性分子は \(N_e = \sum_A Z_A\) 個の電子を持ち、電荷 \(q\) はそれを取り除く（\(q>0\)）または加えます（\(q<0\)）。これが SCF の詰める数 \(N_e\) を決めます。
スピン状態 — それらの電子スピンの並び方。多重度 \(2S+1\) で与えます。

それ以外 — 基底関数、手法、収束設定 — はすべて問題をどう解くかの選択です。構造・電荷・スピンは問題そのもの。 qc-rs ではこの 3 つを qc.chk.new(...) に渡し、中核概念のチェックポイントを作ります。

定義 3 (スピン多重度)

全スピン量子数 \(S\) の状態について、多重度は \(2S+1\) — \(M_S\) 射影（\(-S, \dots, +S\)）の数です。不対電子はそれぞれスピン \(\tfrac12\) を与えるので、不対電子 \(n\) 個では \(S = n/2\)、多重度は \(n+1\):

不対電子 \(n\)	\(S\)	多重度 \(2S+1\)	名称
0	0	1	一重項（singlet）
1	½	2	二重項（doublet）
2	1	3	三重項（triplet）

これは Gaussian の慣習で、qc-rs の spin= が期待するものです — 閉殻一重項なら spin=1、ラジカルなら spin=2、三重項なら spin=3 を渡します。

使い方:`atom=` 文字列#

構造は atom= 引数で記述します — 1 行 1 原子の文字列で、各行は元素ラベルに続けて 3 つのデカルト座標:

Element   x   y   z

原子は改行またはセミコロン ; で区切ります。だから次の 2 つは同じ水分子です:

import qc

# セミコロン区切り（コンパクト、ワンライナー向き）
w1 = qc.chk.new(atom="O 0 0 0; H 0 0.757 0.587; H 0 -0.757 0.587", ao="sto-3g")

# 改行区切り（読みやすい、実際の分子向き）— 三重引用符文字列
w2 = qc.chk.new(atom="""
    O   0.000   0.000   0.117
    H   0.000   0.757  -0.469
    H   0.000  -0.757  -0.469
""", ao="sto-3g")

元素ラベルは、化学記号（H, C, Fe — 大文字小文字を区別しない）、原子番号（6 は炭素）、または同一種の原子を区別するラベル接尾辞つき記号（C1, C2, Fe3）が使えます。ラベルは、後で特定の 1 原子に独自の基底関数を与えたいときに効きます — 基底関数と AO 表現参照。

Tip

Python メモ:三重引用符文字列 """… """ は複数行文字列です。実際の構造を書く自然な方法 — 入力ファイルのように 1 行 1 原子で列を揃える — です。先頭/末尾の空行やインデントは問題ありません;qc-rs は無視します。

単位:オングストロームまたはボーア#

座標は既定でオングストローム（unit="angstrom"）です。原子単位で与えるなら unit="bohr" を設定します:

# これらは*同じ* H2 分子（結合長 ≈ 0.741 Å ≈ 1.400 bohr）
h2_ang  = qc.chk.new(atom="H 0 0 0; H 0 0 0.741", ao="sto-3g", unit="angstrom")
h2_bohr = qc.chk.new(atom="H 0 0 0; H 0 0 1.400", ao="sto-3g", unit="bohr")

内部では qc-rs は座標を常にボーア（原子単位）で保持し、入力時に唯一の定数 \(\text{1 Å} = 1/0.52917721092 \approx 1.8897\) bohr で変換します。だから下の coordinates() アクセサは、あなたがどの単位で打っても常にボーアで返します。

重要

既定はオングストローム — ボーアと思い込まない PDB ファイル・論文・多くの GUI からコピーした構造はオングストロームなので、既定でたいてい正解です。しかし原子単位の座標を貼って unit="bohr" を忘れると、すべての結合が約 1.9 倍長くなり、計算は無意味になります（または収束しません）。迷ったら coordinates() で確認してください。

電荷とスピン#

charge= と spin= は既定で中性一重項（charge=0, spin=1）です。イオンや開殻種では変更します。例えば中性の ヒドロキシルラジカル OH は電子 9 個 — 奇数なので一重項にはなれず、二重項（spin=2）です:

oh = qc.chk.new(atom="O 0 0 0; H 0 0 0.97", ao="sto-3g", unit="angstrom",
                charge=0, spin=2)          # 中性の二重項ラジカル
print(oh.nelectron(), oh.spin)             # 9 2

例 1 (電荷とスピンの選び方)

水 H₂O — 中性、全電子対 → charge=0, spin=1。
水酸化物イオン OH⁻ — 電子 1 個追加 → charge=-1, spin=1（10 電子、閉殻）。
メチルラジカル CH₃ — 中性で不対電子 1 個 → charge=0, spin=2。
二酸素 O₂ — 基底状態は三重項（不対電子 2 個）→ charge=0, spin=3。

手早い確認:電子数と多重度のパリティ（偶奇）は整合していなければなりません。偶数の電子数は二重項になれず、奇数は一重項になれません。矛盾すれば run() がその状態を拒否します。

qc-rs が解釈した内容を確認する#

重い計算を走らせる前に、チェックポイントが何を理解したかを見てください。単位ミス・電子数の誤り・打ち間違いを即座に捕まえられます。関連アクセサ（プロパティのものとメソッドのものがある — () に注意）:

w = qc.chk.new(atom="O 0 0 0; H 0 0.757 0.587; H 0 -0.757 0.587", ao="sto-3g", unit="angstrom")

w.natom               # 3          — （実 + ゴースト）原子数
w.symbols             # ['O','H','H']
w.charge, w.spin      # (0, 1)
w.nelectron()         # 10         — SCF が詰める電子数
w.coordinates()       # (3, 3) 配列、常にボーア
w.nuclear_energy()    # 9.188258   — E_nuc（hartree）、電子エネルギーに加わる定数

これが分子指定の全体で、そのまま返ってきます。nelectron() や nuclear_energy() がおかしければ、今入力を直してください — 下流のすべての数値がこれに依存します。

特殊原子#

通常の原子に加え、atom= 文法は 3 種の特殊原子をサポートします。日常作業では不要ですが、いくつかの重要な手法に不可欠です。

ゴースト原子（`Element-Bq`）— カウンターポイズ / BSSE 用#

ゴースト原子は Element-Bq（例 O-Bq, H-Bq）と書き、その元素の基底関数だけを持ち、核も電子も持たない 点です — 核電荷 0、電子 0、核間反発 0。目的は基底関数重ね合わせ誤差（BSSE）のカウンターポイズ補正です: 2 分子が近づくと、各々が相手の基底関数を借りてエネルギーが不自然に下がります。この人工効果を測るには、片方の単量体を相手の基底をゴーストとして置いて再計算します:

# 単量体を「二量体の完全な基底」で計算:相手の原子をゴーストにする
mono = qc.chk.new(atom="O 0 0 0; H 0 0 1.8; H-Bq 0 0 3.5", ao="sto-3g", unit="bohr")
print(mono.natom, mono.nelectron())     # 3 9   — ゴーストは原子として数えるが、電子は加えない

ゴーストは natom に数えられますが（AO 関数を持つため）、nelectron() や nuclear_energy() には何も寄与しません。ゴーストの基底はまずラベルで、次に基底元素で引かれるので、H-Bq は自動的に H の基底を継承します（詳細は基底関数と AO 表現）。

ダミー原子（`X`, `Xx`）— 幾何学的参照点#

ダミー原子（X, Xx、またはラベル付き X1, X-ref）は純粋な幾何マーカーです:核電荷なし、電子なし、 基底関数もなし。物理的なものを何も加えずに参照できる座標（例:対称軸や測定点の定義）です:

d = qc.chk.new(atom="C 0 0 0; H 0 0 1.089; X 0 0 5.0", ao="sto-3g", unit="angstrom")
print(d.natom, d.nelectron())     # 2 7   — ダミーは数えられず、何も加えない

ゴーストと違い、ダミーは natom に数えられません;dummy_atoms() に格納されます。

並進ベクトル（`TV`）— 周期系#

TV で始まる行は、周期系の格子並進ベクトルを与えます（ポリマーは 1 本、シートは 2 本、結晶は 3 本）。 translation_vectors() に格納され、ダミー同様、原子として数えられません:

poly = qc.chk.new(atom="""
    C  -0.574 -0.143  0.376
    C   0.579  0.022 -0.301
    TV  4.848  0.171  0.511
""", ao="sto-3g", unit="angstrom")

発展:原子ごとの核パラメータ#

原子ラベルの後の括弧内に核パラメータを付けられます — 最もよく使うのは同位体（Iso=）とフラグメント番号（Fragment=、カウンターポイズ計算で原子をグループ化するのに使う）です:

# ¹³C 同位体置換体、原子を 2 フラグメントにタグ付け
m = qc.chk.new(atom="C(Iso=13,Fragment=1) 0 0 0; H(Fragment=1) 0 0 1.089",
               ao="sto-3g", unit="angstrom")

全体（Iso, Spin, ZNuc, ZEff, QMom, NMagM, RadNuclear, Fragment, …）、MM 力場フィールド、オプティマイザのフリーズコードは、リファレンス章分子指定に記載しています。これらメタデータの多くは解析・保存されますが、まだハミルトニアンやオプティマイザには影響しません — どのフィールドが現在有効かはリファレンスに注記があります。

注釈

受け付けられないもの qc-rs は分子構造を上記文法の文字列として受け取ります。（他のパッケージが受け付けるような）Python の タプルのリストは、有効な atom= の値ではありません — 文字列形式を使ってください。また名前による自動構造取得はありません:座標は常に自分で与えます（または後で .opt() で最適化します）。

総合例:作って、確認して、SCF へ渡す#

まとめて — 水を記述し、確認し、計算の準備を整えます:

import qc

water = qc.chk.new(
    atom="""
        O   0.000000   0.000000   0.117300
        H   0.000000   0.757200  -0.469200
        H   0.000000  -0.757200  -0.469200
    """,
    ao="cc-pvdz",
    unit="angstrom",
    charge=0,
    spin=1,
)

# --- 何か作業する前の健全性チェック ---
print("atoms     :", water.natom, water.symbols)   # 3 ['O', 'H', 'H']
print("electrons :", water.nelectron())            # 10
print("E_nuc     :", round(water.nuclear_energy(), 6))

# 計算の準備完了 — チェックポイントがこれを以降の全ステップへ運ぶ
water = water.scf(ref="r").run()
print("E(RHF)    :", round(water.scf.energy, 6))   # -76.026772

これで分子は完全に指定され、チェックポイントがそれを後続の各ステップ — 基底関数（次章）、SCF、物性、勾配 — へ運びます。ここで入力を正しく確認しておけば、下流のすべてが盤石の上に立ちます。

練習 1

次の各化学種について qc.chk.new(...) の呼び出しを書き、nelectron() と spin を確認しなさい:

アンモニウムカチオン NH₄⁺（中性 NH₄ なら電子 11 個）。
基底状態の酸素分子 O₂（結合長 1.208 Å）、正しいスピンを含めて。
ヘリウム二量体の基底で計算したヘリウム原子（2 つ目の He は 3 Å 離れたゴースト）。

解答練習 1

# 1. NH4+ : 中性 NH4（11 e-）から電子 1 個を取り除く -> 10 e-、閉殻一重項
nh4 = qc.chk.new(atom="N 0 0 0; H 0.63 0.63 0.63; H -0.63 -0.63 0.63; "
                      "H -0.63 0.63 -0.63; H 0.63 -0.63 -0.63",
                 ao="sto-3g", unit="angstrom", charge=1, spin=1)
# nh4.nelectron() -> 10 ,  nh4.spin -> 1

# 2. O2 の基底状態は三重項（不対電子 2 個）-> spin=3
o2 = qc.chk.new(atom="O 0 0 0; O 0 0 1.208", ao="sto-3g", unit="angstrom",
                charge=0, spin=3)
# o2.nelectron() -> 16 ,  o2.spin -> 3

# 3. He を He 二量体の基底で:相手はゴースト（Bq）
he = qc.chk.new(atom="He 0 0 0; He-Bq 0 0 3.0", ao="cc-pvdz", unit="angstrom")
# he.nelectron() -> 2   （ゴーストは基底関数を加えるが電子は加えない）

分子が定義できたら、次の選択はその軌道をどう展開するか — 基底関数です。それが次章です。