2013年02月09日

結晶、収束検討、(SCFサイクルでの収束条件、toldfe)

SCFサイクルの収束条件で検討してみます。
収束条件を厳しく(エネルギー差を小さくすれば)、
確かな値になるが、時間がかかるのではないかとの予想です。
以下の部分をいじって計算しました。

toldfe 1.0d-6 #Tolerance on the difference of total energy

以下に結果を示します。
横軸は指数部分でプロットしてあります。

2013020905.png

2013020906.png

全エネルギーは殆ど変化なしです。格子定数は変わります。
取りあえず10d-6以下で計算していれば大丈夫のようです。
posted by イトー at 21:02| セミナーでのチュートリアル | このブログの読者になる | 更新情報をチェックする

結晶、収束検討(k点メッシュ、ngkpt)

k点のメッシュをどれほど取ればいいか検討しました。
基本的には以下の場所を変えます。

ngkpt 2 2 2 #Number of grid points

例えばこれを

ngkpt 4 4 4 #Number of grid points

などに変えます。
その時の全エネルギー、格子定数の変化を以下に示します。

2013020901.png

2013020902.png

2013020903.png

2013020904.png

これらの結果より今まで

ngkpt 2 2 2

で計算していましたが、最低でも

ngkpt 4 4 4

の条件で計算した方がいいことが分かります。
なので、今後の検討は

ngkpt 4 4 4

で行います。
posted by イトー at 20:15| セミナーでのチュートリアル | このブログの読者になる | 更新情報をチェックする

2013年02月04日

結晶、収束検討(カットオフエネルギー、ecut)

では、まずカットオフエネルギーについて検討します。
カットオフエネルギーはどれだけの数の平面波を使うかに関係します。
可能な限りカットオフエネルギーを高くし、用いる平面波の数を大きく
した方が適した計算になるのは想像できると思います。

それではカットオフエネルギーを変化させ、格子定数と全エネルギーを見ていきます。

カットオフエネルギーを検討していると以下のようなエラーメッセージが出てきました。
Abinitのエラーメッセージは結構充実していて読めば分かると思います。
この場合も「初期の格子定数からの変化が大きいからdilatmxを大きくしなさい」との
ことです。

2013020401.png

では、検討に戻ります。abinitin.txtの以下の部分を変えて、abinitexe.batをダブルクリックします。

変更前
#Definition of the planewave basis set
ecut 20.0 # Maximal plane-wave kinetic energy cut-off, in Hartree

変更後
#Definition of the planewave basis set
ecut 2.0 # Maximal plane-wave kinetic energy cut-off, in Hartree

計算語、abinitout.txtを開き格子定数(acell)と全エネルギー(etotal)を確認します。
777行目に
acell 1.0253378556E+01 1.0253378556E+01 1.0253378556E+01 Bohr
782行目に
etotal -3.4833366715E+01
となっていると思います。
値をグラフソフトにコピー&ペーストします。
次の計算をするために、abinitin.txt、abinitexe.bat、files.txt以外は消去するか、
フォルダーを作って移動させてください。
(計算が混乱しないようにするためです。)
ファイルが整理できたら次のカットオフエネルギーで計算してみましょう。

カットオフエネルギーと全エネルギーの関係

Graph2.png

この関係を見るとこの結晶ではカットオフエネルギーは20以上は欲しいところです。
20以上で詳しく見ると

Graph3.png

やはりカットオフエネルギーを大きくした方が全エネルギーが下がります。
ただ、20と40のエネルギー差は0.0006 Hatree (16meV)です。
この差が大きいか否かは、議論したい内容で決まると思います。

カットオフエネルギーと格子定数の関係

Graph0.png

格子定数に関してもカットオフエネルギーは20以上欲しいです。
20以上を詳しく見ると

Graph1.png

20〜40で0.0005 au (0.00025 Å)程度のバラツキとなります。
これもどこまで見たいか、計算時間との兼ね合いだと思います。

このように様々なパラメーターに関して収束検討を行う必要があります。
結晶構造が変わったらもちろん行わなければなりません。
posted by イトー at 21:08| セミナーでのチュートリアル | このブログの読者になる | 更新情報をチェックする

収束検討の重要性

これから可能な限り多くのパラメーターに関して、
収束検討を行おうと思います。

第一原理計算に関しては理想的に言えば、
厳しい条件でやればやるほどいいと思います。
(「いい」が真の値とは違います。そのソフトを用いた場合での「いい計算条件」です。)
例えば、収束条件を厳しくしたり、カットオフエネルギーを高くしたり、
k点メッシュを数多く取ったり。
理想的に言えば厳しいほどいいと思います。
しかし、計算機から言わせるとたまったものではありません。
計算機にとっては楽な計算ほどいいのです。
なので精度が高く、楽な計算を見つける必要があります。
これ以上厳しい計算条件にしても格子定数や全エネルギーが
変わらないことを確認しなければなりません。
多分、第一原理計算をものにするにはこの作業がとても重要になると思います。
この作業を丁寧にやることによって、Abinitがどのような計算でどこまでできるか?
どこまで議論していいのか?分かると思います。
posted by イトー at 21:02| セミナーでのチュートリアル | このブログの読者になる | 更新情報をチェックする

結晶(結果)

結果ですが、以下のファイルに出力されています。

abinitout.txt

いつも使っているテキストエディタで開いてください。
今のところ850行目以下のみ見て頂ければいいと思います。

-outvars: echo values of variables after computation --------
acell 1.0214212321E+01 1.0214212321E+01 1.0214212321E+01 Bohr
amu 2.80855000E+01
dilatmx 1.05000000E+00
ecut 2.00000000E+01 Hartree
ecutsm 5.00000000E-01 Hartree
etotal -3.5482292910E+01

acellが格子定数で、etotalが全エネルギーになります。
格子定数は

10.214212321 au (5.405126233 Å)となります。

全エネルギーは

-35.482292910 Hartree (-965.5299617 eV)

です。

次は計算条件によって格子定数、全エネルギーがどのように変化をするか検討します。
私は第一原理計算で次から行うような収束検討はとても大切だと思います。
収束検討を行うことで、第一原理計算の精度などを身につけられると考えています。
posted by イトー at 20:41| セミナーでのチュートリアル | このブログの読者になる | 更新情報をチェックする

2013年01月23日

結晶(計算)

では、Siの計算の準備ができましたので、実行してみましょう。

abinitのフォルダーの中は以下のようになっています。

2013012301.png

pspのフォルダーの中には14si.pspncの擬ポテンシャルのファイル入っています。
abinitexe.batをダブルクリックすると計算が始まります。
以下のような画面になり。計算が終わるとその下のような画面になります。
何かキーを押すと戻れます。

2013012302.png

2013012303.png

計算後はabinitのフォルダーは以下のようになっています。

2013012304.png

abinitフォルダー内のabinitout.txtを開きます。
最後から3行目が以下のようになっていればOKです。

Calculation completed.
.Delivered 0 WARNINGs and 3 COMMENTs to log file.
+Overall time at end (sec) : cpu= 12.3 wall= 13.0
posted by イトー at 23:20| セミナーでのチュートリアル | このブログの読者になる | 更新情報をチェックする

2013年01月16日

結晶(入力ファイルの確認)

再度表示しますが、入力ファイル(abinitin.txt)は以下のような内容になっています。

#Si

#Definition of the optimization of cell parameter
acell 3*10.0 #Cell parameter
chkprim 0 #Checking primitive cell
optcell 2 #Optimazing cell parameter
ionmov 3 #Control the displacements of ions
ntime 10 #Number of time step
dilatmx 1.05 #Maximal permitted scaling of the lattice parameters
ecutsm 0.5 #Energy cutoff smearing

#Definition of the atoms
ntypat 1 #Number of Types of atoms
znucl 14 #Charge -Z- of the nucleus
natom 8 #Number of atoms
typat 1 1 1 1 1 1 1 1 #Type of atoms
xred #vectors of atom positions
0.125000 0.125000 0.125000
0.875000 0.875000 0.875000
0.625000 0.125000 0.625000
0.375000 0.875000 0.375000
0.125000 0.625000 0.625000
0.875000 0.375000 0.375000
0.625000 0.625000 0.125000
0.375000 0.375000 0.875000

#Definition of the planewave basis set
ecut 20.0 # Maximal plane-wave kinetic energy cut-off, in Hartree

#Definition of the k-point grid
kptopt 1 #k-points option
ngkpt 2 2 2 #Number of grid points

#Definition of the SCF procedure
nstep 10 #Number of SCF
toldfe 1.0d-6 #Tolerance on the difference of total energy

posted by イトー at 23:00| セミナーでのチュートリアル | このブログの読者になる | 更新情報をチェックする

結晶(SCF計算)

次にSCF計算の条件について説明します。

SCF計算については以下のHPが役に立ちます。
(中山先生すみません。使わせてもらいます。)
http://homepage3.nifty.com/mnakayama/research/tips/qm.htm

入力ファイルの以下の部分にSCF計算の条件が記載されています。

#Definition of the SCF procedure
nstep 10 #Number of SCF
toldfe 1.0d-6 #Tolerance on the difference of total energy


nstep
SCF計算の最大サイクル数。この数で計算が止まる場合は、設定を大きくしてください。

詳細説明
http://www.abinit.org/documentation/helpfiles/for-v7.0/input_variables/varbas.html#nstep

toldfe
前の計算と今回の計算のエネルギー差がこの設定値以下になればSCF計算は終了します。大きくすれば計算は直ぐ終わりますが、初心者の方は「1.0d-6」にしておきましょう。単位はHartreeになっています。eVに直すと

 1.0d-6(0.000001) Hartree = 0.0000272 eV = 0.272 meV

となります。

詳細説明
http://www.abinit.org/documentation/helpfiles/for-v7.0/input_variables/varbas.html#toldfe
posted by イトー at 22:52| セミナーでのチュートリアル | このブログの読者になる | 更新情報をチェックする

結晶(k点)

次にk点に関して説明します。
kは波数ベクトルで、これによってどのような平面波か決まります。理想的には全てのkに対して求めなければなりませんが、計算が終わりません。なので適当なkで計算することになります。
結晶は3次元で、逆格子空間であるkも3次元となります。その3次元でどのようにkを取って計算するかが重要になります。つまり、3次元のメッシュでk点を取って計算することになります。

友人のHPを参考にしてください。
http://homepage3.nifty.com/mnakayama/research/tips/qm.htm

入力ファイルでは以下の部分になります。

#Definition of the k-point grid
kptopt 1 #k-points option
ngkpt 2 2 2 #Number of grid points


kptopt
k点の取り方を定義します。Defaultで1になっています。通常1で大丈夫です。多分。詳細は英語を確認してください。確か、逆格子空間の原点から少しずらしたりできるはず。

詳細説明
http://www.abinit.org/documentation/helpfiles/for-v7.0/input_variables/varbas.html#kptopt

ngkpt
逆格子空間の3次元でどれだけのメッシュの交点のkで計算するかを定義します。多ければ多いほど理想に近づきますが、その分、3乗で計算時間が長くなります。k点を大きくしていき、格子定数、全エネルギーが安定した条件を使う方がいいと思います。
また、a,b,c軸で格子定数が違う場合にも注意が必要です。k点は逆格子空間ですから、実空間で格子定数が長いは逆格子空間では短くなります。なので、格子定数の逆数に比例するようにメッシュを取りましょう。例えばa,b軸が5Å、c軸が10Åだった場合は

ngkpt 4 4 2

のように取ります。

詳細説明
http://www.abinit.org/documentation/helpfiles/for-v7.0/input_variables/varbas.html#ngkpt
posted by イトー at 22:12| セミナーでのチュートリアル | このブログの読者になる | 更新情報をチェックする

結晶(カットオフエネルギー)

以下のカットオフエネルギーに関して説明します。

#Definition of the planewave basis set
ecut 20.0 # Maximal plane-wave kinetic energy cut-off, in Hartree


カットオフエネルギーはどれだけの平面波を使って計算するかの入力パラメータです。
カットオフエネルギーが大きければ大きいほど多くの平面波を計算に利用するので、
計算精度が良くなります。しかし、その分、計算量は多くなり、時間がかかります。
カットオフエネルギーによって結果が変わりますので、比較したい計算は必ず、
カットオフエネルギーを同じ値にしましょう。単位はHartreeです。
また、対称となる結晶構造の収束検討を行ってください。
同じ元素、同じ構造の大きさで行うことが大切です。
カットオフエネルギーを大きくしていき、格子定数、全エネルギーの値が変化しなくなったらOKです。

詳細説明
http://www.abinit.org/documentation/helpfiles/for-v7.0/input_variables/varbas.html#ecut
posted by イトー at 21:54| セミナーでのチュートリアル | このブログの読者になる | 更新情報をチェックする
×

この広告は1年以上新しい記事の投稿がないブログに表示されております。