Machikaneyama2000 の使用に関するメモ
大阪大学大学院理学研究科物理学専攻
赤井久純（１）はじめに
Machikaneyama2000 は KKR-CPA 法を用いて置換型合金・混晶の電子状態を密度汎関数法の局所密度近似の枠内で第一原理的に計算するための計算機コードである．CPA(coherent potential approximation)は Soven をはじめとする何人かの人々によって独立にほぼ同時期に考え出された，置換型不規則合金を扱うための近似法であり，統計力学における分子場近似によく似た考え方で不規則系の電子状態を要領よく高い精度で扱うことができる．CPA は最初タイトバインディング模型を念頭に開発されたが，その後，斯波と Soven によって独立にバンド計算の手法である KKR 法（グリーン関数法）に適用可能であることが見いだされ，KKR-CPA 法として発展した．ただし，KKR-CPA 法は普通のバンド計算とは比較にならないほど大規模な計算が必要とされるために，実際にKKR-CPA 法が第一原理計算として使われるようになったのは 1980 年代になってからである．
KKR-CPA 法については講義ノート"Korringa-Kohn-Rostoker Method - Brief summary of the theory ofscattering -"by H. Akai,
を参照いただきたい．教科書としては"Green Functions for Ordered and Disordered Systems"
by A. Gonis
North-Holland, 1994
ISDN 044488009868
をあげておく．
論文としては，歴史的に初めて第一原理計算に応用された例として
H. Akai, J. Phys. Soc. Jpn 51, 468 (1982)
初期の計算手法に関する議論は
H. Akai, J. Phys. Conden. Matter 1, 211 (1989).
全エネルギーや不規則局所磁気モーメントを計算する手法に関しては
H. Akai and P. H. Dederichs, Phys. Rev. B 47, 8793 (1993).
これら以前のものは上記論文の参考文献中にあげられている．
また，最近の磁性半導体に対する応用例としては
H. Akai, Phys, Rev. Lett. 81, 3002 (1998)
をあげておく．

（２）KKR-CPA コード
KKR-CPA コードは PC，ワークステーション，大型汎用計算機，スーパーコンピュターのいずれでも動作可能である．特に最近の PC の発展の結果PC でも余り不自由を感じないで計算が可能になっている．また，パラレル化やベクトル化も効率良く行うことができ，ベクトル化で２０倍程度，パラレル化でもノード数に応じて比較的簡単に高速化をはかれる．物性研のスーパーコンピュター１０ノードを使って２００倍程度の高速化が得られている．
KKR-CPA コードは不規則合金・混晶を想定して作られているが，規則合金や単体の系の計算も問題なく実行することができる．その場合 CPAは実行されず，純粋な系のバンド計算が実行される．ただし，多くのバンド計算で用いられているような対角化の方法にはよらず，グリーン関数を直接計算することよってバンド計算が実行される．したがって，純粋な系のバンド計算を実行すると同時にその系を母体とした不純物系の計算も実行するというようなことも容易に出来る．例えば鉄のバンド計算を実行すると同時に鉄中にニッケル不純物の計算を済ませることができる．
計算機に依存するサブルーチンは以下の２個である．
uclockudate
これらは，ソースファイル中のコメントに従って，適当に計算機に適合するように変更（多くの場合は単にコメントマークの付け直し）する必要がある．いずれも時刻を取得するルーチンであり，変更はほとんど自明である．
パッケージにはソースコード，例題としての input file，コンパイルするための makefile が含まれている．ソースはディレクトリーsource の下にすべておかれている．一つのサブルーティンにつき一つのファイルとなっている．ソースは FORTRAN77 を用いて書かれており，コンパイルは標準的なフォートランコンパイラーを用いて makefile に従ってコンパイルする．
main program は source/specx.f である。オリジナルには main program はsource/spec#.f であって，#のところには 1, 2, 3...等の数字が入り，この数字を扱うことの出来るユニットセル中の最大原子数に対応させていた．
例えば source/spec7.f はユニットセル中に７個の原子が含まれる系までを扱うことができる main program で一般に最大原子数を決めたとき，その系を扱うことの出来る spec#.f は spec1.f をもとにして作成していた．しかし現在使用のパッケージでは一律に specx.f としている．ユーザーは場合に応じてspecx.f を手直しして使うことになる．
specx は main program であるが，主要な役割は dimension の宣言をすることであり，それ以外には input data を読み込むプログラムを起動することと，実質的な main program である source/spmain.f を起動することだけである．

（３）データー入力
input はキーボードなどの標準入力からの入力となるが，ファイルを作っておいてパイプ入力するのが普通であろう．例として，ディレクトリーin に in/ymnal, in/fe3pt, in/fe70ni30, in/fe, in/inmn6mn0as0as等が納められている．in/ymnal にはコメント付きとなっており，一応それだけで理解できる形になっている．input ファイルを作る規則は次の通りである．
1) 自由形式であり，フォーマットはない．どこで行変えがなされてもよく，重要なのはデータの順序だけである．
2) データの区切りは空白またはカンマである．
3) 空白行は読み飛ばされる．
4) 先頭が"C", "c" あるいは"#"で始まる行はコメント行である読み飛ばされる．
このことはしばしばデータ入力の失敗を引き起こすので注意する必要がある．
たとえば炭素に C という名前をつけてデーターとして行の先頭から入力すると，コメントとみなされて読み飛ばされることになる．
5) データは省略することができる．ただし，省略したデータに対応して区切りのカンマを入れる必要がある．省略したデータに対してはデフォルト値が与えられる．
6) １セットの入力データのあとに更に入力データのセットを続けておくことが出来る．それぞれのセットは独立した計算として実行される．

7) 以下で説明するキーワードについては大文字はすべて小文字に変換された上で解釈される．
以下に YMnAl 合金の場合について入力データセットの例をあげて説明する．c で始まる
行はすべてコメント行で実際には読み飛ばされる．

c--------------------------- input data ---------------------------c go/ngo/dos/dsp/spc
file namec----------------------------------------godata/ymnal
c fcc/bcc/hcp/sc/... a b/a c/a alpha betac--------------------------------------------------------fcc
14.0 1.0 1.0, , , ,gammac edelt width nrl/sra mjw/vbh/vwn mag/nmag/-mag/rvrs/kickc-----------------------------------------------------------------------0.001
1.
nrl
mjw
nmag

init/1st/2nd
2nd

c quit/update 0/1/2/../t/l/m/h/u iteration mixingc----------------------------------------------------update
t40
0.024
c number of typec----------------2
c type components rtin field l_max
c-------------------------------------------------------------Y
100.02MnAl
200.02Z

concentration

3925
13

c number of atomsc-------------------6
c positionc----------------------------------0.125a 0.125b 0.125c
0.875a 0.875b 0.875c
0.5a
0.5b
0.5c
0.0a
0.5b
0.5c
0.5a
0.0b
0.5c
0.5a
0.5b
0.0c

typeYY
MnAl
MnAl
MnAl
MnAl

c crystal axis system is also possible, e.g. 0.5a, 0.5b, 0.5c.

100955１）最初のデータは go, ngo, dos, dsp, spc のいずれかで，以下のデータに対する処理の種類を指示する．それぞれの意味は次の通りである．
go 通常の self-consistent な計算を実行する．
ngo 何もしない．すなわち，続く１セットに相当する一連のの入力データは単に読み飛ばされる．
dos self-consistent な反復計算をしないで，最新のポテンシャルデータを用いて状態密度が計算される．したがってあらかじめ通常の self-consistent な計算を完了しておかなければ意味のある結果は得られない．
dsp 途中の結果が出力されす．計算結果がおかしいときなどにチェックのために使用．
spc 通常のバンド計算の dispersion relation に相当するものは，不規則合金ではBloch spectral function である．dispersion relation が与えられた k ベクトルに対するエネルギー固有値に対応したデルタ関数であるのに対して，不規則性による散乱によってエネルギー固有状態が寿命をもつために幅をもつようになる．そのような情報を持ったBloch spectral function を出力する．dos と同様に self-consistent な計算を行わずに最新のポテンシャルデータを用いて計算される．
２）このデータはポテンシャルデータを入出力するファイルを指定する．このファイルは実際に存在すればそのファイルが用いられるし，存在しなければ新たに作られる．
３）１４個のブラベー格子の一つを指定する．(1)fcc (2)bcc (3)hcp (4)sc (5)bct (fct)(6)simple tetragonal (7)face centered orthorhombic (8)body centered orthorhombic(9)base centered orthorhonbic (10)simple orthorhombic(11)base centered monoclinic (12)simple monoclinic(13)triclinic (14)rhombohedral (trigonal) (15)fct (bct)に対してそれぞれ，fcc, bcc, hcp (hex), sc, bct, st, fco, bco, bso, so, bsm, sm,trc, rhb (trg), fct のように指定する．(5)bct と(15)fct は同じものであるが，fcc と bcc のどちらにより近いかということで fct と bct を使い分けることもできる．hcp と hex，rhb と trg はまったく同じものである．
４）～９）格子定数を指定する．a, c/a, b/a, alpha, beta, gamma を指定するが，たとえば立方晶系ならば強制的に c/a=b/a=1, alpha=beta=gamma=90 にとられる．実際に基本格子ベクトルと逆格子ベクトルがどのように作られるかはサブルーチン prmvec をみればわかる．また，c/a，b/a などは多くの場合省略することができ，その場合には省略の規則に従って該当する数だけのカンマが入力される．また，a を省略するかあるいは a=0 とすると適当な理論値がまた，a として１０００００より大きい値を与えると実験値が用いられる．
１０），１１）delta と width が入力される．delta の意味はフェルミエネルギーに付与される微小な虚数部分である．本プログラムパッケージで計算されているのは基本的にグリーン関数である．グリーン関数はその定義から複素エネルギー面上の正の実軸はブランチカットになっており，グリーン関数は定義されない．複素平面状の虚数部分の正の側から実軸に近づいた極限として実エネルギーでのグリーン関数が定義されるが，数値計算上は有限の虚数部分を与えて計算するが，そのときの虚数部分の大きさである．普通は 0.001 程度の値を与えて計算される．
なお，フェルミーエネルギー以外のエネルギーについては一般の複素エネルギーを用いて計算がなされており，当然大きな虚数部分をふくんでいる．また，width の意味はグリーン関数が計算されるエネルギー領域の広さである．たとえば，width=1 であればフェルミエネルギー以下1Ry のところからフェルミエネルギーまでが計算される．width の大きさは価電子帯を完全にカバーするものでなければならない．価電子帯の幅がわからない場合には大きめの width を用いて荒い計算をしておき，状態密度を出力してその結果から価電子帯の幅を見積もり，適切なwidth を選ぶようにすれば良い．

１２）nrl または sra というオプションを指定する．nrl は非相対論的計算を実行，sra はスカラー相対論的な計算（すなわちスピン軌道相互作用を落とした相対論的計算）を実行．
１３）mjw, vbh, vwn の一つを指定．いずれも密度汎関数法の局所近似の parametrization の方法で，mjw は Moruzzi, Janak and Williams によるもの，vbh は von Barth and Hedin によるもの，vwn は Vosko, Wilk and Nusair によるものである．実際にはこれ以外に GGA にも対応しており指定することができるが，その動作については保証していない．詳しくはサブルーチン
potenv のコメントおよび中身をみればわかるようになっている．
１４）mag, nmag, -mag, rvrs, kick のいずれかを指定できる．mag は磁性をもつ解を得たい場合に指定．強磁性や反強磁性やスピングラス，非磁性母体中の磁性不純物などを扱うことができる．nmag は非磁性解を得たいときに使用．非磁性の物質や高温で nonmagnetic になっている状態の計算やテスト計算に用いる．計算時間は半分ですませることができる．
-mag と rvrs は同じものであり，計算の結果や途中で，磁化が負になった場合（それでも一向に差し支えないわけであるが）それを正にしたいときに指定すると磁化の符号を逆転することができる．非磁性状態で計算されたポテンシャルデータを用いて磁性状態を計算しようとする場合，オプションを nmag から mag に変えれば良いわけであるが，それだけでは self-consistent な計算を続けても非磁性状態の結果が求まってしまう（非磁性状態はどんな場合でも解のひとつになり得るから）．スピン空間でも対称性をいったん壊してやらなければ磁性状態の解へと移ってくれない．このようなときに kick というオプションを用いれば人為的に対称性を壊し，磁性状態へと移行することができる．
１５）init, 1st, 2nd のいずれかを指定する．ポテンシャルデータは通常２組保存される．
最初の１組（レコード）は二番目のレコードを作り出すための初期値として用いられたポテンシャルデータが収められている．このように二組のポテンシャルデータを保有することによって，最新のデータと一つ前のデータを参照することができる．これは，たとえば一連の計算によって求まっているポテンシャルデータを用いて，さらに self-consistent計算を行い精度を上げようとしたが，結果が思わしくないので，一つ前の結果に戻りたいなどというときに用いることができる．2nd を指定すると最新のデータを，1st を指定するとそれより一つ前のデータを初期ポテンシャルデータと用いて計算することになる．また，
init を指定すると保存されたポテンシャルデータを用いずに，原子計算から初期ポテンシャルを作り出して，計算をはじめる．省略すると 2nd が仮定されるが，ポテンシャルデータが存在しない，あるいはファイル自身が存在しない場合には強制的に init として計算される．
１６）quit, update のいずれかを指定できる．update を指定すると計算によって得られたポテンシャルデータでファイルのレコードを更新する．更新の仕方は，ファイル上の２番目のレコードで１番目のレコードを置き換え，計算された最新のポテンシャルデータが２番目のデータとなる．quit を指定した場合にはファイルの書き換えは一切起こらない．１）でdos, spc を指定した場合には，たとえここで update を指定していても強制的に quit に置き換えられる．
１７）ブリルアン域内の積分に用いられる k ベクトルの数に関係したオプションを指定する．
整数もしくは t, l, m, h, u のいずれかのデータを与える．整数の場合にはブリルアン域の(0,0,0)から(1,0,0)までの分割数を与える．0 を用いると(0,0,0)の１点のみが用いられる．
また，t は test の意味で 0 を指定したのと同じである．l, m, h, u はそれぞれ，low, medium,high, ultra high quality を意味し，この順番に積分に用いられる k ベクトルの数は増加する．l, m, h, u のそれぞれがどのような整数値に相当するかはブラベー格子によるが，サブルーチン nfqlty に与えられている（かなりいいかげんな値）．

１８）self-consistent な反復計算の最大反復数を与える．たとえば４０とすれば４０回の反復計算の後計算は停止する．それ以前に self-consistent な収束条件を満たせばその時点で停止する．収束条件はメインプログラムの中で与えられており（data 文中の tol の値で 1d-4～1d-6 ぐらいが適当）適当に変更することができる．一般には結果をみながら，収束しているようであるならば，収束条件を満たしていなくても十分であり，収束条件はあくまで単なる目安である．
１９）反復計算のためのいわゆる mixing parameter である．チェビシェフ加速による反復計算を行っているために通常の mixing parameter ではないが，大きいほどポテンシャル更新に程度が大きいと考えて良い．収束する限りにおいては大きい方が収束が早いが通常の反復計算と同じで，大きすぎると容易に発散する．最適な大きさは問題に依存するが，経験上は 0.01～0.03 程度が良いようである．困難な問題の場合（収束が難しい）0.004 程度あるいはそれ以下が用いられる．
２０）これ以降はユニットセルに含まれる原子の情報である．まず最初に与えるデータは含まれるサイトの種類の数である．たとえば YMn_0.95Al_0.05 の場合，Mn と Al は互いに固溶しておりこれは Mn が 95%と Al が 5%からなる一つのサイトとみなせる．また Y は単独で一つのサイトを作っている．従って，この場合にサイトの数は二個である．
２１）～２７）これらのデータはサイトの種類の数だけ繰り返して与えられる．まずサイトに任意の名前を与える．ここでは１番目のサイトに対して Y という名前をつけている．
次に，このサイトはただ一つの原子 Y からなるので，component の数 1 を与えている．
rtin はマフィンティン球の半径であるが，特に指定する必要がなければ 0 を与えておけばtouching radius になるようにプログラムが自動生成する．また field は外部磁場を与える場合で，磁場によるゼーマン分裂の大きさを Ry 単位で与える．l_max はそのサイトのポテンシャルによる散乱が l=0 から l=l_max まで考慮されることを意味する．たとえば
l_max=2 ならば d 散乱までが考慮され，それ以上の部分波に対する散乱はない（すなわち平面波と同じになる）として取り扱うことを意味する．Z はそのサイトを構成する原子の原子番号でその次の concentration とペアで指定する．Y サイトを構成するのはただ一つの原子 Y であるから Y の原子番号 39 と concentration 100 を与えている．ただしconcentration はサイト内で規格化されて使われるので，この場合にはどのような数字を与えても結果はおなじである（すなわち規格化されて最終的に数字 1 に変換される）．
これで１番目のサイトに関するデータはおしまいである．次に二番目のサイトであるが，このサイトには MnAl という名前を与えた．component の数は Mn 原子と Al 原子があるので２である．rtin=0，field=0 としている．また構成原子として Mn である Z=25 と Al である
Z=13 を与えている．それぞれの concentration は 95%と 5%であるので，95 と 5 という数字をあたえているが，これは 0.95 と 0.05 などと書いても結果は規格化されるので同じである．
２８）サイトの種類に関するデータをすべて入力すると次はそれらのサイトの位置に関する情報を与える．まず，ユニットセル内に存在するサイトの数を入力する．同じ種類のサイトであっても位置が違えば異なった位置として勘定しなければならない．今の場合ユニットセルあたり６個の違った位置があるので数字 6 を与える．

２９）～３２）異なった位置のそれぞれについての情報を与える．最初に位置を与える．
３通りの与え方があり，今の場合に様に 0.125a 0.125b 0.125c のように与える方法，
0.5x 0.5y 0.5z のように x, y, z が数字のあとにつく与え方，およびとただ単に数字（たとえば 0.5 0.5 0.5 のように）を与える方法がある．数字の後に a, b, c がついている場合にはそれは基本ベクトル a, b, c の数値倍と解釈される．数字の後に x, y, z がついている場合には x, y, z 方向の単位ベクトルの数値倍と解釈される．また単なる数値の場合には１番目，２番目，３番目の数字に対して x, y, z が省略されているものと解釈される．実際に指定される位置はベクトル 0.125a+0.125b+0.125c あるいは 0.5x+0.5y+0.5zなどとなる．また，これらは混在していてよい．たとえば，0.5 0.5y 0.125c などの指定が可能である．また，a, b, c, x, y, z が単独で現れた場合にはアルファベットの前に数字 1 があるものと解釈される．これらの位置の指定のあとにそのサイトの種類を指定する文字が与えられる．たとえば，いまの場合６個の異なった位置があるが，１番目と２番目の位置はサイト Y であり，位置３～６はサイト MnAl である．一般には同じ種類のサイトがユニットセルの中に複数回現れて良い．また，２１）～２７）で定義していないサイトがここに現れてはならない．
一番最初のデータ１）で spc 以外を指定した場合は以上で１セットのデータの入力が終わりである．このようなデータのセットは任意個続けることができる．つまり，サイトの位置の入力のあとから，再び go などで始まるデータを置いて良い．計算はデータのセットがなくなるまで連続して行なわれる．
spc を指定した場合には例外で，この後に Bloch spectral function を計算したい k 点（k ベクトル）を任意個数指定する．たとえば0000.1 0 00.2 0 00.3 0 00.4 0 00.5 0 0のように入力するとガンマ点から X 点までの k 点について 0.1 間隔で spectral function が計算される．このような k ベクトルのデータがなくなるまで（つまり file の最後が検出されるか次の計算の始まりを示す go などのデータが検出されるまで）計算される．
以上が YMnAl の場合であるが，inmn6mn0as0as の場合も input file の意味を探ることは容易である．ここでは原子番号として Z=0 の指定があるが，これは原子空孔を意味する．
ダイヤモンド構造や zincblende 構造の結晶の場合オープン構造であり，計算の精度をあげるために原子の無い位置にも空孔を入れているが，必ずしも必要なわけではない．

もっとも単純な例として純粋な鉄の場合を以下に示しておく．この場合には不規則合金出ないために，自動的に cpa の部分ははずされて計算するために計算時間は合金の場合に比べてかなり高速になる．コメントがはずされているためにコンパクトになっているが，本質的な部分は変わり無い．
go data/fe bcc 5.27,,,,,,
0.001 1.0 nrl mjw mag 2ndupdate 6 80 0.024
1Fe 1 0 0 26 2 10010 0 0 Fe（３）出力出力の例がディレクトリーout に収められている．上の例 in/fe を用いて計算された結果でそれぞれの機械でコンパイルした実行ファイルを用いて得られた結果と比較できるようになっている．結果は 80 回の反復計算の後 no-convergence で終了しているがこれは問題ない．反復の経過をみればわかるとおり，磁化や全エネルギーは 60 回程度の反復計算で実用上完全に収束している．エネルギー値，磁化等はユニットセルあたりでエネルギーはすべて Ry 単位，磁化はボーアマグネトンを単位にしている．その他，各原子のげマフィンティン球内に含まれる角運動量別の電子数や球内の磁化，超微細磁場やアイソマーシフトに相当する原子核上での電子密度，コア電子のエネルギー準位などが出力される．
（４）メインプログラムの生成について
spec1.f から任意の spec#.f を生成する方法について説明する．specx.f のパラメータ文の変更を行う．変更するパラメータはparameter& (natmmx=1, ncmpmx=1, msizmx=25, mxlmx=4, nk1x=550, nk3x=21,& msex=201, ngmx=15, nrpmx=250, ngpmx=250, npmx=200, msr=400)の内で，natmmx, ncmpmx, msizemx であり，natmmx がユニットセル内の異なった位置の数に対応する．たとえば，先の YMnAl では２８）で指定した 6 が natmmx より小さくなければならない．また，ncmpmx は２１）～２７）に現れる原子の種類の数に対応する．たとえば，
いまの場合は Y と Mn と Al があらわれるので，ncmpmx は３あるいはそれ以上であれば良い．注意しなければならないのはたとえ同じ原子であっても違う種類のサイトに扱われる場合はちがう原子として勘定しなければならない．また，同一の種類のサイトと構成する場合でも異なった原子と考えなければならない場合がある．
たとえばスピングラスで，磁気モーメントが上を向いた Mn 原子と下を向いた Mn 原子を区別して考えるときで，inmn6mn0as0as はその例である．sizemx はユニットセル内の異なった位置のそれぞれで考慮される部分波散乱の数をすべて足し合わせた和である．たとえばすべての位置で d 波散乱（部分波の数は 3^2＝９個）まで考慮する場合には９×６＝５４より大きければ良い．

上記のパラメータの値は最大値を与えるもので，たとえば上記の例で与えられる
main program を用いればユニットセルに６個の異なった位置と３個の異なった成分原子をもつ YMnAl を扱うことも１個の位置と１個の成分原子しかもたない Feを扱うことも可能である．しかし，上記のパラメータが大きいほど必要とするメモリーは大きくなるので，大きなプログラムですべてをまかなうのはリソースの消費につながる場合がある．
これら以外に場合に応じて調整を要するパラメータは nk1x, nk3x, mxl である．
nk1x+nk3x はブリルアン域内の積分に用いられる k 点の数の上限を与えている．
上の例では 581 点まで k 点の数を増やすことができる．一般にはユニットセル内の異なった位置の数が増えれば k 点の数を減らすことができる（また，そうしなければ計算時間が増えて計算が困難になる）．たとえば，ユニットセルに６個の異なった位置が含まれる場合には１０～４０点程度で十分な結果が出せる．テスト計算
ならば数点でもかまわない．nk3x は Bloch spectral function を計算する k 点の数の対応するが，普通 1 程度で十分である．mxl は計算に用いる最大の部分波散乱角運動量+1 である．s, p, d 状態まで計算する場合には mxl=3，f 状態まで計算する場合にはmxl=4 を用いる．希土類やアクチナイドを除いて定性的な議論には mxl=3 で十分であろう．
data 文で変更の可能性があるものは，mse0 と tol である．現在 mse0=35 になっている．
これはエネルギーに関する複素積分を実行する際の複素平面状のメッシュの数に対応するが，高精度の計算が必要な場合には適宜増やすことが可能である．もちろんそれに比例して計算時間も増加する．また，tol は self-consistency の基準を与える
パラメータで 1-d6 は電荷密度の rms (root mean suare) error が 100 万分の 1 という意味である．出力に現れる err はその常用対数である．
上のようなサイトの位置の個数に応じたプログラムの変更を行ったときは，program 名をprogram spec6のように付け替えを行って，それを source/spec6.f という名前をつけて別のファイルとして save しても良い（本来そのようにしていた）．また，そのとき同時に makfile の１行目のprogram = specxを変更してprogram = spec6のように変更して make する．その結果 spec6 のような実行ファイルが生成される．