付録B計算環境の準備とPySCF入門

B.1 まず覚える5つの言葉

B.2 ターミナルとUNIXコマンド

macOSなら「ターミナル」アプリ，WindowsならWSL（Windows Subsystem for Linux）上のUbuntuを起動する．黒い画面に文字を打ち込む，あの画面である．

B.2.1 移動と一覧

$ pwd                       # いまいる場所を表示 (print working directory)
/Users/ymochizuki

$ ls                        # この場所にあるファイルを一覧 (list segments)
Desktop  Documents  Downloads

$ ls -a                     # 隠しファイルも含めて全部表示
$ ls -l                     # 詳細（サイズ・更新日時）つきで表示

$ mkdir 2026_comp_mat_sci   # ディレクトリを作る (make directory)
$ cd 2026_comp_mat_sci      # 移動する (change directory)
$ pwd
/Users/ymochizuki/2026_comp_mat_sci

$ cd ..                     # 1つ上の階層へ戻る
$ cd ~                      # ホームディレクトリへ戻る（cd だけでも同じ）

コードB.1 移動と一覧の基本コマンド．$ はプロンプトなので入力しない．

B.2.2 コピー・削除・表示

$ cp XYZ_H.xyz XYZ_H_backup.xyz    # ファイルをコピー (copy)
$ cp ~/Downloads/XYZ_CH4.xyz .     # 別の場所から今いる場所へコピー
$ cp -r pyscf_tutorial1 backup/    # ディレクトリごとコピー (recursive)

$ mv old.xyz new.xyz               # 名前の変更／移動 (move)

$ rm XYZ_H_backup.xyz              # ファイルを削除 (remove)
$ rm -r backup                     # ディレクトリごと削除

$ cat XYZ_H2O.xyz                  # ファイルの中身を表示 (concatenate)
3
962
  O   6.30121   0.00000  -0.38499
  H   7.63502   0.00000   0.38499
  H   4.96764   0.00000   0.38499

$ head -20 output.txt              # 先頭20行だけ表示
$ tail -20 output.txt              # 末尾20行だけ表示

コードB.2 ファイル操作の基本コマンド．

B.2.3 文字列の置換・検索・切り出し

第8章で $\mathrm{O_2}$ の結合長を系統的に変えるときに使う3つのコマンドである．

# sed: 文字列を置換して別ファイルに保存する
$ sed -e 's/1.20000/2.00000/g' XYZ_O2_1p20000.xyz > XYZ_O2_2p00000.xyz

# grep: 特定の文字列を含む行だけ抜き出す
$ grep "Total Energy" output_1p20000.txt
  Total Energy (eV) : -4069.348

# 複数ファイルをまとめて検索（ワイルドカード）
$ grep "Total Energy" output_*.txt

# awk: 空白区切りの n 列目だけ取り出す
$ grep "Total Energy" output_*.txt | awk '{print $NF}'
-4069.348
-4061.466
...

コードB.3 sed／grep／awk．縦棒 | は「パイプ」で，左のコマンドの出力を右のコマンドの入力に渡す．

awk の $1, $2, … は空白区切りの列番号，$NF は最後の列を意味する．出力の書式に応じて列番号を数え，必要な数字だけを取り出す．

B.3 Python環境を用意する

Pythonを動かす方法は大きく2つある．

B.3.1 ローカル環境

「windows anaconda install」「mac anaconda install」などで検索し，案内に従ってインストールする．インストール後，ターミナルで次を確認する．

$ which python3          # python本体がどこにあるか
/Users/ymochizuki/anaconda3/bin/python3

$ python3 --version
Python 3.11.5

$ echo $SHELL            # シェルの種類
/bin/zsh

$ pip3 list | grep pyscf # PySCFが入っているか

コードB.4 環境の確認．which が何も返さなければ，まだPythonにパスが通っていない．

B.3.2 必要なモジュールのインストール

$ pip3 install pyscf
$ pip3 install numpy matplotlib
$ pip3 install pymatgen        # 第13章で使う
$ pip3 install geometric       # 構造最適化に使われることがある

コードB.5 モジュールのインストール．

Jupyter Notebook（拡張子 .ipynb）を開くには，Visual Studio Code をインストールしてPython拡張機能を入れるのが手軽である．講義配布の pyscf_tutorial1.ipynb などはこの形式になっている．

B.4 Google Colabを使う

環境構築で詰まったら，迷わずGoogle Colabに切り替えてよい．https://colab.research.google.com/ にアクセスし，「ノートブックを新規作成」をクリックするだけで，ブラウザ上にJupyter Notebookが立ち上がる．ファイル名は 2026_マテリアル計算科学.ipynb のように分かりやすく変えておくとよい．

# Colabのセルに書く内容
!pwd
!ls
!pip install pyscf
%cd /content/drive/MyDrive/pyscf_tutorial1
!python3 calc_pyscf.py --xyz XYZ_H.xyz --spin 1

B.4.1 Google Driveと接続する

Colabの作業領域（/content）は，12時間以上放置すると中身がすべて消える．作ったファイルを残すには，自分のGoogle Driveをマウント（接続）しておく．

# Colabのセルで実行
from google.colab import drive
drive.mount('/content/drive')

コードB.6 Google Driveのマウント．認証画面が出るので許可する．

マウントすると /content/drive/MyDrive/ が自分のDriveのルートになる．講義配布の pyscf_tutorial1.zip は必ず解凍してからDriveに置くこと（zipのまま置くと中のファイルにアクセスできない）．

# 12時間放置して環境がリセットされたら，これを順に実行し直すだけでよい
from google.colab import drive
drive.mount('/content/drive')

!pip install pyscf

%cd /content/drive/MyDrive/pyscf_tutorial1

!python3 calc_pyscf.py --xyz XYZ_Be.xyz --spin 0 --basis 6-31g

コードB.7 セッションが切れたあとの復帰手順．

使い終わったらマウントを外しておくとよい．

drive.flush_and_unmount()

コードB.8 マウントの解除．

B.5 PySCFのインストールと最初の計算

準備ができたら，講義配布の pyscf_tutorial1 ディレクトリに移動して，水素原子を計算してみよう．

$ cd ~/Downloads/pyscf_tutorial1
$ ls
XYZ_B.xyz  XYZ_Be.xyz  XYZ_CH4.xyz  XYZ_F.xyz  XYZ_H.xyz  XYZ_He.xyz
XYZ_Li.xyz XYZ_N.xyz   XYZ_Ne.xyz   XYZ_O.xyz  XYZ_O2.xyz  calc_pyscf.py

$ python3 calc_pyscf.py --xyz XYZ_H.xyz --spin 1

コードB.9 水素原子の計算．--spin 1 は「アップスピンとダウンスピンの数の差が1」という意味．

出力には，計算手法（Method : UHF），全エネルギー（Total Energy），分子軌道エネルギー（MO energies）と占有数が並ぶ．STO-3G基底での水素原子の全エネルギーは $-12.696$ eV で，厳密解 $-13.606$ eV に対して 6.7% の誤差である（第6章6.7節）．

B.6 構造ファイルの書式

B.6.1 xyzファイル（分子）

分子の構造を書く最も単純な形式である．

3
962
  O        6.30121      0.00000       -0.38499
  H        7.63502      0.00000        0.38499
  H        4.96764      0.00000        0.38499

コードB.10 XYZ_H2O.xyz．1行目が原子数，2行目はコメント（何を書いてもよい），3行目以降が「元素記号 x y z」．座標の単位はÅ．

B.6.2 SDFファイルからの変換

PubChem（https://pubchem.ncbi.nlm.nih.gov/）で分子名や組成式を検索し，「3D Conformer」の Download Coordinates から SDF 形式で保存できる．ファイル名は SDF_D-glucose.sdf のように「何の分子か」が分かる形に変えておくと，あとで ls したときに一目で分かる．

$ python3 sdf2xyz.py --sdf SDF_D-glucose.sdf --xyz XYZ_D-glucose.xyz
$ head -3 XYZ_D-glucose.xyz
24
converted from SDF_D-glucose.sdf
  C   -0.7810    1.2680    0.4980

コードB.11 SDF → xyz の変換．--xyz を省略すると既定の名前で出力される．

B.6.3 POSCARファイル（結晶）

第13章のpymatgenでは，結晶構造をVASPのPOSCAR形式で与える．Materials Project（https://next-gen.materialsproject.org/）から無償でダウンロードできる．

Na4 Cl4
1.0
   5.6920    0.0000    0.0000
   0.0000    5.6920    0.0000
   0.0000    0.0000    5.6920
Na Cl
4 4
Direct
  0.000  0.000  0.000
  0.000  0.500  0.500
  ...

コードB.12 POSCARの構造．1行目はコメント，2行目はスケール因子，3〜5行目が格子ベクトル，6行目が元素，7行目が各元素の原子数，8行目が座標系（Direct＝分率座標），以降が原子座標．

B.7 計算スクリプト一覧とオプション

講義で配布されるスクリプトは6つある．いずれも --help を付けると使い方が表示される．

$ python3 calc_pyscf_dos.py --help

コードB.13 困ったらまず --help．

B.7.1 共通のオプション

B.7.2 calc_pyscf.py の固有オプション

# 第9章：全エネルギーの内訳を見る
$ python3 calc_pyscf.py --xyz XYZ_He.xyz --spin 0 --basis 6-31g --decompose-total-energy

# 第6章：Li 原子を UHF と ROHF で比べる
$ python3 calc_pyscf.py --xyz XYZ_Li.xyz --spin 1
$ python3 calc_pyscf.py --xyz XYZ_Li.xyz --spin 1 --method rohf

# 第6章：電荷を変える（Li+）
$ python3 calc_pyscf.py --xyz XYZ_Li.xyz --spin 0 --charge 1

# 第7章：炭素原子の低スピンと高スピン
$ python3 calc_pyscf.py --xyz XYZ_C.xyz --spin 2 --basis 6-31g
$ python3 calc_pyscf.py --xyz XYZ_C.xyz --spin 4 --basis 6-31g

コードB.14 calc_pyscf.py の使用例．

B.7.3 構造最適化・DOS・波動関数・振動

# 第8章・第10章：構造緩和．XYZ_○○-finish.xyz が出力される
$ python3 calc_pyscf_relax.py --xyz XYZ_O2_0p900000.xyz --spin 2 --basis 6-31g

# 第10章：状態密度（元素・軌道分解つき，表示範囲を -30〜10 eV に）
$ python3 calc_pyscf_dos.py --xyz XYZ_D-glucose-finish.xyz --spin 0 \
      --basis 6-31g --xrange -30 10 --element-pdos

# 第10章：COHP（結合の強さの解析）
$ python3 calc_pyscf_dos.py --xyz XYZ_CH4-finish.xyz --spin 0 --basis 6-31g \
      --cohp --pair C,H

# 第10章：HOMO と LUMO を cube ファイルに出力
$ python3 calc_pyscf_wf.py --xyz XYZ_CH4-finish.xyz --spin 0 --basis 6-31g \
      --homo --lumo

# 特定の分子軌道を指定して出力（0始まりの番号）
$ python3 calc_pyscf_wf.py --xyz XYZ_CH4-finish.xyz --spin 0 --basis 6-31g \
      --mo 1 2 3 4 7

# 第10章：振動数の計算
$ python3 calc_pyscf_vib.py --xyz XYZ_H2O-finish.xyz --spin 0 --basis 6-31g

コードB.15 各スクリプトの代表的な使い方．行末の \ は「次の行に続く」という意味で，1行で書いてもよい．

calc_pyscf_dos.py には --sigma（Gauss広がりの幅，既定 0.3 eV），--csv／--txt／--output（数値データや図の保存先）といったオプションもある．calc_pyscf_wf.py は --nx --ny --nz でcubeファイルの格子点数，--margin で箱の余白（Bohr単位）を調整できる．

B.8 forループによる系統計算と結果の抽出

第8章で行う「$\mathrm{O_2}$ の結合長を変えて全エネルギーを計算する」という作業を，最初から最後まで通してみよう．手作業なら数十回の繰り返しになるが，ターミナルなら数行で済む．

B.8.1 構造ファイルを量産する

# 元になるファイル（結合長 1.20000 Å）
$ cat XYZ_O2_1p20000.xyz
2
O2 molecule
  O   0.00000   0.00000   0.00000
  O   0.00000   0.00000   1.20000

# 0.90 Å から 1.90 Å まで 0.10 Å 刻みで構造ファイルを作る
$ for d in 0p90000 1p00000 1p10000 1p20000 1p30000 1p40000 1p50000 \
           1p60000 1p70000 1p80000 1p90000
  do
    r=$(echo $d | sed -e 's/p/./')
    sed -e "s/1.20000/$r/" XYZ_O2_1p20000.xyz > XYZ_O2_${d}.xyz
  done

$ ls XYZ_O2_*.xyz
XYZ_O2_0p90000.xyz  XYZ_O2_1p00000.xyz  ...  XYZ_O2_1p90000.xyz

コードB.16 sed と for ループによる構造ファイルの量産．ドル記号は変数を意味する．ファイル名にピリオドを使うと紛らわしいので p で代用している．

B.8.2 まとめて計算する

$ for f in XYZ_O2_*.xyz
  do
    tag=$(basename $f .xyz)
    python3 calc_pyscf.py --xyz $f --spin 2 --basis 6-31g > output_${tag}.txt
  done

コードB.17 全ファイルを順に計算し，結果をテキストファイルに保存する．> は「出力をファイルに書き出す」という意味（リダイレクト）．

B.8.3 結果を抜き出してプロットする

# 全エネルギーの行だけ抜き出す
$ grep "Total Energy" output_*.txt

# 数値だけ取り出す
$ grep "Total Energy" output_*.txt | awk '{print $NF}'

# 1列目に結合長，2列目に全エネルギーを並べたファイルを作る
$ for f in output_XYZ_O2_*.txt
  do
    d=$(echo $f | sed -e 's/.*O2_//' -e 's/.txt//' -e 's/p/./')
    e=$(grep "Total Energy" $f | awk '{print $NF}')
    echo "$d $e"
  done > energy_vs_distance.dat

$ cat energy_vs_distance.dat
0.90000 -4060.401
1.00000 -4066.xxx
...
1.20000 -4069.348
...
2.00000 -4061.466

コードB.18 結果を2列のデータファイルにまとめる．

# plot_energy.py
import numpy as np
import matplotlib.pyplot as plt

d, e = np.loadtxt("energy_vs_distance.dat", unpack=True)

plt.figure(figsize=(5, 4))
plt.plot(d, e, "o-", color="#14684a")
plt.xlabel(r"Bond length $R$ (Å)", fontsize=13)
plt.ylabel("Total energy (eV)", fontsize=13)
plt.grid(alpha=0.4)
plt.tight_layout()
plt.savefig("O2_potential_curve.png", dpi=200)
plt.show()

コードB.19 matplotlibによるプロット．教科書に出てくる原子間ポテンシャル曲線がそのまま得られる．

B.9 VESTAによる可視化

VESTA（https://jp-minerals.org/vesta/）は，結晶構造・分子構造・波動関数を可視化する無償ソフトウェアである．日本で開発され，世界中の第一原理計算ユーザーが使っている．

使い方の要点をまとめる．

• ドラッグで回転，ホイールで拡大縮小．
• 2つの原子を続けてクリックすると結合長，3つ続けてクリックすると結合角が下部に表示される．実験値と比べるときに使う．
• cubeファイルを開くと等値面が表示される．黄色が波動関数の正の領域，青色が負の領域である．等値面の値は Properties → Isosurfaces で変えられる．
• 正と負の領域がどう並んでいるかを見れば，結合性軌道か反結合性軌道かが判定できる．原子間で同符号（節がない）なら結合性，異符号（節がある）なら反結合性である（第11章）．

B.10 pymatgenとMadelungエネルギー

第13章では，イオン結晶の静電エネルギー（Madelungエネルギー）をpymatgenで計算する．pymatgenは，米国が2011年に始めた Materials Genome Initiative の中で開発されたPythonモジュールで，Materials Informatics（材料インフォマティクス）の事実上の標準になっている．

$ pip3 install pymatgen
$ cd ~/Downloads/pymatgen_tutorial
$ ls
POSCAR_NaCl_Fm-3m  POSCAR_MgO_Fm-3m  POSCAR_TiO2_P42mnm_rutile  ...
calc_madelung_constant.py

$ python3 calc_madelung_constant.py POSCAR_NaCl_Fm-3m

コードB.20 NaClのMadelungエネルギーの計算．

出力には，セルあたり（Na 4個・Cl 4個）のMadelungエネルギーと，化学式単位（formula unit，Na 1個・Cl 1個）あたりのMadelungエネルギー，そしてMadelung定数が並ぶ．「セルあたり」と「化学式単位あたり」を取り違えると4倍ずれるので注意してほしい．

イオンの価数が自動推定でうまくいかないときは，明示的に指定できる．

# 価数を明示する
$ python3 calc_madelung_constant.py POSCAR_Mg2Si_Fm-3m --oxidation-state Mg=2 Si=-4

# 特定のカチオン・アニオン対を基準にMadelung定数を出す
$ python3 calc_madelung_constant.py POSCAR_NaCl_Fm-3m --reference-pair Na Cl

コードB.21 価数の指定と参照イオン対の指定．

スクリプトの中身は本質的に数行である．

from pymatgen.analysis.ewald import EwaldSummation
from pymatgen.io.vasp.inputs import Poscar

structure = Poscar.from_file("POSCAR_NaCl_Fm-3m").structure
structure.add_oxidation_state_by_guess()

ewald = EwaldSummation(structure)
print("Madelung energy (cell) :", ewald.total_energy, "eV")

コードB.22 pymatgenの中核部分．EwaldSummation がEwald法による静電エネルギーの計算を担う．

B.11 よくあるエラーと対処

B.12 演習問題

参考文献