分子のおもちゃ箱 bbs 09 No.901...1000 ( 2004.07.04 - 2004.08.31 )
 ( prev | index | next )

●1000 分子のすり抜け Nori - 2004/08/31 16:29 -
書き忘れたのですが、あと上部に配置した原子が結晶を通り過ぎるといった現象がおきました。
相互作用はきいているようなのですが…
最初の結晶構造に隙間があるので、そこから入り込むのでしょうか?
もしくはポテンシャルがおかしいのでしょうか?
たびたび申し訳ありません。

●999 No313no Nori - 2004/08/31 16:27 -
アドバイスありがとうございます。
周期境界条件を入れてみたのですが、やはり結晶構造が崩壊してしまうので、結晶を構成する原子についてはベルレ法を使って位置を求めているところで、
xx[i]=xx[i]と位置を更新しないようにしました…

こうするともちろん結晶構造は維持できるのですが、右側二表示される温度と圧力がすごい数字になってしまいます。
この原因は何でしょうか?
動径分布関数についても同じでグラフが変な形になります。
動径分布関数とは粒子の周りにどれだけの粒子があるかというのを表す関数だったの思うのですが、このアプレットでは、どういうアルゴリズムで、どういった意味があるのでしょうか?
グラフの意味を教えていただきたいです。

●998 アプレット追加 mike - 2004/08/30 17:49 -
アプレット追加しました。

(525) periodic SD3D 水素様原子(3次元-最急降下法-周期的境界条件)
周期的境界条件をつけ、最急降下法とGram-Schmidtの直交化法により水素様原子の固有値と固有関数を求めます。
理論値はε(1s)=-0.5au、ε(2s)=ε(2p)=-0.125auです。箱が有限な場合(アプレットNo.513参照)とは逆に、
2sのエネルギーの方が低くなるようです。

●997 アプレット追加 mike - 2004/08/29 18:17 -
アプレット追加しました。

(523) periodic SD1D 一次元最急降下法(周期的境界条件)
最急降下法とGram-Schmidtの直交化法により、周期的境界条件のもとで一次元ポテンシャル(緑)中の固有値と
固有関数を求めます。最急降下法は、H = -Δ/2 + V(x)として、関数のセット{ ui(t) }から、
 μ{ui(t+dt)-ui(t)}/dt = - (H-εi) ui(t), εi= <ui(t)|H|ui(t)>
から次世代の、{ ui(t+dt) }を求め、これを直交化することをくり返し、収束したεiと{ui}を求める方法です。
V=0(depth=0)とすると、{1,cos(2nπ/L),sin(2nπ/L)}のフーリエ級数の基底関数のセットが得られます。
特筆すべきは、V=0のとき、一様に拡がった(u(x)=1/sqrt(L))解も存在でき、エネルギーは0となります。

(524) periodic SD2D 円形の井戸(2次元最急降下法-周期的境界条件)
周期的境界条件のもとで、最急降下法と直交化により2次元の円形井戸の固有値εiと固有関数{ui}を求めます。
最急降下法は、H = -Δ/2 + V(r),V(r)=-depth(r<radious),0(other)として、関数のセット{ui}から
 μ{ui(t+dt)-ui(t)}/dt = - (H-εi) ui(t), εi= <ui|H|ui>
により次世代の、{ ui(x,y,t+dt) }を求め、直交化することをくり返し、収束したεiと{ui}を求める方法です。
井戸を浅くしていくと、電子は、局在から遍在へと変わり、depth=0のとき、2次元の自由粒子となります。

●996 Re 結晶の分子動力学 mike - 2004/08/29 06:46 -
Noriさん、こんにちは。

>先日お尋ねしました、第3項のパラメータ(D1,B1,D2,B2)なのですが、他の文献を探してもどうしても見つからなかったので、
>河村先生に直接お伺いしました。
>DはkJ/molで、BはÅ^-1のようです。
すばらしいですね。 情報、ありがとうございます。

>そこで、早速ポテンシャルを変更して、実行してみたのですがどうにもうまくいきません。
>初期配置で設定した結晶構造がばらばらになってしまい、No313の
>アプレットのように中央付近にかたまってしまいます。
>・・・(中略)・・・
>なぜタングステンへの吸着のアプレットではあのように結晶構造がくずれないのでしょうか? 力がつりあっているからでしょうか?それとも何か拘束条件のようなものをいれているとか…

思いあたる原因は2つあります。(他にも原因があるかどうかわかりませんが、とりあえず...)
(1)表面張力:はじめの結晶が直方体のような形の場合、頂点の原子が周りから受ける力(表面張力)がものすごく大きくなります。このため、丸くなろうとすることがあります。たとえば、z方向に垂直な表面を作りたい時、x、y方向には周期的境界条件を付けるのがいいと思います。
(2)結晶構造(分極):NaCl構造では周期的境界条件を付けなくても直方体は維持できますが、CaF2構造では周期的境界条件を付けないと結晶構造を維持できません。これは、CaF2構造では、単独では分極ができてしまい、小さな結晶では、これが余分なエネルギーが大きくなるための考えられます。
(CaF2についてはNo.150を参照ください)

>それから最後のForce()で行われていることの意味がよくわかりません。ここではただ値を返すだけではなく、簡単な計算をしているようですが、これにはどういった意味があるのでしょうか?
線形補間を行っています。表を2回引くことで、高い精度の表引きが可能になります。
x[i]<x<x[i+1]のとき、xに対する値y(x)が求めたいとすると、表を2回引き、x[i]のときy[i]、x[i+1]のときy[i+1]として、y=a*y[i]+(1-a)*y[i+1]、a=x[i+1]/(x[i+1]-x[i])から求めます。

●995 Re^3 プログラム借用のお願い mike - 2004/08/29 06:44 -
iotsukaさん、こんにちは。わざわざ、ありがとうございます。
かなり異なった分野ですが、今後とも、よろしくお願いします。

●994 Nori 結晶の分子動力学 - 2004/08/28 22:59 -
遅くなりましたが、ご返答ありがとうございます。

先日お尋ねしました、第3項のパラメータ(D1,B1,D2,B2)なのですが、他の文献を探してもどうしても見つからなかったので、
河村先生に直接お伺いしました。
DはkJ/molで、BはÅ^-1のようです。

そこで、早速ポテンシャルを変更して、実行してみたのですがどうにもうまくいきません。

初期配置で設定した結晶構造がばらばらになってしまい、No313の
アプレットのように中央付近にかたまってしまいます。
私のイメージとしては、タングステンへの吸着のアプレットのように、下の結晶構造は崩れずに上の原子だけが動くようにしたいのですが・・・

原因としては、ポテンシャルパラメータがおかしい、計算アルゴリズムが間違えているなどの可能性がありますが、もしよろしければアドバイスをいただけないでしょうか?

なぜタングステンへの吸着のアプレットではあのように結晶構造がくずれないのでしょうか? 力がつりあっているからでしょうか?それとも何か拘束条件のようなものをいれているとか…

それから最後のForce()で行われていることの意味がよくわかりません。ここではただ値を返すだけではなく、簡単な計算をしているようですが、これにはどういった意味があるのでしょうか?
よろしければ教えていただきたいです。

●993 Re: Re プログラム借用のお願い iotsuka - 2004/08/28 18:28 -
iotsukaと申します。お世話になっております。

>>つきましては、(中略)当方が借用公開いたしてよろしいか、判定していただければと思います。
>OKです。お役にたてれば、幸いです。
借用OKして下さり、どうもありがとうございます。とても助かります。

当方のページの内容が変わって、借用公開の条件が変化した場合は、また適宜見ていただいて、公開して良いか、適宜判定していただければと思います。

それでは、今後ともよろしくお願いいたします。

http://iwao-otsuka.com/dw/indexdw.htm

●992 アプレット追加 mike - 2004/08/28 16:55 -
アプレット追加しました。
(522) alkali empty-core pseudopotential アルカリ金属の空芯-擬ポテンシャル
空芯モデル擬ポテンシャル*) を用いて実空間-密度汎関数法によりアルカリ金属原子の価電子の構造を再現します。
アルカリ金属原子(Na,K,Rb,Cs)は稀ガス閉殻の外側のs軌道に、1個の価電子を持っています。
空芯モデルは、W(r)=-2/r (r>Rc)、W(r)=0(r<Rc)であり、パラメータはRc(およそイオン半径)のみです。
Rcとイオン化エネルギーは、Na:Rc=1.7au, Vi=0.162au(実測値0.189)、K:Rc=2.2au, Vi=0.134au(0.160)、
Rb:Rc=2.45, Vi=0.123au(0.154)、Cs:Rc=2.7,Vi=0.112au(0.143)と少し低めになります。
*) Chekmarev, et.al., J. Chem. Phys. 109, p.768-778 (1998)

●991 Re プログラム借用のお願い mike - 2004/08/28 16:45 -
iotsukaさん、こんにちは。ご連絡、ありがとうございます。

>つきましては、上記の中身を実際に見ていただいて、当方が借用公開いたしてよろしいか、判定していただければと思います。

OKです。お役にたてれば、幸いです。

上記の公開されているページを拝見しました。ドライとウェットという視点は新鮮で、興味深く拝見させていただきました。分子運動と社会科学的な方面はあまり縁がなさそうに思っていたのですが、何か関連があるんですね。

今後とも、よろしくお願いします。

●990 プログラム借用のお願い iotsuka - 2004/08/28 06:47 -
iotsukaと申します。初めまして、よろしくお願いいたします。

この度、管理人様が作られた液体、気体分子運動のプログラムを借用、一部改変する形で、管理人様の名前と共に当方のwebサイトで公開したいと考えております。

具体的には、以下のような形で載せたいと考えております。
(現在、仮公開中。)
http://iwao-otsuka.com/enqb/res/gasliquidcmp1.htm

つきましては、上記の中身を実際に見ていただいて、当方が借用公開いたしてよろしいか、判定していただければと思います。

借用公開不可と判定された場合は、上記ページは直ちに削除します。

以上、よろしくお願いいたします。

http://iwao-otsuka.com

●989 アプレット追加 mike - 2004/08/27 20:26 -
アプレット追加しました。
(521) Hg empty-core pseudopotential 水銀の空芯擬ポテンシャル [ home]
Hgの空芯モデル擬ポテンシャル*) を用いて実空間-密度汎関数法によりHg原子の価電子の構造を再現します。
Z=1.5としてイオン化エネルギーを求めると、0.340au(実測値0.383)となり、少し低めになります。
Hgの電子配置はZ=+80の原子核と電子80個からなる多体系ですが、5d軌道まで閉殻(Au殻)でその外側6sに、
2個の価電子を持っています。空芯モデルでは、閉殻イオンは外から(r>Rc)は、V(r)=-2/r として見え、
殻内(r<Rc)ではV=0とする近似で、粗い近似のわりには半定量的な取扱いができると言われています。
*) Chekmarev, et.al., Phys. Rev. E59, p.479-491 (1999) ; Hgに対して、Rc=0.915auです。

●988 アプレット追加 mike - 2004/08/25 18:57 -
アプレット追加しました。
(520) light atom RSDFT 軽い原子(実空間-密度汎関数法)
実空間-密度汎関数法を用いてH, He, Li, Be, B, C原子を構成します。電子数n=Z-0.5として最外軌道の
エネルギーからイオン化エネルギーを求めると、B:0.259au(実測値0.305)、C:0.709au(0.414)とだんだん
隔たってきます。本アプレットではメッシュはdx=0.25auであり、1sの径方向のrR(r)の最大値は水素様のとき
r=1/Zですので、空間分解能が足りなくなっています。空間分解能を上げないまま、実用的な計算するには、
Bより重い元素では擬ポテンシャルの使用が必須と考えられます。

●987 アプレット追加 mike - 2004/08/22 18:57 -
アプレット追加しました。
(519) Be like atom RSDFT ベリリウム原子(実空間-密度汎関数法)

●986 アプレット追加 mike - 2004/08/21 18:43 -
アプレット追加しました。
(518) H2 molecule RSDFT 水素分子(実空間-密度汎関数法)

核間距離dpp=1.5auのとき、 全エネルギーは最小となり、Etot=-1.12auとなります。
実測値はdpp=1.4au、Etot=-1.17auです。精度はあまり良くないですが、ほぼ、
水素分子を再現できていると思います。No.511と共通なルーチンが多く、比較すると、
No.511の核間距離が合わなかった原因が、小林先生が指摘された「十分に収束していない
電子密度を力の評価に使った」ためだということが検証されたと考えます。

小林先生>適確なアドバイス、ありがとうございます。

●985 アプレット追加 mike - 2004/08/19 18:48 -
アプレット追加しました。
(517) He like atom - fine RSDFT ヘリウム様原子(細かいメッシュ)

●984 アプレット追加 mike - 2004/08/18 20:21 -
アプレット追加しました。
(516) He like atom - large RSDFT ヘリウム様原子(大きい箱)

●983 Re 982 モデルのパラメータ mike - 2004/08/18 18:27 -
Noriさん、こんにちは。

>f0についてですが、f0=1kcal/(mol・Å)という記述がありました。
私は、どこにあったか見つけられなかったのですが、上記の単位とするとSI単位になおすと
f0=kcal/(mol・Å)=4185J/(6.e23個・1.0e-10)=6.975e-11となります。

>ε0=8.8542e-12 でいいのでしょうか?
SI単位ではそれでいいです。

>f0=6.9742e-11でしょうか?
たぶんf0の部分に相当すると思います。

>このポテンシャルの引用先をよろしければ教えていただきたいです。
河村先生の用いられているポテンシャルはクーロン力を近接力と遠くのエバルト和に分解して評価するものですが、この方法は周期的な境界条件のときはいいのですが、孤立系では適用がむずかしいと考えられます。このため、スクリーニングを行う方法(「初心者の..」のp.53とp54参照)にしました(ただ、この方法は電気的中性が保たれているのが必要があります)。 最も一般的なクーロン力の計算方法は多段階の多重極子展開を用いる方法のようですが、アルゴリズムがむずかしく、手がついていません。本アプレットではクーロン項をscreened-Coulombポテンシャルとし、イオン反発項はBorn-Mayer型とし、パラメータはSX-1から採っています。
これで、組んだ結晶の格子定数や膨張率がほぼ合うことを確認しています。

>また、研究報告にあったポテンシャルの第3項のパラメータ(D1,B1,D2,B2)なのですが、単位がわかりません。
たしかに、この単位はわかりませんね。....
Si-Oの結合エネルギーを調べて、これとDの値を比べ、あたりを付けてみてはいかがでしょうか?

●982 SX-モデル、パラメータの単位 Nori - 2004/08/18 17:17 -
返答ありがとうございます。

f0についてですが、f0=1kcal/(mol・Å)という記述がありました。
No313のアプレットのForceTable()で力の計算をしていますが、
eForceConst = 1.0/(4.0*pai*8.8542e-12) になっており、
これは河村先生の研究報告にあった4πε0で割っているところと同じかと思うのですが、ε0=8.8542e-12 でいいのでしょうか?
またf0についても6.9742e-11*Math.exp((a-r)/b)で、
研究報告に記述されている、ポテンシャルの微分の第2項に相当します。
ということは、f0=6.9742e-11でしょうか?
このポテンシャルの引用先をよろしければ教えていただきたいです。

また、研究報告にあったポテンシャルの第3項のパラメータ(D1,B1,D2,B2)なのですが、単位がわかりません。
共有結合を表すとあるので、Morse項として、DをEE、Bを(/AA)としたのですが、うまくいきません。
第1、第2項についてはmikeさんのアプレットの単位を用いさせていただいています。
単位については何も記述がないので、これらのパラメータの単位をご存知でしたら教えていただきたいです。

●981 アプレットNo.313について mike - 2004/08/17 18:07 -
Noriさん、こんにちは。

>・setRadialDistribution()の役割について。
これは原子の配置xx[][][]、yy[][][]、zz[][][]から動径分布関数を計算します。
結果をrdf[0][i]にセットします。

>・initVirial()、setVirial()の役割について。
系の圧力を計算するためのメソッド群です。「初心者の...」のp.29を参照してください。

>本のほうにのっているポテンシャルもカオリナイトのO-H間を表すものとして、適当といえるのでしょうか?
>この二つの違いは、4πε0で割っているところと、Morse項が微妙に違うところなのですが、これらの違いの意味はどういったものでしょうか?
報告書の方がより精密な近似といえると思います。しかしながら、具体的にはわかりません。

>河村先生のポテンシャルを用いようとしているのですが、f0、ε0をどのような値なのか、思案中です。
>本にはf0は単位を整合させるための定数とありますが、これらはどのような値なのでしょうか?
f0は河村先生の本のどこかに書かれていたと思いますが、今探したら見つかりませんでした。

●980 No313のアプレット Nori - 2004/08/17 16:15 -
ご指摘ありがとうございます。
さっそくNo316のアプレットを拝見させていただきました。
わからない箇所が数点ありましたので、質問させていただいてもよろしいでしょうか。
疑問点は以下の点です。
・setRadialDistribution()の役割について。
・initVirial()、setVirial()の役割について。

また、ここで、用いてるポテンシャルは「パソコン分子シミュレーション」のP20〜P25に述べられているSX-1モデルだと思うのですが、これは河村先生の研究報告にあったポテンシャルと非常によく似ていますね。
河村先生の著書なので、当たり前なのですが…。

本のほうにのっているポテンシャルもカオリナイトのO-H間を表すものとして、適当といえるのでしょうか?
この二つの違いは、4πε0で割っているところと、Morse項が微妙に違うところなのですが、これらの違いの意味はどういったものでしょうか?

河村先生のポテンシャルを用いようとしているのですが、f0、ε0をどのような値なのか、思案中です。
本にはf0は単位を整合させるための定数とありますが、これらはどのような値なのでしょうか?

多くの質問ですいません。 わかる範囲でけっこうですので、教えていただきたいです。

●979 表関数について mike - 2004/08/16 19:00 -
Noriさん、こんにちは。

>setForceTable()について詳しく教えていただけないでしょうか。
>表に代入しているということですが、なぜ、100や6000で場合わけをしているのかがちょっと分かりません。
表関数はポテンシャルの関数を計算する代わりに、予めこの関数を計算しておいて配列に記憶しておき、
関数計算の代わりに配列を参照するだけにして、高速化するのがねらいです。
アプレットNo.326のsetForceTable()は配列forceTable[i][j][ir]にi種とj種の間の距離irの時に働く力を記憶します。irは1.0e-4 nmを単位とした整数とします。ir=10000は1nmに相当します。
ir<100はr<0.01nmの場合で、このあたりの数値が使われることはほとんどないのですが、r=0で数値が無限大になるのを防ぐため、f=0.01nmの時の値を返すようにしています。
6000<ir<10000ではカットオフを行います。粒子登録法ではrがある値(本アプレットでは1nm)以上でポテンシャルが0になる必要があります。本例ではir=6000のとき1で、ir=10000に向かって直線的に0となるカットオフ関数scをもとのポテンシャルにかけることでこれを実現しています。

>またregNear()、regProximity()、registrationMorse()についても教えていただけないでしょうか。
これらは似たような処理をしていますが、どのように違うのか、粒子登録法の箇所を読んでもちょっと分かりません。
以前にも書きましたが、No326はFSとMorseが併存するため、かなり見苦しいプログラムとなっています。
河村先生のポテンシャルを用いようとされるなら、むしろNo.313(イオン結合のMDテンプレート)の方がお勧めです。こちらをご覧ください。

●978 ソースについての質問 Nori - 2004/08/15 19:36 -
たびたび申し訳ありませんが、setForceTable()について詳しく教えていただけないでしょうか。
表に代入しているということですが、なぜ、100や6000で場合わけをしているのかがちょっと分かりません。

またregNear()、regProximity()、registrationMorse()についても教えていただけないでしょうか。
これらは似たような処理をしていますが、どのように違うのか、粒子登録法の箇所を読んでもちょっと分かりません。

またご指摘のように、河村先生のポテンシャルを用いようと思うのですが、そうすると、regProximityなどは必要なくなるのでしょうか?
もしよろしければ教えていただきたいです。

●977 アプレット追加 mike - 2004/08/15 15:28 -
アプレット追加しました。
(515) He like atom - real space DFT ヘリウム様原子(実空間-密度汎関数法)

●976 アプレットNo.511について mike - 2004/08/15 15:24 -
小林一昭先生のご指導にしたがい、いくつかの検討を行い、次のような結論になりました。
本アプレットで採用しているlossプロセスにより、定常解に近付けるには、かなりの繰り返しが必要であり、不十分なまま、力の計算をしたため、実際と異なった結果になったと考えられます。

lossプロセスは波束の位相をそろえて動きを止め、運動エネルギーを奪うのですが、このとき、ビリアルを考えると、2<K>と<V>のバランスが崩れるため、定常解と比べると<K>の比率が異常に小さくなります。<K>は<ψ,(-Δ/2)ψ>なのでこれがちいさいことは、ψの広がりが定常解に比べ大きいことを意味し、これによって核にはたらく力も定常解と異なってしまったようです。

したがって、本アプレットは失敗作ということになります。
しかし、削除せず、事例として残そうと思います。

今後、検証の意味でも、SD法で水素分子の力計算、MDを行いたいと思います。

●975 リンク追加 mike - 2004/08/15 07:40 -
neoさんの「OpenGLで分子動力学シミュレーション」を追加しました。
今後の発展がたのしみです。

●974 ありがとうございます neo - 2004/08/13 12:59 -
>もちろん、OKです。こちらから、リンクさせていただいてもよろしいでしょうか?
ありがとうございます。もちろん大丈夫です。
私は量子化学、MDをやっていたのですが、
最近密度汎関数法を勉強しています。
mikeさんと 小林一昭先生の話はとても興味をもって
拝見させていただいています。
今後ともよろしくお願いします。
http://www.geocities.jp/hide_morita2002/

●973 Re はじめまして mike - 2004/08/11 18:21 -
neoさん、はじめまして。「分子のおもちゃ箱」をご覧いただき、ありがとうございます。

ホームページ「OpenGLで分子動力学シミュレーション」拝見しました。
今後予定されている内容は、とても興味深いものですね。
分子動力学を共に学んでいけたらと考えます。

>リンクさせていただきたいのですがよろしいでしょうか?
もちろん、OKです。こちらから、リンクさせていただいてもよろしいでしょうか?
分子動力学の横のつながりができたら素敵ですね。

>できましたらWindowsで作成したMDを見てもらいたいのですが
>送ってもいいでしょうか?
私はMacユーザーなので、あまりwindowsのことがわかりません。すみません。
貴ホームページに公開され、それを拝見できれば、いいのではないかと思います。

●972 はじめまして neo - 2004/08/11 17:09 -
はじめして。いつも拝見させていただいています。
私もMDをWindowsで作成しています。
今度新しくホームページを作りました。そこで、
リンクさせていただきたいのですがよろしいでしょうか?
できましたらWindowsで作成したMDを見てもらいたいのですが
送ってもいいでしょうか?

http://www.geocities.jp/hide_morita2002/

●971 アプレット追加 mike - 2004/08/09 18:03 -
アプレット追加しました。
(514) H like atom coarse SD3D 水素様原子(粗メッシュ 最急降下法)

●970 Re 河村先生の論文 mike - 2004/08/09 18:01 -
Noriさん、こんにちは。

ソースに関しては、我流で書いていて、わかりにくいと思います。なんでもお尋ねください。
河村先生の論文については、いっしょに考えていけたらいいと思いますので、遠慮なく、
なんでも書き込んでください。

●969 ご返答 Nori - 2004/08/08 21:59 -
ありがとうございます。

ソースをさらにじっくり見るとともに河村先生の論文についてもう少し詳しく調べ、自分の研究に活かせるようがんばっていきたいと思います。
またわからない点などありましたら、お尋ねすることがあるかもしれませんがよろしくお願いいたします。

●968 アプレット追加 mike - 2004/08/08 16:40 -
アプレット追加しました。
(512) harmonic oscillator SD2D 2次元調和振動子(最急降下法)
(513) H like atom SD3D 水素様原子(最急降下法)

●967 Re: 511について6 (後半) mike - 2004/08/08 08:42 -
(後半、続きです)

(補足3)ヘルマン-ファインマン力
>時間発展DFT-LDAの計算を見ていると、lossという操作でエネルギーを引いて系を
>収束させているように見えるのですが、収束が十分でない時は、
>上記のようにあまりH-F力は正確でない可能性があります。
ご指摘のように、lossによる収束を十分行うことを試みて見ます。

(補足4)力のチェック
先生のサイト、拝見しました。今、作っているアプレットが完了したら、力の
計算を試みようと思います。

●966 Re: 511について6 (前半 mike - 2004/08/08 08:41 -
小林一昭先生、たいへん詳しい、ご指導、ありがとうございます。
レスが遅れまして、申し訳ありません。

(補足1)井戸型ポテンシャル導入について
>昨日の書き込みでも少し書きましたが、周期的境界条件を課していない実空
>間の計算で、系に強電場を印加して電子を放出させる方法も一つの手かと思い
>ます。この時、基盤側部分の電荷中性を保つようにするかしないによる影響を
>調べることも面白そうです。この場合も励起状態云々の問題はありますが、、、。
今のところ、まだ準備はできていないのですが、渡辺-塚田先生の時間依存のKohn-Sham方程式の
時間発展アルゴリズムを用いて試みようかなと思っています。

>>原子を高温にして熱電子放出が再現できればと考えます。
> この場合、電子に関しても温度の概念が必要となるのですが、密度汎関数法は、
>絶対零度(←基底状態)が前提の計算なので、有限温度への拡張が必要となります。
>DFT+LDAを有限温度へ拡張 するアプローチはありますが、これもまだこれからという分野です。
B-O近似で電子の定常状態を考えることは、絶対零度になるということでしょうか。
まだ、全く知識がありませんので、調べていこうと思います。

(補足2)原子の計算での動径方向のメッシュ
対数メッシュの方法、わかりました。ありがとうございます。
アプレットNo.493は0-32auの等間隔メッシュ、dr=1/16.0において、rRに変数変換し、
2分法で解いているのですが、1sの固有値は-0.4995au とかなり理論値に近い値を与えるようです。

>実空間法でも、”adaptive coordinates”という上記に対応したような手法 があるようですが、
>詳細は筆者もよく知りません。以下に、”adaptive coordinates”というキーワードで把握できた
>文献[1]を示します。実際は、実空間法関連で(検索エンジンや、文献データベースなど)探せば
>この手法に関 しての情報は割と得られるかと思います。
>[1] D. R. Hamann, Phys. Rev. B63, 075107(2001)
ありがとうございます。上記の文献そのものは入手できなかったのですが、googleで
”adaptive coordinates Hamann”の検索し、下記*)のpdfファイルを入手しました。
まずは、これを読んでみようと思います。

*) arXiv: cond- mat/ 9506086 v1 19 Jun 1995
Gil Zumbach, N. A. Modine and Efthimios Kaxiras;
”Adaptive coordinate, real-space electronic structure calculations on parallel computers”

●965 Re: 511について5(後半) 小林一昭 - 2004/08/06 23:12 -
以下、後半部分です。

(補足3)ヘルマン-ファインマン力
ヘルマン-ファインマンの力(以下、H-F力)は、電子状態に関する変分計算
(自己無撞着な計算)が収束していない場合や完全系をなす基底関数が十分で
ない場合には、正しい力とならなくなります(→プーレイ補正)。

時間発展DFT-LDAの計算を見ていると、lossという操作でエネルギーを引い
て系を収束させているように見えるのですが、収束が十分でない時は、上記の
ようにあまりH-F力は正確でない可能性があります。

実空間で上記の基底の数に相当するのが、メッシュの数(細かさ)で、これ
が十分でないとH-F力はそのままでは正しくなくなる可能性があります。

プーレイ補正項に関しては、

http://www.geocities.co.jp/Technopolis/4765/INTRO/yogo.html#DK000010

及び、そこにある参考文献等もご参照下さい。

(補足3)力のチェック
ずれΔrとそれによる全エネルギーの変化ΔEtotによって得られる力、F = -
ΔEtot/Δrと、ヘルマン-ファインマン力の表式から得られる力との比較に関
しては、。

http://www.geocities.co.jp/Technopolis/4765/INTRO/yogo.html#DK000014

もご参照ください。

http://www.bandstructure.jp/

●964 Re: 511について5(前半) 小林一昭 - 2004/08/06 23:11 -
いくつか気付いた点について以下に示します。

(補足1)井戸型ポテンシャル導入について
周期的な系(スラブモデル)で、井戸型のポテンシャルの導入により、スラ
ブ(表面基盤)側の電子が移動、捕獲された場合、スラブ側が正、井戸側が負
(←捕獲、局在した電子による)の電荷に帯電し、双極子を形成します。これ
が(特に周期系において)何か悪さをしないか心配です。

またこれは正に時間依存の問題として扱うべきものですが、井戸に捕獲され
た電子状態を含む系全体は、純然たる励起状態で密度汎関数法+局所密度近似
が、どこまでこれを正しく記述できるかは何とも言えません。おそらく定量的
な精度はそのままではあまり期待できないと思われます。

昨日の書き込みでも少し書きましたが、周期的境界条件を課していない実空
間の計算で、系に強電場を印加して電子を放出させる方法も一つの手かと思い
ます。この時、基盤側部分の電荷中性を保つようにするかしないによる影響を
調べることも面白そうです。この場合も励起状態云々の問題はありますが、、、。
電場の印加に関しては、STM絡みとかで、理論計算が試されている可能性が
あります。この分野は筆者はあまりよく知らないので詳細はご勘弁下さい。

>原子を高温にして熱電子放出が再現できればと考えます。

 この場合、電子に関しても温度の概念が必要となるのですが、
密度汎関数法は、絶対零度(←基底状態)が前提の計算なので、
有限温度への拡張が必要となります。DFT+LDAを有限温度へ拡張
するアプローチはありますが、これもまだこれからという分野
です。

(補足2)原子の計算での動径方向のメッシュ
筆者の使用している原子の計算では、動径方向に関してポテンシャルの原点
近傍はより細かいメッシュ、遠方はより間隔の広いメッシュをとります。これ
を対数メッシュ(以下、log-mesh)と言います。

(定義)
r(log-mesh)=rmax*(exp(h))**N/(exp(h)**421)

r():動径方向座標、exp:指数関数、N:1〜421までの数
rmax=60.0a.u., h=1/32.0

です。N=421(動径方向のメッシュ数とする)のときr(log-mesh)=rmaxとな
ります。この定義から、原点近傍に行くほどメッシュ間の間隔は狭くなり、原
点から遠方になるほど間隔は広くなっていきます。

これは、原点近傍では1/rの変化が急激なので、波動関数の変化も大きなも
のとなるため、それを記述するためにはより細かいメッシュが必要であり、一
方、原点からずっと遠くではポテンシャル、波動関数とも変化が緩やかなので、
広めのメッシュ間隔でも十分記述可能なためです。

実空間法でも、”adaptive coordinates”という上記に対応したような手法
があるようですが、詳細は筆者もよく知りません。以下に、”adaptive
coordinates”というキーワードで把握できた文献[1]を示します。実際は、実
空間法関連で(検索エンジンや、文献データベースなど)探せばこの手法に関
しての情報は割と得られるかと思います。

[1] D. R. Hamann, Phys. Rev. B63, 075107(2001)

http://www.bandstructure.jp/

●963 Re: 511について4 mike - 2004/08/06 18:33 -
小林一昭先生、ご指導、ありがとうございます。

ORBIT ELECTRONS ENERGYについて、
ご紹介いただきました渡辺先生のページの「多電子原子の全電子軌道」を拝見しました。
プログラムを含むため、まだ全部は見れていないのですが、軌道エネルギーの定義は
わかりました。 ありがとうございます。
また、実空間でスピンを入れたLDAが使われており、たいへん興味深いです。

>むしろ実空間法なら、周期的境界条件に縛られる必要がないので、
>この井戸を考える方法が使えるのではと思います。
>井戸というより電場勾配を与えるということでも良いと思います。
電場勾配を与えて電界放出、遠くに井戸型のポテンシャルを設置し、原子を高温にして
熱電子放出が再現できればと考えます。
そのために電子状態の時間発展が必要だと思われます 。

511の核間距離が実測と合わない問題に関し、ご指導いただきましたように、
F = - ∂E_tot/∂Rとの比較から、力の計算部分の検証をやってみようと思います。

●962 Re: 511について3(訂正) 小林一昭 - 2004/08/06 12:28 -
すみません、昨日の書き込みで井戸型のポテンシャルが
正電荷に相当する云々のところは筆者の考え違いでした。
従って、電荷中性を壊すという問題もありません。謹ん
で訂正します。

http://www.bandstructurer.jp/

●961 Re: 511について3 小林一昭 - 2004/08/05 22:45 -
>D. D. Koeling and B. N. Harmonの計算結果に関しまして、
>わからない部分があります。水素原子の各項の意味ですが、
>KINETIC ENERG = <1s|-Δ/2|1s>
>COR-ELE ENERGY = <1s|V(r)|1s>, V(r)=-1/r
>HARTREE ENERGY = <1s|VH|1s>
>EX COR ENERGY = <1s|Vxc|1s>
>TOTAL ENERGY = 上の各項の和(Heについては、上記の各軌道の和)
>と考えてよろしいのでしょうか?

それで良いと思います。

>このとき、ORBIT ELECTRONS ENERGYはどう定義されるのでしょうか?

以上に関しては、以前出てきた渡辺氏のページ
(Yahooの物理学/計算物理学から辿れる)の、
”計算物理のためのC/C++言語入門”、多電子
原子の全電子軌道のページでの詳しい説明が参
考になると思います。

ORBIT ELECTRONS ENERGYは軌道のエネルギーで、
筆者が前の投稿で提示したものは、水素なら軌道
に1個詰まった場合の軌道エネルギー、ヘリウム
なら2個詰まった場合の軌道エネルギーです。
これは、上で紹介したページでのイオン化エネル
ギーとは異なります。イオン化エネルギーの場合
は、水素の場合、軌道に0.5個、ヘリウムでは1.5個
での軌道エネルギーとなります。当該ページ下に、
このイオン化エネルギーと全エネルギーの表があり
ます。単位はa.u.と思われます。全エネルギーの値
は、筆者の提示したものとLDAの表式及び計算方法の
差の範囲内でよく一致しています。表にある各軌道
のエネルギーはイオン化エネルギーで、この値は筆者
の提示したものと異なります。筆者の提示したのは
イオン化させていない(0.5個または1.5個としない)
場合での軌道エネルギーです。一応、水素で、0.5個、
ヘリウムで1.5個にすると当該する表の値と良く整合
する値が得られることは確認しました。

>スラブ近似で電子状態を収束させた後、遠くに(スラブのまん中に)電子を捕獲
>する井戸を作りもう一度収束させ、その場所における電子密度を求めるようなこ
>とは可能でしょうか?

まず周期的なスラブモデルを前提とした場合、スラブ
とスラブの間の真空層を通常よりかなり大きな値をと
る必要があるでしょう。また、電子を捕獲する井戸と
いことは、正電荷に相当するものであり、これはその
ままでは系の電荷中性を壊すのではないかと考えます。
電荷中性でないと電場の勾配ができ、周期系において
は、これは深刻な問題となります。井戸が正電荷に
相当といことにちょっと自信がないのですが、電子を
引き寄せるものとしたら、そう考えるしかないように
思えます。

むしろ実空間法なら、周期的境界条件に縛られる必要
がないので、この井戸を考える方法が使えるのではと
思います。井戸というより電場勾配を与えるというこ
とでも良いと思いますが、具体的にこれ以上のイメージ
は、今の段階では思い付きませ。

>どの部分でしょうか? 電子の識別は、配列の番号で行っていますので、個別
>に扱っているようですが、非個別性は、交換項を入れ、各軌道の密度の和のみ
>をポテンシャルの計算に使うことで保証されると考えているのですが...。

すみません、筆者の計算ではスピンは今の場合up,down
とも完全に縮退していると考えいるので、2電子系で
電子を別々に扱っているのに違和感を感じたのですが、
これは筆者の勘違いと思われます。

http://www.bandstructure.jp/

●960 Re^3 511について mike - 2004/08/05 18:26 -
小林一昭先生、重ね重ねのご指導、詳しい情報、ありがとうございます。

D. D. Koeling and B. N. Harmonの計算結果に関しまして、
わからない部分があります。水素原子の各項の意味ですが、
KINETIC ENERG = <1s|-Δ/2|1s>
COR-ELE ENERGY = <1s|V(r)|1s>, V(r)=-1/r
HARTREE ENERGY = <1s|VH|1s>
EX COR ENERGY = <1s|Vxc|1s>
TOTAL ENERGY = 上の各項の和(Heについては、上記の各軌道の和)
と考えてよろしいのでしょうか?
このとき、ORBIT ELECTRONS ENERGYはどう定義されるのでしょうか?

>ただ、動径方向の外側は、60 a.u.とかなり
>遠いところが境界となります(そこから外は考えない)。
>分子の計算を実空間で解く場合の境界の扱いについての
>検証が必要かもしれません。
SD法による水素原子の固有関数を求めるアプレットを作っているところですが、
ご指摘のように、実空間グリッドでは、境界をかなり遠くにとらないと、
固有値に影響が出るようです。(週末頃にはアップする予定です。)

(電子の放出に関して:補足)に関しまして、
通常のバンド計算で仕事関数を求める方法を解説していただき、ありがとうございます。

>密度汎関数法+局所密度近似(ここでは通常のバンド計算を想定)では、
>電子を個別的に扱うことが困難で、表面から電子が一個(複数個でも)出ていく過程を
>再現するように計算を遂行するのは、そのままでは出来ないように思います。
スラブ近似で電子状態を収束させた後、遠くに(スラブのまん中に)電子を捕獲する井戸を作り
もう一度収束させ、その場所における電子密度を求めるようなことは可能でしょうか?

>時間依存DFT-LDA計算で、電子を個別に扱っているように見える部分があるのですが、
>これは問題ないのでしょうか?。
どの部分でしょうか? 電子の識別は、配列の番号で行っていますので、個別に扱っているようですが、
非個別性は、交換項を入れ、各軌道の密度の和のみをポテンシャルの計算に使うことで保証される
と考えているのですが...。

●959 Re:Re:511について(補足2後半) 小林一昭 - 2004/08/04 23:27 -
すみません、エラーが出たので前半、後半と分けて
投稿しました。

以下に水素とヘリウムの場合の計算結果(固有値、
全エネルギーとそれを各項に分解した値)を示します。

水素原子の場合

ORBIT ELECTRONS ENERGY (RYD.) :リュードベリ単位
1 1S  1.0000     -0.4678477

交換相関項 = Perdew, Zungerの表式

【 HARTREE UNIT 】:原子単位(a.u.)
TOTAL ENERGY = -0.446158598591
KINETIC ENERGY = 0.424958005186
COR-ELE ENERGY = -0.920949909963
HARTREE ENERGY = 0.282836281461
EX COR ENERGY = -0.233002975275

ヘリウム原子の場合

ORBIT ELECTRONS ENERGY (RYD.):リュードベリ単位
1 1S  2.0000     -1.1409423

交換相関項 = Perdew, Zungerの表式

【 HARTREE UNIT 】:原子単位(a.u.)
TOTAL ENERGY = -2.834784654841
KINETIC ENERGY = 2.766404415021
COR-ELE ENERGY = -6.623782412909
HARTREE ENERGY = 1.995396766358
EX COR ENERGY = -0.972803423312

TOTAL ENERGY :全エネルギー
KINETIC ENERGY :運動エネルギー
COR-ELE ENERGY :原子核-電子間の相互作用項
HARTREE ENERGY :ハートリー項
EX COR ENERGY :交換相関項

軌道エネルギーと全エネルギーとで単位が異なるので
ご注意下さい。水素に関しては比較検証する他の計算
値が見つからなかったのですが、ヘリウムには全エネ
ルギーの値が見つかりました。LDAの型は不明(古い
タイプ)ですが、-2.83 a.u.という値で、上記の
ヘリウムの全エネルギーの値と一致します。それから
類推して、水素の値も多分間違いないだろうと考えます。

Wignerの表式を使用した場合、水素の全エネルギーは、
-0.447280129774 a.u.、ヘリウムの場合は、
-2.818996132019 a.u.となります。

以上は、D. D. Koeling and B. N. Harmon, J. Phys.
C10, 3107(1977)による原子の計算によるものです。
この計算手法は原子のみに対応し、分子以上の系には
対応しません。極座標表示におけるディラック方程式
を基に、動径方向の外側と、内側(r = 0 a.u.、原点
近傍の扱いは何か特別な処理をしていたと思います)
両方から解いていき、適当なところで接続して固有値等
を求めるものですが、詳細は筆者もちょっと忘れてしま
いました。ただ、動径方向の外側は、60 a.u.とかなり
遠いところが境界となります(そこから外は考えない)。
分子の計算を実空間で解く場合の境界の扱いについての
検証が必要かもしれません。

http://www.bandstructure.jp/

●958 Re:Re:511について(補足2前半) 小林一昭 - 2004/08/04 23:25 -
まずは訂正があります。先日示した文献、F. Herman and S.Skillman,"Atomic Structure Calculations",
Prentice-Hall Inc.(1963)、の名前が、Hermanでなく
Harmanとなっていました。謹んで訂正します。

この文献は古典的な名著と言えるものですが、筆者も昔
みたことがあるだけで、現在閲覧できる状況にありません。

(水素原子の固有値と水素分子の結合距離)
>DFT-LDAによる1s軌道の固有値は、厳密な理論値と比べて
>正確には一致しない由、軌道エネルギーがズレていても、
>水素分子の核間距離が数%以内で一致するのでしょうか?
>エネルギーの渟留値近くで波動関数の形が変化しても、
>エネルギー変化は小さくなることから逆に考えますと、
>エネルギーがずれると、関数の形が異なってきて、
>重なり積分が変化し、結合距離も変わるように思われ
>るのですが...。
>(的外れであれば、御容赦ください)

最後の方で計算例を提示しますが、水素原子の計算の
固有値は、前回で示した値及び単位で正しいことを確
認しました。その意味で、厳密な理論値と比べて、か
なりD. D. Koeling and B. N. Harmonの計算による
固有値はずれています。この方法は原子のみしか計算
できないので、水素分子に関しては何とも言えません。
ただ、この方法で作成した原子のポテンシャルを基に
擬ポテンシャルを作成し、それを使用したバンド計算
で、以前示したスーパーセルを用いて水素分子の結合
距離を求めると、LDAの計算で、1.46 a.u.程度の値が
得られます(0.01 a.u.程度の変動あり)。GGA(密度
勾配を用いたLDAの拡張)を用いれば、1.41 a.u.程度
の値が得られます。水素分子に関しては、その意味で
かなり実験値に合った結果が得られます。

(電子の放出に関して:補足)
密度汎関数法+局所密度近似(ここでは通常のバンド
計算を想定)では、電子を個別的に扱うことが困難で、
表面から電子が一個(複数個でも)出ていく過程を再
現するように計算を遂行するのは、そのままでは出来
ないように思います。仕事関数は、このような電子を
無限遠に持っていくことは行なわず、フェルミ準位と
真空準位の差から求めます。通常のバンド計算は周期
系を扱うので、スラブ近似を使い、真空層を厚めにと
れば真空層中心付近は大体真空準位となっているとみ
なせる。それでそこそこ合う結果が得られます。

時間依存DFT-LDA計算で、電子を個別に扱っているよ
うに見える部分があるのですが、これは問題ないので
しょうか?。

http://www.bandstructure.jp/

●957 Re 511について(補足) mike - 2004/08/04 17:37 -
小林一昭先生、ご指導、文献のご紹介、ありがとうございます。

DFT-LDAによる1s軌道の固有値は、厳密な理論値と比べて正確には一致しない由、
軌道エネルギーがズレていても、水素分子の核間距離が数%以内で一致するのでしょうか?
エネルギーの渟留値近くで波動関数の形が変化しても、エネルギー変化は小さくなることから
逆に考えますと、エネルギーがずれると、関数の形が異なってきて、重なり積分が変化し、
結合距離も変わるように思われるのですが...。
(的外れであれば、御容赦ください)

再度、アプレットNo.511のプログラムのVHとVxcの部分を見直しましたが、
少なくとも、コーディング上の誤りは見つけられませんでした。
Vxcの部分をスレーターのXα法の
 Vxc(ρ)=-3α[3ρ/4π]^(1/3) ... (0.66<α<1.0)
にしたもので試してみたいと思います。

>>最終的には実時間、実空間で金属表面の電子放出のような
>>現象を扱いたいためです。
>大変興味深い、挑戦的なテーマだと思います。
こう言っていただくと、とても励みになります。また、アドバイスありがとうございます。
まだ道のりは遠そうですが、少しづつ進めていけたらと考えます。
ご指導のほど、よろしくお願いします。

●956 511について(補足) 小林一昭 - 2004/08/03 23:57 -
 ちょっと時間がないのでテストルームの方は見ていないのですが、
少し補足があります。水素原子の固有値に関してですが、密度汎関
数法+局所密度近似による1s軌道の固有値は、厳密な理論値と比べ
てそう正確には一致しません。D. D. KoelingとB. N. Harmonの
方法[1]で解いた場合の水素の1s軌道の固有値は、LDAにPerdew,
Zungerの表式を用いたもので、-0.46785 Ry(a.u.では半分の値
になる)、Wignerの表式で、-0.46941 Ry(これもa.u.では半分)
となります。[1]の原子を解く方法では、相対論効果が考慮されて
いますが、水素の場合その影響はほとんどないと言えます。非相対
論として原子を解くもので参考となる文献として、F. HarmanとS.
Skillmanによるもの[2]があります。少々古めでLDA部分の扱いが、
Xα法かWignerの表式によるものと思われますが、これは論文でなく
書籍で、各原子の計算結果の詳しい表があったと思います。

 筆者もこの場合、どの程度まで合うのが良いのかちょっと計り
かねる部分があります。上記に挙げた数値の単位も間違っていな
いかも含めて、これはもう少し調べてみます。

[1]D. D. Koeling and B. N. Harmon, J. Phys. C10, 3107(1977)
[2]F. Harman and S. Skillman, "Atomic Structure
Calculations",Prentice-Hall Inc.(1963)

>最終的には実時間、実空間で金属表面の電子放出のような
>現象を扱いたいためです。

 大変興味深い、挑戦的なテーマだと思います。LDA(GGAも
含めて)は表面の計算でも仕事関数などは比較的実験値に合う
結果を返します。ただ、以前にも言及したように、LDAでは表面
の鏡像ポテンシャルや、自己相互作用の相殺が不完全なので、
電子を遠方に取り出した場合に、取り出した電子の感じるポテン
シャルが正しくないなどの問題があります。電子放出となると
ここら辺を検討する必要があるかと思います。

http://www.bandstructure.jp/

●955 水素原子の計算による検証 mike - 2004/08/03 20:47 -
水素原子の計算による検証のため、Mikes Room/TEST ROOMに
(TR10) H atom Kohn Sham 3D 水素原子(時間依存-DFT-LDA)をアップしました。

時間依存-局所 密度汎関数法のアプレット(No.504)を用いて、電子数n=1、外部ポテンシャルVext=-1/rとし
H原子を構成します。
この結果、電子のエネルギーは-0.35程度と理論値-0.5と比べ30%程度ズレています。
viewのVHとVxcを比べると、相殺がうまくいっていないようです。
No.480は同じメッシュの水素原子ですが、こちらは-0.474auと理論値から5%程度のズレで済んでいます。
このズレは原点の1/rを有限値で置き換えたことによると考えられます。

VHはPoissonの方程式ΔVH=-4πρ(r)をSOR法で解いています。
どこにも調節パラメータはありませんので、コードを再度見直して見ます。

みなさま>
No.502-No.508のアプレットは当面、ご使用にならないよう、お願いします。

●954 Re 511について mike - 2004/08/03 17:34 -
小林一昭 先生、懇切丁寧なご指導、ありがとうございます。

ご指導いただきました、F = - ∂E_tot/∂Rとの比較から、力の計算部分の検証を
やってみようと思います。

(疑問点1)につきまして、
本アプレットの実空間での時間依存Kohn-Sham方程式は、渡辺-塚田先生のアルゴリズムを
用いてプログラムしています。最終的には実時間、実空間で金属表面の電子放出のような
現象を扱いたいためです。
ですが、ご指摘のように、実空間グリッドで密度汎関数法+局所密度近似による
変分計算を試た後にしようと思います。

(疑問点2)につきまして、
>メッシュなどの条件を変えた場合(507、508)の、エネルギー固有値の値が倍も違っています。
>これは計算 が収束していないことを示しており、収束性に関して検証する必要があります。
No.507は核間距離が1.4auで固定され、No.508は可変でき初期値は3.0auですので、
そのままの比較はできません。No.508の核間距離を1.4auにしてlossボタンを何回も押して、
定常状態に近づけると比較できます。(lossボタンは波束の位相をそろえて運動を止め、
エネルギーを奪います。 時間発展は断熱的ですので、定常解に近付けるにはこの操作が必要です)
この方法で比較すると、軌道エネルギーはそれぞれ、約-0.74(0.25au)と約-0.63(0.5au)
となります。この影響があるかどうかは、わかりません。

(水素原子の計算による検証) につきまして、
>やはり、分子ではなく厳密に解ける水素原子の固有値が 計算で扱う近似の範囲内で
>正しい結果となるか検証する のが良いと思います。
前回も同様のご指摘を受けましたが、試しておくべきでした。 早速やって見ます。
できましたら、TEST ROOMにアップします。

●953 511について 小林一昭 - 2004/08/02 23:44 -
以下、511に関してのコメントです。

まず、力が計算できるということは、全エネルギーE_tot
(F = - ∂E_tot/∂Rなので)も計算可能なはずで、
二原子間の距離に関して、全エネルギーの極小と力がゼロ
になるところが一致するでしょうか?。あと、原子間の距
離のずれによる全エネルギーの差ΔEとそのずれの距離Δrか
ら、F = - ΔE/Δrが計算でき、こと対応するヘルマン-ファ
インマン力として計算した力Fを比較してみると、その力を
求める計算部分の妥当性が検証できると思います。

で、得られた全エネルギーが極小になる原子間距離(=平衡
位置)が、水素分子の実験値の2倍となった場合、これは計
算の前提部分で問題がある可能性があります。通常のLDAの
計算で、水素分子の原子間距離は、伸び目に出ますがせいぜ
い数%のオーダーです。2倍も大きくは出ません。この通常
のLDAの計算は、バンド計算において大きなユニットセルを
想定して、その中に水素分子を置いて計算するというもので
す。セルが十分大きければ、となりのセルからの影響は無視
できます。おそらく、量子化学的な手法で水素分子の計算を
行なった場合でも、LDAを前提にしてれば、ほぼ同様の結果
が出るはずと思われます。

以下、疑問な点を挙げますが、筆者の誤解、理解不足(特に
具体的にどのように計算されているのか把握しきれていない
ので)による勘違いの可能性もあります。その際はご容赦下
さい。

(疑問点1)
そもそも時間依存する形式で解く必要について疑問がありま
す。時間依存する計算はハミルトニアンが時間依存する場合
を扱うためにあり、ポテンシャルが時間とともに変化するよ
うな場合、例えば時間によって変化する動的な外場中や、
非断熱的な反応過程を取り扱う場合などが挙げられます。
水素分子の原子間の平衡距離を求める計算は、時間依存する
問題として取り扱う必然性はなく、ごくまっとうな変分問題
として解いて行くことが可能です。←力を基に水素分子間の
距離を変える毎に、電子状態が常にボルン-オッペンハイマー
面にあるように電子状態を収束させながら平衡位置を求めれ
ばOK。

時間依存問題は、まだ発展途上で方法論もまだ確立したとは
言えない手法で、そもそも、時間依存問題が関わる場合は、
励起状態を扱わなければなりません。密度汎関数法では励起
状態が正しく計算できる保証がありません。経験的には良く
合う場合もありますが、どのような場合でも通用するとは思
えません。

水素分子の動力学という話となると時間依存の話になるので
すが、水素は原子が軽過ぎるため量子効果が無視できないの
で(古典的に扱えない)、かなり厄介な話になります。平衡
位置を求めることや、電子状態の基底状態を求める場合はより
まっとうな通常の手法(密度汎関数法+局所密度近似による
変分計算)を用いるのが第一歩かと思います。

(疑問点2)
メッシュなどの条件を変えた場合(507、508)の、
エネルギー固有値の値が倍も違っています。これは計算
が収束していないことを示しており、収束性に関して検証
する必要があります。この収束性の問題が、今回の水素
分子の平衡位置の問題に関わっている可能性があります。
ただ、これだけが原因であるかは筆者には何とも言えま
せん。

(水素原子の計算による検証)
やはり、分子ではなく厳密に解ける水素原子の固有値が
計算で扱う近似の範囲内で正しい結果となるか検証する
のが良いと思います。既にこの手法で計算、比較検討さ
れていたらご容赦下さい。

http://www.bandstructure.jp/

●952 バグフィックス mike - 2004/08/02 19:40 -
バグフィックス(規格化メソッドの誤りを訂正)しました。
(510) steepest descent method 最急降下法による固有値と固有関数 ver 0.0.2

●951 Re 吸着のアプレット mike - 2004/08/02 18:43 -
Noriさん、こんにちは。

>では、河村先生の研究報告をじっくりとよんでいきたいと思います。
私も量子力学の関連が一段落したら、水をやりたいです。

>吸着のアプレットに関する質問なのですが、vAjustmentメソッドはどこからも呼び出されていないような気がするのですが、
>これは必要なのでしょうか?
No.326のアプレットのことでしょうか?
たしかにどこからも呼ばれていないようですね。作る時に似たようなアプレットを土台として
作りますので、そのときの消し忘れですね。ご指摘ありがとうございます。

>また温度を求めるメソッドとして、tempCont()、temperature()がありますが、これは両方必要なのでしょうか?
>片方だけではだめなのでしょうか?
tempCont()はタングステン原子(MAT0)のみの温度(基板の温度)をコントロールするために
設けました。temperature()は系全体の温度、temp(int k)はある種類だけの温度を返します。

>またzsort()の32767.0という数字はどういう意味があるのでしょうか?
srtz[1][i] = (int)(32767.0*pz[j]/zMax);
で0-32767の整数がsrtz[1][i]にセットされます。ソートできる程度の数であれば
他の数でもいいと思います。
zSort()では奥行の方向にソートし、後ろから表示して前後の重なりを表現しますので、
あまり小さな数でなければ良いと思います。

●950 質問 Nori - 2004/08/02 16:47 -
返答ありがとうございます。
では、河村先生の研究報告をじっくりとよんでいきたいと思います。

吸着のアプレットに関する質問なのですが、vAjustmentメソッドはどこからも呼び出されていないような気がするのですが、
これは必要なのでしょうか?

また温度を求めるメソッドとして、tempCont()、temperature()がありますが、これは両方必要なのでしょうか?
片方だけではだめなのでしょうか?

またzsort()の32767.0という数字はどういう意味があるのでしょうか?
多くの質問申し訳ありません

●949 アプレット追加 mike - 2004/08/01 17:37 -
アプレット追加しました。
(511) H2 molecule QMD 水素分子の動力学(暫定版) ver 0.0.0

このアプレットは、水素原子間の平近距離が実測の2倍程度になってしまいます。
ただ、プログラム上の誤りというより、考え方に原因があると推量したため、
敢て公開させていただきました。
みなさま、ご意見、ご指摘をいただければ幸いです。よろしくお願いします。

●948 アプレット追加 mike - 2004/07/31 15:34 -
アプレット追加しました。
(510) steepest descent method 最急降下法による固有値と固有関数

●947 表計算で拡散方程式 mike - 2004/07/31 05:25 -
わかさん、こんにちは。

拡散方程式を表計算ソフトで計算する場合、横方向に空間的な刻み(たとえば1mmごとのガスの濃度)
縦方向に時間刻みdtを考え、あるセルの濃度c(i,next)は上のセルc(i)とその両隣りの濃度c(i-1)、c(i+1)と
拡散係数D、dtからc(i,next)=c(i)+D*dt{c(i-1)+c(i+1)-2c(i)}/dx^2としたらいいかと思います。
ただし、D*dt/dx^2が0.5未満になるようにdtを選ぶ必要があります。
境界条件は空気と接している側は一定濃度(たとえばc(1)=1.0)、
反対側は濃度勾配なし(たとえばc(N-1)=c(N))とするといいかも知れません。

●946 ありがとうございました わか - 2004/07/30 20:28 -
みけさん、ありがとうございました。返事が遅くなり申し訳ありません。7月17日から、出張やグループ会社会議などのため時間がとれませんでした。おっしゃることわかりました。あるいはわかったつもりです。表計算ソフトで計算に挑戦します。またどうぞよろしくおねがいします。

●945 Re 原子パラメータ mike - 2004/07/30 18:04 -
Noriさん、こんにちは。

>カオリナイトを構成する要素としてO、Si、OHを追加するために、使用するかどうかはわかりませんが、これらのMorseポテンシャルにおけるパラメータを探しています。
>質量はすぐにわかったのですが、D、A、r0がなかなか見つけることができません。
>初心者の〜にあった参考文献も探したのですが、見つけれませんでした。よろしければ参考になさったサイト、資料等ありましたら教えていただけないでしょうか。
>またレナードジョーンズポテンシャルではCOなども要素として使用していましたが、MorseポテンシャルでもOHなどは使用可能なのでしょうか?

Al-O、Si-O、O-Hは共有結合とイオン結合の性格が強く、局在し、方向性があると考えられます。
Morseポテンシャルは、方向性がなく、「初心者の〜」の表にも掲載されているように金属結合は比較的
得意なのですが、上記のような結合を表現するポテンシャルとしては、適しているとは言えません。
おそらく、このような事情のために見つからないのではないでしょうか。

No.936でご紹介いただきました河村先生のポテンシャルの形が適していると思います。
この中の(6)式のようにイオン結合的な部分はクーロンポテンシャルと反発するイオン殻、
共有結合的な部分は相手を特定したMorseポテンシャルで表現できそうです。
また、資料などがありましたら、再度、書き込みします。
>また、パラメータのDは結合エネルギー、r0は原子間距離として、Aは何を意味するパラメータなのでしょうか?

私もハッキリとしたことは知りませんが、Aが小さいほど結合が締まっていて、大きいほど緩やかになります。
Aが小さいほど、結合のバネが硬くなります。遷移金属では小さく、アルカリ金属では大きくなります。

●944 原子パラメータ Nori - 2004/07/30 16:49 -
カオリナイトを構成する要素としてO、Si、OHを追加するために、使用するかどうかはわかりませんが、これらのMorseポテンシャルにおけるパラメータを探しています。
質量はすぐにわかったのですが、D、A、r0がなかなか見つけることができません。
初心者の〜にあった参考文献も探したのですが、見つけれませんでした。よろしければ参考になさったサイト、資料等ありましたら教えていただけないでしょうか。
またレナードジョーンズポテンシャルではCOなども要素として使用していましたが、MorseポテンシャルでもOHなどは使用可能なのでしょうか?
また、パラメータのDは結合エネルギー、r0は原子間距離として、Aは何を意味するパラメータなのでしょうか?

●943 イオン結合 Nori - 2004/07/30 10:54 -
ご指導ありがとうございます。
ぜひ参考にさせていただきます。
現在吸着面のタングステンをカオリナイトの結晶構造に変更するところまでできました。
といってもカオリナイト内の分子間の相互作用は無視して、その計算はカットして静止させ、初期配置を構造通りに並べただけですが…
イオン結合を考慮できるようになるにはもう少し時間がかかりそうです。

●942 Re 分子数 mike - 2004/07/27 18:20 -
Noriさん、こんにちは。

>今、カオリナイトの結晶構造を作成しているのですが、 粒子数を増やすと止まってしまいます。
>これは配列の数が原因だと思い、reg[][]や、rijReg[][]の数を250から1000に変えたところ遅いながら動きました。
>この初期宣言の配列の数の250や他の860、10200といった数には何か意味があるのでしょうか?
reg[][]やrijReg[][]はカットオフと粒子登録法で使う配列です。
あるi番目の粒子に対し、FSポテンシャルのカットオフ距離より少し大きな球を考え、
球内の粒子番号をreg[i][]に登録します。 これによって、高速化しています。
ijReg[][]は一度計算した距離を覚えておいて再度使うためのものです。
250などの数はカットオフ距離内の粒子の最大数より少し大きめの数で、
小さすぎると、ご指摘のように止まってしまいます。
カットオフ距離との兼ね合いで決めました。

int regMorse[][] = new int[NN][860];
はMorseポテンシャルに対する上記と同様のことです。
double forceTable[][][] = new double[3][3][10200];
は力の計算をする代わりに表を引くことによって高速化しています。
forceMorse(double r, int knd0, int knd1)で
forceTable[kindT[knd0]][kindT[knd1]][(int)(r*1.0e13)]としています。

>分子数を増やしても止まらないようにするには他に変更する必要がある箇所などはあるでしょうか?
いま、思い付くかぎりでは、だいじょうぶと思います。

なお、カオリナイトは酸化物系ですので、「No.313 イオン結合のMDテンプレート」を
参考にされるといいかも知れません。酸化物にモース型のポテンシャルを入れる方が
簡単かと思います。大局的に電荷が残らないような構造では、これでいけると思います。
大局的に電荷が残る時は、クーロンポテンシャルについてエバルト法などの方法で
計算する必要があります。

●941 分子数 Nori - 2004/07/27 12:56 -
返答ありがとうございます。
今、カオリナイトの結晶構造を作成しているのですが、
粒子数を増やすと止まってしまいます。
これは配列の数が原因だと思い、reg[][]や、rijReg[][]の数を250から1000に変えたところ遅いながら動きました。
この初期宣言の配列の数の250や他の860、10200といった数には
何か意味があるのでしょうか?
分子数を増やしても止まらないようにするには他に変更する必要がある箇所などはあるでしょうか?

●940 アプレット追加 mike - 2004/07/24 16:33 -
アプレット追加しました。
(509) variational method 変分法による基底状態の近似エネルギーと波動関数

●939 Re 表面吸着 mike - 2004/07/23 18:08 -
Noriさん、こんにちは。

>この研究ではスメクタイトを用いていますが、その代わりにカオリナイトを使えないか、と検討中です。
河村先生のスメクタイトと同じように酸化物系+OHと元素も近いので、有効だと考えます。

>この(No326の)配置の原点はどこになるのでしょうか?
>x軸は左端、y軸は最上部、z軸は一番奥でいいのでしょうか?
>格子でいうと左上の奥ということになりますが・・・
ご指摘の通りです。変則的ですが、左手系(親指x軸、人指し指y軸、中指z軸(前後))ですね。
左下奥の角を原点とする右手系(親指x軸、人指し指y軸、中指z軸(上下))がふつうでしょうか。
鏡像などの問題があると、注意が必要かも知れません。

●938 表面吸着 Nori - 2004/07/23 15:19 -
コメントありがとうございます。
この文献も大いに参考にして研究に役立てたいと思います。
この研究ではスメクタイトを用いていますが、その代わりにカオリナイトを使えないか、と検討中です。

ところでアプレットに関する質問なのですが、No326の吸着のアプレットについて、タングステンを格子状に配列していますが、
この配置の原点はどこになるのでしょうか?
x軸は左端、y軸は最上部、z軸は一番奥でいいのでしょうか?
格子でいうと左上の奥ということになりますが・・・

●937 Re O-H間のポテンシャルについて mike - 2004/07/22 18:09 -
Noriさん、こんにちは。 興味深い情報、ありがとうございます。よく探されましたね。

ご紹介いただきました、河村先生の研究報告を拝見しました。相互作用ポテンシャルとして、2体力(4式)と3体力(5式)が考慮され、「原子間ポテンシャル関数は,すべての原子間の相互作用を記述しており,実際の計算においては分子内,分子間の区別はない。」と書かれていますので、O-H間のポテンシャルとして使えるのではないでしょうか。
また、Oは結合相手によってパラメータを変えているようで、No.933でコメントした「Oに結合する相手によって、O-Hの電子状態が変わる」効果も取り入れておられるようです。
参照している系もバルクというより、層状の構造を持っていますので、表面吸着にも適用できそうですね。
また、水でなければ、3体力は省略できそうです。
「SiOCH系動力学用ポテンシャルパラメータ」に関しては、まだ見れていません。
とりあえず、感想を述べさせていただきました。

●936 O-H間のポテンシャルについて Nori - 2004/07/22 15:58 -
O-H間のポテンシャルについて何か参考になるようなものはないかと調べていましたら以下のようなサイトを見つけました。
利用可能かどうかはわかりませんが、ちょっと調べてみようと思います。
もしお時間のほうありましたら、ご感想などいただけたら幸いです。
SiOCH系動力学用
 ポテンシャルパラメータ
http://www-fact.tokyo.jst.go.jp/content/h10/material/M02/parameter.doc
分子動力学法によるスメクタイトに対するSrの収着構造の研究
http://www.jnc.go.jp/siryou/gihou/pdf2/5983.pdf

後者はパソコン分子シミュレーションの著者である、河村雄行教授の研究報告で、河村教授の研究室のHPにいけばソースをダウンロードできるようです。

●935 アプレット追加 mike - 2004/07/21 17:02 -
アプレット追加しました。
(508) H2 coarse TDKS3D 水素分子(粗い格子DFT-LDA版)

●934 アプレット追加 mike - 2004/07/20 18:23 -
アプレット追加しました。
(507) H2 molecule TDKS3D 水素分子(時間依存-DFT-LDA)

●933 Re O-H間のポテンシャル mike - 2004/07/20 18:21 -
Noriさん、こんにちは。

>ソースを少し変更して2種類の分子を扱えるようになりました。
>おそらくあってると思いますが…
よかったですね。また、癖のあるコードを見ていただき、ありがとうございます。

>さて、話は変わるのですが、シミュレーションによって、O-Hの結合を再現したいのですが、適当なポテンシャル等はご存知でしょうか?
残念ながら、今まで、見たことは、ありません。
O-Hの結合は、私も興味があります。Si-O-Hのように表面の疎水性、親水性を決める重要な因子として
多くの分野で話題になっていますね。ただ、この場合、A-O-HのようにO-HのOに結合する相手によって、
O-Hの電子状態が変わると考えられ、まじめに考えるなら、分子軌道法や第一原理計算が必要と思います。
私も、第一原理計算をやりたくて、現在、量子力学の復習をしているところです。
表面の吸着が扱えるところまでは、まだしばらくかかると思います。

●932 O-H間のポテンシャル Nori - 2004/07/20 17:15 -
返答ありがとうございます。
ソースを少し変更して2種類の分子を扱えるようになりました。
おそらくあってると思いますが…

さて、話は変わるのですが、シミュレーションによって、O-Hの結合を再現したいのですが、適当なポテンシャル等はご存知でしょうか?
水分子のポテンシャルはいくつかあるようですが、それが適応できるかどうか…
水分子ではなくて、水素−酸素間の結合をシミュレーションによって再現したい、というよりはOH基を再現したのですが…
O-H間は共有結合だと思うのですが、そういった結合を再現できるポテンシャルをご存知でしたら教えていただけないでしょうか?

●931 アプレット追加 mike - 2004/07/19 17:52 -
アプレット追加しました。
(506) Gram Schmidt TDKS1D 波動関数の直交化

●930 バージョンアップとアプレット追加 mike - 2004/07/18 17:06 -
バージョンアップ(コントロールパネルの改良)しました。
(33) Brownian motion ブラウン運動 ver 0.0.3

アプレット追加しました。
(505) He atom Kohn Sham 3D ヘリウム原子(時間依存-DFT-LDA)

●929 アプレット追加 mike - 2004/07/17 16:35 -
アプレット追加しました。
(503) time dependent Kohn Sham 2D 二次元の時間依存-密度汎関数法
(504) time dependent Kohn Sham 3D 三次元の時間依存-密度汎関数法

●928 気体の相互拡散 mike - 2004/07/17 05:50 -
わかさん、こんにちは。 「分子のおもちゃ箱」をご覧いただき、ありがとうございます。

さっそくですが、
>このたびの質問は、79番のアルゴン中の各種気体の拡散に関してです。
>あれをもう少し定量的といったらいいのでしょうか、以下のようなことがわかるようにできないでしょうか。
>内径3mm、長さ100mmの閉じた石英ガラス管にアルゴンが1気圧詰まっていて、
>1方を空気中で開けたらどのくらいの時間で、空気中の水分子が他方に到達するか、
>それとどのぐらいの時間で平衡して混ざるか?
おたずねのような(〜mm程度の)スケールは、現在の分子動力学法によって直接シミュレートできる範囲
(<数100nm)より大きいため、難しいと思います。 むしろ、マクロな拡散問題として取り扱うのが良いと
考えます。相互拡散係数を分子動力学法によって計算することはできます。

上記の問題は管内径と長さの比が大きく、しかもマクロな流れが存在しないようなので、
ほぼ理想的な1次元の拡散問題と考えられます。(たとえは、アプレットNo.225を参照)
このとき、必要な情報はアルゴン中の酸素、もしくは窒素の相互拡散係数です。
(あるいは、アルゴン中の水蒸気も考慮する必要があるかもしれません)
アルゴン、窒素、酸素のいずれも質量が近いので、第0近似ではアルゴンの自己拡散係数で代用できそうです。
室温、1気圧のときの拡散係数は、D〜0.1cm2/s(10^-5m2/s)の桁になり、
拡散時間τはL=2sqrt(Dτ), L=0.1m, として
τ=L^2/4D〜0.1^2/(4x10^-5)=250s
約4分ほどで拡散し、ガスの本体が端に到達することになります。

しかし、このようなガスが少しでも到達すると良くない場合は、アプレットNo.225を見てもわかりますが、
これよりはるかに短い時間で、微量ながらすることを留意する必要があります。
この拡散の計算は、表計算ソフトでも簡単に行えますので、挑戦されたらいかがでしょうか。

●927 教えてください わか - 2004/07/16 18:57 -
みけさんこんにちは。ずっと以前に排気と吸着について質問させて頂いて、教えていただきありがとうございました。その後、勉強して自分でシミュレーションを作ろうなどと考えていましたが、あんまり難しいので全く進んでおりません。せっかく教えていただいたのに、以後進めないのは大変失礼と思いますが、良い製品を作ろうという意思に免じてど
うか今後もお世話してください。このたびの質問は、79番のアルゴン中の各種気体の拡散に関してです。あれをもう少し定量的といったらいいのでしょうか、以下のようなことがわかるようにできないでしょうか。内径3mm、長さ100mmの閉じた石英ガラス管にアルゴンが1気圧詰まっていて、1方を空気中で開けたらどのくらいの時間で、空気中の水分子が他方に到達するか、それとどのぐらいの時間で平衡して混ざるか?みけさんには以前教えていただいて、超短時間であるとじぶんでは納得しているのですが、他の人に説明しなければなりません。製品の排気プロセスに関わります。せっかく排気して加熱して内部きれいにしてもいったんブレークするとなんにもならないことを数字で示したいのです。製品そのものでは、そのことを証明しましたが、どうもピンときていないようなので。ご意見お願いします。

●926 バージョンアップ mike - 2004/07/15 18:32 -
バージョンアップ(直交性の表示を追加)しました。
(502) time dependent Kohn Sham 1D 一次元の時間依存-密度汎関数法 ver 0.0.3
時間発展によっても、波動関数の直交性はかなり良い精度で保持されます。

●925 バグフィックス mike - 2004/07/14 18:18 -
バグフィックス(timeEvolution()にバグがあり訂正)しました。
(502) time dependent Kohn Sham 1D 一次元の時間依存-密度汎関数法 ver 0.0.2

>小林先生、わざわざ、ありがとうございます。

●924 リンク確認しました(感謝) 小林一昭 - 2004/07/13 23:10 -
リンクを確認しました。深く感謝です。
http://www.bandstructure.jp/

●923 リンク追加 mike - 2004/07/13 18:02 -
リンク追加しました。
「Band Structure」
小林一昭 先生が運営されているバンド計算の総合情報サイトです。
バンド計算入門、擬ポテンシャルデータCD-R配布、バンド屋さんマップなど

●922 linkの件 mike - 2004/07/12 19:29 -
小林一昭先生、ご指導と励まし、ありがとうございます。
ほかにも誤りや不適切な表現などありましたら、ご指摘ください。

また、リンクのご承諾、ありがとうございます。
現在、HPへのアップロードができない状態になっていますので、
復旧後、リンクをはらせていただきます。

●921 Re Re^4:一次元DFT-LDAの妥当性テスト(補足) 小林一昭 - 2004/07/11 23:25 -
 先ほどの書き込みはちょっと舌足らずだったので、
補足します。

 局所密度近似はかなり大胆な近似で、問題点も多
数ある近似です。それを勘案すれば、一次元DFT-LDA
と492の一致は筆者が思っていた以上に良いです。
 先に挙げたように固有値の問題はちょっとひっかか
りますが。これは今すぐ対処すべき必要もないでしょ
う。

 是非、今後も頑張って下さい。
http://www.bandstructure.jp/

●920 Re Re^4:一次元DFT-LDAの妥当性テスト 小林一昭 - 2004/07/11 23:06 -
>電子密度を評価し、Perdew, Zungerの方法でVxcを求める方法
>で当面やりたいと考えます。

 それで良いと思います。具体的に困難な問題が出てきたら、
その時点で対処すれば良いと思います。そもそも、局所密度
近似自体、以下のページ、

http://www.geocities.co.jp/Technopolis/4765/INTRO/dftldaj.html

の欠点で示したように、問題点が多数あります。

>先生のページへ「分子のおもちゃ箱」からlinkさせて
>いただいてもよろしいでしょうか?

 深く感謝です。リンクは自由に張っていただいて構いません。
http://www.bandstructure.jp/

●919 Re Re^3:一次元DFT-LDAの妥当性テスト mike - 2004/07/11 07:09 -
小林一昭先生、ご指導、ありがとうございます。

一次元DFT-LDAのアプレットNo.498に関し、dxに対し仮想的にdx*lylzの体積を考え、
電子密度を評価し、Perdew, Zungerの方法でVxcを求める方法で当面やりたいと考えます。
よい近似方法があれば取り入れていきたいと考えますので、ご提案いただければ幸いです。

先生の「バンド計算入門」のページを拝見させていただいているのですが、むずかしくて、
なかなか先に進めません。少しづつ読ませていただこうと思います。
先生のページへ「分子のおもちゃ箱」からlinkさせていただいてもよろしいでしょうか?

●918 Re:ReRe:一次元DFT-LDAの妥当性テスト 小林一昭 - 2004/07/10 23:41 -
>変分原理は「どんな試行関数も真の固有関数の<ψ,Hψ>=E
>より小さくない」ことを保証しますが、どんなHの近似つい
>ても、Eより小さくないと言えるのでしょうか?

 これについては、492が真に厳密な解としたら、その固有
値より低い値を与える(492の近似としての)ハミルトニア
ンは、そもそも棄却されるべきでしょう。←つまりハミルト
ニアンおける近似が正しくないか、導入の仕方が妥当でない
と考えられます。

>このあたりのことは、私には糸口さえ見えません。
>なにか考えるヒントがありましたら、おねがいし
>ます。

 これを含めて後半の話に、今の筆者にはご助言できる
ほどのアイデア、情報がないのが実情です。
http://www.bandstructure.jp/

●917 アプレット追加 mike - 2004/07/10 18:21 -
アプレット追加しました。
(502) time dependent Kohn Sham 1D 一次元の時間依存-密度汎関数法

●916 アプレット改名 mike - 2004/07/09 20:19 -
アプレット名を改めました。(小林先生のご指摘により)
(501) Fermi sea DFT LDA 1D 自由電子の海の形成

●915 Re 吸着種 mike - 2004/07/09 20:14 -
Noriさん、こんばんは。

>しかし、2種類にすると現実とはかけ離れたものになるのでしょうか?
>お互いの相互作用もしっかりと働くようプログラムが組まれていると思うのですが…
>現実とは違うとは具体的にはどのようなことなのかよろしければ教えていただきたいです。
一般的に、表面は、バルクの結晶構造が変更を受け、再構成されることが知られています。
一つの原子を取り囲む配位数が表面で変わり、安定状態がバルクと一致しないためと考えられます。
原子間のポテンシャルはバルクのときの構造を再現するように決められているため、これを
表面に適用するとうまくいかないことが多いと思います。
一方、吸着エネルギーは構造にそれほど敏感でなく、Morseポテンシャルでも
うまくいく場合が多いと思います。

なお、WにKとCaを吸着させると高温では必ずW-Ca-Kのように層状に吸着しますが、
これはMorseポテンシャルでも再現します。

●914 吸着種 Nori - 2004/07/09 19:36 -
>吸着種は、比較的簡単な変更で複数にかえることができると思い
>ます。ただ、重ねた時の構造は、Morseポテンシャルですので、
>現実とかけ離れたものになりそうですが...。

吸着種を2種類にするのはたしかに簡単そうに思えます。
実際ソースもMAT2まで用意してありますし…
しかし、2種類にすると現実とはかけ離れたものになるのでしょうか?
お互いの相互作用もしっかりと働くようプログラムが組まれていると思うのですが…
現実とは違うとは具体的にはどのようなことなのかよろしければ教えていただきたいです。

●913 ReRe:一次元DFT-LDAの妥当性テスト mike - 2004/07/09 18:51 -
小林一昭先生、ご指導、ありがとうございます。

>一つ気になるのは、得られた固有値E1の値で、lylz=16の場合、密度汎関数法の方
>のE1の値の方が低くなっています。一次元DFT-LDAの方は変分で解いていることにな
>るので、この場合、E1の値は厳密解よりエネルギーが低くなるのはちょっと問題が
>あるのかもしれません。
「厳密解より低い固有エネルギーが得られるのは、変分原理からみて問題」というご指摘に
関しまして、質問があります。
変分原理は「どんな試行関数も真の固有関数の<ψ,Hψ>=Eより小さくない」ことを保証しますが、
どんなHの近似ついても、Eより小さくないと言えるのでしょうか?

>lylzをより小さい方に変えると、E1の値は、厳密解としての E1_exactより値が高くなり、
>lylzを大きくすると、より低い値になるようです。lylz により擬似的に3次元として与える
>電荷密度 が変分パラメーターとして良くないのか、 そもそもやはり一次元の系に、Perdew,
>Zungerの表式を用いたことに問題があるのか は、今の段階では筆者には何とも言えません。
このあたりのことは、私には糸口さえ見えません。なにか考えるヒントがありましたら、
おねがいします。

lylzがいわゆる調節パラメータとなり、原理計算にとってあまり好ましくないかも知れません。
むしろ3次元の箱を考えて、陽にLy、Lzを考え、変数分離で一次元の問題とするアプローチが
良いのでしょうか...。

●912 Re:一次元DFT-LDAの妥当性テスト 小林一昭 - 2004/07/08 22:44 -
mike's room/TEST ROOM/TR09拝見しました。

>lylz=4x4=16のとき、ほぼVHとVxcは相殺し、
>概ね妥当な結果を与えそうです。

 確かに良く一致していると思います。

 一つ気になるのは、得られた固有値E1の
値で、lylz=16の場合、密度汎関数法の方
のE1の値の方が低くなっています。一次元
DFT-LDAの方は変分で解いていることにな
るので、この場合、E1の値は厳密解より
エネルギーが低くなるのはちょっと問題が
あるのかもしれません。lylzをより小さい
方に変えると、E1の値は、厳密解としての
E1_exactより値が高くなり、lylzを大きく
すると、より低い値になるようです。lylz
により擬似的に3次元として与える電荷密度
が変分パラメーターとして良くないのか、
そもそもやはり一次元の系に、Perdew,
Zungerの表式を用いたことに問題があるのか
は、今の段階では筆者には何とも言えません。
http://www.bandstructure.jp/

●911 一次元DFT-LDAの妥当性テスト mike - 2004/07/07 19:53 -
小林一昭先生のアドバイスをもとに、一次元の密度汎関数法(局所密度近似)のアプレット(No.498)の電子密度ρを3次元的に考え、dx*lxlyとすることの妥当性を評価しました。
mike's room/TEST ROOM/TR09をご覧ください。
上記アプレットの一電子の場合のVeffとEiを一次元の井戸型ポテンシャルのアプレット(No.492)と
比較して見ると、lylz=4x4=16のとき、ほぼVHとVxcは相殺し、概ね妥当な結果を与えそうです。

●910 Re 原子の種類を増やすには mike - 2004/07/07 18:43 -
Noriさん、こんにちは。

>今、吸着面としてタングステンを用いていますが、上から例えばタングステン、鉛、カルシウムなど、違った原子を層状に重ねることはできるでしょうか?
本アプレットはタングステンのW(110)面上にCaなど一種類の原子を吸着させるものです。
基体金属はWに限定されます。これを変えるのは大幅な変更が必要です。吸着種は、比較的簡単な変更で複数にかえることができると思います。ただ、重ねた時の構造は、Morseポテンシャルですので、現実とかけ離れたものになりそうですが...。

>今のところ初期配置だけを変更して、相互作用等は無視してもかまわないのですが…
> ソースを見ると
>....
>となっているので、分子数をNNx*NNy*NNz*4 としてforループを追加していけば上手くいくかと思った>のですが、なかなか上手くできません。
>下の吸着面を層状に違う原子を重ねてみたいのですが、よろしければアドバイスをいただけないでしょうか?
引用されている部分はタングステンのbcc結晶を組み上げる部分です。
(少しずらして、W結晶の立方の部分と体心の部分を配置しています)
このプログラムはver 0.0.1ではないでしょうか?
現在はver 0.0.2で少し変更していますので、こちらをご覧ください。
タングステンのbcc結晶を組み上げた後に
原子間距離を考えて間隔を開けながら組み上げるといいです。
(本アプレットでは遠くから吸着させていますが、直接組み立てないと層構造は作り難いと思います)

>あと原子数をxyz=12,6,12でやりますと、原子数をNNx*NNy*NNz*4
>にすると計算がとまってしまいます。
>これは原子数を増やしたために、どこかの配列数をオーバーしてしまっているのでしょうか? 詳細はcodeを見ないとわかりませんが、ご指摘の配列数のオーバーと複数の粒子に同じ位置を指定した時、止まってしまうことを経験しています。

●909 原子の種類を増やすには Nori - 2004/07/07 17:00 -
早速質問なんですが、N0326のアプレットの吸着シミュレーションについてお聞きしたいことがあります。
今、吸着面としてタングステンを用いていますが、上から例えばタングステン、鉛、カルシウムなど、違った原子を層状に重ねることはできるでしょうか?
今のところ初期配置だけを変更して、相互作用等は無視してもかまわないのですが…
 ソースを見ると
for (i=0; i kind[i] = MAT0;
ix = i/(NNy*NNz);
iy = (i%(NNy*NNz))/NNz;
iz = (i%(NNy*NNz))%NNz;
xx[i] = s*(double)(ix)+0.5*sgm;
yy[i] = s*(double)(iy)+6.5*sgm;
zz[i] = s*(double)(iz)+0.5*sgm;
vx[i] = 100.0*Math.random()-50.0;
vy[i] = 100.0*Math.random()-50.0;
vz[i] = 100.0*Math.random()-50.0;
ffx[i] = 0.0;
ffy[i] = 0.0;
ffz[i] = 0.0;
fbx[i] = 0.0;
fby[i] = 0.0;
fbz[i] = 0.0;
}

for (i=NNx*NNy*NNz; i kind[i] = MAT0;
j = i-NNx*NNy*NNz;
ix = j/(NNy*NNz);
iy = (j%(NNy*NNz))/NNz;
iz = (j%(NNy*NNz))%NNz;
xx[i] = s*(double)(ix)+0.5*s+0.5*sgm;
yy[i] = s*(double)(iy)+0.5*s+6.5*sgm;
zz[i] = s*(double)(iz)+0.5*s+0.5*sgm;
vx[i] = 100.0*Math.random()-50.0+500.0;
vy[i] = 100.0*Math.random()-50.0;
vz[i] = 100.0*Math.random()-50.0;
ffx[i] = 0.0;
ffy[i] = 0.0;
ffz[i] = 0.0;
fbx[i] = 0.0;
fby[i] = 0.0;
fbz[i] = 0.0;
}
となっているので、分子数をNNx*NNy*NNz*4 としてforループを追加していけば上手くいくかと思ったのですが、なかなか上手くできません。

下の吸着面を層状に違う原子を重ねてみたいのですが、よろしければアドバイスをいただけないでしょうか?

あと原子数をxyz=12,6,12でやりますと、原子数をNNx*NNy*NNz*4
にすると計算がとまってしまいます。
これは原子数を増やしたために、どこかの配列数をオーバーしてしまっているのでしょうか?

●908 Re906 2周年 mike - 2004/07/06 17:42 -
Noriさん、お久しぶりです。また、書き込みありがとうございます。

>ついに500個を超えましたね!おめでとうございます!
>これからもますますの更新を楽しみにしております。
ありがとうございます。 新作のペースは落ちていますが、
あきらめず、続けていこうと思っています。

>またアプレットについて、わからない点など教えていただけたら幸いです。
わからないことをいっしょに考えていけたらと思います。
遠慮なく、なんでも書き込んでください。
今後とも、よろしくお願いします。

●907 Re 一次元密度汎関数法コメント3 mike - 2004/07/06 17:40 -
小林一昭 先生、コメント、ご指導ありがとうございます。
No.498の1電子での評価法、試みてみます。

>あと(501)ですが、これは孤立した系での計算でしょうか?
>(境界条件はどうなっているか?)。孤立した系だとすると、
>(501)で示された状況では電子状態がバンドを形成するとは
>思えないのですが、、、。←電子は依然として局在したまま
>ではないか?。
このアプレットの意図は、孤立系で、原子が集まってきて、ある準位が分裂し、
やがてエネルギーの拡がった帯状の準位を作ることをビジュアルに見ることでした。
自由電子モデルの出発点として考えたもので、「Fermi sea」とした方が
良かったかも知れません。(まだ、周期的な境界条件は付けていません)
ご指摘のように「band」は不適切だったと思います。
ただ、この場合でも、(電子数を2個にしてみるとわかりやすいのですが)
電子は局在しているわけではなく、波動関数は全体に拡がっていているので、
自由電子モデルへの橋渡しとしてはいいかなと考えます。

>そもそも自分に馴染みのあるバンドというのは、波数kとエネ
>ルギー固有値(準位)との分散曲線を考えるのですが、(501)
>で言っておられるバンドの意味が良く分からないのですが、、、。
この先、自由電子モデルからバンド理論へと見てわかるアプレットを作っていけたらと
考えています。 ご指導のほど、よろしくお願いします。

●906 2周年おめでとうございます。 Nori - 2004/07/06 16:21 -
ごぶさたしております。
半年以上前になりますが、No326の吸着のアプレットについていろいろとご質問させていただいておりましたNoriです。
この半年間、分子動力学法をさらに学ぶとともに、市販されているソフトなどを用いてなんとか土壌への分子吸着のシミュレーションができないかと模索していましたが、なかなか思うようにできず、もう一度mikeさんのアプレットを参考にと思い、書き込みをさせていただきました。
 アプレットの方はずっと拝見させていただいておりました。
ついに500個を超えましたね!おめでとうございます!
これからもますますの更新を楽しみにしております。
またアプレットについて、わからない点など教えていただけたら幸いです。

●905 一次元密度汎関数法コメント3 小林一昭 - 2004/07/05 23:48 -
>Kohn-Sham方程式の数値解(ポテンシャルVeff、軌道関数φi)
>を比較するということでしょうか?

おっしゃる通りです。

>自己の作る電荷密度によるポテンシャルを取り込んでしまう
>ように思えるのですが

局所密度近似では、電子数に依らずご指摘の自分自身との相互作用
(自己相互作用)が出来てきてしまいます。平均場近似でもハート
リー・フォック近似では、クーロン相互作用の項と交換項がうまく
相殺するので、この自己相互作用部分は消えますが、局所密度近似
では相関項が出てくるため自己相互作用が完全に相殺されずに残っ
てしまいます。これは電子の数に関係しないで出てきます。これが
局所密度近似の欠点の一つです。

あと(501)ですが、これは孤立した系での計算でしょうか?
(境界条件はどうなっているか?)。孤立した系だとすると、
(501)で示された状況では電子状態がバンドを形成するとは
思えないのですが、、、。←電子は依然として局在したまま
ではないか?。

 そもそも自分に馴染みのあるバンドというのは、波数kとエネ
ルギー固有値(準位)との分散曲線を考えるのですが、(501)
で言っておられるバンドの意味が良く分からないのですが、、、。
http://www.bandstructure.jp/

●904 リンク追加 mike - 2004/07/05 19:03 -
リンク追加しました。
Webラーニングプラザ:
科学技術振興機構(JST)が提供する、技術者向けラーニングサービスです。

●903 Re 一次元密度汎関数法コメント2 mike - 2004/07/05 18:29 -
小林一昭 先生、コメント、そしてアドバイスありがとうございます。

>(498)の設定では電子数が1個の場合がないのですが、電子数1個の場合は、条件次第で厳密な解が得られると思われる(?)ので、
>それと密度汎関数法+局所密度近似の結果と比較すると妥当性について判断できるかもしれません。
電子数が1個の場合、物理的には外場のみのSchrodingerの方程式の解になると考えられます。これと形式的にn=1とおいたKohn-Sham方程式の数値解(ポテンシャルVeff、軌道関数φi)を比較するということでしょうか?
(電子数が1個の場合、平均場近似では、自己の作る電荷密度によるポテンシャルを取り込んでしまうように思えるのですが)

●902 一次元密度汎関数法(コメント2) 小林一昭 - 2004/07/04 21:58 -
(498)についてですが、仮想的にly.lzを与えて、3次元としての電荷密度と考えるのは悪くない手だと思います。そもそも局所密度近似はかなり粗い近似なので、次元性による影響を議論する以前に局所密度近似を用いることの妥当性の方が問題になるのかもしれません。

(498)の設定では電子数が1個の場合がないのですが、電子数1個の場合は、条件次第で厳密な解が得られると思われる(?)ので、それと密度汎関数法+局所密度近似の結果と比較すると妥当性について判断できるかもしれません。
http://www.bandstructure.jp/

●901 アプレット追加 mike - 2004/07/04 18:44 -
アプレット追加しました。
(501) band DFT LDA 1D バンドの形成


 ( prev | index | next )
(created 2004.07.04, last updated 2006.07.13)
inserted by FC2 system