(計算プログラム、講議ノートは 英文資料 へ移動しました。)
★ プログラムの達人については、こちらを参照。
★ ちょっとした報告書。 ★ 神戸大学公開講座の予稿。
★ 京都大学
で公演しました。(1998.3.17)
数値計算に「噂」はつきものなのですが、果たして真実やいかに?
Q1. システムサイズ N や、採択状態数 m を増やすと、計算が不安定
になってしまって、信頼性のある結果が得られない?
--- これはウソ。但し、プログラムを作るにあたって注意を払う
べきポイントが幾つかあって、
(a) 全系の波動関数を求める Lanczos 対角化の数値精度は充分か?
また、得られた波動関数を逆反復法で更に改善してあるか?
(手抜き計算では、結果も押して知るべしです。)
(b) 状態数の制限をする時に、「ちょうど m 個で切る」のはダメ。
系が SU(2) 対称性や Z2 対称性を持っていたら、密度行列の
固有値は、その対称性を反映して縮退するので、縮退している
ものは「全部採るか、全部捨てるか二つに一つ」である。これ
は White が原論文で注意してあるポイント。従って、実際に
採択する状態数は「最大値が m の、ある整数」になる。この
注意を怠ると、系のサイズが大きくなった時に数値的不安定性
を招いて、計算が停止する。
(c) バグはないか? --- なにをバカな、と思われるかもしれないが、
意外と良くみかけるケース。例えば、「m 個状態を取る所を、
間違えて m - 1 個にしていた」というケースでは、m が充分
大きければ一見マトモな数値データを得るが、やがて破綻する。
この手の「密度行列の固有値が小さな状態に関するバグ」は、
発見が難しい --- いや、計算の途中で (b) で述べた固有値の
縮重が解けて、やがて計算が破綻するので、慣れて来るとバグ
を見つける為の重要な手がかりになる。
(e) 結構良く見かけるバグが Fortran のメモリー領域違反。m を
大きくすると、色々な変数の占有領域が増えて行くので、バグっ
て変数領域が重なっていたら、一巻の終わり。ここで計算が
破綻してくれたら、まだ発見が容易なのだけど、「何か結果
がおかしい? 」という症状のみの場合もあるので、最終的に
は物理的におかしくないかどうかのチェックは必要不可欠。
なお、間違ったデータが出版されてしまったら、もう手遅れ
なので、東に向かって手を合わせて「投稿料はナンマンダ?」
(f) そもそも、なぜ密度行列の固有値が縮重するのだろうか? それ
は、計算対象の物理系が何らかの対称性を持っているからであ
る。だから、最初から対称性を考慮して、無駄な自由度を落と
してから計算を始めるのが一番良い。
Q2. 励起状態の計算精度が悪い?
これは半分真実。励起エネルギーを求めるには 2 つの方法があって
それぞれ異なる結果を与えます。
(a) 基底状態をターゲットした計算の最中に、Lanczos 対角化の第
二固有値を励起状態のエネルギーとして採択することが可能。
但し、こうして得られた励起状態の精度は悪くて、基底状態の
エネルギーの半分くらいの精度しか出ない。その理由は、ター
ゲットしていない状態のエネルギーを無理に求めている為であ
る。実は、もう一つ理由があって、こうして求めた励起状態は、
「左右のブロックの継ぎ目あたりに弱く局在している」状態で、
真の励起状態よりも高いエネルギーを持っている。
(この方法は、どちらかというと Infinite Algorithm 向け。)
(b) 量子数 (又はパリティー) を色々と変化させて、各量子数で指定
される部分空間の基底状態のエネルギーを DMRG で求める。そ
うして、最後にエネルギーの量子数差分を取って、励起エネル
ギーを評価する。この方法は、Finite Algorithm のみで可能
なので、最後にスケーリング解析をしてサイズ無限大の極限を
取らなければならない。
(a) の方法が簡単なので、まずこちらで解析してみて、より正確な
結果を得たい場合に (b) の方法を利用して励起エネルギーの計算
精度を上げるというのが良い方策だと思います。もちろん、最初
から (b) で攻めても結構です。
(b) の方が精度が良い理由ですが、例えばスピン系のスピン励起の
計算をしてみればわかるのですが、(b) では全 Sz を 1 だけ変化さ
せて、しかもその Sz 編極が系の全体に分布しているのに対して、
(a) の場合では「左右の系の切り口のそば」で擬似的に triplet を
作る励起を拾っているにすぎないというのが、主な要因です。これ
については、 Ostlund らによる離散化ラプラス変換の議論を参照
するのが良いかもしれません。
気をつけるべき点は、いずれの方法を用いるにせよ m を変化させ
てみて、Gap の m 依存性を把握しておくことです。単一の m の
みによる計算結果のみでは、客観的にみて充分に信用できるとは言
えません。
Q3. 相関関数の計算精度が低い?
まず、相関関数の計算法について。これも、2 種類あります。
(a) 弱い外場を「あるサイト」にかけて、その周囲のサイトの磁化
や粒子密度の変化を測定する。系の両端付近の密度変化を測定
しても良い。(静的) 相関関数と、(静的) 局所外場応答の関係は
既知である --- フリーデル振動を参照 --- ので、結果的に相
関長や振動周期を求める事が可能。 (この方法は Finite Algorithm
のみで可能。)
(b) 波動関数は「繰り込み群変換を表わす直交行列の積」で書かれ
る事を応用して、計算の結果得られた行列積型の波動関数から
相関関数を計算する。 (これは、Finite / Infinite の両方で
可能。)
(a) の方が計算が簡単で、精度も高いので、こちらをまず推奨します。
但し、Off-Diagonal な相関などは、どうしても (b) でないと測定で
きない場合がありますので、そういう場合は (b) を用いて下さい。
相関関数も、m に依存しますので、m 外挿は忘れずに行って下さい。
Q. 算出したサイト当たりのエネルギーが、低すぎる?
DMRG は変分法なので、N サイトの系の全エネルギーを計算すると、
それは真の基底エネルギーよりホンの少し高い値になっています。
但し、各サイト間のローカルなエネルギーは、「場所によって真の
ものより高かったり低かったり」します。変分原理は、あくまでも
系全体のエネルギーについてのみ成立しているのです。B・・B と
いうブロック構成を考えた場合には、真ん中の "・・ペア" 間のエ
ネルギーは、通常の場合真のエネルギーより低く出ます。
Q. ランチョス方を何百回も繰り返し使うので、計算時間が長い。(?)
DMRG を加速するテクニックが色々とありますので、この問題は
解決可能です。順番に並べると、
(a) 2-Path Method: ランチョス法を加速する方法として、Lanczos
Vector を作って逆反復法に渡す "2 Path Method" という方法が
知られています。逆反復法をイキナリ任意の初期ベクトルで始める
よりも、基底ベクトルを求める時間が短縮されます。
(b) 原点移動: 逆反復法では、CG 加速ルーチンから渡される基底ベク
トルの候補が、最適な変分ベクトルになっているかどうかチェッ
クします。その際に、変分エネルギーがランチョス法により求め
られた基底エネルギーの近似値を下回れば、逆反復法で使用する
基底エネルギーを「より低い方」に更新します。こうする事によっ
て、逆反復法の収束を改善することができます。
(c) 打ち切りパラメターの調整: ランチョス、逆反復、CG 法はそれぞれ
反復法で、反復打ち切りの内部パラメターを持ちます。これをチュー
ニングすることによって、計算精度を落とすことなくプログラムの
実行速度を数倍することが可能です。
(d) 波動関数の再生: White は 1995 年暮れに Finite DMRG Method
に際して、初期ランチョスベクトルとして「ほぼ基底ベクトルに
等しいものを合成する」方法を発表しました。この方法を用いると、
ランチョス法の反復回数が 10 ランチョス・ステップ以下になり、
実行速度が 10 倍以上高速化されることがわかっています。
(e) 必要なデータだけ計算する: 良く、相図を書こうとして、相境界辺
りの励起エネルギーをバカ丁寧に計算する方がいらっしゃいますが、
これは不毛な努力です。Gap が閉じる相境界付近は、系の揺らぎが
大きく、大規模な数値計算をもってしても、有限サイズ効果その他
が見えてしまいます。まず相境界から外れた所でデータを集めて、
そこから外挿可能かどうかを観察してみて、外挿不可能なら始めて
「君子危うきに近寄る」べきです。勿論、データの信頼性を維持す
る為に、何箇所かは相境界まで丁寧にデータを取らないといけない
のですが、「点 10 個でグラフが書ける所を 100 点計算する」の
は計算費の消化しか頭に無い計算バカのやる事です。
(f) 対称性を使う: 対角化一般に言えることなのですが、系に対称性が
あれば、予めハミルトニアンを対称性で「割って」おいて、ヒル
ベルト空間の次元を下げておくのは当り前です。こうする事によっ
て、計算速度が高速になるとともに、対称性の指標 (量子数) が違
う空間で、それぞれ基底状態を求める事が可能になります。
まあ、以上のような工夫をすることによって、計算を 1000 倍程度加速
することが可能になります。 1 分で済む計算に半日を費やす人を時々み
かけますので、御注意あれ。
Q. Ferromagnetic Bond を含むスピン系に DMRG を使うと計算
が破綻する。(?)
これは、半分ほんとう。Ferromagnetic Bond を含む系では、DMRG
に現われる左右の Block のスピン量子数 (の平均値) が大きくなり、
SU(2) 対称性により S = J の空間は 2J+1 回繰り返しヒルベルト空間
に現われることが理由で、Block を表現する基底の数が増えてしまう。
....直ちにわかる様に、Ferromagnetic Bond が何サイトも続けて出て
来ると、計算が破綻してしまう。
原因は、SU(2) 対称性を用いて Hamiltonian の次元を落としておか
なかったから。IRF 基底またはプラケット基底と呼ばれる SU(2) 対称
性を最初から取り入れた基底を用いて、Hamiltonian の表現行列を
書き表せば、上記の困難は克服される。 (これについての説明は長く
なりますので、e-mail で個別にお問い合わせ下さい。)
Q. Fermi系の「追い越し符号」は、どう取り扱うの?
長距離ホッピングが存在する場合は、演算子に対して Fermion sign
を考慮しなければなりません。これは面倒ではなくて、符号を考える
必要があるのは「新たにサイトを付加して、ブロックを拡大する時」
だけですから、その際に元々の Operator の行列表現に対して加える
サイトの電子数分だけマイナス符号をつければ良いだけです。
Q. 有限温度の方法とは、どんな方法?