やさしい（？）バンド計算プログラムの作り方（各論） [Top]

目次

各論１

各論２（MSD,FORCE,CHAVER,STRNLの説明）

サブルーチンＭＳＤ

サブルーチンＦＯＲＣＥ

サブルーチンＣＨＡＶＥＲ

サブルーチンＳＴＲＮＬ

各論３（その他のサブルーチン１）

各論４（その他のサブルーチン２）

【注意】

各論１

各論２（MSD,FORCE,CHAVER,STRNLの説明）

1. [ まず最初に、コモン変数呼び込みのためのインクルード文、パラメーター変 数、局所的に使う配列（含む作業用配列）の定義、ＶＰＰ用並列化定義（ｋ点 に関する並列化）などが行なわれます。]
   
   
2. [ その後、初期変数値の設定（含むゼロ初期化）を行ない。]
3. [特別な 条件（ストレスによるユニットセルの最適化）の時に非局所擬ポテンシャル部 分の動径方向（実空間）から逆空間へのフーリエ変換を行なうサブルーチンＫ ＢＩＮＴと力の計算（非局所擬ポテンシャル部分、簡易版なので実際の力は計 算していない）サブルーチンＦＯＲＺＦＢを呼び出しています。]
   
   
4. [ そして次に、運動エネルギー、ハートリー、交換相関、局所・非局所擬ポテ ンシャル、各項に対応する行列要素を作ります。]
   
5. [ 波動関数は最急降下法で更新します。但し、収束を加速するため、波動関数 に対する時間に関しての１次微分方程式を解析的に解き、時間刻みをより大き くとれるようにしています（詳しくは筆者のＤ論）][（またはA. R. Williams and J. Solerの文献参照）]
   
   
6. [ 更新された波動関数を直交化（グラムシュミットの方法）します。その後エ ネルギー固有値（固有エネルギー、バンドレベル、エネルギー準位などいろい ろ言い方があります）を求めます。]
7. [求めた エネルギー固有値は必要によって並べ替えを行ないます。]
   
   

1. [ ま ず最初に、ＰＡＣＶＰＰをインクルードし、局所（作業用）変数を定義する。 ]
2. [ そ こでは、ｋ点に関する並列化の設定も行なわれています。]
3. [ 次 に作業用変数等をゼロ初期化する。]
4. [ そ して、ＭＳＤでも使う非局所擬ポテンシャル用行列要素の計算を行ない。 ]
5. [ 最 後に、非局所擬ポテンシャルの力の計算を行ないます。]

1. [ まず最初に、ＰＡＣＶＰＰをインクルードし、局所（作業用）変数を定義する。 ]
2. [ そこでは、ｋ点に関する並列化の設定も行なわれています。]
3. [ 次に作業用変数等をゼロ初期化する。]
4. [ そして、逆空間での波動関数をＦＦＴ用の３次元複素配列に（マッピングしな がら）格納します。]
5. [ これをＦＦＴで逆フーリエ変換（逆空間--->実空間をフーリエ逆変換と定 義しています）して実空間の波動関数（３次元複素配列）を作ります。]
6. [ 次に、畳み込みを行ないます。これは単純に先に求めた実空間波動関数同士を かけているだけです。これによって実空間電荷密度が得られます。この時、バ ンド（エネルギーレベル）とｋ点に関してフェルミ面までの和をとります。こ れにより系の総電荷密度が求まります。]
7. [ そして最後に、求めた実空間電荷密度をフーリエ変換して、逆空間での電荷密 度を求めます。]

1. [ ま ず最初に、ＰＡＣＶＰＰをインクルードし、局所（作業用）変数を定義する。 ]
2. [ そ こでは、ｋ点に関する並列化の設定も行なわれています。]
3. [ 次 に作業用変数等をゼロ初期化する。]
4. [ そ して、各項毎のストレスの計算を行ないます。最初に作業用変数の初期化を行 ない、表式通りに計算を行ない、各配列に値を代入していきます。計算１ ]
5. [ 計 算２ ]
6. [ 計 算３ ]
7. [ 計 算４ ]

各論３

KPMWBZ

これは、ブリュアンゾーン内のｋ点座標の設定を行なうサブルーチ ンで、対称性は考慮していません。従って、第一ブリュアンゾーン全体を等メッ シュで分割したｋ点座標を作ります。この時、ｋ点座標全体を原点（Γ点）を 基準として、ずらすことも可能です（深いピンク色になっている注釈行）。但 し、ずらし方は、慎重に行なわないと思わぬ失敗の原因になります。
またKPMSF は２次元（つまり表面用）の場合でのｋ点を作成します。

KBINT

これは非局所擬ポテンシャル部分の実空間（動径方向）から逆空間 へのフーリエ変換を行なうサブルーチンです。基本的にはKBMATと同じです。
KBINTでは後でス トレス計算で使用するための逆空間用非局所部分を配列変数に格納している部 分があります（深いピンク色で示した部分）。

FORZFB

これはサブルーチンFORCE の簡略版です。つまり、電子状態の計算のみを行なう場合などでは、力の計算 を行なう必要がないので、それらの部分を取り除いた部分からなるのが、この FORZFBです。小規模なルーチンなので、並列化部分の参考例としても良いかと 思います。

STRESS

非局所擬ポテンシャル以外の部分のストレスの計算を行なうサブルー チンです。ストレスの計算は、非常に複雑なので、詳しい説明は避けます。基 本的に、Nielsen and Martinの方法に従っ て、ストレスの各項の計算を行なっています。

FORLOC

力の非局所擬ポテンシャル以外の項（局所擬ポテンシャル項、エバ ルト項）からの寄与部分の計算（エバルト項は別のルーチン〔非公開〕で計算 しています）と、取りまとめ（非局所項も加えて、最終的な力を求めています） をしているサブルーチン。

（１０／１、１９９９）現在、このルーチンでは ZVXCPCという変数を使って力の計算を行なっている部分がありますが、この部 分は不必要であることが判明しました。従って、現在当該部分を注釈行として 扱っています。

MD

FORCE、 FORLOC、 STRESS、 STRNL で求めた力、ストレスをもとに分子動力学計算を行なうサブルーチン。原子は ベレのアルゴリズムで更新されます。また、ユニットセル自身も、ストレスを 使って任意の外部圧力条件下のもとに最適化することが可能です。ユニットセ ルそのものの動力学はベレまたは最急降下法で行なわれます。この部分は後か ら付け足した形になっていて、ルーチンとしての構成はあまり良いものではあ りません。
このサブルーチンは、計算全体の１％以下し か占めません（だからといって疎かにしてはいけません）。

ENERGY

（９／３０、１９９９）このルーチンは各サブルーチンで計算され た結果をもとに、全エネルギーを最終的に求めるサブルーチンです。このルー チンに最近、問題があることが分かりました。それは、変数EPCに関してのものです。
サブルーチンENERGYでは最終的な全エネ ルギーの値ETOTALに変数EPCが引かれています。このEPCはサブルーチンPCCで求 められていますが、系のセル（体積）に依らない単なる定数であり、全エネル ギーにとっては余計なパラメーターであることが判明しています。EPCは系の 体積（格子定数）に全く依存しないので、エネルギーの基準がEPC分シフトし ただけであり、同じ原子構成からなる系での全エネルギーの比較という点では 何ら問題ありません。但し、凝集エネルギーを求める場合（系同士の構成原子 の数が異なる）には問題を生じる可能性があります。
もし、問題があると思われる場合は、当該部分、ETOTAL = ETOTAL-EPCを注釈行にして計算しないようにして下さい。

CONV2

電子状態計算での収束判定をするサブルーチンです。またイタレー ションでの古い電荷密度と新しい電荷密度の混合も行なっています。
混合の仕方には他にも、ブロイデン法 や準ニュートン法によるものなど幾つか存在します。対角化を使用しないこの 手の方法では、波動関数の更新（MSD） の段階で、古い波動関数（１イタレーション前、最急降下法の場合）を新しい 波動関数に混ぜているので、これは電荷混合と同じ効果を既にもたらしていま す。しかし実際には、電荷の混合も合わせて行なうと更に収束を速める場合が 多いです（筆者はそうでない場合を経験したことがありません）。

各論４

XCFFT

これは交換相関ポテンシャル、エネルギーを計算するためのサブルー チンです。具体的には、逆空間の（CHAVER で作った）電荷密度を高速フーリエ逆変換を使って実空間電荷密度にします。 本計算ではWignerの方法を使って、この 実空間電荷密度から交換相関ポテンシャル、エネルギーを求めます。
また、ここでＰＣＣの考慮が必要な場 合は、その計算も行ないます。
そして、求めた実空間での交換相関ポテンシャル、エネルギーを再び高速フー リエ変換を使って、逆空間のものに戻します。
因みに、このサブルーチンではポテンシャルとエネルギーは同時には計算さ れません。（ポテンシャル用、エネルギー用と場合分け〔当然IF文で〕されて います。）

XSTPC

これは交換相関ポテンシャル、エネルギーのストレス版です。特に、 部分内殻補正ＰＣＣがストレス計算で考 慮されています。
部分内殻補正に対するストレスの表式（pdf形 式、159 kb、Appendix参照、CMSMD'96の原稿pdf）。
但し、実を言うとこのサブルーチンには問題があります。深刻なほどではな いのですが、完全に無視できるほどではありません（厳密にはバグとも言えな い）。もう少し新しいバージョンのプログラム（磁性考慮版）では修正されて いますが、本バージョンでは取り敢えずこのまま公開してしまいます。
また、ＰＣＣを使わない場合は、この問題 は顕在化しません。

DIAGON

いよいよ最後（でなくなった。１０／１３）に残った、各論４で最も 重要なサブルーチンの説明を行ないます。このルーチンは本計算プログラムで は第一回目のイタレーションでのみ呼び出されます。
つまり、カー・パリネロ的な計算のためには、初期値としての固有ベクトル （つまり波動関数）と固有値（固有エネルギー）が必要です。そのために最初 のイタレーションで、従来から使われてきた対角化による方法を用いて、初期 固有ベクトル、固有値を求めます。これを第２回目以降の逐次更新型の計算イ タレーションで利用する訳です。
この対角化の計算部分では、最初にマトリックス要素（運動エネルギー項、 ローカル項〔局所擬ポテンシャル項、ハートリー項、交換相関項〕、非局所擬 ポテンシャル項、その他）を計算します。これを対角化することによって、固 有ベクトル、固有値を求めます。

ここでの問題が対角化を行なうサブルーチンですが、 このDIAGON ではCHOBSD、HZES1M、DHEIG2（後の２つは注釈になっています）を呼び出して います。これらはそれぞれNUMPAC（科学技術計算ライブラリ、元は名古屋大学 の研究者が開発したもの、現在は富士通が管理）、日立のメインフレーム用の 科学技術計算ライブラリ、SSL2（富士通のメインフレーム用の科学技術計算ラ イブラリ）の各ライブラリに対応しています。これらは市販されているもので、 当該企業に著作権、所有権があります。従って、これらを公開することはできません。

筆者が一番最初に計算した頃（１９８９年当時）は、ＳＳＬ２のDHEIG（エ ルミート行列の対角化）を使っていました。このサブルーチン（DHEIG）は、 対角化すべき２次元配列の格納方法が特殊だったため、DIAGONでのマトリック ス要素の計算方法が大分複雑になってしまっています。当然、DHEIG以外のルー チンを使う場合は、格納方法が異なるのですが、最初の方法をそのまま踏襲す るような形でルーチンに修正を加えています。従って、このDIAGONサブルーチ ンはあまり可読性が良いと言えません。

FERMI

（１０／１３、１９９７）追加としてサブルーチンFERMI の説明を行ないます。これは、東京理科大学の浜田先生が作成したものです。 これは、バンドが完全に詰まり切っていない場合（部分的にバンドに電子が占 有される状態、つまり金属状態）で、ｋ点とバンド（バンドインデックス）の 足し上げを行なうための簡易方法の一 つです。
具体的な方法は、バンドを二つの放物線関数を張り合わせた幅のある形で表 現し、それが部分的に占有され、金属的なバンド構造の場合でのｋ点、バンド 毎の足し上げを可能とします。

幅の与え方によって、全エネルギーの値が異なるので、系の全エネルギーを 比較する場合は、同じ幅で統一しておかないと比 較する意味が無くなってしまうので注意が必 要です。（本当は、幅をゼロにした場合への外挿なども考慮する必要がありま すが、筆者はそこまではやっていません。ただ極めて厳密な計算が必要でなけ れば、外挿までする必要はないと考えられます。）

尚、サブルーチンWIDTH2 は、幅を持たせたバンドでの、電子の占有数を計算（放物関数を使っているの で厳密に求められる）しています。