バンド計算プログラム＋用語集（拡張編） [Top]

目次 [Top]
はじめに
ストレス編
磁性＋交換相関項＋局所密度近似編
ＧＧＡ編
スピン軌道相互作用（失敗）編

はじめに

最初の思惑や、方針とは裏腹に、どんどん変な方向に向かっています。

ストレス編 [目次]

（１０／２９、２００４）磁性編の前にストレス編を置く。

ストレス計算に関しての覚書（備忘）。

磁性＋交換相関項＋局所密度近似編

導入 [目次]

既にご存知のように本ＷＷＷサイトで公開されている第一原理分子動力学計 算用プログラム、revpe_d.fでは、常 磁性状態でのバンド計算しか行なうことができません。
既に強磁性等、スピンが分極化した場合に対応するプログラムは完成してい ましたが、まだ公開していませんでした。ここでは、この磁性問題対応版プロ グラムの公開と説明（開発方法）を行なうことを第一の目的として話を進めて いきたいと思います。
現時点では、磁性状態としては、単純なFe,Co,Niの強磁性状態の計算テスト しか行なっていません。より複雑な磁気構造（反強磁性、フェリ磁性など）に 関しての計算、経験の蓄積は全くありません。

磁性問題への拡張 [目次]

では実際に現行の、revpe_d.fまたはrevep_d_new.fなどを、磁性を扱えるよ うにする具体的な方法について、ここでは議論したいと思います。磁性を扱え るようにするということは、これまで全く考慮していなかったスピンのアップ とダウンを別々のものとして考えるということです。常磁性のときは１つのバ ンドに２つの電子（アップとダウン）があると考えましたが、磁性を考える場 合は、１つのバンドには１つの（アップ、ダウンいずれか）が対応するように なります。
ではプログラムのどの部分を変更する必要があるのでしょうか。

最も重要なのは交換相関項の部分です。これを磁性考慮版に変更する必要が あります。revpe_d.f等ではWignerの表式 を用いていましたが、残念ながらこの表式では磁性を考慮することができませ ん（そのような表式が用意されていない）。従って、Wigner以外の交換相関項 形式を採用する必要があります。この場合の詳しい表式に関しては、文献「密 度汎関数法とその応用‐分子・クラスターの電子状態」、菅野 暁監修、里子 允敏、大西楢平著、講談社サイエンティフィック刊、４．５．１節、１１２頁 周辺に幾つかの表式が載っていて大変参考になります。

本計算では、ノルム保存型の擬ポテンシャルを使用していますが、基本的に は何ら変更を加える必要はありません。但し、本擬ポテンシャル（ＮＣＰＳ９ ７、２Ｋ）は、その大半が先のWignerの表式を採用しています。従って、バン ド計算部分との整合性を考えれば、擬ポテンシャルもバンド計算側に対応する 交換相関形式でのポテンシャルの再構築の必要性があります。ただここで扱う 擬ポテンシャルは、unscreeningしたものを用いるので、擬ポテンシャル側と バンド計算側の交換相関項の形式が合わなくとも、そう深刻な影響はないと思 われます。
交換相関部分に渡す情報は電荷密度に関するもので、分極（偏極）に関して のパラメーターζ、

ζ = (ρ↑ - ρ↓)/(ρ↑ + ρ↓) = (ρ↑ - ρ↓)/ρ

ρ = ρ↑ + ρ↓

ρ：全電荷密度、ρ↑：up spin、ρ↓：down spin

及び、rs(ρ)パラメータ（無単位量）。無単位量のため本当は、 rsでのρは電荷密度ではなく、電子の（数）密度である。ハート レー原子単位では電荷の単位が、e2 = 1（リュードベリ原子単位 では、= 2）なので、暗黙の内に電荷がかかっていると考える。分極に関して のパラメーターは、常磁性（ρ↑ = ρ↓ = 1/2ρ）での交換相関項の表式と、 スピンが完全に分極した場合（ρ↑ = ρ, ρ↓ = 0またはその逆）の表式を 内挿するために使われる。つまり、交換相関エネルギーεxcはス ピン分極を考慮した場合、

εxc[ρ,ζ] = εxcparr(ρ) + ( εxcferro(ρ) - εxcparr(ρ) )f(ζ)

εxcparr(ρ)：スピン分極していない場合
εxcferro(ρ)：スピンが完全に分極している場合

f(ζ) = { (1 + ζ)4/3 + (1 - ζ)4/3 - 2 }/{ 2(21/3 - 1 ) }

という内挿表式が使用される。常磁性（ρ↑ = ρ↓ = 1/2ρ）なら、ζ = 0であり、f(ζ = 0) = 0となる。従って、

εxc[ρ,ζ] = εxcparr(ρ)

となる。一方、完全に分極した場合（ρ↑ = ρ, ρ↓ = 0またはその逆）、 ζ = 1またはζ = -1となり、f(ζ) = 1となる。従って、

εxc[ρ,ζ] = εxcferro(ρ)

となる。それ以外の分極した場合は、parrとferroが上記内挿式にあるよう に適当に内挿された交換相関エネルギーの値となる。交換相関ポテンシャルも ほぼ同様な内挿が行なわれる。

アップ、ダウン別に擬ポテンシャルを構築する必要があるのは、スピン軌道 相互作用を考慮する場合で、単に（分極した）磁性を扱う場合は、そのような ことを考慮する必要はありません。

交換相関項の部分以外に変更を必要とするのは、波動関数に関しての部分で す。波動関数はアップ用、ダウン用の二つが必要になります。このため、総電 荷密度以外に、アップ用、ダウン用の電荷密度も必要になります。当然、計算 全体でもアップ用の計算、ダウン用の計算の二つに分かれます。このループ （アップとダウンで２回しか回らない）はバンド計算において、一番外側で、 最大のループといえます（ｋ点のループの外側になります）。

一番簡単で、分かりやすいのは、ｋ点ループの外側に、スピンに関しての （アップ、ダウン）のループを作ってしまう方法です。


      DO 100 ISPIN = 1,IFP+1
      Paramagnetic case --> IFP = 0, Ferromagnetic case --> IFP = 1 
      ------
      DO 200 IK    = 1,KNV
      ------
      WAVE_FUNCTION(IK,ISPIN)    = UP or DOWN
      ------
  200 CONTINUE
  100 CONTINUE
ループに関しては別の手段もあります。それはｋ点のループを２倍にして、 それぞれ前半をアップ用、後半をダウン用（逆でも別によい）とする方法です。 


      DO 100 IK    = 1,2*KNV
      ------
      WAVE_FUNCTION(    1,  KNV) = UP
      WAVE_FUNCTION(KNV+1,2*KNV) = DOWN 
      ------
  100 CONTINUE
既存の磁性非考慮のプログラムから、磁性考慮のプログラムへ拡張する場合、 後者の方法が最も簡単に取り掛かることができると想定したのですが、実際に は前者のスピンに関してのループを新たに導入する方法を選びました。

後者の方法を選ばなかった理由は、最初ｋ点ループ数を２倍にしてアップと ダウンの計算を割り振ろうとしたのですが、その制御が思った以上に面倒そう だったので（例えば、常磁性、強磁性の場合分け）、むしろ一番外側にスピン 用のループを作った方が楽だと判断しました。
（８／２、２０１８）ｋ点のループを２倍にする方法は、並列化（ｋ点並 列）にとっては有利に働くと思われます。系が大きくなると、相対的に必要 なｋ点数が少くて済むようになり、並列化にとっては不利になります。ｋ点 が２倍ならば、その分並列化による効果が効くので有利に働きます。

（８／１７、１９９８）最近気付いたのですが、ｋ点のループ以外にバンド のループ数を２倍にする手もあります。バンドはアップスピン用のバンド、ダ ウンスピン用のバンドと分かれているので、ｋ点を２倍にするよりかは拡張し 易いと思われます（ただ、筆者はこの方法は試していません。正直現在思い付 いた段階なので、思わぬ困難があるかもしれません）。

ただ、良く考えると、バンドに関しての配列変数EKO(KEG,KNV)は、先の前者 の方法でもEKO(KEG,KNV,2)となっています。EKO(2*KEG,KNV)として（ｋ点では なく）バンド関係のループを２倍にすることに利点があるのかどうか、現在の ところ不明です（あまりないようにも思えます）。

これら以外の方法としては、スピンアップ部分とダウン部分の計算を独立で 行なってしまう手があります。スピンに関する計算も独立性は高く、ほとんど の計算部分で（アップ、ダウン）互いに影響を及ぼすことはありません。アッ プ、ダウン両方の部分が必要なのは、交換相関項の計算部分、フェルミエネル ギー決定部分と全エネルギー計算部分、力やストレスの計算部分などです（電 子状態を計算するだけなら、力やストレスの計算は必要としません）。力、ス トレスの部分を除けば、計算ルーチンとしての量及び計算時間として占める割 合は大して多くはありません。従って、バンド計算のほとんどの部分をアップ、 ダウンスピン別々で計算できます。

これは並列計算からみても有利で、アップ、ダウン別々のＣＰＵで並列で計 算させることが可能です（但し、２ＣＰＵ並列しかできないのが欠点）。ＦＬ ＡＰＷ法のパッケージで有名なWIEN97 （現在は、WIEN2k）は、アップ、ダウ ンを別々に計算する（場合によっては並列での計算も可能か？）場合のシェル プログラムが標準で用意されています（いるはず）。

スピンのループを作る場合、DO 100 ISPIN = 1,IFP + 1 とかして、変数IFPを、常磁性のときはIFP = 0、強磁性の ときはIFP = 1と場合分けする方法が考えられます。このときIFPは、プログラ ムの最初で指定する（この場合、常磁性、強磁性でいちいち再コンパイルする 必要があります）方法もありますが、筆者は入力データ上（CONTL.DATの磁性 対応版、まだ公開していません）で制御するようにしています。こうすれば再 コンパイルの必要はなくなります（本当は、磁性考慮、非考慮で配列のサイズ が変わりますが、少々無駄を承知で最初からスピンアップ、ダウン分として配 列を宣言しておけば再コンパイルの必要性はなくなります。当然、大規模な系 ではこの手は使えません）。

実際、スピン用のループを作ると、既存の波動関数（固有ベクトル）や電荷 密度などの変数も、スピン用に書き直す必要が生じます。例えば、固有ベクト ル用の配列変数、ZAJ(KNG1,KEG,KNV3)は、ZAJ(KNG1,KEG,KNV3,2)となります。他にも電荷密度関係（波動関数がアッ プ用、ダウン用になるので、自ずと電荷密度もアップ用、ダウン用が必要とな ります。当然全電荷密度用の配列も必要です）、サブルーチンCONV2での電荷 密度混合用の作業配列変数、固有値（EKO）、バンドの占有数（OCCUP）、サブ ルーチンFORCEなどで使用している作業用変数（ZFBB,ZFC）、交換相関項など の配列（ここが一番面倒）がスピン用に書き換えたり、新たに作る必要があり ます。

（７／７、１９９８）磁性を扱う場合、最初のアップ、ダウン別の初期電荷 密度は、どちらかが多くなるように設定します。この時、アップ、ダウンを同 じにするとパラマグネティック（非磁性）の計算と同等になります。もし、Ｌ ＳＤＡ計算が（ＬＳＤ近似の範囲内で）正しく進むなら、最終的にアップ、ダ ウン異なった電荷量になって収束します（分極化した方が安定なら）。磁性を 扱う場合のアップ、ダウンそれぞれの電荷密度の設定はアップ側を全体の７５ ％、ダウン側を２５％にして計算しています（強磁性NI,Co,Feの場合）。この 設定はこれがベストという訳ではなく、５０％、５０％以外ならどのような設 定も可能です（が、あまり極端なものでは、計算が破綻したり、収束しない可 能性があります）。
（１１／３０、２００４）反強磁性やフェリ磁性を扱う場合は、初期電荷の 設定にもう一工夫が必要にあります。それは反強磁性などに向かうように最初 から電荷密度の分布を反強磁性的なものにして置きます。実際にはこれによっ て得られた反強磁性状態と他の磁気状態との全エネルギーの比較で、いずれの 磁気構造が安定か判定することとなります。初期値の設定によっては、反強磁 性的な分布を出発点としても、常磁性に落ち着いてしまう可能性もあります。 また、反強磁性を扱う場合その対称性が常磁性、強磁性と異なる（単位胞も２ 倍以上となる）ものとなります。

（１２／８、１９９９）磁気モーメントの求め方は、用語集をご参照下さい。

（７／７、２０１１）磁気構造の模式的な表現
参考図１ 〔png形式、6 kb〕：摸式的な状態密度で、強磁性、常磁性、反強磁 性、ハーフメタリックを表現したもの。強磁性、常磁性、反強磁性の各状態は 必ずしも金属的である必要はない（＝ギャップが空いてもかまわない）。注： 図とウェブページ上の文章とは必ずしも整合していない。
参考図２ 〔png形式、7 kb〕：↑(Up spin)、↓(Down spin)で、強磁性、フェ リ磁性、反強磁性、ノンコリニアな場合の各磁気構造を２次元的な平面上で表 現したもの（あくまで摸式的な表現であり、厳密性に欠けることに注意）。平 面を区切っている四角が最小の単位（単位胞）に相当。注：図とウェブページ 上の文章とは必ずしも整合していない。

スピン分極：交換エネルギーの 場合 [目次]

（１２／２、２００４）スピン分極する場合の交換エネルギー部分について 以下に説明を行なう。
交換エネルギーEx(ρ)について、これは、

Ex(ρ) = ∫ρεx(ρ)dr ・・・(1)

と表現でき、εx(ρ)を交換エネルギー密度と言い、

εx(ρ) = Aρ1/3 ・・・(2)

とする。Aは適当な係数。（＊）一様な電子ガスでの交換エネルギーの形は、 上記のようなものとなる。
従って、Ex(ρ)は、

Ex(ρ) = A∫ρ4/3dr ・・・(3)

となる。これをスピン分極（偏極）していない場合の交換エネルギーの表式 とする。次に、スピン分極した場合を考える。この場合の交換エネルギーを、 Ex(ρ↑,ρ↓)とすると、

Ex(ρ↑,ρ↓) = ○A∫{ ρ↑4/3 }dr + ○A∫{ ρ↓4/3 }dr
= ○A∫{ ρ↑4/3 + ρ↓4/3 }dr ・・・(4)

と表現される。○は未知の係数とする。常磁性状態の場合（ρ↑ = ρ↓ = ρ/2）を上記分極した交換エネルギーの式に適用すると、

Ex(ρ/2,ρ/2) = ○A∫{ (ρ/2)4/3 + (ρ/2)4/3 }dr
= ○A∫{ 2(1/2)4/3ρ4/3 }dr
= ○2(1/2)4/3A∫ρ4/3 dr ・・・(5)

となる。この(5)の表式は、(3)の表式と一致するはずで、そのためには未知 の係数○が、○ = 21/3でないとならない。従って、スピン分極を 考慮した交換エネルギーの表式(4)は、

Ex(ρ↑,ρ↓) = 21/3A∫{ ρ↑4/3 + ρ↓4/3 }dr ・・・(6)

となる。ここで重要なのは、分極を考慮しない場合と考慮する場合で表式が 異なることである。従って、分極を考慮しない表式をそのままスピン分極する 場合に適用することは磁性状態に関して間違った結果を導く可能性がある。相 関エネルギー部分は表式がより複雑であり、分極の考慮、非考慮の使い分けは より重要となる。

　式(6)は、分極（偏極）度ζを使って、ρ↑=(1+ζ)/2, ρ↓=(1-ζ)/2から、

Ex(ρ↑,ρ↓) = (1/2)A∫ ρ4/3[ (1 + ζ)4/3 + (1 - ζ)4/3]dr ・・・(7)

とも書き表せる。ζ=0なら(3)式に帰着し、ζ=1（またはζ=-1）なら、(6) 式となる。

交換エネルギー、交換ポテンシャ ルの係数 [目次]

（２／１７、２００５）電子ガスにおける交換エネルギー、交換ポテンシャ ル部分の係数に関して以下に説明を行なう。この時、スピンは考えないことに する。単位は、ハートレー原子単位を 用いる。従って、e2 = 1である。交換エネルギーEx (ρ)について、これは(1)式より、

Ex(ρ) = ∫ρεx(ρ)dr ・・・(1)

と表現でき、交換エネルギー密度εx(ρ)は、

εx(ρ) = Aρ1/3 ・・・(2)

ということは既に(2)式で示した。Aは適当な係数であり、局所密度近似下に おける一様な電子ガスでの係数Aは、

A = -3/4 (3/π)1/3 ・・・(8)

である。この場合、Aの値（絶対値）はおおよそ、0.7386...となる。上式の 導出に関しては、バンド計算関連の参考文献、 書籍等を参考にして欲しい。交換エネルギー密度εxは、rsパ ラメータ、フェルミ波数それぞれで表現すると、

εx = -3/4 (3/π)1/3 ρ1/3 ・・・(9a)

= - 0.4582 / rs ・・・(9a)

= - (3/4) kF ・・・(9c)

となる。この時、

rs = (3/4)1/3 (1/π)1/3 (1/ρ)1/3 ・・・(10)

kF = (3 π2 ρ)1/3 ・・・(11)

である。交換ポテンシャルは、

Vx(r) = d(ρεx(ρ))/dρ ・・・(12)

から求められるので、

Vx = d(ρεx(ρ))/dρ = - d( (3/4)(3/π)1/3 ρ4/3)/dρ

= - (3/4)(4/3)(3/π)1/3 ρ1/3

= - (3/π)1/3 ρ1/3(r) ・・・(13)

となる。この時、ρ1/3に付く係数（絶対値）(3/π) 1/3はおおよそ、0.98475...（先のA値の4/3倍）となる。以上は、 ハートレー原子単位で考えている。上 記の式には表示されていないが、ハートレー原子単位では、e2 = 1であるが、リュードベリ単位では、e2 = 2である。また、スレー ターによる交換ポテンシャルは上記交換ポテンシャルの3/2倍となる。Xα近似 では係数αがつき（Vxαとする）、

Vxα = αVx,sl = (3/2)αVx ・・・(14)

となる。Vx,slはスレーターの交換ポテンシャル(= 3/2 Vx)、αは通常、 (2/3)〜1までの範囲をとる。α = 2/3がVxであり、α = 1でVx,slとなる。
スピン偏極を考える場合は、交換エネルギーεxに21/3が付く。 ほとんどの式で省略しているが、Vx = Vx(ρ) =Vx(ρ(r)) = Vx(r)、εx = ε x(ρ) = εx(ρ(r)) = εx(r)、ρ = ρ(r)である。
尚、以上の議論は、分子のおも ちゃ箱の掲示板での議論が元となっている（←係数の問題は重要であり、 覚書として記す必要、意義があると考えた）。

局所密度近似 [目次]

（３／１、２００５）電子ガスにおける交換相関エネルギーExc（交換相関 項、交換相関エネルギー項とも言う）を、

Exc[ρ] = ∫ρεxc(ρ)dr ・・・(15)

とするのが局所密度近似である。ρ は電荷密度（ρ = ρ(r)、rは座標で、r = r）である。上式左辺の括 弧が[]となっていて、密度汎関数法に おいて厳密に求められるべきものであるが、右辺側は近似的に、ρεxc(ρ)と いうように局所的なρ(r)に関しての体積積分で求まるとしている（左辺→右 辺における、[]→()は重要）。勿論、厳密なExc[ρ]を現時点で求めることは 不可能である。求める式も方法も手段も全く未解決である。局所密度近似は、 上式右辺のようにして近似的にExcを求めるものである。交換相関エネルギー 密度εxc（電子１個当たりの交換相関エネルギーと表現する場合もある）の表 式も、近似的に求められたものが使用されている。電子ガスでの低密度、高密 度からの内挿（または外挿）による方法や、RPAを用いて求めたもの、モンテ カルロ法から求めたものなど、εxcにはいくつかの表式が存在する。スピンを考慮しない範囲で は、それぞれの表式で大きな差はみられないが、場合によりあまり良い結果を 与えないこともある。また、磁性が絡む場合 も注意が必要（調査中）。
交換相関エネルギー密度εxc(ρ)が電子１個分に相当すること の（あまり厳密でないかもしれない）説明。電荷密度ρが一様であるとすると、

ρ(r) = ρ = N/Ω ・・・(16)

となる。Nは系（単位胞内）の全電子数。Ωは系（単位胞）の体積である。 こうすると、(15)式右辺は、

∫ρεxc(ρ)dr = ∫(N/Ω)εxc(N/Ω)dr

= (N/Ω)εxc(N/Ω)∫dr = (N/Ω)εxc(N/Ω)Ω = Nεxc(N/Ω) ・・・(17)

となる。以上から、

εxc(N/Ω) = Exc/N ・・・(18)

であり、交換相関エネルギー密度は電子１個分の交換相関エネルギーとなっ ている。以上は、ρが一様（一定）であることを前提としているが、実際の局 所密度近似も電荷密度の空間的な変化が緩やかであることを前提とした近似で ある。空間上で激しく電荷密度が変化するとしたら、局所的な記述が破綻する ことは容易に想像できる。
局所密度近似では不完全な形で相関項を取り入れるため、ハートリー・フォッ ク近似の段階では厳密に相殺される、電子自身が自分と相互作用する自己相互作用部分が残ってしまう。このため無限遠 で感じる電子のポテンシャル（一電子ポテンシャル）が、-1/rにならずに指数 関数的に減衰してしまう。関連用語：自己 相互作用補正(SIC)

（相関項の由来）
系の全エネルギーEtotが厳密に得られているとして、ハートリー・フォック 近似で得られる当該する系の全エネルギーを、EtotHFとする。相 関項Ecは、

Ec = Etot - EtotHF ・・・(19)

として表現される。つまり厳密な全エネルギーの値とハートリー・フォック 近似によるものとの差の部分が相関項（相関エネルギー）となる。相関項は勿 論、（クーロン項＋交換項以外の）電子間の相互作用を記述する項である。 （関連）相互作用していない系での運動エネルギー項←相互作用する部分は、 全て相関項に押し込める。相互作用していない電子の運動エネルギーなので、 ビリアル定理での扱いにも気を付ける必要 がある（調査中）。
（注意）相関項の考え方は、必ずしも上 記の定義の仕方に縛られない（上記以外の説明、定義のされ方がなされる場合 がある。調査中）。

（局所とは）
Exc = ∫ρ(r)εxc(ρ(r))drと、局所的な座標r(= r)のみに依る形 の状態を、”局所”と言う（数学的、物理的には不確かな説明かもしれない）。 もし、Excが、

Exc = ∫f(r,r')drdr' ・・・(20)

のような形で表現される場合、最早座標rのみによる局所的な記述は出来な くなる。これは”非局所”なものとなっている。つまりその場（の座標 r）だけでExを決めることは出来なくなっている。局所、非局所と言う 用語は上記と異なる意味合いで使われることがある（調査中）。関連用語：局 所擬ポテンシャル、非局所擬ポテンシャル

（密度汎関数法〔または密度汎関数理論〕は変分問題）
変分パラメーターは（電子の）電荷密度。密度汎関数法は、絶対零度での変 分問題。従って全エネルギー（＝内部エネルギー）の電荷密度に関しての変分 となる（Kohn, Shamの論文）。有限温度 への拡張は、Merminによるものが最初（論文 ）である。この場合、全エネルギーでなく自由エネルギーに対する変分と なる（変分パラメーターは電荷密度）。ただ、バンド計算で具体的にどのよう な処方を用いればよいかに関しての言及はされていない。実際に有限温度をバ ンド計算で扱うのは、大変難しい。一つの処方としてはバンドにフェルミ分布 関数による広がり（ブローデニングと同じ方法）を持たせて、その広がりに関 するパラメーターを変分の対象として扱うことが行なわれている（有限温度ＬＤＡ、調査中）。

関連ページ：[密度汎関数法、局所密度近似関連 情報]

バンドギャップの過小評価の原因 [目次]

（４／８、２００５）これに関しては、Lecture Notes in Physics, A Primer in Density Functional Theoryの5. R. W. Godby and P. Garcia-Gonzalez, "Density Functional Theories and Self-energy Approach"の196〜199頁（5.3.1節）の部分、福本先生に「よるレビュー（豊田 中央研究所 R & Dレビュー、Vol. 25, No. 4）、浜田先生、大西先生による” 日本物理学会誌”での解説（「固体の電子状態」、第４０巻、第１１号（１９ ８５）、８４２頁の部分）、寺倉先生、浜田先生による雑誌”固体物理”での 解説（「バンド計算法の最近の発展（II） - 密度汎関数法 -」、第２０巻、 第９号（１９８５）、７０６頁の部分〔第５節〕）、里子先生等の著書（「密度汎関数法とその応用‐分子・ク ラスターの電子状態」）などに関連する記述がある。更に参考文献として、

[1] J. P. Perdew, R. G. Parr, M. Levy and J. L. Balduz, Jr., Phys. Rev. Lett., Vol. 49, No. 23, 1691(1982)
[2] L. J. Sham amd M. Schlüter, Phys. Rev. Lett., Vol. 51, No. 20, 1888(1983)

などが挙げられる（福本先生のレビューでの参考文献も有用）。

以下、スピンは考えず、バンドギャップの存在する半導体、絶縁体を対象と する。スピンを考えないので、一つの準位（または軌道、バンドと言ってもよいだろう）に収容でき る電子の数は２である。まず、電子親和力をA、及びイオン化エネルギー（イ オン化ポテンシャル）をIとすると、

A = Etot,N - Etot,N+1 ・・・(21)

I = Etot,N-1 - Etot,N ・・・(22)

である。Etot,NはN個の電子からなる系の基底状態での全エネル ギーを意味する。上記の電子親和力とイオン化エネルギーから、バンドギャッ プEgapは、

Egap = I - A ・・・(23)

と表される。このEgapが、実際の実験等で観測されるバンドギャッ プに相当する。但し、実験手法などによりバンドギャップの解釈は微妙に異な る場合がある。実験が一体何を観測しているのか、そしてその時のバンドギャッ プの定義なり解釈なりの議論は大変重要であるが、ここではその詳細を述べる だけの力量が筆者にないため、これ以上の言及は避けることとする。
一方、このEgapをKohn-Shamの密度汎関数法で表現すると、

Egap = εKS,N+1(N+1) - εKS,N(N) ・・・(24)

となる。上の式で、εKS,N+1(N+1) は、N+1電子系でのN+1個目 の準位（軌道、バンド）、つまり最高占有準位（バンド、軌道またはコー ン・シャム軌道）である。εKS,N(N)は、N電子系でのN個目の 準位（軌道、バンド）である。これはN電子系での最高占有準位（軌道、バン ド）になっている。そして、実際のコーン・シャム形式での密度汎関数法に基 づくバンド計算で得られるバンドギャップは、

εKS,gap = εKS,N(N+1) - εKS,N(N) ・・・(25)

となる。εKS,N(N+1)は、N電子系でのN+1個目の準位（軌道、バ ンド）、つまり空の準位（軌道、バンド）である。以上から、 EgapとεKS,gapとの関係が導かれ、それは、

Egap = (εKS,N(N+1) - εKS,N(N)) + (εKS,N+1(N+1) - εKS,N(N+1))
= εKS,gap + Δxc ・・・(26)

となる。ここで出てくるΔxcは、

Δxc = εKS,N+1(N+1) - εKS,N(N+1) ・・・(27)

であり、これはN+1電子系でのN+1個目の最高占有準位（軌道、バンド）と、 N電子系でのN+1個目の空（非占有）の準位（軌道、バンド）との差となってい る。N電子系に対してバンド計算を行なった場合、電荷中性の要請からそもそ も、N電子系に１個だけ電子を加える（=N+1個系）ような計算が出来ない。周 期的境界条件下での単位胞（またはスーパーセル）で、電子を１個だけ加えれ ば、エバルト項などが発散する。これの対処としてのジェリウムモデルの導入は、元のN電子系との整合性 が崩れ、上記（式26など）にある全エネルギー同士の比較が意味をなさなくな る。
Δxcの起源に関しては、これをN-δからN+δ（Nは整数：電子の 数、δは微小量〔＝無限小〕）における交換相関ポテンシャル部分の不連続な とびによるとする説明が、先の文献[1][2]等でなされている。つまり、

Δxc = vxc(N+δ) - vxc(N-δ) ・・・(28)

である。vxcは、交換相関ポテンシャル。また、δExc/δρ=vxcより、

Δxc = δExc[ρ]/δρ|N+1 - δExc[ρ]/δρ|N + O(1/N) ・・・(29)

となる。ここでNは十分大きな量と考えている（N ≫ 0）。Nが十分大きけれ ば、誤差の部分Oは無視できることとなる。式２８にも本当は誤差部分 が存在する。式２８→式２９に関しては、筆者自身まだ十分理解できていない 部分あり（←調査、勉強中→）。Nが十分大きい（＝無限大）ということと、 δが無限小であるということは、いずれもN+1（N → ∞）、N+δ（δ → 1/∞） での軌道緩和効果が事実上無視できることを意味している。
Δxcが正の有限な量であれば、通常のバンド計算によるバンド ギャップは常に過小評価していることが示せる。交換相関ポテンシャルの非連 続的なとびによる過小評価が、全体としてのバンドギャップのバンド計算によ る誤差の８割方を占めるということである。残り２割の中には、軌道緩和の効 果などが含まれるものと考えられる。系が金属だった場合、Δxc は、ゼロか負の値である必要がある（正ならギャップが生じ得るし、負の値だっ たとしてもそう大きな値にはなり得ない）。但し、Δxc > 0は、 今の段階で筆者にとっては自明ではない（←調査、勉強中）。更に、第一原理 的かつ簡便にΔxcを定量的に評価することは困難である。これを 解決するには、局所密度近似を越える試み における自己相互作用補正や、GW近似などを考える必要がある。

【おまけ】ヤナックの定理より、i番目 のコーン・シャム軌道をεiとし、その占有数をfi（0 < fi < 1とする）となる。この時、イオン化エネルギーIは、

I = Etot,N-1|fi = 0 - Etot,N|fi = 1 = - ∫0 → 1 εi(f) df ・・・(30)

となる。積分部分、∫0 → 1 εi(f) dfにおいて、

∫0 → 1 εi(f) df ≒ εi(0.5) ・・・(31)

とし、I = -εi(0.5)とするのがSlaterの遷移状態理論（近似） である。３１式導出に関しては小口先生の著書「バンド計算」（内田老鶴圃刊） ３．４節参照。"0.5"は占有数が半分であることを意味する。ここで、一つの KS軌道（準位、バンド）の収容電子数（１か２）のことは深く考えないことに する。

ＧＧＡ編 [目次]

現在（６／１、１９９９）、Perdew先生、Burke先生、Ernzerhof先生による ＰＢＥコードを入手、revpe_d.fに導入する作業を遂行中です。

PBEに関しての[参考文献]、Perdew先 生の[ウェブページ]（ＰＢＥコード入 手に関して深く感謝する次第です）。

（７／１６、２００１）現在のPBEコードのある[ところ]（Rutgers大学のKieron Burke 先生のサイト、”past projects”のページにPBEコードへのリンクがある。使 用方法、諸注意は当該ページ、データ等のドキュメント参照のこと）。

ＧＧＡ導入に関して、現在プログラムへの導入作業はほぼ終了し、バルク系、 表面系に対するテストを行なっています。これらの顛末については、現在執筆 準備中（＊）であります。

（２／１、２０００）現在、仮想的なＢＣＣ構造（立方体のユニットセル内 に２個原子が存在）でのＡｌのテスト計算を行ない、原子位置を平衡位置から 少しだけずらせる（いくつかの場合で検証）ことによって生じる力について [検証]しました。
検証の結果、力は（一部の場合、ｋ点数を十分にとる必要あり）十分な精度 で正しく得られていることを確認しました。この場合、力は線形近似（調和近 似つまり”ばね”）の範囲内で十分記述され、全エネルギーの位置変位による 変化分＝１／２×Ｆ×Δxとなります（Ｆ：力、Δx変位長さ）。Ｆ＝ｋΔxな ので変位の長さを１／１０にするとエネルギーの変化量は１／１００になりま す。

少し状況報告（６／２５、１９９９）。ＧＧＡでは電荷密度に対して微分を 行なうため、特に表面系のようなスーパーセル＋スラブ近似を行なっている場 合、真空層での電荷密度がほとんどゼロかそれに近い部分での微分対しての寄 与について考えなければなりません。バルク系の計算では、さして計算に問題 は生じなかったのですが、表面系では、全エネルギーも原子間に働く力も、バ ルク系で扱った手順、条件と同じではまともに計算できないことが判明しまし た。

（２／１９、２００３）現在も少しずつバグ等を潰しつつ作業はゆっくりと 進行中。最近、グラファイトの格子定数をＧＧＡではまともに求められないと いう話を聞いて、こちらでも試してみる（ＧＧＡの関数形は、ＰＢＥ）。結果 としてＧＧＡでは、ｃ軸方向の格子定数が求められなかった（バルクでいるよ り、単独のグラファイトシート〔単層膜〕の方がより安定となってしまう）。
ＬＤＡでは、格子定数は求まるが、採用する関数形によって合う場合と、全 く合わない（Wigner形式では、ｃ軸の格子定数が数％以上伸び目に出る）場合 があり、また精度を上げよう（エネルギーカットオフを上げる）とすると、む しろ結果が悪くなってしまう（筆者の計算での話）。
ＧＧＡでは、van der Waals結合を正 確に表現することが特に困難で（勿論、ＬＤＡでも困難）、他に、窒素の分子 性結晶、グラファイトの面間距離（c軸の格子定数）など、van der Waals力に よって結合している物質の格子定数を正確に求められないことが分かっている。
グラファイト類似の面間の結合が弱いと推察される系（炭素-ホウ素化合物） でも、やはりＧＧＡで平衡面間隔が求められない事例あり。局所密度近似では 最小値、つまり平衡面間隔が存在する。ＧＧＡでは面間隔が20 a.u.以上とい う値にしてもどんどん全エネルギーが下がっていく。→つまり、バルクでなく 単独の一平面である方がより安定と結論される（７／１５、２００４）。
（１１／５、２００３）ＧＧＡ計算に関しては若干の実践的進展あり。HBC （LiBC構造と同じ結晶構造、詳細はLiBC[ページ ]参照）に関して、PBE形式でのＧＧＡ計算結果をIUMRS-ICAM2003にて発表。 まだ、ＧＧＡの基での圧力の導出は出来ていないので、格子定数を少しずつ変 えながら全エネルギーの最小値を探索するという原始的手段を用いている。結 果は、それなりに正しい模様、詳細はIUMRS-ICAM2003のプロシーディングス発 刊（K. Kobayashi, M. Arai and T. Sasaki, "Lattice Anomalies of MBC (M = H, Li, Na) Under Anisotropic Compression", Trans. MRS-J, Vol. 29, No. 8, 3799 - 3802 (2004)[IUMRS-ICAM2003]）以降 に公開。

以上（ＧＧＡ導入）に関しての対処方法について今後も述べる予定（現在も なお調査、検討中、７／１６、２０００→２／１９、２００３→１１／５、２ ００３→７／１５、２００４）。

【参考文献】
[1][PBE]J. P. Perdew, K. Burke and M. Ernzerhof: Phys. Rev. Lett. 77 (1996) 3865[PBE].
[2]J. P. Perdew and W. Yue, Phys. Rev. B33, 8800(1986)[(24)式： GGA exchange potential]
[3]J. P. Perdew, M. Ernzerhof, A. Zupan and K. Burke, J. Chem. Phys. 108, 1522(1998)[Appendix A: GGA correlation potential]
[4]J. A. White and D. M. Bird, Phys. Rev. B50, 4954(1994)[Implementation of gradient-corrected exchange-correlation potentials]

[GGAにおけるストレス関連]
[1]A. D. Corso and R. Resta, Phys. Rev. B50, 4327(1994)
[2]L. C. Balbas, J. L. Martins and J. M. Soler, Phys. Rev. B64, 165110(2001)
[3]N. Nagasako and T. Oguchi, Journal of the Physical Society of Japan, Vol. 82, No. 4, 044701(2013)[Stress formulation][All-electron FLAPW][Accuracy check][GGA]
[4]L. C. Balbas, J. L. Martins and J. M. Soler, Phys. Rev. B64, 165110(2001)
[5]N. A. W. Holzwarth先生(Department of Physics, Wake Forest University)のところに ある、"Notes on GGA implementation"（PDFによるノート）は大変有用。

電荷密度の微分[目次]

一般化された密度勾配近似という名称からも分かるように、ＧＧＡでは電荷 密度ρ(r)の微分（▽ρやΔρなど）が必要となります。数値微分を使って電 荷密度の勾配を求めることも可能ですが、ＧＧＡでは普通フーリエ変換（実際 には高速フーリエ変換〔FFT〕）を使用して勾配を求めます。これは通常、数 値微分よりフーリエ変換を利用する手法の方が、精度的、速度的に有効なため です。この手法は、実空間での電荷密度ρ(r)の逆フーリエ変換（←以下の式 〔G → r〕を逆変換と定義）は、

ρ(r) = ΣGρ(G)exp(iGr)

となります。G：逆格子ベクトル、r：実格子ベクトル、i：虚 数。ρ(G)：実空間での電荷密度をフーリエ変換した（逆格子空間上で の）電荷密度。これに対する勾配の絶対値｜▽ρ(r)｜は、

｜▽ρ(r)｜ = ｜iΣGGρ(G)exp(iGr)｜

となります。上式でのexpの肩の部分でrに関して微分します。 rに依るのはこの肩の部分のみです。つまり、iGρ(G) に対する逆フーリエ変換をすれば、｜▽ρ(r)｜を求めることが出来ます。同 じ要領で、Δρ（ = ▽2ρ(r)）なども求めることができます。

スピン軌道相互作用（失敗）編 [目次]

希土類金属擬ポテンシャルにアプローチするためには、スピン軌道相互作用 を考慮しなければならないということで、摂動形式で導入することを計画しま したが、希土類金属くらい重い原子での計算ではスピン軌道相互作用は摂動的 な取り扱いでは無理ということが判明（これは、金材研での鈴木修吾先生のセ ミナーで強く認識しました）しました。従って、摂動形式での導入だけでは駄 目であることが分かり、現在作業、計画共頓座しています（６／１、１９９９）。

比較的軽い原子である、In、Sb、Snなどでは、スピン軌道相互作用によって バンド構造におけるバンドの分裂が起こりますが、格子定数や体積弾性率など は、スピン軌道相互作用による影響をほとんど受けません。実際、我々の扱っ た、InSb/Sn界面の計算（神奈川工 科大学の山本氏との共同研究）ではスピン軌道相互作用は考慮していませんが、 格子定数、体積弾性率は実験値と十分一致（InはＰＣＣを考慮しても尚２％以上の格子定数 のずれがありましたが、InSbとしてなら、実験値と良く一致）する結果が得ら れています。
このようなIn、Sb、Sn辺りでは摂動でスピン軌道相互作用を導入することに より、より精度の高いバンド構造を得ることが可能と思われます。またこれに よる格子パラメーターへの寄与はほとんどないと思われます。これだけで断定 はできませんが、構造最適化においても、スピン軌道相互作用は（バンド分裂 の影響が出てくるような場合を除いて）ほとんど影響しないだろうと思われま す（確証なし、調査中）。

ただ、他にも幾つか問題があり、実際にバンド計算する段階まで全く達して いません。特に、スピン軌道相互作用用擬ポテンシャルを作成する目処が全く 立っていません（これでは計算にとりかかれるはずもない）。

擬ポテンシャル＋平面波によるスピン軌道相互作用の摂動による表式は、次 節で示してある参考文献から得ることが可能ですが、これをさらに具体的に表 したものを示します（png画像、30KBほど です。方向を修正し見易くしました。６／２１、１９９９）。
（注意）この表式は、計算による確認等を全く行なって いません。従って、本当に正しいかどうかは不明です（筆者及び物質・材料研 究機構は一切責任を負えません）。

スピン軌道相互作用関連文献[ 目次]

筆者が現在、把握しているものを以下に示します（当然、全てを網羅してい る訳ではありません。擬ポテンシャルが関連しているものが中心です）。 以上の他に、（旧）電総研（産総研）の織田先生の博士論文（高圧化でのヨ ウ素の計算で、スピン軌道相互作用を考慮した計算も行なっている）も参考に なりました。擬ポテンシャルの計算は摂動によるものです。

[先頭][総目次 ][最初に戻る][使ってみよう編][バンド屋さんマップ][密度汎関数法、局所密度近似関連情報]