5  入力ファイルを一行ずつ読む

一枚の入力ファイルは、計算への命令であると同時に、ひとつの物理の宣言文だ。
そこに書かれた各行は、どんな仮定を置き、何を出力させているのか。

この章のゴール:Si SCF の入力ファイルを一行ずつ「物理の供述書」として読み、各設定が何を仮定し何を出力させるかを述べられるようになる。

1 入力ファイルは命令書であり、供述書でもある

前章で、全エネルギーはそれ単独では何も語らず、差にだけ意味があると分かった。そして差をとるには、二つの計算が同じ土俵――同じ近似、同じ刻み――に乗っていなければならなかった。では、その「土俵」は、いったいどこに書いてあるのか。

答えは、数十行のテキストファイルの中である。シリコンの全エネルギーを一度求めるために計算機へ手渡すのは、たとえば次のような一枚だ。慣例で si.scf.in と名づけておく。

&control
    calculation = 'scf'
    prefix      = 'si'
    outdir      = './out/'
    pseudo_dir  = './pseudo/'
/
&system
    ibrav     = 2
    celldm(1) = 10.26
    nat       = 2
    ntyp      = 1
    ecutwfc   = 60
    ecutrho   = 240
/
&electrons
    conv_thr    = 1.0d-8
    mixing_beta = 0.7
/
ATOMIC_SPECIES
  Si  28.0855 Si.pbe-n-kjpaw_psl.1.0.0.UPF
ATOMIC_POSITIONS (alat)
  Si  0.00  0.00  0.00
  Si  0.25  0.25  0.25
K_POINTS (automatic)
  8 8 8  0 0 0

このファイルを pw.x という計算プログラムに食わせる(pw.x < si.scf.in のように打つ)と、しばらくして画面に全エネルギーが返ってくる。一見すると、機械に向けた事務的な命令の列だ。だが読み方を変えれば、これは「私はこういうシリコンを、こういう近似と刻みで調べました」という供述書でもある。前章でいう「同じ土俵」とは、煎じ詰めれば、二枚の供述書のこの中身が一致している、ということに他ならない。

そこで本章は、この一枚を一行ずつ読む。各行を命令としてではなく、結果を読むときに添えるべき前提として読む。最後にはあなた自身に、どの行が「どんなシリコンか」を決め、どの行が「どう調べるか」を決めているかを、仕分けてもらう。

2 &control:何を計算し、何を出力させるか

ファイルの頭にある &control から / までが最初のひと区切り――Fortran でいうネームリストだ。ここは「何を計算し、何を出力させるか」を宣言する。

calculation = 'scf' が、いちばん大事な一語である。scf は self-consistent field、つまり「原子は今の配置に固定したまま、その配置での電子状態と全エネルギーだけを解け」という指示だ。原子を動かして安定形を探す relax や、胞の形まで緩める vc-relax もあるが、それは第8章の話。いまは配置を一つ固定し、第1章でいう変換器を一度だけ回す、と宣言している。

ここには書いていないが、力や応力(ストレス)も数字として書き出させたいなら、tprnfor = .true.tstress = .true. の二行を加えればよい。第1章で、変換器が返すのは全エネルギーと各原子にはたらく力だと言った。その力を、内部で使うだけでなく外へ書き出させるのが tprnfor である。tstress は胞全体を縮めよう/広げようとする圧力で、第8章で格子定数を詰めるときの手がかりになる。だがこの実行はエネルギー一つに絞ったので、どちらも付けていない。力そのものは、第8章で必要になったときに改めて呼び出せばよい。

残る prefixoutdirpseudo_dir は、物理ではなく場所の指定――計算結果をどこに書き、擬ポテンシャルをどこから読むか――の事務連絡だ。同じ &control の中でも、物理を語る行と帳簿の行が同居している。一行ずつ読む、とはこの違いに気づくことでもある。

3 &system:結晶と電子の前提を宣言する

&system は、計算する結晶そのものを宣言する区画だ。ここの数行が、第2章までに手に入れた「結晶を数件のデータで置く」作業の、そのままの行き先になる。

最初の二行 ibrav = 2celldm(1) = 10.26 が、結晶の骨格と寸法を決める。ibrav=2 は面心立方(FCC)の格子を選ぶ番号で、第2章で見たとおり、シリコンの結晶はこの FCC の格子に2原子の基底を載せた形をしている。celldm(1) はその格子の寸法で、単位は Å ではなく bohr だ。

ここで検算をしておこう。celldm(1)=10.26 が、本当に私たちの知るシリコンを指しているか。\(1\ \mathrm{bohr}=0.529\ \mathrm{\text{Å}}\) で Å に戻すと

\[ a = 10.26\ \mathrm{bohr} \times 0.529\ \mathrm{\text{Å}/bohr} \approx 5.43\ \mathrm{\text{Å}}. \]

実験で知られるシリコンの格子定数は \(5.431\ \mathrm{\text{Å}}\)。小数第2位までぴたりと合う。残る差は換算係数を丸めたことによるもので、0.1 % に満たない。第2章で \(5.431/0.529 \approx 10.26\) と換算したあの数が、いま入力の中で私たちを出迎えている。ただし、この celldm(1) は計算の結果ではなく、私たちがあらかじめ置いた前提だという点は外せない。計算機がこの格子を自分でも選ぶかどうかは、格子定数を少しずつ振って全エネルギーが最小になる点を探して初めて言える話で、第8章と第11章で実際にやってみる。一般に PBE はこの平衡格子をやや大きめに見積もる傾向が文献で知られているが、私たちのシリコンでどれだけずれるかも、最小を探すまでは言えない。いまの celldm(1) は、あくまで前提の側に立っている。

次の nat = 2 は「この胞に原子が2個」。ここで第1章を思い出してほしい。立方の単位胞には8個入ると数えたはずだ。なのに2個とは、どういうことか。ibrav=2 が FCC を選んでいるからである。QE は立方体ではなく、その \(1/4\) の体積をもつ最小の繰り返し単位(原始胞)をそのまま胞として使う。だから原子は2個。立方体で数えれば8個、原始胞で数えれば2個。同じ結晶を、別の窓から数えているだけだ。

体積でも確かめておこう。FCC の原始胞の体積は立方体の \(1/4\)、すなわち \(a^3/4 \approx 5.431^3/4 \approx 40\ \mathrm{\text{Å}^3}\)。これを2原子で割ると、1原子あたり \(20\ \mathrm{\text{Å}^3}\)。第1章で密度から見積もった「原子1個 約 \(20\ \mathrm{\text{Å}^3}\)」と重なる。nat=2ibrav=2 の二行は、確かにあのシリコンを置いている。

残る ntyp = 1 は「原子の種類は1つ」、シリコンだけ、の宣言だ。そして ecutwfc = 60 は平面波のカットオフ――電子の波をどこまで細かい波まで使って表すかを決める目盛りで、単位は Ry である。ここでは 60 Ry を置いたが、この数字をいくつにすべきかは、まだ宙に浮いている。次章で正面から決める。すぐ下の ecutrho = 240 は、電子の波ではなく電荷密度の側のカットオフで、ここでは波動関数の4倍に取った(PAW を使うときの目安だ)。この比をどう選ぶかにも立ち入らず、ecutwfc と合わせて次章で扱う。

4 &electrons:つじつま合わせの作法

&system が「何を」なら、&electrons は「どう解くか」を決める。

電子は、自分たちがつくる電場の中を動く。動けば電場が変わり、電場が変われば動きも変わる。この鶏と卵を、計算機は反復で解く。ある電子配置を仮に置いて電場を計算し、その電場で電子を置き直し……を繰り返し、入口と出口がほとんど変わらなくなったところで「つじつまが合った(self-consistent)」とみなして止める。scf の名は、ここから来ている。

conv_thr = 1.0d-8 は、その「ほとんど変わらない」の基準だ。反復の前後でエネルギーの変化がこの値を下回ったら、ループを閉じてよい、という合図である。単位は Ry で、\(10^{-8}\ \mathrm{Ry} \approx 1.4\times10^{-7}\ \mathrm{eV}\)。前章で比べたかった全エネルギーの差が meV の桁だったことを思えば、ループはその千分の一よりさらに細かいところまで閉じている。mixing_beta = 0.7 は、新旧の電場をどんな割合で混ぜて次の一手に進むかの加減だ。

この区画は、いまは作法として受け取っておけばよい。反復がなぜ必要で、どんなときに収束をしくじり、その目盛りをどう選ぶか――それは第7章で、つじつま合わせそのものを主題にするときに正面から立ち返る。ここでは「電子状態は一発では決まらず、反復で詰める」という一点だけ、頭に置いておこう。

5 擬ポテンシャルという約束

ATOMIC_SPECIES の行には、元素記号 Si、原子量、そしてもう一つ、長い名前のファイルが並ぶ。

Si  28.0855 Si.pbe-n-kjpaw_psl.1.0.0.UPF

末尾の .UPF ファイルが擬ポテンシャルである。これが何のための約束かを、一段だけ説明しておく。

シリコン原子は14個の電子をもつ。だが結合に与るのは外側の4個(\(3s^2 3p^2\))だけで、内側の10個(\(1s^2 2s^2 2p^6\))は原子核にきつく縛られ、結晶になっても素知らぬ顔で芯に居続ける。ならば、その動かない10個を原子核と一括りにして一つの「実効的な原子核」とみなし、計算では外側の4個だけを相手にすればよい。これが擬ポテンシャルという約束だ。1原子あたり14個を4個に減らせるのだから、節約は大きい。芯の電子を律儀に解いても、結合についてはほとんど何も新しいことを教えてくれない――その見切りが、この約束を支えている。

便利だが、あくまで約束であることは頭の隅に留めておく。「芯は動かない」と決めてかかっているのだから、もしそれが破れる状況(極端な高圧など)に踏み込めば、約束の方が先に音を上げる。なお、ファイル名の pbe は、その擬ポテンシャルをどの近似のもとで作ったかの刻印だ。この pbe という名前は、本書の最後(第9章)でもう一度、別の顔をして効いてくる。いまは名前だけ覚えておけばよい。

6 原子を置く行、刻みを決める行 ― ATOMIC_POSITIONSK_POINTS

ネームリストの後ろに、残り二枚のカードがある。配置の仕上げと、周期の刻みだ。

ATOMIC_POSITIONS (alat) は、胞の中で2個の原子をどこに置くかを言う。(alat) は座標の単位が celldm(1)(=alat)であることを示す合図で、Si 0 0 0Si 0.25 0.25 0.25 の二行は、一方を原点に、もう一方を対角へ \(1/4\) ずらして置く、と読む。これが第2章で見たダイヤモンド型の基底そのものだ。この二行こそ、文字どおり「原子をどこに置くか」を書いた、配置のなかの配置である。

K_POINTS (automatic)8 8 8 0 0 0 は、無限に続く周期構造をどれだけ細かく標本化するかを指定する。前半の 8 8 8 は、ブリルアンゾーンを各方向に8分割した \(8\times8\times8\) の格子点で電子状態を拾う、という意味(Monkhorst–Pack のメッシュ)。後半の 0 0 0 は、その刻みをずらさない指定で、原点すなわち Γ 点を含む対称な格子の上で標本点を取る。これも ecutwfc と同じく「どこまで細かく刻むか」の目盛りで、ただし刻む対象が違う。ecutwfc波の細かさを、K_POINTS周期の標本化の細かさを決める。いくつにすべきかは第6章で決める。

こうして一枚を読み終えると、行がきれいに二色に分かれてくることに気づく。次の節で、それを自分の手で塗り分けてみよう。

7 二色に塗り分ける

flowchart TB
    F["si.scf.in<br>(一枚の供述書)"]
    F --> A["配置を決める行<br>どんなシリコンか"]
    F --> B["近似・規則を決める行<br>どう調べるか"]
    A --> A1["ibrav=2 / celldm(1)=10.26<br>nat=2 / ntyp=1<br>ATOMIC_POSITIONS"]
    B --> B1["calculation='scf' / ecutrho<br>ecutwfc / K_POINTS<br>conv_thr / 擬ポテンシャル"]
    classDef placement fill:#d4edda,stroke:#28a745,color:#000
    classDef knob fill:#d1ecf1,stroke:#17a2b8,color:#000
    class A,A1 placement
    class B,B1 knob
図 1: 一枚の入力ファイルを二色で塗り分ける。緑は結晶の配置そのものを置く行(どんなシリコンか)、青はその配置をどんな近似・刻みで調べるかを決める行(どう調べるか)。

この二色――結晶を置く行と、その結晶をどう調べるかを決める行――の見分けが、供述書を読む目の核心だ。次の小課題で、一行ずつ自分の手で振り分けてほしい。

小課題 4.1

本章冒頭の si.scf.in から、次の各行を「配置を決める行(どんなシリコンか)」と「近似・規則を決める行(どう調べるか)」のどちらかに仕分けよ。さらに、それぞれが結晶の何を決めているか、物理の言葉で一言ずつ添えよ。

  1. ibrav = 2 (b) celldm(1) = 10.26 (c) ecutwfc = 60 (d) ATOMIC_POSITIONS の2行 (e) K_POINTS 8 8 8 (f) conv_thr = 1.0d-8 (g) Si ... Si.pbe-n-kjpaw_psl.1.0.0.UPF(擬ポテンシャル)

解答例

配置を決める行(どんなシリコンか)

    1. ibrav = 2:格子の型を FCC に決める。結晶の骨格そのもの。
    1. celldm(1) = 10.26:格子の寸法(bohr)。結晶の大きさ。Å に戻せば 5.43 Å で、実験値 5.431 Å と合う。
    1. ATOMIC_POSITIONS の2行:胞の中の原子の座標。ダイヤモンド型の2原子基底。

近似・規則を決める行(どう調べるか)

    1. ecutwfc = 60:電子の波をどこまで細かい波まで使うか(波の解像度)。
    1. K_POINTS 8 8 8:無限の周期をどれだけ細かく標本化するか(周期の刻み)。
    1. conv_thr = 1.0d-8:つじつま合わせをどこまで詰めて止めるか。

境界線の上にあるもの

    1. 擬ポテンシャル:一面では「Si という元素」を指すから配置の側に見える。だがその中身は「芯電子を一括りにする」という近似――約束の側だ。一枚のカードが両側に足をかけている。この二面性に気づくこと自体が、供述書を読む練習になる。

よくある誤解

入力ファイルのフラグは、いったん計算が走り出せば用済みになる事務書類だ、と片づけたくなる。だが順序はむしろ逆である。ecutwfc=60K_POINTS 8 8 8 も、シリコンが生まれつき持っている性質ではない。それらは、私たちがシリコンに質問を投げかけたときの「質問のしかた」だ。だから出てきた数字は、いつでもこの質問のしかたと一対で読む。前章の言葉でいえば、二つのエネルギーを比べてよいのは、二枚の供述書のこの部分が揃っているときだけだった。入力は、結果に添えて初めて意味をなす前提のラベルであり、結果が出たあとにこそ効いてくる。

8 この章のまとめ

一枚のテキストを、物理の供述書として読み通した。

  • 入力ファイルは計算への命令書であると同時に、どんなシリコンを、どんな近似・刻みで調べたかを記した供述書であること。前章でいう「同じ土俵」は、この供述書の一致のことだった。
  • 行は二色に分かれること。配置を決める行ibravcelldm(1)natATOMIC_POSITIONS)と、近似・規則を決める行calculationecutwfcecutrhoK_POINTSconv_thr・擬ポテンシャル)。
  • 入力の数字が、まだ前提の側にあること。celldm(1)=10.26(bohr)を Å に戻すと 5.43 Å、nat=2ibrav=2 から1原子 20 ų。どちらも本物のシリコン(5.431 Å、密度 2.33 g/cm³)と合い、まだ計算の結果ではない。
  • 擬ポテンシャルは「芯電子を一括りにする」便利な約束であり、その約束の名前(pbe)は第9章で再会すること。
  • 出てきた数字は、この供述書と一対で読むこと。入力は結果に添える前提のラベルである。

二色のうち、配置の側は第2章までに置けるようになった。残る規則の側で、いちばん最初に宙に浮いたまま残っているのが ecutwfc だった。電子の波を、有限個の波の重ね合わせとして表す――そう決めた瞬間に、ひとつの問いが立ち上がる。では、どこまで細かい波まで数えれば足りるのか。粗すぎれば本物を取り逃し、細かすぎれば日が暮れる。次章は、この一つの目盛りをどこに合わせるかに、正面から取り組む。