3.2. 計算条件と再現手順

このページでは、水素原子を並べると電子バンドは生まれるのかの数値と図を得た計算の条件、スクリプト、数値データ、および再実行の手順をまとめます。

3.2.1. 計算環境

計算は2026年10月に、WSL上のDockerで行いました。使用したのは、Psi4 1.10を含むRadonPy用のローカルDockerイメージ（ radonpy-local:v1.0-beta.1 ）です。Quloud上では実行していません。このイメージは公開していないため、以下ではPsi4 1.10、NumPy、Matplotlibが使えるPython環境で各スクリプトを実行する形で手順を示します。他のPsi4のバージョンでの再実行は確認していません。

全計算に共通する条件は次のとおりです。

中性・直線・等間隔・開放端の水素鎖。原子は \(x\) 軸上に1.8 Å間隔で配置（ symmetry c1 、 no_reorient 、 no_com ）
基底関数STO-3G、電荷0、スピン多重度1（ \(N_\alpha = N_\beta = N/2\) ）
SCF収束条件はエネルギー・密度ともに \(10^{-9}\) 、最大500反復
PBEの積分グリッドはPsi4の既定値（動径75点、角度302点）
スレッド数2

3.2.2. 配布ファイル

3.2.2.1. Psi4入力スクリプト

hydrogen_chain_psi4.py : 原子数を変えながらRHFまたはRKSを計算し、軌道エネルギーと擬似波数をCSVに出力します（DF-SCF）。
hydrogen_chain_bond_order.py : RHFとRKS-PBEのLöwdin結合次数を求め、RHF解の安定性解析を行います（PK積分）。
hydrogen_chain_spin.py : RHF、RKS-PBE、UHF、UKS-PBEを順に計算し、エネルギー、軌道ギャップ、結合次数、スピン密度、 \(\langle S^2\rangle\) を出力します（PK積分）。

3.2.2.2. 図の作成スクリプト

plot_article_figures.py : 下記のCSVからモデル図、軌道エネルギー、ギャップ、スピンと結合次数の比較図を作成します。Psi4の計算は行いません。記事中のその他の図は配布済みのPNGを使用しています。

3.2.2.3. 数値データ（CSV）

ファイル	内容
`rhf_orbitals.csv`	RHFの軌道エネルギー、擬似波数、占有数（ \(N\) = 2–64、DF）
`pbe_core_orbitals.csv`	RKS-PBEの同上（ \(N\) = 2–32、コアハミルトニアン初期推定、DF）
`pbe_sad_orbitals.csv`	RKS-PBEの同上（ \(N\) = 32, 48, 64、SAD初期推定、DF）
`rhf_pbe_comparison.csv`	上の3ファイルから抜き出したHOMO、LUMO、軌道ギャップ
`hydrogen_chain_bond_order.csv`	RHFとRKS-PBEの隣接原子間Löwdin結合次数（ \(N\) = 16, 32, 64、PK）。原子番号は0始まり
`hydrogen_chain_spin_summary.csv`	4手法の全エネルギー、軌道ギャップ、中央の交替振幅、スピン密度、 \(\langle S^2\rangle\) （ \(N\) = 16, 32、PK）
`hydrogen_chain_spin_sites.csv`	4手法の原子ごとのLöwdin電子数とスピン密度。原子番号は1始まり
`hydrogen_chain_spin_bonds.csv`	4手法の隣接原子間Löwdin結合次数。原子番号は1始まり

軌道エネルギーの単位はeV、全エネルギーの単位はHartree（ \(E_\mathrm{h}\) ）です。

3.2.3. 計算手順

以下のコマンドは、Psi4 1.10の環境で、スクリプトを置いたディレクトリから実行します。各スクリプトはPsi4のログ（ *_psi4.log ）も出力します。

3.2.3.1. 1. RHFの原子数依存性

python hydrogen_chain_psi4.py --method rhf --a 1.8 --ns 2 4 8 16 24 32 48 64 80 100 --prefix rhf

rhf_orbitals.csv が得られます。CSVには収束した系だけが書き出されます。この条件では \(N\) = 80, 100が500反復以内に収束せず、結果から除外されました。

3.2.3.2. 2. RKS-PBEの原子数依存性

python hydrogen_chain_psi4.py --method rks --xc pbe --a 1.8 --ns 2 4 8 16 24 32 --prefix pbe_core

\(N\) = 48, 64では、初期推定をSADに変えて計算しました。 hydrogen_chain_psi4.py の set_options にある 'guess': 'core' を 'guess': 'sad' に書き換えたコピー（ここでは pbe_sad_input.py とします）を使います。

python pbe_sad_input.py --method rks --xc pbe --a 1.8 --ns 32 48 64 --prefix pbe_sad

コアハミルトニアン初期推定では、 \(N\) = 64が500反復以内に収束しませんでした。 \(N\) = 32では2つの初期推定のギャップが表示桁数で一致し、図と表にはコアハミルトニアン初期推定の結果を使っています。

3.2.3.3. 3. 結合次数とRHFの安定性解析

python hydrogen_chain_bond_order.py --a 1.8 --ns 16 32 64

hydrogen_chain_bond_order.csv が得られます。Psi4 1.10の安定性解析は密度フィッティングのSCFに対応していないため、このスクリプトはPK積分を使います。RHF→RHFとRHF→UHFの安定性固有値は、ログの Lowest singlet (RHF->RHF) stability eigenvalues と Lowest triplet (RHF->UHF) stability eigenvalues に出力されます。

3.2.3.4. 4. スピン制限解と非制限解の比較

python hydrogen_chain_spin.py --ns 16 32

hydrogen_chain_spin_summary.csv 、 _sites.csv 、 _bonds.csv 、 _orbitals.csv （α・β軌道のエネルギー）が得られます。

非制限計算の初期軌道は次のようにつくります。まず同じ手法の制限解（RHFまたはRKS-PBE）のFock行列またはKohn–Sham行列をLöwdin直交化し、対角成分に \(\pm 0.5\ E_\mathrm{h}\) の交互の場をαとβで逆符号に加えて対角化します。得られた軌道をPsi4の restart_file として読み込ませます。この場は初期軌道の生成だけに使い、SCF本計算には加えていません。ログに Reading orbitals が出力されていることで、初期軌道の読み込みを確認できます。単純なHOMO–LUMO混合による初期推定では、 \(N\) = 16のUHFが500反復以内に収束しませんでした。

3.2.3.5. 5. 図の作成

python plot_article_figures.py

scripts/ と同じ階層の data/ にあるCSVを読み、 figures/ に上記の図を出力します。あわせて、中央の結合次数の交替振幅を表示します。

3.2.4. 数値を扱う上での注意

原子数依存性（手順1・2）はDF-SCF、結合次数とスピンの比較（手順3・4）はPK積分で計算しています。同じ \(N\) でも、軌道ギャップは末尾の桁（0.003 eV以下）で一致しません。例えば \(N\) = 64のRHFギャップはDFで4.764 eV、PKで4.762 eVです。
中央の交替振幅は、鎖中央の8本の結合を1本おきの2組に分け、2組の平均値の差の絶対値を2で割った値です。
非制限計算の軌道ギャップは、両スピンを通じた最高被占有軌道と最低空軌道の差です。スピンごとのギャップも hydrogen_chain_spin_summary.csv に記録しています。
\(\langle S^2\rangle\) は密度行列から計算した一行列式に対する値で、Psi4のログの値と照合しました。UKSではKohn–Sham行列式の診断値として扱います。
SCFが収束しても、得られた解が唯一の解や大域的な最低エネルギー解であることは保証されません。
擬似波数 \(ka = j\pi/(N+1)\) はエネルギー順の軌道番号から割り当てたもので、周期系のk点計算や、軌道係数から求めた波数ではありません。