Your first calculation: CO molecule in a box

input file

 0  0  0  0  0  0                      : dummy line (6 integers)
 5.50 20.00  2  2  2                   : GMAX, GMAXP, NTYP, NATM, NATM2
 1  0                                  : space group number, bravis lattice type
 6.00  4.00  4.00  90.00  90.00  90.00 : a, b, c, alpha, beta, gamma
 1  1  1  1  1  1                      : knx, kny, knz, k-point shift
 1  0                                  : NCORD, NINV
 0.0000  0.0000  0.0000  1  1  1       : cps, iwei, imdtyp, ityp
 2.2000  0.0000  0.0000  1  1  2       : cps, iwei, imdtyp, ityp
 6  0.1500  51577.50 3 1 0.d0          : IATOMN, ALFA, AMION, ILOC, IVAN, ZETA1
 8  0.1500  51577.50 3 1 0.d0          : IATOMN, ALFA, AMION, ILOC, IVAN, ZETA1
 0  0  0  0  0                         : ICOND, INIPOS, INIVEL, ININOSE, INIACC
 0  1                                  : IPRE, IPRI
 200  200   0    57200.00  0           : NMD1, NMD2, iter_last, CPUMAX, ifstop
 3   1                                 : way_mix, mix_what
 0    8  0.8                           : starting mixing, kbxmix,alpha
 0.60  0.50  0.60  0.70  1.00          : DTIM1, DTIM2, DTIM3, DTIM4, dtim_last
 30.00    2     1  0.10D-08 1.d-06     : DTIO, IMDALG, IEXPL, EDELTA
  0.0010  0.10D+02    0                : WIDTH, FORCCR, ISTRESS
ggapbe          1                      : XCTYPE, nspin
  1.00     3                           : destm, n_stm
102                                    : NBZTYP
   0   0   0                           : NKX,  NKY,  NKZ  (dummy)
   0   0   0                           : NKX2, NKY2, NKZ2 (dummy)
   8                                   : NEG (# of bands)
       1                               : NEXTST (1: G-space, 0: R-space)
       0                               : 0; random numbers, 1; matrix diagon
       2                               : imsd (2: Davidson, 1: RMM)
       0                               : eval. eko diff.: .0 = no ,1 = yes
       0                               : npdosao
       0    0.0                        : SM_dopping

ダミー行

 0  0  0  0  0  0                      : dummy line (6 integers)

歴史的経緯でそのままになっている。過去のバージョンでは1カラム目の変数が1の場合は※基本並進ベクトルをカーテシアン座標で与え、0の場合はa, b, c, alpha, beta, gammaで与えていたようである。
※バージョン539tでは活きているようです。

カットオフと原子数

 5.50 20.00  2  2  2                   : GMAX, GMAXP, NTYP, NATM, NATM2
  • GMAX: 波動関数を展開する平面波のカットオフエネルギーを指定。E_{cut}^{wf} = GMAX**2
  • GMAXP: 補強電荷を展開する平面波のカットオフエネルギー。E_{cut}^{dens} = GMAXP**2

GMAX=5.00 GMAXP=15.0 の条件で計算を行う場合は、一般的にpbe3 pbe1sのpseudo potentialのセットを用いるとよい
GMAX=5.50-6.0 GMAXP=20.0 の条件で計算を行う場合はpbe1のpseudo potentialセットを用いるとよい
※これは一般的なルールなのでエネルギーと力のカットオフエネルギーについての収束性は必ずチェックすべきである。

  • NTYP: 原子の種類
  • NATM: 空間反転対称性(NINV=1)がある場合の非等価な原子の数
  • NATM2:全原子数、空間反転対称性が無い場合はNATM=NATM2である。

格子の対称性と種類

 1  0                                  : space group number, bravais lattice type
  • NUM_SPACE_GROUP: 空間群の番号
  • TYPE: ブラベー格子の種類を番号で指定:
    #bravais lattice type
    0simple
    1body-centered
    2face-centered
    3a-face-centered
    4b-face-centered
    5c-face-centered
    6rhombohedral

 6.00  4.00  4.00  90.00  90.00  90.00 : a, b, c, alpha, beta, gamma
  • a, b, c: 格子ベクトルの長さ
  • alpha, beta, gamma: 格子ベクトルのなす角度

これらの代わりに以下のように"Cartesian"と記述するとそれ以下の行で格子ベクトルを直接カーテシアン座標で指定することができる。

Cartesian
 6.00  0.00  0.00
 0.00  4.00  0.00
 0.00  0.00  4.00

 1  1  1  1  1  1                      : knx, kny, knz, k-point shift

Brillouine zone積分で用いるk点のメッシュ(knx, kny, knz)を最初の3カラムで、k点のオフセットを次の3カラムで指定する。 4-6カラムでのk点のオフセットの指定方法は、ガンマ点を含む場合には1を、ガンマ点を含まない場合(Monkhorst-Pack)は2を指定する。 六方晶の場合にはk点はオフセットしない方が良いようである。

原子位置の指定方法と空間反転対称性

 1  0                                  : NCORD, NINV
  • NCORD: 1: カーテシアン座標で原子座標を入力(単位はBohr)、0: 基本並進ベクトルを単位として入力、2: conventional unit cellの格子ベクトルを単位とする
  • NINV: 0: 系に空間反転対称性がない、1: 空間反転対称性がある。
    このオプションを使用すると、空間反転対称性が存在する場合には明示的に考慮する原子数をおよそ半分にすることが可能で、計算時間の短縮をすることが可能である。

原子位置と原子種

 0.0000  0.0000  0.0000  1  1  1       : cps, iwei, imdtyp, ityp
 2.2000  0.0000  0.0000  1  1  2       : cps, iwei, imdtyp, ityp
  • 1-3カラム: cps(katm,3) (pos(katm,3)): 原子座標。入力方法はNCORDに依存。
  • 4カラム: iwei(katm): 空間反転対称性について等価な原子数
  • 5カラム: imdtyp(katm): Hellmann-Feynman力に従って原子を動かす場合には1、そうではない場合には0を指定。有限温度分子動力学(FTMD)の場合、1000以上の数を指定するとその原子は熱浴に接することを指定。
  • 6カラム: ityp(katm): 原子の種類を番号で指定。原子番号等は後に続く列に記述する。

原子(擬イオン)に関する情報

 6  0.1500  51577.50 3 1 0.d0          : IATOMN, ALFA, AMION, ILOC, IVAN, ZETA1
 8  0.1500  51577.50 3 1 0.d0          : IATOMN, ALFA, AMION, ILOC, IVAN, ZETA1
  • 1カラム: 原子番号
  • 3カラム: 原子の質量(Da(旧のa.m.u.)(=水素原子が約1))。後述するが、上の値は現実のものとは異なる。Quenched MD(IMDALG=2)を行うときには原子質量は同じ値、例えばCの12.0107やSiの28.0855などにしておくと計算を安定して行うことができる。 有限温度の分子動力学を実行するときには現実の質量に設定するべきである。
  • 6カラム: zeta=(rho(up)-rho(down))/(rho(up)+rho(down))の初期値を与える。スピン分極の計算を行う際はnspin=2としてzetaとしては0でない有限の値を設定する。zeta=0とすると非磁性の解に収束するだろう。
  • 6カラム: (光田くんによる説明)spin分極の初期値(一般的に分極してそうなものは0.5、そうでないものは0.0を入力する)
  • その他: ilocはローカルポテンシャルの軌道量子数l+1、ivanはVanderbiltのultrasoft pseudopotentialを使用するかどうか、といったパラメータであるが、入力ファイルからではなく擬ポテンシャルのファイルから読むので実際にはダミーパラメータである。ちなみにAlFAという変数は現在は使用されていない。電荷密度を乱数から発生させるときに使用していたようである。

計算の継続(初期化)に関するオプション

 0  0  0  0  0                         : ICOND, INIPOS, INIVEL, ININOSE, INIACC
ICOND
0初期波動関数は乱数から発生させる。scratchからSCF計算を開始する時に使用する。
1初期波動関数と電荷密度は前の計算のものを使用する。SCF計算を継続する際に使用する。zaj.dataとpotential.dataが必要
2初期波動関数は乱数から、電荷密度は前の計算のものを使用するが、電荷密度(ポテンシャル)はfixしたまま計算を行う。potential.dataが必要
3初期波動関数とポテンシャルは古いものから出発するが、電荷密度(ポテンシャル)はfixしたまま計算を行う. ポテンシャルをfixした計算(ICOND=2)を継続する際に用いる。potential.dataとzaj.dataが必要。
4電荷密度をfixしたままSCF計算を行う。ICOND=2と同じ。
9, 11実空間での電荷密度を出力. (9は全電荷密度、11はsoft partのみ)詳細
10STM
12DOS_CAL 電子状態を読み込んで修正計算せずにDOSを計算する
14, 24Partial density of state (PDOS)の計算(24はK点分解で出力)
15 (115)WFN_FFT 実空間での波動関数の出力(ソースファイルの変更が必要. ソースの変更が不要なバージョンは濱田から入手可能)(115はK点分解出力用と思われるが詳細不明)
17, 117Crystal Orbital Overlap Population (COOP Analysis)(117はK点分解で出力,16は実験的旧バージョン?)
18CHG_DIFF_ANA 電荷密度差の表示(全系と、それを2つに分割した系の和との差、3系準備して入力)
20INITIAL_MG 詳細不明
22, 23BAND STRUCTURE バンド計算用(23は継続計算)
33, 133 (13)ALDOS XY面に平行な面でSuperCellをスライスした断片各々のDOS(133はK点分解で出力、13は旧バージョン(使用しない))
40規約ブリルアンゾーン内で波動関数(zag.data)と電荷密度(potential.data)をGWST用に書き換える(STATE 5.3.8b用)
4141: nfkpt.dataで指定した経路上で波動関数と電荷密度をGWST用に書き換える(STATE 5.3.8b用)


INIPOS原子位置
0原子位置を入力ファイルから読み込む
1restart.dataから原子座標を読み込んでMDまたは構造最適化を継続する
2GEOMETRYから原子座標を読み込んで計算を継続する(restart.dataが必要)。restart.dataはバイナリーのデータで扱いにくいが、GEOMETRYはascii形式であるので、比較的扱いやすい。計算を現在のgeometryと異なるものから始めるときに使用すると便利である。
  • INIVEL: 原子の速度
    • INIPOSと同様、1の場合はrestart.dataから、2の場合はGEOMETRYから速度を読み込む。INIPOSの時と同様にrestart.dataが必要。
  • ININOS: 能勢の温度制御
  • INIACC: accumulator:物理量の統計平均の計算の継続に関するオプション

計算の継続する例 (1): SCF計算の継続

  • SCF計算のみを継続したいとき(geometryは入力ファイルから読み込む)
     1  0  0  0  0                         : ICOND,INIPOS,INIVEL,ININOS, INIACC
  • 構造最適化を継続
     1  1  0  0  0                         : ICOND,INIPOS,INIVEL,ININOS, INIACC
  • MDの継続(IMDALG <=0)
     1  1  1  0  0                         : ICOND,INIPOS,INIVEL,ININOS, INIACC

 0  1                                  : IPRE, IPRI
  • 圧力(ストレス)の計算(現在は実装されていないのでダミー変数である)
  • 標準出力のレベル(通常は1、1よりも大きいと多くの情報が出力される。デバッグ用のオプションではあるが、出力される情報はあまり参考にならないようである)

繰り返しの回数と計算時間に関するオプション

 200  200   0    57200.00  0           : NMD1, NMD2, iter_last, CPUMAX, ifstop
  • NMD1: 電子系のSCF計算の回数の上限を指定する。
  • NMD2: 分子動力学のステップ数の上限を指定する。NMD2+1が実際のMD(構造最適化)の回数の上限である。
  • iter_last:
  • CPUMAX: CPU時間の上限を秒単位で指定
  • ifstop:

ミキシング

 3   1                                 : way_mix, mix_what
way_mixミキシングの方法
1simple mixing
2Broyden
3Broyden2
4DFP
5Pulay
6Blugel


mix_whatミキシングの対象
1電荷密度
2ポテンシャル

ミキシングに関するパラメータ

 0    8  0.8                           : starting mixing, kbxmix,alpha
  • start mixing: dummy
  • kbmix: ミキシングに用いる履歴の回数
  • alpha=mix_alpha: ミキシングパラメータ

ミキシングに関するパラメータ(?)

 0.60  0.50  0.60  0.70  1.00          : DTIM1, DTIM2, DTIM3, DTIM4, dtim_last 
  • DTIM1: RMMで使用、その他はダミー変数.
  • dtim_last: ダミー変数

構造最適化または分子動力学(MD)に関するオプション

 30.00    2     1  0.10D-08 1.d-06     : DTIO,IMDALG,IEXPL,EDELTA
  • DTIO: イオン系のMDで用いる時間刻み(タイムステップ)(a.u.)。Quenched MD (IMDALG=2)またはMDで有効。 1a.u.=0.024188fs, 41.3428a.u.=1fs
  • IMDALG: MDまたは構造最適化の方法を指定。正の値は構造最適化、0または負の場合は有限温度MDを実行
  • Molecular dynamicsも参考にしてください。
IMDALG構造最適化の方法
1Newtonian dynamics
2Quenched molecular dynamics
3Vibrational mode analysis(nfvibrate.dataというファイルが必要)Vibratinal 詳細はこちら
4GDIIS
5TS search by GDIIS
6NEB Nudged Elastic Band method
7CINEB Climbing Image NEB method
0Newtonian dynamics
-1有限温度のNewtonian dynamics(MVELSC=0:Velocity scaling =10:Nose-Hoover)
-2Langevin MD
  • IEXPL: dummy
  • EDELTA: 電子系の収束の閾値を規定。収束は一原子あたりの全エネルギーで評価している。

  0.0010  0.10D+02    0                : WIDTH, FORCCR, ISTRESS
  • WIDTH: Fermi面のぼかしの幅. 正の場合はデルタ関数を多項式で、負の場合にはHermite-Gaussianの一次式で近似する。WIDTH<-10.0のときにはTetrahedron法によるBrillouine zone積分を行う。この時はNBZTYPが有効になる。
  • FORCCR: Forceの閾値。Forceがこの値よりも小さくなると構造最適化を終了する。IMDALG>0のみ有効である。
  • ISTRESS: ストレスの計算用オプション。現在は使用されていない。

ggapbe          1                      : XCTYPE, nspin
  • XCTYPE: 交換相関汎関数を指定.
    • ggapbe: Perdew, Burke, Ernzerhof GGA (1996)
    • ggapw91: Perdew-Wang GGA (1991)
    • ldapw91: Perdew-Wang L(S)DA (1991)
    • その他RPBEとrevPBEが使用可能. 最新版にはWu-Cohen GGAも実装(Experimental)
  • NSPIN: 1ならば非磁性、2でcollinearな磁性体の計算を行う。

STM用オプション

  1.00     3                           : destm, n_stm

102                                    : NBZTYP
  • WIDTH < -10.0の時のみ有効。101が推奨値。

ダミー行

   0   0   0                           : NKX,  NKY,  NKZ  (dummy)
   0   0   0                           : NKX2, NKY2, NKZ2 (dummy)

計算に用いるバンドの数

   8                                   : NEG (# of bands)

STATEの計算では、常に非占有状態を含めて計算を行う必要がある。非磁性のSiならば、例えば4 (電子数/2)+4などとするのが良いだろう。

nonlocal pseudopotentialの計算方法

       1                               : NEXTST (1: G-space, 0: R-space)
  • 0: R-space (Davidsonの時は使用しない)
  • 1: G-space

ダミー行

       0                               : 0; random numbers, 1; matrix diagon

対角化法の選択

       2                               : imsd (2: Davidson, 1: RMM)
  • 1: RMM-DIIS
  • 2: Davidson
  • 3: ???

大規模な計算をするときにはRMM-DIISとReal space projectionの組み合わせ(NEXTST=0、IMSD=1)が推奨されますが、初期波動関数は乱数から発生させるものではなく、Davidsonである程度収束させたものを使うのが良いでしょう。そうしなければ最安定な電子状態に収束しない可能性がある。

       0                               : eval. eko diff.: .0 = no ,1 = yes

前回のiterationにおける擬固有値との比較を行うかどうかを指定。ほぼダミー変数

PDOSの計算オプション

       0                                                   : npdosao

0でない場合にPDOSの計算を行う。但しnpdosaoは原子の数より小さくなければならない。具体的な計算方法は、柳澤氏による以下の文書を参照のこと。

f電子用の半経験的パラメータ用のオプション

       0    0.0                                          : SM_dopping

現在使用している人はほとんどいないため十分テストされていません。

PDOS/AOLDOSの設定

(ここのDOSの記述はAOLDOS。PDOSはたしかLayer-resolved DOSだったような気がする。Layer-resolved DOSの計算は井関氏の卒論研究データに詳しい。)

原子軌道に射影したPDOS(Partial Density Of States、もしくはAO_LDOS)の計算方法

(以下、pdos用のオプション。仮にnpdosao = 5とし、通し番号が1,3,4,9,10の原子の軌道に投影したDOSを計算、原子の種類は2種類。ただし、lda+uでない場合)

   5                                   : npdosao
   1                                   (以下、npdosao個の原子の通し番号)
   3
   4
   9
  10
-15.00 5.00 0.20 501                   :EPDOS(1), EPDOS(2), EPDOS(3), NPDOSE
2.2   0.3                              :rad, width for ityp 1
1.0   0.3                              :rad, width for ityp 2
0.2   12                               :DR, NR
  • npdosao: DOSを投影させる原子数。その次の行にnpdosao個の原子の通し番号を続ける
  • EPDOS(1,2): pdosを表示するエネルギー領域(Fermi準位との相対値、eV)、EPDOS(3): 未使用
  • NPDOSE: EPDOS(1)~EPDOS(2)の間のエネルギーメッシュ数
  • rad: DOSを投影させる原子中心の波動関数(擬波動関数)の、動径方向のカットオフ値を指定
    (bohr単位。通常、原子の価電子半径程度を指定。e.g. H:0.6 程度, C,N,O:1.4 程度)
  • width: 上記のrad前後で、原子の擬波動関数をFermi分布関数にしたがってぼかす際のぼかし幅 (0.3程度でよい模様)
  • DR,NR: (おそらく)未使用
    (元々は、擬波動関数同士やプロジェクターbetaとの重なり積分を実空間で計算する際の刻み幅・積分範囲のパラメータを意図していた模様)

Revision

  • 2007-05-14 (月) 17:24:59 I. Hamada
  • 2006-12-21 (木) 12:40:00 I. Hamada (created)
トップ   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS