LaTeXで図が真っ白になる問題について(texlive/2013)(Windows)
論文を書いていて、TeXworksでpLaTex(ptex2pdf)でpdfを作製。
どれどれと思って見たら、epsの図のところだけポッカリ真っ白。
LaTeXを書いていてeps画像がPDFファイルに出力されない問題を解決 - 透明の練習 第2刷
パソコン雑記: TeXでEPSファイルがPDF出力できない現象について(2)
を参考にして見たが、結局上手く行かない。
というか、そもそもdviファイルが出て来てなかった。
「pLaTexはダイレクトにtex -> pdfに行くから、dviは出ない」という先入観が良くなかった。
あと、コンソールを常時表示させていなかったのも良くなかった。
コンソールを見ると、Warningとあり、どうやらGhostscriptを「gs」のコマンドで呼んでいてコケていた。
ん?と思って、自分のC:/gs/gs9.20/binを見ると、gswin64.exeがいらっしゃった。
つまり、gs.exeをLaTexは探していたが、あるのはgswin64.exeだからシカトされていたわけです。
gswin64.exeをコピーしてgs.exeを用意したら、ちゃんと図が出るようになりました。
logファイルを見ればいいやと思ってたら甘かった。。。コンソール大事ですね。
VASPでEMAXを増やして計算したらcore dumpする件について(Ubuntu12.0.4)
VASPのコンパイルの問題ではなく、Linux側でデフォルトの設定だとstackの上限が8MB程度と小さいことが原因(あくまで自分のケースでは)。オーバーフローとか別のエラーを吐いて欲しかった。。。
ulimit -a
でいろんな上限が見れる。
一時的には
ulimit -s unlimited
でstackの上限を外せるが、ターミナルを閉じたりするとデフォルトに戻る。めんどい。
/etc/security/limits.conf
を弄ってstackの上限を外す。自分は適当に、
* soft stack unlimited
* hard stack unlimited
と付け加えた。*は全てのユーザーで、soft/hardは一般ユーザー/rootの上限を設定することを表す。ちゃんと理解していないので、ユーザーの指定とsoft/hardの役割が被っている気がしてならないが、VASPが動けばどうでも良い。
変更を反映させるにはrebootする必要がある。うちのワークステーションはSDDなのであまり気にならないが、立ち上がりが遅いと面倒な作業ではある。
変更されたかを見るには
ulimit -a
をターミナルに打ち込む。ulimitをunlimitだと思ってターミナルが反応せず一瞬頭を抱えたのは内緒。
VASPの問題なのかCPU側の問題なのかを判別するのは非常に難しいと痛感したトピックスではあった。
=== 追記 ===
stackの上限が問題になることももちろんあったが、追記として「Intel MKLとopenmpiを併用する場合」を補足しておく。
Intel MKLでBLACSのライブラリーを用いる場合、-lmkl_blacs_intelmpiか-lmkl_blacs_openmpiが選べる。ここでのコンパイルはopenmpiのコンパイラーであるmpifortもしくはmpif90で用いるのだが、intelmpiのライブラリを使うと、コンパイルは通るがコアダンプになることがわかった。
intelmpiはちゃんとIntel MPIを買って使いましょう。
散乱理論:罠1
散乱状態の漸近形が
とあるとき、散乱粒子が単位時間あたり表面要素を通過する確率は、
...ん?入射波はどこ行ったん??
っての確率密度だよね??
とずっとモヤモヤしていたが、ランダウの量子力学の脚注に何か書いてあるのを発見した。
「この議論では、散乱に関する実際の実験ではいつも満たされているように、入射粒子束は(回折効果を避けるために)幅広く、しかも有限のスリットによって仕切られていることを暗に意味している。このような訳で、両項間の干渉は存在しない;絶対値の自乗は入射波の存在しない点で取ることができる.」
スッキリしたようなしないような。。。
つまり、であれば、で入射波のないところに取れるから、散乱波だけでOKということか。
でものときどうするの?と思ってしまう。
こんな散乱理論の初っ端で一般性を落とさないといけないのだろうか?
っていうか、影散乱って散乱波と入射波の干渉で影が出来まっせ!って話だから、干渉を無くせますって言われても何かしっくり来ない。
あれか、中途半端に先を知ってしまったから影散乱とか気にしてしまうのかも知れない。
もう少し様子を見て、ランダウを読んでみることにした。
ちなみに持っているランダウの量子力学は恩師から譲り受けたものなので、いろいろ書き込みがあって面白い。
正直、恩師でも全部は読んでなかったことが判明してちょっと安心したりもしたりしている。
では、より散乱ライフを。
ブリルアンゾーンの対称性の高い点
ブリルアンゾーンなのかブリユアンゾーンなのかブリュアンゾーンなのか知りませんが、いっくら探してもG点とかX点とかM点とかK点とかR点とかの説明が無い。
みんなband計算でどうしてるの?と思うくらいに全く記述が無いに等しい。
逆格子ベクトルの定義とBZの立体図は死ぬほど見た。特にscとかbccとかfccとかhcpとか。
でも残念ながら計算したいのがmonoclinicなんですが?私に死ねと申すのですか?
Wikipediaの英語のBZのpageにはほぼ全部の構造が網羅されていると思われるBZの絵がある。
が、点の座標は無い。
図には以下のDOIのリンクへ行けとある。
10.1016/J.Commatsci.2010.05.010
行ってみたところ、Computational Material ScienceのW. Setyawana, S. Curtaroloの論文にぶち当たる。
残念ながら私の所属機関ではこの論文は落とせない。
と諦めかけていたら、arXivに落ちてた。
[1004.2974] High-throughput electronic band structure calculations: challenges and tools
神ってる、とは正にこのことだと悟った。
それでは、良いbandライフをお過ごし下さい。
バンド違いですが、昔ブルジュ・ハリファと砂漠の映像と北野武の語りが入るポカリスエットのCMがありましたが、そのときに流れてたtoeの"i dance alone"のアレンジがカッコ良過ぎたので、いつか音源化されないかなぁとボンヤリ思っていたりします。
研究という名の"i dance alone"に勤しみます。
多電子原子の波動関数(途中)
今日、朝起きてからずっと多電子原子波動関数の作成プログラムやっているけど、なかなかうまく行かない。
コードは引き続き下記のサイトを参照。
http://www-cms.phys.s.u-tokyo.ac.jp/~naoki/CIPINTRO/CIP/atom.html
とりあえずs軌道は大丈夫。
でもp軌道になった途端に、束縛状態がなくなってしまう。遠心力ポテンシャルに負けちゃう。
Liで結果を比較したけれども、イオン化エネルギーはバッチリ合っている。
でもトータルが一緒にならない。
全く同じコードを書いてるわけじゃないから、全く同じじゃないのは良いにしても、イオン化エネルギー合ってるのにトータルエネルギー合わないのはめっちゃ気持ち悪い。
てかp, d軌道を計算したい。。。
休みの間にケリがつかないのは気持ち悪いがしょうがない。
以下、コードです。
import numpy as np import matplotlib.pyplot as plt import prettyplotlib as ppl """ calculation condition """ dx = 0.0625 r_num = 255 + 1 r_far = r_num - 1 r_o = 160 output_r_max = 64.0 # 1s, 2s, 2p, 3s, 3p z = 3. deg = [ 2., 1. ] # degeneracy = occupancy n = [ 1, 2 ] # principal quantum number l = [ 0, 0 ] # partial wave eps = 1e-6 r, u = [], [] for i in range( r_num ): r.append( np.exp( dx * ( i - r_o ) ) / z ) u.append( - z ) """ functions """ # simga = r^2 | \phi |^2 # total sigma def Calc_sigma( ): # Gset and deg are global values sig = np.zeros( r_num ) for i_d, d in enumerate( deg ): # in this case, G[ r, ni ] is needed. NormCheck( Gset[ i_d ] ) for i_r, gi in enumerate( Gset[ i_d ] ): sig[ i_r ] += d * ( gi**2 ) sum = 0. for i_r, s in enumerate( sig ): sum += s * r[ i_r ] * dx print 'total electron', sum return sig def TotalEnergy( TE0 ): corre = 0. for i, sig in enumerate( sigma ): corre += sig * r[i] * ( 2. * Calc_uh( i ) + Calc_uxc( i ) ) * dx return TE0 - 0.25 * corre def Calc_uc( i ): return - z def Calc_uh( i ): uh = 0. for j in range( r_num ): if j < i : uh += sigma[j] * r[j] * dx # dr = \delta e^{\delta * x } dx = \delta r dx else: # i < j uh += sigma[j] * r[i] * dx return uh def Calc_rs( i ): return ( 6. * r[i] * r[i] / sigma[i] )**( 1. / 3. ) def Calc_beta( rs ): return 1. + 0.0368 * rs * np.log( 1. + 21. / rs ) def Calc_uxc( i ): rs = Calc_rs( i ) beta = Calc_beta( rs ) return - 2. * beta * ( 3. / ( 32. * np.pi * np.pi ) * r[i] * sigma[i] )**( 1. / 3. ) def Calc_u( i ): uc = Calc_uc( i ) uh = Calc_uh( i ) uxc = Calc_uxc( i ) return uc + uh + uxc def P( i, i_n ): return 2. * r[i] * u[i] - 2. * ( r[i]**2 ) * E + l[ i_n ] * ( l[ i_n ] + 1 ) def Pzero( i_n ): for i in range( r_far, -1, -1 ): if P( i, i_n ) < 0.: break if i == 0: i = r_o return i def Calc_dG( i, Gi, Fi, i_n ): return Fi def Calc_dF( i, Gi, Fi, i_n ): return P( i, i_n ) * Gi + Fi def HammingA( i_fit, i_n ): node = 0 Gpc, Fpc = 0., 0. G, F = np.zeros( r_num ), np.zeros( r_num ) dG, dF = np.zeros( r_num ), np.zeros( r_num ) for i in range( 0, 4 ): G[i] = r[i] ** ( l[ i_n ] + 1 ) F[i] = ( l[ i_n ] + 1 ) * G[i] dG[i] = Calc_dG( i, G[i], F[i], i_n ) dF[i] = Calc_dF( i, G[i], F[i], i_n ) for i in range( 4, i_fit + 1 ): Gp = G[i-4] + ( 4.0 * dx / 3.0 ) * ( 2. * dG[i-1] - dG[i-2] + 2. * dG[i-3] ) Fp = F[i-4] + ( 4.0 * dx / 3.0 ) * ( 2. * dF[i-1] - dF[i-2] + 2. * dF[i-3] ) Gm = Gp - ( 112. / 121. ) * Gpc Fm = Fp - ( 112. / 121. ) * Fpc dGm = Calc_dG( i, Gm, Fm, i_n ) dFm = Calc_dF( i, Gm, Fm, i_n ) Gc = ( 9. / 8. ) * G[i-1] - ( 1. / 8. ) * G[i-3] + ( 3. * dx / 8. ) * ( dGm + 2. * dG[i-1] - dG[i-2] ) Fc = ( 9. / 8. ) * F[i-1] - ( 1. / 8. ) * F[i-3] + ( 3. * dx / 8. ) * ( dFm + 2. * dF[i-1] - dF[i-2] ) Gpc = Gp - Gc Fpc = Fp - Fc G[i] = Gc + ( 9. / 121. ) * Gpc F[i] = Fc + ( 9. / 121. ) * Fpc dG[i] = Calc_dG( i, G[i], F[i], i_n ) dF[i] = Calc_dF( i, G[i], F[i], i_n ) if G[i] * G[i-1] < 0. : node += 1 return node, G, F def HammingB( i_fit, i_n ): node = 0 Gpc, Fpc = 0., 0. G, F = np.zeros( r_num ), np.zeros( r_num ) dG, dF = np.zeros( r_num ), np.zeros( r_num ) kapper = np.sqrt( np.abs( 2. * E ) ) for i in range( r_far, r_far - 4, -1 ): G[i] = np.exp( - kapper * r[i] ) F[i] = - kapper * r[i] * G[i] dG[i] = Calc_dG( i, G[i], F[i], i_n ) dF[i] = Calc_dF( i, G[i], F[i], i_n ) for i in range( r_far - 4, i_fit - 1, -1 ): Gp = G[ i + 4 ] - ( 4.0 * dx / 3.0 ) * ( 2. * dG[ i + 1 ] - dG[ i + 2 ] + 2. * dG[ i + 3 ] ) Fp = F[ i + 4 ] - ( 4.0 * dx / 3.0 ) * ( 2. * dF[ i + 1 ] - dF[ i + 2 ] + 2. * dF[ i + 3 ] ) Gm = Gp - ( 112. / 121. ) * Gpc Fm = Fp - ( 112. / 121. ) * Fpc dGm = Calc_dG( i, Gm, Fm, i_n ) dFm = Calc_dF( i, Gm, Fm, i_n ) Gc = ( 9. / 8. ) * G[ i + 1 ] - ( 1. / 8. ) * G[ i + 3 ] - ( 3. * dx / 8. ) * ( dGm + 2. * dG[ i + 1 ] - dG[ i + 2 ] ) Fc = ( 9. / 8. ) * F[ i + 1 ] - ( 1. / 8. ) * F[ i + 3 ] - ( 3. * dx / 8. ) * ( dFm + 2. * dF[ i + 1 ] - dF[ i + 2 ] ) Gpc = Gp - Gc Fpc = Fp - Fc G[i] = Gc + ( 9. / 121. ) * Gpc F[i] = Fc + ( 9. / 121. ) * Fpc dG[i] = Calc_dG( i, G[i], F[i], i_n ) dF[i] = Calc_dF( i, G[i], F[i], i_n ) if G[i] * G[ i + 1 ] < 0. : node += 1 return node, G, F def Normalization( GA, FA, GB, FB, i_fit ): for i in range( i_fit ): G[i] = GA[i] F[i] = FA[i] coeff = GA[ i_fit ] / GB[ i_fit ] for i in range( i_fit, r_num ): G[i] = GB[i] * coeff F[i] = FB[i] * coeff sum = 0. for i in range( r_num ): sum += ( G[i]**2 ) * r[i] * dx # dr = e^x dx = r dx norm = 1. / np.sqrt( sum ) for i in range( r_num ): G[i] *= norm F[i] *= norm return G, F, norm def NormCheck( vec ): sum = 0. for i_v, vi in enumerate( vec ): sum += r[ i_v ] * ( vi**2 ) * dx print 'Normalization check', sum def Solver( i_n ): true_node = n[i_n] - l[i_n] - 1 i_fit = Pzero( i_n ) nodeA, GA, FA = HammingA( i_fit, i_n ) LogA = FA[ i_fit ] / ( r[ i_fit ] * GA[ i_fit ] ) nodeB, GB, FB = HammingB( i_fit, i_n ) LogB = FB[ i_fit ] / ( r[ i_fit ] * GB[ i_fit ] ) node = nodeA + nodeB G, F, norm = Normalization( GA, FA, GB, FB, i_fit ) #NormCheck( G ) dE0 = - 0.5 * ( G[ i_fit ] ** 2 ) * ( LogB - LogA ) if not true_node == node: dE = float( true_node - node ) / ( n[i_n] ** 3 ) elif dE0 > 0.01: dE = dE0 * 0.5 else: dE = dE0 if not E < 0: print 'E', E + dE print ' not a bound state' exit() return node, G, F, dE """ program main """ dTE = 1e-5 Gset = [ np.zeros( r_num ) for i_n, ni in enumerate( n ) ] sigma = [ 1e-12 for i in range( r_num ) ] # initialize the global sigma u = np.array([ Calc_u( i ) for i in range( r_num ) ]) TE, TE_old = 0., 0. while abs( dTE ) > eps: for i_n, ni in enumerate( n ): E = - 0.5 * ( z**2 ) / ( ni**2 ) - 0.05 print '-' * 20 print ' ' * 5, 'n=', ni, 'l=', l[i_n] print '-' * 20 G, F = np.zeros( r_num ), np.zeros( r_num ) dE = 1e-5 if l[ i_n ] == 1: ppl.plot( r, [ P( i_r, i_n ) for i_r in range( r_num ) ] ) plt.ylim( -5, 5 ) plt.show() while abs( dE ) > eps: node, G, F, dE = Solver( i_n ) E = E + dE Gset[i_n] = G sigma = Calc_sigma() """ merge old and new versions of u """ u = np.array([ 0.6 * ui + 0.4 * Calc_u(i) for i, ui in enumerate( u ) ]) print deg[ i_n ], E TE += deg[ i_n ] * E dTE = TE - TE_old TE_old = TE TE = 0. print 'dTE', dTE print 'Total Energy0', TE_old TE = TotalEnergy( TE_old ) print 'Total Energy', TE for i_gs, gs in enumerate( Gset ): ppl.plot( r, gs ) plt.show() ppl.plot( r, map( Calc_u, range( r_num ) ) ) plt.show()
水素原子の波動関数の数値計算
束縛状態の波動関数の計算に抵抗がありまくりんぐなので、作ってみました。
普段は連続状態を計算していて、境界条件って言っても解析解につなぐだけなので、「正しいエネルギー(固有値)を見つける」というプロセスに心のハードルが無限の井戸型ポテンシャル状態でした。
コードは以下のサイトを参考に、pythonで書きました。
http://www-cms.phys.s.u-tokyo.ac.jp/~naoki/CIPINTRO/CIP/radschro.html#S_Orbital
自分で分かっていなかったポイントを整理すると、
- 原点周辺も無限遠周辺も、漸近解を初期値に使う。厳密な初期値かどうか気にするよりも、大雑把でも物理的要請を満たしている方がよっぽど大事。
- 漸近解を用いるのであれば、ルンゲクッタだろうが予測子修正子だろうが、初期値はただ漸近解の値を使えば良いから、運用上の違いはそんなに気にしなくて良い。もちろん誤差は違うだろうが。
- 最初に与える試行固有値は、やっぱり解析解をある程度は参照する必要があると思われる。例えば節の数が全然合わないときには真のエネルギー近傍に無理やり持って行ったりとかは知らないと出来ないなぁと。
- 試行固有値が真の固有値かどうかの判定として、原点近傍と無限遠近傍の両方から計算して、両者の解があるところで一次微分までで繋がるかどうかを見ることが使える。で、そのようなある点で微分不連続を与えるようなポテンシャルはデルタ関数に違いないということを考え、デルタ関数の分を一次摂動としてエネルギーから差っ引く。これで真のエネルギーに近づくでしょうという発想。
- 本来は真の波動関数があって、そこにデルタ関数があるから波動関数が滑らかにならず、真の波動関数で一次摂動を計算したものを差っ引くべきなんだが、ここで現在得られている不連続な波動関数を真の波動関数と読みかえ、それで近似的に一次摂動を計算している。なかなかアクロバットやなぁと感心する。
- 節を数えるのも結構苦手意識があったけど、単純に前後の波動関数の値を掛けて負になったら一個節があると思えば十分だった。
- 波動関数の規格化が長方形のやつを使っていて、普段スプライン積分で中身が半分ブラックボックスで使っている身としてはカルチャーショックだった。多分自分で一からやろうとしたら、「これで良いの?」と自問自答して前に進めなかったと思う。数値計算には漢気が大事だと学んだ。
自分で作ると、やっぱり感慨深いものがある。
試しに3s軌道が初期値からどのように修正されるかを見てみる。
以下、コードを示す。
import numpy as np import matplotlib.pyplot as plt import prettyplotlib as ppl z = 1. dx = 0.0625 r_num = 255 + 1 r_far = r_num - 1 r_o = 160 output_r_max = 64.0 r, u = [], [] for i in range( r_num ): r.append( np.exp( dx * ( i - r_o ) ) / z ) u.append( - z ) def P( i, l, E ): return 2. * r[i] * u[i] - 2. * ( r[i]**2 ) * E + l * ( l + 1 ) def Pzero( l, E ): for i in range( r_far, -1, -1 ): if P( i, l, E ) < 0.: break return i def Calc_dG( i, G, F, l, E ): return F def Calc_dF( i, G, F, l, E ): return P( i, l, E ) * G + F def HammingA( i_fit, l, E ): node = 0 Gpc, Fpc = 0., 0. G, F = np.zeros( r_num ), np.zeros( r_num ) dG, dF = np.zeros( r_num ), np.zeros( r_num ) for i in range( 0, 4 ): G[i] = r[i] ** ( l + 1 ) F[i] = ( l + 1 ) * G[i] dG[i] = Calc_dG( i, G[i], F[i], l, E ) dF[i] = Calc_dF( i, G[i], F[i], l, E ) for i in range( 4, i_fit + 1 ): Gp = G[i-4] + ( 4.0 * dx / 3.0 ) * ( 2. * dG[i-1] - dG[i-2] + 2. * dG[i-3] ) Fp = F[i-4] + ( 4.0 * dx / 3.0 ) * ( 2. * dF[i-1] - dF[i-2] + 2. * dF[i-3] ) Gm = Gp - ( 112. / 121. ) * Gpc Fm = Fp - ( 112. / 121. ) * Fpc dGm = Calc_dG( i, Gm, Fm, l, E ) dFm = Calc_dF( i, Gm, Fm, l, E ) Gc = ( 9. / 8. ) * G[i-1] - ( 1. / 8. ) * G[i-3] + ( 3. * dx / 8. ) * ( dGm + 2. * dG[i-1] - dG[i-2] ) Fc = ( 9. / 8. ) * F[i-1] - ( 1. / 8. ) * F[i-3] + ( 3. * dx / 8. ) * ( dFm + 2. * dF[i-1] - dF[i-2] ) Gpc = Gp - Gc Fpc = Fp - Fc G[i] = Gc + ( 9. / 121. ) * Gpc F[i] = Fc + ( 9. / 121. ) * Fpc dG[i] = Calc_dG( i, G[i], F[i], l, E ) dF[i] = Calc_dF( i, G[i], F[i], l, E ) if G[i] * G[i-1] < 0. : node += 1 return node, G, F def HammingB( i_fit, l, E ): node = 0 Gpc, Fpc = 0., 0. G, F = np.zeros( r_num ), np.zeros( r_num ) dG, dF = np.zeros( r_num ), np.zeros( r_num ) kapper = np.sqrt( np.abs( 2. * E ) ) for i in range( r_far, r_far - 4, -1 ): G[i] = np.exp( - kapper * r[i] ) F[i] = - kapper * r[i] * G[i] dG[i] = Calc_dG( i, G[i], F[i], l, E ) dF[i] = Calc_dF( i, G[i], F[i], l, E ) for i in range( r_far - 4, i_fit - 1, -1 ): Gp = G[ i + 4 ] - ( 4.0 * dx / 3.0 ) * ( 2. * dG[ i + 1 ] - dG[ i + 2 ] + 2. * dG[ i + 3 ] ) Fp = F[ i + 4 ] - ( 4.0 * dx / 3.0 ) * ( 2. * dF[ i + 1 ] - dF[ i + 2 ] + 2. * dF[ i + 3 ] ) Gm = Gp - ( 112. / 121. ) * Gpc Fm = Fp - ( 112. / 121. ) * Fpc dGm = Calc_dG( i, Gm, Fm, l, E ) dFm = Calc_dF( i, Gm, Fm, l, E ) Gc = ( 9. / 8. ) * G[ i + 1 ] - ( 1. / 8. ) * G[ i + 3 ] - ( 3. * dx / 8. ) * ( dGm + 2. * dG[ i + 1 ] - dG[ i + 2 ] ) Fc = ( 9. / 8. ) * F[ i + 1 ] - ( 1. / 8. ) * F[ i + 3 ] - ( 3. * dx / 8. ) * ( dFm + 2. * dF[ i + 1 ] - dF[ i + 2 ] ) Gpc = Gp - Gc Fpc = Fp - Fc G[i] = Gc + ( 9. / 121. ) * Gpc F[i] = Fc + ( 9. / 121. ) * Fpc dG[i] = Calc_dG( i, G[i], F[i], l, E ) dF[i] = Calc_dF( i, G[i], F[i], l, E ) if G[i] * G[ i + 1 ] < 0. : node += 1 return node, G, F def Normalization( GA, FA, GB, FB, i_fit ): for i in range( i_fit ): G[i] = GA[i] F[i] = FA[i] coeff = GA[ i_fit ] / GB[ i_fit ] for i in range( i_fit, r_num ): G[i] = GB[i] * coeff F[i] = FB[i] * coeff sum = 0. for i in range( r_num ): sum += ( G[i]**2 ) * r[i] * dx # dr = e^x dx = r dx norm = 1. / np.sqrt( sum ) for i in range( r_num ): G[i] *= norm F[i] *= norm return G, F def Solver( n, l, E ): true_node = n - l - 1 i_fit = Pzero( l, E ) nodeA, GA, FA = HammingA( i_fit, l, E ) LogA = FA[ i_fit ] / ( r[ i_fit ] * GA[ i_fit ] ) nodeB, GB, FB = HammingB( i_fit, l, E ) LogB = FB[ i_fit ] / ( r[ i_fit ] * GB[ i_fit ] ) node = nodeA + nodeB G, F = Normalization( GA, FA, GB, FB, i_fit ) dE0 = - 0.5 * ( G[ i_fit ] ** 2 ) * ( LogB - LogA ) if not true_node == node: dE = float( true_node - node ) / ( n * n * n ) elif dE0 > 0.01: dE = dE0 * 0.5 else: dE = dE0 E = E + dE print 'true_node', true_node, 'node', node print 'dE0', dE0, 'dE', dE print 'E', E ppl.plot( r, G ) if not E < 0: print 'E', E print ' not a bound state' exit() return node, G, F, dE eps = 1e-6 for n in range( 1, 5 ): E = - 0.5 * ( z**2 ) / ( n**2 ) - 0.05 for l in range( 4 ): if l >= n: continue print 'n=', n, 'l=', l G, F = np.zeros( r_num ), np.zeros( r_num ) dE = 1e-5 while abs( dE ) > eps: node, G, F, dE = Solver( n, l, E ) E = E + dE plt.xlim( 0, 50 ) plt.ylim( -0.8, 0.8 ) plt.show()
ミクロカノニカルとカノニカル
ただの愚痴であるが、位相空間が苦手である。
大抵、ミクロカノニカル分布から入ると思うが、その説明に
が使われる。
その中で一番状態数の多い状態が平衡に対応するとし、変分をかまし、最後にマクスウェル-ボルツマン分布を得る。
正直、これだけよくわからんものをホイホイ導入されて、これをとにかく納得しろという方が無理。そもそもイメージ出来ないし。
あと、これで持って来られるとカノニカル分布との対応が非常にわかりにくい。
てかカノニカル分布と同じノリじゃダメなの?と思う。
カノニカル分布は
- エネルギーで状態が指定される箱(系)を用意。
- 箱を複数連結させる。これを集団(アンサンブル)と呼ぶ。
- 箱の間でエネルギーのやり取り可能。
- 集団全体でエネルギーは固定。
- 集団のエネルギー = 箱1 + 箱2 + 箱3 +...
- 全部の箱のエネルギーが一定とは全く限らないから、箱同士のエネルギーの奪い合いが、今、始まる!
的な感じで、
全体のエネルギーが保存する様に各箱にエネルギーを割り振ったときに、どんな風にエネルギーを割り振るのが一番可能性が高いか?
となる。
「箱の間のエネルギーのやり取り可能」と書きましたが、正直ここは僕はこの表現嫌いです。
「やり取りがあるからエネルギーが割り振れる」という立場からすると、必要と感じますが、やり取りを真面目に考えると箱のエネルギーを独立に足せないし、割り振っている間のエネルギー移動のプロセスは全く気にしないので、個人的には、「箱同士は独立に扱い、それらのエネルギー状態は確率的に発生する」とした方がしっくりくる。下手に物理的に扱うと変なことになる気がする。
ちなみに、全部の箱にエネルギーを均等に割り振る場合の数は「1通り」しかないから、ほとんど発生しない。
箱1のエネルギーがちょっと高くて、その代わりに箱2のエネルギーがちょっと低いのは、その逆もあるし、箱2じゃなくて箱3が低くても良いから、こういう方が場合の数が多くなる。そして結果的にエネルギー分布が出来るのがわかる。
それでエネルギーの保存とか粒子数の保存とかの制限が入って、カノニカル分布と言われるものが得られる。
これをミクロカノニカルにそのまま転用すれば、
- エネルギーで状態が指定される粒子を用意。
- 粒子を箱の中に入れる。
- 粒子間でエネルギーのやり取り可能
- 箱でエネルギーは固定。
- 箱のエネルギー = 粒子1 + 粒子2 + 粒子3 +...
- 粒子の箱のエネルギーが一定とは全く限らないから、粒子同士のエネルギーの奪い合いが、今、始まる!
で良くない?
ここでもカノニカルのときと同じで、粒子は独立だけど確率的にエネルギー状態が発現して、箱でエネルギーが保存していなければならないから、全ての粒子のエネルギーが等しいのではなく、マクスウェル-ボルツマン分布に従う形になる。
逆に、粒子が独立じゃなくて箱のエネルギーを各粒子の和で書けないから、いっそ箱単位の扱いに拡張しようというのがカノニカル分布だと理解した。
例えば「クーロン相互作用のある電子の集団」を考えたとき、ハートリー・フォック近似で独立粒子近似をしたとしても、全エネルギーは各粒子の和にはならない。そのまま足すと相互作用を二重に数えてしまう。
何となく、バンド計算とかを齧っていると、
ミクロカノニカル → カノニカル
unit cell → supercell
な対応に近いなぁとも最近思う。
例えばunit cellに原子一個だけとかだと、原子一個解いてあと並進対称性でズルすればよいけど、不純物が入っているとかなんとかでsupercellにすると、その大きいセルの中で真面目に解かないといけなくなる。
一応、言及しておくが、等重率の原理は使っている。場合の数を数えるときにはすべて同じ重みで扱う。
更に言及しておくが、多分、田崎先生の統計力学に似たようなことが書いてあったと思う。
今、手元にないので、割と適当なことを言ってしまっていると思う。
でも、阿部大先生の本の1章読んでて、やっぱり位相空間から入るのは性に合わないと感じた。