How to make the FPMD code easily(?)

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

目次: 準備（バンド計算プログラムを作成するのに必要な予備知識、道具）; バンド計算プログラム作成前に注意しておくべきこと; 実際のコーディング; プログラム全体の構成; プログラムのメイン部分; メインルーチンの内容（各論）

【注意】：今は、”作る”より”使う”[時代 ]（１１／２６、２００２）

本当は、バンド計算（第一原理電子状態計算）プログラムを１つでも作ることは、大変な時間と労力と知識を必要とし、決してやさしいものではありません。

ただそれを言ってしまっただけで話しを完結させてしまうと身もふたもないので、可能な限りわかりやすく説明を行ないたいと思います。
このようなものを書く理由の一つに、自分の持っているプログラムの整備を行ないたいという意味も込められています。ご存知のように擬ポテンシャルデータベースＮＣＰＳ９５を構築中（現在、NCPS2Kを[配布中]：１１／２６、２００２）ですが、これはこのポテンシャルデータを利用できるプログラムがないことには、このデータベースは何の役にも立たないことになります。

こういう意味では、自分のプログラムも一緒に公開すれば良いのですが、なかなかそうは行かないわけです。その理由として、

（１）プログラムに汎用性が（全く）ない。
（２）プログラムの可読性が非常に悪い。
（３）プログラム内に他の共同研究者との共作の部分があり、著作権の問題が絡んでくる。
（いろいろと公開するための事前の手続きが面倒）

特に、（１）、（２）の問題が深刻で、公開することを想定していないでバンド計算プログラム（第一原理分子動力学法）を作ったので、非常に難解で、動かしにくいものになっています。またプログラムに関してのドキュメント、マニュアル等が全くないことも、公開を妨げる原因の一つになっています。

そこで、公開することを想定（最近一部公開しました、３／１８）して、プログラムに関してのドキュメントを充実させるために、このような企画を立てました。これによって、自分のプログラム作成に関するノウハウを公開し、かつこれによってドキュメントも充実させることができます。

準備

バンド計算プログラムを作成するのに必要な予備知識、道具

まず最低線、量子力学の知識が必要です。成積はどうでもよい（というわけでもないが）ですが、どのような参考書や、どの本のどこらへんの部分を調べれば良いか位はわかる程度の知識は必要です。特に縮退とか固有値問題、固有ベクトル（線形代数の知識でもある）、バンド理論（物性論の知識でもある）、変分原理、ハートリー・フォック近似（平均場近似）、摂動等の用語を知っているか、すぐに調べられる（そして理解できる）程度でなければなりません。

量子力学以外にも、固体物理の知識も必要です。ブロッホの定理やバンド理論、逆格子空間の概念、ブリュアンゾーン、結晶対称性や群論など結晶学的な知識などを理解しておく必要があります。（群論などは絶対必須というわけではありませんが、知らないと後々苦労することがあります。）

バンド計算と言うからには計算機や計算機言語に対しても深い知識が必要です。普通、第一原理電子状態計算のプログラムのソースコードは、全部で１万行（コメント等を含めて）以上にもなります。人によってこれが大きいか、小さいかは異なると思いますが、普通の計算物理未経験の研究者にとってみればこれは非常に巨大な部類に入ります。できれば千行以上のソースコード（用途は問わない）を作成、運用、保守の経験があるのが望ましいと思います。

使用言語としては、大抵のバンド計算プログラムはＦＯＲＴＲＡＮ（それも７７）で書いてあります。古いものではＦＯＲＴＲＡＮ６６以前の形式で書かれているものもあります。最近、ＦＯＲＴＲＡＮ９０も大部普及してきました。もし、全くゼロの状態から作り始めるのなら、ＦＯＲＴＲＡＮ９０の仕様でプログラムをコードした方が良いかもしれません。（少なくとも、７７仕様でもすぐにかつ楽にＦＯＲＴＲＡＮ９０に移行できるようにしておくと良いと思います。）

最近、ＣやＣ＋＋でバンド計算プログラムを書いたり、ＦＯＲＴＲＡＮから移植する例もあります。私個人としてはＣ、Ｃ＋＋が必ずしも向いているとは思いませんが、どのような言語で書くかは当事者の自由です。（ＪＡＶＡで書こうとする者も出てくるかもしれません。）

バンド計算プログラム作成前に注意しておくべきこと

プログラム作成一般に言えることと思いますが、以下のようなものが挙げられます。

プログラムに汎用性があること

可読性が良いこと

いろいろな機種、ＯＳで動くこと（を想定して作る→汎用性）

計算対象マシンがベクトル計算機ならベクトル化について留意する

著作権等の留意（特に複数のユーザーで作る時や、公開版として作る場合）

改造、拡張、他機種への移植が容易なこと（汎用性、可読性）

以上を実現するために、プログラムコーディング前に十分計画を練る

最初の段階で、コーディングの方針と作業計画を十分練る必要があります。まず全て全くゼロから作るのか、誰かのプログラムを拡張させるのかで、大分方針が異なってきます。特に元になるプログラムがある場合、その元プログラムの汎用性や可読性が、新たに作っていくプログラムの運命を決定付けてしまいます。

また元プログラムが存在する場合、その設計方針（使用言語、稼働対象マシン、ＯＳ、プログラム仕様等）に従わなければなりません。もし従わないと、結局ゼロから始めた方がマシになります。

言語をＦＯＲＴＲＡＮ７７にするか９０にするかは、大変重要な問題です。ＨＰＦ（High performance FORTRAN）という選択もありますが、これはＦＯＲＴＲＡＮ９０程（日本では９０でさえまだ普及しているとは言えない）には普及していません。

ＦＯＲＴＲＡＮ９０は（１）動的なメモリの使用が可能、（２）モジュール文の存在（７７でのコモン文はいらなくなる）、（３）構造体の存在、（４）強力な行列演算機能、（５）行列表現の簡略化、（６）１行に複数の命令が可能、（６）再帰呼び出し可能、（７） SELECT---CASE文、DO---WHILE文、DO---END DO文等新しい命令の存在、（８）ポインターの概念の存在などなどが挙げられます。

特に、（１）の機能は大変魅力あるもので、これは実行中に必要な配列サイズを、事前に考慮しなくても計算を遂行できるようにすることが可能で（機種とＯＳ、計算機運用環境に依存し、バッチで処理して複数のユーザーが使用するシステムでは動的メモリ割付けは実際上不可能）、非常に柔軟な計算が可能になります。例えば、系のサイズの変化や、計算条件の変更で、いちいち配列のサイズ変更による再コンパイルを繰り返す必要がなくなります。

その他に、構造体の概念や行列演算部分の簡略表示などは、うまく利用すれば、プログラムを非常に使い易いものにし、汎用性や可読性をも向上させることができます。一方、バンド計算は基本的に、たして、ひいて、かけて、わって、だけで計算できてしまうもので、このようなものにポインターのような概念が本当に（バンド計算エンドユーザーにとって）必要かどうか甚だ疑問ではあります。（少なくとも必須のものではないし、これまで、ポインターの概念がなくてもバンド計算プログラムはきちんと動いてきました。（当然、非常に計算に達者な者はこの限りではありません。）

もし既存のプログラム（ＦＯＲＴＲＡＮ７７仕様）を改良、拡張するなら、無理をしないで、ＦＯＲＴＲＡＮ７７でプログラムをコーディングしていったらいいと思います。ただ近い将来ＦＯＲＴＲＡＮ９０に移行することを考えてコーディングしておいた方が良いと思います。（なるべくコモン文は使用しない、ＤＯ文はＥＮＤＤＯで終るようにする、ちゃんと汎用性、可読性を考慮してプログラムを書けば、９０への移行はそう困難なことではないと思います。）

既存のコードを流用する場合、元のプログラムが非常に、汎用性や可読性に欠ける場合はどうすればよいでしょうか？
結論としてはどうしようもないと思います。自分で、修正するか、汎用性や可読性のなさに目をつむってそのまま拡張、改良していくしかないと思います。なるべく元になるコードも、汎用性、可読性の良いプログラムを使用することを勧めます。（ただ最初に見ただけでは、なかなかそれがわからないとは思います。）

実際のコーディング

では次に、実際のコーディングについて説明を行ないたいと思います。

自分の所有するプログラムは第一原理分子動力学手法を使って、系の安定構造とその電子状態を計算するもので、これは自分が東大物性研の寺倉研にいた時代に作成を開始しました。（寺倉先生は現在工業技術院産業技術融合領域研究所に所属しています。）

これには元になるプログラムがあり、これは当時寺倉研の助手であった石田さん（石田先生は現在、日大文理に所属しています。）が作ったものでした。この段階で、ＢＨＳの論文（G. B. Bachlet, D. R. Hamann and M .Schluter, Phys. Rev. B26, 4199(1982)）の表のパラメーターをもとに平面波基底を用いたバンド計算が可能でした。これに力の計算部分を加え、Kleinman-Bylanderの分離形（ L. Kleinman and D. M. Bylander, Phys. Rev. Lett., 48, 1425(1982)）を導入し、１９９０年春にはカー・パリネロ計算（第一原理分子動力学計算）が可能となりました。

この１９９０年春の段階では、プログラムとしては、最急降下法による対角化を用いない電子状態の計算と、ユニットセル内の力の計算ができ、それをもとにした初歩的な分子動力学計算ができるだけでした。

ここで、現在のプログラムの先頭部分を示したいと思います。
これはプログラムの最初の注釈行です。御覧のように１９８９年５月４日から本格的なカー・パリネロプログラムの作成を開始し、いろいろな改良（部分内殻補正や対称性の部分は現在産業技術融合領域研究所の森川先生が作成した。）、拡張を行ない現在に至っています。現段階で、このプログラムは以下の計算が可能です。

ＢＣＣ、ＦＣＣ、ＨＣＰ、ＳＣ、正方晶、斜方晶等の結晶対称性を考慮した計算が可能

ｓ、ｐ、ｄ、（ｆ）を非局所ポテンシャルとした計算が可能

磁性を考慮した計算が可能（ここで示しているプログラムのバージョンでは不可）

ユニットセル内の原子に働く力以外に、系そのものに懸かるストレス（圧力）の計算が可能

完全ではないが、上記ストレスを使用して、定圧条件での構造最適化が可能

部分内殻補正を考慮した計算が可能（力、ストレスに対する補正も考慮）

富士通のベクトル並列計算機ＶＰＰシリーズでの並列計算が可能

直接プログラム内容とは関わりないが、ＮＣＰＳ９５を使って、非常に広範な範囲の電子状態計算が可能

逆に、現在できない機能として以下ののものが挙げられます。

ＬＤＡを越えた計算はできない。（ＧＧＡ、ＳＩＣなどはできない。
但し、現在ＳＩＣ擬ポテンシャルを検討中、３／１８、１９９７）

スピン軌道相互作用は考慮されていない。

電子状態計算部分は基本的に絶対零度の計算（有限温度の計算は不可）

非常に複雑な対称性を持つ系に対しては、対称性なしで計算するしかない。

ＶＰＰ以外での並列計算（ＰＶＭ、ＭＰＩ等）には対応していない。

プログラム全体の構成

では次に、プログラム全体の構成について話したいと思います。
まずプログラムのサブルーチン構造ツリーを示します。

次に、プログラムのおおまかな流れ（ループ）を示します。
LOOP0は、最も外側のループで、一定圧力下での構造最適化を行なう場合のものです。この場合、ユニットセルを規定する格子パラメーターそのものが最適化されるので、ループの範囲はこのようになります。

LOOP1は、通常のカー・パリネロ法による、系の内部構造とその電子状態を最適化するループです。
LOOP2は、電子状態部分のみ計算する場合のループで、サブルーチンＭＤは計算されません。

更に、各サブルーチンのおおまかな説明を示します。

次に、このプログラムで使用しているインクルードファイルと入力データを示します。
インクルードファイル”PACVPP”
インクルードファイル”REINIT”
入力データ”I.DAT”
入力データ”CONTL.DAT”

他にも、擬ポテンシャルデータ入力ファイルがありますが、これは擬ポテンシャルの説明の記事で言及されているので略します。

インクルードファイル”PACVPP”は、主に変数（含む配列変数）の定義された部分です。これをプログラム本体（revpe_d.f）中のinclude文で読み込みます（REINITも同じ）。PACVPPではパラメーターの設定、主にコモン文による配列変数の定義が行なわれています。最後の方にある、注釈行!XOCLで始まる部分は富士通のベクトル並列スーパーコンピューターＶＰＰ５００用の並列化制御のための指示命令行です。

これを見てわかるように、コモン（大域定義）文によって定義された変数や配列変数が多数あります。前述のように、可読性や汎用性、移植性を考えると、このような変数定義の仕方はお勧めすることができません。あまり良くない例と思って見てやってください。

インクルードファイル”REINIT”は、ファイル入出力のためのOPEN文に始まり、計算制御用入力データCONTL.DAT入力部分（READ文使用）とからなります。

入力データ"I.DAT"は、この場合ヨウ素（Ｉ）の入力データを意味しています。データの内容は主にユニットセルの格子パラメーターと、内部原子座標に関するものです。但し、一部計算制御用のパラメーターが混じっています。また最後の行に局所擬ポテンシャル計算で使われるＶcore用のパラメーターがあります。

入力データ”CONTL.DAT”は、計算制御用のパラメーターが入っています。外部から懸かる圧力、結晶構造の種類、ＰＣＣを考慮するかしないか、イタレーションの数などが主なパラメーターです。

ここでの重要な問題は、二つの入力用データが用途として完全に独立していないことです。座標定義用のI.DATに制御用のパラメーターが混じっています。これはデーター内容の理解や、プログラム運用を甚だ面倒なものにします。

新たに、バンド計算プログラムを作ろうという方は、入力データの設計には十分な配慮（データの割り振りと定義をきちんとする）が必要と思います。行き当たりばったり的な入力データの構築をすると、後で後悔することになります。

プログラムのメイン部分

では次に、プログラムのメイン部分を示したいと思います。

メインルーチンだけでも４００行もあります。ここでは第一原理分子動力学計算の制御と、サブルーチンの呼び出しが行なわれているだけです。
このルーチンの問題点は、第一原理分子動力学計算の制御にあります。この計算は３つに場合分けされています。

（１）第一原理電子状態計算のみ
（２）第一原理電子状態計算＋ユニットセル内原子の構造最適化
（３）（２）＋ストレスを使った系のユニットセルの最適化（定圧問題用）

最初は（１）と（２）の場合のみの計算に対応していました。この段階でも、（１）のみの計算と（２）の計算を、最初の設定条件（CONTL.DATで設定）によって計算を切替えるようになっていて、メインルーチンのプログラムの流れを大変分り難いものにしていました。その後、ストレスの計算を導入し、系のユニットセルそのものを最適化できるようにするため、（１）、（２）、（３）の各計算を切替えて計算できるようにプログラムを書き替えたのですが、（３）の付加は、ほとんど場当たり的に、無理に付け足したものとなり、メインルーチンの構造は、複雑怪奇なものになってしまいました。

プログラム作成においてメインルーチンは、そのプログラムの顔とも言える部分であり、ここが分りやすくすっきりしているか、そうでないかで、プログラムの見ためや、可読性を大きく左右します。教訓としては、先の入力データ等のところでも述べたように、事前の設計や計算方針の決定は十分に時間をかけて、練りに練った上で、コーディングを開始すべきだと思います。

尚、メイン文がどのようなサブルーチンを呼び出しているかは、前述のサブルーチンツリーを参照して下さい。第一段目にコール（呼び出し）されているものがメインルーチンで呼び出されているサブルーチンです。

では具体的に良いメインルーチンの書き方を考えて見ましょう。まず、

（１）どのようなバンド計算をするか、方針をきちんと決める。
（２）どのようにコードするかじっくり考える。
（３）計画、方針に則って、メインルーチンを書く。
（４）メインルーチンでは、具体的な計算を行なわない。
（５）基本的には、サブルーチンコールだけにする。
（６）呼び出すサブルーチン毎に何をするサブルーチンか注釈をつける。
（７）各サブルーチン毎に、時間を計れるようにしておくと、実行時での各サブルーチンの処要時間を割出すことができる。
（８）サブルーチン名も、なるべく誰が（少なくとも他のバンド屋が）みてもわかるようなものにしておく。（例、DIAGON、EIGEN、FERMI、VELRET、SD等）
（９）計算が電子状態だけか、構造最適化も行なうかなどの場合分けは、すっきり、分り易くする。（なるべく最初の段階で、方針を決めたら、後で安易に変更しない。）

などが考えられます。

先ほど言及しましたが、このバンド計算プログラムでは、２つのインクルードファイルを読み込んでいます。[PACVPP][REINIT]

PACVPPは主に、変数（配列変数）等の定義、REINITはオープン文によるファイルのＩ／Ｏ定義を行なっています。PACVPPは、メインルーチンを含めのほとんど全てのサブルーチンでも取り込む（インクルードする）ようになっています。REINITはメインルーチンで１回だけ取り込みます。

OPEN文に関しては、I／O用のファイル名を変更する毎に、プログラムをコンパイルし直さなければならないので、少々不便です。ＯＳやベンダーに依存しますが、ＯＳ上の環境設定段階で、ファイル入出力を制御することも可能です。例えば、昔のメインフレーム大型コンピューターでは、ＪＣＬを使ってOPEN文を使わないで、ファイル入出力が制御でき、この場合、プログラムの再コンパイルをする必要はありません。UNIX環境でも、これに似たことをすることが可能です（確か。setenvを使う）が、ここではこれ以上の説明はしません。各自で、自分の環境に対応した方法を探してください。

メインルーチンの内容（各論）

ここで、メインルーチンにおける、計算の大まかな流れを説明したいと思います。
（１）まず最初に行なうのが初期設定部分で、計算するべき系の格子パラメーターや、擬ポテンシャルデータ、計算制御パラメーターの入力。これらの入力データをもとにした、逆格子空間座標の計算。必要（計算制御パラメータで指示してあれば）ならば、空間群による対称性の計算。初期電荷密度の計算（ガウス関数を用いた初期電荷密度を使用）。エバルト和の計算（いわゆるマーデルングエネルギーの計算）。フォームファクターの計算。などが挙げられます。

（２）次に、読み込んだ擬ポテンシャル（局所、非局所）のフーリエ変換を行ないます。以下に、局所ポテンシャルの座標変換の式を示します。

V_loc(G)=int{j_0(G*r)*V_loc(r)*exp(-i*G*r)}/Ω

ここで、Ωは体積、V_loc(G)は逆格子空間での局所擬ポテンシャル、 V_loc(r)は実空間（動径方向）の局所擬ポテンシャル、j_0(r)は０次球ベッセル関数、intは積分の意、G、rはそれぞれ逆格子空間座標、実空間（動径）座標です。

非局所擬ポテンシャルは、局所の場合より複雑になります。（従って、表式による説明は行ないません。）特に、このプログラムではＫＢ近似による分離形を使っているので、計算手順がいろいろと複雑になっています。これは後のプログラムの説明を見ていくと次第に明らかになっていくものと思います。

さて、久しぶりに続きを書いていきたいと思います。（１１月７日、１９９６年）

（３）次に、交換相関エネルギーまたは交換相関ポテンシャルの計算を行なっています。どちらの計算をするかは、状況によって異なりますが、メインルーチン上では、交換相関ポテンシャルの計算を行なっています。（どちらを計算するかわ、サブルーチンの引数によって選択できるようになっています。）
これは後で、電子状態計算のための行列要素の計算を行なう時に使います。

（４）一回目のイタレーションでは、行列の対角化を使って電子状態を解きます。これは初期固有エネルギーと初期固有ベクトル（つまり波動関数）を用意するためです。これら固有エネルギーと固有ベクトルを使って、次のイタレーションからカー・パリネロ的な逐次的な実行により波動関数を求めていきます。

対角化の際、解くべき行列の要素数は、実際の平面波の数より少なめにとります。本計算では普通、ＧＭＡＸ／２（つまり通常の平面波の１／８の空間）に対応する行列の対角化を行なっています。これだけでも、初期の固有ベクトル（波動関数）としては十分です。非常に大きな系を扱う場合は、更に、初期の行列要素数を少なくしますが、その分、電子状態計算の収束が遅くなります。（当然、極端に少な過ぎると、計算が収束せずに、破綻する可能性があります。）

２回目以降のイタレーションでは、カー・パリネロ的な手法（実際は時間の１次微分方程式として扱う、最急降下法〔ＳＤ法〕）を用いて、電子状態を解いていきます（固有値、固有ベクトルを求める）。解き方の内容の詳細は、私の書いたＤ論（東京大学、１９９０年度）に書いてあります。現在も方法論的にはほとんど進歩していません。カー・パリネロ法に関しては、カー・パリネロの最初の論文（R. Car and M. Parrinello, Phys. Rev. Lett., 55, 2741(1985)）をはじめ、多数の論文や詳しい解説（日本語によるものもある）が出ていて、探せば容易に入手することができます。

メインルーチンでは、一回目と二回目以降で、対角化とＳＤ法による場合分けを行なっています。以下参照。


      IF (ITER.EQ.1) THEN                                               
          CALL TIME(ITIME1,'      ',1)                                  
          CALL KBMAT 
          CALL DIAGON(IREC8)                                            
          CALL KBINT 
          CALL TIME(ITIME1,'KBMSD ',2)                                  
        ELSE                                                            
          CALL TIME(ITIME1,'      ',1)                                  
          CALL MSD(IREC8,IREC9,IMD,ICAR,ISTR) 
          CALL TIME(ITIME1,'MSD   ',2)                                  
      END IF

これは、比較的単純なIF文による場合分けです。メインルーチンにはもっと複雑怪奇な場合分けをしている部分もあります。当然、条件判定を変更して、最初から最後まで、行列の対角化のみで、電子状態の計算を遂行していくことも可能です。但し、最初から（対角化による初期固有ベクトルなしに）ＳＤ法で始められるようにはなっていません。

バンド計算グループによっては、初期固有ベクトルを乱数を使って、求めておくという手法を用いているところもあります。この方法では対角化という比較的面倒な手法を用いないので、労力や計算時間が稼げますが、その分、電子状態の収束のスピードが遅くなります。（曲がりなりにも、対角化で求めた固有ベクトル（波動関数）の方が真の解により近いため。）

（５）次に、電子状態（固有値、固有ベクトル）からフェルミエネルギーを求めます。フェルミエネルギーは、実際は固有値（固有エネルギー）つまり、バンドの底から、扱っているユニットセル内の総電子数分を数えていき、１つのバンドに２個（up-spinとdown-spinで２個分）電子が詰まるとして、その詰まり切ったところ（位置）を、フェルミエネルギーとしています。
詰めていく場合、各バンドとｋ点について足していきます。
Σ_i,k[ε(i,k)]=N

ε(E,k)が固有値で、これをi（バンド数）とk（ｋ点）で足していき、足した数Nが、ユニットセル内の総電子数と一致した場合のε(i,k)の値がフェルミエネルギーとなります。但し、バンドとｋ点はとびとびの値で、固有値自身も連続ではないので、金属的なバンドで、フェルミエネルギーを決める場合は、なんらかの内挿を行なうか（例えばテトラヘドロン法〔四面体法ともいう〕など）、バンドに適当な幅を与える（大抵は、フェルミ分布関数か、ガウス関数によってバンドにエネルギー方向に関しての広がり〔ブロードニング〕を与える）などの方法をとります。これは一種の温度の概念を導入することでありますが、本当の意味で密度汎関数計算で温度の計算が可能になったわけではありません。我々の計算でも後者の。バンドに広がりを与える方法を使いますが、方法はもっと簡単で、放物線関数を使った非常に簡単な関数を使って、バンドの広がりに関しての積分を解析的に行なってしまいます。（詳しくは私のＤ論参照）

（６）フェルミエネルギーが求まると、これによって力や電荷密度、全エネルギー、ストレスなどの計算が可能となります。これは、これらの物理量がバンドとｋ点に関しての足し上げによって求まり、その足し上げがフェルミエネルギーのところまでだからです。つまり（５）でフェルミエネルギーが求まっていなければ、力、電荷密度、全エネルギー、ストレス等は求めることができません。フェルミエネルギーを求めるサブルーチン、FERMIを実行した後は、まず力の計算が行なわれています。但し、第一原理分子動力学計算では、最初電子状態のみを断熱ポテンシャル面に落すため、ユニットセル内の原子に働く力の計算を行なう必要がなく、電子状態の計算のみが必要となります。また、分子動力学計算を行なっている時も、毎回（毎イタレーション）力の計算を行なわず、数イタレーションに一回力の計算を行ない、原子をその時だけ動かすとした方が、計算全体としての収束が速い場合があります（と言うより、ほとんどの場合で収束は速くなります）。このプログラムでは、力の計算を行なうサブルーチンで、非局所擬ポテンシャル用の行列要素の計算で必要な部分の計算も一緒にやらせており、力の計算を必要としない場合は、この行列要素のための計算のみを行なうサブルーチンが用意されています。そして必要に応じて切替えてサブルーチンが呼び出されるようになっています。以下に、そのプログラム部分を示します。


      IF (MOD(ITER,IMDI).EQ.IRES) THEN
        IF (ITER.EQ.1 .OR. ITER.GT.IMD) THEN                            
          WRITE (6,*) 'CALL FORCE'
          CALL FORCE(IREC8) 
        ELSE
          CALL FORZFB 
        END IF
      ELSE
        WRITE (6,*) 'CALL FORZFB'
        CALL FORZFB 
      END IF

先の、DIAGON、MSD と比べて、場合分けがかなり複雑になっています。まず最初のＩＦ文では、何回かに１回の割で、力の計算を行なうように設定しています。内側にあるＩＦ文では、外側のＩＦ文の条件を満たし、かつ第一回目のイタレーションか、またはＩＭＤ回より多くイタレーションが進んでいる場合に力の計算が行なわれます。そうでない時は力を除いた部分のみの計算（CALL FORZFB）が行なわれます。制御構造が大分複雑になり、プログラム作者でさえ、ちょっと考えないとわからないくらいになっています。

（７）次に、電荷密度の計算を行なっています。これはサブルーチンＭＳＤで計算される波動関数（固有ベクトル）を使って求められます。この時、波動関数同士のかけ算が、畳み込みの形になっているので、ＦＦＴ（高速フーリエ変換）を使って、計算量をＮの自乗のオーダーをＮｌｏｇＮのオーダーまで少なくすることができます。

（８）その後、必要ならばストレスの計算を行ないます。一応、一定圧力の条件でのユニットセルの形そのものの構造最適化も可能です。サブルーチンを呼び出している部分を以下に示します。


          IF (ISTR.EQ.1) THEN
              CALL TIME(ITIME1,'      ',1)                              
              CALL STRESS(IPCC,SCHGPC
     &       ,ETOT1,TOTCH,EPC,PCM,KOPR,KBZTYP) 
              CALL TIME(ITIME1,'STRESS',2)                              
          END IF

ここで、制御用の変数ＩＳＴＲがクセ者で、一応ＩＳＴＲの値が１の場合に、ストレスの計算を行なうようにさせていますが、この場合、ユニットセルの形が変わることを想定したプログラム制御を行なう必要があり、そのための制御をするための部分を、後から（場当たり的に）付け足したため、メインルーチンの構造が複雑怪奇なものになってしまいました。

（９）次に、非局所部分以外の力の寄与の部分の計算と、場合によっては対称性を考慮するための計算を行ない全ての力を足し上げます。その後、この力を使って分子動力学計算（ユニットセル内の原子の最適化、及び、必要ならストレスを使って、ユニットセル自身の構造最適化）を行ないます。分子動力学そのものは、単純なベレのアルゴリズムを使った方法を用いています。
第一原理分子動力学計算では、この分子動力学計算部分は計算全体の１％にも満たないので、あまりプログラムに関してテクニカルな方法は用いていません。

（１０）そして、最終的に系の全エネルギーを求めます。ここでは交換相関エネルギーが必要なので、全エネルギー計算のためのサブルーチンで、交換相関ポテンシャル、エネルギーを計算するためのサブルーチンを呼び出しています。

（１１）最後に、全エネルギーに関しての収束判定と、収束していない場合、次のイタレーションに進むために新しい（更新された）電荷密度と、古い電荷密度の混合を行ないます。本当は、ＭＳＤでの波動関数の更新処理の段階で、電荷密度の混合と類似（厳密には等価ではない）な処理を施しているのですが、実際の計算では、ここ（収束判定のところ）でも新しい電荷と古い電荷の混合を行なっています。混合の仕方は単純混合です。この方法を用いた方が、時間の刻み幅を小さくするより、経験的に収束がずっと速いことがわかっています。
（注意）”電荷密度の混合と類似”：これは波動関数の混合操作を意味しません。また、電荷密度の混合は各イタレーション毎の最後（収束判定）のところで一度だけ行ないます。

単純な混合より、ブロイデンの方法など、より収束を速める混合手法がありますが、現在このプログラムでは使用していないので、これ以上の説明は行ないません。

１２月１６日（１９９６年頃？）、久々に話しを進めて行きたいと思います。一応、以上の説明でメインプログラム部分の説明はおしまいです。これでプログラム全体の概要がある程度説明できたと思います。次の[各論編]へ進む。

（寄り道）
不肖この私目も物性研や原研のＶＰＰを利用して、プログラムの並列化、及びその実行を行なっています（既に、過去の話：１１／２６、２００２）。この並列化作業の内容についての話は、[並列化（ＶＰＰ）]を参照。

[先頭][総目次 ][最初に戻る][各論編][デバッグ編 ]